c3d645eecb283bf77ec6602203352644967f1223 — axect 2 months ago 9f47f90
Ver 0.15.3 - SIMD meet GL4
M Cargo.toml => Cargo.toml +1 -1
@@ 1,6 1,6 @@
 [package]
 name = "peroxide"
-version = "0.15.2"
+version = "0.15.3"
 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 RELEASES.md => RELEASES.md +2 -1
@@ 1,7 1,8 @@
-# Release 0.15.3 (2019-09-08) (Candidates)
+# Release 0.15.3 (2019-09-08)
 
 * Smart pointer of Vector - `RedoxVector`
 * SIMD integrated - `packed_simd`
+* Revive Gauss-Legendre 4th order ODE solver
 
 # Release 0.15.2 (2019-09-04)
 

A example/ode_gl4.rs => example/ode_gl4.rs +39 -0
@@ 0,0 1,39 @@
+extern crate peroxide;
+use peroxide::*;
+
+fn main() {
+    let mut im_test = ImplicitODE::new(test_fn);
+
+    let init_state = State::<f64>::new(0f64, c!(1), c!(0));
+    im_test
+        .set_initial_condition(init_state)
+        .set_method(ImMethod::GL4)
+        .set_step_size(0.01f64)
+        .set_rtol(1e-6)
+        .set_times(1000);
+
+    let result = im_test.integrate();
+
+    let x = result.col(0);
+    let y = result.col(1);
+
+    #[cfg(feature = "plot")]
+    {
+        let mut plt = Plot2D::new();
+        plt.set_domain(x)
+            .insert_image(y)
+            .set_title("Test Figure")
+            .set_fig_size((10, 6))
+            .set_dpi(300)
+            .set_legend(vec!["GL4"])
+            .set_path("example_data/gl4_plot.png");
+        plt.savefig();
+    }
+}
+
+fn test_fn(st: &mut State<Dual>) {
+    let x = st.param;
+    let y = &st.value;
+    let dy = &mut st.deriv;
+    dy[0] = (5f64 * x.powi(2) - y[0]) / (x + y[0]).exp();
+}

A example_data/gl4_plot.png => example_data/gl4_plot.png +0 -0

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

D src/bin/dot.rs => src/bin/dot.rs +0 -9
@@ 1,9 0,0 @@
-extern crate peroxide;
-use peroxide::*;
-
-fn main() {
-    let a = vec![0f64; 1000_0000];
-    let b = vec![0f64; 1000_0000];
-
-    a.dot(&b).print();
-}>
\ No newline at end of file

D src/bin/s_add.rs => src/bin/s_add.rs +0 -7
@@ 1,7 0,0 @@
-extern crate peroxide;
-use peroxide::*;
-
-fn main() {
-    let a = vec![0f64; 1000_0000];
-    a.s_add(1f64).print();
-}>
\ No newline at end of file

D src/bin/scal.rs => src/bin/scal.rs +0 -8
@@ 1,8 0,0 @@
-extern crate peroxide;
-use peroxide::*;
-
-fn main() {
-    let a = vec![0f64; 1000_0000];
-    let c = 1f64;
-    a.s_mul(c).print();
-}>
\ No newline at end of file

M src/numerical/gauss_legendre.rs => src/numerical/gauss_legendre.rs +125 -125
@@ 1,125 1,125 @@
-//use numerical::utils::jacobian;
-//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.fmap(|x| x * GL4[0][1] * h)
-//                        .add(&k2.fmap(|x| x * GL4[0][2] * h)),
-//                ),
-//            ),
-//            f(
-//                t2,
-//                yn.add(
-//                    &k1.fmap(|x| x * GL4[1][1] * h)
-//                        .add(&k2.fmap(|x| x * GL4[1][2] * h)),
-//                ),
-//            ),
-//        )
-//    };
-//
-//    // 2. Obtain Jacobian
-//    let I = eye(2 * n);
-//
-//    let mut Dg = jacobian(k_curr.clone(), g.clone());
-//    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(k_curr.clone(), g.clone());
-//        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))
-//}
+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 +2 -2
@@ 1,9 1,9 @@
 //! Differential equations & Numerical Analysis tools
 
 //pub mod bdf;
-//pub mod gauss_legendre;
+pub mod gauss_legendre;
 pub mod interp;
-//pub mod newton;
+pub mod newton;
 pub mod ode;
 pub mod optimize;
 pub mod spline;

