LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Spline.cc
Go to the documentation of this file.
1 
7 #include <limits>
8 
9 #include "boost/format.hpp"
12 #include "lsst/afw/geom/Angle.h"
13 
14 namespace afwGeom = lsst::afw::geom;
15 
16 namespace lsst { namespace afw { namespace math { namespace detail {
17 
18  static int search_array(double z, double const *x, int n, int i);
19 
20 /*****************************************************************************/
21 /*
22  * Allocate the storage a Spline needs
23  */
24 void
25 Spline::_allocateSpline(int const nknot)
26 {
27  _knots.resize(nknot);
28  _coeffs.resize(4);
29  for (unsigned int i = 0; i != _coeffs.size(); ++i) {
30  _coeffs[i].reserve(nknot);
31  }
32 }
33 
37 void
38 Spline::interpolate(std::vector<double> const& x,
39  std::vector<double> & y
40  ) const
41 {
42  int const nknot = _knots.size();
43  int const n = x.size();
44 
45  y.resize(n); // may default-construct elements which is a little inefficient
46  /*
47  * For _knots[i] <= x <= _knots[i+1], the interpolant
48  * has the form
49  * val = _coeff[0][i] +dx*(_coeff[1][i] + dx*(_coeff[2][i]/2 + dx*_coeff[3][i]/6))
50  * with
51  * dx = x - knots[i]
52  */
53  int ind = -1; // no idea initially
54  for (int i = 0; i != n; ++i) {
55  ind = search_array(x[i], &_knots[0], nknot, ind);
56 
57  if(ind < 0) { // off bottom
58  ind = 0;
59  } else if(ind >= nknot) { // off top
60  ind = nknot - 1;
61  }
62 
63  double const dx = x[i] - _knots[ind];
64  y[i] = _coeffs[0][ind] + dx*(_coeffs[1][ind] + dx*(_coeffs[2][ind]/2 + dx*_coeffs[3][ind]/6));
65  }
66 }
67 
68 /*****************************************************************************/
72 void
73 Spline::derivative(std::vector<double> const& x,
74  std::vector<double> &dydx
75  ) const
76 {
77  int const nknot = _knots.size();
78  int const n = x.size();
79 
80  dydx.resize(n); // may default-construct elements which is a little inefficient
81  /*
82  * For _knots[i] <= x <= _knots[i+1], the * interpolant has the form
83  * val = _coeff[0][i] +dx*(_coeff[1][i] + dx*(_coeff[2][i]/2 + dx*_coeff[3][i]/6))
84  * with
85  * dx = x - knots[i]
86  * so the derivative is
87  * val = _coeff[1][i] + dx*(_coeff[2][i] + dx*_coeff[3][i]/2))
88  */
89 
90  int ind = -1; // no idea initially
91  for (int i = 0; i != n; ++i) {
92  ind = search_array(x[i], &_knots[0], nknot, ind);
93 
94  if(ind < 0) { // off bottom
95  ind = 0;
96  } else if(ind >= nknot) { // off top
97  ind = nknot - 1;
98  }
99 
100  double const dx = x[i] - _knots[ind];
101  dydx[i] = _coeffs[1][ind] + dx*(_coeffs[2][ind] + dx*_coeffs[3][ind]/2);
102  }
103 }
104 
105 /*****************************************************************************/
178  TautSpline::TautSpline(std::vector<double> const& x,
179  std::vector<double> const& y,
180  double const gamma0,
181  Symmetry type
182  )
183  {
184  if(x.size() != y.size()) {
185  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
186  (boost::format("TautSpline: x and y must have the same size; saw %d %d\n")
187  % x.size() % y.size()).str());
188  }
189 
190  int const ntau = x.size(); /* size of tau and gtau, must be >= 2*/
191  if(ntau < 2) {
192  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
193  (boost::format("TautSpline: ntau = %d, should be >= 2\n")
194  % ntau).str());
195  }
196 
197  switch (type) {
198  case Unknown:
199  calculateTautSpline(x, y, gamma0);
200  break;
201  case Even:
202  case Odd:
203  calculateTautSplineEvenOdd(x, y, gamma0, type == Even);
204  break;
205  }
206  }
207 
208 /************************************************************************************************************/
212  void
213  TautSpline::calculateTautSpline(std::vector<double> const& x,
214  std::vector<double> const& y,
215  double const gamma0
216  )
217  {
218  const double *tau = &x[0];
219  const double *gtau = &y[0];
220  int const ntau = x.size(); // size of tau and gtau, must be >= 2
221 
222  if(ntau < 4) { // use a single quadratic
223  int const nknot = ntau;
224 
225  _allocateSpline(nknot);
226 
227  _knots[0] = tau[0];
228  for (int i = 1; i < nknot;i++) {
229  _knots[i] = tau[i];
230  if(tau[i - 1] >= tau[i]) {
231  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
232  (boost::format("point %d and the next, %f %f, are out of order")
233  % (i - 1) % tau[i - 1] % tau[i]).str());
234  }
235  }
236 
237  if(ntau == 2) {
238  _coeffs[0][0] = gtau[0];
239  _coeffs[1][0] = (gtau[1] - gtau[0])/(tau[1] - tau[0]);
240  _coeffs[2][0] = _coeffs[3][0] = 0;
241 
242  _coeffs[0][1] = gtau[1];
243  _coeffs[1][1] = (gtau[1] - gtau[0])/(tau[1] - tau[0]);
244  _coeffs[2][1] = _coeffs[3][1] = 0;
245  } else { /* must be 3 */
246  double tmp = (tau[2]-tau[0])*(tau[2]-tau[1])*(tau[1]-tau[0]);
247  _coeffs[0][0] = gtau[0];
248  _coeffs[1][0] = ((gtau[1] - gtau[0])*pow(tau[2] - tau[0],2) -
249  (gtau[2] - gtau[0])*pow(tau[1] - tau[0],2))/tmp;
250  _coeffs[2][0] = -2*((gtau[1] - gtau[0])*(tau[2] - tau[0]) -
251  (gtau[2] - gtau[0])*(tau[1] - tau[0]))/tmp;
252  _coeffs[3][0] = 0;
253 
254  _coeffs[0][1] = gtau[1];
255  _coeffs[1][1] = _coeffs[1][0] + (tau[1] - tau[0])*_coeffs[2][0];
256  _coeffs[2][1] = _coeffs[2][0];
257  _coeffs[3][1] = 0;
258 
259  _coeffs[0][2] = gtau[2];
260  _coeffs[1][2] = _coeffs[1][0] + (tau[2] - tau[0])*_coeffs[2][0];
261  _coeffs[2][2] = _coeffs[2][0];
262  _coeffs[3][2] = 0;
263  }
264 
265  return;
266  }
267 /*
268  * Allocate scratch space
269  * s[0][...] = dtau = tau(.+1) - tau
270  * s[1][...] = diag = diagonal in linear system
271  * s[2][...] = u = upper diagonal in linear system
272  * s[3][...] = r = right side for linear system (initially)
273  * = fsecnd = solution of linear system , namely the second
274  * derivatives of interpolant at tau
275  * s[4][...] = z = indicator of additional knots
276  * s[5][...] = 1/hsecnd(1,x) with x = z or = 1-z. see below.
277  */
278  std::vector<std::vector<double> > s(6); // scratch space
279 
280  for (int i = 0; i != 6; i++) {
281  s[i].resize(ntau);
282  }
283 /*
284  * Construct delta tau and first and second (divided) differences of data
285  */
286 
287  for (int i = 0; i < ntau - 1; i++) {
288  s[0][i] = tau[i + 1] - tau[i];
289  if(s[0][i] <= 0.) {
290  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
291  (boost::format("point %d and the next, %f %f, are out of order")
292  % i % tau[i] % tau[i+1]).str());
293  }
294  s[3][i + 1] = (gtau[i + 1] - gtau[i])/s[0][i];
295  }
296  for (int i = 1; i < ntau - 1; ++i) {
297  s[3][i] = s[3][i + 1] - s[3][i];
298  }
299 /*
300  * Construct system of equations for second derivatives at tau. At each
301  * interior data point, there is one continuity equation, at the first
302  * and the last interior data point there is an additional one for a
303  * total of ntau equations in ntau unknowns.
304  */
305  s[1][1] = s[0][0]/3;
306 
307  int method;
308  double gamma = gamma0; // control the smoothing
309  if(gamma <= 0) {
310  method = 1;
311  } else if(gamma > 3) {
312  gamma -= 3;
313  method = 3;
314  } else {
315  method = 2;
316  }
317  double const onemg3 = 1 - gamma/3;
318 
319  int nknot = ntau; // count number of knots
320 /*
321  * Some compilers don't realise that the flow of control always initialises
322  * these variables; so initialise them to placate e.g. gcc
323  */
324  double zeta = 0.0;
325  double alpha = 0.0;
326  double ratio = 0.0;
327 
328  //double c, d;
329  //double z, denom, factr2;
330  //double onemzt, zt2, del;
331 
332 
333  double entry3 = 0.0;
334  double factor = 0.0;
335  double onemzt = 0;
336  double zt2 = 0;
337  double z_half = 0;
338 
339  for (int i = 1; i < ntau - 2; ++i) {
340  /*
341  * construct z[i] and zeta[i]
342  */
343  double z = .5;
344  if((method == 2 && s[3][i]*s[3][i + 1] >= 0) || method == 3) {
345  double const temp = fabs(s[3][i + 1]);
346  double const denom = fabs(s[3][i]) + temp;
347  if(denom != 0) {
348  z = temp/denom;
349  if (fabs(z - 0.5) <= 1.0/6.0) {
350  z = 0.5;
351  }
352  }
353  }
354 
355  s[4][i] = z;
356 /*
357  set up part of the i-th equation which depends on the i-th interval
358 */
359  z_half = z - 0.5;
360  if(z_half < 0) {
361  zeta = gamma*z;
362  onemzt = 1 - zeta;
363  zt2 = zeta*zeta;
364  double temp = onemg3/onemzt;
365  alpha = (temp < 1 ? temp : 1);
366  factor = zeta/(alpha*(zt2 - 1) + 1);
367  s[5][i] = zeta*factor/6;
368  s[1][i] += s[0][i]*((1 - alpha*onemzt)*factor/2 - s[5][i]);
369 /*
370  * if z = 0 and the previous z = 1, then d[i] = 0. Since then
371  * also u[i-1] = l[i+1] = 0, its value does not matter. Reset
372  * d[i] = 1 to insure nonzero pivot in elimination.
373  */
374  if (s[1][i] <= 0.) {
375  s[1][i] = 1;
376  }
377  s[2][i] = s[0][i]/6;
378 
379  if(z != 0) { /* we'll get a new knot */
380  nknot++;
381  }
382  } else if(z_half == 0) {
383  s[1][i] += s[0][i]/3;
384  s[2][i] = s[0][i]/6;
385  } else {
386  onemzt = gamma*(1 - z);
387  zeta = 1 - onemzt;
388  double const temp = onemg3/zeta;
389  alpha = (temp < 1 ? temp : 1);
390  factor = onemzt/(1 - alpha*zeta*(onemzt + 1));
391  s[5][i] = onemzt*factor/6;
392  s[1][i] += s[0][i]/3;
393  s[2][i] = s[5][i]*s[0][i];
394 
395  if(onemzt != 0) { /* we'll get a new knot */
396  nknot++;
397  }
398  }
399  if (i == 1) {
400  s[4][0] = 0.5;
401 /*
402  * the first two equations enforce continuity of the first and of
403  * the third derivative across tau[1].
404  */
405  s[1][0] = s[0][0]/6;
406  s[2][0] = s[1][1];
407  entry3 = s[2][1];
408  if(z_half < 0) {
409  const double factr2 = zeta*(alpha*(zt2 - 1.) + 1.)/(alpha*(zeta*zt2 - 1.) + 1.);
410  ratio = factr2*s[0][1]/s[1][0];
411  s[1][1] = factr2*s[0][1] + s[0][0];
412  s[2][1] = -factr2*s[0][0];
413  } else if (z_half == 0) {
414  ratio = s[0][1]/s[1][0];
415  s[1][1] = s[0][1] + s[0][0];
416  s[2][1] = -s[0][0];
417  } else {
418  ratio = s[0][1]/s[1][0];
419  s[1][1] = s[0][1] + s[0][0];
420  s[2][1] = -s[0][0]*6*alpha*s[5][1];
421  }
422 /*
423  * at this point, the first two equations read
424  * diag[0]*x0 + u[0]*x1 + entry3*x2 = r[1]
425  * -ratio*diag[0]*x0 + diag[1]*x1 + u[1]*x2 = 0
426  * set r[0] = r[1] and eliminate x1 from the second equation
427  */
428  s[3][0] = s[3][1];
429 
430  s[1][1] += ratio*s[2][0];
431  s[2][1] += ratio*entry3;
432  s[3][1] = ratio*s[3][1];
433  } else {
434 /*
435  * the i-th equation enforces continuity of the first derivative
436  * across tau[i]; it reads
437  * -ratio*diag[i-1]*x_{i-1} + diag[i]*x_i + u[i]*x_{i+1} = r[i]
438  * eliminate x_{i-1} from this equation
439  */
440  s[1][i] += ratio*s[2][i - 1];
441  s[3][i] += ratio*s[3][i - 1];
442  }
443 /*
444  * Set up the part of the next equation which depends on the i-th interval.
445  */
446  if(z_half < 0) {
447  ratio = -s[5][i]*s[0][i]/s[1][i];
448  s[1][i + 1] = s[0][i]/3;
449  } else if(z_half == 0) {
450  ratio = -(s[0][i]/6)/s[1][i];
451  s[1][i + 1] = s[0][i]/3;
452  } else {
453  ratio = -(s[0][i]/6)/s[1][i];
454  s[1][i + 1] = s[0][i]*((1 - zeta*alpha)*factor/2 - s[5][i]);
455  }
456  }
457 
458  s[4][ntau - 2] = 0.5;
459 /*
460  * last two equations, which enforce continuity of third derivative and
461  * of first derivative across tau[ntau - 2]
462  */
463  double const entry_ = ratio*s[2][ntau - 3] + s[1][ntau - 2] + s[0][ntau - 2]/3;
464  s[1][ntau - 1] = s[0][ntau - 2]/6;
465  s[3][ntau - 1] = ratio*s[3][ntau - 3] + s[3][ntau - 2];
466  if (z_half < 0) {
467  ratio = s[0][ntau - 2]*6*s[5][ntau - 3]*alpha/s[1][ntau - 3];
468  s[1][ntau - 2] = ratio*s[2][ntau - 3] + s[0][ntau - 2] + s[0][ntau - 3];
469  s[2][ntau - 2] = -s[0][ntau - 3];
470  } else if (z_half == 0) {
471  ratio = s[0][ntau - 2]/s[1][ntau - 3];
472  s[1][ntau - 2] = ratio*s[2][ntau - 3] + s[0][ntau - 2] + s[0][ntau - 3];
473  s[2][ntau - 2] = -s[0][ntau - 3];
474  } else {
475  const double factr2 = onemzt*(alpha*(onemzt*onemzt - 1) + 1)/(alpha*(onemzt*onemzt*onemzt - 1) + 1);
476  ratio = factr2*s[0][ntau - 2]/s[1][ntau - 3];
477  s[1][ntau - 2] = ratio*s[2][ntau - 3] + factr2*s[0][ntau - 3] + s[0][ntau - 2];
478  s[2][ntau - 2] = -factr2*s[0][ntau - 3];
479  }
480 /*
481  * at this point, the last two equations read
482  * diag[i]*x_i + u[i]*x_{i+1} = r[i]
483  * -ratio*diag[i]*x_i + diag[i+1]*x_{i+1} = r[i+1]
484  * eliminate x_i from last equation
485  */
486  s[3][ntau - 2] = ratio*s[3][ntau - 3];
487  ratio = -entry_/s[1][ntau - 2];
488  s[1][ntau - 1] += ratio*s[2][ntau - 2];
489  s[3][ntau - 1] += ratio*s[3][ntau - 2];
490 
491 /*
492  * back substitution
493  */
494  s[3][ntau - 1] /= s[1][ntau - 1];
495  for (int i = ntau - 2; i > 0; --i) {
496  s[3][i] = (s[3][i] - s[2][i]*s[3][i + 1])/s[1][i];
497  }
498 
499  s[3][0] = (s[3][0] - s[2][0]*s[3][1] - entry3*s[3][2])/s[1][0];
500 /*
501  * construct polynomial pieces; first allocate space for the coefficients
502  */
503 #if 1
504 /*
505  * Start by counting the knots
506  */
507  {
508  int const nknot0 = nknot;
509  int nknot = ntau;
510 
511  for (int i = 0; i < ntau - 1; ++i) {
512  double const z = s[4][i];
513  if((z < 0.5 && z != 0.0) || (z > 0.5 && (1 - z) != 0.0)) {
514  nknot++;
515  }
516  }
517  assert (nknot == nknot0);
518  }
519 #endif
520  _allocateSpline(nknot);
521 
522  _knots[0] = tau[0];
523  int j = 0;
524  for (int i = 0; i < ntau - 1; ++i) {
525  _coeffs[0][j] = gtau[i];
526  _coeffs[2][j] = s[3][i];
527  double const divdif = (gtau[i + 1] - gtau[i])/s[0][i];
528  double z = s[4][i];
529  double const z_half = z - 0.5;
530  if (z_half < 0) {
531  if (z == 0) {
532  _coeffs[1][j] = divdif;
533  _coeffs[2][j] = 0;
534  _coeffs[3][j] = 0;
535  } else {
536  zeta = gamma*z;
537  onemzt = 1 - zeta;
538  double const c = s[3][i + 1]/6;
539  double const d = s[3][i]*s[5][i];
540  j++;
541 
542  double const del = zeta*s[0][i];
543  _knots[j] = tau[i] + del;
544  zt2 = zeta*zeta;
545  double temp = onemg3/onemzt;
546  alpha = (temp < 1 ? temp : 1);
547  factor = onemzt*onemzt*alpha;
548  temp = s[0][i];
549  _coeffs[0][j] = gtau[i] + divdif*del + temp*temp*(d*onemzt*(factor - 1) + c*zeta*(zt2 - 1));
550  _coeffs[1][j] = divdif + s[0][i]*(d*(1 - 3*factor) + c*(3*zt2 - 1));
551  _coeffs[2][j] = (d*alpha*onemzt + c*zeta)*6;
552  _coeffs[3][j] = (c - d*alpha)*6/s[0][i];
553  if(del*zt2 == 0) {
554  _coeffs[1][j - 1] = 0; /* would be NaN in an */
555  _coeffs[3][j - 1] = 0; /* 0-length interval */
556  } else {
557  _coeffs[3][j - 1] = _coeffs[3][j] - d*6*(1 - alpha)/(del*zt2);
558  _coeffs[1][j - 1] = _coeffs[1][j] -
559  del*(_coeffs[2][j] - del/2*_coeffs[3][j - 1]);
560  }
561  }
562  } else if (z_half == 0) {
563  _coeffs[1][j] = divdif - s[0][i]*(s[3][i]*2 + s[3][i + 1])/6;
564  _coeffs[3][j] = (s[3][i + 1] - s[3][i])/s[0][i];
565  } else {
566  onemzt = gamma*(1 - z);
567  if (onemzt == 0) {
568  _coeffs[1][j] = divdif;
569  _coeffs[2][j] = 0;
570  _coeffs[3][j] = 0;
571  } else {
572  zeta = 1 - onemzt;
573  double const temp = onemg3/zeta;
574  alpha = (temp < 1 ? temp : 1);
575  double const c = s[3][i + 1]*s[5][i];
576  double const d = s[3][i]/6;
577  double const del = zeta*s[0][i];
578  _knots[j + 1] = tau[i] + del;
579  _coeffs[1][j] = divdif - s[0][i]*(2*d + c);
580  _coeffs[3][j] = (c*alpha - d)*6/s[0][i];
581  j++;
582 
583  _coeffs[3][j] = _coeffs[3][j - 1] +
584  (1 - alpha)*6*c/(s[0][i]*(onemzt*onemzt*onemzt));
585  _coeffs[2][j] = _coeffs[2][j - 1] + del*_coeffs[3][j - 1];
586  _coeffs[1][j] = _coeffs[1][j - 1] +
587  del*(_coeffs[2][j - 1] + del/2*_coeffs[3][j - 1]);
588  _coeffs[0][j] = _coeffs[0][j - 1] +
589  del*(_coeffs[1][j - 1] +
590  del/2*(_coeffs[2][j - 1] + del/3*_coeffs[3][j - 1]));
591  }
592  }
593 
594  j++;
595  _knots[j] = tau[i + 1];
596  }
597 /*
598  * If there are discontinuities some of the knots may be at the same
599  * position; in this case we generated some NaNs above. As they only
600  * occur for 0-length segments, it's safe to replace them by 0s
601  *
602  * Due to the not-a-knot condition, the last set of coefficients isn't
603  * needed (the last-but-one is equivalent), but it makes the book-keeping
604  * easier if we _do_ generate them
605  */
606  double const del = tau[ntau - 1] - _knots[nknot - 2];
607 
608  _coeffs[0][nknot - 1] = _coeffs[0][nknot - 2] +
609  del*(_coeffs[1][nknot - 2] + del*(_coeffs[2][nknot - 2]/2 +
610  del*_coeffs[3][nknot - 2]/6));
611  _coeffs[1][nknot - 1] = _coeffs[1][nknot - 2] +
612  del*(_coeffs[2][nknot - 2] + del*_coeffs[3][nknot - 2]/2);
613  _coeffs[2][nknot - 1] = _coeffs[2][nknot - 2] + del*_coeffs[3][nknot - 2];
614  _coeffs[3][nknot - 1] = _coeffs[3][nknot - 2];
615 
616  assert (j + 1 == nknot);
617 }
618 
619 /*****************************************************************************/
620 /*
621  * Here's the code to fit smoothing splines through data points
622  */
623 static void
624 spcof1(const double x[], double avh, const double y[], const double dy[],
625  int n, double p, double q, double a[], double *c[3], double u[],
626  const double v[]);
627 
628 static void
629 sperr1(const double x[], double avh, const double dy[], int n,
630  double *r[3], double p, double var, std::vector<double> *se);
631 
632 static double
633 spfit1(const double x[], double avh, const double dy[], int n,
634  double rho, double *p, double *q, double var, double stat[],
635  const double a[], double *c[3], double *r[3], double *t[2],
636  double u[], double v[]);
637 
638 static double
639 spint1(const double x[], double *avh, const double y[], double dy[], int n,
640  double a[], double *c[3], double *r[3], double *t[2]);
641 
711  std::vector<double> const& x,
712  std::vector<double> const& f,
713  std::vector<double> const& df,
714  double s,
715  double *chisq,
716  std::vector<double> *errs
717  )
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 }
864 
865 /*****************************************************************************/
866 /*
867  * initializes the arrays c, r and t for one dimensional cubic
868  * smoothing spline fitting by subroutine spfit1. The values
869  * df[i] are scaled so that the sum of their squares is n.
870  * The average of the differences x[i+1] - x[i] is calculated
871  * in avh in order to avoid underflow and overflow problems in spfit1.
872  *
873  * Return the initial rms value of dy.
874  */
875 static double
876 spint1(const double x[],
877  double *avh,
878  const double f[],
879  double df[],
880  int n,
881  double a[],
882  double *c[3],
883  double *r[3],
884  double *t[2])
885 {
886  double avdf;
887  double e, ff, g, h;
888  int i;
889 
890  assert (n >= 3);
891 /*
892  * Get average x spacing in avh
893  */
894  g = 0;
895  for (i = 0; i < n - 1; ++i) {
896  h = x[i + 1] - x[i];
897  assert (h > 0);
898 
899  g += h;
900  }
901  *avh = g/(n - 1);
902 /*
903  * Scale relative weights
904  */
905  g = 0;
906  for (int i = 0; i < n; ++i) {
907  assert(df[i] > 0);
908 
909  g += df[i]*df[i];
910  }
911  avdf = sqrt(g / n);
912 
913  for (i = 0; i < n; ++i) {
914  df[i] /= avdf;
915  }
916 /*
917  * Initialize h,f
918  */
919  h = (x[1] - x[0])/(*avh);
920  ff = (f[1] - f[0])/h;
921 /*
922  * Calculate a,t,r
923  */
924  for (i = 1; i < n - 1; ++i) {
925  g = h;
926  h = (x[i + 1] - x[i])/(*avh);
927  e = ff;
928  ff = (f[i + 1] - f[i])/h;
929  a[i] = ff - e;
930  t[0][i] = (g + h) * 2./3.;
931  t[1][i] = h / 3.;
932  r[2][i] = df[i - 1]/g;
933  r[0][i] = df[i + 1]/h;
934  r[1][i] = -df[i]/g - df[i]/h;
935  }
936 /*
937  * Calculate c = r'*r
938  */
939  r[1][n-1] = 0;
940  r[2][n-1] = 0;
941  r[2][n] = 0;
942 
943  for (i = 1; i < n - 1; i++) {
944  c[0][i] = r[0][i]*r[0][i] + r[1][i]*r[1][i] + r[2][i]*r[2][i];
945  c[1][i] = r[0][i]*r[1][i+1] + r[1][i]*r[2][i+1];
946  c[2][i] = r[0][i]*r[2][i+2];
947  }
948 
949  return(avdf);
950 }
951 
952 /*****************************************************************************/
953 /*
954  * Fits a cubic smoothing spline to data with relative
955  * weighting dy for a given value of the smoothing parameter
956  * rho using an algorithm based on that of C.H. Reinsch (1967),
957  * Numer. Math. 10, 177-183.
958  *
959  * The trace of the influence matrix is calculated using an
960  * algorithm developed by M.F.Hutchinson and F.R.De Hoog (numer.
961  * math., 47 p.99 (1985)), enabling the generalized cross validation
962  * and related statistics to be calculated in order n operations.
963  *
964  * The arrays a, c, r and t are assumed to have been initialized
965  * by the routine spint1. overflow and underflow problems are
966  * avoided by using p=rho/(1 + rho) and q=1/(1 + rho) instead of
967  * rho and by scaling the differences x[i+1] - x[i] by avh.
968  *
969  * The values in dy are assumed to have been scaled so that the
970  * sum of their squared values is n. the value in var, when it is
971  * non-negative, is assumed to have been scaled to compensate for
972  * the scaling of the values in dy.
973  *
974  * the value returned in fun is an estimate of the true mean square
975  * when var is non-negative, and is the generalized cross validation
976  * when var is negative.
977  */
978 static double
979 spfit1(const double x[],
980  double avh,
981  const double dy[],
982  int n,
983  double rho,
984  double *pp,
985  double *pq,
986  double var,
987  double stat[],
988  const double a[],
989  double *c[3],
990  double *r[3],
991  double *t[2],
992  double u[],
993  double v[])
994 {
995  double const eps = std::numeric_limits<double>::epsilon();
996  double e, f, g, h;
997  double fun;
998  int i;
999  double p, q; /* == *pp, *pq */
1000 /*
1001  * Use p and q instead of rho to prevent overflow or underflow
1002  */
1003  if(fabs(rho) < eps) {
1004  p = 0; q = 1;
1005  } else if(fabs(1/rho) < eps) {
1006  p = 1; q = 0;
1007  } else {
1008  p = rho/(1 + rho);
1009  q = 1/(1 + rho);
1010  }
1011 /*
1012  * Rational Cholesky decomposition of p*c + q*t
1013  */
1014  f = 0;
1015  g = 0;
1016  h = 0;
1017  r[0][-1] = r[0][0] = 0;
1018 
1019  for (int i = 1; i < n - 1; ++i) {
1020  r[2][i-2] = g*r[0][i - 2];
1021  r[1][i-1] = f*r[0][i - 1];
1022  {
1023  double tmp = p*c[0][i] + q*t[0][i] - f*r[1][i-1] - g*r[2][i-2];
1024  if(tmp == 0.0) {
1025  r[0][i] = 1e30;
1026  } else {
1027  r[0][i] = 1/tmp;
1028  }
1029  }
1030  f = p*c[1][i] + q*t[1][i] - h*r[1][i-1];
1031  g = h;
1032  h = p*c[2][i];
1033  }
1034 /*
1035  * Solve for u
1036  */
1037  u[-1] = u[0] = 0;
1038  for (int i = 1; i < n - 1; i++) {
1039  u[i] = a[i] - r[1][i-1]*u[i-1] - r[2][i-2]*u[i-2];
1040  }
1041  u[n-1] = u[n] = 0;
1042  for (int i = n - 2; i > 0; i--) {
1043  u[i] = r[0][i]*u[i] - r[1][i]*u[i+1] - r[2][i]*u[i+2];
1044  }
1045 /*
1046  * Calculate residual vector v
1047  */
1048  e = h = 0;
1049  for (int i = 0; i < n - 1; i++) {
1050  g = h;
1051  h = (u[i + 1] - u[i])/((x[i+1] - x[i])/avh);
1052  v[i] = dy[i]*(h - g);
1053  e += v[i]*v[i];
1054  }
1055  v[n-1] = -h*dy[n-1];
1056  e += v[n-1]*v[n-1];
1057 /*
1058  * Calculate upper three bands of inverse matrix
1059  */
1060  r[0][n-1] = r[1][n-1] = r[0][n] = 0;
1061  for (i = n - 2; i > 0; i--) {
1062  g = r[1][i];
1063  h = r[2][i];
1064  r[1][i] = -g*r[0][i+1] - h*r[1][i+1];
1065  r[2][i] = -g*r[1][i+1] - h*r[0][i+2];
1066  r[0][i] = r[0][i] - g*r[1][i] - h*r[2][i];
1067  }
1068 /*
1069  * Calculate trace
1070  */
1071  f = g = h = 0;
1072  for (i = 1; i < n - 1; ++i) {
1073  f += r[0][i]*c[0][i];
1074  g += r[1][i]*c[1][i];
1075  h += r[2][i]*c[2][i];
1076  }
1077  f += 2*(g + h);
1078 /*
1079  * Calculate statistics
1080  */
1081  stat[0] = p;
1082  stat[1] = f*p;
1083  stat[2] = n*e/(f*f + 1e-20);
1084  stat[3] = e*p*p/n;
1085  stat[5] = e*p/((f < 0) ? f - 1e-10 : f + 1e-10);
1086  if (var >= 0) {
1087  fun = stat[3] - 2*var*stat[1]/n + var;
1088 
1089  stat[4] = fun;
1090  } else {
1091  stat[4] = stat[5] - stat[3];
1092  fun = stat[2];
1093  }
1094 
1095  *pp = p; *pq = q;
1096 
1097  return(fun);
1098 }
1099 
1100 /*
1101  * Calculates coefficients of a cubic smoothing spline from
1102  * parameters calculated by spfit1()
1103  */
1104 static void
1105 spcof1(const double x[],
1106  double avh,
1107  const double y[],
1108  const double dy[],
1109  int n,
1110  double p,
1111  double q,
1112  double a[],
1113  double *c[3],
1114  double u[],
1115  const double v[])
1116 {
1117  double h;
1118  int i;
1119  double qh;
1120 
1121  qh = q/(avh*avh);
1122 
1123  for (i = 0; i < n; ++i) {
1124  a[i] = y[i] - p * dy[i] * v[i];
1125  u[i] = qh*u[i];
1126  }
1127 /*
1128  * calculate c
1129  */
1130  for (i = 0; i < n - 1; ++i) {
1131  h = x[i+1] - x[i];
1132  c[2][i] = (u[i + 1] - u[i])/(3*h);
1133  c[0][i] = (a[i + 1] - a[i])/h - (h*c[2][i] + u[i])*h;
1134  c[1][i] = u[i];
1135  }
1136 }
1137 
1138 /*****************************************************************************/
1139 /*
1140  * Calculates Bayesian estimates of the standard errors of the fitted
1141  * values of a cubic smoothing spline by calculating the diagonal elements
1142  * of the influence matrix.
1143  */
1144 static void
1145 sperr1(const double x[],
1146  double avh,
1147  const double dy[],
1148  int n,
1149  double *r[3],
1150  double p,
1151  double var,
1152  std::vector<double> *se)
1153 {
1154  double f, g, h;
1155  int i;
1156  double f1, g1, h1;
1157 /*
1158  * Initialize
1159  */
1160  h = avh/(x[1] - x[0]);
1161  (*se)[0] = 1 - p*dy[0]*dy[0]*h*h*r[0][1];
1162  r[0][0] = r[1][0] = r[2][0] = 0;
1163 /*
1164  * Calculate diagonal elements
1165  */
1166  for (i = 1; i < n - 1; ++i) {
1167  f = h;
1168  h = avh/(x[i+1] - x[i]);
1169  g = -(f + h);
1170  f1 = f*r[0][i-1] + g*r[1][i-1] + h*r[2][i-1];
1171  g1 = f*r[1][i-1] + g*r[0][i] + h*r[1][i];
1172  h1 = f*r[2][i-1] + g*r[1][i] + h*r[0][i+1];
1173  (*se)[i] = 1 - p*dy[i]*dy[i]*(f*f1 + g*g1 + h*h1);
1174  }
1175  (*se)[n-1] = 1 - p*dy[n-1]*dy[n-1]*h*h*r[0][n-2];
1176 /*
1177  * Calculate standard error estimates
1178  */
1179  for (int i = 0; i < n; ++i) {
1180  double const tmp = (*se)[i]*var;
1181  (*se)[i] = (tmp >= 0) ? sqrt(tmp)*dy[i] : 0;
1182  }
1183 }
1184 
1185 /*****************************************************************************/
1186 /*
1187  * <EXTRACT AUTO>
1188  *
1189  * Fit a taut spline to a set of data, forcing the resulting spline to
1190  * obey S(x) = +-S(-x). The input points must have tau[] >= 0.
1191  *
1192  * See TautSpline::TautSpline() for a discussion of the algorithm, and
1193  * the meaning of the parameter gamma
1194  *
1195  * This is done by duplicating the input data for -ve x, so consider
1196  * carefully before using this function on many-thousand-point datasets
1197  */
1198 void
1199 TautSpline::calculateTautSplineEvenOdd(std::vector<double> const& _tau,
1200  std::vector<double> const& _gtau,
1201  double const gamma,
1202  bool const even // ensure Even symmetry
1203  )
1204 {
1205  const double *tau = &_tau[0];
1206  const double *gtau = &_gtau[0];
1207  int const ntau = _tau.size(); // size of tau and gtau, must be >= 2
1208  std::vector<double> x, y; // tau and gtau, extended to -ve tau
1209 
1210  if(tau[0] == 0.0f) {
1211  int const np = 2*ntau - 1;
1212  x.resize(np);
1213  y.resize(np);
1214 
1215  x[ntau - 1] = tau[0]; y[ntau - 1] = gtau[0];
1216  for (int i = 1; i != ntau; ++i) {
1217  if (even) {
1218  x[ntau - 1 + i] = tau[i]; y[ntau - 1 + i] = gtau[i];
1219  x[ntau - 1 - i] = -tau[i]; y[ntau - 1 - i] = gtau[i];
1220  } else {
1221  x[ntau - 1 + i] = tau[i]; y[ntau - 1 + i] = gtau[i];
1222  x[ntau - 1 - i] = -tau[i]; y[ntau - 1 - i] = -gtau[i];
1223  }
1224  }
1225  } else {
1226  int const np = 2*ntau;
1227  x.resize(np);
1228  y.resize(np);
1229 
1230  for (int i = 0; i != ntau; ++i) {
1231  if (even) {
1232  x[ntau + i] = tau[i]; y[ntau + i] = gtau[i];
1233  x[ntau - 1 - i] = -tau[i]; y[ntau - 1 - i] = gtau[i];
1234  } else {
1235  x[ntau + i] = tau[i]; y[ntau + i] = gtau[i];
1236  x[ntau - 1 - i] = -tau[i]; y[ntau - 1 - i] = -gtau[i];
1237  }
1238  }
1239  }
1240 
1241  TautSpline sp(x, y, gamma); // fit a taut spline to x, y
1242 /*
1243  * Now repackage that spline to reflect the original points
1244  */
1245  int ii;
1246  for (ii = sp._knots.size() - 1; ii >= 0; --ii) {
1247  if(sp._knots[ii] < 0.0f) {
1248  break;
1249  }
1250  }
1251  int const i0 = ii + 1;
1252  int const nknot = sp._knots.size() - i0;
1253 
1254  _allocateSpline(nknot);
1255 
1256  for (int i = i0; i != static_cast<int>(sp._knots.size()); ++i) {
1257  _knots[i - i0] = sp._knots[i];
1258  for (int j = 0; j != 4; ++j) {
1259  _coeffs[j][i - i0] = sp._coeffs[j][i];
1260  }
1261  }
1262 }
1263 
1264 /*****************************************************************************/
1265 /*
1266  * returns index i of first element of x >= z; the input i is an initial guess
1267  *
1268  * N.b. we could use std::lower_bound except that we use i as an initial hint
1269  */
1270 static int
1271 search_array(double z, double const *x, int n, int i)
1272 {
1273  int lo, hi, mid;
1274  double xm;
1275 
1276  if(i < 0 || i >= n) { /* initial guess is useless */
1277  lo = -1;
1278  hi = n;
1279  } else {
1280  unsigned int step = 1; /* how much to step up/down */
1281 
1282  if(z > x[i]) { /* expand search upwards */
1283  if(i == n - 1) { /* off top of array */
1284  return(n - 1);
1285  }
1286 
1287  lo = i; hi = lo + 1;
1288  while(z >= x[hi]) {
1289  lo = hi;
1290  step += step; /* double step size */
1291  hi = lo + step;
1292  if(hi >= n) { /* reached top of array */
1293  hi = n - 1;
1294  break;
1295  }
1296  }
1297  } else { /* expand it downwards */
1298  if(i == 0) { /* off bottom of array */
1299  return(-1);
1300  }
1301 
1302  hi = i; lo = i - 1;
1303  while(z < x[lo]) {
1304  hi = lo;
1305  step += step; /* double step size */
1306  lo = hi - step;
1307  if(lo < 0) { /* off bottom of array */
1308  lo = -1;
1309  break;
1310  }
1311  }
1312  }
1313  }
1314 /*
1315  * perform bisection
1316  */
1317  while(hi - lo > 1) {
1318  mid = (lo + hi)/2;
1319  xm = x[mid];
1320  if(z <= xm) {
1321  hi = mid;
1322  } else {
1323  lo = mid;
1324  }
1325  }
1326 
1327  if(lo == -1) { /* off the bottom */
1328  return(lo);
1329  }
1330 /*
1331  * If there's a discontinuity (many knots at same x-value), choose the
1332  * largest
1333  */
1334  xm = x[lo];
1335  while(lo < n - 1 && x[lo + 1] == xm) lo++;
1336 
1337  return(lo);
1338 }
1339 
1340 /*****************************************************************************/
1341 /*
1342  * Move the roots that lie in the specified range [x0,x1) from newRoots to roots
1343  */
1344 static void
1345 keep_valid_roots(std::vector<double>& roots,
1346  std::vector<double>& newRoots, double x0, double x1)
1347 {
1348  for (unsigned int i = 0; i != newRoots.size(); ++i) {
1349  if(newRoots[i] >= x0 && newRoots[i] < x1) { // keep this root
1350  roots.push_back(newRoots[i]);
1351  }
1352  }
1353 
1354  newRoots.clear();
1355 }
1356 
1357 /*****************************************************************************/
1358 /*
1359  * find the real roots of a quadratic ax^2 + bx + c = 0
1360  */
1361 namespace {
1362 void
1363 do_quadratic(double a, double b, double c, std::vector<double> & roots)
1364 {
1365  if(::fabs(a) < std::numeric_limits<double>::epsilon()) {
1366  if(::fabs(b) >= std::numeric_limits<double>::epsilon()) {
1367  roots.push_back(-c/b);
1368  }
1369  } else {
1370  double const tmp = b*b - 4*a*c;
1371 
1372  if(tmp >= 0) {
1373  if (b >= 0) {
1374  roots.push_back((-b - sqrt(tmp))/(2*a));
1375  } else {
1376  roots.push_back((-b + sqrt(tmp))/(2*a));
1377  }
1378  roots.push_back(c/(a*roots[0]));
1379 /*
1380  * sort roots
1381  */
1382  if(roots[0] > roots[1]) {
1383  double const tmp2 = roots[0]; roots[0] = roots[1]; roots[1] = tmp2;
1384  }
1385  }
1386  }
1387 }
1388 
1389 /*****************************************************************************/
1390 /*
1391  * find the real roots of a cubic ax^3 + bx^2 + cx + d = 0
1392  */
1393 void
1394 do_cubic(double a, double b, double c, double d, std::vector<double> & roots)
1395 {
1396  if (::fabs(a) < std::numeric_limits<double>::epsilon()) {
1397  do_quadratic(b, c, d, roots);
1398  return;
1399  }
1400  b /= a; c /= a; d /= a;
1401 
1402  double const q = (b*b - 3*c)/9;
1403  double const r = (2*b*b*b - 9*b*c + 27*d)/54;
1404  /*
1405  * n.b. note that the test for the number of roots is carried out on the
1406  * same variables as are used in (e.g.) the acos, as it is possible for
1407  * r*r < q*q*q && r > sq*sq*sq due to rounding.
1408  */
1409  double const sq = (q >= 0) ? sqrt(q) : -sqrt(-q);
1410  double const sq3 = sq*sq*sq;
1411  if(::fabs(r) < sq3) { // three real roots
1412  double const theta = ::acos(r/sq3); // sq3 cannot be zero
1413 
1414  roots.push_back(-2*sq*cos(theta/3) - b/3);
1415  roots.push_back(-2*sq*cos((theta + afwGeom::TWOPI)/3) - b/3);
1416  roots.push_back(-2*sq*cos((theta - afwGeom::TWOPI)/3) - b/3);
1417  /*
1418  * sort roots
1419  */
1420  if(roots[0] > roots[1]) {
1421  std::swap(roots[0], roots[1]);
1422  }
1423  if(roots[1] > roots[2]) {
1424  std::swap(roots[1], roots[2]);
1425  }
1426  if(roots[0] > roots[1]) {
1427  std::swap(roots[0], roots[1]);
1428  }
1429 
1430  return;
1431  } else if(::fabs(r) == sq3) { /* no more than two real roots */
1432  double const aa = -((r < 0) ? -::pow(-r,1.0/3.0) : ::pow(r,1.0/3.0));
1433 
1434  if(::fabs(aa) < std::numeric_limits<double>::epsilon()) { /* degenerate case; one real root */
1435  roots.push_back(-b/3);
1436  return;
1437  } else {
1438  roots.push_back(2*aa - b/3);
1439  roots.push_back(-aa - b/3);
1440 
1441  if(roots[0] > roots[1]) {
1442  std::swap(roots[0], roots[1]);
1443  }
1444 
1445  return;
1446  }
1447  } else { /* only one real root */
1448  double tmp = ::sqrt(r*r - (q > 0 ? sq3*sq3 : -sq3*sq3));
1449  tmp = r + (r < 0 ? -tmp : tmp);
1450  double const aa = -((tmp < 0) ? -::pow(-tmp,1.0/3.0) : ::pow(tmp,1.0/3.0));
1451  double const bb = (fabs(aa) < std::numeric_limits<double>::epsilon()) ? 0 : q/aa;
1452 
1453  roots.push_back((aa + bb) - b/3);
1454 #if 0
1455  roots.push_back(-(aa + bb)/2 - b/3); // the real
1456  roots.push_back(::sqrt(3)/2*(aa - bb)); // and imaginary parts of the complex roots
1457 #endif
1458  return;
1459  }
1460 }
1461 }
1462 
1463 /*****************************************************************************/
1464 /*
1465  * <AUTO EXTRACT>
1466  *
1467  * Find the roots of
1468  * Spline - val = 0
1469  * in the range [x0, x1). Return a vector of all the roots found
1470  */
1471 std::vector<double>
1472 Spline::roots(double const value,
1473  double a,
1474  double const b
1475  ) const
1476 {
1477  /*
1478  * Strategy: we know that the interpolant has the form
1479  * val = coef[0][i] +dx*(coef[1][i] + dx*(coef[2][i]/2 + dx*coef[3][i]/6))
1480  * so we can use the usual analytic solution for a cubic. Note that the
1481  * cubic quoted above returns dx, the distance from the previous knot,
1482  * rather than x itself
1483  */
1484  std::vector<double> roots; /* the roots found */
1485  double x0 = a; // lower end of current range
1486  double const x1 = b;
1487  int const nknot = _knots.size();
1488 
1489  int i0 = search_array(x0, &_knots[0], nknot, -1);
1490  int const i1 = search_array(x1, &_knots[0], nknot, i0);
1491  assert (i1 >= i0 && i1 <= nknot - 1);
1492 
1493  std::vector<double> newRoots; // the roots we find in some interval
1494 /*
1495  * Deal with special case that x0 may be off one end or the other of
1496  * the array of knots.
1497  */
1498  if(i0 < 0) { /* off bottom */
1499  i0 = 0;
1500  do_cubic(_coeffs[3][i0]/6, _coeffs[2][i0]/2, _coeffs[1][i0], _coeffs[0][i0] - value, newRoots);
1501  //
1502  // Could use
1503  // std::transform(newRoots.begin(), newRoots.end(), newRoots.begin(),
1504  // std::tr1::bind(std::plus<double>(), _1, _knots[i0]));
1505  // but let's not
1506  //
1507  for (unsigned int j = 0; j != newRoots.size(); ++j) {
1508  newRoots[j] += _knots[i0];
1509  }
1510  keep_valid_roots(roots, newRoots, x0, _knots[i0]);
1511 
1512  x0 = _knots[i0];
1513  } else if(i0 >= nknot) { /* off top */
1514  i0 = nknot - 1;
1515  assert (i0 >= 0);
1516  do_cubic(_coeffs[3][i0]/6, _coeffs[2][i0]/2, _coeffs[1][i0], _coeffs[0][i0] - value, newRoots);
1517 
1518  for (unsigned int j = 0; j != newRoots.size(); ++j) {
1519  newRoots[j] += _knots[i0];
1520  }
1521  keep_valid_roots(roots, newRoots, x0, x1);
1522 
1523  return roots;
1524  }
1525 /*
1526  * OK, now search in main body of spline. Note that i1 may be nknot - 1, and
1527  * in any case the right hand limit of the last segment is at x1, not a knot
1528  */
1529  for (int i = i0;i <= i1;i++) {
1530  do_cubic(_coeffs[3][i]/6, _coeffs[2][i]/2, _coeffs[1][i], _coeffs[0][i] - value, newRoots);
1531 
1532  for (unsigned int j = 0; j != newRoots.size(); ++j) {
1533  newRoots[j] += _knots[i];
1534  }
1535  keep_valid_roots(roots, newRoots, ((i == i0) ? x0 : _knots[i]), ((i == i1) ? x1 : _knots[i + 1]));
1536  }
1537 
1538  return roots;
1539 }
1540 }}}}
int y
std::vector< std::vector< double > > _coeffs
Definition: Spline.h:33
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)
Definition: Spline.cc:710
void swap(ImageBase< PixelT > &a, ImageBase< PixelT > &b)
Definition: Image.cc:291
std::vector< double > roots(double const value, double const x0, double const x1) const
Definition: Spline.cc:1472
void interpolate(std::vector< double > const &x, std::vector< double > &y) const
Definition: Spline.cc:38
double const TWOPI
Definition: Angle.h:19
int const x0
Definition: saturated.cc:45
Matrix alpha
void derivative(std::vector< double > const &x, std::vector< double > &dydx) const
Definition: Spline.cc:73
TautSpline(std::vector< double > const &x, std::vector< double > const &y, double const gamma=0, Symmetry type=Unknown)
Definition: Spline.cc:178
std::vector< double > _knots
Definition: Spline.h:32
int x
double chisq
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
void _allocateSpline(int const nknot)
Definition: Spline.cc:25
afw::table::Key< double > b
void calculateTautSpline(std::vector< double > const &x, std::vector< double > const &y, double const gamma0)
Definition: Spline.cc:213
void calculateTautSplineEvenOdd(std::vector< double > const &x, std::vector< double > const &y, double const gamma0, bool even)
Definition: Spline.cc:1199