lsst::afw::math::detail::SmoothedSpline Class Reference

#include <Spline.h>

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

Public Member Functions
	SmoothedSpline (std::vector< double > const &x, std::vector< double > const &y, std::vector< double > const &dy, double s, double *chisq=nullptr, std::vector< double > *errs=nullptr)
	Cubic spline data smoother. 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 159 of file Spline.h.

Constructor & Destructor Documentation

SmoothedSpline()

lsst::afw::math::detail::SmoothedSpline::SmoothedSpline	(	std::vector< double > const &	x,
		std::vector< double > const &	y,
		std::vector< double > const &	dy,
		double	s,
		double *	chisq = nullptr,
		std::vector< double > *	errs = nullptr
	)

Cubic spline data smoother.

Algorithm 642 collected algorithms from ACM. Algorithm appeared in Acm-Trans. Math. Software, vol.12, no. 2, Jun., 1986, p. 150.

Translated from fortran by a combination of f2c and RHL.

    Author              - M.F.Hutchinson
                          CSIRO Division of Mathematics and Statistics
                          P.O. Box 1965
                          Canberra, ACT 2601
                          Australia

latest revision - 15 August 1985

Parameters

[in]	x	array of length n containing the abscissae of the n data points (x(i),f(i)) i=0..n-1. x must be ordered so that x(i) < x(i+1)
[in]	y	vector of length >= 3 containing the ordinates (or function values) of the data points
[in]	dy	vector of standard deviations of y the error associated with the data point; each dy[] must be positive.
[in]	s	desired chisq
[out]	chisq	final chisq (if non-NULL)
[out]	errs	error estimates, (if non-NULL). You'll need to delete it

Note

y,c: spline coefficients (output). y is an array of length n; c is an n-1 by 3 matrix. The value of the spline approximation at t is s(t) = c[2][i]*d^3 + c[1][i]*d^2 + c[0][i]*d + y[i] where x[i] <= t < x[i+1] and d = t - x[i].

var: error variance. If var is negative (i.e. unknown) then the smoothing parameter is determined by minimizing the generalized cross validation and an estimate of the error variance is returned. If var is non-negative (i.e. known) then the smoothing parameter is determined to minimize an estimate, which depends on var, of the true mean square error. In particular, if var is zero, then an interpolating natural cubic spline is calculated. Set var to 1 if absolute standard deviations have been provided in dy (see above).

Additional information on the fit is available in the stat array. on normal exit the values are assigned as follows: stat[0] = smoothing parameter (= rho/(rho + 1)) stat[1] = estimate of the number of degrees of freedom of the residual sum of squares; this reduces to the usual value of n-2 when a least squares regression line is calculated. stat[2] = generalized cross validation stat[3] = mean square residual stat[4] = estimate of the true mean square error at the data points stat[5] = estimate of the error variance; chi^2/nu in the case of linear regression

If stat[0]==0 (rho==0) an interpolating natural cubic spline has been calculated; if stat[0]==1 (rho==infinite) a least squares regression line has been calculated.

Returns stat[4], an estimate of the true rms error

precision/hardware - double (originally VAX double)

the number of arithmetic operations required by the subroutine is proportional to n. The subroutine uses an algorithm developed by M.F. Hutchinson and F.R. de Hoog, 'Smoothing Noisy Data with Spline Functions', Numer. Math. 47 p.99 (1985)

