LSSTApplications  10.0+286,10.0+36,10.0+46,10.0-2-g4f67435,10.1+152,10.1+37,11.0,11.0+1,11.0-1-g47edd16,11.0-1-g60db491,11.0-1-g7418c06,11.0-2-g04d2804,11.0-2-g68503cd,11.0-2-g818369d,11.0-2-gb8b8ce7
LSSTDataManagementBasePackage
Public Member Functions | List of all members
lsst::afw::math::detail::SmoothedSpline Class Reference

#include <Spline.h>

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

Public Member 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)
 
- Public Member Functions inherited from lsst::afw::math::detail::Spline
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
 

Additional Inherited Members

- Protected Member Functions inherited from lsst::afw::math::detail::Spline
 Spline ()
 
void _allocateSpline (int const nknot)
 
- Protected Attributes inherited from lsst::afw::math::detail::Spline
std::vector< double > _knots
 
std::vector< std::vector
< double > > 		_coeffs
 

Detailed Description

Definition at line 57 of file Spline.h.

Constructor & Destructor Documentation

lsst::afw::math::detail::SmoothedSpline::SmoothedSpline ( std::vector< double > const &  x,
std::vector< double > const &  f,
std::vector< double > const &  df,
double  s,
double *  chisq = NULL,
std::vector< double > *  errs = NULL 
)

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

purpose - cubic spline data smoother

arguments:

Parameters
xarray 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) .lt. x(i+1)
fvector of length >= 3 containing the ordinates (or function values) of the data points
dfvector of standard deviations of the error associated with the data point; each df[] must be positive.

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 df (see above).

Notes:

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)

Parameters
xpoints where function's specified; monotonic increasing
fvalues of function at x
dferror in function at x
sdesired chisq
chisqfinal chisq (if non-NULL)
errserror estimates, (if non-NULL). You'll need to delete it

Definition at line 710 of file Spline.cc.