M src/numerical/newton.rs => src/numerical/newton.rs +33 -34
@@ 1,34 1,33 @@
-//use numerical::utils::jacobian;
-//use structure::dual::*;
-///// Newton-Raphson Method
-//use structure::matrix::*;
-//use structure::vector::*;
-//
-//pub fn newton<F>(init_cond: Vec<f64>, f: F, rtol: f64) -> Vec<f64>
-//where
-//    F: Fn(Vec<Dual>) -> Vec<Dual> + Copy,
-//{
-//    let mut x = init_cond;
-//    let mut x_next = update(x.clone(), f);
-//    let mut err = (x_next.sub(&x)).norm();
-//
-//    while err >= rtol {
-//        x = x_next.clone();
-//        x_next = update(x_next.clone(), f);
-//        err = (x_next.sub(&x)).norm();
-//    }
-//
-//    x_next
-//}
-//
-//fn update<F>(xs: Vec<f64>, f: F) -> Vec<f64>
-//where
-//    F: Fn(Vec<Dual>) -> Vec<Dual> + Copy,
-//{
-//    let j = jacobian(xs.clone(), f);
-//    let pinv_j = j.pseudo_inv().unwrap();
-//    let xs_dual = xs.conv_dual();
-//    let fx = f(xs_dual).values();
-//
-//    xs.sub(&(pinv_j * fx).col(0))
-//}
+use numerical::utils::jacobian;
+use operation::number::{Number, NumberVector};
+/// Newton-Raphson Method
+use structure::matrix::*;
+use structure::vector::*;
+
+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 err = (x_next.sub(&x)).norm();
+
+    while err >= rtol {
+        x = x_next.clone();
+        x_next = update(&x_next, f);
+        err = (x_next.sub(&x)).norm();
+    }
+
+    x_next
+}
+
+fn update<F>(xs: &Vec<f64>, f: F) -> Vec<f64>
+where
+    F: Fn(Vec<Number>) -> Vec<Number> + Copy,
+{
+    let j = jacobian(f, &xs);
+    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))
+}

M src/numerical/ode.rs => src/numerical/ode.rs +306 -36
@@ 205,15 205,17 @@
 use self::BoundaryCondition::Dirichlet;
 use self::ExMethod::{Euler, RK4};
 use self::ODEOptions::{BoundCond, InitCond, Method, StepSize, StopCond, Times};
+use self::ImMethod::{BDF1, GL4};
 use operation::extra_ops::Real;
 use operation::mut_ops::MutFP;
 use std::collections::HashMap;
 use structure::dual::Dual;
-use structure::matrix::{Matrix, Row, FP};
+use structure::matrix::{Matrix, Row, FP, LinearAlgebra};
 use structure::vector::FPVector;
-use util::non_macro::{cat, zeros};
+use numerical::utils::jacobian_dual;
+use util::non_macro::{cat, concat, zeros, eye};
 use util::print::Printable;
-use VecOps;
+use {VecOps, VecWithDual, Dualist};
 #[cfg(feature = "oxidize")]
 use {blas_daxpy, blas_daxpy_return};
 


@@ 227,6 229,12 @@ pub enum ExMethod {
     RK4,
 }
 
+#[derive(Debug, Copy, Clone, Hash, PartialOrd, PartialEq, Eq)]
+pub enum ImMethod {
+    BDF1,
+    GL4,
+}
+
 /// Kinds of Boundary Conditions
 ///
 /// * Dirichlet