Definition at line 534 of file Spline.cc.

                                                        {
    float var = 1;  // i.e. df is the absolute s.d.  N.B. ADD GCV Variant with var=-1
    int const n = x.size();
    double const ratio = 2.0;
    double const tau = 1.618033989; /* golden ratio */
    double avdf, avar, stat[6];
    double p, q, delta, r1, r2, r3, r4;
    double gf1, gf2, gf3, gf4, avh, err;
    /*
     * allocate scratch space
     */
    _allocateSpline(n);
 
    double *y = &_coeffs[0][0];
    double *c[3];
    c[0] = &_coeffs[1][0];
    c[1] = &_coeffs[2][0];
    c[2] = &_coeffs[3][0];
 
    std::vector<double> scratch(7 * (n + 2));  // scratch space
 
    double *r[3];
    r[0] = &scratch[0] + 1;  // we want indices -1..n
    r[1] = r[0] + (n + 2);
    r[2] = r[1] + (n + 2);
    double *t[2];
    t[0] = r[2] + (n + 2);
    t[1] = t[0] + (n + 2);
    double *u = t[1] + (n + 2);
    double *v = u + (n + 2);
    /*
     * and so to work.
     */
    std::vector<double> sdf = df;  // scaled values of df
 
    avdf = spint1(&x[0], &avh, &f[0], &sdf[0], n, y, c, r, t);
    avar = var;
    if (var > 0) {
        avar *= avdf * avdf;
    }
 
    if (var == 0) { /* simply find a natural cubic spline*/
        r1 = 0;
 
        gf1 = spfit1(&x[0], avh, &sdf[0], n, r1, &p, &q, avar, stat, y, c, r, t, u, v);
    } else { /* Find local minimum of gcv or the
                expected mean square error */
        r1 = 1;
        r2 = ratio * r1;
        gf2 = spfit1(&x[0], avh, &sdf[0], n, r2, &p, &q, avar, stat, y, c, r, t, u, v);
        bool set_r3 = false;  // was r3 set?
        for (;;) {
            gf1 = spfit1(&x[0], avh, &sdf[0], n, r1, &p, &q, avar, stat, y, c, r, t, u, v);
            if (gf1 > gf2) {
                break;
            }
 
            if (p <= 0) {
                break;
            }
            r2 = r1;
            gf2 = gf1;
            r1 /= ratio;
        }
 
        if (p <= 0) {
            set_r3 = false;
            r3 = 0; /* placate compiler */
        } else {
            r3 = ratio * r2;
            set_r3 = true;
 
            for (;;) {
                gf3 = spfit1(&x[0], avh, &sdf[0], n, r3, &p, &q, avar, stat, y, c, r, t, u, v);
                if (gf3 >= gf2) {
                    break;
                }
 
                if (q <= 0) {
                    break;
                }
                r2 = r3;
                gf2 = gf3;
                r3 = ratio * r3;
            }
        }
 
        if (p > 0 && q > 0) {
            assert(set_r3);
            r2 = r3;
            gf2 = gf3;
            delta = (r2 - r1) / tau;
            r4 = r1 + delta;
            r3 = r2 - delta;
            gf3 = spfit1(&x[0], avh, &sdf[0], n, r3, &p, &q, avar, stat, y, c, r, t, u, v);
            gf4 = spfit1(&x[0], avh, &sdf[0], n, r4, &p, &q, avar, stat, y, c, r, t, u, v);
            /*
             * Golden section search for local minimum
             */
            do {
                if (gf3 <= gf4) {
                    r2 = r4;
                    gf2 = gf4;
                    r4 = r3;
                    gf4 = gf3;
                    delta /= tau;
                    r3 = r2 - delta;
                    gf3 = spfit1(&x[0], avh, &sdf[0], n, r3, &p, &q, avar, stat, y, c, r, t, u, v);
                } else {
                    r1 = r3;
                    gf1 = gf3;
                    r3 = r4;
                    gf3 = gf4;
                    delta /= tau;
                    r4 = r1 + delta;
                    gf4 = spfit1(&x[0], avh, &sdf[0], n, r4, &p, &q, avar, stat, y, c, r, t, u, v);
                }
 
                err = (r2 - r1) / (r1 + r2);
            } while (err * err + 1 > 1 && err > 1e-6);
 
            r1 = (r1 + r2) * .5;
            gf1 = spfit1(&x[0], avh, &sdf[0], n, r1, &p, &q, avar, stat, y, c, r, t, u, v);
        }
    }
    /*
     * Calculate spline coefficients
     */
    spcof1(&x[0], avh, &f[0], &sdf[0], n, p, q, y, c, u, v);
 
    stat[2] /= avdf * avdf; /* undo scaling */
    stat[3] /= avdf * avdf;
    stat[4] /= avdf * avdf;
    /*
     * Optionally calculate standard error estimates
     */
    if (errs != nullptr) {
        sperr1(&x[0], avh, &sdf[0], n, r, p, avar, errs);
    }
    /*
     * clean up
     */
    if (chisq != nullptr) {
        *chisq = n * stat[4];
    }
}

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