5bd569873c3bf01cbc9e88cc30a107469adf5fb7 — axect 3 months ago be3134a
Ver 0.16.0 - Peroxide meets Low level
M Cargo.toml => Cargo.toml +1 -1
@@ 1,6 1,6 @@
 [package]
 name = "peroxide"
-version = "0.15.3"
+version = "0.16.0"
 authors = ["axect <edeftg@gmail.com>"]
 description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"
 repository = "https://github.com/Axect/Peroxide"

M README.md => README.md +1 -1
@@ 189,7 189,7 @@ Corresponding to `0.15.3`
       - [reg.rs](src/ml/reg.rs) : Regression tools
   - __numerical__ : To do numerical things
     - [bdf.rs](src/numerical/bdf.rs) : Backward Differentiation Formula
-    - [gauss_legendre.rs](src/numerical/gauss_legendre.rs) : Gauss-Legendre 4th order
+    - [gauss_legendre.rs](src/grave/gauss_legendre.rs) : Gauss-Legendre 4th order
     - [interp.rs](src/numerical/interp.rs) : Interpolation
     - [mod.rs](src/numerical/mod.rs)
     - [newton.rs](src/numerical/newton.rs) : Newton's Method

M RELEASES.md => RELEASES.md +53 -0
@@ 1,3 1,56 @@
+# Release 0.16.0 (2019-09-14)
+
+## Now, peroxide meets low level
+
+### Added
+* `operation::mut_ops::MutMatrix`
+    * `fn col_mut(&mut self, idx: usize) -> Vec<*mut f64>`
+    * `fn row_mut(&mut self, idx: usize) -> Vec<*mut f64>`
+* `structure::Matrix::FP::{col_mut_map, row_mut_map}`
+    
+### Modified
+* `structure::matrix::Matrix`
+    * `fn ptr(&self) -> *const f64`
+    * `fn mut_ptr(&self) -> *mut f64`
+    * `Index for Matrix`
+    * `IndexMut for Matrix`
+
+### Example
+```rust
+fn main() {
+    // ===================================
+    // Low Level
+    // ===================================
+    let mut a = ml_matrix("1 2; 3 4");
+    a.print();
+    //      c[0] c[1]
+    // r[0]    1    2
+    // r[1]    3    4
+
+    unsafe {
+        let mut p: Vec<*mut f64> = a.col_mut(1); // Mutable second column
+        for i in 0 .. p.len() {
+            *p[i] = i as f64;
+        }
+    }
+
+    a.print();
+    //      c[0] c[1]
+    // r[0]    1    0
+    // r[1]    3    1
+    
+    // ===================================
+    // High Level
+    // ===================================
+    let mut b = ml_matrix("1 2 3; 4 5 6");
+    b.col_mut_map(|x| x.normalize());
+    b.print();
+    //        c[0]   c[1]   c[2]
+    // r[0] 0.2425 0.3714 0.4472
+    // r[1] 0.9701 0.9285 0.8944
+}
+```
+
 # Release 0.15.3 (2019-09-12)
 
 * Smart pointer of Vector - `RedoxVector`

M example_data/lm_test.png => example_data/lm_test.png +0 -0