@@ 421,48 429,18 @@ impl ODE for ExplicitODE {
                 let k1 = self.state.deriv.clone();
                 let k1_add = k1.s_mul(h2);
                 self.state.param += h2;
-
-                match () {
-                    #[cfg(feature = "oxidize")]
-                    () => {
-                        blas_daxpy(1f64, &k1_add, &mut self.state.value);
-                    }
-                    _ => {
-                        self.state.value.mut_zip_with(|x, y| x + y, &k1_add);
-                    }
-                }
-
+                self.state.value.mut_zip_with(|x, y| x + y, &k1_add);
                 (self.func)(&mut self.state);
 
                 let k2 = self.state.deriv.clone();
                 let k2_add = k2.zip_with(|x, y| h2 * x - y, &k1_add);
-
-                match () {
-                    #[cfg(feature = "oxidize")]
-                    () => {
-                        blas_daxpy(1f64, &k2_add, &mut self.state.value);
-                    }
-                    _ => {
-                        self.state.value.mut_zip_with(|x, y| x + y, &k2_add);
-                    }
-                }
-
+                self.state.value.mut_zip_with(|x, y| x + y, &k2_add);
                 (self.func)(&mut self.state);
 
                 let k3 = self.state.deriv.clone();
                 let k3_add = k3.zip_with(|x, y| h * x - y, &k2_add);
                 self.state.param += h2;
-
-                match () {
-                    #[cfg(feature = "oxidize")]
-                    () => {
-                        blas_daxpy(1f64, &k3_add, &mut self.state.value);
-                    }
-                    _ => {
-                        self.state.value.mut_zip_with(|x, y| x + y, &k3_add);
-                    }
-                }
-
+                self.state.value.mut_zip_with(|x, y| x + y, &k3_add);
                 (self.func)(&mut self.state);
 
                 let k4 = self.state.deriv.clone();


@@ 618,3 596,295 @@ impl ODE for ExplicitODE {
         true
     }
 }
