lsst::afw::math::detail::Spline Class Reference

#include <Spline.h>

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

Public Member Functions
virtual	~Spline ()
void	interpolate (std::vector< double > const &x, std::vector< double > &y) const
void	derivative (std::vector< double > const &x, std::vector< double > &dydx) const
std::vector< double >	roots (double const value, double const x0, double const x1) const

Protected Member Functions
	Spline ()
void	_allocateSpline (int const nknot)

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

Detailed Description

Definition at line 12 of file Spline.h.

Constructor & Destructor Documentation

virtual lsst::afw::math::detail::Spline::~Spline ( )

inlinevirtual

Definition at line 14 of file Spline.h.

{}

lsst::afw::math::detail::Spline::Spline ( )

inlineprotected

Definition at line 29 of file Spline.h.

{}

Member Function Documentation

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

protected

Definition at line 25 of file Spline.cc.

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

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

Find the derivative of a Spline.

Parameters

x	points to evaluate derivative at
dydx	derivatives at x

Definition at line 73 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);
    }
}

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

Interpolate a Spline.

Parameters

x	points to interpolate at
y	interpolated values at x

Definition at line 38 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));
    }
}

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

Parameters

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

Definition at line 1472 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::tr1::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

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

protected

Definition at line 33 of file Spline.h.

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

protected

Definition at line 32 of file Spline.h.

---

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

• /home/lsstsw/stack/Linux64/afw/11.0-2-g04d2804/include/lsst/afw/math/detail/Spline.h
• /home/lsstsw/stack/Linux64/afw/11.0-2-g04d2804/src/math/Spline.cc