Ver 0.36.1

- Generic Butcher Tableau
- Fix all warnings
- More & safe non-macro

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "peroxide"
-version = "0.36.0"
+version = "0.36.1"
 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"

RELEASES.md

@@ -1,3 +1,39 @@
+# Release 0.36.1 (2024-04-09)
+
+- Fix all warnings in peroxide
+- Change redundant method
+  - `Vec<f64>::resize` -> `Vec<f64>::reshape`
+- Error handling for concatenation
+  - `cbind` & `rbind` now returns `Result<Matrix, ConcatenateError>`
+- New non-macro utils
+  - `column_stack(&[Vec<f64>]) -> Result<Matrix, ConcatenateError>`
+  - `row_stack(&[Vec<f64>]) -> Result<Matrix, ConcatenateError>`
+  - `rand_with_rng(usize, usize, &mut Rng) -> Matrix`
+- Generic Butcher tableau trait (now for embedded Runge-Kutta methods)
+
+  ```rust
+  pub trait ButcherTableau {
+      const C: &'static [f64];
+      const A: &'static [&'static [f64]];
+      const BH: &'static [f64];
+      const BL: &'static [f64];
+
+      fn tol(&self) -> f64;
+      fn safety_factor(&self) -> f64;
+      fn max_step_size(&self) -> f64;
+      fn min_step_size(&self) -> f64;
+      fn max_step_iter(&self) -> usize;
+  }
+  ```
+
+  - Implement `ODEIntegrator` for `ButcherTableau`
+    - Just declare `ButcherTableau` then `step` is free
+
+  - Three available embedded Runge-Kutta methods
+    - `RKF45`: Runge-Kutta-Fehlberg 4/5th order
+    - `DP45`: Dormand-Prince 4/5th order
+    - `TSIT45`: Tsitouras 4/5th order
+
 # Release 0.36.0 (2024-04-08)
 
 ## Huge Update - Error handling & Whole new ODE

examples/lorenz_dp45.rs

@@ -0,0 +1,47 @@
+use peroxide::fuga::*;
+
+#[allow(unused_variables)]
+fn main() -> Result<(), Box<dyn Error>> {
+    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_mat = py_matrix(y_vec);
+    let y0 = y_mat.col(0);
+    let y2 = y_mat.col(2);
+
+    #[cfg(feature = "plot")]
+    {
+        let mut plt = Plot2D::new();
+        plt
+            .set_domain(y0)
+            .insert_image(y2)
+            .set_xlabel(r"$y_0$")
+            .set_ylabel(r"$y_2$")
+            .set_style(PlotStyle::Nature)
+            .tight_layout()
+            .set_dpi(600)
+            .set_path("example_data/lorenz_dp45.png")
+            .savefig()?;
+    }
+
+    Ok(())
+}
+
+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]) -> Result<(), ODEError> {
+        dy[0] = 10f64 * (y[1] - y[0]);
+        dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];
+        dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1];
+        Ok(())
+    }
+}

examples/lorenz_gl4.rs

@@ -1,7 +1,8 @@
 use peroxide::fuga::*;
 
+#[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
-    let gl4 = GL4::new(ImplicitSolver::FixedPoint, 100, 1e-6);
+    let gl4 = GL4::new(ImplicitSolver::FixedPoint, 1e-6, 100);
     let basic_ode_solver = BasicODESolver::new(gl4);
     let (_, y_vec) = basic_ode_solver.solve(
         &Lorenz,

examples/lorenz_rk4.rs

@@ -1,5 +1,6 @@
 use peroxide::fuga::*;
 
+#[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
     let basic_ode_solver = BasicODESolver::new(RK4);
     let (_, y_vec) = basic_ode_solver.solve(
@@ -36,7 +37,7 @@ impl ODEProblem for Lorenz {
         vec![10f64, 1f64, 1f64]
     }
 
-    fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> {
+    fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> {
         dy[0] = 10f64 * (y[1] - y[0]);
         dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];
         dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1];

examples/lorenz_rkf45.rs

@@ -1,5 +1,6 @@
 use peroxide::fuga::*;
 
+#[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
     let rkf45 = RKF45::new(1e-4, 0.9, 1e-6, 1e-1, 100);
     let basic_ode_solver = BasicODESolver::new(rkf45);
@@ -37,7 +38,7 @@ impl ODEProblem for Lorenz {
         vec![10f64, 1f64, 1f64]
     }
 
-    fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> {
+    fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> {
         dy[0] = 10f64 * (y[1] - y[0]);
         dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];
         dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1];

examples/lorenz_tsit45.rs

@@ -0,0 +1,47 @@
+use peroxide::fuga::*;
+
+#[allow(unused_variables)]
+fn main() -> Result<(), Box<dyn Error>> {
+    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_mat = py_matrix(y_vec);
+    let y0 = y_mat.col(0);
+    let y2 = y_mat.col(2);
+
+    #[cfg(feature = "plot")]
+    {
+        let mut plt = Plot2D::new();
+        plt
+            .set_domain(y0)
+            .insert_image(y2)
+            .set_xlabel(r"$y_0$")
+            .set_ylabel(r"$y_2$")
+            .set_style(PlotStyle::Nature)
+            .tight_layout()
+            .set_dpi(600)
+            .set_path("example_data/lorenz_tsit45.png")
+            .savefig()?;
+    }
+
+    Ok(())
+}
+
+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]) -> Result<(), ODEError> {
+        dy[0] = 10f64 * (y[1] - y[0]);
+        dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];
+        dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1];
+        Ok(())
+    }
+}

examples/ode_env.rs