+
+#[derive(Clone)]
+pub struct ImplicitODE {
+    state: State<Dual>,
+    func: fn(&mut State<Dual>),
+    step_size: f64,
+    rtol: f64,
+    method: ImMethod,
+    init_cond: State<f64>,
+    bound_cond1: (State<f64>, BoundaryCondition),
+    bound_cond2: (State<f64>, BoundaryCondition),
+    stop_cond: fn(&Self) -> bool,
+    times: usize,
+    options: HashMap<ODEOptions, bool>,
+}
+
+impl ImplicitODE {
+    pub fn new(f: ImUpdater) -> Self {
+        let mut default_to_use: HashMap<ODEOptions, bool> = HashMap::new();
+        default_to_use.insert(InitCond, false);
+        default_to_use.insert(StepSize, false);
+        default_to_use.insert(BoundCond, false);
+        default_to_use.insert(Method, false);
+        default_to_use.insert(StopCond, false);
+        default_to_use.insert(Times, false);
+
+        ImplicitODE {
+            state: Default::default(),
+            func: f,
+            step_size: 0.0,
+            rtol: 1e-6,
+            method: GL4,
+            init_cond: Default::default(),
+            bound_cond1: (Default::default(), Dirichlet),
+            bound_cond2: (Default::default(), Dirichlet),
+            stop_cond: |_x| false,
+            times: 0,
+            options: default_to_use,
+        }
+    }
+
+    pub fn get_state(&self) -> &State<Dual> {
+        &self.state
+    }
+
+    pub fn set_rtol(&mut self, rtol: f64) -> &mut Self {
+        self.rtol = rtol;
+        self
+    }
+}
+
+/// Value of 3f64.sqrt()
+const SQRT3: f64 = 1.7320508075688772;
+
+/// Butcher tableau for Gauss_legendre 4th order
+const GL4_TAB: [[f64; 3]; 2] = [
+    [0.5 - SQRT3 / 6f64, 0.25, 0.25 - SQRT3 / 6f64],
+    [0.5 + SQRT3 / 6f64, 0.25 + SQRT3 / 6f64, 0.25],
+];
+
+#[allow(non_snake_case)]
+impl ODE for ImplicitODE {
+    type Records = Matrix;
+    type Param = Dual;
+    type ODEMethod = ImMethod;
+    fn mut_update(&mut self) { 
+        match self.method {
+            BDF1 => {
+                unimplemented!()
+            }
+            GL4 => {
+                let f = |t: Dual, y: Vec<Dual>| {
+                    let mut st = State::new(t, y.clone(), y);
+                    (self.func)(&mut st);
+                    st.deriv
+                };
+
+                let h = self.step_size;
+                let t = self.state.param;
+                let t1: Dual = t + GL4_TAB[0][0] * h;
+                let t2: Dual = t + GL4_TAB[1][0] * h;
+                let yn = &self.state.value;
+                let n = yn.len();
+
+                // 0. Initial Guess
+                let k1_init: Vec<f64> = f(t, yn.clone()).values();
+                let k2_init: Vec<f64> = f(t, yn.clone()).values();
+                let mut k_curr: Vec<f64> = 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_TAB[0][1] * h)
+                                    .add(&k2.s_mul(GL4_TAB[0][2]*h)),
+                            ),
+                        ),
+                        f(
+                            t2,
+                            yn.add(
+                                &k1.s_mul(GL4_TAB[1][1] * h)
+                                    .add(&k2.s_mul(GL4_TAB[1][2] * h)),
+                            ),
+                        ),
+                    )
+                };
+
+                // 2. Obtain Jacobian
+                let I = eye(2 * n);
+
+                let mut Dg = jacobian_dual(g, &k_curr);
+                let mut DG = &I - &Dg;
+                let mut DG_inv = DG.inv().unwrap();
+                let mut G = k_curr.sub(&g(k_curr.conv_dual()).values());
+                let mut num_iter: usize = 0;
+
+                // 3. Iteration
+                while err >= self.rtol && num_iter <= 10 {
+                    let DGG = &DG_inv * &G;
+                    let k_prev = k_curr.clone();
+                    k_curr.mut_zip_with(|x, y| x - y, &DGG.col(0));
+                    Dg = jacobian_dual(g, &k_curr);
+                    DG = &I - &Dg;
+                    DG_inv = DG.inv().unwrap();
+                    G = k_curr.sub(&g(k_curr.conv_dual()).values());
+                    err = k_curr.sub(&k_prev).norm();
+                    num_iter += 1;
+                }
+
+                // 4. Find k1, k2
+                let (k1, k2) = (k_curr.take(n), k_curr.skip(n));
+
+                // Set Derivative from state
+                (self.func)(&mut self.state);
+
+                // 5. Merge k1, k2
+                let mut y_curr = self.state.value.values();
+                y_curr = y_curr.add(&k1.s_mul(0.5 * h).add(&k2.s_mul(0.5 * h)));
+                self.state.value = y_curr.conv_dual();
+                self.state.param = self.state.param + h;
+            }
+        }
+    }
+
+    fn integrate(&mut self) -> Self::Records { 
+        assert!(self.check_enough(), "Not enough fields!");
+
+        let mut result = zeros(self.times + 1, self.state.value.len() + 1);
+
+        result.subs_row(0, &cat(self.state.param.to_f64(), self.state.value.values()));
+
+        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.to_f64(), self.state.value.values()));
+                    key += 1;
+                    if (self.stop_cond)(&self) {
+                        println!("Reach the stop condition!");
+                        print!("Current values are: ");
+                        cat(self.state.param.to_f64(), self.state.value.values()).print();
+                        break;
+                    }
+                }
+                return result.take(key, Row);
+            }
+            _ => {
+                for i in 1..self.times + 1 {
+                    self.mut_update();
+                    result.subs_row(i, &cat(self.state.param.to_f64(), self.state.value.values()));
+                }
+                return result;
+            }
+        }
+    }
+
+    fn set_initial_condition<T: Real>(&mut self, init: State<T>) -> &mut Self { 
+        if let Some(x) = self.options.get_mut(&InitCond) {
+            *x = true
+        }
+        self.init_cond = init.to_f64();
+        self.state = init.to_dual();
+        self
+    }
+
+    fn set_boundary_condition<T: Real>(
+        &mut self,
+        bound1: (State<T>, BoundaryCondition),
+        bound2: (State<T>, BoundaryCondition),
+    ) -> &mut Self {
+        if let Some(x) = self.options.get_mut(&BoundCond) {
+            *x = true
+        }
+        self.bound_cond1 = (bound1.0.to_f64(), bound1.1);
+        self.bound_cond2 = (bound2.0.to_f64(), bound2.1);
+        self
+    }
+
+    fn set_step_size(&mut self, dt: f64) -> &mut Self { 
+        if let Some(x) = self.options.get_mut(&StepSize) {
+            *x = true
+        }
+        self.step_size = dt;
+        self
+    }
+
+    fn set_method(&mut self, method: Self::ODEMethod) -> &mut Self { 
+        if let Some(x) = self.options.get_mut(&Method) {
+            *x = true
+        }
+        self.method = method;
+        self
+    }
+
+    fn set_stop_condition(&mut self, f: fn(&Self) -> bool) -> &mut Self {
+        if let Some(x) = self.options.get_mut(&StopCond) {
+            *x = true
+        }
+        self.stop_cond = f;
+        self
+    }
+
+    fn set_times(&mut self, n: usize) -> &mut Self {
+        if let Some(x) = self.options.get_mut(&Times) {
+            *x = true
+        }
+        self.times = n;
+        self
+    }
+
+    fn check_enough(&self) -> bool {
+        // Method
+        match self.options.get(&Method) {
+            Some(x) => {
+                if !*x {
+                    return false;
+                }
+            }
+            None => {
+                return false;
+            }
+        }
+
+        // Step size
+        match self.options.get(&StepSize) {
+            Some(x) => {
+                if !*x {
+                    return false;
+                }
+            }
+            None => {
+                return false;
+            }
+        }
+
+        // Initial or Boundary
+        match self.options.get(&InitCond) {
+            None => {
+                return false;
+            }
+            Some(x) => {
+                if !*x {
+                    match self.options.get(&BoundCond) {
+                        None => {
+                            return false;
+                        }
+                        Some(_) => (),
+                    }
+                }
+            }
+        }
+
+        // Set Time?
+        match self.options.get(&Times) {
+            None => {
+                return false;
+            }
+            Some(x) => {
+                if !*x {
+                    return false;
+                }
+            }
+        }
+        true
+    }    
+}<
\ No newline at end of file

