rsl_interpolation/
spline.rs

1use crate::Accelerator;
2use crate::DomainError;
3use crate::DynInterpType;
4use crate::InterpType;
5use crate::Interpolation;
6use crate::InterpolationError;
7
8/// 1D Higher level interface.
9///
10/// A Spline owns the data it is constructed with, and provides the same evaluation methods as the
11/// lower-level Interpolator object, without needing to provide the data arrays in every call.
12///
13/// # Example
14/// ```
15/// # use rsl_interpolation::*;
16/// #
17/// # fn main() -> Result<(), InterpolationError>{
18/// let mut acc = Accelerator::new();
19///
20/// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
21/// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
22///
23/// let interp = Cubic.build(&xa, &ya)?;
24///
25/// let typ = Cubic;
26/// let spline = Spline::new(typ, &xa, &ya)?;
27///
28/// let x = 1.5;
29/// let y_interp = interp.eval(&xa, &ya, x, &mut acc)?;
30/// let y_spline = spline.eval(x, &mut acc)?;
31///
32/// assert_eq!(y_interp, y_spline);
33/// #
34/// # Ok(())
35/// # }
36/// ```
37pub struct Spline<I, T>
38where
39    I: InterpType<T> + Send + Sync + 'static,
40{
41    /// The lower-level [`Interpolator`].
42    ///
43    /// [`Interpolator`]: Interpolation#implementors
44    pub interp: I::Interpolation,
45    /// The owned x data.
46    pub xa: Box<[T]>,
47    /// The owned y data.
48    pub ya: Box<[T]>,
49    name: Box<str>,
50    min_size: usize,
51}
52
53impl<I, T> Spline<I, T>
54where
55    I: InterpType<T> + Send + Sync + 'static,
56{
57    /// Constructs a Spline of an Interpolation type `typ` from the data arrays `xa` and `ya`.
58    ///
59    /// # Example
60    /// ```
61    /// # use rsl_interpolation::*;
62    /// #
63    /// # fn main() -> Result<(), InterpolationError>{
64    /// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
65    /// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
66    /// let typ = Cubic;
67    ///
68    /// let spline = Spline::new(typ, &xa, &ya)?;
69    /// #
70    /// # Ok(())
71    /// # }
72    #[doc(alias = "gsl_spline_init")]
73    pub fn new(typ: I, xa: &[T], ya: &[T]) -> Result<Self, InterpolationError>
74    where
75        T: Clone,
76    {
77        Ok(Self {
78            interp: typ.build(xa, ya)?,
79            xa: xa.into(),
80            ya: ya.into(),
81            name: typ.name().into(),
82            min_size: typ.min_size(),
83        })
84    }
85
86    /// Returns the interpolated value `y` for a given point `x`, using the [`Accelerator`] `acc`.
87    ///
88    /// # Example
89    ///
90    /// ```
91    /// # use rsl_interpolation::*;
92    /// #
93    /// # fn main() -> Result<(), InterpolationError>{
94    /// let mut acc = Accelerator::new();
95    ///
96    /// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
97    /// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
98    /// let typ = Cubic;
99    /// let spline = Spline::new(typ, &xa, &ya)?;
100    /// #
101    /// let y = spline.eval(1.5, &mut acc)?;
102    ///
103    /// assert_eq!(y, 3.0);
104    /// # Ok(())
105    /// # }
106    /// ```
107    ///
108    /// # Errors
109    ///
110    /// Returns a [`DomainError`] if `x` is outside the range of `xa`.
111    #[doc(alias = "gsl_spline_eval")]
112    #[doc(alias = "gsl_spline_eval_e")]
113    pub fn eval(&self, x: T, acc: &mut Accelerator) -> Result<T, DomainError> {
114        self.interp.eval(&self.xa, &self.ya, x, acc)
115    }
116
117    /// Returns the derivative `dy/dx` of an interpolated function for a given point `x`, using the
118    /// [`Accelerator`] `acc`.
119    ///
120    /// # Example
121    ///
122    /// ```
123    /// # use rsl_interpolation::*;
124    /// #
125    /// # fn main() -> Result<(), InterpolationError>{
126    /// let mut acc = Accelerator::new();
127    ///
128    /// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
129    /// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
130    /// let typ = Cubic;
131    /// let spline = Spline::new(typ, &xa, &ya)?;
132    ///
133    /// let dydx = spline.eval_deriv(1.5, &mut acc)?;
134    ///
135    /// assert_eq!(dydx, 2.0);
136    /// # Ok(())
137    /// # }
138    /// ```
139    ///
140    /// # Errors
141    ///
142    /// Returns a [`DomainError`] if `x` is outside the range of `xa`.
143    #[doc(alias = "gsl_spline_eval_deriv")]
144    #[doc(alias = "gsl_spline_eval_deriv_e")]
145    pub fn eval_deriv(&self, x: T, acc: &mut Accelerator) -> Result<T, DomainError> {
146        self.interp.eval_deriv(&self.xa, &self.ya, x, acc)
147    }
148
149    /// Returns the second derivative `d²y/dx²` of an interpolated function for a given point `x`, using the
150    /// [`Accelerator`] `acc`.
151    ///
152    /// # Example
153    /// ```
154    /// # use rsl_interpolation::*;
155    /// #
156    /// # fn main() -> Result<(), InterpolationError>{
157    /// let mut acc = Accelerator::new();
158    ///
159    /// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
160    /// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
161    /// let typ = Cubic;
162    /// let spline = Spline::new(typ, &xa, &ya)?;
163    ///
164    /// let dydx = spline.eval_deriv2(1.5, &mut acc)?;
165    ///
166    /// assert_eq!(dydx, 0.0);
167    /// # Ok(())
168    /// # }
169    /// ```
170    ///
171    /// # Errors
172    ///
173    /// Returns a [`DomainError`] if `x` is outside the range of `xa`.
174    #[doc(alias = "gsl_spline_eval_deriv2")]
175    #[doc(alias = "gsl_spline_eval_deriv2_e")]
176    pub fn eval_deriv2(&self, x: T, acc: &mut Accelerator) -> Result<T, DomainError> {
177        self.interp.eval_deriv2(&self.xa, &self.ya, x, acc)
178    }
179
180    #[allow(rustdoc::broken_intra_doc_links)]
181    /// Returns the numerical integral of an interpolated function over the range [`a` ,`b`], using the
182    /// [`Accelerator`] `acc`.
183    ///
184    /// # Example
185    /// ```
186    /// # use rsl_interpolation::*;
187    /// #
188    /// # fn main() -> Result<(), InterpolationError>{
189    /// let mut acc = Accelerator::new();
190    ///
191    /// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
192    /// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
193    /// let typ = Cubic;
194    /// let spline = Spline::new(typ, &xa, &ya)?;
195    ///
196    /// let int = spline.eval_integ(0.0, 2.0, &mut acc)?;
197    ///
198    /// assert_eq!(int, 4.0);
199    /// # Ok(())
200    /// # }
201    /// ```
202    ///
203    /// # Errors
204    ///
205    /// Returns a [`DomainError`] if `a` or `b` is outside the range of xa.
206    #[doc(alias = "gsl_spline_eval_integ")]
207    #[doc(alias = "gsl_spline_eval_integ_e")]
208    pub fn eval_integ(&self, a: T, b: T, acc: &mut Accelerator) -> Result<T, DomainError> {
209        self.interp.eval_integ(&self.xa, &self.ya, a, b, acc)
210    }
211
212    /// Returns the name of the Interpolator.
213    #[doc(alias = "gsl_spline_name")]
214    pub fn name(&self) -> &str {
215        &self.name
216    }
217
218    /// Returns the minimum number of points required by the Interpolator.
219    #[doc(alias = "gsl_spline_min_size")]
220    pub fn min_size(&self) -> usize {
221        self.min_size
222    }
223}
224
225/// 2D Spline with runtime-determined Interpolation Type.
226pub type DynSpline<T> = Spline<DynInterpType<T>, T>;
227
228impl<T> DynSpline<T> {
229    /// Constructs a Spline of a dynamic Interpolation type `typ` from the data arrays `xa` and
230    /// `ya`.
231    ///
232    /// # Example
233    /// ```
234    /// # use rsl_interpolation::*;
235    /// #
236    /// # fn main() -> Result<(), InterpolationError> {
237    /// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
238    /// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
239    /// let typ = "cubic";
240    ///
241    /// let spline = match typ {
242    ///     "cubic" => Spline::new_dyn(Cubic, &xa, &ya)?,
243    ///     "akima" => Spline::new_dyn(Akima, &xa, &ya)?,
244    ///     // ...
245    ///     _ => unreachable!()
246    /// };
247    /// # Ok(())
248    /// # }
249    /// ```
250    #[doc(alias = "gsl_spline_init")]
251    pub fn new_dyn<I>(typ: I, xa: &[T], ya: &[T]) -> Result<Self, InterpolationError>
252    where
253        T: Clone,
254        I: InterpType<T> + Send + Sync + 'static,
255        I::Interpolation: Send + Sync + 'static,
256    {
257        Self::new(DynInterpType::new(typ), xa, ya)
258    }
259}
260
261/// Creates a [`DynSpline`] of `typ` type.
262///
263/// Useful when `typ` is not known at compile time.
264///
265/// # Example
266/// ```
267/// # use rsl_interpolation::*;
268/// #
269/// # fn main() -> Result<(), InterpolationError> {
270/// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
271/// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
272/// let typ = "cubic";
273///
274/// let spline = make_spline(typ, &xa, &ya)?;
275/// # Ok(())
276/// # }
277/// ```
278pub fn make_spline<T>(typ: &str, xa: &[T], ya: &[T]) -> Result<DynSpline<T>, InterpolationError>
279where
280    T: crate::Num + ndarray_linalg::Lapack,
281{
282    use crate::*;
283
284    match typ.to_lowercase().as_str() {
285        "linear" => Ok(DynSpline::new_dyn(Linear, xa, ya)?),
286        "cubic" => Ok(DynSpline::new_dyn(Cubic, xa, ya)?),
287        "akima" => Ok(DynSpline::new_dyn(Akima, xa, ya)?),
288        "cubicperiodic" | "cubic periodic" => Ok(DynSpline::new_dyn(CubicPeriodic, xa, ya)?),
289        "akimaperiodic" | "akima periodic" => Ok(DynSpline::new_dyn(AkimaPeriodic, xa, ya)?),
290        "steffen" => Ok(DynSpline::new_dyn(Steffen, xa, ya)?),
291        _ => Err(InterpolationError::InvalidType(typ.into())),
292    }
293}
294
295#[cfg(test)]
296mod test {
297    use super::*;
298    use crate::*;
299
300    #[test]
301    fn test_spline_creation() {
302        let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
303        let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
304
305        let spline = Spline::new(Cubic, &xa, &ya).unwrap();
306        let _: &str = spline.name();
307        let _: usize = spline.min_size();
308    }
309
310    #[test]
311    fn test_spline_eval() {
312        let xa = [0.0, 1.0, 2.0];
313        let ya = [0.0, 1.0, 2.0];
314        let spline = Spline::new(Cubic, &xa, &ya).unwrap();
315        let mut acc = Accelerator::new();
316
317        let x = 0.5;
318        let y = spline.eval(x, &mut acc).unwrap();
319        let dy = spline.eval_deriv(x, &mut acc).unwrap();
320        let dy2 = spline.eval_deriv2(x, &mut acc).unwrap();
321        let int = spline.eval_integ(0.0, x, &mut acc).unwrap();
322
323        assert_eq!(y, 0.5);
324        assert_eq!(dy, 1.0);
325        assert_eq!(dy2, 0.0);
326        assert_eq!(int, 0.125);
327    }
328
329    #[test]
330    fn test_dyn_spline() {
331        let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
332        let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
333
334        Spline::new(Cubic, &xa, &ya).unwrap();
335        Spline::new_dyn(Cubic, &xa, &ya).unwrap();
336    }
337}
rsl_interpolation/spline.rs

rsl_interpolation/
spline.rs