@@ -1,5 +1,6 @@
 use peroxide::fuga::*;
 
+#[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
     let t = seq(0, 10, 1);
     let y = t

examples/ode_test_gl4.rs

@@ -1,9 +1,8 @@
-#[macro_use]
-extern crate peroxide;
 use peroxide::fuga::*;
 
+#[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
-    let gl4 = GL4::new(ImplicitSolver::FixedPoint, 100, 1e-6);
+    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,

examples/ode_test_rk4.rs

@@ -1,5 +1,6 @@
 use peroxide::fuga::*;
 
+#[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
     let basic_ode_solver = BasicODESolver::new(RK4);
     let (t_vec, y_vec) = basic_ode_solver.solve(

examples/ode_test_rkf45.rs

@@ -1,5 +1,6 @@
 use peroxide::fuga::*;
 
+#[allow(unused_variables)]
 fn main() -> Result<(), Box<dyn Error>> {
     let rkf = RKF45::new(1e-4, 0.9, 1e-6, 1e-1, 100);
     let basic_ode_solver = BasicODESolver::new(rkf);

src/lib.rs

@@ -147,14 +147,17 @@
 //!
 //! ## Useful tips for features
 //!
-//! * After `0.28.0`, `dataframe` feature is replaced by `nc` feature.
-//! * If you want to use *QR* or *SVD* or *Cholesky Decomposition* then should use `O3` feature (there are no implementations for these decompositions in `default`)
-//! * If you want to write your numerical results, then use `parquet` or `nc` features (corresponding to `parquet` or `netcdf` format. (It is much more effective than `csv` and `json`.)
-//! * After `0.23.0`, there are two options - `fuga`, `prelude`. Choose proper option for your computations.
-//! * To plot, use `parquet` or `nc` feature to export data as parquet or netcdf format and use python to draw plot.
-//!     * `plot` feature has limited plot abilities.
-//!     * To read parquet file in python, use `pandas` & `pyarrow` libraries.
-//!     * There is a template of python code for netcdf. - [Socialst](https://github.com/Axect/Socialst/blob/master/Templates/PyPlot_Template/nc_plot.py)
+//! * If you want to use _QR_, _SVD_, or _Cholesky Decomposition_, you should use the `O3` feature. These decompositions are not implemented in the `default` feature.
+//!
+//! * If you want to save your numerical results, consider using the `parquet` or `nc` features, which correspond to the `parquet` and `netcdf` file formats, respectively. These formats are much more efficient than `csv` and `json`.
+//!
+//! * For plotting, it is recommended to use the `plot` feature. However, if you require more customization, you can use the `parquet` or `nc` feature to export your data in the parquet or netcdf format and then use Python to create the plots.
+//!
+//!     * To read parquet files in Python, you can use the `pandas` and `pyarrow` libraries.
+//!
+//!     * A template for Python code that works with netcdf files can be found in the [Socialst](https://github.com/Axect/Socialst/blob/master/Templates/PyPlot_Template/nc_plot.py) repository.
+
+//!
 #![cfg_attr(docsrs, feature(doc_auto_cfg))]
 
 #[cfg(feature = "O3")]

src/numerical/eigen.rs

@@ -2,13 +2,9 @@
 //!
 //! * Reference : Press, William H., and William T. Vetterling. *Numerical Recipes.* Cambridge: Cambridge Univ. Press, 2007.
 
-#[cfg(feature = "O3")]
-use lapack::{dorgtr, dsteqr, dsytrd};
-
 pub use self::EigenMethod::*;
 use crate::structure::matrix::Matrix;
 use crate::util::non_macro::eye_shape;
-use std::f64::EPSILON;
 
 #[derive(Debug, Copy, Clone)]
 pub enum EigenMethod {
@@ -114,11 +110,11 @@ impl JacobiTemp {
             for ip in 0..n - 1 {
                 for iq in ip + 1..n {
                     let g = 100f64 * a[(ip, iq)].abs();
-                    if i > 4 && g <= EPSILON * d[ip].abs() && g <= EPSILON * d[iq].abs() {
+                    if i > 4 && g <= f64::EPSILON * d[ip].abs() && g <= f64::EPSILON * d[iq].abs() {
                         a[(ip, iq)] = 0f64;
                     } else if a[(ip, iq)].abs() > tresh {
                         h = d[iq] - d[ip];
-                        let t = if g <= EPSILON * h.abs() {
+                        let t = if g <= f64::EPSILON * h.abs() {
                             a[(ip, iq)] / h
                         } else {
                             let theta = 0.5 * h / a[(ip, iq)];

src/numerical/mod.rs

@@ -1,7 +1,5 @@
 //! Differential equations & Numerical Analysis tools
 
-//pub mod bdf;
-//pub mod gauss_legendre;
 pub mod eigen;
 pub mod integral;
 pub mod interp;

src/numerical/newton.rs

@@ -12,12 +12,12 @@ pub fn newton<F: Fn(&Vec<AD>) -> Vec<AD> + Copy>(init_cond: Vec<f64>, f: F, rtol
     let mut x_next = init_cond;
     let mut x = x_next.clone();
     update(&mut x_next, f);
-    let mut err = (x_next.sub_vec(&x)).norm(Norm::L2);
+    let mut err = x_next.sub_vec(&x).norm(Norm::L2);
 
     while err >= rtol {
         x = x_next.clone();
         update(&mut x_next, f);
-        err = (x_next.sub_vec(&x)).norm(Norm::L2);
+        err = x_next.sub_vec(&x).norm(Norm::L2);
     }
 
     x_next