rgsl/types/
basis_spline.rs

1//
2// A rust binding for the GSL library by Guillaume Gomez (guillaume1.gomez@gmail.com)
3//
4
5/*!
6B-splines are commonly used as basis functions to fit smoothing curves to large data sets.
7To do this, the abscissa axis is broken up into some number of intervals, where the endpoints of
8each interval are called breakpoints.
9
10These breakpoints are then converted to knots by imposing various continuity and smoothness
11conditions at each interface. Given a nondecreasing knot vector t = {t_0, t_1, …, t_{n+k-1}}, the n
12basis splines of order k are defined by
13
14```latex
15B_(i,1)(x) = (1, t_i <= x < t_(i+1)
16
17             (0, else
18
19B_(i,k)(x) = [(x - t_i)/(t_(i+k-1) - t_i)] B_(i,k-1)(x)
20
21              + [(t_(i+k) - x)/(t_(i+k) - t_(i+1))] B_(i+1,k-1)(x)
22```
23
24for i = 0, …, n-1. The common case of cubic B-splines is given by k = 4. The above recurrence
25relation can be evaluated in a numerically stable way by the de Boor algorithm.
26
27If we define appropriate knots on an interval [a,b] then the B-spline basis functions form a
28complete set on that interval. Therefore we can expand a smoothing function as
29
30f(x) = \sum_i c_i B_(i,k)(x)
31
32given enough (x_j, f(x_j)) data pairs. The coefficients c_i can be readily obtained from a
33least-squares fit.
34
35### References and Further Reading
36
37Further information on the algorithms described in this section can be found in the following book,
38
39C. de Boor, A Practical Guide to Splines (1978), Springer-Verlag, ISBN 0-387-90356-9.
40Further information of Greville abscissae and B-spline collocation can be found in the following
41paper,
42
43Richard W. Johnson, Higher order B-spline collocation at the Greville abscissae. Applied Numerical
44Mathematics. vol. 52, 2005, 63–75.
45
46A large collection of B-spline routines is available in the PPPACK library available at
47http://www.netlib.org/pppack, which is also part of SLATEC.
48!*/
49
50use crate::Value;
51use ffi::FFI;
52use types::VectorF64;
53
54ffi_wrapper!(
55    BSpLineWorkspace,
56    *mut sys::gsl_bspline_workspace,
57    gsl_bspline_free
58);
59
60impl BSpLineWorkspace {
61    /// This function allocates a workspace for computing B-splines of order k.
62    ///
63    /// The number of breakpoints is given by nbreak. This leads to n = nbreak + k - 2 basis
64    /// functions.
65    ///
66    /// Cubic B-splines are specified by k = 4. The size of the workspace is O(5k + nbreak).
67    #[doc(alias = "gsl_bspline_alloc")]
68    pub fn new(k: usize, nbreak: usize) -> Option<Self> {
69        let tmp = unsafe { sys::gsl_bspline_alloc(k, nbreak) };
70
71        if tmp.is_null() {
72            None
73        } else {
74            Some(Self::wrap(tmp))
75        }
76    }
77
78    /// This function computes the knots associated with the given breakpoints and stores them
79    /// internally in w->knots.
80    #[doc(alias = "gsl_bspline_knots")]
81    pub fn knots(&mut self, breakpts: &VectorF64) -> Result<(), Value> {
82        let ret = unsafe { sys::gsl_bspline_knots(breakpts.unwrap_shared(), self.unwrap_unique()) };
83        result_handler!(ret, ())
84    }
85
86    /// This function assumes uniformly spaced breakpoints on [a,b] and constructs the corresponding
87    /// knot vector using the previously specified nbreak parameter.
88    /// The knots are stored in w->knots.
89    #[doc(alias = "gsl_bspline_knots_uniform")]
90    pub fn knots_uniform(&mut self, a: f64, b: f64) -> Result<(), Value> {
91        let ret = unsafe { sys::gsl_bspline_knots_uniform(a, b, self.unwrap_unique()) };
92        result_handler!(ret, ())
93    }
94
95    /// This function evaluates all B-spline basis functions at the position x and stores them in
96    /// the vector B, so that the i-th element is B_i(x).
97    ///
98    /// The vector B must be of length n = nbreak + k - 2. This value may also be obtained by
99    /// calling gsl_bspline_ncoeffs.
100    ///
101    /// Computing all the basis functions at once is more efficient than computing them
102    /// individually, due to the nature of the defining recurrence relation.
103    #[doc(alias = "gsl_bspline_eval")]
104    pub fn eval(&mut self, x: f64, B: &mut VectorF64) -> Result<(), Value> {
105        let ret = unsafe { sys::gsl_bspline_eval(x, B.unwrap_unique(), self.unwrap_unique()) };
106        result_handler!(ret, ())
107    }
108
109    /// This function evaluates all potentially nonzero B-spline basis functions at the position x
110    /// and stores them in the vector Bk, so that the i-th element is B_(istart+i)(x).
111    ///
112    /// The last element of Bk is B_(iend)(x). The vector Bk must be of length k.
113    /// By returning only the nonzero basis functions, this function allows quantities involving
114    /// linear combinations of the B_i(x) to be computed without unnecessary terms (such linear
115    /// combinations occur, for example, when evaluating an interpolated function).
116    ///
117    /// Returns `(Value, istart, iend)`.
118    // checker:ignore
119    #[doc(alias = "gsl_bspline_eval_nonzero")]
120    pub fn eval_non_zero(&mut self, x: f64, Bk: &mut VectorF64) -> Result<(usize, usize), Value> {
121        let mut istart = 0;
122        let mut iend = 0;
123        let ret = unsafe {
124            sys::gsl_bspline_eval_nonzero(
125                x,
126                Bk.unwrap_unique(),
127                &mut istart,
128                &mut iend,
129                self.unwrap_unique(),
130            )
131        };
132        result_handler!(ret, (istart, iend))
133    }
134
135    /// This function returns the number of B-spline coefficients given by n = nbreak + k - 2.
136    #[doc(alias = "gsl_bspline_ncoeffs")]
137    pub fn ncoeffs(&mut self) -> usize {
138        unsafe { sys::gsl_bspline_ncoeffs(self.unwrap_unique()) }
139    }
140
141    /// The Greville abscissae are defined to be the mean location of k-1 consecutive knots in the
142    /// knot vector for each basis spline function of order k.
143    ///
144    /// With the first and last knots in the gsl_bspline_workspace knot vector excluded, there are
145    /// gsl_bspline_ncoeffs Greville abscissae for any given B-spline basis.
146    ///
147    /// These values are often used in B-spline collocation applications and may also be called
148    /// Marsden-Schoenberg points.
149    ///
150    /// Returns the location of the i-th Greville abscissa for the given B-spline basis.
151    /// For the ill-defined case when k=1, the implementation chooses to return breakpoint interval
152    /// midpoints.
153    #[doc(alias = "gsl_bspline_greville_abscissa")]
154    pub fn greville_abscissa(&mut self, i: usize) -> f64 {
155        unsafe { sys::gsl_bspline_greville_abscissa(i, self.unwrap_unique()) }
156    }
157}
rgsl/types/basis_spline.rs

rgsl/types/
basis_spline.rs