A src/grave/gauss_legendre.rs => src/grave/gauss_legendre.rs +125 -0
@@ 0,0 1,125 @@
+//use numerical::utils::{jacobian_dual};
+//use structure::dual::*;
+//use structure::matrix::*;
+//use structure::vector::*;
+//use util::non_macro::{cat, concat, eye, zeros};
+//
+///// Value of 3f64.sqrt()
+//const SQRT3: f64 = 1.7320508075688772;
+//
+///// Butcher tableau for Gauss_legendre 4th order
+//const GL4: [[f64; 3]; 2] = [
+//    [0.5 - SQRT3 / 6f64, 0.25, 0.25 - SQRT3 / 6f64],
+//    [0.5 + SQRT3 / 6f64, 0.25 + SQRT3 / 6f64, 0.25],
+//];
+//
+///// Gauss-Legendre 4th order method
+/////
+///// # Description
+/////
+///// ## 1. Find `k1, k2`
+///// * `k1 = f(t + p0*h, y + h(p1*k1 + p2*k2))`
+///// * `k2 = f(t + q0*h, y + h(q1*k1 + q2*k2))`
+/////
+///// ## 2. Iteration
+///// * `y_{n+1} = y_n + 0.5*h*(k1 + k2)`
+//pub fn gl4<F>(f: F, t_init: f64, y_init: Vec<f64>, h: f64, rtol: f64, num: usize) -> Matrix
+//where
+//    F: Fn(Dual, Vec<Dual>) -> Vec<Dual> + Copy,
+//{
+//    let mut t = t_init;
+//    let mut y_curr = y_init.clone();
+//    let mut records = zeros(num + 1, y_curr.len() + 1);
+//    records.subs_row(0, &cat(t, y_curr.clone()));
+//
+//    for i in 0..num {
+//        let (k1, k2) = k_newton(f, t, y_curr.clone(), h, rtol);
+//        y_curr = y_curr.add(&k1.fmap(|x| 0.5 * x * h).add(&k2.fmap(|x| 0.5 * x * h)));
+//        t += h;
+//        records.subs_row(i + 1, &cat(t, y_curr.clone()));
+//    }
+//
+//    records
+//}
+//
+///// Newton's Method for find k in GL4
+/////
+///// # Description
+/////
+///// ## 0. Initial Guess by Euler method
+///// * `k1 = f(t, y)`
+///// * `k2 = f(t, y)`
+/////
+///// ## 1. Combine below two equations to one equation
+///// * `k1 = f(t1, y + h(p1*k1 + p2*k2))`
+///// * `k2 = f(t2, y + h(q1*k1 + q2*k2))`
+///// * `k = g(k)`
+/////
+///// ## 2. Obtain Jacobian
+///// * `DG(k^l) = I - Dg(k^l)`
+/////
+///// ## 3. Iteration by Newton's Method
+///// * `k^{l+1} = k^l - DG^{-1}G(k^l)`
+//#[allow(non_snake_case)]
+//pub fn k_newton<F>(f: F, t: f64, y: Vec<f64>, h: f64, rtol: f64) -> (Vec<f64>, Vec<f64>)
+//where
+//    F: Fn(Dual, Vec<Dual>) -> Vec<Dual> + Copy,
+//{
+//    let t1 = dual(t + GL4[0][0] * h, 0.);
+//    let t2 = dual(t + GL4[1][0] * h, 0.);
+//    let tn = dual(t, 0.);
+//    let yn = y.conv_dual();
+//    let n = y.len();
+//
+//    // 0. Initial Guess
+//    let k1_init = f(tn, yn.clone()).values();
+//    let k2_init = f(tn, yn.clone()).values();
+//    let mut k_curr = concat(k1_init.clone(), k2_init.clone());
+//    let mut err = 1f64;
+//
+//    // 1. Combine two functions to one function
+//    let g = |k: Vec<Dual>| -> Vec<Dual> {
+//        let k1 = k.take(n);
+//        let k2 = k.skip(n);
+//        concat(
+//            f(
+//                t1,
+//                yn.add(
+//                    &k1.s_mul(GL4[0][1] * h)
+//                        .add(&k2.s_mul(GL4[0][2]*h)),
+//                ),
+//            ),
+//            f(
+//                t2,
+//                yn.add(
+//                    &k1.s_mul(GL4[1][1] * h)
+//                        .add(&k2.s_mul(GL4[1][2] * h)),
+//                ),
+//            ),
+//        )
+//    };
+//
+//    // 2. Obtain Jacobian
+//    let I = eye(2 * n);
+//
+//    let mut Dg = jacobian_dual(g, &k_curr);
+//    let mut DG = I.clone() - Dg.clone();
+//    let mut DG_inv = DG.inv().unwrap();
+//    let mut G = k_curr.sub(&g.clone()(k_curr.conv_dual()).values());
+//    let mut num_iter: usize = 0;
+//
+//    // 3. Iteration
+//    while err >= rtol && num_iter <= 10 {
+//        let k_prev = k_curr.clone();
+//        let DGG = DG_inv.clone() * G.clone();
+//        k_curr = k_curr.sub(&DGG.col(0));
+//        Dg = jacobian_dual(g, &k_curr);
+//        DG = I.clone() - Dg.clone();
+//        DG_inv = DG.inv().unwrap();
+//        G = k_curr.sub(&g.clone()(k_curr.conv_dual()).values());
+//        err = k_curr.sub(&k_prev).norm();
+//        num_iter += 1;
+//    }
+//
+//    (k_curr.take(n), k_curr.skip(n))
+//}

