Ver 0.37.3

Axect · Axect · commit d6665b7d7697 · 2024-05-01T18:14:16.000+09:00
- Add Nan/infinite guard to `gauss_kronrod_quadrature` (early exit) (#59) - Add complex feature & complex module (#35) - Implement Cubic B-Spline basis function
diff --git a/Cargo.toml b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "peroxide"
-version = "0.37.2"
+version = "0.37.3"
 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"
@@ -38,6 +38,7 @@ 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 }
+num-complex = { version = "0.4", optional = true }
 
 [package.metadata.docs.rs]
 rustdoc-args = [ "--html-in-header", "katex-header.html", "--cfg", "docsrs"]
@@ -48,3 +49,4 @@ O3 = ["blas", "lapack"]
 plot = ["pyo3"]
 nc = ["netcdf"]
 parquet = ["arrow2"]
+complex = ["num-complex", "matrixmultiply/cgemm"]
diff --git a/README.md b/README.md
@@ -173,6 +173,8 @@ Peroxide can do many things.
     - Cubic Hermite Spline
       - Estimate slope via Akima
       - Estimate slope via Quadratic interpolation
+    - B-Spline (incomplete)
+      - Uniform Cubic B-Spline basis function
   - Non-linear regression
     - Gradient Descent
     - Levenberg Marquardt
