lsst::afw::math::detail::TautSpline Class Reference

#include <Spline.h>

Inheritance diagram for lsst::afw::math::detail::TautSpline:

Public Types
enum	Symmetry { Unknown , Odd , Even }

Public Member Functions
	TautSpline (std::vector< double > const &x, std::vector< double > const &y, double const gamma=0, Symmetry type=Unknown)
	Construct cubic spline interpolant to given data. More...
void	interpolate (std::vector< double > const &x, std::vector< double > &y) const
	Interpolate a Spline. More...
void	derivative (std::vector< double > const &x, std::vector< double > &dydx) const
	Find the derivative of a Spline. More...
std::vector< double >	roots (double const value, double const x0, double const x1) const
	Find the roots of Spline - val = 0 in the range [x0, x1). More...

Protected Member Functions
void	_allocateSpline (int const nknot)
	Allocate the storage a Spline needs. More...

Protected Attributes
std::vector< double >	_knots
std::vector< std::vector< double > >	_coeffs

Detailed Description

Definition at line 59 of file Spline.h.

Member Enumeration Documentation

Symmetry

enum lsst::afw::math::detail::TautSpline::Symmetry

Enumerator
Unknown
Odd
Even

Definition at line 61 of file Spline.h.

{ Unknown, Odd, Even };

Constructor & Destructor Documentation

TautSpline()

lsst::afw::math::detail::TautSpline::TautSpline	(	std::vector< double > const &	x,
		std::vector< double > const &	y,
		double const	gamma = 0,
		Symmetry	type = Unknown
	)

Construct cubic spline interpolant to given data.

Adapted from A Practical Guide to Splines by C. de Boor (N.Y. : Springer-Verlag, 1978). (His routine tautsp converted to C by Robert Lupton and then to C++ by an older and grayer Robert Lupton)

If gamma > 0, additional knots are introduced where needed to make the interpolant more flexible locally. This avoids extraneous inflection points typical of cubic spline interpolation at knots to rapidly changing data. Values for gamma are:

• = 0, no additional knots
• in (0, 3), under certain conditions on the given data at points i-1, i, i+1, and i+2, a knot is added in the i-th interval, i=1,...,ntau-3. See notes below. The interpolant gets rounded with increasing gamma. A value of 2.5 for gamma is typical.
• in (3, 6), same as above, except that knots might also be added in intervals in which an inflection point would be permitted. A value of 5.5 for gamma is typical.

Parameters

x	points where function's specified
y	values of function at tau[]
gamma	control extra knots. See main description for details.
type	specify the desired symmetry (e.g. Even)

Exceptions

pex::exceptions::InvalidParameterError

Thrown if x and y do not have the same length or do not have at least two points

Note

on the i-th interval, (tau[i], tau[i+1]), the interpolant is of the form (*) f(u(x)) = a + b*u + c*h(u,z) + d*h(1-u,1-z) , with u = u(x) = (x - tau[i])/dtau[i]. Here, z = z(i) = addg(i+1)/(addg(i) + addg(i+1)) (= .5, in case the denominator vanishes). with ddg(j) = dg(j+1) - dg(j), addg(j) = abs(ddg(j)), dg(j) = divdif(j) = (gtau[j+1] - gtau[j])/dtau[j] and h(u,z) = alpha*u**3 + (1 - alpha)*(max(((u-zeta)/(1-zeta)),0)**3 with alpha(z) = (1-gamma/3)/zeta zeta(z) = 1 - gamma*min((1 - z), 1/3) thus, for 1/3 .le. z .le. 2/3, f is just a cubic polynomial on the interval i. otherwise, it has one additional knot, at tau[i] + zeta*dtau[i] . as z approaches 1, h(.,z) has an increasingly sharp bend near 1, thus allowing f to turn rapidly near the additional knot. in terms of f(j) = gtau[j] and fsecnd[j] = 2*derivative of f at tau[j], the coefficients for (*) are given as a = f(i) - d b = (f(i+1) - f(i)) - (c - d) c = fsecnd[i+1]*dtau[i]**2/hsecnd(1,z) d = fsecnd[i]*dtau[i]**2/hsecnd(1,1-z) hence can be computed once fsecnd[i],i=0,...,ntau-1, is fixed.

f is automatically continuous and has a continuous second derivat- ive (except when z = 0 or 1 for some i). we determine fscnd(.) from the requirement that also the first derivative of f be continuous. in addition, we require that the third derivative be continuous across tau[1] and across tau[ntau-2] . this leads to a strictly diagonally dominant tridiagonal linear system for the fsecnd[i] which we solve by gauss elimination without pivoting.

There must be at least 4 interpolation points for us to fit a taut cubic spline, but if you provide fewer we'll fit a quadratic or linear polynomial (but you must provide at least 2)

Definition at line 86 of file Spline.cc.

                                      {
    if (x.size() != y.size()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                          (boost::format("TautSpline: x and y must have the same size; saw %d %d\n") %
                           x.size() % y.size())
                                  .str());
    }
 
    int const ntau = x.size(); /* size of tau and gtau, must be >= 2*/
    if (ntau < 2) {
        throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                          (boost::format("TautSpline: ntau = %d, should be >= 2\n") % ntau).str());
    }
 
    switch (type) {
        case Unknown:
            calculateTautSpline(x, y, gamma0);
            break;
        case Even:
        case Odd:
            calculateTautSplineEvenOdd(x, y, gamma0, type == Even);
            break;
    }
}