D src/numerical/gauss_legendre.rs => src/numerical/gauss_legendre.rs +0 -125
@@ 1,125 0,0 @@
-use numerical::utils::{jacobian_dual};
-use structure::dual::*;
-use structure::matrix::*;
-use structure::vector::*;
-use util::non_macro::{cat, concat, eye, zeros};
-
-/// Value of 3f64.sqrt()
-const SQRT3: f64 = 1.7320508075688772;
-
-/// Butcher tableau for Gauss_legendre 4th order
-const GL4: [[f64; 3]; 2] = [
-    [0.5 - SQRT3 / 6f64, 0.25, 0.25 - SQRT3 / 6f64],
-    [0.5 + SQRT3 / 6f64, 0.25 + SQRT3 / 6f64, 0.25],
-];
-
-/// Gauss-Legendre 4th order method
-///
-/// # Description
-///
-/// ## 1. Find `k1, k2`
-/// * `k1 = f(t + p0*h, y + h(p1*k1 + p2*k2))`
-/// * `k2 = f(t + q0*h, y + h(q1*k1 + q2*k2))`
-///
-/// ## 2. Iteration
-/// * `y_{n+1} = y_n + 0.5*h*(k1 + k2)`
-pub fn gl4<F>(f: F, t_init: f64, y_init: Vec<f64>, h: f64, rtol: f64, num: usize) -> Matrix
-where
-    F: Fn(Dual, Vec<Dual>) -> Vec<Dual> + Copy,
-{
-    let mut t = t_init;
-    let mut y_curr = y_init.clone();
-    let mut records = zeros(num + 1, y_curr.len() + 1);
-    records.subs_row(0, &cat(t, y_curr.clone()));
-
-    for i in 0..num {
-        let (k1, k2) = k_newton(f, t, y_curr.clone(), h, rtol);
-        y_curr = y_curr.add(&k1.fmap(|x| 0.5 * x * h).add(&k2.fmap(|x| 0.5 * x * h)));
-        t += h;
-        records.subs_row(i + 1, &cat(t, y_curr.clone()));
-    }
-
-    records
-}
-
-/// Newton's Method for find k in GL4
-///
-/// # Description
-///
-/// ## 0. Initial Guess by Euler method
-/// * `k1 = f(t, y)`
-/// * `k2 = f(t, y)`
-///
-/// ## 1. Combine below two equations to one equation
-/// * `k1 = f(t1, y + h(p1*k1 + p2*k2))`
-/// * `k2 = f(t2, y + h(q1*k1 + q2*k2))`
-/// * `k = g(k)`
-///
-/// ## 2. Obtain Jacobian
-/// * `DG(k^l) = I - Dg(k^l)`
-///
-/// ## 3. Iteration by Newton's Method
-/// * `k^{l+1} = k^l - DG^{-1}G(k^l)`
-#[allow(non_snake_case)]
-pub fn k_newton<F>(f: F, t: f64, y: Vec<f64>, h: f64, rtol: f64) -> (Vec<f64>, Vec<f64>)
-where
-    F: Fn(Dual, Vec<Dual>) -> Vec<Dual> + Copy,
-{
-    let t1 = dual(t + GL4[0][0] * h, 0.);
-    let t2 = dual(t + GL4[1][0] * h, 0.);
-    let tn = dual(t, 0.);
-    let yn = y.conv_dual();
-    let n = y.len();
-
-    // 0. Initial Guess
-    let k1_init = f(tn, yn.clone()).values();
-    let k2_init = f(tn, yn.clone()).values();
-    let mut k_curr = concat(k1_init.clone(), k2_init.clone());
-    let mut err = 1f64;
-
-    // 1. Combine two functions to one function
-    let g = |k: Vec<Dual>| -> Vec<Dual> {
-        let k1 = k.take(n);
-        let k2 = k.skip(n);
-        concat(
-            f(
-                t1,
-                yn.add(
-                    &k1.s_mul(GL4[0][1] * h)
-                        .add(&k2.s_mul(GL4[0][2]*h)),
-                ),
-            ),
-            f(
-                t2,
-                yn.add(
-                    &k1.s_mul(GL4[1][1] * h)
-                        .add(&k2.s_mul(GL4[1][2] * h)),
-                ),
-            ),
-        )
-    };
-
-    // 2. Obtain Jacobian
-    let I = eye(2 * n);
-
-    let mut Dg = jacobian_dual(g, &k_curr);
-    let mut DG = I.clone() - Dg.clone();
-    let mut DG_inv = DG.inv().unwrap();
-    let mut G = k_curr.sub(&g.clone()(k_curr.conv_dual()).values());
-    let mut num_iter: usize = 0;
-
-    // 3. Iteration
-    while err >= rtol && num_iter <= 10 {
-        let k_prev = k_curr.clone();
-        let DGG = DG_inv.clone() * G.clone();
-        k_curr = k_curr.sub(&DGG.col(0));
-        Dg = jacobian_dual(g, &k_curr);
-        DG = I.clone() - Dg.clone();
-        DG_inv = DG.inv().unwrap();
-        G = k_curr.sub(&g.clone()(k_curr.conv_dual()).values());
-        err = k_curr.sub(&k_prev).norm();
-        num_iter += 1;
-    }
-
-    (k_curr.take(n), k_curr.skip(n))
-}

