B-Spline Interpolation Implementation

[yfnaji]

This pull request is to address issue #5 and is part of the roadmap to v1 (#187).

B-Splines

B-splines are constructed by joining polynomials of degree $d$, resulting in a curve with continuity $C^{d-1}$.

For a B-spline of degree $d$, we define a set of control points:

$$C= [c_0, c_1, \ldots, c_{k-1}]$$

We also define a non-decreasing sequence of knots:

$$T = [t_0, t_1, \ldots, t_{k+d}]$$

where we require a total of $d + k + 1$ knots.

The basis functions are defined as follows:

$$ N_{i, 0}(t) = \begin{cases} 1 & t\in[t_i, t_{i+1}) \\ 0 & \text{otherwise} \end{cases} $$

$$ N_{i, d}(t) = \frac{t-t_i}{t_{i + d} - t_{i}} B_{i, d-1}(t) + \frac{t_{i+d+1} - t}{t_{i+d+1} - t_{i+1}} B_{i+1, d-1}(t) $$

for $i\in\set{0,1,\ldots,k+d-1}$.

The basis functions $B_{i, d}(t)$ for $d>0$ can be evaluated via the Cox de Boor algorithm.

The splines are then constructed as follows (mathematical notation could not be rendered, so an image has been included instead):

The B-spline defined above can only be evaluated within the interval $t\in[t_{d}, t_{k}]$.

Implementation

A B-spline configuration can be set up as follows:

use RustQuant::math::interpolation::BSplineInterpolator;

let knots = vec![0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let control_points = vec![-1.0, 2.0, 0.0, -1.0];

let mut interpolator = BSplineInterpolator::new(knots, control_points.clone(), 2).unwrap();
let _ = interpolator.fit();

An approximation can then be performed:

println!("{}", interpolator.interpolate(2.5).unwrap()); // 1.375

If the length of the knot vector does not align with the number of control points and the degree, an error will be raised:

use RustQuant::math::interpolation::BSplineInterpolator;

let knots = vec![0.0, 1.0, 2.0, 3.0, 4.0,];
let control_points = vec![-1.0, 2.0, 0.0, -1.0];
match BSplineInterpolator::new(knots.clone(), control_points.clone(), 2) {
            Ok(_) => {}
            Err(e) => println!("{}", e.to_string()), // For 4 control points and degree 2, we need 4 + 2 + 1 (7) knots.
}

Additionally, the domain over which you can perform approximations is limited, as mentioned earlier:

let knots = vec![0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let control_points = vec![-1.0, 2.0, 0.0, -1.0];
let mut interpolator = BSplineInterpolator::new(knots, control_points, 2).unwrap();
let _ = interpolator.fit();

match interpolator.interpolate(5.5) {
            Ok(_) => {},
            Err(e) => println!("{}", e.to_string()) // Point 5.5 is outside of the interpolation range [2, 4]
            )
        }
    }

Notes:

• The add_point() method places the "x-axis" coordinate into its appropriate position within the knots vector, maintaining its sorted order. However, the corresponding "y-axis" value is appended to the end of the control_points vector, since the ordering of control points is independent of the knot values and is instead determined by the fitting process.

• Some of the unit tests were benchmarked against SciPy’s BSpline() implementation to verify that the functionality behaves as expected.

[avhz]

awesome work, cheers ! :)