diff --git a/CHANGELOG.md b/CHANGELOG.md index bfad977..800fbde 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,17 @@ +# 0.5.0 + - update `BuilderError` + - combine variants `DimensionError` and `AxisLenght` into `ShapeError` + - add `ValueError` variant + - update `CubicSpline` stragegie + - Move `CubicSpline` interpolator to `interp1d::cubic_spline` module + - add extrapolation + - add not-a-knot boundary condition + - add clamped boundary condition + - add periodic boundary condition + - make not-a-knot boundary condition the default + - allow any first or second derivative as boundary condition + - fix typo `Biliniar` -> `Bilinear` + # 0.4.1 - major performance improvement for `interp_scalar()` methods ~-50% - keywords and categorys in crate metadata @@ -37,7 +51,7 @@ logarithmic spaced values. - updated package structure - replaced Interp1DStrategy enum with individual structs - added Strategy and StrategyBuilder trait - - added QubicSpline strategy + - added CubicSpline strategy - added traits for custom strategies # 0.1.1 diff --git a/Cargo.toml b/Cargo.toml index 63c0ccc..2a78ec9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "ndarray-interp" -version = "0.4.1" +version = "0.5.0" edition = "2021" license = "MIT" description = "Interpolation package for ndarray" diff --git a/src/interp1d.rs b/src/interp1d.rs index 95e0939..75c3473 100644 --- a/src/interp1d.rs +++ b/src/interp1d.rs @@ -10,7 +10,7 @@ //! //! # Strategies //! - [`Linear`] Linear interpolation strategy -//! - [`CubicSpline`] Cubic spline interpolation strategy +//! - [`cubic_spline`] Cubic spline interpolation strategy use std::{any::TypeId, fmt::Debug, ops::Sub}; @@ -30,7 +30,14 @@ use crate::{ mod aliases; mod strategies; pub use aliases::*; -pub use strategies::{CubicSpline, Interp1DStrategy, Interp1DStrategyBuilder, Linear}; +pub use strategies::linear::Linear; +pub use strategies::{Interp1DStrategy, Interp1DStrategyBuilder}; + +pub mod cubic_spline { + pub use super::strategies::cubic_spline::{ + BoundaryCondition, CubicSpline, RowBoundary, SingleBoundary, + }; +} /// One dimensional interpolator #[derive(Debug)] @@ -444,7 +451,7 @@ where let Interp1DBuilder { x, data, strategy } = self; if data.ndim() < 1 { - return Err(DimensionError( + return Err(ShapeError( "data dimension is 0, needs to be at least 1".into(), )); } @@ -460,7 +467,7 @@ where )); } if x.len() != data.shape()[0] { - return Err(BuilderError::AxisLenght(format!( + return Err(BuilderError::ShapeError(format!( "Lengths of x and data axis need to match. Got x: {:}, data: {:}", x.len(), data.shape()[0], diff --git a/src/interp1d/strategies.rs b/src/interp1d/strategies.rs index b873125..73f4c11 100644 --- a/src/interp1d/strategies.rs +++ b/src/interp1d/strategies.rs @@ -6,11 +6,8 @@ use num_traits::Num; use super::Interp1D; use crate::{BuilderError, InterpolateError}; -mod cubic_spline; -mod linear; - -pub use cubic_spline::CubicSpline; -pub use linear::Linear; +pub mod cubic_spline; +pub mod linear; pub trait Interp1DStrategyBuilder where diff --git a/src/interp1d/strategies/cubic_spline.rs b/src/interp1d/strategies/cubic_spline.rs index ebdaae4..7a4f0a2 100644 --- a/src/interp1d/strategies/cubic_spline.rs +++ b/src/interp1d/strategies/cubic_spline.rs @@ -1,12 +1,13 @@ use std::{ fmt::Debug, - ops::{Add, Sub, SubAssign}, + ops::{Add, Neg, Sub, SubAssign}, }; use ndarray::{ - s, Array, ArrayBase, ArrayViewMut, Axis, Data, Dimension, Ix1, RemoveAxis, ScalarOperand, Zip, + s, Array, Array1, ArrayBase, ArrayView, ArrayViewMut, Axis, Data, Dimension, FoldWhile, Ix1, + IxDyn, RemoveAxis, ScalarOperand, Slice, Zip, }; -use num_traits::{cast, Num, NumCast, Pow}; +use num_traits::{cast, Euclid, Num, NumCast, Pow}; use crate::{interp1d::Interp1D, BuilderError, InterpolateError}; @@ -14,6 +15,43 @@ use super::{Interp1DStrategy, Interp1DStrategyBuilder}; const AX0: Axis = Axis(0); +/// Marker trait that is implemented for anithing that satisfies +/// the trait bounds required to be used as an element in the QubicSpline +/// strategy. +pub trait SplineNum: + Debug + + Num + + Copy + + PartialOrd + + Sub + + SubAssign + + Neg + + NumCast + + Add + + Pow + + ScalarOperand + + Euclid + + Send +{ +} + +impl SplineNum for T where + T: Debug + + Num + + Copy + + PartialOrd + + Sub + + SubAssign + + Neg + + NumCast + + Add + + Pow + + ScalarOperand + + Euclid + + Send +{ +} + /// The CubicSpline 1d interpolation Strategy /// /// # Example @@ -21,6 +59,7 @@ const AX0: Axis = Axis(0); /// ``` /// # use ndarray_interp::*; /// # use ndarray_interp::interp1d::*; +/// # use ndarray_interp::interp1d::cubic_spline::*; /// # use ndarray::*; /// # use approx::*; /// @@ -28,41 +67,209 @@ const AX0: Axis = Axis(0); /// let x = array![-1.0, 0.0, 3.0]; /// let query = Array::linspace(-1.0, 3.0, 10); /// let interpolator = Interp1DBuilder::new(y) -/// .strategy(CubicSpline) +/// .strategy(CubicSpline::new()) /// .x(x) /// .build().unwrap(); /// /// let result = interpolator.interp_array(&query).unwrap(); /// let expect = array![ /// 0.5, -/// 0.2109053497942387, -/// 0.020576131687242816, -/// 0.01851851851851849, -/// 0.21364883401920443, -/// 0.5733882030178327, -/// 1.0648148148148144, -/// 1.6550068587105617, -/// 2.3110425240054866, -/// 3.0 +/// 0.1851851851851852, +/// 0.01851851851851853, +/// -5.551115123125783e-17, +/// 0.12962962962962965, +/// 0.40740740740740755, +/// 0.8333333333333331, +/// 1.407407407407407, +/// 2.1296296296296293, 3.0 /// ]; /// # assert_abs_diff_eq!(result, expect, epsilon=f64::EPSILON); /// ``` #[derive(Debug)] -pub struct CubicSpline; -impl Interp1DStrategyBuilder for CubicSpline +pub struct CubicSpline { + extrapolate: bool, + boundary: BoundaryCondition, +} + +/// Boundary conditions for the whole dataset +/// +/// The boundary condition is structured in three hirarchic enum's: +/// - [`BoundaryCondition`] The toplevel boundary applys to the whole dataset +/// - [`RowBoundary`] applys to a single row in the dataset +/// - [`SingleBoundary`] applys to an individual boundary of a single row +/// +/// the default is the [`NotAKnot`](BoundaryCondition::NotAKnot) boundary in each level +/// +/// There are different possibilities for the boundary condition in each level: +/// - [`NotAKnot`](BoundaryCondition::NotAKnot) - all levels +/// - [`Natural`](BoundaryCondition::Natural) - all levels (same as `SecondDeriv(0.0)`) +/// - [`Clamped`](BoundaryCondition::Clamped) - all levels (same as `FirstDeriv(0.0)`) +/// - [`Periodic`](BoundaryCondition::Periodic) - not in [`SingleBoundary`] +/// - [`FirstDeriv`](SingleBoundary::FirstDeriv) - only in [`SingleBoundary`] +/// - [`SecondDeriv`](SingleBoundary::SecondDeriv) - only in [`SingleBoundary`] +/// +/// ## Example +/// In a complex case all boundaries can be set individually: +/// ``` rust +/// # use ndarray_interp::*; +/// # use ndarray_interp::interp1d::*; +/// # use ndarray_interp::interp1d::cubic_spline::*; +/// # use ndarray::*; +/// # use approx::*; +/// +/// let y = array![ +/// [0.5, 1.0], +/// [0.0, 1.5], +/// [3.0, 0.5], +/// ]; +/// let x = array![-1.0, 0.0, 3.0]; +/// +/// // first data column: natural +/// // second data column top: NotAKnot +/// // second data column bottom: first derivative == 0.5 +/// let boundaries = array![ +/// [ +/// RowBoundary::Natural, +/// RowBoundary::Mixed { left: SingleBoundary::NotAKnot, right: SingleBoundary::FirstDeriv(0.5)} +/// ], +/// ]; +/// let strat = CubicSpline::new().boundary(BoundaryCondition::Individual(boundaries)); +/// let interpolator = Interp1DBuilder::new(y) +/// .x(x) +/// .strategy(strat) +/// .build().unwrap(); +/// +/// ``` +#[derive(Debug, PartialEq, Eq)] +pub enum BoundaryCondition { + /// Not a knot boundary. The first and second segment at a curve end are the same polynomial. + NotAKnot, + /// Natural boundary. The second derivative at the curve end is 0 + Natural, + /// Clamped boundary. The first derivative at the curve end is 0 + Clamped, + /// Periodic spline. + /// The interpolated functions is assumed to be periodic. + /// The first and last element in the data must be equal. + Periodic, + /// Set individual boundary conditions for each row in the data + /// and/or individual conditions for the left and right boundary + Individual(Array, D>), +} + +impl Default for BoundaryCondition { + fn default() -> Self { + Self::NotAKnot + } +} + +/// Boundary condition for a single data row +#[derive(Debug, PartialEq, Eq, Clone)] +pub enum RowBoundary { + /// ![`BoundaryCondition::NotAKnot`] + NotAKnot, + /// ![`BoundaryCondition::Natural`] + Natural, + /// ![`BoundaryCondition::Clamped`] + Clamped, + /// Set individual boundary conditions at the left and right end of the curve + Mixed { + left: SingleBoundary, + right: SingleBoundary, + }, +} + +impl Default for RowBoundary { + fn default() -> Self { + Self::NotAKnot + } +} + +impl From> for InternalBoundary { + fn from(val: RowBoundary) -> Self { + match val { + RowBoundary::NotAKnot => InternalBoundary::NotAKnot, + RowBoundary::Natural => InternalBoundary::Natural, + RowBoundary::Clamped => InternalBoundary::Clamped, + RowBoundary::Mixed { left, right } => InternalBoundary::Mixed { left, right }, + } + } +} + +/// This is essentially [`RowBoundary`] but including the Periodic variant. +/// The periodic variant can not be applied to a single row only all or nothing. +/// But we still need it for calculating the coefficients, which may or may not be done +/// for each row individually. +enum InternalBoundary { + NotAKnot, + Natural, + Clamped, + Periodic, + Mixed { + left: SingleBoundary, + right: SingleBoundary, + }, +} + +impl InternalBoundary { + fn specialize(self) -> Self { + use SingleBoundary::*; + match self { + InternalBoundary::Natural => Self::Mixed { + left: Natural, + right: Natural, + }, + InternalBoundary::NotAKnot => Self::Mixed { + left: NotAKnot, + right: NotAKnot, + }, + InternalBoundary::Clamped => Self::Mixed { + left: Clamped, + right: Clamped, + }, + _ => self, + } + } +} + +/// Boundary condition for a single boundary (one side of one data row) +#[derive(Debug, PartialEq, Eq, Clone)] +pub enum SingleBoundary { + /// ![`BoundaryCondition::NotAKnot`] + NotAKnot, + /// This ist the same as `SingleBoundary::SecondDeriv(0.0)` + /// ![`BoundaryCondition::Natural`] + Natural, + /// This ist the same as `SingleBoundary::FirstDeriv(0.0)` + /// ![`BoundaryCondition::Clamped`] + Clamped, + /// Set a value for the first derivative at the curve end + FirstDeriv(T), + /// Set a value for the second derivative at the curve end + SecondDeriv(T), +} + +impl SingleBoundary { + fn specialize(self) -> Self { + use SingleBoundary::*; + match self { + SingleBoundary::Natural => SecondDeriv(cast(0.0).unwrap_or_else(|| unimplemented!())), + SingleBoundary::Clamped => FirstDeriv(cast(0.0).unwrap_or_else(|| unimplemented!())), + _ => self, + } + } +} + +impl Default for SingleBoundary { + fn default() -> Self { + Self::NotAKnot + } +} + +impl Interp1DStrategyBuilder for CubicSpline where Sd: Data, - Sd::Elem: Debug - + Num - + Copy - + PartialOrd - + Sub - + SubAssign - + NumCast - + Add - + Pow - + ScalarOperand - + Send, + Sd::Elem: SplineNum, Sx: Data, D: Dimension + RemoveAxis, { @@ -77,43 +284,142 @@ where where Sx2: Data, { - let (a, b) = self.calc_coefficients(x, data); - Ok(CubicSplineStrategy { a, b }) + let (a, b) = self.calc_coefficients(x, data)?; + let extrapolate = if !self.extrapolate { + Extrapolate::No + } else if matches!(self.boundary, BoundaryCondition::Periodic) { + Extrapolate::Periodic + } else { + Extrapolate::Yes + }; + Ok(CubicSplineStrategy { a, b, extrapolate }) } } -impl CubicSpline { - fn calc_coefficients( - self, +impl CubicSpline +where + D: Dimension + RemoveAxis, + T: SplineNum, +{ + /// Calculate the coefficients `a` and `b` + fn calc_coefficients( + &self, x: &ArrayBase, data: &ArrayBase, - ) -> (Array, Array) + ) -> Result<(Array, Array), BuilderError> where - Sd: Data, - Sd::Elem: Num - + Copy - + Sub - + SubAssign - + NumCast - + Add - + Pow - + ScalarOperand - + Debug, - Sx: Data, - D: Dimension + RemoveAxis, + Sd: Data, + Sx: Data, { let dim = data.raw_dim(); let len = dim[0]; + let mut k = Array::zeros(dim.clone()); + let kv = k.view_mut(); + match self.boundary { + BoundaryCondition::Periodic => { + Self::solve_for_k(kv, x, data, InternalBoundary::Periodic) + } + BoundaryCondition::Natural => Self::solve_for_k(kv, x, data, InternalBoundary::Natural), + BoundaryCondition::Clamped => Self::solve_for_k(kv, x, data, InternalBoundary::Clamped), + BoundaryCondition::NotAKnot => { + Self::solve_for_k(kv, x, data, InternalBoundary::NotAKnot) + } + BoundaryCondition::Individual(ref bounds) => { + let mut bounds_shape = kv.raw_dim(); + bounds_shape[0] = 1; + if bounds_shape != bounds.raw_dim() { + return Err(BuilderError::ShapeError(format!( + "Boundary conditions array has wrong shape. Expected: {bounds_shape:?}, got: {:?}", + bounds.raw_dim() + ))); + } + Self::solve_for_k_individual( + kv.into_dyn(), + x, + data.view().into_dyn(), + bounds.view().into_dyn(), + ) + } + }?; + let mut a_b_dim = data.raw_dim(); a_b_dim[0] -= 1; + let mut c_a = Array::zeros(a_b_dim.clone()); + let mut c_b = Array::zeros(a_b_dim); + for index in 0..len - 1 { + Zip::from(c_a.index_axis_mut(AX0, index)) + .and(c_b.index_axis_mut(AX0, index)) + .and(k.index_axis(AX0, index)) + .and(k.index_axis(AX0, index + 1)) + .and(data.index_axis(AX0, index)) + .and(data.index_axis(AX0, index + 1)) + .for_each(|c_a, c_b, &k, &k_right, &y, &y_right| { + *c_a = k * (x[index + 1] - x[index]) - (y_right - y); + *c_b = (y_right - y) - k_right * (x[index + 1] - x[index]); + }) + } + + Ok((c_a, c_b)) + } + + fn solve_for_k_individual( + mut k: ArrayViewMut, + x: &ArrayBase, + data: ArrayView, + boundary: ArrayView, IxDyn>, + ) -> Result<(), BuilderError> + where + Sx: Data, + { + if k.ndim() > 1 { + let ax = Axis(k.ndim() - 1); + Zip::from(k.axis_iter_mut(ax)) + .and(data.axis_iter(ax)) + .and(boundary.axis_iter(ax)) + .fold_while(Ok(()), |_, k, data, boundary| { + Self::solve_for_k_individual(k, x, data, boundary).map_or_else( + |err| FoldWhile::Done(Err(err)), + |_| FoldWhile::Continue(Ok(())), + ) + }) + .into_inner() + } else { + Self::solve_for_k( + k, + x, + &data, + boundary + .first() + .cloned() + .unwrap_or_else(|| unreachable!()) + .into(), + ) + } + } + + /// solves the linear equation `A * k = rhs` with the [`RowBoundary`] used for + /// each row in the data + /// + /// **returns** k + fn solve_for_k( + mut k: ArrayViewMut, + x: &ArrayBase, + data: &ArrayBase, + boundary: InternalBoundary, + ) -> Result<(), BuilderError> + where + _D: Dimension + RemoveAxis, + Sd: Data, + Sx: Data, + { + let dim = data.raw_dim(); + let len = dim[0]; /* * Calculate the coefficients c_a and c_b for the cubic spline the method is outlined on * https://en.wikipedia.org/wiki/Spline_interpolation#Example * * This requires solving the Linear equation A * k = rhs - * The Thomas algorithm is used, because the matrix A will be tridiagonal and diagonally dominant. - * (https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm) */ // upper, middle and lower diagonal of A @@ -121,79 +427,260 @@ impl CubicSpline { let mut a_mid = Array::zeros(len); let mut a_low = Array::zeros(len); - let one: Sd::Elem = cast(1.0).unwrap_or_else(|| unimplemented!()); - let two: Sd::Elem = cast(2.0).unwrap_or_else(|| unimplemented!()); - let three: Sd::Elem = cast(3.0).unwrap_or_else(|| unimplemented!()); + let zero: T = cast(0.0).unwrap_or_else(|| unimplemented!()); + let one: T = cast(1.0).unwrap_or_else(|| unimplemented!()); + let two: T = cast(2.0).unwrap_or_else(|| unimplemented!()); + let three: T = cast(3.0).unwrap_or_else(|| unimplemented!()); Zip::from(a_up.slice_mut(s![1..-1])) .and(a_mid.slice_mut(s![1..-1])) .and(a_low.slice_mut(s![1..-1])) .and(x.windows(3)) .for_each(|a_up, a_mid, a_low, x| { - let x_left = x[0]; - let x_mid = x[1]; - let x_right = x[2]; + let dxn = x[2] - x[1]; + let dxn_1 = x[1] - x[0]; - *a_up = one / (x_right - x_mid); - *a_mid = two / (x_mid - x_left) + two / (x_right - x_mid); - *a_low = one / (x_mid - x_left); + *a_up = dxn_1; + *a_mid = two * (dxn + dxn_1); + *a_low = dxn; }); - let x_0 = x[0]; - let x_1 = x[1]; - - a_up[0] = one / (x_1 - x_0); - a_mid[0] = two / (x_1 - x_0); + // RHS vector + let mut rhs = Array::zeros(dim.clone()); - // x_n and xn-1 - let x_n = x[len - 1]; - let x_n1 = x[len - 2]; - a_mid[len - 1] = two / (x_n - x_n1); - a_low[len - 1] = one / (x_n - x_n1); + for n in 1..len - 1 { + let rhs = rhs.index_axis_mut(AX0, n); + let y_left = data.index_axis(AX0, n - 1); + let y_mid = data.index_axis(AX0, n); + let y_right = data.index_axis(AX0, n + 1); - // RHS vector - let mut rhs: Array = Array::zeros(dim.clone()); + let dxn = x[n + 1] - x[n]; // dx(n) + let dxn_1 = x[n] - x[n - 1]; // dx(n-1) - for i in 1..len - 1 { - let rhs = rhs.index_axis_mut(AX0, i); - let y_left = data.index_axis(AX0, i - 1); - let y_mid = data.index_axis(AX0, i); - let y_right = data.index_axis(AX0, i + 1); - let x_left = x[i - 1]; - let x_mid = x[i]; - let x_right = x[i + 1]; Zip::from(y_left).and(y_mid).and(y_right).map_assign_into( rhs, |&y_left, &y_mid, &y_right| { - three - * ((y_mid - y_left) / (x_mid - x_left).pow(two) - + (y_right - y_mid) / (x_right - x_mid).pow(two)) + three * (dxn * (y_mid - y_left) / dxn_1 + dxn_1 * (y_right - y_mid) / dxn) }, ); } - let rhs_0 = rhs.index_axis_mut(AX0, 0); - let data_0 = data.index_axis(AX0, 0); - let data_1 = data.index_axis(AX0, 1); - Zip::from(rhs_0) - .and(data_0) - .and(data_1) - .for_each(|rhs_0, &y_0, &y_1| { - *rhs_0 = three * (y_1 - y_0) / (x_1 - x_0).pow(two); - }); - - let rhs_n = rhs.index_axis_mut(AX0, len - 1); - let data_n = data.index_axis(AX0, len - 1); - let data_n1 = data.index_axis(AX0, len - 2); - Zip::from(rhs_n) - .and(data_n) - .and(data_n1) - .for_each(|rhs_n, &y_n, &y_n1| { - *rhs_n = three * (y_n - y_n1) / (x_n - x_n1).pow(two); - }); - - // now solving With Thomas algorithm + let dx0 = x[1] - x[0]; + let dx1 = x[2] - x[1]; + let dx_1 = x[len - 1] - x[len - 2]; + let dx_2 = x[len - 2] - x[len - 3]; + + // apply boundary conditions + match (boundary.specialize(), len) { + (InternalBoundary::Periodic, 3) => { + let y0 = data.index_axis(AX0, 0); + let y2 = data.index_axis(AX0, 2); + if y0 != y2 { + if data.ndim() == 1 { + return Err(BuilderError::ValueError(format!("for periodic boundary condition the first and last value must be equal. First: {:?}, last: {:?}", data.first().unwrap_or_else(||unreachable!()), data.last().unwrap_or_else(||unreachable!())))); + } else { + return Err(BuilderError::ValueError(format!("for periodic boundary condition the first and last value must be equal. First: {y0:?}, last: {y2:?}"))); + } + } + + let y1 = data.index_axis(AX0, 1); + let slope0: Array = (&y1 - &y0) / dx0; + let slope1: Array = (&y2 - &y1) / dx1; + k.assign(&((slope0 / dx0 + slope1 / dx1) / (one / dx0 + one / dx1))); + return Ok(()); + } + + (InternalBoundary::Periodic, _) => { + let y0 = data.index_axis(AX0, 0); + let y_1 = data.index_axis(AX0, len - 1); + if y0 != y_1 { + if data.ndim() == 1 { + return Err(BuilderError::ValueError(format!("for periodic boundary condition the first and last value must be equal. First: {:?}, last: {:?}", data.first().unwrap_or_else(||unreachable!()), data.last().unwrap_or_else(||unreachable!())))); + } else { + return Err(BuilderError::ValueError(format!("for periodic boundary condition the first and last value must be equal. First: {y0:?}, last: {y_1:?}"))); + } + } + + // due to the preriodicity we need to solve one less equation + // the system matrix a is also condensed + // https://web.archive.org/web/20151220180652/http://www.cfm.brown.edu/people/gk/chap6/node14.html + a_up.slice_axis_inplace(AX0, Slice::from(0..-2)); + a_mid.slice_axis_inplace(AX0, Slice::from(0..-2)); + a_low.slice_axis_inplace(AX0, Slice::from(0..-2)); + rhs.slice_axis_inplace(AX0, Slice::from(0..-1)); + + a_mid[0] = two * (dx_1 + dx0); + a_up[0] = dx_1; + + let y1 = data.index_axis(AX0, 1); + let slope0: Array = (&y1 - &y0) / dx0; + + let y_1 = data.index_axis(AX0, len - 1); + let y_2 = data.index_axis(AX0, len - 2); + let y_3 = data.index_axis(AX0, len - 3); + let slope_1: Array = (&y_1 - &y_2) / dx_1; + let slope_2: Array = (&y_2 - &y_3) / dx_2; + + rhs.index_axis_mut(AX0, 0) + .assign(&((&slope_1 * dx0 + &slope0 * dx_1) * three)); + rhs.index_axis_mut(AX0, len - 1 - 1) + .assign(&((slope_2 * dx_1 + slope_1 * dx_2) * three)); + + let rhs1 = rhs.slice_axis(AX0, Slice::from(0..-1)).to_owned(); + let mut rhs2 = Array::zeros(rhs1.raw_dim()); + rhs2.index_axis_mut(AX0, 0).fill(-dx0); // = -dx0; + let dx_3 = x[len - 3] - x[len - 4]; + rhs2.index_axis_mut(AX0, len - 3).fill(-dx_3); + + let mut k1 = Array::zeros(rhs1.raw_dim()); + let mut k2 = Array::zeros(rhs1.raw_dim()); + + Self::thomas( + k1.view_mut(), + a_up.clone(), + a_mid.clone(), + a_low.clone(), + rhs1, + ); + Self::thomas(k2.view_mut(), a_up, a_mid, a_low, rhs2); + + let k_m1 = (&rhs.index_axis(AX0, len - 2) + - &k1.index_axis(AX0, 0) * dx_2 + - &k1.index_axis(AX0, len - 3) * dx_1) + / (&k2.index_axis(AX0, 0) * dx_2 + + &k2.index_axis(AX0, len - 3) * dx_1 + + two * (dx_1 + dx_2)); + + k.slice_axis_mut(AX0, Slice::from(0..-2)) + .assign(&(k1 + &k_m1 * k2)); + k.index_axis_mut(AX0, len - 2).assign(&k_m1); + let k0 = k.index_axis(AX0, 0).to_owned(); + k.index_axis_mut(AX0, len - 1).assign(&k0); + return Ok(()); + } + (InternalBoundary::Clamped, _) => unreachable!(), + (InternalBoundary::Natural, _) => unreachable!(), + (InternalBoundary::NotAKnot, _) => unreachable!(), + ( + InternalBoundary::Mixed { + left: SingleBoundary::NotAKnot, + right: SingleBoundary::NotAKnot, + }, + 3, + ) => { + // We handle this case by constructing a parabola passing through given points. + + let y0 = data.index_axis(AX0, 0); + let y1 = data.index_axis(AX0, 1); + let y2 = data.index_axis(AX0, 2); + let slope0 = (y1.to_owned() - y0) / dx0; + let slope1 = (y2.to_owned() - y1) / dx1; + + a_mid[0] = one; // [0, 0] + a_up[0] = one; // [0, 1] + a_low[1] = dx1; // [1, 0] + a_mid[1] = two * (dx0 + dx1); // [1, 1] + a_up[1] = dx0; // [1, 2] + a_low[2] = one; // [2, 1] + a_mid[2] = one; // [2, 2] + + rhs.index_axis_mut(AX0, 0).assign(&(&slope0 * two)); + rhs.index_axis_mut(AX0, 1) + .assign(&((&slope1 * dx0 + &slope0 * dx1) * three)); + rhs.index_axis_mut(AX0, 2).assign(&(slope1 * two)); + } + (InternalBoundary::Mixed { left, right }, _) => { + match left.specialize() { + SingleBoundary::NotAKnot => { + a_mid[0] = dx1; + let d = x[2] - x[0]; + a_up[0] = d; + let tmp1 = (dx0 + two * d) * dx1; + Zip::from(rhs.index_axis_mut(AX0, 0)) + .and(data.index_axis(AX0, 0)) + .and(data.index_axis(AX0, 1)) + .and(data.index_axis(AX0, 2)) + .for_each(|b, &y0, &y1, &y2| { + *b = (tmp1 * (y1 - y0) / dx0 + dx0.pow(two) * (y2 - y1) / dx1) / d; + }); + } + SingleBoundary::Natural => unreachable!(), + SingleBoundary::Clamped => unreachable!(), + SingleBoundary::FirstDeriv(deriv) => { + a_mid[0] = one; + a_up[0] = zero; + rhs.index_axis_mut(AX0, 0).fill(deriv); + } + SingleBoundary::SecondDeriv(deriv) => { + a_up[0] = dx0; + a_mid[0] = two * dx0; + let rhs_0 = rhs.index_axis_mut(AX0, 0); + let data_0 = data.index_axis(AX0, 0); + let data_1 = data.index_axis(AX0, 1); + Zip::from(rhs_0) + .and(data_0) + .and(data_1) + .for_each(|rhs_0, &y_0, &y_1| { + *rhs_0 = three * (y_1 - y_0) - deriv * dx0.pow(two) / two; + }); + } + }; + match right.specialize() { + SingleBoundary::NotAKnot => { + a_mid[len - 1] = dx_1; + let d = x[len - 1] - x[len - 3]; + a_low[len - 1] = d; + let tmp1 = (two * d + dx_1) * dx_2; + Zip::from(rhs.index_axis_mut(AX0, len - 1)) + .and(data.index_axis(AX0, len - 1)) + .and(data.index_axis(AX0, len - 2)) + .and(data.index_axis(AX0, len - 3)) + .for_each(|b, &y_1, &y_2, &y_3| { + *b = (dx_1.pow(two) * (y_2 - y_3) / dx_2 + + tmp1 * (y_1 - y_2) / dx_1) + / d; + }); + } + SingleBoundary::Natural => unreachable!(), + SingleBoundary::Clamped => unreachable!(), + SingleBoundary::FirstDeriv(deriv) => { + a_mid[len - 1] = one; + a_low[len - 1] = zero; + rhs.index_axis_mut(AX0, len - 1).fill(deriv); + } + SingleBoundary::SecondDeriv(deriv) => { + a_mid[len - 1] = two * dx_1; + a_low[len - 1] = dx_1; + let rhs_n = rhs.index_axis_mut(AX0, len - 1); + let data_n = data.index_axis(AX0, len - 1); + let data_n1 = data.index_axis(AX0, len - 2); + Zip::from(rhs_n) + .and(data_n) + .and(data_n1) + .for_each(|rhs_n, &y_n, &y_n1| { + *rhs_n = three * (y_n - y_n1) + deriv * dx_1.pow(two) / two; + }); + } + }; + } + } + Self::thomas(k, a_up, a_mid, a_low, rhs); + Ok(()) + } + /// The Thomas algorithm is used, because the matrix A will be tridiagonal and diagonally dominant + /// [https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm] + fn thomas<_D>( + mut k: ArrayViewMut, + a_up: Array1, + mut a_mid: Array1, + a_low: Array1, + mut rhs: Array, + ) where + _D: Dimension + RemoveAxis, + { + let dim = rhs.raw_dim(); + let len = dim[0]; let mut rhs_left = rhs.index_axis(AX0, 0).into_owned(); for i in 1..len { let w = a_low[i] / a_mid[i - 1]; @@ -209,7 +696,6 @@ impl CubicSpline { }); } - let mut k = Array::zeros(dim); Zip::from(k.index_axis_mut(AX0, len - 1)) .and(rhs.index_axis(AX0, len - 1)) .for_each(|k, &rhs| { @@ -227,36 +713,46 @@ impl CubicSpline { *k_right = new_k; }) } + } - let mut c_a = Array::zeros(a_b_dim.clone()); - let mut c_b = Array::zeros(a_b_dim); - for index in 0..len - 1 { - Zip::from(c_a.index_axis_mut(AX0, index)) - .and(c_b.index_axis_mut(AX0, index)) - .and(k.index_axis(AX0, index)) - .and(k.index_axis(AX0, index + 1)) - .and(data.index_axis(AX0, index)) - .and(data.index_axis(AX0, index + 1)) - .for_each(|c_a, c_b, &k, &k_right, &y, &y_right| { - *c_a = k * (x[index + 1] - x[index]) - (y_right - y); - *c_b = (y_right - y) - k_right * (x[index + 1] - x[index]); - }) + /// create a cubic-spline interpolation stratgy + pub fn new() -> Self { + Self { + extrapolate: false, + boundary: BoundaryCondition::NotAKnot, } + } - (c_a, c_b) + /// does the strategy extrapolate? Default is `false` + pub fn extrapolate(mut self, extrapolate: bool) -> Self { + self.extrapolate = extrapolate; + self } - pub fn new() -> Self { - Self + /// set the boundary condition. default is [`BoundaryCondition::Natural`] + pub fn boundary(mut self, boundary: BoundaryCondition) -> Self { + self.boundary = boundary; + self } } -impl Default for CubicSpline { +impl Default for CubicSpline +where + D: Dimension + RemoveAxis, + T: SplineNum, +{ fn default() -> Self { Self::new() } } +#[derive(Debug)] +enum Extrapolate { + Yes, + No, + Periodic, +} + #[derive(Debug)] pub struct CubicSplineStrategy where @@ -265,12 +761,13 @@ where { a: Array, b: Array, + extrapolate: Extrapolate, } impl Interp1DStrategy for CubicSplineStrategy where Sd: Data, - Sd::Elem: Num + PartialOrd + NumCast + Copy + Debug + Sub + Send, + Sd::Elem: SplineNum, Sx: Data, D: Dimension + RemoveAxis, { @@ -280,12 +777,20 @@ where target: ArrayViewMut<'_, ::Elem, ::Smaller>, x: ::Elem, ) -> Result<(), InterpolateError> { - if !interp.is_in_range(x) { + let in_range = interp.is_in_range(x); + if matches!(self.extrapolate, Extrapolate::No) && !in_range { return Err(InterpolateError::OutOfBounds(format!( "x = {x:#?} is not in range", ))); } + let mut x = x; + if matches!(self.extrapolate, Extrapolate::Periodic) && !in_range { + let x0 = interp.x[0]; + let xn = interp.x[interp.x.len() - 1]; + x = ((x - x0).rem_euclid(&(xn - x0))) + x0; + } + let idx = interp.get_index_left_of(x); let (x_left, data_left) = interp.index_point(idx); let (x_right, data_right) = interp.index_point(idx + 1); diff --git a/src/interp1d/strategies/linear.rs b/src/interp1d/strategies/linear.rs index c218ec3..78a957a 100644 --- a/src/interp1d/strategies/linear.rs +++ b/src/interp1d/strategies/linear.rs @@ -75,7 +75,7 @@ impl Linear { Self { extrapolate: false } } - /// set the extrapolate property, default is `false` + /// does the strategy extrapolate? Default is `false` pub fn extrapolate(mut self, extrapolate: bool) -> Self { self.extrapolate = extrapolate; self diff --git a/src/interp2d.rs b/src/interp2d.rs index 62bf309..c77817a 100644 --- a/src/interp2d.rs +++ b/src/interp2d.rs @@ -475,9 +475,7 @@ where strategy: stratgy_builder, } = self; if data.ndim() < 2 { - return Err(DimensionError( - "data dimension needs to be at least 2".into(), - )); + return Err(ShapeError("data dimension needs to be at least 2".into())); } if data.shape()[0] < Strat::MINIMUM_DATA_LENGHT { return Err(NotEnoughData(format!("The 0-dimension has not enough data for the chosen interpolation strategy. Provided: {}, Reqired: {}", data.shape()[0], Strat::MINIMUM_DATA_LENGHT))); @@ -486,14 +484,14 @@ where return Err(NotEnoughData(format!("The 1-dimension has not enough data for the chosen interpolation strategy. Provided: {}, Reqired: {}", data.shape()[1], Strat::MINIMUM_DATA_LENGHT))); } if x.len() != data.shape()[0] { - return Err(AxisLenght(format!( + return Err(ShapeError(format!( "Lenghts of x-axis and data-0-axis need to match. Got x: {}, data-0: {}", x.len(), data.shape()[0] ))); } if y.len() != data.shape()[1] { - return Err(AxisLenght(format!( + return Err(ShapeError(format!( "Lenghts of y-axis and data-1-axis need to match. Got y: {}, data-1: {}", y.len(), data.shape()[1] diff --git a/src/lib.rs b/src/lib.rs index 725633e..b8e20a6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -5,13 +5,114 @@ //! The ndarray-interp crate provides interpolation algorithms //! for interpolating _n_-dimesional data. //! -//! 1D and 2D interpolation is supported. See the modules [interp1d] and [interp2d] +//! # 1D Interpolation +//! The [interp1d] module provides the [`Interp1D`](interp1d::Interp1D) interpolator +//! and different interpolation strategies +//! +//! **1D Strategies** +//! - [`interp1d::Linear`] - Linear interpolation and extrapolation +//! - [`interp1d::cubic_spline`] - Cubic Spline interpolation with different boundary conditions. +//! +//! # 2D Interpolation +//! The [interp2d] module provides the [`Interp2D`](interp2d::Interp2D) interpolator +//! and different interpolation strategies +//! +//! **2D Strategies** +//! - [`interp2d::Bilinear`] - Bilinear interpolation and extrapolation //! //! # Custom interpolation strategy //! This crate defines traits to allow implementation of user //! defined interpolation algorithms. -//! see the `custom_strategy.rs` example. +//! A 1D interpolation strategy can be created by implementing the +//! [`Interp1DStrategy`](interp1d::Interp1DStrategy) and +//! [`Interp1DStrategyBuilder`](interp1d::Interp1DStrategyBuilder) traits. +//! A 2D interpolation strategy can be created by implementing the +//! [`Interp2DStrategy`](interp2d::Interp2DStrategy) and +//! [`Interp2DStrategyBuilder`](interp2d::Interp2DStrategyBuilder) traits. +//! +//! See also the `custom_strategy.rs` example. +//! +//! # Examples +//! **1D Example** +//! ``` rust +//! use ndarray_interp::interp1d::*; +//! use ndarray::*; +//! +//! let data = array![0.0, 1.0, 1.5, 1.0, 0.0 ]; +//! let interp = Interp1DBuilder::new(data).build().unwrap(); +//! +//! let result = interp.interp_scalar(3.5).unwrap(); +//! assert!(result == 0.5); +//! let result = interp.interp_array(&array![0.0, 0.5, 1.5]).unwrap(); +//! assert!(result == array![0.0, 0.5, 1.25]) +//! ``` +//! +//! **1D Example with multidimensional data** +//! ```rust +//! use ndarray_interp::interp1d::*; +//! use ndarray::*; +//! +//! let data = array![ +//! [0.0, 1.0], +//! [1.0, 2.0], +//! [1.5, 2.5], +//! [1.0, 2.0], +//! ]; +//! let x = array![1.0, 2.0, 3.0, 4.0]; +//! +//! let interp = Interp1D::builder(data) +//! .strategy(Linear::new().extrapolate(true)) +//! .x(x) +//! .build().unwrap(); +//! +//! let result = interp.interp(0.5).unwrap(); +//! assert!(result == array![-0.5, 0.5]); +//! let result = interp.interp_array(&array![0.5, 4.0]).unwrap(); +//! assert!(result == array![[-0.5, 0.5], [1.0, 2.0]]); +//! ``` +//! +//! **2D Example** +//! ```rust +//! use ndarray_interp::interp2d::*; +//! use ndarray::*; +//! +//! let data = array![ +//! [1.0, 2.0, 2.5], +//! [3.0, 4.0, 3.5], +//! ]; +//! let interp = Interp2D::builder(data).build().unwrap(); +//! +//! let result = interp.interp_scalar(0.0, 0.5).unwrap(); +//! assert!(result == 1.5); +//! let result = interp.interp_array(&array![0.0, 1.0], &array![0.5, 2.0]).unwrap(); +//! assert!(result == array![1.5, 3.5]); +//! ``` +//! +//! **1D Example with multidimensional data** +//! ``` rust +//! use ndarray_interp::interp2d::*; +//! use ndarray::*; +//! +//! let data = array![ +//! // ---------------------------------> y +//! [[1.0, -1.0], [2.0, -2.0], [3.0, -3.0]], // | +//! [[4.0, -4.0], [5.0, -5.0], [6.0, -6.0]], // | +//! [[7.0, -7.0], [8.0, -8.0], [9.0, -9.0]], // V +//! [[7.5, -7.5], [8.5, -8.5], [9.5, -9.5]], // x +//! ]; +//! let x = array![1.0, 2.0, 3.0, 4.0]; +//! let y = array![1.0, 2.0, 3.0]; +//! +//! let interp = Interp2D::builder(data) +//! .x(x) +//! .y(y) +//! .build().unwrap(); //! +//! let result = interp.interp(1.5, 2.0).unwrap(); +//! assert!(result == array![3.5, -3.5]); +//! let result = interp.interp_array(&array![1.5, 1.5], &array![2.0, 2.5]).unwrap(); +//! assert!(result == array![[3.5, -3.5],[4.0, -4.0]]); +//! ``` use std::mem::ManuallyDrop; @@ -31,12 +132,10 @@ pub enum BuilderError { /// A interpolation axis is not strict monotonic rising #[error("{0}")] Monotonic(String), - /// The lengths of interpolation axis and the - /// corresponding data axis do not match #[error("{0}")] - AxisLenght(String), + ShapeError(String), #[error("{0}")] - DimensionError(String), + ValueError(String), } /// Errors during Interpolation diff --git a/tests/cubic_spline_strat.rs b/tests/cubic_spline_strat.rs index 794f4bd..3a22ff2 100644 --- a/tests/cubic_spline_strat.rs +++ b/tests/cubic_spline_strat.rs @@ -1,52 +1,29 @@ -use approx::assert_abs_diff_eq; -use ndarray::{array, Array1}; -use ndarray_interp::interp1d::{CubicSpline, Interp1D}; +use approx::assert_relative_eq; +use ndarray::{array, stack, Array1, Axis}; +use ndarray_interp::interp1d::cubic_spline::{ + BoundaryCondition, CubicSpline, RowBoundary, SingleBoundary, +}; +use ndarray_interp::interp1d::{Interp1D, Interp1DBuilder}; use ndarray_interp::{BuilderError, InterpolateError}; #[test] -fn interp() { +fn interp_natural() { let data = array![1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0, 0.0, 2.0, 4.0, 6.0, 8.0]; let interp = Interp1D::builder(data) - .strategy(CubicSpline::new()) + .strategy(CubicSpline::new().boundary(BoundaryCondition::Natural)) .build() .unwrap(); let q = Array1::linspace(0.0, 11.0, 30); let res = interp.interp_array(&q).unwrap(); - // values from scipy.interpolate.QubicSpline + // values from scipy.interpolate.QubicSpline with bc_type="natural" let expect = array![ - 1.0, - 1.3917082281418252, - 1.7709152572751259, - 2.125720997885402, - 2.4735200559559645, - 2.873596855334901, - 3.3692218872560726, - 3.822919531969092, - 3.998240261438613, - 3.75923077015136, - 3.279709933108678, - 2.7881342665115523, - 2.390891499049402, - 2.0569231634621636, - 1.744119027809967, - 1.38442936840091, - 0.8991930736934348, - 0.327385578986533, - -0.01567970348252848, - 0.2056442153017282, - 0.9653909358084248, - 1.9164377865351583, - 2.757368677491977, - 3.485961877773172, - 4.197630493134489, - 4.947868508672761, - 5.71920917646552, - 6.487721497405632, - 7.246383891907155, - 8.0 + 1., 1.39170823, 1.77091526, 2.125721, 2.47352006, 2.87359686, 3.36922189, 3.82291953, + 3.99824026, 3.75923077, 3.27970993, 2.78813427, 2.3908915, 2.05692316, 1.74411903, + 1.38442937, 0.89919307, 0.32738558, -0.0156797, 0.20564422, 0.96539094, 1.91643779, + 2.75736868, 3.48596188, 4.19763049, 4.94786851, 5.71920918, 6.4877215, 7.24638389, 8. ]; - assert_abs_diff_eq!(res, expect, epsilon = f64::EPSILON); + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); } #[test] @@ -66,7 +43,7 @@ fn enough_data() { } #[test] -fn extrapolate() { +fn extrapolate_false() { let interp = Interp1D::builder(array![1.0, 2.0, 1.0]) .strategy(CubicSpline::new()) .build() @@ -76,3 +53,557 @@ fn extrapolate() { let err = interp.interp(3.5); assert!(matches!(err, Err(InterpolateError::OutOfBounds(_)))); } + +#[test] +fn extrapolate_natural() { + let data = array![1.0, 2.0, 2.5, 2.5, 3.0, 2.0, 1.0, -2.0, 3.0, 5.0, 6.3, 8.0]; + let interp = Interp1D::builder(data) + .strategy( + CubicSpline::new() + .extrapolate(true) + .boundary(BoundaryCondition::Natural), + ) + .build() + .unwrap(); + let q = Array1::linspace(-3.0, 15.0, 30); + let res = interp.interp_array(&q).unwrap(); + + // values from scipy.interpolate.QubicSpline with bc_type="natural" + let expect = array![ + -0.10117811, + -0.50187696, + -0.46744049, + -0.11138225, + 0.45278419, + 1.11154527, + 1.75138741, + 2.25775994, + 2.49749363, + 2.442418, + 2.62405156, + 3.00988064, + 2.60389947, + 1.96187505, + 1.6459892, + -0.21920517, + -2.0380548, + 0.35839389, + 3.69754559, + 4.82435282, + 5.45047974, + 6.35498498, + 7.39691304, + 8.48312564, + 9.5339106, + 10.46955574, + 11.21034887, + 11.67657779, + 11.78853034, + 11.46649431 + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} + +#[test] +fn extrapolate_not_a_knot() { + let data = array![1f32, 2.0, 2.5, 2.5, 3.0, 2.0, 1.0, -2.0, 3.0, 5.0, 6.3, 8.0]; + let interp = Interp1D::builder(data) + .strategy( + CubicSpline::new() + .extrapolate(true) + .boundary(BoundaryCondition::NotAKnot), + ) + .build() + .unwrap(); + let q = Array1::linspace(-3.0, 15.0, 30); + let res = interp.interp_array(&q).unwrap(); + // values from scipy.interpolate.QubicSpline with bc_type="not-a-knot" + let expect = array![ + 0.94398816f32, + 0.09886458, + -0.16503997, + 0.01013939, + 0.48226752, + 1.109_209_3, + 1.748_829_5, + 2.258_993_1, + 2.497_564_8, + 2.4421474, + 2.624_124, + 3.009_909_2, + 2.603_880_2, + 1.961_875_3, + 1.645_976_8, + -0.21916762, + -2.038_032_5, + 0.35816476, + 3.697_835_4, + 4.825_070_4, + 5.447_815_4, + 6.3556859, + 7.409_040_5, + 8.452_749, + 9.331_685, + 9.890_717, + 9.974_716, + 9.428_555, + 8.097_102, + 5.825_231 + ]; + assert_relative_eq!(res, expect, epsilon = f32::EPSILON, max_relative = 0.001); +} + +#[test] +fn not_a_knot_3_values() { + let interp = Interp1D::builder(array![1.0, 2.0, 0.0]) + .strategy( + CubicSpline::new() + .boundary(BoundaryCondition::NotAKnot) + .extrapolate(true), + ) + .build() + .unwrap(); + + let q = Array1::linspace(-1.0, 3.0, 15); + let res = interp.interp_array(&q).unwrap(); + + let expect = array![ + -3., + -1.55102041, + -0.34693878, + 0.6122449, + 1.32653061, + 1.79591837, + 2.02040816, + 2., + 1.73469388, + 1.2244898, + 0.46938776, + -0.53061224, + -1.7755102, + -3.26530612, + -5. + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} + +#[test] +fn multidim_multi_bounds() { + let y = array![[0.5, 1.0], [0.0, 1.5], [3.0, 0.5],]; + let x = array![-1.0, 0.0, 3.0]; + + // first data column: natural + // second data column top: NotAKnot + // second data column bottom: first derivative == 0.5 + let boundaries = array![[ + RowBoundary::Natural, + RowBoundary::Mixed { + left: SingleBoundary::NotAKnot, + right: SingleBoundary::FirstDeriv(0.5) + } + ],]; + let strat = CubicSpline::new() + .boundary(BoundaryCondition::Individual(boundaries)) + .extrapolate(true); + let interp = Interp1DBuilder::new(y) + .x(x) + .strategy(strat) + .build() + .unwrap(); + + let query = Array1::linspace(-2.0, 4.0, 15); + let res = interp.interp_array(&query).unwrap(); + + let expect = stack![ + Axis(1), + [ + 1., + 0.85787172, + 0.59766764, + 0.30794461, + 0.07725948, + -0.00655977, + 0.10058309, + 0.375, + 0.78717201, + 1.30758017, + 1.90670554, + 2.55502915, + 3.22303207, + 3.88119534, + 4.5 + ], + [ + -1.13194444, + 0.02834467, + 0.81235828, + 1.27749433, + 1.48115079, + 1.48072562, + 1.33361678, + 1.09722222, + 0.82893991, + 0.5861678, + 0.42630385, + 0.40674603, + 0.58489229, + 1.01814059, + 1.76388889, + ] + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} + +#[test] +fn extrapolate_clamped() { + let data = array![1.0, 2.0, 2.5, 2.5, 3.0, 2.0, 1.0, -2.0, 3.0, 5.0, 6.3, 8.0]; + let interp = Interp1D::builder(data) + .strategy( + CubicSpline::new() + .extrapolate(true) + .boundary(BoundaryCondition::Clamped), + ) + .build() + .unwrap(); + let q = Array1::linspace(-3.0, 15.0, 30); + let res = interp.interp_array(&q).unwrap(); + + // values from scipy.interpolate.QubicSpline with bc_type="clamped" + let expect = array![ + 41.28722497, + 23.28738691, + 11.50757146, + 4.70085655, + 1.6203201, + 1.01904002, + 1.65009422, + 2.30659337, + 2.50031574, + 2.43169729, + 2.62693014, + 3.01102652, + 2.60307096, + 1.96191635, + 1.64574608, + -0.21831221, + -2.03751124, + 0.35279783, + 3.70463099, + 4.84190082, + 5.38534268, + 6.37212173, + 7.69341241, + 7.7404559, + 4.5896631, + -3.68255511, + -18.99978784, + -43.28562421, + -78.46365334, + -126.45746433 + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} + +#[test] +fn extrapolate_deriv1() { + let data = array![1.0, 2.0, 2.5, 2.5, 3.0, 2.0, 1.0, -2.0, 3.0, 5.0, 6.3, 8.0]; + let interp = Interp1D::builder(data) + .strategy( + CubicSpline::new() + .extrapolate(true) + .boundary(BoundaryCondition::Individual(array![RowBoundary::Mixed { + left: SingleBoundary::FirstDeriv(-0.1), + right: SingleBoundary::FirstDeriv(-0.5) + },])), + ) + .build() + .unwrap(); + let q = Array1::linspace(-3.0, 15.0, 30); + let res = interp.interp_array(&q).unwrap(); + + // values from scipy.interpolate.QubicSpline with bc_type=((1,-0.1),(1,-0.5)) + let expect = array![ + 45.12263976, + 25.49190916, + 12.61728065, + 5.14680023, + 1.72851392, + 1.01046772, + 1.64070764, + 2.31111841, + 2.50057718, + 2.43070534, + 2.62719459, + 3.01112854, + 2.60301259, + 1.96191065, + 1.64564649, + -0.2180452, + -2.03735486, + 0.35120098, + 3.70664967, + 4.84689904, + 5.36679077, + 6.37700245, + 7.77785832, + 7.52893643, + 3.18149421, + -7.71321086, + -27.60392136, + -58.93937981, + -104.16832878, + -165.7395108 + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} + +#[test] +fn extrapolate_deriv2() { + let data = array![1.0, 2.0, 2.5, 2.5, 3.0, 2.0, 1.0, -2.0, 3.0, 5.0, 6.3, 8.0]; + let interp = Interp1D::builder(data) + .strategy( + CubicSpline::new() + .extrapolate(true) + .boundary(BoundaryCondition::Individual(array![RowBoundary::Mixed { + left: SingleBoundary::SecondDeriv(-0.1), + right: SingleBoundary::SecondDeriv(-0.5) + },])), + ) + .build() + .unwrap(); + let q = Array1::linspace(-3.0, 15.0, 30); + let res = interp.interp_array(&q).unwrap(); + + // values from scipy.interpolate.QubicSpline with bc_type=((2,-0.1),(2,-0.5)) + let expect = array![ + -1.20835424, + -1.1382612, + -0.78778322, + -0.24011435, + 0.42155137, + 1.11401989, + 1.75409718, + 2.25645344, + 2.49741809, + 2.44270565, + 2.62397325, + 3.00984762, + 2.60393207, + 1.96186855, + 1.645952, + -0.21912456, + -2.03800922, + 0.35793208, + 3.69812853, + 4.82579579, + 5.4451242, + 6.35639393, + 7.42129049, + 8.42206522, + 9.12740733, + 9.306006, + 8.72655042, + 7.15772979, + 4.36823329, + 0.12675012 + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} + +#[test] +#[should_panic(expected = "Expected: [1, 2], got: [1, 3]")] +fn bounds_shape_error1() { + let y = array![[0.5, 1.0], [0.0, 1.5], [3.0, 0.5],]; + let boundaries = BoundaryCondition::Individual(array![[ + RowBoundary::Natural, + RowBoundary::Clamped, + RowBoundary::NotAKnot + ],]); + Interp1DBuilder::new(y) + .strategy(CubicSpline::new().boundary(boundaries)) + .build() + .unwrap(); +} + +#[test] +#[should_panic(expected = "Expected: [1, 2], got: [2, 2]")] +fn bounds_shape_error2() { + let y = array![[0.5, 1.0], [0.0, 1.5], [3.0, 0.5],]; + let boundaries = BoundaryCondition::Individual(array![ + [RowBoundary::Natural, RowBoundary::NotAKnot], + [RowBoundary::Natural, RowBoundary::NotAKnot], + ]); + Interp1DBuilder::new(y) + .strategy(CubicSpline::new().boundary(boundaries)) + .build() + .unwrap(); +} + +#[test] +#[should_panic( + expected = "First: [0.5, 1.0], shape=[2], strides=[1], layout=CFcf (0xf), const ndim=1, last: [0.5, 1.1]" +)] +fn periodic_wrong_values() { + let y = array![[0.5, 1.0], [0.0, 1.5], [0.5, 1.1],]; + Interp1DBuilder::new(y) + .strategy(CubicSpline::new().boundary(BoundaryCondition::Periodic)) + .build() + .unwrap(); +} + +#[test] +fn extrapolate_periodic() { + let data = array![1.0, 2.0, 2.5, 2.5, 3.0, 2.0, 1.0, -2.0, 3.0, 5.0, 6.3, 1.0]; + let interp = Interp1D::builder(data) + .strategy( + CubicSpline::new() + .extrapolate(true) + .boundary(BoundaryCondition::Periodic), + ) + .build() + .unwrap(); + + let q = Array1::linspace(-3.0, 15.0, 30); + let res = interp.interp_array(&q).unwrap(); + let expect = array![ + 3., + 4.45171164, + 5.5978812, + 6.54905092, + 3.79486808, + 0.76011398, + 1.36656494, + 2.4432986, + 2.50822019, + 2.40158688, + 2.63514361, + 3.01451693, + 2.59950279, + 1.96267846, + 1.65029582, + -0.22831889, + -2.04318459, + 0.41031552, + 3.63201944, + 4.66215778, + 6.05245899, + 6.19632834, + 2.68818585, + 0.64246067, + 1.77979077, + 2.52789822, + 2.46676892, + 2.41681682, + 2.76866398, + 3. + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} + +#[test] +fn extrapolate_periodic_multidim() { + let y = array![[0.5, 1.0], [0.0, 1.5], [0.0, 1.5], [0.5, 1.0],]; + let x = array![-1.0, 0.0, 2.0, 3.0]; + let interp = Interp1D::builder(y) + .x(x) + .strategy( + CubicSpline::new() + .extrapolate(true) + .boundary(BoundaryCondition::Periodic), + ) + .build() + .unwrap(); + + let q = Array1::linspace(-1.5, 3.5, 15); + let res = interp.interp_array(&q).unwrap(); + let expect = array![ + [0.325, 1.175], + [0.48279883, 1.01720117], + [0.46260933, 1.03739067], + [0.28075802, 1.21924198], + [0.04424198, 1.45575802], + [-0.14693878, 1.64693878], + [-0.26173469, 1.76173469], + [-0.3, 1.8], + [-0.26173469, 1.76173469], + [-0.14693878, 1.64693878], + [0.04424198, 1.45575802], + [0.28075802, 1.21924198], + [0.46260933, 1.03739067], + [0.48279883, 1.01720117], + [0.325, 1.175] + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} + +#[test] +fn extrapolate_periodic_len3() { + let y = array![0.5, 0.0, 0.5]; + let x = array![-1.0, 0.0, 3.0]; + let interp = Interp1D::builder(y) + .x(x) + .strategy( + CubicSpline::new() + .extrapolate(true) + .boundary(BoundaryCondition::Periodic), + ) + .build() + .unwrap(); + + let q = Array1::linspace(-1.5, 3.5, 15); + let res = interp.interp_array(&q).unwrap(); + let expect = array![ + 0.55555556, + 0.53773891, + 0.40889213, + 0.20845481, + 0.02623907, + -0.05701328, + -0.03717201, + 0.05555556, + 0.19080013, + 0.33819242, + 0.46736314, + 0.54794299, + 0.54956268, + 0.44314869, + 0.25 + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} + +#[test] +fn extrapolate_periodic_len3_multidim() { + let y = array![[0.5, 1.0], [0.0, 2.5], [0.5, 1.0],]; + let x = array![-1.0, 0.0, 3.0]; + let interp = Interp1D::builder(y) + .x(x) + .strategy( + CubicSpline::new() + .extrapolate(true) + .boundary(BoundaryCondition::Periodic), + ) + .build() + .unwrap(); + + let q = Array1::linspace(-1.5, 3.5, 15); + let res = interp.interp_array(&q).unwrap(); + let expect = array![ + [0.55555556, 0.83333333], + [0.53773891, 0.88678328], + [0.40889213, 1.27332362], + [0.20845481, 1.87463557], + [0.02623907, 2.4212828], + [-0.05701328, 2.67103984], + [-0.03717201, 2.61151603], + [0.05555556, 2.33333333], + [0.19080013, 1.92759961], + [0.33819242, 1.48542274], + [0.46736314, 1.09791059], + [0.54794299, 0.85617104], + [0.54956268, 0.85131195], + [0.44314869, 1.17055394], + [0.25, 1.75] + ]; + assert_relative_eq!(res, expect, epsilon = f64::EPSILON, max_relative = 0.001); +} diff --git a/tests/interp1d.rs b/tests/interp1d.rs index 888886e..4154e2d 100644 --- a/tests/interp1d.rs +++ b/tests/interp1d.rs @@ -129,7 +129,7 @@ fn interp_builder_errors() { Interp1DBuilder::new(array![1, 2]) .x(array![1, 2, 3]) .build(), - Err(BuilderError::AxisLenght(_)) + Err(BuilderError::ShapeError(_)) )); assert!(matches!( Interp1DBuilder::new(array![1, 2, 3]) diff --git a/tests/interp2d.rs b/tests/interp2d.rs index f7aa34a..4da56e8 100644 --- a/tests/interp2d.rs +++ b/tests/interp2d.rs @@ -290,25 +290,25 @@ fn builder_errors() { Interp2D::builder(array![[1, 2], [3, 4]]) .x(array![1]) .build(), - Err(BuilderError::AxisLenght(_)) + Err(BuilderError::ShapeError(_)) )); assert!(matches!( Interp2D::builder(array![[1, 2], [3, 4]]) .x(array![1, 2, 3]) .build(), - Err(BuilderError::AxisLenght(_)) + Err(BuilderError::ShapeError(_)) )); assert!(matches!( Interp2D::builder(array![[1, 2], [3, 4]]) .y(array![1]) .build(), - Err(BuilderError::AxisLenght(_)) + Err(BuilderError::ShapeError(_)) )); assert!(matches!( Interp2D::builder(array![[1, 2], [3, 4]]) .y(array![1, 2, 3]) .build(), - Err(BuilderError::AxisLenght(_)) + Err(BuilderError::ShapeError(_)) )); assert!(matches!( Interp2D::builder(array![[1, 2], [3, 4]])