~axect/Peroxide

e57c0ccac329b2a41efa59b8cb4bf809bc88def4 — axect 7 months ago 4cedd6b
Ver 0.13.0 - Little bit optimization
M Cargo.toml => Cargo.toml +2 -4
@@ 1,6 1,6 @@
[package]
name = "peroxide"
version = "0.12.4"
version = "0.13.0"
authors = ["axect <edeftg@gmail.com>"]
description = "Pure rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"
repository = "https://github.com/Axect/Peroxide"


@@ 9,7 9,7 @@ categories = ["science"]
readme = "README.md"
documentation = "https://docs.rs/peroxide/"
keywords = ["Numeric", "R", "MATLAB", "Python", "Scientific"]
exclude = ["example_data/", "src/bin/"]
exclude = ["example_data/", "src/bin/", "benches/"]

[badges]
maintenance = { status = "actively-developed" }


@@ 18,8 18,6 @@ maintenance = { status = "actively-developed" }
csv = "1"
rand = "0.7"
special = "0.8"
serde = "1"
serde-pickle = "0.5"
pyo3 = "0.7"

[package.metadata.docs.rs]

M RELEASES.md => RELEASES.md +14 -0
@@ 1,3 1,17 @@
# Release 0.13.0 (2019-08-11)

* **[Important]** Remove inefficient things
    * Remove `util/pickle.rs`
    * Remove `serde`, `serde-pickle` dependencies
* Optimize original code
    * Matrix
        * `change_shape`
        * `col(usize), row(usize)`
        * `diag`
        * `subs_col, subs_row`
        * `Add<Matrix>`
        * `Mul<Matrix>` - Block matrix multiplication

# Release 0.12.4 (2019-08-09)

* Add `normalize` for `Vec<f64>`

A benches/data/add_test_12.md => benches/data/add_test_12.md +3 -0
@@ 0,0 1,3 @@
| Command | Mean [s] | Min [s] | Max [s] | Relative |
|:---|---:|---:|---:|---:|
| `./target/release/add_test` | 3.752 ± 0.219 | 3.494 | 4.176 | 1.0 |

A benches/data/add_test_13.md => benches/data/add_test_13.md +3 -0
@@ 0,0 1,3 @@
| Command | Mean [s] | Min [s] | Max [s] | Relative |
|:---|---:|---:|---:|---:|
| `./target/release/add_test` | 2.581 ± 0.107 | 2.475 | 2.784 | 1.0 |

A benches/data/col_test_12_col.md => benches/data/col_test_12_col.md +3 -0
@@ 0,0 1,3 @@
| Command | Mean [s] | Min [s] | Max [s] | Relative |
|:---|---:|---:|---:|---:|
| `./target/release/col_test` | 1.009 ± 0.003 | 1.005 | 1.016 | 1.0 |

A benches/data/col_test_13_col.md => benches/data/col_test_13_col.md +3 -0
@@ 0,0 1,3 @@
| Command | Mean [ms] | Min [ms] | Max [ms] | Relative |
|:---|---:|---:|---:|---:|
| `./target/release/col_test` | 354.2 ± 1.3 | 352.9 | 356.3 | 1.0 |

A benches/data/col_test_13_idx.md => benches/data/col_test_13_idx.md +3 -0
@@ 0,0 1,3 @@
| Command | Mean [ms] | Min [ms] | Max [ms] | Relative |
|:---|---:|---:|---:|---:|
| `./target/release/col_test` | 352.7 ± 1.1 | 351.3 | 354.5 | 1.0 |

A benches/data/idx_test_12.md => benches/data/idx_test_12.md +3 -0
@@ 0,0 1,3 @@
| Command | Mean [ms] | Min [ms] | Max [ms] | Relative |
|:---|---:|---:|---:|---:|
| `./target/release/idx_test` | 352.0 ± 1.8 | 350.6 | 356.4 | 1.0 |

A benches/data/row_test_12_row.md => benches/data/row_test_12_row.md +3 -0
@@ 0,0 1,3 @@
| Command | Mean [ms] | Min [ms] | Max [ms] | Relative |
|:---|---:|---:|---:|---:|
| `./target/release/row_test` | 677.8 ± 4.0 | 674.6 | 688.4 | 1.0 |

A benches/data/row_test_13_idx.md => benches/data/row_test_13_idx.md +3 -0
@@ 0,0 1,3 @@
| Command | Mean [ms] | Min [ms] | Max [ms] | Relative |
|:---|---:|---:|---:|---:|
| `./target/release/row_test` | 352.0 ± 1.0 | 350.9 | 354.2 | 1.0 |

A benches/data/row_test_13_row.md => benches/data/row_test_13_row.md +3 -0
@@ 0,0 1,3 @@
| Command | Mean [ms] | Min [ms] | Max [ms] | Relative |
|:---|---:|---:|---:|---:|
| `./target/release/row_test` | 357.5 ± 4.2 | 353.6 | 365.5 | 1.0 |

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

M src/bin/jacobian.rs => src/bin/jacobian.rs +2 -2
@@ 3,7 3,7 @@ use peroxide::*;
use Number::D;

fn main() {
    let x = vec![D(dual(1f64, 1f64)), D(dual(2,1))];
    let x = vec![D(dual(1f64, 1f64)), D(dual(2, 1))];
    jacobian(f, x.to_f64_vec()).print();
}



@@ 18,4 18,4 @@ fn g(x: Number) -> Number {
    } else {
        x.powi(3)
    }
}
\ No newline at end of file
}

M src/bin/lda_ls.rs => src/bin/lda_ls.rs +2 -2
@@ 37,8 37,8 @@ fn main() {
    let X_tilde = cbind(X_bias, X.clone());

    let mut T = zeros(2 * N, 2);
    T.subs_col(0, concat(vec![0f64; N], vec![1f64; N]));
    T.subs_col(1, concat(vec![1f64; N], vec![0f64; N]));
    T.subs_col(0, &concat(vec![0f64; N], vec![1f64; N]));
    T.subs_col(1, &concat(vec![1f64; N], vec![0f64; N]));

    x1.print();
    x2.print();

M src/bin/lm.rs => src/bin/lm.rs +5 -2
@@ 1,3 1,4 @@
#[macro_use]
extern crate peroxide;
use peroxide::*;



@@ 21,8 22,10 @@ fn main() {
}

fn f(domain: &Vec<f64>, p: Vec<Number>) -> Vec<Number> {
    domain.clone().into_iter()
    domain
        .clone()
        .into_iter()
        .map(|t| Number::from_f64(t))
        .map(|t| p[0] * (-t / p[1]).exp() + p[2] * t * (-t / p[3]).exp())
        .collect()
}
\ No newline at end of file
}

