File Spline.h

namespace lsst

Class for a simple mapping implementing a generic AstrometryTransform.

Remove all non-astronomical counts from the Chunk Exposure’s pixels.

Forward declarations for lsst::utils::Cache

For details on the Cache class, see the Cache.h file.

It uses a template rather than a pointer so that the derived classes can use the specifics of the transform. The class simplePolyMapping overloads a few routines.

A base class for image defects

Numeric constants used by the Integrate.h integrator routines.

Compute Image Statistics

Note

Gauss-Kronrod-Patterson quadrature coefficients for use in quadpack routine qng. These coefficients were calculated with 101 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov 1981.

Note

The Statistics class itself can only handle lsst::afw::image::MaskedImage() types. The philosophy has been to handle other types by making them look like lsst::afw::image::MaskedImage() and reusing that code. Users should have no need to instantiate a Statistics object directly, but should use the overloaded makeStatistics() factory functions.

namespace afw
namespace math
namespace detail
class SmoothedSpline : public lsst::afw::math::detail::Spline

Public Functions

SmoothedSpline(std::vector<double> const &x, std::vector<double> const &y, std::vector<double> const &dy, double s, double *chisq = NULL, std::vector<double> *errs = NULL)

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

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].

Note

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).

Note

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

Note

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.

Note

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

Note

precision/hardware - double (originally VAX double)

Note

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)

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

class Spline

Subclassed by lsst::afw::math::detail::SmoothedSpline, lsst::afw::math::detail::TautSpline

Public Functions

virtual ~Spline()
Spline(Spline const&)
Spline(Spline&&)
Spline &operator=(Spline const&)
Spline &operator=(Spline&&)
void interpolate(std::vector<double> const &x, std::vector<double> &y) const

Interpolate a Spline.

Parameters

  • [in] x: points to interpolate at

  • [out] y: values of spline interpolation at x

void derivative(std::vector<double> const &x, std::vector<double> &dydx) const

Find the derivative of a Spline.

Parameters

  • [in] x: points to evaluate derivative at

  • [out] dydx: derivatives at x

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). Return a vector of all the roots found

Parameters

  • value: desired value

  • x0x1: specify desired range is [x0,x1)

Protected Functions

Spline()
void _allocateSpline(int const nknot)

Allocate the storage a Spline needs

Protected Attributes

std::vector<double> _knots
std::vector<std::vector<double>> _coeffs
class TautSpline : public lsst::afw::math::detail::Spline

Public Types

enum Symmetry

Values:

Unknown
Odd
Even

Public 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.

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.

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.

Note

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.

Note

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)

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

Private Functions

void calculateTautSpline(std::vector<double> const &x, std::vector<double> const &y, double const gamma0)

Here’s the worker routine for the TautSpline ctor

Parameters

  • x: points where function’s specified

  • y: values of function at tau[]

  • gamma0: control extra knots

void calculateTautSplineEvenOdd(std::vector<double> const &x, std::vector<double> const &y, double const gamma0, bool even)

Fit a taut spline to a set of data, forcing the resulting spline to obey S(x) = +-S(-x). The input points must have tau[] >= 0.

See TautSpline::TautSpline() for a discussion of the algorithm, and the meaning of the parameter gamma

This is done by duplicating the input data for -ve x, so consider carefully before using this function on many-thousand-point datasets