From 345f764ea5dc1e8c1d8d0ba1774e14a199c3b3b1 Mon Sep 17 00:00:00 2001 From: Axect Date: Fri, 22 Nov 2024 16:53:52 +0900 Subject: [PATCH 1/3] FMT --- src/complex/matrix.rs | 2 +- src/prelude/simpler.rs | 2 +- src/structure/matrix.rs | 2 +- src/structure/sparse.rs | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/complex/matrix.rs b/src/complex/matrix.rs index 37cb06e..1d05969 100644 --- a/src/complex/matrix.rs +++ b/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, diff --git a/src/prelude/simpler.rs b/src/prelude/simpler.rs index 71aa4a8..f2e70d0 100644 --- a/src/prelude/simpler.rs +++ b/src/prelude/simpler.rs @@ -12,7 +12,7 @@ use crate::structure::matrix::Matrix; use crate::structure::polynomial; use crate::traits::math::{Norm, Normed}; #[allow(unused_imports)] -use crate::traits::matrix::{Form, LinearAlgebra, MatrixTrait, SolveKind, PQLU, QR, WAZD, UPLO}; +use crate::traits::matrix::{Form, LinearAlgebra, MatrixTrait, SolveKind, PQLU, QR, UPLO, WAZD}; #[cfg(feature = "parquet")] use arrow2::io::parquet::write::CompressionOptions; #[cfg(feature = "parquet")] diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs index 5dead26..539fe75 100644 --- a/src/structure/matrix.rs +++ b/src/structure/matrix.rs @@ -624,7 +624,7 @@ use crate::traits::{ fp::{FPMatrix, FPVector}, general::Algorithm, math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector}, - matrix::{Form, LinearAlgebra, MatrixTrait, SolveKind, PQLU, QR, SVD, WAZD, UPLO}, + matrix::{Form, LinearAlgebra, MatrixTrait, SolveKind, PQLU, QR, SVD, UPLO, WAZD}, mutable::MutMatrix, }; #[allow(unused_imports)] diff --git a/src/structure/sparse.rs b/src/structure/sparse.rs index 6ce51f3..35b3576 100644 --- a/src/structure/sparse.rs +++ b/src/structure/sparse.rs @@ -5,7 +5,7 @@ use crate::structure::matrix::Matrix; use crate::traits::math::LinearOp; #[allow(unused_imports)] -use crate::traits::matrix::{Form, LinearAlgebra, SolveKind, PQLU, QR, SVD, WAZD, UPLO}; +use crate::traits::matrix::{Form, LinearAlgebra, SolveKind, PQLU, QR, SVD, UPLO, WAZD}; //use crate::traits::math::{InnerProduct, LinearOp, Norm, Normed, Vector}; use crate::util::non_macro::zeros; use std::ops::Mul; From f2490f2983bb0c565f8c7597c5370525c194efe3 Mon Sep 17 00:00:00 2001 From: Axect Date: Fri, 22 Nov 2024 17:08:14 +0900 Subject: [PATCH 2/3] CHGE: Decouple initial conditions with ODEProblem --- README.md | 6 ++---- examples/lorenz_dp45.rs | 7 ++----- examples/lorenz_gl4.rs | 7 ++----- examples/lorenz_rk4.rs | 7 ++----- examples/lorenz_rkf45.rs | 7 ++----- examples/lorenz_tsit45.rs | 7 ++----- examples/ode_env.rs | 8 +++----- examples/ode_test_gl4.rs | 7 ++----- examples/ode_test_rk4.rs | 7 ++----- examples/ode_test_rkf45.rs | 7 ++----- src/numerical/ode.rs | 22 +++++++--------------- 11 files changed, 28 insertions(+), 64 deletions(-) diff --git a/README.md b/README.md index abbc8cc..b244dd0 100644 --- a/README.md +++ b/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> { + 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> { struct Lorenz; impl ODEProblem for Lorenz { - fn initial_conditions(&self) -> Vec { - 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]; diff --git a/examples/lorenz_dp45.rs b/examples/lorenz_dp45.rs index ce8c8a7..87a204f 100644 --- a/examples/lorenz_dp45.rs +++ b/examples/lorenz_dp45.rs @@ -2,9 +2,10 @@ use peroxide::fuga::*; #[allow(unused_variables)] fn main() -> Result<(), Box> { + 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> { struct Lorenz; impl ODEProblem for Lorenz { - fn initial_conditions(&self) -> Vec { - 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]; diff --git a/examples/lorenz_gl4.rs b/examples/lorenz_gl4.rs index 9a64b57..96893e3 100644 --- a/examples/lorenz_gl4.rs +++ b/examples/lorenz_gl4.rs @@ -2,9 +2,10 @@ use peroxide::fuga::*; #[allow(unused_variables)] fn main() -> Result<(), Box> { + 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> { struct Lorenz; impl ODEProblem for Lorenz { - fn initial_conditions(&self) -> Vec { - 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]; diff --git a/examples/lorenz_rk4.rs b/examples/lorenz_rk4.rs index 8ebc458..8af7876 100644 --- a/examples/lorenz_rk4.rs +++ b/examples/lorenz_rk4.rs @@ -2,8 +2,9 @@ use peroxide::fuga::*; #[allow(unused_variables)] fn main() -> Result<(), Box> { + 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> { struct Lorenz; impl ODEProblem for Lorenz { - fn initial_conditions(&self) -> Vec { - 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]; diff --git a/examples/lorenz_rkf45.rs b/examples/lorenz_rkf45.rs index da0453b..3ccbff1 100644 --- a/examples/lorenz_rkf45.rs +++ b/examples/lorenz_rkf45.rs @@ -2,9 +2,10 @@ use peroxide::fuga::*; #[allow(unused_variables)] fn main() -> Result<(), Box> { + 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> { struct Lorenz; impl ODEProblem for Lorenz { - fn initial_conditions(&self) -> Vec { - 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]; diff --git a/examples/lorenz_tsit45.rs b/examples/lorenz_tsit45.rs index 5bc147f..5076d0d 100644 --- a/examples/lorenz_tsit45.rs +++ b/examples/lorenz_tsit45.rs @@ -2,9 +2,10 @@ use peroxide::fuga::*; #[allow(unused_variables)] fn main() -> Result<(), Box> { + 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> { struct Lorenz; impl ODEProblem for Lorenz { - fn initial_conditions(&self) -> Vec { - 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]; diff --git a/examples/ode_env.rs b/examples/ode_env.rs index e2836bb..cbcf3b9 100644 --- a/examples/ode_env.rs +++ b/examples/ode_env.rs @@ -11,9 +11,11 @@ fn main() -> Result<(), Box> { 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 = y_vec.into_iter().flatten().collect(); #[cfg(feature = "plot")] @@ -38,10 +40,6 @@ struct Test { } impl ODEProblem for Test { - fn initial_conditions(&self) -> Vec { - vec![1f64] - } - fn rhs(&self, t: f64, _y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> { Ok(dy[0] = self.cs.eval(t)) } diff --git a/examples/ode_test_gl4.rs b/examples/ode_test_gl4.rs index bf31dee..69b83c9 100644 --- a/examples/ode_test_gl4.rs +++ b/examples/ode_test_gl4.rs @@ -2,9 +2,10 @@ use peroxide::fuga::*; #[allow(unused_variables)] fn main() -> Result<(), Box> { + 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 = y_vec.into_iter().flatten().collect(); #[cfg(feature = "plot")] @@ -27,10 +28,6 @@ fn main() -> Result<(), Box> { struct Test; impl ODEProblem for Test { - fn initial_conditions(&self) -> Vec { - 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()) } diff --git a/examples/ode_test_rk4.rs b/examples/ode_test_rk4.rs index 6b90f5c..2bc10f2 100644 --- a/examples/ode_test_rk4.rs +++ b/examples/ode_test_rk4.rs @@ -2,8 +2,9 @@ use peroxide::fuga::*; #[allow(unused_variables)] fn main() -> Result<(), Box> { + 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 = y_vec.into_iter().flatten().collect(); #[cfg(feature = "plot")] @@ -26,10 +27,6 @@ fn main() -> Result<(), Box> { struct Test; impl ODEProblem for Test { - fn initial_conditions(&self) -> Vec { - 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()) } diff --git a/examples/ode_test_rkf45.rs b/examples/ode_test_rkf45.rs index 04ea96c..b08240e 100644 --- a/examples/ode_test_rkf45.rs +++ b/examples/ode_test_rkf45.rs @@ -2,9 +2,10 @@ use peroxide::fuga::*; #[allow(unused_variables)] fn main() -> Result<(), Box> { + 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 = y_vec.into_iter().flatten().collect(); #[cfg(feature = "plot")] @@ -27,10 +28,6 @@ fn main() -> Result<(), Box> { struct Test; impl ODEProblem for Test { - fn initial_conditions(&self) -> Vec { - 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()) } diff --git a/src/numerical/ode.rs b/src/numerical/ode.rs index 4025bbe..d5b477e 100644 --- a/src/numerical/ode.rs +++ b/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 = 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 { -//! 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 { -/// 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; 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 { 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, Vec>)>; } @@ -193,12 +186,14 @@ pub trait ODESolver { /// use peroxide::fuga::*; /// /// fn main() -> Result<(), Box> { +/// 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 = y_vec.into_iter().flatten().collect(); /// @@ -208,10 +203,6 @@ pub trait ODESolver { /// struct Test; /// /// impl ODEProblem for Test { -/// fn initial_conditions(&self) -> Vec { -/// 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 ODESolver for BasicODESolver { problem: &P, t_span: (f64, f64), dt: f64, + initial_conditions: &[f64], ) -> Result<(Vec, Vec>)> { 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()]; From db310b88d686d0675a2e40eb09e068b7445a8632 Mon Sep 17 00:00:00 2001 From: Axect Date: Fri, 22 Nov 2024 17:12:08 +0900 Subject: [PATCH 3/3] RLSE: Ver 0.39.0 - Decouple initial conditions from ODEProblem --- Cargo.toml | 24 +++++++++++++++--------- RELEASES.md | 5 +++++ 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index e940f8e..de83e83 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "peroxide" -version = "0.38.1" +version = "0.39.0" authors = ["axect "] 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 } diff --git a/RELEASES.md b/RELEASES.md index bf3a848..8e15eb9 100644 --- a/RELEASES.md +++ b/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