Ver 0.39.0

Decouple initial conditions from ODEProblem

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "peroxide"
-version = "0.38.1"
+version = "0.39.0"
 authors = ["axect <[email protected]>"]
 edition = "2018"
 description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"
@@ -11,12 +11,12 @@ readme = "README.md"
 documentation = "https://axect.github.io/Peroxide_Doc"
 keywords = ["Numeric", "Science", "Dataframe", "Plot", "LinearAlgebra"]
 exclude = [
-    "example_data/",
-    "src/bin/",
-    "benches/",
-    "example/",
-    "test_data/",
-    "peroxide-ad2",
+  "example_data/",
+  "src/bin/",
+  "benches/",
+  "example/",
+  "test_data/",
+  "peroxide-ad2",
 ]
 
 [badges]
@@ -40,12 +40,18 @@ anyhow = "1.0"
 paste = "1.0"
 #num-complex = "0.3"
 netcdf = { version = "0.7", optional = true, default-features = false }
-pyo3 = { version = "0.22", optional = true, features = ["auto-initialize", "gil-refs"] }
+pyo3 = { version = "0.23", optional = true, features = [
+  "auto-initialize",
+  "gil-refs",
+] }
 blas = { version = "0.22", optional = true }
 lapack = { version = "0.19", optional = true }
 serde = { version = "1.0", features = ["derive"], optional = true }
 json = { version = "0.12", optional = true }
-arrow2 = { version = "0.18", features = ["io_parquet", "io_parquet_compression"], optional = true }
+arrow2 = { version = "0.18", features = [
+  "io_parquet",
+  "io_parquet_compression",
+], optional = true }
 num-complex = { version = "0.4", optional = true }
 rayon = { version = "1.10", optional = true }
 

README.md

@@ -260,12 +260,14 @@ The example code demonstrates how Peroxide can be used to simulate the Lorenz at
 use peroxide::fuga::*;
 
 fn main() -> Result<(), Box<dyn Error>> {
+    let initial_conditions = vec![10f64, 1f64, 1f64];
     let rkf45 = RKF45::new(1e-4, 0.9, 1e-6, 1e-2, 100);
     let basic_ode_solver = BasicODESolver::new(rkf45);
     let (_, y_vec) = basic_ode_solver.solve(
         &Lorenz,
         (0f64, 100f64),
         1e-2,
+        &initial_conditions,
     )?; // Error handling with `?` - can check constraint violation and etc.
     let y_mat = py_matrix(y_vec);
     let y0 = y_mat.col(0);
@@ -290,10 +292,6 @@ fn main() -> Result<(), Box<dyn Error>> {
 struct Lorenz;
 
 impl ODEProblem for Lorenz {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![10f64, 1f64, 1f64]
-    }
-
     fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         dy[0] = 10f64 * (y[1] - y[0]);
         dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];

RELEASES.md

@@ -1,3 +1,8 @@
+# Release 0.39.0 (2024-11-22)
+
+- Decouple `initial_conditions` from `ODEProblem`
+  - Now, we can define `initial_conditions` in solving phase 
+
 # Release 0.38.1 (2024-11-06)
 
 - Fix error in `O3` feature

examples/lorenz_dp45.rs

@@ -2,9 +2,10 @@ use peroxide::fuga::*;
 
 #[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
+    let initial_conditions = vec![10f64, 1f64, 1f64];
     let dp45 = DP45::new(1e-4, 0.9, 1e-6, 1e-1, 100);
     let basic_ode_solver = BasicODESolver::new(dp45);
-    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2)?;
+    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2, &initial_conditions)?;
     let y_mat = py_matrix(y_vec);
     let y0 = y_mat.col(0);
     let y2 = y_mat.col(2);
@@ -29,10 +30,6 @@ fn main() -> Result<(), Box<dyn Error>> {
 struct Lorenz;
 
 impl ODEProblem for Lorenz {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![10f64, 1f64, 1f64]
-    }
-
     fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         dy[0] = 10f64 * (y[1] - y[0]);
         dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];

