rateslib/splines/spline_py.rs
1//! Wrapper to export spline functionality to Python
2
3use crate::dual::{Dual, Dual2, Number};
4use crate::json::json_py::DeserializedObj;
5use crate::json::JSON;
6use crate::splines::spline::{
7 bspldnev_single_f64, bsplev_single_f64, PPSpline, PPSplineDual, PPSplineDual2, PPSplineF64,
8};
9use std::cmp::PartialEq;
10
11use numpy::{PyArray2, ToPyArray};
12use pyo3::exceptions::{PyTypeError, PyValueError};
13use pyo3::prelude::*;
14
15macro_rules! create_interface {
16 ($name: ident, $type: ident) => {
17 #[pymethods]
18 impl $name {
19 #[new]
20 #[pyo3(signature = (k, t, c=None))]
21 fn new(k: usize, t: Vec<f64>, c: Option<Vec<$type>>) -> Self {
22 Self {
23 inner: PPSpline::new(k, t, c),
24 }
25 }
26
27 fn __getnewargs__(&self, _py: Python) -> PyResult<(usize, Vec<f64>, Option<Vec<$type>>)> {
28 Ok((self.k()?, self.t()?, self.c()?))
29 }
30
31 /// The dimension of the pp spline.
32 #[getter]
33 fn n(&self) -> PyResult<usize> {
34 Ok(*self.inner.n())
35 }
36
37 /// The order of the pp spline.
38 #[getter]
39 fn k(&self) -> PyResult<usize> {
40 Ok(*self.inner.k())
41 }
42
43 /// The knot sequence of the pp spline, of length ``n+k``.
44 #[getter]
45 fn t(&self) -> PyResult<Vec<f64>> {
46 Ok(self.inner.t().clone())
47 }
48
49 /// The spline coefficients of length ``n``.
50 #[getter]
51 fn c(&self) -> PyResult<Option<Vec<$type>>> {
52 match self.inner.c() {
53 Some(val) => Ok(Some(val.clone().into_raw_vec_and_offset().0)),
54 None => Ok(None)
55 }
56 }
57
58 /// Solve the spline coefficients given the data sites.
59 ///
60 /// Parameters
61 /// ----------
62 /// tau: list[f64]
63 /// The data site `x`-coordinates.
64 /// y: list[type]
65 /// The data site `y`-coordinates in appropriate type (float, *Dual* or *Dual2*)
66 /// for *self*.
67 /// left_n: int
68 /// The number of derivatives to evaluate at the left side of the data sites,
69 /// i.e. defining an endpoint constraint.
70 /// right_n: int
71 /// The number of derivatives to evaluate at the right side of the datasites,
72 /// i.e. defining an endpoint constraint.
73 /// allow_lsq: bool
74 /// Whether to permit least squares solving using non-square matrices.
75 ///
76 /// Returns
77 /// -------
78 /// None
79 fn csolve(
80 &mut self,
81 tau: Vec<f64>,
82 y: Vec<$type>,
83 left_n: usize,
84 right_n: usize,
85 allow_lsq: bool
86 ) -> PyResult<()> {
87 self.inner.csolve(&tau, &y, left_n, right_n, allow_lsq)
88 }
89
90 /// Evaluate a single *x* coordinate value on the pp spline.
91 ///
92 /// Parameters
93 /// ----------
94 /// x: float
95 /// The x-axis value at which to evaluate value.
96 ///
97 /// Returns
98 /// -------
99 /// float, Dual or Dual2 based on self
100 ///
101 /// Notes
102 /// -----
103 /// The value of the spline at *x* is the sum of the value of each b-spline
104 /// evaluated at *x* multiplied by the spline coefficients, *c*.
105 ///
106 /// .. math::
107 ///
108 /// \$(x) = \sum_{i=1}^n c_i B_{(i,k,\mathbf{t})}(x)
109 ///
110 fn ppev_single(&self, x: Number) -> PyResult<$type> {
111 match x {
112 Number::F64(f) => self.inner.ppdnev_single(&f, 0),
113 Number::Dual(_) => Err(PyTypeError::new_err(
114 "Cannot index PPSpline with `Dual`, use either `ppev_single(float(x))` or `ppev_single_dual(x)`."
115 )),
116 Number::Dual2(_) => Err(PyTypeError::new_err(
117 "Cannot index PPSpline with `Dual2`, use either `ppev_single(float(x))` or `ppev_single_dual2(x)`.")),
118 }
119 }
120
121 /// Evaluate a single *x* coordinate value on the pp spline.
122 ///
123 /// Parameters
124 /// ----------
125 /// x: Dual
126 /// The x-axis value at which to evaluate value.
127 ///
128 /// Returns
129 /// -------
130 /// Dual
131 ///
132 /// Notes
133 /// -----
134 /// This function guarantees preservation of accurate AD :class:`~rateslib.dual.Dual`
135 /// sensitivities. It also prohibits type mixing and will raise if *Dual2* data types
136 /// are encountered.
137 fn ppev_single_dual(&self, x: Number) -> PyResult<Dual> {
138 match x {
139 Number::F64(f) => self.inner.ppdnev_single_dual(&Dual::new(f, vec![]), 0),
140 Number::Dual(d) => self.inner.ppdnev_single_dual(&d, 0),
141 Number::Dual2(_) => Err(PyTypeError::new_err("Cannot mix `Dual2` and `Dual` types, use `ppev_single_dual2(x)`.")),
142 }
143 }
144
145 /// Evaluate a single *x* coordinate value on the pp spline.
146 ///
147 /// Parameters
148 /// ----------
149 /// x: Dual2
150 /// The x-axis value at which to evaluate value.
151 ///
152 /// Returns
153 /// -------
154 /// Dual2
155 ///
156 /// Notes
157 /// -----
158 /// This function guarantees preservation of accurate AD :class:`~rateslib.dual.Dual2`
159 /// sensitivities. It also prohibits type mixing and will raise if *Dual* data types
160 /// are encountered.
161 fn ppev_single_dual2(&self, x: Number) -> PyResult<Dual2> {
162 match x {
163 Number::F64(f) => self.inner.ppdnev_single_dual2(&Dual2::new(f, vec![]), 0),
164 Number::Dual(_) => Err(PyTypeError::new_err("Cannot mix `Dual2` and `Dual` types, use `ppev_single_dual(x)`.")),
165 Number::Dual2(d) => self.inner.ppdnev_single_dual2(&d, 0),
166 }
167 }
168
169 /// Evaluate an array of *x* coordinates derivatives on the pp spline.
170 ///
171 /// Repeatedly applies :meth:`~rateslib.splines.PPSplineF64.ppev_single`, and
172 /// is typically used for minor performance gains in chart plotting.
173 ///
174 /// .. warning::
175 ///
176 /// The *x* coordinates supplied to this function are treated as *float*, or are
177 /// **converted** to *float*. Therefore it does not guarantee the preservation of AD
178 /// sensitivities. If you need to index by *x* values which are
179 /// :class:`~rateslib.dual.Dual` or :class:`~rateslib.dual.Dual2`, then
180 /// you should choose to iteratively map the
181 /// provided methods :meth:`~rateslib.splines.PPSplineF64.ppev_single_dual` or
182 /// :meth:`~rateslib.splines.PPSplineF64.ppev_single_dual2` respectively.
183 ///
184 /// Returns
185 /// -------
186 /// 1-d array of float
187 fn ppev(&self, x: Vec<f64>) -> PyResult<Vec<$type>> {
188 let out: Vec<$type> = x.iter().map(|v| self.inner.ppdnev_single(&v, 0)).collect::<Result<Vec<$type>, _>>()?;
189 Ok(out)
190 }
191
192 /// Evaluate a single *x* coordinate derivative from the right on the pp spline.
193 ///
194 /// Parameters
195 /// ----------
196 /// x: float
197 /// The x-axis value at which to evaluate value.
198 /// m: int
199 /// The order of derivative to calculate value for (0 is function value).
200 ///
201 /// Returns
202 /// -------
203 /// float, Dual, or Dual2, based on self
204 ///
205 /// Notes
206 /// -----
207 /// The value of derivatives of the spline at *x* is the sum of the value of each
208 /// b-spline derivatives evaluated at *x* multiplied by the spline
209 /// coefficients, *c*.
210 ///
211 /// Due to the definition of the splines this derivative will return the value
212 /// from the right at points where derivatives are discontinuous.
213 ///
214 /// .. math::
215 ///
216 /// \frac{d^m\$(x)}{d x^m} = \sum_{i=1}^n c_i \frac{d^m B_{(i,k,\mathbf{t})}(x)}{d x^m}
217 fn ppdnev_single(&self, x: Number, m: usize) -> PyResult<$type> {
218 match x {
219 Number::Dual(_) => Err(PyTypeError::new_err("Splines cannot be indexed with Duals use `float(x)`.")),
220 Number::F64(f) => self.inner.ppdnev_single(&f, m),
221 Number::Dual2(_) => Err(PyTypeError::new_err("Splines cannot be indexed with Duals use `float(x)`.")),
222 }
223 }
224
225 /// Evaluate a single *x* coordinate derivative from the right on the pp spline.
226 ///
227 /// Parameters
228 /// ----------
229 /// x: Dual
230 /// The x-axis value at which to evaluate value.
231 /// m: int
232 /// The order of derivative to calculate value for (0 is function value).
233 ///
234 /// Returns
235 /// -------
236 /// Dual
237 ///
238 /// Notes
239 /// -----
240 /// This function guarantees preservation of accurate AD :class:`~rateslib.dual.Dual`
241 /// sensitivities. It also prohibits type mixing and will raise if any *Dual2*
242 /// data types are encountered.
243 fn ppdnev_single_dual(&self, x: Number, m: usize) -> PyResult<Dual> {
244 match x {
245 Number::F64(f) => self.inner.ppdnev_single_dual(&Dual::new(f, vec![]), m),
246 Number::Dual(d) => self.inner.ppdnev_single_dual(&d, m),
247 Number::Dual2(_) => Err(PyTypeError::new_err("Cannot mix `Dual2` and `Dual` types, use `ppdnev_single_dual2(x)`.")),
248 }
249 }
250
251 /// Evaluate a single *x* coordinate derivative from the right on the pp spline.
252 ///
253 /// Parameters
254 /// ----------
255 /// x: Dual2
256 /// The x-axis value at which to evaluate value.
257 /// m: int
258 /// The order of derivative to calculate value for (0 is function value).
259 ///
260 /// Returns
261 /// -------
262 /// Dual2
263 ///
264 /// Notes
265 /// -----
266 /// This function guarantees preservation of accurate AD :class:`~rateslib.dual.Dual2`
267 /// sensitivities. It also prohibits type mixing and will raise if any *Dual*
268 /// data types are encountered.
269 fn ppdnev_single_dual2(&self, x: Number, m: usize) -> PyResult<Dual2> {
270 match x {
271 Number::F64(f) => self.inner.ppdnev_single_dual2(&Dual2::new(f, vec![]), m),
272 Number::Dual(_) => Err(PyTypeError::new_err("Cannot mix `Dual2` and `Dual` types, use `ppdnev_single_dual(x)`.")),
273 Number::Dual2(d) => self.inner.ppdnev_single_dual2(&d, m),
274 }
275 }
276
277 /// Evaluate an array of x coordinates derivatives on the pp spline.
278 ///
279 /// Repeatedly applies :meth:`~rateslib.splines.PPSplineF64.ppdnev_single`.
280 ///
281 /// .. warning::
282 ///
283 /// The *x* coordinates supplied to this function are treated as *float*, or are
284 /// **converted** to *float*. Therefore it does not guarantee the preservation of AD
285 /// sensitivities.
286 ///
287 /// Parameters
288 /// ----------
289 /// x: 1-d array of float
290 /// x-axis coordinates.
291 /// m: int
292 /// The order of derivative to calculate value for.
293 ///
294 /// Returns
295 /// -------
296 /// 1-d array of float
297 fn ppdnev(&self, x: Vec<f64>, m: usize) -> PyResult<Vec<$type>> {
298 let out: Vec<$type> = x.iter().map(|v| self.inner.ppdnev_single(&v, m)).collect::<Result<Vec<$type>, _>>()?;
299 Ok(out)
300 }
301
302 /// Evaluate value of the *i* th b-spline at x coordinates.
303 ///
304 /// Repeatedly applies :meth:`~rateslib.splines.bsplev_single`.
305 ///
306 /// .. warning::
307 ///
308 /// The *x* coordinates supplied to this function are treated as *float*, or are
309 /// **converted** to *float*. Therefore it does not guarantee the preservation of AD
310 /// sensitivities.
311 ///
312 /// Parameters
313 /// ----------
314 /// x: 1-d array of float
315 /// x-axis coordinates
316 /// i: int
317 /// Index of the B-spline to evaluate.
318 ///
319 /// Returns
320 /// -------
321 /// 1-d array of float
322 fn bsplev(&self, x: Vec<f64>, i: usize) -> PyResult<Vec<f64>> {
323 Ok(self.inner.bspldnev(&x, &i, &0))
324 }
325
326 /// Evaluate *m* order derivative on the *i* th b-spline at *x* coordinates.
327 ///
328 /// Repeatedly applies :meth:`~rateslib.splines.bspldnev_single`.
329 ///
330 /// .. warning::
331 ///
332 /// The *x* coordinates supplied to this function are treated as *float*, or are
333 /// **converted** to *float*. Therefore it does not guarantee the preservation of AD
334 /// sensitivities.
335 ///
336 /// Parameters
337 /// ----------
338 /// x: 1-d array of float
339 /// x-axis coordinates.
340 /// i: int
341 /// The index of the B-spline to evaluate.
342 /// m: int
343 /// The order of derivative to calculate value for.
344 ///
345 /// Returns
346 /// -------
347 /// 1-d array
348 fn bspldnev(&self, x: Vec<f64>, i: usize, m: usize) -> PyResult<Vec<f64>> {
349 Ok(self.inner.bspldnev(&x, &i, &m))
350 }
351
352 /// Evaluate the 2d spline collocation matrix at each data site.
353 ///
354 /// Parameters
355 /// ----------
356 /// tau: 1-d array of float
357 /// The data sites `x`-axis values which will instruct the pp spline.
358 /// left_n: int
359 /// The order of derivative to use for the left most data site and top row
360 /// of the spline collocation matrix.
361 /// right_n: int
362 /// The order of derivative to use for the right most data site and bottom row
363 /// of the spline collocation matrix.
364 ///
365 /// Returns
366 /// -------
367 /// 2-d array of float
368 ///
369 /// Notes
370 /// -----
371 /// The spline collocation matrix is defined as,
372 ///
373 /// .. math::
374 ///
375 /// [\mathbf{B}_{k, \mathbf{t}}(\mathbf{\tau})]_{j,i} = B_{i,k,\mathbf{t}}(\tau_j)
376 ///
377 /// where each row is a call to :meth:`~rateslib.splines.PPSplineF64.bsplev`, except the top and bottom rows
378 /// which can be specifically adjusted to account for
379 /// ``left_n`` and ``right_n`` such that, for example, the first row might be,
380 ///
381 /// .. math::
382 ///
383 /// [\mathbf{B}_{k, \mathbf{t}}(\mathbf{\tau})]_{1,i} = \frac{d^n}{dx}B_{i,k,\mathbf{t}}(\tau_1)
384 fn bsplmatrix<'py>(
385 &'py self,
386 py: Python<'py>,
387 tau: Vec<f64>,
388 left_n: usize,
389 right_n: usize
390 ) -> PyResult<Bound<'py, PyArray2<f64>>> {
391 Ok(self.inner.bsplmatrix(&tau, left_n, right_n).to_pyarray(py))
392 }
393
394 fn __eq__(&self, other: &Self) -> PyResult<bool> {
395 Ok(self.inner.eq(&other.inner))
396 }
397
398 fn __copy__(&self) -> Self {
399 $name { inner: self.inner.clone() }
400 }
401
402 fn __repr__(&self) -> String {
403 format!("<rl.{} at {:p}>", stringify!($name) ,self)
404 }
405
406 // JSON
407 /// Return a JSON representation of the object.
408 ///
409 /// Returns
410 /// -------
411 /// str
412 #[pyo3(name = "to_json")]
413 fn to_json_py(&self) -> PyResult<String> {
414 match DeserializedObj::$name(self.clone()).to_json() {
415 Ok(v) => Ok(v),
416 Err(_) => Err(PyValueError::new_err("Failed to serialize `PPSpline` to JSON.")),
417 }
418 }
419 }
420 };
421}
422
423create_interface!(PPSplineF64, f64);
424create_interface!(PPSplineDual, Dual);
425create_interface!(PPSplineDual2, Dual2);
426
427/// Calculate the value of an indexed b-spline at *x*.
428///
429/// Parameters
430/// ----------
431/// x: float
432/// The *x* value at which to evaluate the b-spline.
433/// i: int
434/// The index of the b-spline to evaluate.
435/// k: int
436/// The order of the b-spline (note that k=4 is a cubic spline).
437/// t: sequence of float
438/// The knot sequence of the pp spline.
439/// org_k: int, optional
440/// The original k input. Used only internally when recursively calculating
441/// successive b-splines. Users will not typically use this parameters.
442///
443/// Notes
444/// -----
445/// B-splines can be recursively defined as:
446///
447/// .. math::
448///
449/// B_{i,k,\mathbf{t}}(x) = \frac{x-t_i}{t_{i+k-1}-t_i}B_{i,k-1,\mathbf{t}}(x) + \frac{t_{i+k}-x}{t_{i+k}-t_{i+1}}B_{i+1,k-1,\mathbf{t}}(x)
450///
451/// and such that the basic, stepwise, b-spline or order 1 are:
452///
453/// .. math::
454///
455/// B_{i,1,\mathbf{t}}(x) = \left \{ \begin{matrix} 1, & t_i \leq x < t_{i+1} \\ 0, & \text{otherwise} \end{matrix} \right .
456///
457/// For continuity on the right boundary the rightmost basic b-spline is also set equal
458/// to 1 there: :math:`B_{n,1,\mathbf{t}}(t_{n+k})=1`.
459#[pyfunction]
460#[pyo3(signature = (x, i, k, t, org_k=None))]
461pub(crate) fn bsplev_single(
462 x: f64,
463 i: usize,
464 k: usize,
465 t: Vec<f64>,
466 org_k: Option<usize>,
467) -> PyResult<f64> {
468 Ok(bsplev_single_f64(&x, i, &k, &t, org_k))
469}
470
471/// Calculate the *m* th order derivative (from the right) of an indexed b-spline at *x*.
472///
473/// Parameters
474/// ----------
475/// x: float
476/// The *x* value at which to evaluate the b-spline.
477/// i: int
478/// The index of the b-spline to evaluate.
479/// k: int
480/// The order of the b-spline (note that k=4 is a cubic spline).
481/// t: sequence of float
482/// The knot sequence of the pp spline.
483/// m: int
484/// The order of the derivative of the b-spline to evaluate.
485/// org_k: int, optional
486/// The original k input. Used only internally when recursively calculating
487/// successive b-splines. Users will not typically use this parameter.
488///
489/// Notes
490/// -----
491/// B-splines derivatives can be recursively defined as:
492///
493/// .. math::
494///
495/// \frac{d}{dx}B_{i,k,\mathbf{t}}(x) = (k-1) \left ( \frac{B_{i,k-1,\mathbf{t}}(x)}{t_{i+k-1}-t_i} - \frac{B_{i+1,k-1,\mathbf{t}}(x)}{t_{i+k}-t_{i+1}} \right )
496///
497/// and such that the basic, stepwise, b-spline derivative is:
498///
499/// .. math::
500///
501/// \frac{d}{dx}B_{i,1,\mathbf{t}}(x) = 0
502///
503/// During this recursion the original order of the spline is registered so that under
504/// the given knot sequence, :math:`\mathbf{t}`, lower order b-splines which are not
505/// the rightmost will register a unit value. For example, the 4'th order knot sequence
506/// [1,1,1,1,2,2,2,3,4,4,4,4] defines 8 b-splines. The rightmost is measured
507/// across the knots [3,4,4,4,4]. When the knot sequence remains constant and the
508/// order is lowered to 3 the rightmost, 9'th, b-spline is measured across [4,4,4,4],
509/// which is effectively redundant since its domain has zero width. The 8'th b-spline
510/// which is measured across the knots [3,4,4,4] is that which will impact calculations
511/// and is therefore given the value 1 at the right boundary. This is controlled by
512/// the information provided by ``org_k``.
513///
514/// Examples
515/// --------
516/// The derivative of the 4th b-spline of the following knot sequence
517/// is discontinuous at `x` = 2.0.
518///
519/// .. ipython:: python
520/// :suppress:
521///
522/// from rateslib import bspldnev_single
523///
524/// .. ipython:: python
525///
526/// t = [1,1,1,1,2,2,2,3,4,4,4,4]
527/// bspldnev_single(x=2.0, i=3, k=4, t=t, m=1)
528/// bspldnev_single(x=1.99999999, i=3, k=4, t=t, m=1)
529///
530/// .. plot::
531///
532/// from rateslib.splines import *
533/// import matplotlib.pyplot as plt
534/// from datetime import datetime as dt
535/// import numpy as np
536/// t = [1,1,1,1,2,2,2,3,4,4,4,4]
537/// spline = PPSplineF64(k=4, t=t)
538/// x = np.linspace(1, 4, 76)
539/// fix, ax = plt.subplots(1,1)
540/// ax.plot(x, spline.bspldnev(x, 3, 0))
541/// plt.show()
542#[pyfunction]
543#[pyo3(signature = (x, i, k, t, m, org_k=None))]
544pub(crate) fn bspldnev_single(
545 x: f64,
546 i: usize,
547 k: usize,
548 t: Vec<f64>,
549 m: usize,
550 org_k: Option<usize>,
551) -> PyResult<f64> {
552 Ok(bspldnev_single_f64(&x, i, &k, &t, m, org_k))
553}