diff --git a/RELEASES.md b/RELEASES.md
@@ -1,3 +1,11 @@
+# Release 0.37.3 (2024-05-01)
+
+- Add Nan/infinite guard to `gauss_kronrod_quadrature` (early exit) ([#59](https://github.com/Axect/Peroxide/pull/59)) (Thanks to [@GComitini](https://github.com/GComitini))
+- Add complex feature & complex module ([#35](https://github.com/Axect/Peroxide/issues/35))
+- Implement Cubic B-Spline basis functions
+  - `UnitCubicBasis`
+  - `CubicBSplineBases`
+
 # Release 0.37.2 (2024-04-16)
 
 - Do not include legend box if there is no legend ([#58](https://github.com/Axect/Peroxide/pull/58)) (Thanks to [@GComitini](https://github.com/GComitini))
diff --git a/example_data/cubic_b_spline.png b/example_data/cubic_b_spline.png
diff --git a/example_data/rkf45_test.png b/example_data/rkf45_test.png
diff --git a/example_data/test_plot.png b/example_data/test_plot.png
diff --git a/src/complex/mod.rs b/src/complex/mod.rs
@@ -0,0 +1,5 @@
+use num_complex::Complex;
+
+pub type C64 = Complex<f64>;
+
+pub mod vector;
diff --git a/src/complex/vector.rs b/src/complex/vector.rs
@@ -3,16 +3,14 @@ use crate::traits::fp::FPVector;
 use crate::traits::math::{Vector, Normed, Norm, InnerProduct};
 use crate::traits::sugar::VecOps;
 
-pub type C64 = Complex<f64>;
-
 impl Vector for Complex<f64> {
     type Scalar = Self;
 
-    fn add_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self {
+    fn add_vec(&self, rhs: &Self) -> Self {
         self + rhs
     }
 
-    fn sub_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self {
+    fn sub_vec(&self, rhs: &Self) -> Self {
         self - rhs
     }
 
@@ -70,15 +68,15 @@ impl FPVector for Vec<Complex<f64>> {
     fn filter<F>(&self, f: F) -> Self
     where
             F: Fn(Self::Scalar) -> bool {
-        self.into_iter().filter(|&x| f(*x)).map(|&t| t).collect()
+        self.iter().filter(|&x| f(*x)).cloned().collect()
     }
 
     fn take(&self, n: usize) -> Self {
-        self.iter().take(n).map(|&x| x).collect()
+        self.iter().take(n).cloned().collect()
     }
 
     fn skip(&self, n: usize) -> Self {
-        self.iter().skip(n).map(|&x| x).collect()
+        self.iter().skip(n).cloned().collect()
     }
 
     fn sum(&self) -> Self::Scalar {
@@ -93,11 +91,11 @@ impl FPVector for Vec<Complex<f64>> {
 impl Vector for Vec<Complex<f64>> {
     type Scalar = Complex<f64>;
 
-    fn add_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self {
+    fn add_vec(&self, rhs: &Self) -> Self {
         self.zip_with(|x, y| x + y, rhs)
     }
 
-    fn sub_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self {
+    fn sub_vec(&self, rhs: &Self) -> Self {
         self.zip_with(|x, y| x - y, rhs)
     }
 
diff --git a/src/lib.rs b/src/lib.rs
@@ -205,3 +205,6 @@ pub mod statistics;
 pub mod structure;
 pub mod traits;
 pub mod util;
+
+#[cfg(feature = "complex")]
+pub mod complex;
diff --git a/src/numerical/integral.rs b/src/numerical/integral.rs
@@ -240,6 +240,9 @@ where
                     tol
                 };
                 if (G - K).abs() < tol_curr || a == b || max_iter == 0 {
+                    if ! G.is_finite() {
+                        return G;
+                    }
                     I += G;
                 } else {
                     S.push((a, c, tol / 2f64, max_iter - 1));
diff --git a/src/numerical/spline.rs b/src/numerical/spline.rs
@@ -194,9 +194,15 @@
 //! }
 //! ```
 //!
+//! # B-Spline (incomplete)
+//!
+//! - `UnitCubicBasis`: Single cubic B-Spline basis function
+//! - `CubicBSplineBases`: Uniform Cubic B-Spline basis functions
+//!
 //! # References
 //!
-//! * Gary D. Knott, *Interpolating Splines*, Birkhäuser Boston, MA, (2000).
+//! - Gary D. Knott, *Interpolating Splines*, Birkhäuser Boston, MA, (2000).
+/// - [Wikipedia - Irwin-Hall distribution](https://en.wikipedia.org/wiki/Irwin%E2%80%93Hall_distribution#Special_cases)
 
 use self::SplineError::{NotEnoughNodes, NotEqualNodes, NotEqualSlopes, RedundantNodeX};
 #[allow(unused_imports)]
@@ -843,3 +849,141 @@ fn quadratic_slopes(x: &[f64], y: &[f64]) -> Result<Vec<f64>> {
 
     Ok(m)
 }
+
+// =============================================================================
+// B-Spline
+// =============================================================================
+/// Unit Cubic Basis Function
+///
+/// # Description
+/// Unit cubic basis function from Irwin-Hall distribution (n=4).
+/// For general interval, we substitute t = 4 * (x - a) / (b - a).
+///
+/// # Reference
+/// [Wikipedia](https://en.wikipedia.org/wiki/Irwin%E2%80%93Hall_distribution#Special_cases)
+#[derive(Debug, Copy, Clone)]
+#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
+pub struct UnitCubicBasis {
+    pub x_min: f64,
+    pub x_max: f64,
+    pub scale: f64,
+}
+
+impl UnitCubicBasis {
+    pub fn new(x_min: f64, x_max: f64, scale: f64) -> Self {
+        Self { x_min, x_max, scale }
+    }
+
+    pub fn eval(&self, x: f64) -> f64 {
+        let t = 4f64 * (x - self.x_min) / (self.x_max - self.x_min);
+
+        let result = if (0f64..1f64).contains(&t) {
+            t.powi(3) / 6f64
+        } else if (1f64..2f64).contains(&t) {
+            (-3f64 * t.powi(3) + 12f64 * t.powi(2) - 12f64 * t + 4f64) / 6f64
+        } else if (2f64..3f64).contains(&t) {
+            (3f64 * t.powi(3) - 24f64 * t.powi(2) + 60f64 * t - 44f64) / 6f64
+        } else if (3f64..4f64).contains(&t) {
+            (4f64 - t).powi(3) / 6f64
+        } else {
+            0f64
+        };
+
+        self.scale * result
+    }
+
+    pub fn eval_vec(&self, x: &[f64]) -> Vec<f64> {
+        x.iter().map(|x| self.eval(*x)).collect()
+    }
+}
+
+/// Uniform Cubic B-Spline basis functions
+///
+/// # Example
+///
+/// ```rust
+/// use peroxide::fuga::*;
+/// use core::ops::Range;
+/// 
+/// # #[allow(unused_variables)]
+/// fn main() -> anyhow::Result<()> {
+///     let cubic_b_spline = CubicBSplineBases::from_interval((0f64, 1f64), 5);
+///     let x = linspace(0f64, 1f64, 1000);
+///     let y = cubic_b_spline.eval_vec(&x);
+/// 
+/// #   #[cfg(feature = "plot")] {
+///     let mut plt = Plot2D::new();
+///     plt.set_domain(x.clone());
+/// 
+///     for basis in &cubic_b_spline.bases {
+///         plt.insert_image(basis.eval_vec(&x));
+///     }
+///
+///     plt
+///         .insert_image(y)
+///         .set_xlabel(r"$x$")
+///         .set_ylabel(r"$y$")
+///         .set_style(PlotStyle::Nature)
+///         .tight_layout()
+///         .set_dpi(600)
+///         .set_path("example_data/cubic_b_spline.png")
+///         .savefig()?;
+/// #   }
+///     Ok(())
+/// }
+/// ```
+#[derive(Debug, Clone)]
+#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
+pub struct CubicBSplineBases {
+    pub ranges: Vec<Range<f64>>,
+    pub bases: Vec<UnitCubicBasis>,
+}
+
+impl CubicBSplineBases {
+    /// Create new Cubic B-Spline basis functions
+    pub fn new(ranges: Vec<Range<f64>>, bases: Vec<UnitCubicBasis>) -> Self {
+        Self { ranges, bases }
+    }
+
+    /// Create new Cubic B-Spline basis functions for `[a, b]`
+    pub fn from_interval((a, b): (f64, f64), num_bases: usize) -> Self {
+        let nodes = linspace(a, b, num_bases + 4);
+        let (ranges, bases) = nodes
+            .iter()
+            .zip(nodes.iter().skip(4))
+            .map(|(a, b)| (Range { start: *a, end: *b }, UnitCubicBasis::new(*a, *b, 1f64)))
+            .unzip();
+        
+        Self::new(ranges, bases)
+    }
+
+    /// Rescale all basis functions
+    ///
+    /// # Arguments
+    /// - `scale_vec` - scale vector
+    pub fn rescale(&mut self, scale_vec: &[f64]) -> Result<()> {
+        if scale_vec.len() != self.bases.len() {
+            bail!("The number of scales should be equal to the number of basis functions");
+        }
+
+        for (basis, scale) in self.bases.iter_mut().zip(scale_vec) {
+            basis.scale = *scale;
+        }
+
+        Ok(())
+    }
+
+    pub fn eval(&self, x: f64) -> f64 {
+        self.ranges.iter()
+            .enumerate()
+            .filter(|(_, range)| range.contains(&x))
+            .fold(0f64, |acc, (i, _)| {
+                let basis = &self.bases[i];
+                acc + basis.eval(x)
+            })
+    }
+
+    pub fn eval_vec(&self, x: &[f64]) -> Vec<f64> {
+        x.iter().map(|x| self.eval(*x)).collect()
+    }
+}
diff --git a/src/structure/vector.rs b/src/structure/vector.rs