Member Function Documentation

_allocateSpline()

void lsst::afw::math::detail::Spline::_allocateSpline ( int const  nknot )

protectedinherited

Allocate the storage a Spline needs.

Definition at line 21 of file Spline.cc.

                                            {
    _knots.resize(nknot);
    _coeffs.resize(4);
    for (unsigned int i = 0; i != _coeffs.size(); ++i) {
        _coeffs[i].reserve(nknot);
    }
}

derivative()

void lsst::afw::math::detail::Spline::derivative ( std::vector< double > const &  x, std::vector< double > &  dydx  ) const

inherited

Find the derivative of a Spline.

Parameters

[in]	x	points to evaluate derivative at
[out]	dydx	derivatives at x

Definition at line 57 of file Spline.cc.

                                                                                 {
    int const nknot = _knots.size();
    int const n = x.size();
 
    dydx.resize(n);  // may default-construct elements which is a little inefficient
                     /*
                      * For _knots[i] <= x <= _knots[i+1], the * interpolant has the form
                      *    val = _coeff[0][i] +dx*(_coeff[1][i] + dx*(_coeff[2][i]/2 + dx*_coeff[3][i]/6))
                      * with
                      *    dx = x - knots[i]
                      * so the derivative is
                      *    val = _coeff[1][i] + dx*(_coeff[2][i] + dx*_coeff[3][i]/2))
                      */
 
    int ind = -1;  // no idea initially
    for (int i = 0; i != n; ++i) {
        ind = search_array(x[i], &_knots[0], nknot, ind);
 
        if (ind < 0) {  // off bottom
            ind = 0;
        } else if (ind >= nknot) {  // off top
            ind = nknot - 1;
        }
 
        double const dx = x[i] - _knots[ind];
        dydx[i] = _coeffs[1][ind] + dx * (_coeffs[2][ind] + dx * _coeffs[3][ind] / 2);
    }
}

interpolate()

void lsst::afw::math::detail::Spline::interpolate ( std::vector< double > const &  x, std::vector< double > &  y  ) const

inherited

Interpolate a Spline.

Parameters

[in]	x	points to interpolate at
[out]	y	values of spline interpolation at x

Definition at line 29 of file Spline.cc.

                                                                               {
    int const nknot = _knots.size();
    int const n = x.size();
 
    y.resize(n);   // may default-construct elements which is a little inefficient
                   /*
                    * For _knots[i] <= x <= _knots[i+1], the interpolant
                    * has the form
                    *    val = _coeff[0][i] +dx*(_coeff[1][i] + dx*(_coeff[2][i]/2 + dx*_coeff[3][i]/6))
                    * with
                    *    dx = x - knots[i]
                    */
    int ind = -1;  // no idea initially
    for (int i = 0; i != n; ++i) {
        ind = search_array(x[i], &_knots[0], nknot, ind);
 
        if (ind < 0) {  // off bottom
            ind = 0;
        } else if (ind >= nknot) {  // off top
            ind = nknot - 1;
        }
 
        double const dx = x[i] - _knots[ind];
        y[i] = _coeffs[0][ind] +
               dx * (_coeffs[1][ind] + dx * (_coeffs[2][ind] / 2 + dx * _coeffs[3][ind] / 6));
    }
}

