From b9710ed47be2c797bc3734cffc3ed36425aea032 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= Date: Fri, 31 Oct 2025 21:06:43 +0100 Subject: [PATCH 1/7] Update the `pyo3` dependency to version 0.27.1. This way `peroxied` compiles in CI on the latest images. --- Cargo.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index c480002..4114419 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -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 } From 76bb265ec1b86b78a2dc01e13f225fa4f6db053b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= Date: Fri, 31 Oct 2025 21:34:00 +0100 Subject: [PATCH 2/7] Update `plot.rs` to use the latest version of `pyo3` --- src/util/plot.rs | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/src/util/plot.rs b/src/util/plot.rs index 978ea7d..377a293 100644 --- a/src/util/plot.rs +++ b/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; @@ -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(()) }) } From b5d298901373da84b4fda21a6629b90997756b14 Mon Sep 17 00:00:00 2001 From: Developing Human Date: Fri, 19 Dec 2025 22:13:13 -0500 Subject: [PATCH 3/7] FIX: compare to epsilon during rref --- src/structure/matrix.rs | 5 ++-- tests/matrix.rs | 51 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 2 deletions(-) diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs index 72220e2..be93717 100644 --- a/src/structure/matrix.rs +++ b/src/structure/matrix.rs @@ -3293,6 +3293,7 @@ impl LinearAlgebra for Matrix { /// /// Implementation of [RosettaCode](https://rosettacode.org/wiki/Reduced_row_echelon_form) fn rref(&self) -> Matrix { + let epsilon = 1e-10; let mut lead = 0usize; let mut result = self.clone(); 'outer: for r in 0..self.row { @@ -3300,7 +3301,7 @@ impl LinearAlgebra for Matrix { 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 +3315,7 @@ impl LinearAlgebra 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); } diff --git a/tests/matrix.rs b/tests/matrix.rs index defb914..d3a2110 100644 --- a/tests/matrix.rs +++ b/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 has a tendency 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)); + } +} From e4518f4e0332b746a27b0bb1e63875526fb00b0f Mon Sep 17 00:00:00 2001 From: Developing Human Date: Fri, 19 Dec 2025 22:26:11 -0500 Subject: [PATCH 4/7] cleanup comment --- tests/matrix.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/matrix.rs b/tests/matrix.rs index d3a2110..3ce5d3c 100644 --- a/tests/matrix.rs +++ b/tests/matrix.rs @@ -137,7 +137,7 @@ fn test_rref() { fn test_rref_unstable() { let epsilon = 1e-10; - // this matrix has a tendency to become unstable during rref, + // this matrix used to become unstable during rref let a = ml_matrix( r#" 1 1 0 0 0 1 0 1 31; From f987e345a1808ecebd5514addb15f475848affbc Mon Sep 17 00:00:00 2001 From: Developing Human Date: Sat, 20 Dec 2025 12:43:04 -0500 Subject: [PATCH 5/7] compute epsilon based on the scale of matrix data --- src/structure/matrix.rs | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs index be93717..7e0241d 100644 --- a/src/structure/matrix.rs +++ b/src/structure/matrix.rs @@ -3293,7 +3293,9 @@ impl LinearAlgebra for Matrix { /// /// Implementation of [RosettaCode](https://rosettacode.org/wiki/Reduced_row_echelon_form) fn rref(&self) -> Matrix { - let epsilon = 1e-10; + 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 { From 36c7089fd5d52ac05110cdc5d8e72e054e2430cc Mon Sep 17 00:00:00 2001 From: Axect Date: Fri, 6 Feb 2026 19:00:00 +0800 Subject: [PATCH 6/7] FIX: Correct adaptive step size exponent for embedded RK --- src/numerical/ode.rs | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/src/numerical/ode.rs b/src/numerical/ode.rs index 079aa15..5834a31 100644 --- a/src/numerical/ode.rs +++ b/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 ODEIntegrator for BU { @@ -324,7 +330,7 @@ impl 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 + } } // ┌─────────────────────────────────────────────────────────┐ From 0692f998a95a7773e0812c36f6d6b7e105d68d36 Mon Sep 17 00:00:00 2001 From: Axect Date: Fri, 6 Feb 2026 19:06:46 +0800 Subject: [PATCH 7/7] RLSE: Ver 0.40.1 - 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 --- Cargo.toml | 2 +- README.md | 3 ++- RELEASES.md | 11 +++++++++++ 3 files changed, 14 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 4114419..786aab8 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "peroxide" -version = "0.40.0" +version = "0.40.1" authors = ["axect "] edition = "2018" description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax" diff --git a/README.md b/README.md index 8dcb78d..bfa44fc 100644 --- a/README.md +++ b/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 "" * `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 diff --git a/RELEASES.md b/RELEASES.md index b36fa07..a21d354 100644 --- a/RELEASES.md +++ b/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`