diff --git a/Cargo.toml b/Cargo.toml index 60be93ec..50b4a0e1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "peroxide" -version = "0.37.5" +version = "0.37.6" 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/RELEASES.md b/RELEASES.md index 8aec81ed..f918b169 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -1,3 +1,18 @@ +# Release 0.37.6 (2024-06-19) + +## Huge Spline Change + +- Generic Spline trait + - `Spline`: desired output type is `T` +- Split `PolynomialSpline` from `Spline` + - `CubicSpline` & `CubicHermiteSpline` are now `PolynomialSpline` + - Implement `Spline` for `PolynomialSpline` +- Implement B-Spline + - `BSpline { degree: usize, knots: Vec, control_points: Vec> }` + - `BSpline::open(degree, knots, control_points)` : Open B-Spline + - `BSpline::clamped(degree, knots, control_points)` : Clamped B-Spline +- Implement `Spline<(f64, f64)>` for `BSpline` + # Release 0.37.5 (2024-06-10) - More generic & stable root finding macros (except `Newton`) diff --git a/example_data/b_spline_test.png b/example_data/b_spline_test.png new file mode 100644 index 00000000..c1641f93 Binary files /dev/null and b/example_data/b_spline_test.png differ diff --git a/examples/b_spline_test.rs b/examples/b_spline_test.rs new file mode 100644 index 00000000..46cdda59 --- /dev/null +++ b/examples/b_spline_test.rs @@ -0,0 +1,46 @@ +use peroxide::fuga::*; + +fn main() -> Result<(), Box> { + let knots = vec![0f64, 1f64, 2f64, 3f64]; + let degree = 3; + let control_points = vec![ + vec![0f64, 2f64], + vec![0.2, -1f64], + vec![0.4, 1f64], + vec![0.6, -1f64], + vec![0.8, 1f64], + vec![1f64, 2f64], + ]; + + + let spline = BSpline::clamped(degree, knots, control_points.clone())?; + let t = linspace(0f64, 3f64, 200); + let (x, y): (Vec, Vec) = spline.eval_vec(&t).into_iter().unzip(); + + #[cfg(feature = "plot")] + { + let control_x = control_points.iter().map(|v| v[0]).collect::>(); + let control_y = control_points.iter().map(|v| v[1]).collect::>(); + + let mut plt = Plot2D::new(); + plt + .insert_pair((x.clone(), y.clone())) + .insert_pair((control_x.clone(), control_y.clone())) + .set_plot_type(vec![(1, PlotType::Scatter)]) + .set_color(vec![(0, "darkblue"), (1, "red")]) + .set_xlabel(r"$x$") + .set_ylabel(r"$y$") + .set_style(PlotStyle::Nature) + .set_dpi(600) + .set_path("example_data/b_spline_test.png") + .savefig()?; + } + + let mut df = DataFrame::new(vec![]); + df.push("t", Series::new(t)); + df.push("x", Series::new(x)); + df.push("y", Series::new(y)); + df.print(); + + Ok(()) +} diff --git a/src/numerical/spline.rs b/src/numerical/spline.rs index f6be15b7..c1da0df1 100644 --- a/src/numerical/spline.rs +++ b/src/numerical/spline.rs @@ -4,35 +4,48 @@ //! //! * Cubic spline //! * Cubic Hermite spline +//! * B-spline //! -//! # `Spline` trait +//! # `Spline` trait //! //! ## Methods //! -//! Let `T: Into + Copy` -//! * `fn eval(&self, x: T) -> f64` : Evaluate the spline at x -//! * `fn eval_vec(&self, v: &[T]) -> Vec` : Evaluate spline values for an array v -//! * `fn polynomial_at(&self, x: T) -> &Polynomial` : Get the polynomial at x +//! * `fn eval(&self, t: f64) -> T` : Evaluate the spline at t +//! * For Cubic Spline, t means x (domain) and `T = f64` +//! * For B-Spline, t means parameter of curve and `T = (f64, f64)` +//! * `fn eval_vec(&self, v: &[f64]) -> Vec` : Evaluate spline values for an array v +//! * `fn eval_with_cond(&self, t: f64, cond: F) -> T` : Evaluate the spline at t, with a condition +//! * `fn eval_vec_with_cond(&self, v: &[f64], cond: F) -> Vec` : Evaluate spline values for an array v, with a condition +//! +//! # `PolynomialSpline` trait +//! +//! ## Methods +//! +//! * `fn polynomial_at(&self, x: f64) -> &Polynomial` : Get the polynomial at x //! * `fn number_of_polynomials(&self) -> usize` : Get the number of polynomials //! * `fn get_ranged_polynomials(&self) -> &Vec<(Range, Polynomial)>` : Get the polynomials -//! * `fn eval_with_cond(&self, x: f64, cond: F) -> f64` : Evaluate the spline at x, with a condition -//! * `fn eval_vec_with_cond(&self, v: &[f64], cond: F) -> Vec` : Evaluate spline values for an array v, with a condition //! //! # Low-level interface //! //! ## Members //! //! * `CubicSpline`: Structure for cubic spline -//! * `fn from_nodes(node_x: &[f64], node_y: &[f64]) -> Result` : Create a cubic spline from nodes -//! * `fn extend_with_nodes(&mut self, node_x: Vec, node_y: Vec) -> Result<()>` : Extend the spline with nodes +//! * `fn from_nodes(node_x: &[f64], node_y: &[f64]) -> Result` : Create a cubic spline from nodes +//! * `fn extend_with_nodes(&mut self, node_x: Vec, node_y: Vec) -> Result<()>` : Extend the spline with nodes //! * `CubicHermiteSpline`: Structure for cubic Hermite spline -//! * `fn from_nodes_with_slopes(node_x: &[f64], node_y: &[f64], m: &[f64]) -> Result` : Create a Cubic Hermite spline from nodes with slopes -//! * `fn from_nodes(node_x: &[f64], node_y: &[f64], slope_method: SlopeMethod) -> Result` : Create a Cubic Hermite spline from nodes with slope estimation methods -//! * `SlopeMethod`: Enum for slope estimation methods +//! * `fn from_nodes_with_slopes(node_x: &[f64], node_y: &[f64], m: &[f64]) -> Result` : Create a Cubic Hermite spline from nodes with slopes +//! * `fn from_nodes(node_x: &[f64], node_y: &[f64], slope_method: SlopeMethod) -> Result` : Create a Cubic Hermite spline from nodes with slope estimation methods +//! * `SlopeMethod`: Enum for slope estimation methods //! * `Akima`: Akima's method to estimate slopes ([Akima (1970)](https://dl.acm.org/doi/abs/10.1145/321607.321609)) //! * `Quadratic`: Using quadratic interpolation to estimate slopes +//! * `BSpline`: Structure for B-Spline +//! * `fn open(degree: usize, knots: Vec, control_points: Vec>) -> Result` : Create an open B-Spline +//! * `fn clamped(degree: usize, knots: Vec, control_points: Vec>) -> Result` +//! : Create a clamped B-Spline +//! * `fn cox_de_boor(t: f64, i: f64)` : Cox-de Boor recursion formula (Here, use iteration +//! instead of recursion) //! -//! ## Usage +//! ## Usage (Cubic Spline Family) //! //! ```rust //! use peroxide::fuga::*; @@ -79,6 +92,60 @@ //! } //! ``` //! +//! ## Usage (B-Spline) +//! +//! ```rust +//! use peroxide::fuga::*; +//! +//! fn main() -> Result<(), Box> { +//! let knots = vec![0f64, 1f64, 2f64, 3f64]; +//! let degree = 3; +//! let control_points = vec![ +//! vec![0f64, 2f64], +//! vec![0.2, -1f64], +//! vec![0.4, 1f64], +//! vec![0.6, -1f64], +//! vec![0.8, 1f64], +//! vec![1f64, 2f64], +//! ]; +//! +//! +//! let spline = BSpline::clamped(degree, knots, control_points.clone())?; +//! let t = linspace(0f64, 3f64, 200); +//! let (x, y): (Vec, Vec) = spline.eval_vec(&t).into_iter().unzip(); +//! +//! # #[cfg(feature = "plot")] +//! # { +//! let control_x = control_points.iter().map(|v| v[0]).collect::>(); +//! let control_y = control_points.iter().map(|v| v[1]).collect::>(); +//! +//! let mut plt = Plot2D::new(); +//! plt +//! .insert_pair((x.clone(), y.clone())) +//! .insert_pair((control_x.clone(), control_y.clone())) +//! .set_plot_type(vec![(1, PlotType::Scatter)]) +//! .set_color(vec![(0, "darkblue"), (1, "red")]) +//! .set_xlabel(r"$x$") +//! .set_ylabel(r"$y$") +//! .set_style(PlotStyle::Nature) +//! .set_dpi(600) +//! .set_path("example_data/b_spline_test.png") +//! .savefig()?; +//! # } +//! +//! let mut df = DataFrame::new(vec![]); +//! df.push("t", Series::new(t)); +//! df.push("x", Series::new(x)); +//! df.push("y", Series::new(y)); +//! df.print(); +//! +//! Ok(()) +//! } +//! ``` +//! +//! - Result for above code is: +//! ![b_spline_test](https://github.com/Axect/Peroxide/blob/master/example_data/b_spline_test.png?raw=true) +//! //! # High-level interface //! //! ## Functions @@ -133,7 +200,7 @@ //! } //! ``` //! -//! # Calculus with splines +//! # Calculus with polynomial splines //! //! ## Usage //! @@ -194,7 +261,7 @@ //! } //! ``` //! -//! # B-Spline (incomplete) +//! # B-Spline utils //! //! - `UnitCubicBasis`: Single cubic B-Spline basis function //! - `CubicBSplineBases`: Uniform Cubic B-Spline basis functions @@ -222,29 +289,32 @@ use std::convert::From; use std::ops::{Index, Range}; use anyhow::{Result, bail}; +pub trait Spline { + fn eval(&self, t: f64) -> T; + fn eval_vec(&self, v: &[f64]) -> Vec { + v.iter().map(|&t| self.eval(t)).collect() + } + fn eval_with_cond T>(&self, t: f64, cond: F) -> T { + cond(self.eval(t)) + } + fn eval_vec_with_cond T + Copy>(&self, t: &[f64], cond: F) -> Vec { + t.iter().map(|&x| self.eval_with_cond(x, cond)).collect() + } +} + /// Trait for spline interpolation /// /// # Available Splines /// /// - `CubicSpline` /// - `CubicHermiteSpline` -pub trait Spline { - fn eval + Copy>(&self, x: T) -> f64 { - let x = x.into(); - +impl Spline for P { + fn eval(&self, x: f64) -> f64 { self.polynomial_at(x).eval(x) } +} - fn eval_vec + Copy>(&self, v: &[T]) -> Vec { - let mut result = vec![0f64; v.len()]; - - for (i, x) in v.iter().enumerate() { - result[i] = self.eval(*x); - } - - result - } - +pub trait PolynomialSpline { fn polynomial_at + Copy>(&self, x: T) -> &Polynomial { let x = x.into(); @@ -271,14 +341,6 @@ pub trait Spline { } fn get_ranged_polynomials(&self) -> &Vec<(Range, Polynomial)>; - - fn eval_with_cond f64>(&self, x: f64, cond: F) -> f64 { - cond(self.eval(x)) - } - - fn eval_vec_with_cond f64 + Copy>(&self, x: &[f64], cond: F) -> Vec { - x.iter().map(|&x| self.eval_with_cond(x, cond)).collect() - } } // ============================================================================= @@ -383,7 +445,7 @@ pub struct CubicSpline { polynomials: Vec<(Range, Polynomial)>, } -impl Spline for CubicSpline { +impl PolynomialSpline for CubicSpline { fn get_ranged_polynomials(&self) -> &Vec<(Range, Polynomial)> { &self.polynomials } @@ -633,7 +695,7 @@ pub struct CubicHermiteSpline { polynomials: Vec<(Range, Polynomial)>, } -impl Spline for CubicHermiteSpline { +impl PolynomialSpline for CubicHermiteSpline { fn get_ranged_polynomials(&self) -> &Vec<(Range, Polynomial)> { &self.polynomials } @@ -806,8 +868,8 @@ fn akima_slopes(x: &[f64], y: &[f64]) -> Result> { let y_f = l_f.eval(x_f); let y_ff = l_f.eval(x_ff); - let new_x = concat(&concat(&vec![x_ii, x_i], &x.to_vec()), &vec![x_f, x_ff]); - let new_y = concat(&concat(&vec![y_ii, y_i], &y.to_vec()), &vec![y_f, y_ff]); + let new_x = concat(&concat(&[x_ii, x_i], x), &[x_f, x_ff]); + let new_y = concat(&concat(&[y_ii, y_i], y), &[y_f, y_ff]); for i in 0..new_x.len() - 1 { let dx = new_x[i + 1] - new_x[i]; @@ -817,12 +879,12 @@ fn akima_slopes(x: &[f64], y: &[f64]) -> Result> { s[i] = (new_y[i + 1] - new_y[i]) / dx; } - for i in 0..x.len() { + for (i, m_i) in m.iter_mut().enumerate() { let j = i + 2; let ds_f = (s[j + 1] - s[j]).abs(); let ds_i = (s[j - 1] - s[j - 2]).abs(); - m[i] = if ds_f == 0f64 && ds_i == 0f64 { + *m_i = if ds_f == 0f64 && ds_i == 0f64 { (s[j - 1] + s[j]) / 2f64 } else { (ds_f * s[j - 1] + ds_i * s[j]) / (ds_f + ds_i) @@ -987,3 +1049,145 @@ impl CubicBSplineBases { x.iter().map(|x| self.eval(*x)).collect() } } + +/// B-Spline +/// +/// # Description +/// : B-Spline is a generalization of Bezier curve. +/// +/// # Caution +/// - Let K = the number of knots, C = the number of control points +/// - For open, K = C + degree + 1 (C = K - degree - 1) +/// - For clamped, K + 2 * degree = C + degree + 1 (C = K + degree - 1) +/// +/// # Example +/// ``` +/// use peroxide::fuga::*; +/// use anyhow::Result; +/// +/// fn main() -> Result<()> { +/// let degree = 3; +/// let knots = vec![0f64, 0.5, 1f64]; +/// let control_points = vec![ +/// vec![0f64, 0f64], +/// vec![0.3, 0.5], +/// vec![0.5, 0.7], +/// vec![0.7, 0.3], +/// vec![1f64, 1f64], +/// ]; +/// let spline = BSpline::clamped(degree, knots, control_points)?; +/// +/// let t = linspace(0f64, 1f64, 100); +/// let (x, y): (Vec, Vec) = spline.eval_vec(&t).into_iter().unzip(); +/// assert_eq!(x.len(), 100); +/// assert_eq!(y.len(), 100); +/// +/// Ok(()) +/// } +/// ``` +#[derive(Debug, Clone, Default)] +#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] +pub struct BSpline { + pub degree: usize, + pub knots: Vec, + pub control_points: Vec>, +} + +impl BSpline { + /// Create new open B-Spline + /// + /// # Arguments + /// - `degree` - Degree of B-Spline + /// - `knots` - Knots (length = K) + /// - `control_points` - Control points (length = C) + /// + /// # Caution + /// - K = C + degree + 1 (C = K - degree - 1) + pub fn open(degree: usize, knots: Vec, control_points: Vec>) -> Result { + if knots.len() != control_points.len() + degree + 1 { + bail!("The number of knots ({}) should be equal to the number of control points ({}) + degree ({}) + 1", knots.len(), control_points.len(), degree); + } + + Ok(Self { degree, knots, control_points }) + } + + /// Create new clamped B-Spline + /// + /// # Arguments + /// - `degree` - Degree of B-Spline + /// - `knots` - Knots (length = K) + /// - `control_points` - Control points (length = C) + /// + /// # Caution + /// - K + 2 * degree = C + degree + 1 (C = K + degree - 1) + pub fn clamped(degree: usize, knots: Vec, control_points: Vec>) -> Result { + let mut knots = knots.clone(); + + if knots.len() != control_points.len() + 1 - degree { + bail!("For clamped, the number of knots ({}) should be equal to the number of control points ({}) + 1 - degree ({})", knots.len(), control_points.len(), degree); + } + for _ in 0..degree { + knots.insert(0, knots[0]); + knots.push(knots[knots.len() - 1]); + } + if knots.len() != control_points.len() + degree + 1 { + bail!("The number of knots ({}) should be equal to the number of control points ({}) + degree ({}) + 1", knots.len(), control_points.len(), degree); + } + + Ok(Self { degree, knots, control_points }) + } + + /// Obtain basis function via Cox-de Boor algorithm + #[allow(non_snake_case)] + pub fn cox_de_boor(&self, t: f64, i: usize) -> f64 { + let p = self.degree; + let mut B = vec![vec![0f64; p + 1]; p + 1]; + + // Initialize 0th degree basis functions + for (j, B_j) in B.iter_mut().enumerate() { + if (self.knots[i + j] <= t && t < self.knots[i + j + 1]) || (i + j == self.knots.len() - (p + 2) && t == self.knots[i + j + 1]) { + B_j[0] = 1f64; + } else { + B_j[0] = 0f64; + } + } + + // Compute the basis functions for higher degree + for k in 1..=p { + for j in 0..=(p-k) { + let a = if self.knots[i + j + k] == self.knots[i+j] { + 0f64 + } else { + (t - self.knots[i+j]) / (self.knots[i+j+k] - self.knots[i+j]) + }; + + let b = if self.knots[i + j + k + 1] == self.knots[i+j+1] { + 0f64 + } else { + (self.knots[i+j+k+1] - t) / (self.knots[i+j+k+1] - self.knots[i+j+1]) + }; + + B[j][k] = a * B[j][k-1] + b * B[j+1][k-1]; + } + } + + B[0][p] + } +} + +impl Spline<(f64, f64)> for BSpline { + #[allow(non_snake_case)] + fn eval(&self, t: f64) -> (f64, f64) { + let n = self.control_points.len(); + + let mut x = 0f64; + let mut y = 0f64; + for i in 0 .. n { + let B = self.cox_de_boor(t, i); + x += B * self.control_points[i][0]; + y += B * self.control_points[i][1]; + } + + (x, y) + } +}