examples/lorenz_gl4.rs

@@ -2,9 +2,10 @@ use peroxide::fuga::*;
 
 #[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
+    let initial_conditions = vec![10f64, 1f64, 1f64];
     let gl4 = GL4::new(ImplicitSolver::FixedPoint, 1e-6, 100);
     let basic_ode_solver = BasicODESolver::new(gl4);
-    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2)?;
+    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2, &initial_conditions)?;
     let y_mat = py_matrix(y_vec);
     let y0 = y_mat.col(0);
     let y2 = y_mat.col(2);
@@ -29,10 +30,6 @@ fn main() -> Result<(), Box<dyn Error>> {
 struct Lorenz;
 
 impl ODEProblem for Lorenz {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![10f64, 1f64, 1f64]
-    }
-
     fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         dy[0] = 10f64 * (y[1] - y[0]);
         dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];

examples/lorenz_rk4.rs

@@ -2,8 +2,9 @@ use peroxide::fuga::*;
 
 #[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
+    let initial_conditions = vec![10f64, 1f64, 1f64];
     let basic_ode_solver = BasicODESolver::new(RK4);
-    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2)?;
+    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2, &initial_conditions)?;
     let y_mat = py_matrix(y_vec);
     let y0 = y_mat.col(0);
     let y2 = y_mat.col(2);
@@ -28,10 +29,6 @@ fn main() -> Result<(), Box<dyn Error>> {
 struct Lorenz;
 
 impl ODEProblem for Lorenz {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![10f64, 1f64, 1f64]
-    }
-
     fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         dy[0] = 10f64 * (y[1] - y[0]);
         dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];

examples/lorenz_rkf45.rs

@@ -2,9 +2,10 @@ use peroxide::fuga::*;
 
 #[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
+    let initial_conditions = vec![10f64, 1f64, 1f64];
     let rkf45 = RKF45::new(1e-4, 0.9, 1e-6, 1e-2, 100);
     let basic_ode_solver = BasicODESolver::new(rkf45);
-    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2)?;
+    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2, &initial_conditions)?;
     let y_mat = py_matrix(y_vec);
     let y0 = y_mat.col(0);
     let y2 = y_mat.col(2);
@@ -29,10 +30,6 @@ fn main() -> Result<(), Box<dyn Error>> {
 struct Lorenz;
 
 impl ODEProblem for Lorenz {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![10f64, 1f64, 1f64]
-    }
-
     fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         dy[0] = 10f64 * (y[1] - y[0]);
         dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];

examples/lorenz_tsit45.rs

@@ -2,9 +2,10 @@ use peroxide::fuga::*;
 
 #[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
+    let initial_conditions = vec![10f64, 1f64, 1f64];
     let tsit45 = TSIT45::new(1e-2, 0.9, 1e-6, 1e-1, 100);
     let basic_ode_solver = BasicODESolver::new(tsit45);
-    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2)?;
+    let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2, &initial_conditions)?;
     let y_mat = py_matrix(y_vec);
     let y0 = y_mat.col(0);
     let y2 = y_mat.col(2);
@@ -29,10 +30,6 @@ fn main() -> Result<(), Box<dyn Error>> {
 struct Lorenz;
 
 impl ODEProblem for Lorenz {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![10f64, 1f64, 1f64]
-    }
-
     fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         dy[0] = 10f64 * (y[1] - y[0]);
         dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];

examples/ode_env.rs