roots()

std::vector< double > lsst::afw::math::detail::Spline::roots ( double const  value, double const  x0, double const  x1  ) const

inherited

Find the roots of Spline - val = 0 in the range [x0, x1).

Return a vector of all the roots found

Parameters

value	desired value
x0,x1	specify desired range is [x0,x1)

Definition at line 1226 of file Spline.cc.

                                                                                  {
    /*
     * Strategy: we know that the interpolant has the form
     *    val = coef[0][i] +dx*(coef[1][i] + dx*(coef[2][i]/2 + dx*coef[3][i]/6))
     * so we can use the usual analytic solution for a cubic. Note that the
     * cubic quoted above returns dx, the distance from the previous knot,
     * rather than x itself
     */
    std::vector<double> roots; /* the roots found */
    double x0 = a;             // lower end of current range
    double const x1 = b;
    int const nknot = _knots.size();
 
    int i0 = search_array(x0, &_knots[0], nknot, -1);
    int const i1 = search_array(x1, &_knots[0], nknot, i0);
    assert(i1 >= i0 && i1 <= nknot - 1);
 
    std::vector<double> newRoots;  // the roots we find in some interval
                                   /*
                                    * Deal with special case that x0 may be off one end or the other of
                                    * the array of knots.
                                    */
    if (i0 < 0) {                  /* off bottom */
        i0 = 0;
        do_cubic(_coeffs[3][i0] / 6, _coeffs[2][i0] / 2, _coeffs[1][i0], _coeffs[0][i0] - value, newRoots);
        //
        // Could use
        //    std::transform(newRoots.begin(), newRoots.end(), newRoots.begin(),
        //                   std::bind(std::plus<double>(), _1, _knots[i0]));
        // but let's not
        //
        for (unsigned int j = 0; j != newRoots.size(); ++j) {
            newRoots[j] += _knots[i0];
        }
        keep_valid_roots(roots, newRoots, x0, _knots[i0]);
 
        x0 = _knots[i0];
    } else if (i0 >= nknot) { /* off top */
        i0 = nknot - 1;
        assert(i0 >= 0);
        do_cubic(_coeffs[3][i0] / 6, _coeffs[2][i0] / 2, _coeffs[1][i0], _coeffs[0][i0] - value, newRoots);
 
        for (unsigned int j = 0; j != newRoots.size(); ++j) {
            newRoots[j] += _knots[i0];
        }
        keep_valid_roots(roots, newRoots, x0, x1);
 
        return roots;
    }
    /*
     * OK, now search in main body of spline. Note that i1 may be nknot - 1, and
     * in any case the right hand limit of the last segment is at x1, not a knot
     */
    for (int i = i0; i <= i1; i++) {
        do_cubic(_coeffs[3][i] / 6, _coeffs[2][i] / 2, _coeffs[1][i], _coeffs[0][i] - value, newRoots);
 
        for (unsigned int j = 0; j != newRoots.size(); ++j) {
            newRoots[j] += _knots[i];
        }
        keep_valid_roots(roots, newRoots, ((i == i0) ? x0 : _knots[i]), ((i == i1) ? x1 : _knots[i + 1]));
    }
 
    return roots;
}

Member Data Documentation

_coeffs

std::vector<std::vector<double> > lsst::afw::math::detail::Spline::_coeffs

protectedinherited

Definition at line 56 of file Spline.h.

_knots

std::vector<double> lsst::afw::math::detail::Spline::_knots

protectedinherited

Definition at line 55 of file Spline.h.

---

The documentation for this class was generated from the following files:

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-0.7.0/Linux64/afw/22.0.1-42-gca6935d93+ba5e5ca3eb/include/lsst/afw/math/detail/Spline.h
• /j/snowflake/release/lsstsw/stack/lsst-scipipe-0.7.0/Linux64/afw/22.0.1-42-gca6935d93+ba5e5ca3eb/src/math/Spline.cc