M src/bin/lm2.rs => src/bin/lm2.rs +6 -4
@@ 30,7 30,7 @@ fn main() {
    // =========================================================================
    // Update
    // =========================================================================
    for _i in 0 .. 30 {
    for _i in 0..30 {
        let b = &jtj + &(lambda * jtj.to_diag());
        let h: Matrix;
        match b.inv() {


@@ 42,14 42,15 @@ fn main() {
        let p_temp = &p_mat + &h;
        let p_temp_num: Vec<Number> = NumberVector::from_f64_vec(p_temp.col(0));
        let y_hat_temp = f(&t, p_temp_num.clone()).to_f64_vec().to_matrix();
        let rho = ((&y - &y_hat_temp).t() * (&y - &y_hat_temp))[(0,0)] / (&h.t() * &(lambda * (&jtj.to_diag() * &h) + j.t() * (&y - &y_hat)))[(0,0)];
        let rho = ((&y - &y_hat_temp).t() * (&y - &y_hat_temp))[(0, 0)]
            / (&h.t() * &(lambda * (&jtj.to_diag() * &h) + j.t() * (&y - &y_hat)))[(0, 0)];
        if rho > 0f64 {
            p_mat = p_temp;
            p = NumberVector::from_f64_vec(p_mat.col(0));
            y_hat = y_hat_temp;
            j = jacobian(|p| f(&t, p), p.to_f64_vec());
            jtj = &j.t() * &j;
            lambda = lambda * max(vec![1f64/3f64, 1f64 - (2f64*rho - 1f64).powi(3)]);
            lambda = lambda * max(vec![1f64 / 3f64, 1f64 - (2f64 * rho - 1f64).powi(3)]);
            nu = 2f64;
        } else {
            lambda = lambda * nu;


@@ 60,7 61,8 @@ fn main() {
}

fn f(t: &Domain, p: Parameter) -> Parameter {
    t.clone().into_iter()
    t.clone()
        .into_iter()
        .map(|x| p[0] * (-x / p[1]).exp() + p[2] * x * (-x / p[3]).exp())
        .collect()
}

M src/bin/lm_test.rs => src/bin/lm_test.rs +8 -4
@@ 1,3 1,4 @@
#[macro_use]
extern crate peroxide;
use peroxide::*;



@@ 15,7 16,8 @@ fn main() {
    let data = hstack!(x.clone(), y.clone());

    let mut opt = Optimizer::new(data, quad);
    let p = opt.set_init_param(n_init)
    let p = opt
        .set_init_param(n_init)
        .set_max_iter(50)
        .set_method(LevenbergMarquardt)
        .optimize();


@@ 32,12 34,14 @@ fn main() {
        .set_title("Test ($y=x^2$ with noise)")
        .set_path("example_data/lm_test.png")
        .set_marker(vec![Point, Line])
        .savefig().expect("Can't draw a plot");
        .savefig()
        .expect("Can't draw a plot");
}

fn quad(x: &Vec<f64>, n: Vec<Number>) -> Vec<Number> {
    x.clone().into_iter()
    x.clone()
        .into_iter()
        .map(|t| Number::from_f64(t))
        .map(|t| t.powf(n[0]))
        .collect()
}
\ No newline at end of file
}

A src/bin/lu.rs => src/bin/lu.rs +7 -0
@@ 0,0 1,7 @@
extern crate peroxide;
use peroxide::*;

fn main() {
    let a = rand(1000, 1000);
    let pqlu = a.lu().unwrap();
}

A src/bin/matmul1.rs => src/bin/matmul1.rs +7 -0
@@ 0,0 1,7 @@
extern crate peroxide;
use peroxide::*;

fn main() {
    let a = matrix(c!(1,2,3,4,5,6), 3, 2, Row);
    (&a * &a.t()).print();
}
\ No newline at end of file

A src/bin/matmul2.rs => src/bin/matmul2.rs +24 -0
@@ 0,0 1,24 @@
extern crate peroxide;
use peroxide::*;

fn main() {
    let a = zeros(1200, 1200);
    let _b = matmul(&a, &a);
}

fn matmul(a: &Matrix, b: &Matrix) -> Matrix {
    match (a.row, a.col) {
        (p, q) if p <= 100 && q <= 100 => a * b,
        _ => {
            let (a1, a2, a3, a4) = a.block();
            let (b1, b2, b3, b4) = b.block();

            let m1 = matmul(&a1,&b1) + matmul(&a2,&b3);
            let m2 = matmul(&a1,&b2) + matmul(&a2,&b4);
            let m3 = matmul(&a3,&b1) + matmul(&a4,&b3);
            let m4 = matmul(&a3,&b2) + matmul(&a4,&b4);

            combine(m1, m2, m3, m4)
        }
    }
}
\ No newline at end of file

A src/bin/matmul3.rs => src/bin/matmul3.rs +41 -0
@@ 0,0 1,41 @@
extern crate peroxide;
use peroxide::*;

fn main() {
    let a = zeros(2048, 2048);
    let _b = strassen(&a, &a);
}

/// Strassen's algorithm
/// # Assumption for convenience
/// * All matrices are square matrix
/// * Except base case, size of all matrices are power of 2
fn strassen(a: &Matrix, b: &Matrix) -> Matrix {
    // Guarantee same length
    assert_eq!(a.col, b.row);

    // Base case
    match (a.row, a.col) {
        (p, q) if p <= 100 && q <= 100 => a * b,
        _ => {
            let (a1, a2, a3, a4) = a.block(); // Split a
            let (b1, b2, b3, b4) = b.block(); // Split b

            // Do Strassen's algorithm
            let p1 = strassen(&a1, &(&b2 - &b4));          // A (F - H)
            let p2 = strassen(&(&a1 + &a2), &b4);          // (A + B) H
            let p3 = strassen(&(&a3 + &a4), &b1);          // (C + D) E
            let p4 = strassen(&a4, &(&b3 - &b1));          // D (G - E)
            let p5 = strassen(&(&a1 + &a4), &(&b1 + &b4)); // (A + D)(E + H)
            let p6 = strassen(&(&a2 - &a4), &(&b3 + &b4)); // (B - D)(G + H)
            let p7 = strassen(&(&a1 - &a3), &(&b1 + &b2)); // (A - C)(E + F)

            let part1 = &(&p5 + &p4) + &(&p6 - &p2); // P5 + P4 - (P2 - P6)
            let part2 = &p1 + &p2;                       // P1 + P2
            let part3 = &p3 + &p4;                       // P3 + P4
            let part4 = &(&p1 - &p3) + &(&p5 - &p7); // P1 + P5 - (P3 + P7)

            combine(part1, part2, part3, part4)
        }
    }
}
\ No newline at end of file

M src/bin/ode_test_lorenz.rs => src/bin/ode_test_lorenz.rs +5 -5
@@ 19,11 19,11 @@ fn main() {

    let results2 = ex_test2.integrate();

    let mut wt = SimpleWriter::new();
    wt.set_path("example_data/lorenz.pickle")
        .insert_matrix(results)
        .insert_matrix(results2)
        .write_pickle();
    // let mut wt = SimpleWriter::new();
    // wt.set_path("example_data/lorenz.pickle")
    //     .insert_matrix(results)
    //     .insert_matrix(results2)
    //     .write_pickle();
}

fn f(st: &mut State<f64>) {

M src/bin/ode_test_rk4.rs => src/bin/ode_test_rk4.rs +4 -4
@@ 14,10 14,10 @@ fn main() {

    let result = ode_solver.integrate();

    let mut st = SimpleWriter::new();
    st.set_path("example_data/rk4_test.pickle")
        .insert_matrix(result)
        .write_pickle();
    // let mut st = SimpleWriter::new();
    // st.set_path("example_data/rk4_test.pickle")
    //     .insert_matrix(result)
    //     .write_pickle();
}

fn test_fn(st: &mut State<f64>) {

M src/bin/optimize.rs => src/bin/optimize.rs +55 -22
@@ 11,7 11,12 @@ fn main() {
    let noise = Normal(0., 0.5);
    let p_true: Vec<Number> = NumberVector::from_f64_vec(vec![20f64, 10f64, 1f64, 50f64]);
    let real = f(p_true.clone()).to_f64_vec();
    let y = matrix(zip_with(|x, y| x + y, &real, &noise.sample(size)), size, 1, Col);
    let y = matrix(
        zip_with(|x, y| x + y, &real, &noise.sample(size)),
        size,
        1,
        Col,
    );
    y.print();
    real.print();



@@ 22,22 27,37 @@ fn main() {
    // Gradient Descent
    let mut p_gd = matrix(p_init.clone(), 4, 1, Col);
    let mut j_gd = j_init.clone();
    let mut y_hat_gd = matrix(f(NumberVector::from_f64_vec(p_init.clone())).to_f64_vec(), size, 1, Col);

    for _i in 0 .. 30 {
    let mut y_hat_gd = matrix(
        f(NumberVector::from_f64_vec(p_init.clone())).to_f64_vec(),
        size,
        1,
        Col,
    );

    for _i in 0..30 {
        let h = 0.02 * j_gd.t() * (&y - &y_hat_gd);
        p_gd = &p_gd + &h;
        j_gd = jacobian(f, p_gd.data.clone());
        y_hat_gd = matrix(f(NumberVector::from_f64_vec(p_gd.data.clone())).to_f64_vec(), size, 1, Col);
        y_hat_gd = matrix(
            f(NumberVector::from_f64_vec(p_gd.data.clone())).to_f64_vec(),
            size,
            1,
            Col,
        );
    }

    p_gd.print();

    // Gauss_Newton
    let mut p_gn = matrix(p_init.clone(), 4, 1, Col);
    let mut y_hat_gn = matrix(f(NumberVector::from_f64_vec(p_init.clone())).to_f64_vec(), size, 1, Col);
    let mut y_hat_gn = matrix(
        f(NumberVector::from_f64_vec(p_init.clone())).to_f64_vec(),
        size,
        1,
        Col,
    );
    let mut j_gn = j_init.clone();
    for _i in 0 .. 10 {
    for _i in 0..10 {
        let h_gn: Matrix;
        match j_gn.pseudo_inv() {
            Some(W) => h_gn = W * (&y - &y_hat_gn),


@@ 45,7 65,12 @@ fn main() {
        }
        p_gn = &p_gn + &h_gn;
        j_gn = jacobian(f, p_gn.data.clone());
        y_hat_gn = matrix(f(NumberVector::from_f64_vec(p_gn.data.clone())).to_f64_vec(), size, 1, Col);
        y_hat_gn = matrix(
            f(NumberVector::from_f64_vec(p_gn.data.clone())).to_f64_vec(),
            size,
            1,
            Col,
        );
    }

    p_gn.print();


@@ 53,19 78,21 @@ fn main() {
    // Levenberg-Marquardt
    let (lambda_0, _eps1, _eps2) = (1e-2, 1e-6, 1e-6);
    let mut p_lm = p_init.to_matrix();
    let mut y_hat_lm = f(NumberVector::from_f64_vec(p_lm.data.clone())).to_f64_vec().to_matrix();
    let mut y_hat_lm = f(NumberVector::from_f64_vec(p_lm.data.clone()))
        .to_f64_vec()
        .to_matrix();
    let mut j_lm = j_init.clone();
    let mut A = &j_lm.t() * &j_lm;
    let mut lambda = lambda_0 * max(A.diag());
    let mut chi2 = ((&y - &y_hat_lm).t() * (&y - &y_hat_lm))[(0,0)];
    let mut chi2 = ((&y - &y_hat_lm).t() * (&y - &y_hat_lm))[(0, 0)];
    let mut nu = 2f64;
    
    for _i in 0 .. 30 {

    for _i in 0..30 {
        let h_lm: Matrix;
        let mut A_diag = eye(A.row);
        let A_ref = A.diag();
        for i in 0 .. A.row {
            A_diag[(i,i)] = A_ref[i];
        for i in 0..A.row {
            A_diag[(i, i)] = A_ref[i];
        }
        match (A.clone() + lambda * A_diag.clone()).inv() {
            Some(B) => h_lm = B * j_lm.t() * (&y - &y_hat_lm),


@@ 75,10 102,13 @@ fn main() {
        let p_lm_temp = &p_lm + &h_lm;
        let j_lm_temp = jacobian(f, p_lm.data.clone());
        let A_temp = &j_lm.t() * &j_lm;
        let y_hat_temp = f(NumberVector::from_f64_vec(p_lm_temp.data.clone())).to_f64_vec().to_matrix();
        let chi2_temp = ((&y - &y_hat_temp).t() * (&y - &y_hat_temp))[(0,0)];
        let y_hat_temp = f(NumberVector::from_f64_vec(p_lm_temp.data.clone()))
            .to_f64_vec()
            .to_matrix();
        let chi2_temp = ((&y - &y_hat_temp).t() * (&y - &y_hat_temp))[(0, 0)];

        let rho = (chi2 - chi2_temp) / (h_lm.t() * (lambda * A_diag * h_lm.clone() + j_lm.t() * (&y - &y_hat_lm)))[(0,0)];
        let rho = (chi2 - chi2_temp)
            / (h_lm.t() * (lambda * A_diag * h_lm.clone() + j_lm.t() * (&y - &y_hat_lm)))[(0, 0)];

        rho.print();



@@ 89,7 119,7 @@ fn main() {
            y_hat_lm = y_hat_temp;
            chi2 = chi2_temp;

            lambda = lambda * max(c!(1f64/3f64, 1f64 - (2f64*rho - 1f64).powi(3)));
            lambda = lambda * max(c!(1f64 / 3f64, 1f64 - (2f64 * rho - 1f64).powi(3)));
            nu = 2f64;
        } else {
            lambda = lambda * nu;


@@ 100,7 130,7 @@ fn main() {
    p_lm.print();

    // Plot
    let p_x = seq(0, (SIZE-1) as i32, 1);
    let p_x = seq(0, (SIZE - 1) as i32, 1);
    let p_y = y.data;
    let f_y = f(NumberVector::from_f64_vec(p_lm.data)).to_f64_vec();



@@ 116,11 146,14 @@ fn main() {
        .set_xlabel("$x$")
        .set_ylabel("$y$")
        .set_marker(vec![Point, Line])
        .savefig().expect("Can't draw plot");
        .savefig()
        .expect("Can't draw plot");
}

// Non autonomous case
fn f(p: Vec<Number>) -> Vec<Number> {
    (0 .. SIZE).map(|t| Number::from_f64(t as f64)).map(|t| p[0] * (-t / p[1]).exp() + p[2] * t * (-t / p[3]).exp()).collect()
    (0..SIZE)
        .map(|t| Number::from_f64(t as f64))
        .map(|t| p[0] * (-t / p[1]).exp() + p[2] * t * (-t / p[3]).exp())
        .collect()
}


D src/bin/pickle.rs => src/bin/pickle.rs +0 -26
@@ 1,26 0,0 @@
extern crate peroxide;
use peroxide::*;

fn main() {
    let a = seq(1, 100, 1);
    //a.write_single_pickle("example_data/pickle_example_vec.pickle").expect("Can't write pickle");

    let mut b = Matrix::from_index(|i, j| (i + j) as f64, (100, 2));
    match b.shape {
        Row => {
            b = b.change_shape(); // To column shape
        }
        Col => (),
    }

    let c = MATLAB::new("1 2; 3 4");

    //b.write_single_pickle("example_data/pickle_example_mat.pickle").expect("Can't write pickle");
    let mut w = SimpleWriter::new();
    w.insert_header(vec!["x", "y"])
        .insert_matrix(b)
        .insert_matrix(c)
        .insert_vector(a)
        .set_path("example_data/pickle_example.pickle")
        .write_pickle();
}

M src/bin/plot.rs => src/bin/plot.rs +1 -1
@@ 41,4 41,4 @@ fn test_fn(st: &mut State<f64>) {
fn stop(st: &ExplicitODE) -> bool {
    let y = &st.get_state().value[0];
    (*y - 2.4).abs() < 0.0001
}
\ No newline at end of file
}

M src/bin/read_csv.rs => src/bin/read_csv.rs +1 -1
@@ 2,7 2,7 @@ extern crate peroxide;
use peroxide::*;

fn main() {
    let a = Matrix::read("test.csv", false, ',');
    let a = Matrix::read("example_data/test.csv", false, ',');
    match a {
        Ok(m) => m.print(),
        Err(err) => println!("{}", err),

M src/bin/thomas.rs => src/bin/thomas.rs +2 -2
@@ 14,8 14,8 @@ fn main() {

    let a = tdma(prev, center, post, y);
    a.print();
    a.write_single_pickle("example_data/tdma.pickle")
        .expect("Can't write pickle file");
    // a.write_single_pickle("example_data/tdma.pickle")
    //     .expect("Can't write pickle file");
}

/// RHS with Dirichlet Boundary Condition

M src/lib.rs => src/lib.rs +5 -5
@@ 66,8 66,8 @@
//!     }
//!     ```

extern crate rand;
extern crate pyo3;
extern crate rand;

pub mod statistics;
pub mod structure;


@@ 80,6 80,9 @@ pub mod special;
pub mod util;

#[allow(unused_imports)]
pub use macros::{julia_macro::*, matlab_macro::*, r_macro::*};

#[allow(unused_imports)]
pub use structure::matrix::*;

#[allow(unused_imports)]


@@ 155,9 158,6 @@ pub use special::function::*;
pub use statistics::ops::*;

#[allow(unused_imports)]
pub use util::pickle::*;

#[allow(unused_imports)]
pub use structure::hyper_dual::*;

#[allow(unused_imports)]


@@ 176,4 176,4 @@ pub use operation::number::*;
pub use util::plot::*;

#[allow(unused_imports)]
pub use numerical::optimize::*;
\ No newline at end of file
pub use numerical::optimize::*;

M src/macros/julia_macro.rs => src/macros/julia_macro.rs +1 -1
@@ 60,4 60,4 @@ macro_rules! vstack {
            matrix(temp0, r, c, Row)
        }
    };
}
\ No newline at end of file
}

M src/macros/mod.rs => src/macros/mod.rs +1 -1
@@ 1,5 1,5 @@
//! Useful macros

pub mod julia_macro;
pub mod matlab_macro;
pub mod r_macro;
pub mod julia_macro;
\ No newline at end of file

M src/macros/r_macro.rs => src/macros/r_macro.rs +1 -0
@@ 273,6 273,7 @@ macro_rules! runif {
/// ```
/// extern crate peroxide;
/// use peroxide::*;
/// use peroxide::*;
///
/// let a = matrix!(1;5;1, 5, 1, Col);
/// let b = matrix(c!(3.7, 4.2, 4.9, 5.7, 6.0), 5, 1, Col);

M src/numerical/mod.rs => src/numerical/mod.rs +1 -1
@@ 5,6 5,6 @@
pub mod interp;
//pub mod newton;
pub mod ode;
pub mod optimize;
pub mod spline;
pub mod utils;
pub mod optimize;
\ No newline at end of file

M src/numerical/ode.rs => src/numerical/ode.rs +15 -26
@@ 152,16 152,7 @@
//!     let results = ex_test.integrate();
//!     let results2 = ex_test2.integrate();
//!
//!     // =========================================
//!     //  Write results to pickle
//!     // =========================================
//!     let mut wt = SimpleWriter::new();
//!
//!     wt
//!         .set_path("example_data/lorenz.pickle")
//!         .insert_matrix(results)
//!         .insert_matrix(results2)
//!         .write_pickle();
//!     // Plot or extract
//! }
//!
//! fn f(st: &mut State<f64>) {


@@ 200,10 191,7 @@
//!
//!     let result = ode_solver.integrate();
//!
//!     let mut st = SimpleWriter::new();
//!     st.set_path("example_data/rk4_test.pickle")
//!         .insert_matrix(result)
//!         .write_pickle();
//!     // Plot or Extract..
//! }
//!
//! fn test_fn(st: &mut State<f64>) {


@@ 214,16 202,17 @@
//! }
//! ```

use self::BoundaryCondition::Dirichlet;
use self::ExMethod::{Euler, RK4};
use self::ODEOptions::{BoundCond, InitCond, Method, StepSize, StopCond, Times};
use operation::extra_ops::Real;
use operation::mut_ops::MutFP;
use std::collections::HashMap;
use BoundaryCondition::Dirichlet;
use ExMethod::{Euler, RK4};
use ODEOptions::{BoundCond, InitCond, Method, StepSize, StopCond, Times};
use {cat, zeros};
use {Dual, Real};
use {FPVector, Matrix, MutFP};
use structure::dual::Dual;
use structure::matrix::{Matrix, Row, FP};
use structure::vector::FPVector;
use util::non_macro::{cat, zeros};
use util::print::Printable;
use FP;
use ::Shape::Row;

/// Explicit ODE Methods
///


@@ 449,14 438,14 @@ impl ODE for ExplicitODE {

        let mut result = zeros(self.times + 1, self.state.value.len() + 1);

        result.subs_row(0, cat(self.state.param, self.state.value.clone()));
        result.subs_row(0, &cat(self.state.param, self.state.value.clone()));

        match self.options.get(&StopCond) {
            Some(stop) if *stop => {
                let mut key = 1usize;
                for i in 1..self.times + 1 {
                    self.mut_update();
                    result.subs_row(i, cat(self.state.param, self.state.value.clone()));
                    result.subs_row(i, &cat(self.state.param, self.state.value.clone()));
                    key += 1;
                    if (self.stop_cond)(&self) {
                        println!("Reach the stop condition!");


@@ 466,11 455,11 @@ impl ODE for ExplicitODE {
                    }
                }
                return result.take(key, Row);
            },
            }
            _ => {
                for i in 1..self.times + 1 {
                    self.mut_update();
                    result.subs_row(i, cat(self.state.param, self.state.value.clone()));
                    result.subs_row(i, &cat(self.state.param, self.state.value.clone()));
                }
                return result;
            }

M src/numerical/optimize.rs => src/numerical/optimize.rs +45 -41
@@ 1,14 1,14 @@
//! To optimize parametric model (non-linear regression)
//! 
//!
//! ## `Optimizer` structure
//! 
//!
//! ### Declaration
//! 
//!
//! ```rust
//! extern crate peroxide;
//! use peroxide::{Number, OptMethod, OptOption};
//! use std::collections::HashMap;
//! 
//!
//! pub struct Optimizer {
//!     domain: Vec<f64>,
//!     observed: Vec<f64>,


@@ 20,48 20,48 @@
//!     option: HashMap<OptOption, bool>,
//! }
//! ```
//! 
//!
//! ### Method (Should fill)
//! 
//!
//! * `new` : Declare new Optimizer. **Should be mutable**
//! * `set_init_param` : Input initial parameter
//! * `set_max_iter` : Set maximum number of iterations
//! * `set_method` : Set method to optimization
//! 
//!
//! ### Method (Optional)
//! 
//!
//! * `get_domain` : Get domain
//! * `get_error` : Root mean square error
//! 
//!
//! ### Method (Generate result)
//! 
//!
//! * `optimize` : Optimize
//! 
//!
//! ## Example
//! 
//!
//! * Optimize $y = x^n$ model with $y = x^2$ with gaussian noise.
//! 
//!
//! ```rust
//! extern crate peroxide;
//! use peroxide::*;
//! 
//!
//! fn main() {
//!     // To prepare noise
//!     let normal = Normal(0f64, 0.1f64);
//!     let normal2 = Normal(0f64, 100f64);
//! 
//!
//!     // Noise to domain
//!     let mut x = seq(0., 99., 1f64);
//!     x = zip_with(|a, b| (a + b).abs(), &x, &normal.sample(x.len()));
//! 
//!
//!     // Noise to image
//!     let mut y = x.fmap(|t| t.powi(2));
//!     y = zip_with(|a, b| a + b, &y, &normal2.sample(y.len()));
//! 
//!
//!     // Initial paramter
//!     let n_init = vec![1f64];
//!     let data = hstack!(x.clone(), y.clone());
//! 
//!
//!     // Optimizer setting
//!     let mut opt = Optimizer::new(data, quad);
//!     let p = opt.set_init_param(n_init)


@@ 70,10 70,10 @@
//!         .optimize();
//!     p.print();                  // Optimized parameter
//!     opt.get_error().print();    // Optimized RMSE
//! 
//!
//!     // To prepare plotting
//!     let z = quad(&x, NumberVector::from_f64_vec(p)).to_f64_vec();
//! 
//!
//!     // Plot
//!     let mut plt = Plot2D::new();
//!     plt.set_domain(x)


@@ 85,7 85,7 @@
//!         .set_marker(vec![Point, Line])
//!         .savefig().expect("Can't draw a plot");
//! }
//! 
//!
//! fn quad(x: &Vec<f64>, n: Vec<Number>) -> Vec<Number> {
//!     x.clone().into_iter()
//!         .map(|t| Number::from_f64(t))


@@ 93,15 93,16 @@
//!         .collect()
//! }
//! ```
//! 
//!
//! ![LM test](https://raw.githubusercontent.com/Axect/Peroxide/master/example_data/lm_test.png)


pub use self::OptMethod::{GaussNewton, GradientDescent, LevenbergMarquardt};
use self::OptOption::{InitParam, MaxIter};
use numerical::utils::jacobian;
use operation::number::{Number, NumberVector};
use std::collections::HashMap;
use ::{Matrix, Number, NumberVector, LinearAlgebra, LinearOps};
use ::{max, jacobian};
use ::OptOption::{InitParam, MaxIter};
pub use ::OptMethod::{GradientDescent, GaussNewton, LevenbergMarquardt};
use structure::matrix::{LinearAlgebra, LinearOps, Matrix};
use util::useful::max;

#[derive(Debug, Clone, Copy)]
pub enum OptMethod {


@@ 117,7 118,7 @@ pub enum OptOption {
}

/// Optimizer for optimization (non-linear regression)
/// 
///
/// # Methods
/// * Gradient Descent
/// * Gauss Newton


@@ 200,25 201,25 @@ impl Optimizer {

        match self.method {
            GradientDescent(alpha) => {
                for _i in 0 .. max_iter {
                for _i in 0..max_iter {
                    let h = alpha * j.t() * (&y - &y_hat);
                    p = &p + &h;
                    j = jacobian(f, p.data.clone());
                    y_hat = f(NumberVector::from_f64_vec(p.data.clone())).to_f64_vec().to_matrix();
                    y_hat = f(NumberVector::from_f64_vec(p.data.clone()))
                        .to_f64_vec()
                        .to_matrix();
                }
            }

            GaussNewton => {
                unimplemented!()
            }
            GaussNewton => unimplemented!(),

            LevenbergMarquardt => {
                let mut chi2 = ((&y - &y_hat).t() * (&y - &y_hat))[(0,0)];
                let mut chi2 = ((&y - &y_hat).t() * (&y - &y_hat))[(0, 0)];
                let mut nu = 2f64;
                let lambda_0 = 1e-3;
                let mut lambda = lambda_0 * max(jtj.diag());

                for _i in 0 .. max_iter {
                for _i in 0..max_iter {
                    let h: Matrix;

                    match (jtj.clone() + lambda * jtj.to_diag()).inv() {


@@ 228,10 229,14 @@ impl Optimizer {

                    let p_temp = &p + &h;
                    let j_temp = jacobian(f, p.data.clone());
                    let y_hat_temp = f(NumberVector::from_f64_vec(p_temp.data.clone())).to_f64_vec().to_matrix();
                    let chi2_temp = ((&y - &y_hat_temp).t() * (&y - &y_hat_temp))[(0,0)];
                    let y_hat_temp = f(NumberVector::from_f64_vec(p_temp.data.clone()))
                        .to_f64_vec()
                        .to_matrix();
                    let chi2_temp = ((&y - &y_hat_temp).t() * (&y - &y_hat_temp))[(0, 0)];

                    let rho = (chi2 - chi2_temp) / (h.t() * (lambda * jtj.to_diag() * h.clone() + j.t() * (&y - &y_hat)))[(0,0)];
                    let rho = (chi2 - chi2_temp)
                        / (h.t() * (lambda * jtj.to_diag() * h.clone() + j.t() * (&y - &y_hat)))
                            [(0, 0)];

                    if rho > 0f64 {
                        p = p_temp;


@@ 239,7 244,7 @@ impl Optimizer {
                        jtj = &j.t() * &j;
                        y_hat = y_hat_temp;
                        chi2 = chi2_temp;
                        lambda *= max(vec![1f64/3f64, 1f64 - (2f64*rho - 1f64).powi(3)]);
                        lambda *= max(vec![1f64 / 3f64, 1f64 - (2f64 * rho - 1f64).powi(3)]);
                        nu = 2f64;
                    } else {
                        lambda *= nu;


@@ 249,8 254,7 @@ impl Optimizer {
            }
        }
        let error_temp = &y - &y_hat;
        self.error = ((error_temp.t() * error_temp)[(0,0)] / y.row as f64).sqrt();
        self.error = ((error_temp.t() * error_temp)[(0, 0)] / y.row as f64).sqrt();
        p.data
    }
}


M src/numerical/utils.rs => src/numerical/utils.rs +1 -1
@@ 1,7 1,7 @@
use operation::number::{Number, NumberVector};
use structure::dual::*;
use structure::matrix::*;
use util::non_macro::{cat, zeros};
use ::{Number, NumberVector};

/// Jacobian Matrix
///

M src/operation/extra_ops.rs => src/operation/extra_ops.rs +2 -1
@@ 25,7 25,8 @@
//!     ```

use std::ops::{Add, Div, Mul, Sub};
use {Dual, HyperDual};
use structure::dual::Dual;
use structure::hyper_dual::HyperDual;

pub trait PowOps: Sized {
    fn powi(&self, n: i32) -> Self;

M src/operation/number.rs => src/operation/number.rs +9 -6
@@ 1,9 1,9 @@
use operation::extra_ops::{ExpLogOps, PowOps, Real, TrigOps};
use operation::number::Number::{D, E, F};
use std::ops::{Add, Div, Mul, Neg, Sub};
use std::process::exit;
use {Dual, Real};
use {ExpLogOps, TrigOps};
use {HyperDual, PowOps};
use structure::dual::Dual;
use structure::hyper_dual::HyperDual;

#[derive(Debug, Copy, Clone, PartialEq)]
pub enum NumError {


@@ 325,7 325,7 @@ impl PowOps for Number {
            (D(x), D(y)) => D(x.powf(*y)),
            (E(x), _) => E(x.to_owned()),
            (_, E(y)) => E(y.to_owned()),
            _ => unreachable!()
            _ => unreachable!(),
        }
    }



@@ 433,7 433,10 @@ impl NumberVector for Vec<Number> {
    }

    fn to_hyper_vec(&self) -> Vec<HyperDual> {
        self.clone().into_iter().map(|x| x.to_hyper_dual()).collect()
        self.clone()
            .into_iter()
            .map(|x| x.to_hyper_dual())
            .collect()
    }

    fn from_dual_vec(v: Vec<Dual>) -> Self {


@@ 447,4 450,4 @@ impl NumberVector for Vec<Number> {
    fn from_hyper_vec(v: Vec<HyperDual>) -> Self {
        v.into_iter().map(|x| Number::from_hyper_dual(x)).collect()
    }
}
\ No newline at end of file
}

M src/statistics/stat.rs => src/statistics/stat.rs +2 -2
@@ 289,8 289,8 @@ pub fn cov(v1: &Vector, v2: &Vector) -> f64 {

    for (x, y) in v1.into_iter().zip(v2) {
        ss += x * y;
        sx += x;
        sy += y;
        sx += *x;
        sy += *y;
        l += 1f64;
    }
    assert_ne!(l, 1f64);

M src/structure/hyper_dual.rs => src/structure/hyper_dual.rs +0 -2
@@ 415,8 415,6 @@ impl PowOps for HyperDual {
        } else {
            unimplemented!()
        }


    }

    fn sqrt(&self) -> Self {

M src/structure/matrix.rs => src/structure/matrix.rs +143 -221
@@ 125,8 125,8 @@
//!         // r[1]    1    2
//!
//!         let mut b = ml_matrix("1 2;3 4");
//!         b.subs_col(0, c!(5, 6));
//!         b.subs_row(1, c!(7, 8));
//!         b.subs_col(0, &c!(5, 6));
//!         b.subs_row(1, &c!(7, 8));
//!         b.print();
//!         //       c[0] c[1]
//!         // r[0]    5    2


@@ 136,7 136,7 @@
//!
//! ## Read & Write
//!
//! In peroxide, we can write matrix to `csv` or `pickle`.
//! In peroxide, we can write matrix to `csv`
//!
//! ### CSV (Legacy)
//!


@@ 149,10 149,10 @@
//!
//!     fn main() {
//!         let a = ml_matrix("1 2;3 4");
//!         a.write("matrix.csv").expect("Can't write file");
//!         a.write("example_data/matrix.csv").expect("Can't write file");
//!
//!         let b = ml_matrix("1 2; 3 4; 5 6");
//!         b.write_with_header("header.csv", vec!["odd", "even"])
//!         b.write_with_header("example_data/header.csv", vec!["odd", "even"])
//!             .expect("Can't write header file");
//!     }
//!     ```


@@ 176,34 176,6 @@
//!     }
//!     ```
//!
//! ### Pickle (Export as python object)
//!
//! * `SimpleWriter` : Struct to write pickle
//!     * Necessary method
//!         * `set_path` : Set path
//!         * `insert_matrix` or `insert_vector`
//!         * `write_pickle` : Must be at last
//!     * Optional method
//!         * `set_round_level` : Set round-off level
//!         * `insert_header` : Insert header
//!
//!     ```rust
//!     extern crate peroxide;
//!     use peroxide::*;
//!
//!     fn main() {
//!         let a = ml_matrix("1 2;3 4");
//!         let v = c!(1,2,3,4);
//!
//!         let mut wrt = SimpleWriter::new();
//!         wrt.set_path("example_data/sample.pickle")
//!             .insert_matrix(a)
//!             .insert_vector(v)
//!             .set_round_level(4)
//!             .write_pickle();
//!     }
//!     ```
//!
//! ## Concatenation
//!
//! There are two options to concatenate matrices.


@@ 367,7 339,7 @@
//!
//!     ```rust
//!     extern crate peroxide;
//!     use peroxide::{Matrix};
//!     use peroxide::*;
//!
//!     #[derive(Debug, Clone)]
//!     pub struct PQLU {


@@ 453,7 425,6 @@
//!     }
//!     ```


extern crate csv;

use self::csv::{ReaderBuilder, StringRecord, WriterBuilder};


@@ 661,7 632,7 @@ impl Matrix {
        assert_eq!(r * c, self.data.len());
        let l = r * c - 1;
        let mut data: Vec<f64> = self.data.clone();
        let ref_data: Vec<f64> = self.data.clone();
        let ref_data = &self.data;

        match self.shape {
            Row => {


@@ 793,26 764,9 @@ impl Matrix {
    /// ```
    pub fn col(&self, index: usize) -> Vector {
        assert!(index < self.col);
        let mut container: Vec<f64> = Vec::new();
        match self.shape {
            Row => {
                let l: usize = self.row * self.col;
                for i in 0..l {
                    if i % self.col == index {
                        container.push(self.data[i]);
                    }
                }
            }
            Col => {
                let s: usize = self.row * index;
                container = self
                    .data
                    .clone()
                    .into_iter()
                    .skip(s)
                    .take(self.row)
                    .collect::<Vec<f64>>();
            }
        let mut container: Vec<f64> = vec![0f64; self.row];
        for i in 0..self.row {
            container[i] = self[(i, index)];
        }
        container
    }


@@ 829,26 783,9 @@ impl Matrix {
    /// ```
    pub fn row(&self, index: usize) -> Vector {
        assert!(index < self.row);
        let mut container: Vec<f64> = Vec::new();
        match self.shape {
            Row => {
                let s: usize = self.col * index;
                container = self
                    .data
                    .clone()
                    .into_iter()
                    .skip(s)
                    .take(self.col)
                    .collect::<Vec<f64>>();
            }
            Col => {
                let l: usize = self.row * self.col;
                for i in 0..l {
                    if i % self.row == index {
                        container.push(self.data[i]);
                    }
                }
            }
        let mut container: Vec<f64> = vec![0f64; self.col];
        for i in 0..self.col {
            container[i] = self[(index, i)];
        }
        container
    }


@@ 864,12 801,13 @@ impl Matrix {
    /// assert_eq!(a.diag(), c!(1,4));
    /// ```
    pub fn diag(&self) -> Vector {
        let mut container: Vector = Vec::new();
        let mut container = vec![0f64; self.row];
        let r = self.row;
        let c = self.col;
        assert_eq!(r, c);
        let c2 = c + 1;
        for i in 0..r {
            container.push(self.data[i * (c + 1)]);
            container[i] = self.data[i * c2];
        }
        container
    }


@@ 924,7 862,7 @@ impl Matrix {
    /// use peroxide::*;
    ///
    /// let a = matrix(c!(1,2,3,3,2,1), 3, 2, Col);
    /// a.write("test.csv");
    /// a.write("example_data/test.csv");
    /// ```
    pub fn write(&self, file_path: &str) -> Result<(), Box<dyn Error>> {
        let mut wtr = WriterBuilder::new().from_path(file_path)?;


@@ 949,7 887,7 @@ impl Matrix {
    /// use peroxide::*;
    ///
    /// let a = matrix(c!(1,2,3,3,2,1), 3, 2, Col);
    /// a.write_round("test.csv", 0);
    /// a.write_round("example_data/test.csv", 0);
    /// ```
    pub fn write_round(&self, file_path: &str, round: usize) -> Result<(), Box<dyn Error>> {
        let mut wtr = WriterBuilder::new().from_path(file_path)?;


@@ 1018,9 956,9 @@ impl Matrix {
    /// use std::process;
    ///
    /// let a = matrix(c!(1,2,3,3,2,1), 3, 2, Col);
    /// a.write_round("test.csv", 0);
    /// a.write_round("example_data/test.csv", 0);
    ///
    /// let b = Matrix::read("test.csv", false, ','); // header = false, delimiter = ','
    /// let b = Matrix::read("example_data/test.csv", false, ','); // header = false, delimiter = ','
    /// match b {
    ///     Ok(mat) => println!("{}", mat),
    ///     Err(err) => {


@@ 1055,14 993,14 @@ impl Matrix {
    }

    /// Substitute Col
    pub fn subs_col(&mut self, idx: usize, v: Vec<f64>) {
    pub fn subs_col(&mut self, idx: usize, v: &Vec<f64>) {
        for i in 0..self.row {
            self[(i, idx)] = v[i];
        }
    }

    /// Substitute Row
    pub fn subs_row(&mut self, idx: usize, v: Vec<f64>) {
    pub fn subs_row(&mut self, idx: usize, v: &Vec<f64>) {
        for j in 0..self.col {
            self[(idx, j)] = v[j];
        }


@@ 1092,7 1030,7 @@ impl Matrix {
    /// To send `Matrix` to `inline-python`
    pub fn to_vec(&self) -> Vec<Vec<f64>> {
        let mut result = vec![vec![0f64; self.col]; self.row];
        for i in 0 .. self.row {
        for i in 0..self.row {
            result[i] = self.row(i);
        }
        result


@@ 1102,8 1040,8 @@ impl Matrix {
        assert!(self.row == self.col, "Should be square matrix");
        let mut result = matrix(vec![0f64; self.row * self.col], self.row, self.col, Row);
        let diag = self.diag();
        for i in 0 .. self.row {
            result[(i,i)] = diag[i];
        for i in 0..self.row {
            result[(i, i)] = diag[i];
        }
        result
    }


@@ 1192,20 1130,13 @@ impl Add<Matrix> for Matrix {
    fn add(self, other: Self) -> Self {
        assert_eq!(&self.row, &other.row);
        assert_eq!(&self.col, &other.col);
        if self.shape == other.shape {
            matrix(
                self.data
                    .into_iter()
                    .zip(&other.data)
                    .map(|(x, y)| x + y)
                    .collect::<Vec<f64>>(),
                self.row,
                self.col,
                self.shape,
            )
        } else {
            self.change_shape().add(other)
        let mut result = matrix(self.data.clone(), self.row, self.col, self.shape);
        for i in 0..self.row {
            for j in 0..self.col {
                result[(i, j)] += other[(i, j)];
            }
        }
        result
    }
}



@@ 1216,22 1147,13 @@ impl<'a, 'b> Add<&'b Matrix> for &'a Matrix {
    fn add(self, other: &'b Matrix) -> Self::Output {
        assert_eq!(self.row, other.row);
        assert_eq!(self.col, other.col);
        if self.shape == other.shape {
            Matrix {
                data: self
                    .data
                    .clone()
                    .into_iter()
                    .zip(&other.data)
                    .map(|(x, y)| x + y)
                    .collect::<Vec<f64>>(),
                row: self.row,
                col: self.col,
                shape: self.shape,
        let mut result = matrix(self.data.clone(), self.row, self.col, self.shape);
        for i in 0..self.row {
            for j in 0..self.col {
                result[(i, j)] += other[(i, j)];
            }
        } else {
            (&self.change_shape()).add(other)
        }
        result
    }
}



@@ 1412,20 1334,13 @@ impl Sub<Matrix> for Matrix {
    fn sub(self, other: Self) -> Self {
        assert_eq!(&self.row, &other.row);
        assert_eq!(&self.col, &other.col);
        if self.shape == other.shape {
            matrix(
                self.data
                    .into_iter()
                    .zip(&other.data)
                    .map(|(x, y)| x - y)
                    .collect::<Vec<f64>>(),
                self.row,
                self.col,
                self.shape,
            )
        } else {
            self.change_shape().sub(other)
        let mut result = self.clone();
        for i in 0 .. self.row {
            for j in 0 .. self.col {
                result[(i,j)] -= other[(i,j)];
            }
        }
        result
    }
}



@@ 1436,22 1351,13 @@ impl<'a, 'b> Sub<&'b Matrix> for &'a Matrix {
    fn sub(self, other: &'b Matrix) -> Self::Output {
        assert_eq!(self.row, other.row);
        assert_eq!(self.col, other.col);
        if self.shape == other.shape {
            Matrix {
                data: self
                    .data
                    .clone()
                    .into_iter()
                    .zip(&other.data)
                    .map(|(x, y)| x - y)
                    .collect::<Vec<f64>>(),
                row: self.row,
                col: self.col,
                shape: self.shape,
        let mut result = self.clone();
        for i in 0 .. self.row {
            for j in 0 .. self.col {
                result[(i,j)] -= other[(i,j)];
            }
        } else {
            (&self.change_shape()).sub(other)
        }
        result
    }
}



@@ 1691,29 1597,7 @@ where
    type Output = Self;

    fn mul(self, other: T) -> Self {
        let r_self = self.row;
        let c_self = self.col;
        let new_other = other.to_matrix();
        let r_other = new_other.row;
        let c_other = new_other.col;

        assert_eq!(c_self, r_other);

        let r_new = r_self;
        let c_new = c_other;

        let mut result = matrix(vec![0f64; r_new * c_new], r_new, c_new, self.shape);

        for i in 0..r_new {
            for j in 0..c_new {
                let mut s = 0f64;
                for k in 0..c_self {
                    s += self[(i, k)] * new_other[(k, j)];
                }
                result[(i, j)] = s;
            }
        }
        result
        matmul(&self, &other.to_matrix())
    }
}



@@ 1724,27 1608,7 @@ where
    type Output = Matrix;

    fn mul(self, other: &'b T) -> Self::Output {
        let new_other = other.to_matrix();

        assert_eq!(self.col, new_other.row);

        let mut result = matrix(
            vec![0f64; self.row * new_other.col],
            self.row,
            new_other.col,
            self.shape,
        );

        for i in 0..self.row {
            for j in 0..new_other.col {
                let mut s = 0f64;
                for k in 0..self.col {
                    s += self[(i, k)] * new_other[(k, j)];
                }
                result[(i, j)] = s;
            }
        }
        result
        matmul(self, &other.to_matrix())
    }
}



@@ 1875,9 1739,14 @@ impl FP for Matrix {
                }
                match self.shape {
                    Row => {
                        let new_data = self.data.clone().into_iter().take(n * self.col).collect::<Vec<f64>>();
                        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 {


@@ 1893,9 1762,14 @@ impl FP for Matrix {
                }
                match self.shape {
                    Col => {
                        let new_data = self.data.clone().into_iter().take(n * self.row).collect::<Vec<f64>>();
                        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 {


@@ 1952,11 1826,11 @@ impl FP for Matrix {
            data: vec![0f64; self.row * self.col],
            row: self.row,
            col: self.col,
            shape: Col
            shape: Col,
        };

        for i in 0 .. self.row {
            result.subs_col(i, f(self.col(i)));
        for i in 0..self.row {
            result.subs_col(i, &f(self.col(i)));
        }

        result


@@ 1964,17 1838,17 @@ impl FP for Matrix {

    fn row_map<F>(&self, f: F) -> Matrix
    where
        F: Fn(Vec<f64>) -> Vec<f64>
        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
            shape: Row,
        };

        for i in 0 .. self.col {
            result.subs_row(i, f(self.row(i)));
        for i in 0..self.col {
            result.subs_row(i, &f(self.row(i)));
        }

        result


@@ 2068,7 1942,12 @@ pub struct PQLU {

impl PQLU {
    pub fn extract(&self) -> (Perms, Perms, Matrix, Matrix) {
        (self.p.clone(), self.q.clone(), self.l.clone(), self.u.clone())
        (
            self.p.clone(),
            self.q.clone(),
            self.l.clone(),
            self.u.clone(),
        )
    }
}



@@ 2289,29 2168,30 @@ impl LinearAlgebra for Matrix {
    fn block(&self) -> (Self, Self, Self, Self) {
        let r = self.row;
        let c = self.col;
        let l = min(self.row / 2, self.col / 2);
        let r_l = r - l;
        let c_l = c - l;

        let mut m1 = matrix(vec![0f64; l * l], l, l, self.shape);
        let mut m2 = matrix(vec![0f64; l * c_l], l, c_l, self.shape);
        let mut m3 = matrix(vec![0f64; r_l * l], r_l, l, self.shape);
        let l_r = self.row / 2;
        let l_c = self.col / 2;
        let r_l = r - l_r;
        let c_l = c - l_c;

        let mut m1 = matrix(vec![0f64; l_r * l_c], l_r, l_c, self.shape);
        let mut m2 = matrix(vec![0f64; l_r * c_l], l_r, c_l, self.shape);
        let mut m3 = matrix(vec![0f64; r_l * l_c], r_l, l_c, self.shape);
        let mut m4 = matrix(vec![0f64; r_l * c_l], r_l, c_l, self.shape);

        for idx_row in 0..r {
            for idx_col in 0..c {
                match (idx_row, idx_col) {
                    (i, j) if (i < l) && (j < l) => {
                    (i, j) if (i < l_r) && (j < l_c) => {
                        m1[(i, j)] = self[(i, j)];
                    }
                    (i, j) if (i < l) && (j >= l) => {
                        m2[(i, j - l)] = self[(i, j)];
                    (i, j) if (i < l_r) && (j >= l_c) => {
                        m2[(i, j - l_c)] = self[(i, j)];
                    }
                    (i, j) if (i >= l) && (j < l) => {
                        m3[(i - l, j)] = self[(i, j)];
                    (i, j) if (i >= l_r) && (j < l_c) => {
                        m3[(i - l_r, j)] = self[(i, j)];
                    }
                    (i, j) if (i >= l) && (j >= l) => {
                        m4[(i - l, j - l)] = self[(i, j)];
                    (i, j) if (i >= l_r) && (j >= l_c) => {
                        m4[(i - l_r, j - l_c)] = self[(i, j)];
                    }
                    _ => (),
                }


@@ 2409,29 2289,30 @@ impl LinearAlgebra for Matrix {
/// assert_eq!(n, b);
/// ```
pub fn combine(m1: Matrix, m2: Matrix, m3: Matrix, m4: Matrix) -> Matrix {
    let l = m1.col;
    let l_r = m1.row;
    let l_c = m1.col;
    let c_l = m2.col;
    let r_l = m3.row;

    let r = l + r_l;
    let c = l + c_l;
    let r = l_r + r_l;
    let c = l_c + c_l;

    let mut m = matrix(vec![0f64; r * c], r, c, m1.shape);

    for idx_row in 0..r {
        for idx_col in 0..c {
            match (idx_row, idx_col) {
                (i, j) if (i < l) && (j < l) => {
                (i, j) if (i < l_r) && (j < l_c) => {
                    m[(i, j)] = m1[(i, j)];
                }
                (i, j) if (i < l) && (j >= l) => {
                    m[(i, j)] = m2[(i, j - l)];
                (i, j) if (i < l_r) && (j >= l_c) => {
                    m[(i, j)] = m2[(i, j - l_c)];
                }
                (i, j) if (i >= l) && (j < l) => {
                    m[(i, j)] = m3[(i - l, j)];
                (i, j) if (i >= l_r) && (j < l_c) => {
                    m[(i, j)] = m3[(i - l_r, j)];
                }
                (i, j) if (i >= l) && (j >= l) => {
                    m[(i, j)] = m4[(i - l, j - l)];
                (i, j) if (i >= l_r) && (j >= l_c) => {
                    m[(i, j)] = m4[(i - l_r, j - l_c)];
                }
                _ => (),
            }


@@ 2515,3 2396,44 @@ pub fn inv_u(u: Matrix) -> Matrix {
        }
    }
}

fn matmul(a: &Matrix, b: &Matrix) -> Matrix {
    match (a.row, a.col) {
        (p, q) if p <= 100 && q <= 100 => {
            let r_self = a.row;
            let c_self = a.col;
            let new_other = b;
            let r_other = new_other.row;
            let c_other = new_other.col;

            assert_eq!(c_self, r_other);

            let r_new = r_self;
            let c_new = c_other;

            let mut result = matrix(vec![0f64; r_new * c_new], r_new, c_new, a.shape);

            for i in 0..r_new {
                for j in 0..c_new {
                    let mut s = 0f64;
                    for k in 0..c_self {
                        s += a[(i, k)] * new_other[(k, j)];
                    }
                    result[(i, j)] = s;
                }
            }
            result
        },
        _ => {
            let (a1, a2, a3, a4) = a.block();
            let (b1, b2, b3, b4) = b.block();

            let m1 = matmul(&a1,&b1) + matmul(&a2,&b3);
            let m2 = matmul(&a1,&b2) + matmul(&a2,&b4);
            let m3 = matmul(&a3,&b1) + matmul(&a4,&b3);
            let m4 = matmul(&a3,&b2) + matmul(&a4,&b4);

            combine(m1, m2, m3, m4)
        }
    }
}
\ No newline at end of file

M src/structure/vector.rs => src/structure/vector.rs +16 -10
@@ 262,10 262,10 @@
//!     }
//!     ```

use operation::extra_ops::Real;
use std::cmp::min;
use std::convert;
use std::f64::MIN;
use Real;
use std::cmp::min;

pub type Vector = Vec<f64>;



@@ 400,37 400,43 @@ impl FPVector for Vector {
}

pub fn map<F, T>(f: F, xs: &Vec<T>) -> Vec<T>
    where F: Fn(T) -> T, T: Real + Default
where
    F: Fn(T) -> T,
    T: Real + Default,
{
    let l = xs.len();
    let mut result = vec![T::default(); l];
    for i in 0 .. l {
    for i in 0..l {
        result[i] = f(xs[i]);
    }
    result
}

pub fn reduce<F, T>(f: F, init: T, xs: &Vec<T>) -> T
where F: Fn(T, T) -> T, T: Real {
where
    F: Fn(T, T) -> T,
    T: Real,
{
    let mut s = init;
    for i in 0 .. xs.len() {
    for i in 0..xs.len() {
        s = f(s, xs[i]);
    }
    s
}

pub fn zip_with<F, T>(f: F, xs: &Vec<T>, ys: &Vec<T>) -> Vec<T>
    where F: Fn(T, T) -> T, T: Real + Default
where
    F: Fn(T, T) -> T,
    T: Real + Default,
{
    let l = min(xs.len(), ys.len());
    let mut result = vec![T::default(); l];
    for i in 0 .. l {
    for i in 0..l {
        result[i] = f(xs[i], ys[i]);
    }
    result
}


/// Some algorithms for Vector
pub trait Algorithm {
    fn rank(&self) -> Vec<usize>;


@@ 554,7 560,7 @@ impl VecOps for Vector {

    /// Dot product
    fn dot(&self, other: &Self) -> f64 {
        zip_with(|x,y| x * y, &self, other).reduce(0, |x, y| x + y)
        zip_with(|x, y| x * y, &self, other).reduce(0, |x, y| x + y)
    }

    /// Norm

M src/util/mod.rs => src/util/mod.rs +0 -1
@@ 2,7 2,6 @@

pub mod api;
pub mod non_macro;
pub mod pickle;
pub mod plot;
pub mod print;
pub mod useful;

D src/util/pickle.rs => src/util/pickle.rs +0 -90
@@ 1,90 0,0 @@
//! To deal with python pickle data structure

extern crate serde;
extern crate serde_pickle;

use std::fs::File;
use std::io::Write;
use std::process::exit;
use structure::matrix::*;
use structure::vector::*;

/// Pickle trait
///
/// # Description
///
/// Use python pickle to export vector or matrix
pub trait Pickle {
    fn write_single_pickle(&self, path: &str) -> serde_pickle::Result<()>;
    fn write_pickle<W: Write>(&self, writer: &mut W) -> serde_pickle::Result<()>;
}

impl Pickle for Vector {
    fn write_single_pickle(&self, path: &str) -> serde_pickle::Result<()> {
        let mut writer: Box<dyn Write>;

        match File::create(path) {
            Ok(p) => writer = Box::new(p),
            Err(e) => {
                println!("{:?}", e);
                exit(1);
            }
        }

        serde_pickle::to_writer(&mut writer, &self, true)
    }

    fn write_pickle<W: Write>(&self, writer: &mut W) -> serde_pickle::Result<()> {
        serde_pickle::to_writer(writer, &self, true)
    }
}

impl Pickle for Matrix {
    fn write_single_pickle(&self, path: &str) -> serde_pickle::Result<()> {
        let mut writer: Box<dyn Write>;

        match File::create(path) {
            Ok(p) => writer = Box::new(p),
            Err(e) => {
                println!("{:?}", e);
                exit(1);
            }
        }

        let mut container: Vec<Vec<f64>> = Vec::new();;

        match self.shape {
            Row => {
                for i in 0..self.row {
                    container.push(self.row(i));
                }
            }
            Col => {
                for i in 0..self.col {
                    container.push(self.col(i));
                }
            }
        }

        serde_pickle::to_writer(&mut writer, &container, true)
    }

    fn write_pickle<W: Write>(&self, writer: &mut W) -> serde_pickle::Result<()> {
        let mut container: Vec<Vec<f64>> = Vec::new();;

        match self.shape {
            Row => {
                for i in 0..self.row {
                    container.push(self.row(i));
                }
            }
            Col => {
                for i in 0..self.col {
                    container.push(self.col(i));
                }
            }
        }

        serde_pickle::to_writer(writer, &container, true)
    }
}

M src/util/plot.rs => src/util/plot.rs +19 -16
@@ 21,11 21,7 @@
//!         .set_times(1000);
//!
//!     let result = ode_solver.integrate();
//!
//!     let mut st = SimpleWriter::new();
//!     st.set_path("example_data/rk4_test.pickle")
//!         .insert_matrix(result)
//!         .write_pickle();
//!     result.write("example_data/test.csv");
//! }
//!
//! fn test_fn(st: &mut State<f64>) {


@@ 82,14 78,13 @@
//!
//! ![test_plot](https://raw.githubusercontent.com/Axect/Peroxide/master/example_data/test_plot.png)


extern crate pyo3;
use self::pyo3::types::IntoPyDict;
use self::pyo3::{PyResult, Python};
pub use self::Grid::{Off, On};
pub use self::Markers::{Circle, Line, Point};
use self::PlotOptions::{Domain, Images, Legends, Path};
use std::collections::HashMap;
pub use Grid::*;
use PlotOptions::{Domain, Images, Legends, Path};
pub use Markers::{Point, Line, Circle};

type Vector = Vec<f64>;



@@ 211,7 206,10 @@ impl Plot for Plot2D {
        if let Some(x) = self.options.get_mut(&Legends) {
            *x = true
        }
        self.legends = legends.into_iter().map(|x| x.to_owned()).collect::<Vec<String>>();
        self.legends = legends
            .into_iter()
            .map(|x| x.to_owned())
            .collect::<Vec<String>>();
        self
    }



@@ 324,16 322,21 @@ impl Plot for Plot2D {

        if self.markers.len() == 0 {
            for i in 0..y_length {
                plot_string.push_str(&format!("plt.plot(x,y[{}],label=r\"{}\")\n", i, legends[i])[..])
                plot_string
                    .push_str(&format!("plt.plot(x,y[{}],label=r\"{}\")\n", i, legends[i])[..])
            }
        } else {
            for i in 0 .. y_length {
            for i in 0..y_length {
                match self.markers[i] {
                    Line => plot_string.push_str(&format!("plt.plot(x,y[{}],label=r\"{}\")\n", i, legends[i])[..]),
                    Point => plot_string.push_str(&format!("plt.plot(x,y[{}],\".\",label=r\"{}\")\n", i, legends[i])[..]),
                    Circle => plot_string.push_str(&format!("plt.plot(x,y[{}],\"o\",label=r\"{}\")\n", i, legends[i])[..]),
                    Line => plot_string
                        .push_str(&format!("plt.plot(x,y[{}],label=r\"{}\")\n", i, legends[i])[..]),
                    Point => plot_string.push_str(
                        &format!("plt.plot(x,y[{}],\".\",label=r\"{}\")\n", i, legends[i])[..],
                    ),
                    Circle => plot_string.push_str(
                        &format!("plt.plot(x,y[{}],\"o\",label=r\"{}\")\n", i, legends[i])[..],
                    ),
                }

            }
        }


M src/util/useful.rs => src/util/useful.rs +11 -5
@@ 75,13 75,16 @@ pub fn choose_longer_vec(x1: &Vector, x2: &Vector) -> Vector {
    }
}

pub fn max<T>(v: Vec<T>) -> T where T: PartialOrd + Copy + Clone {
pub fn max<T>(v: Vec<T>) -> T
where
    T: PartialOrd + Copy + Clone,
{
    let l = v.len();
    if l == 1 {
        v[0]
    } else {
        let mut t = if v[0] >= v[1] { v[0] } else { v[1] };
        for i in 2 .. v.len() {
        for i in 2..v.len() {
            if v[i] > t {
                t = v[i];
            }


@@ 90,17 93,20 @@ pub fn max<T>(v: Vec<T>) -> T where T: PartialOrd + Copy + Clone {
    }
}

pub fn min<T>(v: Vec<T>) -> T where T: PartialOrd + Copy + Clone {
pub fn min<T>(v: Vec<T>) -> T
where
    T: PartialOrd + Copy + Clone,
{
    let l = v.len();
    if l == 1 {
        v[0]
    } else {
        let mut t = if v[0] <= v[1] { v[0] } else { v[1] };
        for i in 2 .. v.len() {
        for i in 2..v.len() {
            if v[i] < t {
                t = v[i];
            }
        }
        t
    }
}
\ No newline at end of file
}

M src/util/writer.rs => src/util/writer.rs +49 -50
@@ 1,12 1,11 @@
//! More convenient matrix writer

pub use self::ToWriter::{Data, Header, Path, Round};
use std::collections::HashMap;
use std::fs::File;
use std::io::Write;
use std::process::exit;
use Matrix;
use Pickle;
pub use ToWriter::{Data, Header, Path, Round};
use structure::matrix::Matrix;

#[derive(Debug, Clone, Copy, Hash, PartialOrd, PartialEq, Eq)]
pub enum ToWriter {


@@ 110,51 109,51 @@ impl SimpleWriter {
        unimplemented!()
    }

    pub fn write_pickle(&self) {
        let mut writer: Box<dyn Write>;

        // Error handling - Path
        if let Some(p) = self.to_write.get(&Path) {
            assert!(*p, "No determined path!");
        }

        // Error handling - Data
        if let Some(dat) = self.to_write.get(&Data) {
            assert!(*dat, "No inserted data!");
        }

        match File::create(self.path.clone()) {
            Ok(p) => writer = Box::new(p),
            Err(e) => {
                println!("{:?}", e);
                exit(1);
            }
        }

        if let Some(head) = self.to_write.get(&Header) {
            if *head {
                serde_pickle::to_writer(&mut writer, &self.header, true)
                    .expect("Can't write header to pickle");
            }
        }

        let mut queued = self.queue.clone().into_iter();
        let mut matrices = self.matrices.clone().into_iter();
        let mut vectors = self.vectors.clone().into_iter();

        loop {
            match queued.next() {
                Some(Queue::Matrix) => {
                    let mat = matrices.next().unwrap();
                    mat.write_pickle(&mut writer)
                        .expect("Can't insert matrices");
                }
                Some(Queue::Vector) => {
                    let vec = vectors.next().unwrap();
                    vec.write_pickle(&mut writer).expect("Can't insert vectors");
                }
                None => return,
            }
        }
    }
    // pub fn write_pickle(&self) {
    //     let mut writer: Box<dyn Write>;

    //     // Error handling - Path
    //     if let Some(p) = self.to_write.get(&Path) {
    //         assert!(*p, "No determined path!");
    //     }

    //     // Error handling - Data
    //     if let Some(dat) = self.to_write.get(&Data) {
    //         assert!(*dat, "No inserted data!");
    //     }

    //     match File::create(self.path.clone()) {
    //         Ok(p) => writer = Box::new(p),
    //         Err(e) => {
    //             println!("{:?}", e);
    //             exit(1);
    //         }
    //     }

    //     if let Some(head) = self.to_write.get(&Header) {
    //         if *head {
    //             serde_pickle::to_writer(&mut writer, &self.header, true)
    //                 .expect("Can't write header to pickle");
    //         }
    //     }

    //     let mut queued = self.queue.clone().into_iter();
    //     let mut matrices = self.matrices.clone().into_iter();
    //     let mut vectors = self.vectors.clone().into_iter();

    //     loop {
    //         match queued.next() {
    //             Some(Queue::Matrix) => {
    //                 let mat = matrices.next().unwrap();
    //                 mat.write_pickle(&mut writer)
    //                     .expect("Can't insert matrices");
    //             }
    //             Some(Queue::Vector) => {
    //                 let vec = vectors.next().unwrap();
    //                 vec.write_pickle(&mut writer).expect("Can't insert vectors");
    //             }
    //             None => return,
    //         }
    //     }
    // }
}

M tests/dual.rs => tests/dual.rs +1 -0
@@ 1,3 1,4 @@
#[macro_use]
extern crate peroxide;
use peroxide::*;


M tests/jacobian.rs => tests/jacobian.rs +3 -5
@@ 1,3 1,4 @@
#[macro_use]
extern crate peroxide;
use peroxide::*;



@@ 18,8 19,5 @@ fn test_jacobian() {
fn f(xs: Vec<Number>) -> Vec<Number> {
    let x = xs[0];
    let y = xs[1];
    vec![
        x.powi(2) * y,
        5f64 * x + y.sin(),
    ]
}
\ No newline at end of file
    vec![x.powi(2) * y, 5f64 * x + y.sin()]
}

M tests/matrix.rs => tests/matrix.rs +26 -6
@@ 1,3 1,4 @@
#[macro_use]
extern crate peroxide;
use peroxide::*;



@@ 30,9 31,8 @@ fn test_accumulation() {

#[test]
fn test_add_matrix() {
    let a = matrix!(1;101;1, 10, 10, Row);
    let b = a.fmap(|x| 2.0 * x);
    assert_eq!(a.clone() + a, b);
    let a = matrix!(1;100;1, 10, 10, Row);
    assert_eq!(&a + &a, 2f64 * a);
}

#[test]


@@ 64,7 64,17 @@ fn test_print() {
fn test_row_map() {
    let m = ml_matrix("1 2;3 4");
    let n = m.row_map(|v| v.normalize());
    let o = matrix(vec![1f64/5f64.sqrt(), 2f64/5f64.sqrt(), 3f64 / 5f64, 4f64 / 5f64], 2, 2, Row);
    let o = matrix(
        vec![
            1f64 / 5f64.sqrt(),
            2f64 / 5f64.sqrt(),
            3f64 / 5f64,
            4f64 / 5f64,
        ],
        2,
        2,
        Row,
    );
    assert_eq!(n, o);
}



@@ 72,6 82,16 @@ fn test_row_map() {
fn test_col_map() {
    let m = ml_matrix("1 3;2 4");
    let n = m.col_map(|v| v.normalize());
    let o = matrix(vec![1f64/5f64.sqrt(), 2f64/5f64.sqrt(), 3f64 / 5f64, 4f64 / 5f64], 2, 2, Col);
    let o = matrix(
        vec![
            1f64 / 5f64.sqrt(),
            2f64 / 5f64.sqrt(),
            3f64 / 5f64,
            4f64 / 5f64,
        ],
        2,
        2,
        Col,
    );
    assert_eq!(n, o);
}
\ No newline at end of file
}

M tests/number.rs => tests/number.rs +1 -1
@@ 1,5 1,5 @@
#[macro_use]
extern crate peroxide;
use peroxide::operation::number::Number;
use peroxide::operation::number::Number::{D, F};
use peroxide::*;