718 {
719  float var = 1; // i.e. df is the absolute s.d. N.B. ADD GCV Variant with var=-1
720  int const n = x.size();
721  double const ratio = 2.0;
722  double const tau = 1.618033989; /* golden ratio */
723  double avdf, avar, stat[6];
724  double p, q, delta, r1, r2, r3, r4;
725  double gf1, gf2, gf3, gf4, avh, err;
726 /*
727  * allocate scratch space
728  */
729  _allocateSpline(n);
730 
731  double *y = &_coeffs[0][0];
732  double *c[3];
733  c[0] = &_coeffs[1][0];
734  c[1] = &_coeffs[2][0];
735  c[2] = &_coeffs[3][0];
736 
737  std::vector<double> scratch(7*(n+2)); // scratch space
738 
739  double *r[3];
740  r[0] = &scratch[0] + 1; // we want indices -1..n
741  r[1] = r[0] + (n + 2);
742  r[2] = r[1] + (n + 2);
743  double *t[2];
744  t[0] = r[2] + (n + 2);
745  t[1] = t[0] + (n + 2);
746  double *u = t[1] + (n + 2);
747  double *v = u + (n + 2);
748 /*
749  * and so to work.
750  */
751  std::vector<double> sdf = df; // scaled values of df
752 
753  avdf = spint1(&x[0], &avh, &f[0], &sdf[0], n, y, c, r, t);
754  avar = var;
755  if (var > 0) {
756  avar *= avdf*avdf;
757  }
758 
759  if (var == 0) { /* simply find a natural cubic spline*/
760  r1 = 0;
761 
762  gf1 = spfit1(&x[0], avh, &sdf[0], n, r1, &p, &q, avar, stat, y, c, r, t, u, v);
763  } else { /* Find local minimum of gcv or the
764  expected mean square error */
765  r1 = 1;
766  r2 = ratio*r1;
767  gf2 = spfit1(&x[0], avh, &sdf[0], n, r2, &p, &q, avar, stat, y, c, r, t, u, v);
768  bool set_r3 = false; // was r3 set?
769  for (;;) {
770  gf1 = spfit1(&x[0], avh, &sdf[0], n, r1, &p, &q, avar, stat, y, c, r, t, u,v);
771  if (gf1 > gf2) {
772  break;
773  }
774 
775  if (p <= 0) {
776  break;
777  }
778  r2 = r1;
779  gf2 = gf1;
780  r1 /= ratio;
781  }
782 
783  if(p <= 0) {
784  set_r3 = false;
785  r3 = 0; /* placate compiler */
786  } else {
787  r3 = ratio * r2;
788  set_r3 = true;
789 
790  for (;;) {
791  gf3 = spfit1(&x[0], avh, &sdf[0], n, r3, &p, &q, avar, stat, y, c, r, t, u, v);
792  if (gf3 >= gf2) {
793  break;
794  }
795 
796  if (q <= 0) {
797  break;
798  }
799  r2 = r3;
800  gf2 = gf3;
801  r3 = ratio*r3;
802  }
803  }
804 
805  if(p > 0 && q > 0) {
806  assert (set_r3);
807  r2 = r3;
808  gf2 = gf3;
809  delta = (r2 - r1) / tau;
810  r4 = r1 + delta;
811  r3 = r2 - delta;
812  gf3 = spfit1(&x[0], avh, &sdf[0], n, r3, &p, &q, avar, stat, y, c, r, t, u, v);
813  gf4 = spfit1(&x[0], avh, &sdf[0], n, r4, &p, &q, avar, stat, y, c, r, t, u, v);
814 /*
815  * Golden section search for local minimum
816  */
817  do {
818  if (gf3 <= gf4) {
819  r2 = r4;
820  gf2 = gf4;
821  r4 = r3;
822  gf4 = gf3;
823  delta /= tau;
824  r3 = r2 - delta;
825  gf3 = spfit1(&x[0], avh, &sdf[0], n, r3, &p, &q, avar, stat, y, c, r, t, u, v);
826  } else {
827  r1 = r3;
828  gf1 = gf3;
829  r3 = r4;
830  gf3 = gf4;
831  delta /= tau;
832  r4 = r1 + delta;
833  gf4 = spfit1(&x[0], avh, &sdf[0], n, r4, &p, &q, avar, stat, y, c, r, t, u, v);
834  }
835 
836  err = (r2 - r1) / (r1 + r2);
837  } while(err*err + 1 > 1 && err > 1e-6);
838 
839  r1 = (r1 + r2) * .5;
840  gf1 = spfit1(&x[0], avh, &sdf[0], n, r1, &p, &q, avar, stat, y, c, r, t, u,v);
841  }
842  }
843 /*
844  * Calculate spline coefficients
845  */
846  spcof1(&x[0], avh, &f[0], &sdf[0], n, p, q, y, c, u, v);
847 
848  stat[2] /= avdf*avdf; /* undo scaling */
849  stat[3] /= avdf*avdf;
850  stat[4] /= avdf*avdf;
851 /*
852  * Optionally calculate standard error estimates
853  */
854  if(errs != NULL) {
855  sperr1(&x[0], avh, &sdf[0], n, r, p, avar, errs);
856  }
857 /*
858  * clean up
859  */
860  if(chisq != NULL) {
861  *chisq = n*stat[4];
862  }
863 }
y
int y
Definition: GaussianCentroid.cc:36
lsst::afw::math::detail::Spline::_coeffs
std::vector< std::vector< double > > _coeffs
Definition: Spline.h:33
lsst::afw::math::detail::Spline::_allocateSpline
void _allocateSpline(int const nknot)
Definition: Spline.cc:25
chisq
double chisq
Definition: GaussianCentroid.cc:86
The documentation for this class was generated from the following files:
Generated on Thu Sep 24 2015 02:29:32 for LSSTApplications by   doxygen 1.8.5