M src/numerical/mod.rs => src/numerical/mod.rs +1 -1
@@ 1,7 1,7 @@
 //! Differential equations & Numerical Analysis tools
 
 //pub mod bdf;
-pub mod gauss_legendre;
+//pub mod gauss_legendre;
 pub mod interp;
 pub mod newton;
 pub mod ode;

M src/numerical/newton.rs => src/numerical/newton.rs +8 -6
@@ 1,27 1,29 @@
 use numerical::utils::jacobian;
 use operation::number::{Number, NumberVector};
-/// Newton-Raphson Method
 use structure::matrix::*;
 use structure::vector::*;
+use MutFP;
 
+/// Newton-Raphson Method
 pub fn newton<F>(init_cond: Vec<f64>, f: F, rtol: f64) -> Vec<f64>
 where
     F: Fn(Vec<Number>) -> Vec<Number> + Copy,
 {
-    let mut x = init_cond;
-    let mut x_next = update(&x, f);
+    let mut x_next = init_cond;
+    let mut x = x_next.clone();
+    update(&mut x_next, f);
     let mut err = (x_next.sub(&x)).norm();
 
     while err >= rtol {
         x = x_next.clone();
-        x_next = update(&x_next, f);
+        update(&mut x_next, f);
         err = (x_next.sub(&x)).norm();
     }
 
     x_next
 }
 
-fn update<F>(xs: &Vec<f64>, f: F) -> Vec<f64>
+fn update<F>(xs: &mut Vec<f64>, f: F)
 where
     F: Fn(Vec<Number>) -> Vec<Number> + Copy,
 {


@@ 29,5 31,5 @@ where
     let pinv_j = j.pseudo_inv().unwrap();
     let fx = f(NumberVector::from_f64_vec(xs.clone())).to_f64_vec();
 
-    xs.sub(&(pinv_j * fx).col(0))
+    xs.mut_zip_with(|x,y| x - y, &(pinv_j * fx).col(0))
 }

M src/numerical/ode.rs => src/numerical/ode.rs +3 -3
@@ 210,7 210,7 @@ use operation::extra_ops::Real;
 use operation::mut_ops::MutFP;
 use std::collections::HashMap;
 use structure::dual::Dual;
-use structure::matrix::{Matrix, Row, FP, LinearAlgebra};
+use structure::matrix::{Matrix, FP, LinearAlgebra};
 use structure::vector::FPVector;
 use numerical::utils::jacobian_dual;
 use util::non_macro::{cat, concat, zeros, eye};


