Merge pull request #92 from Axect/release

Ver 0.40.1

Author: Axect

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "peroxide"
-version = "0.40.0"
+version = "0.40.1"
 authors = ["axect <axect@outlook.kr>"]
 edition = "2018"
 description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"
@@ -39,9 +39,8 @@ peroxide-num = "0.1"
 anyhow = "1.0"
 paste = "1.0"
 netcdf = { version = "0.7", optional = true, default-features = false }
-pyo3 = { version = "0.22", optional = true, features = [
+pyo3 = { version = "0.27.1", optional = true, features = [
   "auto-initialize",
-  "gil-refs",
 ] }
 blas = { version = "0.22", optional = true }
 lapack = { version = "0.19", optional = true }

README.md

@@ -321,7 +321,7 @@ How's that? Let me know if there's anything else you'd like me to improve!
 
 ## Latest README version
 
-Corresponding to `0.38.0`
+Corresponding to `0.40.1`
 
 ## Pre-requisite
 
@@ -351,6 +351,7 @@ cargo add peroxide --features "<FEATURES>"
 * `csv`: Adds CSV support for DataFrame
 * `parquet`: Adds Parquet support for DataFrame
 * `serde`: Enables serialization/deserialization for Matrix and polynomial
+* `rkyv`: Enables zero-copy serialization/deserialization with [rkyv](https://rkyv.org)
 
 ### Install Examples
 

RELEASES.md

@@ -1,3 +1,14 @@
+# Release 0.40.1 (2026-02-06)
+
+## Bug Fixes & Improvements
+
+- Fix numerical instability in RREF (Reduced Row Echelon Form) by comparing to epsilon instead of zero ([#90](https://github.com/Axect/Peroxide/pull/90)) (Thanks to [@developing-human](https://github.com/developing-human))
+- Update `pyo3` dependency to 0.27.1 for `plot` feature compatibility ([#89](https://github.com/Axect/Peroxide/pull/89)) (Thanks to [@JSorngard](https://github.com/JSorngard))
+- Fix adaptive step size control exponent for embedded Runge-Kutta methods
+  - Add `order()` method to `ButcherTableau` trait for correct exponent `1/(p+1)`
+  - BS23: `1/3`, RKF45/DP45/TSIT45: `1/5`, RKF78: `1/8`
+- Fix misleading comments on RKF78 BU/BE coefficients
+
 # Release 0.40.0 (2025-07-24)
 
 ## Move from `arrow2` to `arrow` & `parquet`

src/numerical/ode.rs

@@ -290,6 +290,12 @@ pub trait ButcherTableau {
     fn max_step_iter(&self) -> usize {
         unimplemented!()
     }
+
+    /// Order of the lower-order solution for adaptive step size control.
+    /// The exponent `1/(order+1)` is used in the step size formula.
+    fn order(&self) -> usize {
+        4
+    }
 }
 
 impl<BU: ButcherTableau> ODEIntegrator for BU {
@@ -324,7 +330,7 @@ impl<BU: ButcherTableau> ODEIntegrator for BU {
                     error = error.max(dt * s.abs())
                 }
 
-                let factor = (self.tol() / error).powf(0.2);
+                let factor = (self.tol() / error).powf(1.0 / (self.order() as f64 + 1.0));
                 let new_dt = self.safety_factor() * dt * factor;
                 let new_dt = new_dt.clamp(self.min_step_size(), self.max_step_size());
 
@@ -558,6 +564,9 @@ impl ButcherTableau for BS23 {
     fn max_step_iter(&self) -> usize {
         self.max_step_iter
     }
+    fn order(&self) -> usize {
+        2
+    }
 }
 
 /// Runge-Kutta-Fehlberg 4/5th order integrator.
@@ -1098,9 +1107,7 @@ impl ButcherTableau for RKF78 {
         ],
     ];
 
-    // Coefficients for the 7th order solution (propagated solution)
-    // BU_i = BE_i (8th order) - ErrorCoeff_i
-    // ErrorCoeff_i = [-41/840, 0, ..., 0, -41/840 (for k11), 41/840 (for k12), 41/840 (for k13)]
+    // Coefficients for the 8th order solution (propagated via local extrapolation)
     const BU: &'static [f64] = &[
         0.0,
         0.0,
@@ -1117,8 +1124,8 @@ impl ButcherTableau for RKF78 {
         41.0 / 840.0,
     ];
 
-    // Coefficients for the 8th order solution (used for error estimation)
-    // These are from the y[i+1] formula in the MATLAB description
+    // Synthetic coefficients for error estimation
+    // BU - BE yields the Fehlberg error formula: 41/840 * (k1 + k11 - k12 - k13)
     const BE: &'static [f64] = &[
         41.0 / 840.0,
         0.0,
@@ -1150,6 +1157,9 @@ impl ButcherTableau for RKF78 {
     fn max_step_iter(&self) -> usize {
         self.max_step_iter
     }
+    fn order(&self) -> usize {
+        7
+    }
 }
 
 // ┌─────────────────────────────────────────────────────────┐

src/structure/matrix.rs

@@ -3293,14 +3293,17 @@ impl LinearAlgebra<Matrix> for Matrix {
     ///
     /// Implementation of [RosettaCode](https://rosettacode.org/wiki/Reduced_row_echelon_form)
     fn rref(&self) -> Matrix {
+        let max_abs = self.data.iter().fold(0f64, |acc, &x| acc.max(x.abs()));
+        let epsilon = (max_abs * 1e-12).max(1e-15);
+
         let mut lead = 0usize;
         let mut result = self.clone();
         'outer: for r in 0..self.row {
             if self.col <= lead {
                 break;
             }
             let mut i = r;
-            while result[(i, lead)] == 0f64 {
+            while result[(i, lead)].abs() < epsilon {
                 i += 1;
                 if self.row == i {
                     i = r;
@@ -3314,7 +3317,7 @@ impl LinearAlgebra<Matrix> for Matrix {
                 result.swap(i, r, Row);
             }
             let tmp = result[(r, lead)];
-            if tmp != 0f64 {
+            if tmp.abs() > epsilon {
                 unsafe {
                     result.row_mut(r).iter_mut().for_each(|t| *(*t) /= tmp);
                 }

src/util/plot.rs

@@ -77,12 +77,14 @@
 //! - `savefig` : Save plot with given path
 
 extern crate pyo3;
-use self::pyo3::types::IntoPyDict;
+use self::pyo3::types::{IntoPyDict, PyDictMethods};
 use self::pyo3::{PyResult, Python};
 pub use self::Grid::{Off, On};
 use self::PlotOptions::{Domain, Images, Pairs, Path};
 use std::collections::HashMap;
 use std::fmt::Display;
+use std::borrow::BorrowMut;
+use std::ffi::CString;
 
 type Vector = Vec<f64>;
 
@@ -463,7 +465,7 @@ impl Plot for Plot2D {
         }
 
         // Plot
-        Python::with_gil(|py| {
+        Python::attach(|py| {
             // Input data
             let x = self.domain.clone();
             let ys = self.images.clone();
@@ -502,24 +504,24 @@ impl Plot for Plot2D {
             let plot_type = self.plot_type.clone();
 
             // Global variables to plot
-            let globals =
-                vec![("plt", py.import_bound("matplotlib.pyplot")?)].into_py_dict_bound(py);
-            globals.as_gil_ref().set_item("x", x)?;
-            globals.as_gil_ref().set_item("y", ys)?;
-            globals.as_gil_ref().set_item("pair", pairs)?;
-            globals.as_gil_ref().set_item("n", y_length)?;
-            globals.as_gil_ref().set_item("p", pair_length)?;
+            let mut globals =
+                vec![("plt", py.import("matplotlib.pyplot")?)].into_py_dict(py)?;
+            globals.borrow_mut().set_item("x", x)?;
+            globals.borrow_mut().set_item("y", ys)?;
+            globals.borrow_mut().set_item("pair", pairs)?;
+            globals.borrow_mut().set_item("n", y_length)?;
+            globals.borrow_mut().set_item("p", pair_length)?;
             if let Some(fs) = fig_size {
-                globals.as_gil_ref().set_item("fs", fs)?;
+                globals.borrow_mut().set_item("fs", fs)?;
             }
-            globals.as_gil_ref().set_item("dp", dpi)?;
-            globals.as_gil_ref().set_item("gr", grid)?;
-            globals.as_gil_ref().set_item("pa", path)?;
+            globals.borrow_mut().set_item("dp", dpi)?;
+            globals.borrow_mut().set_item("gr", grid)?;
+            globals.borrow_mut().set_item("pa", path)?;
             if let Some(xl) = self.xlim {
-                globals.as_gil_ref().set_item("xl", xl)?;
+                globals.borrow_mut().set_item("xl", xl)?;
             }
             if let Some(yl) = self.ylim {
-                globals.as_gil_ref().set_item("yl", yl)?;
+                globals.borrow_mut().set_item("yl", yl)?;
             }
 
             // Plot Code
@@ -702,7 +704,7 @@ impl Plot for Plot2D {
                 plot_string.push_str(&format!("plt.savefig(pa, dpi={})", dpi)[..]);
             }
 
-            py.run_bound(&plot_string[..], Some(&globals), None)?;
+            py.run(&CString::new(plot_string)?, Some(&globals), None)?;
             Ok(())
         })
     }

tests/matrix.rs

@@ -111,3 +111,54 @@ fn test_kronecker() {
     let c1 = a1.kronecker(&b1);
     assert_eq!(c1, ml_matrix("0 5 0 10;6 7 12 14;0 15 0 20;18 21 24 28"));
 }
+
+#[test]
+fn test_rref() {
+    let a = ml_matrix(
+        r#"
+        -3 2 -1 -1;
+        6 -6 7 -7;
+        3 -4 4 -6"#,
+    );
+    let b = a.rref();
+
+    assert_eq!(
+        b,
+        ml_matrix(
+            r#"
+    1 0 0 2;
+    0 1 0 2;
+    0 0 1 -1"#
+        )
+    );
+}
+
+#[test]
+fn test_rref_unstable() {
+    let epsilon = 1e-10;
+
+    // this matrix used to become unstable during rref
+    let a = ml_matrix(
+        r#"
+    1 1 0 0 0 1 0 1 31;
+    1 1 1 1 0 0 1 1 185;
+    0 0 1 0 0 1 1 1 165;
+    1 0 1 0 1 1 0 1 32;
+    1 0 1 0 0 0 1 1 174;
+    0 0 1 0 1 1 1 1 171;
+    0 1 1 0 1 1 0 1 27;
+    1 0 0 1 0 1 0 0 20;
+    1 0 1 1 0 1 0 0 23"#,
+    );
+
+    let b = a.rref();
+
+    // creating a row like "0 0 0 0 0 0 0 0 1" will "prove" 0 == 1
+    // which is a tell of numeric instability
+    for row in 0..b.row {
+        let ends_in_1 = (b[(row, b.col - 1)] - 1.0).abs() < epsilon;
+        let rest_zeroes = (0..b.col - 1).all(|col| b[(row, col)].abs() < epsilon);
+
+        assert!(!(ends_in_1 && rest_zeroes));
+    }
+}