[/
Copyright (c) 2017 Nick Thompson
Use, modification and distribution are subject to the
Boost Software License, Version 1.0. (See accompanying file
LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
]

[section:cubic_b Cubic B-spline interpolation]

[heading Synopsis]
``
  #include <boost/math/interpolators/cubic_b_spline.hpp>
``

  namespace boost { namespace math {

    template <class Real>
    class cubic_b_spline
    {
    public:

	    template <class BidiIterator>
        cubic_b_spline(BidiIterator a, BidiIterator b, Real left_endpoint, Real step_size,
                       Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
                       Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
        cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
                       Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
                       Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());

        Real operator()(Real x) const;

        Real prime(Real x) const;
    };

  }} // namespaces


[heading Cubic B-Spline Interpolation]

The cubic B-spline class provided by boost allows fast and accurate interpolation of a function which is known at equally spaced points.
The cubic B-spline interpolation is numerically stable as it uses compactly supported basis functions constructed via iterative convolution.
This is to be contrasted to traditional cubic spline interpolation is ill-conditioned as the global support of cubic polynomials causes small changes far from the evaluation point exert a large influence on the calculated value.

There are many use cases for interpolating a function at equally spaced points.
One particularly important example is solving ODE's whose coefficients depend on data determined from experiment or numerical simulation.
Since most ODE steppers are adaptive, they must be able to sample the coefficients at arbitrary points; not just at the points we know the values of our function.

The first two arguments to the constructor are either:

* A pair of bidirectional iterators into the data, or
* A pointer to the data, and a length of the data array.

These are then followed by:

* The start of the functions domain
* The step size

Optionally, you may provide two additional arguments to the constructor, namely the derivative of the function at the left endpoint, and the derivative at the right endpoint.
If you do not provide these arguments, they will be estimated using one-sided finite-difference formulas.
An example of a valid call to the constructor is

    std::vector<double> f{0.01, -0.02, 0.3, 0.8, 1.9, -8.78, -22.6};
    double t0 = 0;
    double h = 0.01;
    boost::math::cubic_b_spline<double> spline(f.begin(), f.end(), t0, h);

The endpoints are estimated using a one-sided finite-difference formula. If you know the derivative at the endpoint, you may pass it to the constructor via

    boost::math::cubic_b_spline<double> spline(f.begin(), f.end(), t0, h, a_prime, b_prime);


To evaluate the interpolant at a point, we simply use

    double y = spline(x);

and to evaluate the derivative of the interpolant we use

    double yp = spline.prime(x);

Be aware that the accuracy guarantees on the derivative of the spline are an order lower than the guarantees on the original function, see [@http://www.springer.com/us/book/9780387984087 Numerical Analysis, Graduate Texts in Mathematics, 181, Rainer Kress] for details.

Finally, note that this is an interpolator, not an extrapolator.
Therefore, you should strenuously avoid evaluating the spline outside the endpoints.
However, it is not an error if you do, as often you cannot control where (say) an ODE stepper will evaluate your function.
As such the interpolant tries to do something reasonable when the passed a value outside the endpoints.
For evaluation within one stepsize of the interval, you can assume something somewhat reasonable was returned.
As you move further away from the endpoints, the interpolant decays to its average on the interval.


[heading Complexity and Performance]

The call to the constructor requires [bigo](/n/) operations, where /n/ is the number of points to interpolate.
Each call the the interpolant is [bigo](1) (constant time).
On the author's Intel Xeon E3-1230, this takes 21ns as long as the vector is small enough to fit in cache.

[heading Accuracy]

Let /h/ be the stepsize. If /f/ is four-times continuously differentiable, then the interpolant is ['[bigo](h[super 4])] accurate and the derivative is ['[bigo](h[super 3])] accurate.

[heading Testing]

Since the interpolant obeys ['s(x[sub j]) = f(x[sub j])] at all interpolation points,
the tests generate random data and evaluate the interpolant at the interpolation points,
validating that equality with the data holds.

In addition, constant, linear, and quadratic functions are interpolated to ensure that the interpolant behaves as expected.

[heading Example]

[import ../../example/cubic_b_spline_example.cpp]

[cubic_b_spline_example]

[cubic_b_spline_example_out]

[endsect]