@@ -11,9 +11,11 @@ fn main() -> Result<(), Box<dyn Error>> {
 
     let c = cubic_hermite_spline(&t, &y, Quadratic)?;
 
+    let initial_conditions = vec![1f64];
     let test_problem = Test { cs: c };
     let basic_ode_solver = BasicODESolver::new(RK4);
-    let (t_vec, y_vec) = basic_ode_solver.solve(&test_problem, (0f64, 10f64), 0.01)?;
+    let (t_vec, y_vec) =
+        basic_ode_solver.solve(&test_problem, (0f64, 10f64), 0.01, &initial_conditions)?;
     let y_vec: Vec<f64> = y_vec.into_iter().flatten().collect();
 
     #[cfg(feature = "plot")]
@@ -38,10 +40,6 @@ struct Test {
 }
 
 impl ODEProblem for Test {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![1f64]
-    }
-
     fn rhs(&self, t: f64, _y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         Ok(dy[0] = self.cs.eval(t))
     }

examples/ode_test_gl4.rs

@@ -2,9 +2,10 @@ use peroxide::fuga::*;
 
 #[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
+    let initial_conditions = vec![1f64];
     let gl4 = GL4::new(ImplicitSolver::FixedPoint, 1e-6, 100);
     let basic_ode_solver = BasicODESolver::new(gl4);
-    let (t_vec, y_vec) = basic_ode_solver.solve(&Test, (0f64, 10f64), 0.01)?;
+    let (t_vec, y_vec) = basic_ode_solver.solve(&Test, (0f64, 10f64), 0.01, &initial_conditions)?;
     let y_vec: Vec<f64> = y_vec.into_iter().flatten().collect();
 
     #[cfg(feature = "plot")]
@@ -27,10 +28,6 @@ fn main() -> Result<(), Box<dyn Error>> {
 struct Test;
 
 impl ODEProblem for Test {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![1f64]
-    }
-
     fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         Ok(dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp())
     }

examples/ode_test_rk4.rs

@@ -2,8 +2,9 @@ use peroxide::fuga::*;
 
 #[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
+    let initial_conditions = vec![1f64];
     let basic_ode_solver = BasicODESolver::new(RK4);
-    let (t_vec, y_vec) = basic_ode_solver.solve(&Test, (0f64, 10f64), 1e-3)?;
+    let (t_vec, y_vec) = basic_ode_solver.solve(&Test, (0f64, 10f64), 1e-3, &initial_conditions)?;
     let y_vec: Vec<f64> = y_vec.into_iter().flatten().collect();
 
     #[cfg(feature = "plot")]
@@ -26,10 +27,6 @@ fn main() -> Result<(), Box<dyn Error>> {
 struct Test;
 
 impl ODEProblem for Test {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![1f64]
-    }
-
     fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         Ok(dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp())
     }

examples/ode_test_rkf45.rs

@@ -2,9 +2,10 @@ use peroxide::fuga::*;
 
 #[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
+    let initial_conditions = vec![1f64];
     let rkf = RKF45::new(1e-4, 0.9, 1e-6, 1e-1, 100);
     let basic_ode_solver = BasicODESolver::new(rkf);
-    let (t_vec, y_vec) = basic_ode_solver.solve(&Test, (0f64, 10f64), 0.01)?;
+    let (t_vec, y_vec) = basic_ode_solver.solve(&Test, (0f64, 10f64), 0.01, &initial_conditions)?;
     let y_vec: Vec<f64> = y_vec.into_iter().flatten().collect();
 
     #[cfg(feature = "plot")]
@@ -27,10 +28,6 @@ fn main() -> Result<(), Box<dyn Error>> {
 struct Test;
 
 impl ODEProblem for Test {
-    fn initial_conditions(&self) -> Vec<f64> {
-        vec![1f64]
-    }
-
     fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
         Ok(dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp())
     }

src/complex/matrix.rs

@@ -16,7 +16,7 @@ use crate::{
     traits::fp::{FPMatrix, FPVector},
     traits::general::Algorithm,
     traits::math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector},
-    traits::matrix::{Form, LinearAlgebra, MatrixTrait, SolveKind, PQLU, QR, SVD, WAZD, UPLO},
+    traits::matrix::{Form, LinearAlgebra, MatrixTrait, SolveKind, PQLU, QR, SVD, UPLO, WAZD},
     traits::mutable::MutMatrix,
     util::low_level::{copy_vec_ptr, swap_vec_ptr},
     util::non_macro::ConcatenateError,

src/numerical/ode.rs

@@ -48,10 +48,12 @@
 //!         max_step_iter: 100,
 //!     };
 //!     let basic_ode_solver = BasicODESolver::new(rkf);
+//!     let initial_conditions = vec![1f64];
 //!     let (t_vec, y_vec) = basic_ode_solver.solve(
 //!         &Test,
 //!         (0f64, 10f64),
 //!         0.01,
+//!         &initial_conditions,
 //!     )?;
 //!     let y_vec: Vec<f64> = y_vec.into_iter().flatten().collect();
 //!     println!("{}", y_vec.len());
@@ -77,10 +79,6 @@
 //! struct Test;
 //!
 //! impl ODEProblem for Test {
-//!     fn initial_conditions(&self) -> Vec<f64> {
-//!         vec![1f64]
-//!     }
-//!
 //!     fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
 //!         Ok(dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp())
 //!     }
@@ -101,10 +99,6 @@ use anyhow::{bail, Result};
 /// struct MyODEProblem;
 ///
 /// impl ODEProblem for MyODEProblem {
-///     fn initial_conditions(&self) -> Vec<f64> {
-///         vec![1.0, 2.0]
-///     }
-///
 ///     fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
 ///         dy[0] = -0.5 * y[0];
 ///         dy[1] = y[0] - y[1];
@@ -113,7 +107,6 @@ use anyhow::{bail, Result};
 /// }
 /// ```
 pub trait ODEProblem {
-    fn initial_conditions(&self) -> Vec<f64>;
     fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<()>;
 }
 
@@ -143,7 +136,6 @@ pub trait ODEIntegrator {
 /// }
 ///
 /// impl ODEProblem for ConstrainedProblem {
-///     fn initial_conditions(&self) -> Vec<f64> { vec![0.0] } // y_0 = 0
 ///     fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
 ///         if y[0] < self.y_constraint {
 ///             anyhow::bail!(ODEError::ConstraintViolation(t, y.to_vec(), dy.to_vec()));
@@ -182,6 +174,7 @@ pub trait ODESolver {
         problem: &P,
         t_span: (f64, f64),
         dt: f64,
+        initial_conditions: &[f64],
     ) -> Result<(Vec<f64>, Vec<Vec<f64>>)>;
 }
 
@@ -193,12 +186,14 @@ pub trait ODESolver {
 /// use peroxide::fuga::*;
 ///
 /// fn main() -> Result<(), Box<dyn Error>> {
+///     let initial_conditions = vec![1f64];
 ///     let rkf = RKF45::new(1e-4, 0.9, 1e-6, 1e-1, 100);
 ///     let basic_ode_solver = BasicODESolver::new(rkf);
 ///     let (t_vec, y_vec) = basic_ode_solver.solve(
 ///         &Test,
 ///         (0f64, 10f64),
 ///         0.01,
+///         &initial_conditions,
 ///     )?;
 ///     let y_vec: Vec<f64> = y_vec.into_iter().flatten().collect();
 ///
@@ -208,10 +203,6 @@ pub trait ODESolver {
 /// struct Test;
 ///
 /// impl ODEProblem for Test {
-///     fn initial_conditions(&self) -> Vec<f64> {
-///         vec![1f64]
-///     }
-///
 ///     fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
 ///         dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp();
 ///         Ok(())
@@ -234,10 +225,11 @@ impl<I: ODEIntegrator> ODESolver for BasicODESolver<I> {
         problem: &P,
         t_span: (f64, f64),
         dt: f64,
+        initial_conditions: &[f64],
     ) -> Result<(Vec<f64>, Vec<Vec<f64>>)> {
         let mut t = t_span.0;
         let mut dt = dt;
-        let mut y = problem.initial_conditions();
+        let mut y = initial_conditions.to_vec();
         let mut t_vec = vec![t];
         let mut y_vec = vec![y.clone()];