@@ 474,7 474,7 @@ impl ODE for ExplicitODE {
                         break;
                     }
                 }
-                return result.take(key, Row);
+                return result.take_row(key);
             }
             _ => {
                 for i in 1..self.times + 1 {


@@ 766,7 766,7 @@ impl ODE for ImplicitODE {
                         break;
                     }
                 }
-                return result.take(key, Row);
+                return result.take_row(key);
             }
             _ => {
                 for i in 1..self.times + 1 {

M src/operation/mut_ops.rs => src/operation/mut_ops.rs +59 -0
@@ 1,3 1,5 @@
+use ::{Matrix, Shape};
+
 pub trait MutFP {
     type Scalar;
     fn mut_map<F>(&mut self, f: F)


@@ 29,3 31,60 @@ impl MutFP for Vec<f64> {
         }
     }
 }
+
+pub trait MutMatrix {
+    unsafe fn col_mut(&mut self, idx: usize) -> Vec<*mut f64>;
+    unsafe fn row_mut(&mut self, idx: usize) -> Vec<*mut f64>;
+}
+
+impl MutMatrix for Matrix {
+    unsafe fn col_mut(&mut self, idx: usize) -> Vec<*mut f64> {
+        assert!(idx < self.col, "Index out of range");
+        match self.shape {
+            Shape::Col => {
+                let mut v: Vec<*mut f64> = Vec::with_capacity(self.row);
+                v.set_len(self.row);
+                let start_idx = idx * self.row;
+                let p = self.mut_ptr();
+                for (i, j) in (start_idx .. start_idx + v.len()).enumerate() {
+                    v[i] = p.add(j);
+                }
+                v
+            }
+            Shape::Row => {
+                let mut v: Vec<*mut f64> = Vec::with_capacity(self.row);
+                v.set_len(self.row);
+                let p = self.mut_ptr();
+                for i in 0 .. v.len() {
+                    v[i] = p.add(idx + i * self.col);
+                }
+                v
+            }
+        }
+    }
+
+    unsafe fn row_mut(&mut self, idx: usize) -> Vec<*mut f64> {
+        assert!(idx < self.row, "Index out of range");
+        match self.shape {
+            Shape::Row => {
+                let mut v: Vec<*mut f64> = Vec::with_capacity(self.col);
+                v.set_len(self.col);
+                let start_idx = idx * self.col;
+                let p = self.mut_ptr();
+                for (i, j) in (start_idx .. start_idx + v.len()).enumerate() {
+                    v[i] = p.add(j);
+                }
+                v
+            }
+            Shape::Col => {
+                let mut v: Vec<*mut f64> = Vec::with_capacity(self.col);
+                v.set_len(self.col);
+                let p = self.mut_ptr();
+                for i in 0 .. v.len() {
+                    v[i] = p.add(idx + i * self.row);
+                }
+                v
+            }
+        }
+    }
+}<
\ No newline at end of file

M src/structure/matrix.rs => src/structure/matrix.rs +150 -98
@@ 447,9 447,12 @@ use std::convert;
 pub use std::error::Error;
 use std::fmt;
 use std::ops::{Add, Index, IndexMut, Mul, Neg, Sub};
+use std::iter;
+use std::slice::{Iter, IterMut};
 #[allow(unused_imports)]
 use structure::vector::*;
 use util::useful::*;
+use MutMatrix;
 
 pub type Perms = Vec<(usize, usize)>;
 


@@ 624,6 627,14 @@ impl PartialEq for Matrix {
 /// Main matrix structure
 #[allow(dead_code)]
 impl Matrix {
+    pub fn ptr(&self) -> *const f64 {
+        &self.data[0] as *const f64
+    }
+
+    pub fn mut_ptr(&mut self) -> *mut f64 {
+        &mut self.data[0] as *mut f64
+    }
+
     /// Change Bindings
     ///
     /// `Row` -> `Col` or `Col` -> `Row`


@@ 653,12 664,7 @@ impl Matrix {
                     data[i] = ref_data[s];
                 }
                 data[l] = ref_data[l];
-                Self {
-                    data: data.clone(),
-                    row: r,
-                    col: c,
-                    shape: Col,
-                }
+                matrix(data.clone(),r,c,Col)
             }
             Col => {
                 for i in 0..l {


@@ 666,12 672,7 @@ impl Matrix {
                     data[i] = ref_data[s];
                 }
                 data[l] = ref_data[l];
-                Self {
-                    data: data.clone(),
-                    row: r,
-                    col: c,
-                    shape: Row,
-                }
+                matrix(data.clone(),r,c,Row)
             }
         }
     }


@@ 701,13 702,13 @@ impl Matrix {
             let part = if r <= 10 {
                 key_row = r;
                 key_col = 100;
-                self.take(100, Col)
+                self.take_col(100)
             } else if c <= 10 {
                 key_row = 100;
                 key_col = c;
-                self.take(100, Row)
+                self.take_row(100)
             } else {
-                self.take(20, Row).take(20, Col)
+                self.take_row(20).take_col(20)
             };
             return format!(
                 "Result is too large to print - {}x{}\nonly print {}x{} parts:\n{}",


@@ 1004,6 1005,25 @@ impl Matrix {
         Ok(m)
     }
 
+    /// Should check shape
+    pub fn subs(&mut self, idx: usize, v: &Vec<f64>) {
+        let p = &mut self.mut_ptr();
+        match self.shape {
+            Row => {
+                let c = self.col;
+                unsafe {
+                    p.add(idx * c).copy_from(v.as_ptr(), c);
+                }
+            }
+            Col => {
+                let r = self.row;
+                unsafe {
+                    p.add(idx * r).copy_from(v.as_ptr(), r);
+                }
+            }
+        }
+    }
+
     /// Substitute Col
     pub fn subs_col(&mut self, idx: usize, v: &Vec<f64>) {
         for i in 0..self.row {


@@ 1888,13 1908,7 @@ impl<'a, 'b> Mul<&'b Matrix> for &'a Vector {
     fn mul(self, other: &'b Matrix) -> Self::Output {
         assert_eq!(self.len(), other.row);
         let l = self.len();
-        Matrix {
-            data: self.clone(),
-            row: 1,
-            col: l,
-            shape: Col,
-        }
-        .mul(other.clone())
+        matrix(self.clone(), 1, l, Col).mul(other.clone())
     }
 }
 


@@ 1914,11 1928,21 @@ impl Index<(usize, usize)> for Matrix {
     type Output = f64;
 
     fn index(&self, pair: (usize, usize)) -> &f64 {
+        let p = self.ptr();
         let i = pair.0;
         let j = pair.1;
+        assert!(i < self.row && j < self.col, "Index out of range");
         match self.shape {
-            Row => &self.data[i * self.col + j],
-            Col => &self.data[i + j * self.row],
+            Row => {
+                unsafe {
+                    &*p.add(i * self.col + j)
+                }
+            }   
+            Col => {
+                unsafe {
+                    &*p.add(i + j * self.row)
+                }  
+            },
         }
     }
 }


@@ 1942,14 1966,20 @@ impl IndexMut<(usize, usize)> for Matrix {
         let j = pair.1;
         let r = self.row;
         let c = self.col;
+        assert!(i < self.row && j < self.col, "Index out of range");
+        let mut p = self.mut_ptr();
         match self.shape {
             Row => {
                 let idx = i * c + j;
-                &mut self.data[idx]
+                unsafe {
+                    &mut *p.add(idx)
+                }
             }
             Col => {
                 let idx = i + j * r;
-                &mut self.data[idx]
+                unsafe {
+                    &mut *p.add(idx)
+                }
             }
         }
     }


@@ 1959,8 1989,10 @@ impl IndexMut<(usize, usize)> for Matrix {
 // Functional Programming Tools (Hand-written)
 // =============================================================================
 pub trait FP {
-    fn take(&self, n: usize, shape: Shape) -> Matrix;
-    fn skip(&self, n: usize, shape: Shape) -> Matrix;
+    fn take_row(&self, n: usize) -> Matrix;
+    fn take_col(&self, n: usize) -> Matrix;
+    fn skip_row(&self, n: usize) -> Matrix;
+    fn skip_col(&self, n: usize) -> Matrix;
     fn fmap<F>(&self, f: F) -> Matrix
     where
         F: Fn(f64) -> f64;


@@ 1970,6 2002,12 @@ pub trait FP {
     fn row_map<F>(&self, f: F) -> Matrix
     where
         F: Fn(Vec<f64>) -> Vec<f64>;
+    fn col_mut_map<F>(&mut self, f: F)
+    where
+        F: Fn(Vec<f64>) -> Vec<f64>;
+    fn row_mut_map<F>(&mut self, f: F)
+    where
+        F: Fn(Vec<f64>) -> Vec<f64>;
     fn reduce<F, T>(&self, init: T, f: F) -> f64
     where
         F: Fn(f64, f64) -> f64,


@@ 1980,80 2018,74 @@ pub trait FP {
 }
 
 impl FP for Matrix {
-    fn take(&self, n: usize, shape: Shape) -> Self {
-        match shape {
+    fn take_row(&self, n: usize) -> Self {
+        if n >= self.row {
+            return self.clone();
+        }
+        match self.shape {
             Row => {
-                if n >= self.row {
-                    return self.clone();
-                }
-                match self.shape {
-                    Row => {
-                        let new_data = self
-                            .data
-                            .clone()
-                            .into_iter()
-                            .take(n * self.col)
-                            .collect::<Vec<f64>>();
-                        matrix(new_data, n, self.col, Row)
-                    }
-                    Col => {
-                        let mut temp_data: Vec<f64> = Vec::new();
-                        for i in 0..n {
-                            temp_data.extend(self.row(i));
-                        }
-                        matrix(temp_data, n, self.col, Row)
-                    }
-                }
+                let new_data = self
+                    .data
+                    .clone()
+                    .into_iter()
+                    .take(n * self.col)
+                    .collect::<Vec<f64>>();
+                matrix(new_data, n, self.col, Row)
             }
             Col => {
-                if n >= self.col {
-                    return self.clone();
-                }
-                match self.shape {
-                    Col => {
-                        let new_data = self
-                            .data
-                            .clone()
-                            .into_iter()
-                            .take(n * self.row)
-                            .collect::<Vec<f64>>();
-                        matrix(new_data, self.row, n, Col)
-                    }
-                    Row => {
-                        let mut temp_data: Vec<f64> = Vec::new();
-                        for i in 0..n {
-                            temp_data.extend(self.col(i));
-                        }
-                        matrix(temp_data, self.row, n, Col)
-                    }
+                let mut temp_data: Vec<f64> = Vec::new();
+                for i in 0..n {
+                    temp_data.extend(self.row(i));
                 }
+                matrix(temp_data, n, self.col, Row)
             }
         }
     }
 
-    fn skip(&self, n: usize, shape: Shape) -> Self {
-        match shape {
-            Row => {
-                assert!(n < self.row, "Skip range is larger than row of matrix");
-
-                let mut temp_data: Vec<f64> = Vec::new();
-                for i in n..self.row {
-                    temp_data.extend(self.row(i));
-                }
-                matrix(temp_data, self.row - n, self.col, Row)
-            }
+    fn take_col(&self, n: usize) -> Self {
+        if n >= self.col {
+            return self.clone();
+        }
+        match self.shape {
             Col => {
-                assert!(n < self.col, "Skip range is larger than col of matrix");
-
+                let new_data = self
+                    .data
+                    .clone()
+                    .into_iter()
+                    .take(n * self.row)
+                    .collect::<Vec<f64>>();
+                matrix(new_data, self.row, n, Col)
+            }
+            Row => {
                 let mut temp_data: Vec<f64> = Vec::new();
-                for i in n..self.col {
+                for i in 0..n {
                     temp_data.extend(self.col(i));
                 }
-                matrix(temp_data, self.row, self.col - n, Col)
+                matrix(temp_data, self.row, n, Col)
             }
         }
     }
 
+    fn skip_row(&self, n: usize) -> Self {
+        assert!(n < self.row, "Skip range is larger than row of matrix");
+
+        let mut temp_data: Vec<f64> = Vec::new();
+        for i in n..self.row {
+            temp_data.extend(self.row(i));
+        }
+        matrix(temp_data, self.row - n, self.col, Row)
+    }
+
+    fn skip_col(&self, n: usize) -> Self {
+        assert!(n < self.col, "Skip range is larger than col of matrix");
+
+        let mut temp_data: Vec<f64> = Vec::new();
+        for i in n..self.col {
+            temp_data.extend(self.col(i));
+        }
+        matrix(temp_data, self.row, self.col - n, Col)
+    }
+
     fn fmap<F>(&self, f: F) -> Matrix
     where
         F: Fn(f64) -> f64,


@@ 2071,12 2103,7 @@ impl FP for Matrix {
     where
         F: Fn(Vec<f64>) -> Vec<f64>,
     {
-        let mut result = Matrix {
-            data: vec![0f64; self.row * self.col],
-            row: self.row,
-            col: self.col,
-            shape: Col,
-        };
+        let mut result = matrix(vec![0f64; self.row * self.col], self.row, self.col, Col);
 
         for i in 0..self.row {
             result.subs_col(i, &f(self.col(i)));


@@ 2089,12 2116,7 @@ impl FP for Matrix {
     where
         F: Fn(Vec<f64>) -> Vec<f64>,
     {
-        let mut result = Matrix {
-            data: vec![0f64; self.row * self.col],
-            row: self.row,
-            col: self.col,
-            shape: Row,
-        };
+        let mut result = matrix(vec![0f64; self.row * self.col], self.row, self.col, Row);
 
         for i in 0..self.col {
             result.subs_row(i, &f(self.row(i)));


@@ 2103,6 2125,36 @@ impl FP for Matrix {
         result
     }
 
+    fn col_mut_map<F>(&mut self, f: F)
+    where
+        F: Fn(Vec<f64>) -> Vec<f64>,
+    {
+        for i in 0 .. self.col {
+            unsafe {
+                let mut p = self.col_mut(i);
+                let fv = f(self.col(i));
+                for j in 0..p.len() {
+                    *p[j] = fv[j];
+                }
+            }
+        }
+    }
+
+    fn row_mut_map<F>(&mut self, f: F)
+    where
+        F: Fn(Vec<f64>) -> Vec<f64>,
+    {
+        for i in 0 .. self.col {
+            unsafe {
+                let mut p = self.row_mut(i);
+                let fv = f(self.row(i));
+                for j in 0 .. p.len() {
+                    *p[j] = fv[j];
+                }
+            }
+        }
+    }
+
     fn reduce<F, T>(&self, init: T, f: F) -> f64
     where
         F: Fn(f64, f64) -> f64,

M src/structure/vector.rs => src/structure/vector.rs +2 -2
@@ 402,8 402,8 @@ impl FPVector for Vector {
     fn skip(&self, n: usize) -> Vector {
         let l = self.len();
         let mut v = vec![0f64; l - n];
-        for i in n..l {
-            v[i - n] = self[i];
+        for (i, j) in (n..l).enumerate() {
+            v[i] = self[j];
         }
         return v;
     }

M src/util/api.rs => src/util/api.rs +3 -18
@@ 40,12 40,7 @@ impl MATLAB for Matrix {
                     .collect::<Vec<f64>>()
             })
             .collect::<Vec<f64>>();
-        Matrix {
-            data: data,
-            row: r,
-            col: c,
-            shape: Row,
-        }
+        matrix(data, r, c, Row)
     }
 }
 


@@ 63,12 58,7 @@ impl PYTHON for Matrix {
                 data[idx] = v[i][j].into();
             }
         }
-        Matrix {
-            data: data,
-            row: r,
-            col: c,
-            shape: Row,
-        }
+        matrix(data, r, c, Row)
     }
 }
 


@@ 77,11 67,6 @@ impl R for Matrix {
     where
         T: convert::Into<f64>,
     {
-        Matrix {
-            data: v.into_iter().map(|t| t.into()).collect::<Vec<f64>>(),
-            row: x,
-            col: y,
-            shape: shape,
-        }
+        matrix(v.into_iter().map(|t| t.into()).collect::<Vec<f64>>(), x, y, shape)
     }
 }