M src/numerical/optimize.rs => src/numerical/optimize.rs +3 -3
@@ 198,7 198,7 @@ impl Optimizer {
 
         // Declare mutable values
         let mut p = p_init_vec.to_matrix();
-        let mut j = jacobian(f, p_init_vec.clone());
+        let mut j = jacobian(f, &p_init_vec);
         let mut y_hat = f(p_init.clone()).to_f64_vec().to_matrix();
         let mut jtj = &j.t() * &j;
 


@@ 207,7 207,7 @@ impl Optimizer {
                 for _i in 0..max_iter {
                     let h = alpha * j.t() * (&y - &y_hat);
                     p = &p + &h;
-                    j = jacobian(f, p.data.clone());
+                    j = jacobian(f, &p.data);
                     y_hat = f(NumberVector::from_f64_vec(p.data.clone()))
                         .to_f64_vec()
                         .to_matrix();


@@ 231,7 231,7 @@ impl Optimizer {
                     }
 
                     let p_temp = &p + &h;
-                    let j_temp = jacobian(f, p.data.clone());
+                    let j_temp = jacobian(f, &p.data);
                     let y_hat_temp = f(NumberVector::from_f64_vec(p_temp.data.clone()))
                         .to_f64_vec()
                         .to_matrix();

M src/numerical/utils.rs => src/numerical/utils.rs +32 -4
@@ 17,7 17,7 @@ use util::non_macro::{cat, zeros};
 /// use peroxide::*;
 ///
 /// let x = c!(1, 1);
-/// let j = jacobian(f, x);
+/// let j = jacobian(f, &x);
 /// j.print();
 ///
 /// //      c[0] c[1]


@@ 35,7 35,7 @@ use util::non_macro::{cat, zeros};
 /// }
 /// ```
 #[allow(non_snake_case)]
-pub fn jacobian<F>(g: F, x: Vec<f64>) -> Matrix
+pub fn jacobian<F>(g: F, x: &Vec<f64>) -> Matrix
 where
     F: Fn(Vec<Number>) -> Vec<Number>,
 {


@@ 43,8 43,36 @@ where
 
     let f = |x: Vec<Dual>| g(NumberVector::from_dual_vec(x)).to_dual_vec();
 
-    let x_var: Vec<Dual> = merge_dual(x.clone(), vec![1f64; l]);
-    let x_const = x.clone().conv_dual();
+    let x_var: Vec<Dual> = merge_dual(&x, &vec![1f64; l]);
+    let x_const = x.conv_dual();
+
+    let l2 = f(x_const.clone()).len();
+
+    let mut J = zeros(l2, l);
+
+    let mut x_temp = x_const.clone();
+
+    for i in 0..l {
+        x_temp[i] = x_var[i];
+        let dual_temp = f(x_temp.clone());
+        let slope_temp = dual_temp.slopes();
+        for j in 0..l2 {
+            J[(j, i)] = slope_temp[j];
+        }
+        x_temp = x_const.clone();
+    }
+    J
+}
+
+#[allow(non_snake_case)]
+pub fn jacobian_dual<F>(f: F, x: &Vec<f64>) -> Matrix
+where
+    F: Fn(Vec<Dual>) -> Vec<Dual>,
+{
+    let l = x.len();
+
+    let x_var: Vec<Dual> = merge_dual(&x, &vec![1f64; l]);
+    let x_const = x.conv_dual();
 
     let l2 = f(x_const.clone()).len();
 

M src/redox/redoxable.rs => src/redox/redoxable.rs +19 -30
@@ 1,70 1,59 @@
 use std::ops::{Add, Deref, Index, IndexMut};
+use structure::vector::VecOps;
 
 /// Smart pointer of Vector
-pub struct RedoxVector<T> {
-    pub data: Vec<T>
+pub struct RedoxVector {
+    pub data: Vec<f64>
 }
 
-impl<T: Default + Clone> RedoxVector<T> {
+impl RedoxVector {
     pub fn new(n: usize) -> Self {
         RedoxVector {
-            data: vec![T::default(); n]
+            data: vec![f64::default(); n]
         }
     }
 
-    pub fn from_vec(v: Vec<T>) -> Self {
+    pub fn from_vec(v: Vec<f64>) -> Self {
         RedoxVector {
             data: v
         }
     }
 }
 
-impl<T> Deref for RedoxVector<T> {
-    type Target = Vec<T>;
+impl Deref for RedoxVector {
+    type Target = Vec<f64>;
 
     fn deref(&self) -> &Self::Target {
         &self.data
     }
 }
 
-impl<T> Index<usize> for RedoxVector<T> {
-    type Output = T;
+impl Index<usize> for RedoxVector {
+    type Output = f64;
 
     fn index(&self, index: usize) -> &Self::Output {
         self.data.index(index)
     }
 }
 
-impl<T> IndexMut<usize> for RedoxVector<T> {
-    fn index_mut(&mut self, index: usize) -> &mut T {
+impl IndexMut<usize> for RedoxVector {
+    fn index_mut(&mut self, index: usize) -> &mut f64 {
         self.data.index_mut(index)
     }
 }
 
-impl<T> Add<RedoxVector<T>> for RedoxVector<T>
-where T: Default + Add<Output=T> + Copy + Clone
-{
-    type Output = RedoxVector<T>;
+impl Add<RedoxVector> for RedoxVector {
+    type Output = RedoxVector;
 
     fn add(self, rhs: Self) -> Self::Output {
-        let mut rv: RedoxVector<T> = RedoxVector::new((*self).len());
-        for i in 0 .. rv.len() {
-            rv[i] = self[i] + rhs[i];
-        }
-        rv
+        RedoxVector::from_vec(self.data.add(&rhs.data))
     }
 }
 
-impl<'a, 'b, T> Add<&'b RedoxVector<T>> for &'a RedoxVector<T>
-where T: Default + Add<Output=T> + Copy + Clone
-{
-    type Output = RedoxVector<T>;
+impl<'a, 'b> Add<&'b RedoxVector> for &'a RedoxVector {
+    type Output = RedoxVector;
 
-    fn add(self, rhs: &'b RedoxVector<T>) -> Self::Output {
-        let mut rv: RedoxVector<T> = RedoxVector::new((*self).len());
-        for i in 0 .. rv.len() {
-            rv[i] = self[i] + rhs[i];
-        }
-        rv
+    fn add(self, rhs: &'b RedoxVector) -> Self::Output {
+        RedoxVector::from_vec(self.data.add(&rhs.data))
     }
 }=
\ No newline at end of file

M src/structure/dual.rs => src/structure/dual.rs +7 -5
@@ 539,11 539,13 @@ impl Dualist for Vec<Dual> {
     }
 }
 
-pub fn merge_dual(y: Vec<f64>, dy: Vec<f64>) -> Vec<Dual> {
-    y.into_iter()
-        .zip(&dy)
-        .map(|(t, &dt)| Dual::new(t, dt))
-        .collect::<Vec<Dual>>()
+pub fn merge_dual(y: &Vec<f64>, dy: &Vec<f64>) -> Vec<Dual> {
+    let mut result = vec![dual(0, 0); y.len()];
+
+    for i in 0 .. y.len() {
+        result[i] = dual(y[i], dy[i]);
+    }
+    result
 }
 
 impl FPVector for Vec<Dual> {

M tests/jacobian.rs => tests/jacobian.rs +1 -1
@@ 4,7 4,7 @@ use peroxide::*;
 #[test]
 fn test_jacobian() {
     let x = c!(1, 0);
-    let j = jacobian(f, x);
+    let j = jacobian(f, &x);
     assert_eq!(j, ml_matrix("0 1; 5 1"));
 }