rateslib/splines/
spline.rs

1use crate::dual::linalg::{dmul11_, fdmul11_, fdsolve, fouter11_};
2use crate::dual::{Dual, Dual2, Gradient1, Gradient2, Number, NumberMapping};
3use ndarray::{Array1, Array2};
4use num_traits::{Signed, Zero};
5use pyo3::exceptions::{PyTypeError, PyValueError};
6use pyo3::{pyclass, PyErr};
7use serde::{Deserialize, Serialize};
8use std::{
9    cmp::PartialEq,
10    iter::{zip, Sum},
11    ops::{Mul, Sub},
12};
13
14/// Evaluate the `x` value on the `i`'th B-spline with order `k` and knot sequence `t`.
15///
16/// Note `org_k` should be input as None, it is used internally for recursively calculating
17/// spline derivatives, where it is set to the original `k` value from the outer scope.
18pub fn bsplev_single_f64(x: &f64, i: usize, k: &usize, t: &Vec<f64>, org_k: Option<usize>) -> f64 {
19    let org_k: usize = org_k.unwrap_or(*k);
20
21    // Short circuit (positivity and support property)
22    if *x < t[i] || *x > t[i + k] {
23        return 0.0_f64;
24    }
25
26    // Right side end point support
27    if *x == t[t.len() - 1] && i >= (t.len() - org_k - 1) {
28        return 1.0_f64;
29    }
30
31    // Recursion
32    if *k == 1_usize {
33        if t[i] <= *x && *x < t[i + 1] {
34            1.0_f64
35        } else {
36            0.0_f64
37        }
38    } else {
39        let mut left: f64 = 0.0_f64;
40        let mut right: f64 = 0.0_f64;
41        if t[i] != t[i + k - 1] {
42            left = (x - t[i]) / (t[i + k - 1] - t[i]) * bsplev_single_f64(x, i, &(k - 1), t, None);
43        }
44        if t[i + 1] != t[i + k] {
45            right = (t[i + k] - x) / (t[i + k] - t[i + 1])
46                * bsplev_single_f64(x, i + 1, &(k - 1), t, None);
47        }
48        left + right
49    }
50}
51
52/// Evaluate the `x` value on the `i`'th B-spline with order `k` and knot sequence `t`.
53///
54/// Note `org_k` should be input as None, it is used internally for recursively calculating
55/// spline derivatives, where it is set to the original `k` value from the outer scope.
56pub fn bsplev_single_dual(
57    x: &Dual,
58    i: usize,
59    k: &usize,
60    t: &Vec<f64>,
61    org_k: Option<usize>,
62) -> Dual {
63    let b_f64 = bsplev_single_f64(&x.real(), i, k, t, org_k);
64    let dbdx_f64 = bspldnev_single_f64(&x.real(), i, k, t, 1, org_k);
65    Dual::clone_from(x, b_f64, dbdx_f64 * x.dual())
66}
67
68/// Evaluate the `x` value on the `i`'th B-spline with order `k` and knot sequence `t`.
69///
70/// Note `org_k` should be input as None, it is used internally for recursively calculating
71/// spline derivatives, where it is set to the original `k` value from the outer scope.
72pub fn bsplev_single_dual2(
73    x: &Dual2,
74    i: usize,
75    k: &usize,
76    t: &Vec<f64>,
77    org_k: Option<usize>,
78) -> Dual2 {
79    let b_f64 = bsplev_single_f64(&x.real(), i, k, t, org_k);
80    let dbdx_f64 = bspldnev_single_f64(&x.real(), i, k, t, 1, org_k);
81    let d2bdx2_f64 = bspldnev_single_f64(&x.real(), i, k, t, 2, org_k);
82    let dual2 =
83        dbdx_f64 * x.dual2() + 0.5 * d2bdx2_f64 * fouter11_(&x.dual().view(), &x.dual().view());
84    Dual2::clone_from(x, b_f64, dbdx_f64 * x.dual(), dual2)
85}
86
87/// Evaluate the `m`'th order derivative of the `x` value on the `i`'th B-spline with
88/// order `k` and knot sequence `t`.
89///
90/// Note `org_k` should be input as None, it is used internally for recursively calculating
91/// spline derivatives, where it is set to the original `k` value from the outer scope.
92pub fn bspldnev_single_f64(
93    x: &f64,
94    i: usize,
95    k: &usize,
96    t: &Vec<f64>,
97    m: usize,
98    org_k: Option<usize>,
99) -> f64 {
100    if m == 0 {
101        return bsplev_single_f64(x, i, k, t, None);
102    } else if *k == 1 || m >= *k {
103        return 0.0_f64;
104    }
105
106    let org_k: usize = org_k.unwrap_or(*k);
107    let mut r: f64 = 0.0;
108    let div1: f64 = t[i + k - 1] - t[i];
109    let div2: f64 = t[i + k] - t[i + 1];
110
111    if m == 1 {
112        if div1 != 0_f64 {
113            r += bsplev_single_f64(x, i, &(k - 1), t, Some(org_k)) / div1;
114        }
115        if div2 != 0_f64 {
116            r -= bsplev_single_f64(x, i + 1, &(k - 1), t, Some(org_k)) / div2;
117        }
118        r *= (k - 1) as f64;
119    } else {
120        if div1 != 0_f64 {
121            r += bspldnev_single_f64(x, i, &(k - 1), t, m - 1, Some(org_k)) / div1;
122        }
123        if div2 != 0_f64 {
124            r -= bspldnev_single_f64(x, i + 1, &(k - 1), t, m - 1, Some(org_k)) / div2;
125        }
126        r *= (k - 1) as f64
127    }
128    r
129}
130
131/// Evaluate the `m`'th order derivative of the `x` value on the `i`'th B-spline with
132/// order `k` and knot sequence `t`.
133///
134/// Note `org_k` should be input as None, it is used internally for recursively calculating
135/// spline derivatives, where it is set to the original `k` value from the outer scope.
136pub fn bspldnev_single_dual(
137    x: &Dual,
138    i: usize,
139    k: &usize,
140    t: &Vec<f64>,
141    m: usize,
142    org_k: Option<usize>,
143) -> Dual {
144    let b_f64 = bspldnev_single_f64(&x.real(), i, k, t, m, org_k);
145    let dbdx_f64 = bspldnev_single_f64(&x.real(), i, k, t, m + 1, org_k);
146    Dual::clone_from(x, b_f64, dbdx_f64 * x.dual())
147}
148
149/// Evaluate the `m`'th order derivative of the `x` value on the `i`'th B-spline with
150/// order `k` and knot sequence `t`.
151///
152/// Note `org_k` should be input as None, it is used internally for recursively calculating
153/// spline derivatives, where it is set to the original `k` value from the outer scope.
154pub fn bspldnev_single_dual2(
155    x: &Dual2,
156    i: usize,
157    k: &usize,
158    t: &Vec<f64>,
159    m: usize,
160    org_k: Option<usize>,
161) -> Dual2 {
162    let b_f64 = bspldnev_single_f64(&x.real(), i, k, t, m, org_k);
163    let dbdx_f64 = bspldnev_single_f64(&x.real(), i, k, t, m + 1, org_k);
164    let d2bdx2_f64 = bspldnev_single_f64(&x.real(), i, k, t, m + 2, org_k);
165    let dual2 =
166        dbdx_f64 * x.dual2() + 0.5 * d2bdx2_f64 * fouter11_(&x.dual().view(), &x.dual().view());
167    Dual2::clone_from(x, b_f64, dbdx_f64 * x.dual(), dual2)
168}
169
170/// A piecewise polynomial spline of given order and knot sequence.
171#[derive(Clone, Debug, Deserialize, Serialize)]
172pub struct PPSpline<T> {
173    k: usize,
174    t: Vec<f64>,
175    c: Option<Array1<T>>,
176    n: usize,
177}
178
179impl<T> PPSpline<T> {
180    pub fn k(&self) -> &usize {
181        &self.k
182    }
183
184    pub fn t(&self) -> &Vec<f64> {
185        &self.t
186    }
187
188    pub fn n(&self) -> &usize {
189        &self.n
190    }
191
192    pub fn c(&self) -> &Option<Array1<T>> {
193        &self.c
194    }
195}
196
197impl<T> PPSpline<T>
198where
199    T: PartialOrd + Signed + Clone + Sum + Zero,
200    for<'a> &'a T: Sub<&'a T, Output = T>,
201    for<'a> &'a f64: Mul<&'a T, Output = T>,
202{
203    /// Create a PPSpline from its order `k`, knot sequence `t` and optional spline coefficents `c`.
204    pub fn new(k: usize, t: Vec<f64>, c: Option<Vec<T>>) -> Self {
205        // t is given and is non-decreasing
206        assert!(t.len() > 1);
207        assert!(zip(&t[1..], &t[..(t.len() - 1)]).all(|(a, b)| a >= b));
208        let n = t.len() - k;
209        let c_ = c.map(Array1::from_vec);
210        PPSpline { k, t, n, c: c_ }
211    }
212
213    pub fn ppdnev_single(&self, x: &f64, m: usize) -> Result<T, PyErr> {
214        let b: Array1<f64> = Array1::from_vec(
215            (0..self.n)
216                .map(|i| bspldnev_single_f64(x, i, &self.k, &self.t, m, None))
217                .collect(),
218        );
219        match &self.c {
220            Some(c) => Ok(fdmul11_(&b.view(), &c.view())),
221            None => Err(PyValueError::new_err(
222                "Must call `csolve` before evaluating PPSpline.",
223            )),
224        }
225    }
226
227    pub fn csolve(
228        &mut self,
229        tau: &[f64],
230        y: &[T],
231        left_n: usize,
232        right_n: usize,
233        allow_lsq: bool,
234    ) -> Result<(), PyErr> {
235        if tau.len() != self.n && !(allow_lsq && tau.len() > self.n) {
236            return Err(PyValueError::new_err(
237                "`csolve` cannot complete if length of `tau` < n or `allow_lsq` is false.",
238            ));
239        }
240        if tau.len() != y.len() {
241            return Err(PyValueError::new_err(
242                "`tau` and `y` must have the same length.",
243            ));
244        }
245        let b: Array2<f64> = self.bsplmatrix(tau, left_n, right_n);
246        let ya: Array1<T> = Array1::from_vec(y.to_owned());
247        let c: Array1<T> = fdsolve(&b.view(), &ya.view(), allow_lsq);
248        self.c = Some(c);
249        Ok(())
250    }
251
252    // pub fn bsplev(&self, x: &Vec<f64>, i: &usize) -> Vec<f64> {
253    //     x.iter().map(|v| bsplev_single_f64(v, *i, self.k(), self.t(), None)).collect()
254    // }
255
256    pub fn bspldnev(&self, x: &[f64], i: &usize, m: &usize) -> Vec<f64> {
257        x.iter()
258            .map(|v| bspldnev_single_f64(v, *i, self.k(), self.t(), *m, None))
259            .collect()
260    }
261
262    pub fn bsplmatrix(&self, tau: &[f64], left_n: usize, right_n: usize) -> Array2<f64> {
263        let mut b = Array2::zeros((tau.len(), self.n));
264        for i in 0..self.n {
265            b[[0, i]] = bspldnev_single_f64(&tau[0], i, &self.k, &self.t, left_n, None);
266            b[[tau.len() - 1, i]] =
267                bspldnev_single_f64(&tau[tau.len() - 1], i, &self.k, &self.t, right_n, None);
268            for j in 1..(tau.len() - 1) {
269                b[[j, i]] = bsplev_single_f64(&tau[j], i, &self.k, &self.t, None)
270            }
271        }
272        b
273    }
274}
275
276impl NumberMapping for PPSpline<f64> {
277    fn mapped_value(&self, x: &Number) -> Result<Number, PyErr> {
278        match x {
279            Number::F64(f) => Ok(Number::F64(self.ppdnev_single(f, 0_usize)?)),
280            Number::Dual(d) => Ok(Number::Dual(self.ppdnev_single_dual(d, 0_usize)?)),
281            Number::Dual2(d) => Ok(Number::Dual2(self.ppdnev_single_dual2(d, 0_usize)?)),
282        }
283    }
284}
285
286impl PPSpline<f64> {
287    pub fn ppdnev_single_dual(&self, x: &Dual, m: usize) -> Result<Dual, PyErr> {
288        let b: Array1<Dual> = Array1::from_vec(
289            (0..self.n)
290                .map(|i| bspldnev_single_dual(x, i, &self.k, &self.t, m, None))
291                .collect(),
292        );
293        match &self.c {
294            Some(c) => Ok(fdmul11_(&c.view(), &b.view())),
295            None => Err(PyValueError::new_err(
296                "Must call `csolve` before evaluating PPSpline.",
297            )),
298        }
299    }
300
301    pub fn ppdnev_single_dual2(&self, x: &Dual2, m: usize) -> Result<Dual2, PyErr> {
302        let b: Array1<Dual2> = Array1::from_vec(
303            (0..self.n)
304                .map(|i| bspldnev_single_dual2(x, i, &self.k, &self.t, m, None))
305                .collect(),
306        );
307        match &self.c {
308            Some(c) => Ok(fdmul11_(&c.view(), &b.view())),
309            None => Err(PyValueError::new_err(
310                "Must call `csolve` before evaluating PPSpline.",
311            )),
312        }
313    }
314}
315
316impl NumberMapping for PPSpline<Dual> {
317    fn mapped_value(&self, x: &Number) -> Result<Number, PyErr> {
318        match x {
319            Number::F64(f) => Ok(Number::Dual(self.ppdnev_single(f, 0_usize)?)),
320            Number::Dual(d) => Ok(Number::Dual(self.ppdnev_single_dual(d, 0_usize)?)),
321            Number::Dual2(d) => Ok(Number::Dual2(self.ppdnev_single_dual2(d, 0_usize)?)),
322        }
323    }
324}
325
326impl PPSpline<Dual> {
327    pub fn ppdnev_single_dual2(&self, _x: &Dual2, _m: usize) -> Result<Dual2, PyErr> {
328        Err(PyTypeError::new_err(
329            "Cannot index with type `Dual2` on PPSpline<Dual>`.",
330        ))
331    }
332
333    pub fn ppdnev_single_dual(&self, x: &Dual, m: usize) -> Result<Dual, PyErr> {
334        let b: Array1<Dual> = Array1::from_vec(
335            (0..self.n)
336                .map(|i| bspldnev_single_dual(x, i, &self.k, &self.t, m, None))
337                .collect(),
338        );
339        match &self.c {
340            Some(c) => Ok(dmul11_(&c.view(), &b.view())),
341            None => Err(PyValueError::new_err(
342                "Must call `csolve` before evaluating PPSpline.",
343            )),
344        }
345    }
346}
347
348impl NumberMapping for PPSpline<Dual2> {
349    fn mapped_value(&self, x: &Number) -> Result<Number, PyErr> {
350        match x {
351            Number::F64(f) => Ok(Number::Dual2(self.ppdnev_single(f, 0_usize)?)),
352            Number::Dual(d) => Ok(Number::Dual(self.ppdnev_single_dual(d, 0_usize)?)),
353            Number::Dual2(d) => Ok(Number::Dual2(self.ppdnev_single_dual2(d, 0_usize)?)),
354        }
355    }
356}
357
358impl PPSpline<Dual2> {
359    pub fn ppdnev_single_dual(&self, _x: &Dual, _m: usize) -> Result<Dual, PyErr> {
360        Err(PyTypeError::new_err(
361            "Cannot index with type `Dual` on PPSpline<Dual2>.",
362        ))
363    }
364
365    pub fn ppdnev_single_dual2(&self, x: &Dual2, m: usize) -> Result<Dual2, PyErr> {
366        let b: Array1<Dual2> = Array1::from_vec(
367            (0..self.n)
368                .map(|i| bspldnev_single_dual2(x, i, &self.k, &self.t, m, None))
369                .collect(),
370        );
371        match &self.c {
372            Some(c) => Ok(dmul11_(&c.view(), &b.view())),
373            None => Err(PyValueError::new_err(
374                "Must call `csolve` before evaluating PPSpline.",
375            )),
376        }
377    }
378}
379
380impl<T> PartialEq for PPSpline<T>
381where
382    T: PartialEq,
383{
384    /// Equality of `PPSpline` if
385
386    fn eq(&self, other: &Self) -> bool {
387        if self.k != other.k || self.n != other.n {
388            return false;
389        }
390        if !self.t.eq(&other.t) {
391            return false;
392        }
393        match (&self.c, &other.c) {
394            (Some(c1), Some(c2)) => c1.eq(&c2),
395            (Some(_c), None) => false,
396            (None, Some(_c)) => false,
397            (None, None) => true,
398        }
399    }
400}
401
402/// Definitive [f64] type variant of a [PPSpline].
403#[pyclass(module = "rateslib.rs")]
404#[derive(Clone, Deserialize, Serialize)]
405pub struct PPSplineF64 {
406    pub(crate) inner: PPSpline<f64>,
407}
408
409/// Definitive [Dual] type variant of a [PPSpline].
410#[pyclass(module = "rateslib.rs")]
411#[derive(Clone, Deserialize, Serialize)]
412pub struct PPSplineDual {
413    pub(crate) inner: PPSpline<Dual>,
414}
415
416/// Definitive [Dual2] type variant of a [PPSpline].
417#[pyclass(module = "rateslib.rs")]
418#[derive(Clone, Deserialize, Serialize)]
419pub struct PPSplineDual2 {
420    pub(crate) inner: PPSpline<Dual2>,
421}
422
423impl PartialEq for PPSplineF64 {
424    /// Equality of `PPSplineF64` if
425    fn eq(&self, other: &Self) -> bool {
426        self.inner.eq(&other.inner)
427    }
428}
429
430impl PartialEq for PPSplineDual {
431    /// Equality of `PPSplineDual` if
432    fn eq(&self, other: &Self) -> bool {
433        self.inner.eq(&other.inner)
434    }
435}
436
437impl PartialEq for PPSplineDual2 {
438    /// Equality of `PPSplineDual2` if
439    fn eq(&self, other: &Self) -> bool {
440        self.inner.eq(&other.inner)
441    }
442}
443
444// UNIT TESTS
445
446//
447
448//
449
450#[cfg(test)]
451mod tests {
452    use super::*;
453    use crate::dual::Dual;
454    use ndarray::{arr1, arr2};
455    use num_traits::One;
456
457    fn is_close(a: &f64, b: &f64, abs_tol: Option<f64>) -> bool {
458        // used rather than equality for float numbers
459        (a - b).abs() < abs_tol.unwrap_or(1e-8)
460    }
461
462    #[test]
463    fn bsplev_single_f64_() {
464        let x: f64 = 1.5_f64;
465        let t: Vec<f64> = Vec::from(&[1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.]);
466        let k: usize = 4;
467        let result: Vec<f64> = (0..8)
468            .map(|i| bsplev_single_f64(&x, i as usize, &k, &t, None))
469            .collect();
470        let expected: Vec<f64> = Vec::from(&[0.125, 0.375, 0.375, 0.125, 0., 0., 0., 0.]);
471        assert_eq!(result, expected)
472    }
473
474    #[test]
475    fn bsplev_single_dual_() {
476        let x: Dual = Dual::new(1.5, vec!["x".to_string()]);
477        let t: Vec<f64> = Vec::from(&[1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.]);
478        let k: usize = 4;
479        let result: Vec<Dual> = (0..8)
480            .map(|i| bsplev_single_dual(&x, i as usize, &k, &t, None))
481            .collect();
482        let expected: Vec<f64> = Vec::from(&[0.125, 0.375, 0.375, 0.125, 0., 0., 0., 0.]);
483        for i in 0..8 {
484            assert_eq!(result[i].real(), expected[i])
485        }
486        // These are values from the bspldnev_single evaluation test
487        let dual_expected: Vec<f64> = Vec::from(&[-0.75, -0.75, 0.75, 0.75, 0., 0., 0., 0.]);
488        for i in 0..8 {
489            assert_eq!(result[i].dual()[0], dual_expected[i])
490        }
491    }
492
493    #[test]
494    fn bsplev_single_f64_right() {
495        let x: f64 = 4.0_f64;
496        let t: Vec<f64> = Vec::from(&[1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.]);
497        let k: usize = 4;
498        let result: Vec<f64> = (0..8)
499            .map(|i| bsplev_single_f64(&x, i as usize, &k, &t, None))
500            .collect();
501        let expected: Vec<f64> = Vec::from(&[0., 0., 0., 0., 0., 0., 0., 1.0]);
502        assert_eq!(result, expected)
503    }
504
505    #[test]
506    fn bspldnev_single_f64_() {
507        let x: f64 = 1.5_f64;
508        let t: Vec<f64> = Vec::from(&[1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.]);
509        let k: usize = 4;
510        let result: Vec<f64> = (0..8)
511            .map(|i| bspldnev_single_f64(&x, i as usize, &k, &t, 1_usize, None))
512            .collect();
513        let expected: Vec<f64> = Vec::from(&[-0.75, -0.75, 0.75, 0.75, 0., 0., 0., 0.]);
514        assert_eq!(result, expected)
515    }
516
517    #[test]
518    fn bspldnev_single_f64_right() {
519        let x: f64 = 4.0_f64;
520        let t: Vec<f64> = Vec::from(&[1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.]);
521        let k: usize = 4;
522        let result: Vec<f64> = (0..8)
523            .map(|i| bspldnev_single_f64(&x, i as usize, &k, &t, 1_usize, None))
524            .collect();
525        let expected: Vec<f64> = Vec::from(&[0., 0., 0., 0., 0., 0., -3., 3.]);
526        assert_eq!(result, expected)
527    }
528
529    #[test]
530    fn bspldnev_single_shortcut() {
531        let x: f64 = 1.5_f64;
532        let t: Vec<f64> = Vec::from(&[1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.]);
533        let k: usize = 4;
534        let result: Vec<f64> = (0..8)
535            .map(|i| bspldnev_single_f64(&x, i as usize, &k, &t, 6_usize, None))
536            .collect();
537        let expected: Vec<f64> = Vec::from(&[0., 0., 0., 0., 0., 0., 0., 0.]);
538        assert_eq!(result, expected)
539    }
540
541    #[test]
542    fn bspldnev_single_m() {
543        let x: f64 = 4.0_f64;
544        let t: Vec<f64> = Vec::from(&[1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.]);
545        let k: usize = 4;
546        let result: Vec<f64> = (0..8)
547            .map(|i| bspldnev_single_f64(&x, i as usize, &k, &t, 2_usize, None))
548            .collect();
549        let expected: Vec<f64> = Vec::from(&[0., 0., 0., 0., 0., 3., -9., 6.]);
550        assert_eq!(result, expected)
551    }
552
553    #[test]
554    fn bspldnev_single_m_zero() {
555        let x: f64 = 1.5_f64;
556        let t: Vec<f64> = Vec::from(&[1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.]);
557        let k: usize = 4;
558        let result: Vec<f64> = (0..8)
559            .map(|i| bspldnev_single_f64(&x, i as usize, &k, &t, 0_usize, None))
560            .collect();
561        let expected: Vec<f64> = Vec::from(&[0.125, 0.375, 0.375, 0.125, 0., 0., 0., 0.]);
562        assert_eq!(result, expected)
563    }
564
565    #[test]
566    fn ppspline_new() {
567        let _pps: PPSpline<f64> = PPSpline::new(
568            4,
569            vec![1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.],
570            None,
571        );
572    }
573
574    #[test]
575    fn ppspline_bsplmatrix() {
576        let pps: PPSpline<f64> = PPSpline::new(4, vec![1., 1., 1., 1., 2., 3., 3., 3., 3.], None);
577        let result = pps.bsplmatrix(&vec![1., 1., 2., 3., 3.], 2_usize, 2_usize);
578        let expected: Array2<f64> = arr2(&[
579            [6., -9., 3., 0., 0.],
580            [1., 0., 0., 0., 0.],
581            [0., 0.25, 0.5, 0.25, 0.],
582            [0., 0., 0., 0., 1.],
583            [0., 0., 3., -9., 6.],
584        ]);
585        assert_eq!(result, expected)
586    }
587
588    #[test]
589    fn csolve_() {
590        let t = vec![0., 0., 0., 0., 4., 4., 4., 4.];
591        let tau = vec![0., 1., 3., 4.];
592        let val = vec![0., 0., 2., 2.];
593        let mut pps: PPSpline<f64> = PPSpline::new(4, t, None);
594        let _ = pps.csolve(&tau, &val, 0, 0, false);
595        let expected = vec![0., -1.11111111, 3.111111111111, 2.0];
596        let v: Vec<bool> = pps
597            .c
598            .expect("csolve")
599            .into_raw_vec_and_offset()
600            .0
601            .iter()
602            .zip(expected.iter())
603            .map(|(x, y)| is_close(&x, &y, None))
604            .collect();
605
606        assert!(v.iter().all(|x| *x));
607    }
608
609    #[test]
610    fn csolve_dual() {
611        let t = vec![0., 0., 0., 0., 4., 4., 4., 4.];
612        let tau = vec![0., 1., 3., 4.];
613        let d1 = Dual::one();
614        let val = vec![0. * &d1, 0. * &d1, 2. * &d1, 2. * &d1];
615        let mut pps = PPSpline::new(4, t, None);
616        let _ = pps.csolve(&tau, &val, 0, 0, false);
617        let expected = vec![0. * &d1, -1.11111111 * &d1, 3.111111111111 * &d1, 2.0 * &d1];
618        let v: Vec<bool> = pps
619            .c
620            .expect("csolve")
621            .into_raw_vec_and_offset()
622            .0
623            .iter()
624            .zip(expected.iter())
625            .map(|(x, y)| is_close(&x.real(), &y.real(), None))
626            .collect();
627
628        assert!(v.iter().all(|x| *x));
629    }
630
631    #[test]
632    fn ppev_single_() {
633        let t = vec![1., 1., 1., 1., 2., 2., 2., 3., 4., 4., 4., 4.];
634        let mut pps = PPSpline::new(4, t, None);
635        pps.c = Some(arr1(&[1., 2., -1., 2., 1., 1., 2., 2.]));
636        let r1 = pps.ppdnev_single(&1.1, 0).unwrap();
637        assert!(is_close(&r1, &1.19, None));
638        let r2 = pps.ppdnev_single(&1.8, 0).unwrap();
639        assert!(is_close(&r2, &0.84, None));
640        let r3 = pps.ppdnev_single(&2.8, 0).unwrap();
641        assert!(is_close(&r3, &1.136, None));
642    }
643
644    #[test]
645    fn partialeq_() {
646        let pp1 = PPSpline::<f64>::new(2, vec![1., 1., 2., 2.], None);
647        let pp2 = PPSpline::<f64>::new(2, vec![1., 1., 2., 2.], None);
648        assert!(pp1 == pp2);
649        let pp3 = PPSpline::new(2, vec![1., 1., 2., 2.], Some(vec![1.5, 0.2]));
650        let pp4 = PPSpline::new(2, vec![1., 1., 2., 2.], Some(vec![1.5, 0.2]));
651        assert!(pp3 == pp4);
652        assert!(pp3 != pp2);
653        assert!(pp1 != pp4);
654    }
655
656    #[test]
657    #[should_panic]
658    fn backwards_definition() {
659        let _pp1 = PPSpline::<f64>::new(4, vec![3., 3., 3., 3., 2., 1., 1., 1., 1.], None);
660    }
661}
rateslib/splines/spline.rs

rateslib/splines/
spline.rs
