LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
Spline.cc
Go to the documentation of this file.
1 /*
2  * A number of non-trivial splines
3  *
4  * @note These should be merged into lsst::afw::math::Interpolate, but its current implementation
5  * (and to some degree interface) uses gsl explicitly
6  */
7 #include <limits>
8 
9 #include "boost/format.hpp"
12 #include "lsst/geom/Angle.h"
13 
14 namespace lsst {
15 namespace afw {
16 namespace math {
17 namespace detail {
18 
19 static int search_array(double z, double const *x, int n, int i);
20 
21 void Spline::_allocateSpline(int const nknot) {
22  _knots.resize(nknot);
23  _coeffs.resize(4);
24  for (unsigned int i = 0; i != _coeffs.size(); ++i) {
25  _coeffs[i].reserve(nknot);
26  }
27 }
28 
30  int const nknot = _knots.size();
31  int const n = x.size();
32 
33  y.resize(n); // may default-construct elements which is a little inefficient
34  /*
35  * For _knots[i] <= x <= _knots[i+1], the interpolant
36  * has the form
37  * val = _coeff[0][i] +dx*(_coeff[1][i] + dx*(_coeff[2][i]/2 + dx*_coeff[3][i]/6))
38  * with
39  * dx = x - knots[i]
40  */
41  int ind = -1; // no idea initially
42  for (int i = 0; i != n; ++i) {
43  ind = search_array(x[i], &_knots[0], nknot, ind);
44 
45  if (ind < 0) { // off bottom
46  ind = 0;
47  } else if (ind >= nknot) { // off top
48  ind = nknot - 1;
49  }
50 
51  double const dx = x[i] - _knots[ind];
52  y[i] = _coeffs[0][ind] +
53  dx * (_coeffs[1][ind] + dx * (_coeffs[2][ind] / 2 + dx * _coeffs[3][ind] / 6));
54  }
55 }
56 
58  int const nknot = _knots.size();
59  int const n = x.size();
60 
61  dydx.resize(n); // may default-construct elements which is a little inefficient
62  /*
63  * For _knots[i] <= x <= _knots[i+1], the * interpolant has the form
64  * val = _coeff[0][i] +dx*(_coeff[1][i] + dx*(_coeff[2][i]/2 + dx*_coeff[3][i]/6))
65  * with
66  * dx = x - knots[i]
67  * so the derivative is
68  * val = _coeff[1][i] + dx*(_coeff[2][i] + dx*_coeff[3][i]/2))
69  */
70 
71  int ind = -1; // no idea initially
72  for (int i = 0; i != n; ++i) {
73  ind = search_array(x[i], &_knots[0], nknot, ind);
74 
75  if (ind < 0) { // off bottom
76  ind = 0;
77  } else if (ind >= nknot) { // off top
78  ind = nknot - 1;
79  }
80 
81  double const dx = x[i] - _knots[ind];
82  dydx[i] = _coeffs[1][ind] + dx * (_coeffs[2][ind] + dx * _coeffs[3][ind] / 2);
83  }
84 }
85 
87  Symmetry type) {
88  if (x.size() != y.size()) {
90  (boost::format("TautSpline: x and y must have the same size; saw %d %d\n") %
91  x.size() % y.size())
92  .str());
93  }
94 
95  int const ntau = x.size(); /* size of tau and gtau, must be >= 2*/
96  if (ntau < 2) {
98  (boost::format("TautSpline: ntau = %d, should be >= 2\n") % ntau).str());
99  }
100 
101  switch (type) {
102  case Unknown:
103  calculateTautSpline(x, y, gamma0);
104  break;
105  case Even:
106  case Odd:
107  calculateTautSplineEvenOdd(x, y, gamma0, type == Even);
108  break;
109  }
110 }
111 
112 void TautSpline::calculateTautSpline(std::vector<double> const &x, std::vector<double> const &y,
113  double const gamma0) {
114  const double *tau = &x[0];
115  const double *gtau = &y[0];
116  int const ntau = x.size(); // size of tau and gtau, must be >= 2
117 
118  if (ntau < 4) { // use a single quadratic
119  int const nknot = ntau;
120 
121  _allocateSpline(nknot);
122 
123  _knots[0] = tau[0];
124  for (int i = 1; i < nknot; i++) {
125  _knots[i] = tau[i];
126  if (tau[i - 1] >= tau[i]) {
128  (boost::format("point %d and the next, %f %f, are out of order") % (i - 1) %
129  tau[i - 1] % tau[i])
130  .str());
131  }
132  }
133 
134  if (ntau == 2) {
135  _coeffs[0][0] = gtau[0];
136  _coeffs[1][0] = (gtau[1] - gtau[0]) / (tau[1] - tau[0]);
137  _coeffs[2][0] = _coeffs[3][0] = 0;
138 
139  _coeffs[0][1] = gtau[1];
140  _coeffs[1][1] = (gtau[1] - gtau[0]) / (tau[1] - tau[0]);
141  _coeffs[2][1] = _coeffs[3][1] = 0;
142  } else { /* must be 3 */
143  double tmp = (tau[2] - tau[0]) * (tau[2] - tau[1]) * (tau[1] - tau[0]);
144  _coeffs[0][0] = gtau[0];
145  _coeffs[1][0] = ((gtau[1] - gtau[0]) * pow(tau[2] - tau[0], 2) -
146  (gtau[2] - gtau[0]) * pow(tau[1] - tau[0], 2)) /
147  tmp;
148  _coeffs[2][0] =
149  -2 * ((gtau[1] - gtau[0]) * (tau[2] - tau[0]) - (gtau[2] - gtau[0]) * (tau[1] - tau[0])) /
150  tmp;
151  _coeffs[3][0] = 0;
152 
153  _coeffs[0][1] = gtau[1];
154  _coeffs[1][1] = _coeffs[1][0] + (tau[1] - tau[0]) * _coeffs[2][0];
155  _coeffs[2][1] = _coeffs[2][0];
156  _coeffs[3][1] = 0;
157 
158  _coeffs[0][2] = gtau[2];
159  _coeffs[1][2] = _coeffs[1][0] + (tau[2] - tau[0]) * _coeffs[2][0];
160  _coeffs[2][2] = _coeffs[2][0];
161  _coeffs[3][2] = 0;
162  }
163 
164  return;
165  }
166  /*
167  * Allocate scratch space
168  * s[0][...] = dtau = tau(.+1) - tau
169  * s[1][...] = diag = diagonal in linear system
170  * s[2][...] = u = upper diagonal in linear system
171  * s[3][...] = r = right side for linear system (initially)
172  * = fsecnd = solution of linear system , namely the second
173  * derivatives of interpolant at tau
174  * s[4][...] = z = indicator of additional knots
175  * s[5][...] = 1/hsecnd(1,x) with x = z or = 1-z. see below.
176  */
177  std::vector<std::vector<double> > s(6); // scratch space
178 
179  for (int i = 0; i != 6; i++) {
180  s[i].resize(ntau);
181  }
182  /*
183  * Construct delta tau and first and second (divided) differences of data
184  */
185 
186  for (int i = 0; i < ntau - 1; i++) {
187  s[0][i] = tau[i + 1] - tau[i];
188  if (s[0][i] <= 0.) {
190  (boost::format("point %d and the next, %f %f, are out of order") % i % tau[i] %
191  tau[i + 1])
192  .str());
193  }
194  s[3][i + 1] = (gtau[i + 1] - gtau[i]) / s[0][i];
195  }
196  for (int i = 1; i < ntau - 1; ++i) {
197  s[3][i] = s[3][i + 1] - s[3][i];
198  }
199  /*
200  * Construct system of equations for second derivatives at tau. At each
201  * interior data point, there is one continuity equation, at the first
202  * and the last interior data point there is an additional one for a
203  * total of ntau equations in ntau unknowns.
204  */
205  s[1][1] = s[0][0] / 3;
206 
207  int method;
208  double gamma = gamma0; // control the smoothing
209  if (gamma <= 0) {
210  method = 1;
211  } else if (gamma > 3) {
212  gamma -= 3;
213  method = 3;
214  } else {
215  method = 2;
216  }
217  double const onemg3 = 1 - gamma / 3;
218 
219  int nknot = ntau; // count number of knots
220  /*
221  * Some compilers don't realise that the flow of control always initialises
222  * these variables; so initialise them to placate e.g. gcc
223  */
224  double zeta = 0.0;
225  double alpha = 0.0;
226  double ratio = 0.0;
227 
228  // double c, d;
229  // double z, denom, factr2;
230  // double onemzt, zt2, del;
231 
232  double entry3 = 0.0;
233  double factor = 0.0;
234  double onemzt = 0;
235  double zt2 = 0;
236  double z_half = 0;
237 
238  for (int i = 1; i < ntau - 2; ++i) {
239  /*
240  * construct z[i] and zeta[i]
241  */
242  double z = .5;
243  if ((method == 2 && s[3][i] * s[3][i + 1] >= 0) || method == 3) {
244  double const temp = fabs(s[3][i + 1]);
245  double const denom = fabs(s[3][i]) + temp;
246  if (denom != 0) {
247  z = temp / denom;
248  if (fabs(z - 0.5) <= 1.0 / 6.0) {
249  z = 0.5;
250  }
251  }
252  }
253 
254  s[4][i] = z;
255  /*
256  set up part of the i-th equation which depends on the i-th interval
257  */
258  z_half = z - 0.5;
259  if (z_half < 0) {
260  zeta = gamma * z;
261  onemzt = 1 - zeta;
262  zt2 = zeta * zeta;
263  double temp = onemg3 / onemzt;
264  alpha = (temp < 1 ? temp : 1);
265  factor = zeta / (alpha * (zt2 - 1) + 1);
266  s[5][i] = zeta * factor / 6;
267  s[1][i] += s[0][i] * ((1 - alpha * onemzt) * factor / 2 - s[5][i]);
268  /*
269  * if z = 0 and the previous z = 1, then d[i] = 0. Since then
270  * also u[i-1] = l[i+1] = 0, its value does not matter. Reset
271  * d[i] = 1 to insure nonzero pivot in elimination.
272  */
273  if (s[1][i] <= 0.) {
274  s[1][i] = 1;
275  }
276  s[2][i] = s[0][i] / 6;
277 
278  if (z != 0) { /* we'll get a new knot */
279  nknot++;
280  }
281  } else if (z_half == 0) {
282  s[1][i] += s[0][i] / 3;
283  s[2][i] = s[0][i] / 6;
284  } else {
285  onemzt = gamma * (1 - z);
286  zeta = 1 - onemzt;
287  double const temp = onemg3 / zeta;
288  alpha = (temp < 1 ? temp : 1);
289  factor = onemzt / (1 - alpha * zeta * (onemzt + 1));
290  s[5][i] = onemzt * factor / 6;
291  s[1][i] += s[0][i] / 3;
292  s[2][i] = s[5][i] * s[0][i];
293 
294  if (onemzt != 0) { /* we'll get a new knot */
295  nknot++;
296  }
297  }
298  if (i == 1) {
299  s[4][0] = 0.5;
300  /*
301  * the first two equations enforce continuity of the first and of
302  * the third derivative across tau[1].
303  */
304  s[1][0] = s[0][0] / 6;
305  s[2][0] = s[1][1];
306  entry3 = s[2][1];
307  if (z_half < 0) {
308  const double factr2 = zeta * (alpha * (zt2 - 1.) + 1.) / (alpha * (zeta * zt2 - 1.) + 1.);
309  ratio = factr2 * s[0][1] / s[1][0];
310  s[1][1] = factr2 * s[0][1] + s[0][0];
311  s[2][1] = -factr2 * s[0][0];
312  } else if (z_half == 0) {
313  ratio = s[0][1] / s[1][0];
314  s[1][1] = s[0][1] + s[0][0];
315  s[2][1] = -s[0][0];
316  } else {
317  ratio = s[0][1] / s[1][0];
318  s[1][1] = s[0][1] + s[0][0];
319  s[2][1] = -s[0][0] * 6 * alpha * s[5][1];
320  }
321  /*
322  * at this point, the first two equations read
323  * diag[0]*x0 + u[0]*x1 + entry3*x2 = r[1]
324  * -ratio*diag[0]*x0 + diag[1]*x1 + u[1]*x2 = 0
325  * set r[0] = r[1] and eliminate x1 from the second equation
326  */
327  s[3][0] = s[3][1];
328 
329  s[1][1] += ratio * s[2][0];
330  s[2][1] += ratio * entry3;
331  s[3][1] = ratio * s[3][1];
332  } else {
333  /*
334  * the i-th equation enforces continuity of the first derivative
335  * across tau[i]; it reads
336  * -ratio*diag[i-1]*x_{i-1} + diag[i]*x_i + u[i]*x_{i+1} = r[i]
337  * eliminate x_{i-1} from this equation
338  */
339  s[1][i] += ratio * s[2][i - 1];
340  s[3][i] += ratio * s[3][i - 1];
341  }
342  /*
343  * Set up the part of the next equation which depends on the i-th interval.
344  */
345  if (z_half < 0) {
346  ratio = -s[5][i] * s[0][i] / s[1][i];
347  s[1][i + 1] = s[0][i] / 3;
348  } else if (z_half == 0) {
349  ratio = -(s[0][i] / 6) / s[1][i];
350  s[1][i + 1] = s[0][i] / 3;
351  } else {
352  ratio = -(s[0][i] / 6) / s[1][i];
353  s[1][i + 1] = s[0][i] * ((1 - zeta * alpha) * factor / 2 - s[5][i]);
354  }
355  }
356 
357  s[4][ntau - 2] = 0.5;
358  /*
359  * last two equations, which enforce continuity of third derivative and
360  * of first derivative across tau[ntau - 2]
361  */
362  double const entry_ = ratio * s[2][ntau - 3] + s[1][ntau - 2] + s[0][ntau - 2] / 3;
363  s[1][ntau - 1] = s[0][ntau - 2] / 6;
364  s[3][ntau - 1] = ratio * s[3][ntau - 3] + s[3][ntau - 2];
365  if (z_half < 0) {
366  ratio = s[0][ntau - 2] * 6 * s[5][ntau - 3] * alpha / s[1][ntau - 3];
367  s[1][ntau - 2] = ratio * s[2][ntau - 3] + s[0][ntau - 2] + s[0][ntau - 3];
368  s[2][ntau - 2] = -s[0][ntau - 3];
369  } else if (z_half == 0) {
370  ratio = s[0][ntau - 2] / s[1][ntau - 3];
371  s[1][ntau - 2] = ratio * s[2][ntau - 3] + s[0][ntau - 2] + s[0][ntau - 3];
372  s[2][ntau - 2] = -s[0][ntau - 3];
373  } else {
374  const double factr2 =
375  onemzt * (alpha * (onemzt * onemzt - 1) + 1) / (alpha * (onemzt * onemzt * onemzt - 1) + 1);
376  ratio = factr2 * s[0][ntau - 2] / s[1][ntau - 3];
377  s[1][ntau - 2] = ratio * s[2][ntau - 3] + factr2 * s[0][ntau - 3] + s[0][ntau - 2];
378  s[2][ntau - 2] = -factr2 * s[0][ntau - 3];
379  }
380  /*
381  * at this point, the last two equations read
382  * diag[i]*x_i + u[i]*x_{i+1} = r[i]
383  * -ratio*diag[i]*x_i + diag[i+1]*x_{i+1} = r[i+1]
384  * eliminate x_i from last equation
385  */
386  s[3][ntau - 2] = ratio * s[3][ntau - 3];
387  ratio = -entry_ / s[1][ntau - 2];
388  s[1][ntau - 1] += ratio * s[2][ntau - 2];
389  s[3][ntau - 1] += ratio * s[3][ntau - 2];
390 
391  /*
392  * back substitution
393  */
394  s[3][ntau - 1] /= s[1][ntau - 1];
395  for (int i = ntau - 2; i > 0; --i) {
396  s[3][i] = (s[3][i] - s[2][i] * s[3][i + 1]) / s[1][i];
397  }
398 
399  s[3][0] = (s[3][0] - s[2][0] * s[3][1] - entry3 * s[3][2]) / s[1][0];
400 /*
401  * construct polynomial pieces; first allocate space for the coefficients
402  */
403 #if 1
404  /*
405  * Start by counting the knots
406  */
407  {
408  int const nknot0 = nknot;
409  int nknot = ntau;
410 
411  for (int i = 0; i < ntau - 1; ++i) {
412  double const z = s[4][i];
413  if ((z < 0.5 && z != 0.0) || (z > 0.5 && (1 - z) != 0.0)) {
414  nknot++;
415  }
416  }
417  assert(nknot == nknot0);
418  }
419 #endif
420  _allocateSpline(nknot);
421 
422  _knots[0] = tau[0];
423  int j = 0;
424  for (int i = 0; i < ntau - 1; ++i) {
425  _coeffs[0][j] = gtau[i];
426  _coeffs[2][j] = s[3][i];
427  double const divdif = (gtau[i + 1] - gtau[i]) / s[0][i];
428  double z = s[4][i];
429  double const z_half = z - 0.5;
430  if (z_half < 0) {
431  if (z == 0) {
432  _coeffs[1][j] = divdif;
433  _coeffs[2][j] = 0;
434  _coeffs[3][j] = 0;
435  } else {
436  zeta = gamma * z;
437  onemzt = 1 - zeta;
438  double const c = s[3][i + 1] / 6;
439  double const d = s[3][i] * s[5][i];
440  j++;
441 
442  double const del = zeta * s[0][i];
443  _knots[j] = tau[i] + del;
444  zt2 = zeta * zeta;
445  double temp = onemg3 / onemzt;
446  alpha = (temp < 1 ? temp : 1);
447  factor = onemzt * onemzt * alpha;
448  temp = s[0][i];
449  _coeffs[0][j] = gtau[i] + divdif * del +
450  temp * temp * (d * onemzt * (factor - 1) + c * zeta * (zt2 - 1));
451  _coeffs[1][j] = divdif + s[0][i] * (d * (1 - 3 * factor) + c * (3 * zt2 - 1));
452  _coeffs[2][j] = (d * alpha * onemzt + c * zeta) * 6;
453  _coeffs[3][j] = (c - d * alpha) * 6 / s[0][i];
454  if (del * zt2 == 0) {
455  _coeffs[1][j - 1] = 0; /* would be NaN in an */
456  _coeffs[3][j - 1] = 0; /* 0-length interval */
457  } else {
458  _coeffs[3][j - 1] = _coeffs[3][j] - d * 6 * (1 - alpha) / (del * zt2);
459  _coeffs[1][j - 1] = _coeffs[1][j] - del * (_coeffs[2][j] - del / 2 * _coeffs[3][j - 1]);
460  }
461  }
462  } else if (z_half == 0) {
463  _coeffs[1][j] = divdif - s[0][i] * (s[3][i] * 2 + s[3][i + 1]) / 6;
464  _coeffs[3][j] = (s[3][i + 1] - s[3][i]) / s[0][i];
465  } else {
466  onemzt = gamma * (1 - z);
467  if (onemzt == 0) {
468  _coeffs[1][j] = divdif;
469  _coeffs[2][j] = 0;
470  _coeffs[3][j] = 0;
471  } else {
472  zeta = 1 - onemzt;
473  double const temp = onemg3 / zeta;
474  alpha = (temp < 1 ? temp : 1);
475  double const c = s[3][i + 1] * s[5][i];
476  double const d = s[3][i] / 6;
477  double const del = zeta * s[0][i];
478  _knots[j + 1] = tau[i] + del;
479  _coeffs[1][j] = divdif - s[0][i] * (2 * d + c);
480  _coeffs[3][j] = (c * alpha - d) * 6 / s[0][i];
481  j++;
482 
483  _coeffs[3][j] =
484  _coeffs[3][j - 1] + (1 - alpha) * 6 * c / (s[0][i] * (onemzt * onemzt * onemzt));
485  _coeffs[2][j] = _coeffs[2][j - 1] + del * _coeffs[3][j - 1];
486  _coeffs[1][j] = _coeffs[1][j - 1] + del * (_coeffs[2][j - 1] + del / 2 * _coeffs[3][j - 1]);
487  _coeffs[0][j] = _coeffs[0][j - 1] +
488  del * (_coeffs[1][j - 1] +
489  del / 2 * (_coeffs[2][j - 1] + del / 3 * _coeffs[3][j - 1]));
490  }
491  }
492 
493  j++;
494  _knots[j] = tau[i + 1];
495  }
496  /*
497  * If there are discontinuities some of the knots may be at the same
498  * position; in this case we generated some NaNs above. As they only
499  * occur for 0-length segments, it's safe to replace them by 0s
500  *
501  * Due to the not-a-knot condition, the last set of coefficients isn't
502  * needed (the last-but-one is equivalent), but it makes the book-keeping
503  * easier if we _do_ generate them
504  */
505  double const del = tau[ntau - 1] - _knots[nknot - 2];
506 
507  _coeffs[0][nknot - 1] = _coeffs[0][nknot - 2] +
508  del * (_coeffs[1][nknot - 2] +
509  del * (_coeffs[2][nknot - 2] / 2 + del * _coeffs[3][nknot - 2] / 6));
510  _coeffs[1][nknot - 1] =
511  _coeffs[1][nknot - 2] + del * (_coeffs[2][nknot - 2] + del * _coeffs[3][nknot - 2] / 2);
512  _coeffs[2][nknot - 1] = _coeffs[2][nknot - 2] + del * _coeffs[3][nknot - 2];
513  _coeffs[3][nknot - 1] = _coeffs[3][nknot - 2];
514 
515  assert(j + 1 == nknot);
516 }
517 
518 /*
519  * Here's the code to fit smoothing splines through data points
520  */
521 static void spcof1(const double x[], double avh, const double y[], const double dy[], int n, double p,
522  double q, double a[], double *c[3], double u[], const double v[]);
523 
524 static void sperr1(const double x[], double avh, const double dy[], int n, double *r[3], double p, double var,
525  std::vector<double> *se);
526 
527 static double spfit1(const double x[], double avh, const double dy[], int n, double rho, double *p, double *q,
528  double var, double stat[], const double a[], double *c[3], double *r[3], double *t[2],
529  double u[], double v[]);
530 
531 static double spint1(const double x[], double *avh, const double y[], double dy[], int n, double a[],
532  double *c[3], double *r[3], double *t[2]);
533 
535  std::vector<double> const &df, double s, double *chisq,
536  std::vector<double> *errs) {
537  float var = 1; // i.e. df is the absolute s.d. N.B. ADD GCV Variant with var=-1
538  int const n = x.size();
539  double const ratio = 2.0;
540  double const tau = 1.618033989; /* golden ratio */
541  double avdf, avar, stat[6];
542  double p, q, delta, r1, r2, r3, r4;
543  double gf1, gf2, gf3, gf4, avh, err;
544  /*
545  * allocate scratch space
546  */
547  _allocateSpline(n);
548 
549  double *y = &_coeffs[0][0];
550  double *c[3];
551  c[0] = &_coeffs[1][0];
552  c[1] = &_coeffs[2][0];
553  c[2] = &_coeffs[3][0];
554 
555  std::vector<double> scratch(7 * (n + 2)); // scratch space
556 
557  double *r[3];
558  r[0] = &scratch[0] + 1; // we want indices -1..n
559  r[1] = r[0] + (n + 2);
560  r[2] = r[1] + (n + 2);
561  double *t[2];
562  t[0] = r[2] + (n + 2);
563  t[1] = t[0] + (n + 2);
564  double *u = t[1] + (n + 2);
565  double *v = u + (n + 2);
566  /*
567  * and so to work.
568  */
569  std::vector<double> sdf = df; // scaled values of df
570 
571  avdf = spint1(&x[0], &avh, &f[0], &sdf[0], n, y, c, r, t);
572  avar = var;
573  if (var > 0) {
574  avar *= avdf * avdf;
575  }
576 
577  if (var == 0) { /* simply find a natural cubic spline*/
578  r1 = 0;
579 
580  gf1 = spfit1(&x[0], avh, &sdf[0], n, r1, &p, &q, avar, stat, y, c, r, t, u, v);
581  } else { /* Find local minimum of gcv or the
582  expected mean square error */
583  r1 = 1;
584  r2 = ratio * r1;
585  gf2 = spfit1(&x[0], avh, &sdf[0], n, r2, &p, &q, avar, stat, y, c, r, t, u, v);
586  bool set_r3 = false; // was r3 set?
587  for (;;) {
588  gf1 = spfit1(&x[0], avh, &sdf[0], n, r1, &p, &q, avar, stat, y, c, r, t, u, v);
589  if (gf1 > gf2) {
590  break;
591  }
592 
593  if (p <= 0) {
594  break;
595  }
596  r2 = r1;
597  gf2 = gf1;
598  r1 /= ratio;
599  }
600 
601  if (p <= 0) {
602  set_r3 = false;
603  r3 = 0; /* placate compiler */
604  } else {
605  r3 = ratio * r2;
606  set_r3 = true;
607 
608  for (;;) {
609  gf3 = spfit1(&x[0], avh, &sdf[0], n, r3, &p, &q, avar, stat, y, c, r, t, u, v);
610  if (gf3 >= gf2) {
611  break;
612  }
613 
614  if (q <= 0) {
615  break;
616  }
617  r2 = r3;
618  gf2 = gf3;
619  r3 = ratio * r3;
620  }
621  }
622 
623  if (p > 0 && q > 0) {
624  assert(set_r3);
625  r2 = r3;
626  gf2 = gf3;
627  delta = (r2 - r1) / tau;
628  r4 = r1 + delta;
629  r3 = r2 - delta;
630  gf3 = spfit1(&x[0], avh, &sdf[0], n, r3, &p, &q, avar, stat, y, c, r, t, u, v);
631  gf4 = spfit1(&x[0], avh, &sdf[0], n, r4, &p, &q, avar, stat, y, c, r, t, u, v);
632  /*
633  * Golden section search for local minimum
634  */
635  do {
636  if (gf3 <= gf4) {
637  r2 = r4;
638  gf2 = gf4;
639  r4 = r3;
640  gf4 = gf3;
641  delta /= tau;
642  r3 = r2 - delta;
643  gf3 = spfit1(&x[0], avh, &sdf[0], n, r3, &p, &q, avar, stat, y, c, r, t, u, v);
644  } else {
645  r1 = r3;
646  gf1 = gf3;
647  r3 = r4;
648  gf3 = gf4;
649  delta /= tau;
650  r4 = r1 + delta;
651  gf4 = spfit1(&x[0], avh, &sdf[0], n, r4, &p, &q, avar, stat, y, c, r, t, u, v);
652  }
653 
654  err = (r2 - r1) / (r1 + r2);
655  } while (err * err + 1 > 1 && err > 1e-6);
656 
657  r1 = (r1 + r2) * .5;
658  gf1 = spfit1(&x[0], avh, &sdf[0], n, r1, &p, &q, avar, stat, y, c, r, t, u, v);
659  }
660  }
661  /*
662  * Calculate spline coefficients
663  */
664  spcof1(&x[0], avh, &f[0], &sdf[0], n, p, q, y, c, u, v);
665 
666  stat[2] /= avdf * avdf; /* undo scaling */
667  stat[3] /= avdf * avdf;
668  stat[4] /= avdf * avdf;
669  /*
670  * Optionally calculate standard error estimates
671  */
672  if (errs != NULL) {
673  sperr1(&x[0], avh, &sdf[0], n, r, p, avar, errs);
674  }
675  /*
676  * clean up
677  */
678  if (chisq != NULL) {
679  *chisq = n * stat[4];
680  }
681 }
682 
692 static double spint1(const double x[], double *avh, const double f[], double df[], int n, double a[],
693  double *c[3], double *r[3], double *t[2]) {
694  double avdf;
695  double e, ff, g, h;
696  int i;
697 
698  assert(n >= 3);
699  /*
700  * Get average x spacing in avh
701  */
702  g = 0;
703  for (i = 0; i < n - 1; ++i) {
704  h = x[i + 1] - x[i];
705  assert(h > 0);
706 
707  g += h;
708  }
709  *avh = g / (n - 1);
710  /*
711  * Scale relative weights
712  */
713  g = 0;
714  for (int i = 0; i < n; ++i) {
715  assert(df[i] > 0);
716 
717  g += df[i] * df[i];
718  }
719  avdf = sqrt(g / n);
720 
721  for (i = 0; i < n; ++i) {
722  df[i] /= avdf;
723  }
724  /*
725  * Initialize h,f
726  */
727  h = (x[1] - x[0]) / (*avh);
728  ff = (f[1] - f[0]) / h;
729  /*
730  * Calculate a,t,r
731  */
732  for (i = 1; i < n - 1; ++i) {
733  g = h;
734  h = (x[i + 1] - x[i]) / (*avh);
735  e = ff;
736  ff = (f[i + 1] - f[i]) / h;
737  a[i] = ff - e;
738  t[0][i] = (g + h) * 2. / 3.;
739  t[1][i] = h / 3.;
740  r[2][i] = df[i - 1] / g;
741  r[0][i] = df[i + 1] / h;
742  r[1][i] = -df[i] / g - df[i] / h;
743  }
744  /*
745  * Calculate c = r'*r
746  */
747  r[1][n - 1] = 0;
748  r[2][n - 1] = 0;
749  r[2][n] = 0;
750 
751  for (i = 1; i < n - 1; i++) {
752  c[0][i] = r[0][i] * r[0][i] + r[1][i] * r[1][i] + r[2][i] * r[2][i];
753  c[1][i] = r[0][i] * r[1][i + 1] + r[1][i] * r[2][i + 1];
754  c[2][i] = r[0][i] * r[2][i + 2];
755  }
756 
757  return (avdf);
758 }
759 
785 static double spfit1(const double x[], double avh, const double dy[], int n, double rho, double *pp,
786  double *pq, double var, double stat[], const double a[], double *c[3], double *r[3],
787  double *t[2], double u[], double v[]) {
788  double const eps = std::numeric_limits<double>::epsilon();
789  double e, f, g, h;
790  double fun;
791  int i;
792  double p, q; /* == *pp, *pq */
793  /*
794  * Use p and q instead of rho to prevent overflow or underflow
795  */
796  if (fabs(rho) < eps) {
797  p = 0;
798  q = 1;
799  } else if (fabs(1 / rho) < eps) {
800  p = 1;
801  q = 0;
802  } else {
803  p = rho / (1 + rho);
804  q = 1 / (1 + rho);
805  }
806  /*
807  * Rational Cholesky decomposition of p*c + q*t
808  */
809  f = 0;
810  g = 0;
811  h = 0;
812  r[0][-1] = r[0][0] = 0;
813 
814  for (int i = 1; i < n - 1; ++i) {
815  r[2][i - 2] = g * r[0][i - 2];
816  r[1][i - 1] = f * r[0][i - 1];
817  {
818  double tmp = p * c[0][i] + q * t[0][i] - f * r[1][i - 1] - g * r[2][i - 2];
819  if (tmp == 0.0) {
820  r[0][i] = 1e30;
821  } else {
822  r[0][i] = 1 / tmp;
823  }
824  }
825  f = p * c[1][i] + q * t[1][i] - h * r[1][i - 1];
826  g = h;
827  h = p * c[2][i];
828  }
829  /*
830  * Solve for u
831  */
832  u[-1] = u[0] = 0;
833  for (int i = 1; i < n - 1; i++) {
834  u[i] = a[i] - r[1][i - 1] * u[i - 1] - r[2][i - 2] * u[i - 2];
835  }
836  u[n - 1] = u[n] = 0;
837  for (int i = n - 2; i > 0; i--) {
838  u[i] = r[0][i] * u[i] - r[1][i] * u[i + 1] - r[2][i] * u[i + 2];
839  }
840  /*
841  * Calculate residual vector v
842  */
843  e = h = 0;
844  for (int i = 0; i < n - 1; i++) {
845  g = h;
846  h = (u[i + 1] - u[i]) / ((x[i + 1] - x[i]) / avh);
847  v[i] = dy[i] * (h - g);
848  e += v[i] * v[i];
849  }
850  v[n - 1] = -h * dy[n - 1];
851  e += v[n - 1] * v[n - 1];
852  /*
853  * Calculate upper three bands of inverse matrix
854  */
855  r[0][n - 1] = r[1][n - 1] = r[0][n] = 0;
856  for (i = n - 2; i > 0; i--) {
857  g = r[1][i];
858  h = r[2][i];
859  r[1][i] = -g * r[0][i + 1] - h * r[1][i + 1];
860  r[2][i] = -g * r[1][i + 1] - h * r[0][i + 2];
861  r[0][i] = r[0][i] - g * r[1][i] - h * r[2][i];
862  }
863  /*
864  * Calculate trace
865  */
866  f = g = h = 0;
867  for (i = 1; i < n - 1; ++i) {
868  f += r[0][i] * c[0][i];
869  g += r[1][i] * c[1][i];
870  h += r[2][i] * c[2][i];
871  }
872  f += 2 * (g + h);
873  /*
874  * Calculate statistics
875  */
876  stat[0] = p;
877  stat[1] = f * p;
878  stat[2] = n * e / (f * f + 1e-20);
879  stat[3] = e * p * p / n;
880  stat[5] = e * p / ((f < 0) ? f - 1e-10 : f + 1e-10);
881  if (var >= 0) {
882  fun = stat[3] - 2 * var * stat[1] / n + var;
883 
884  stat[4] = fun;
885  } else {
886  stat[4] = stat[5] - stat[3];
887  fun = stat[2];
888  }
889 
890  *pp = p;
891  *pq = q;
892 
893  return (fun);
894 }
895 
900 static void spcof1(const double x[], double avh, const double y[], const double dy[], int n, double p,
901  double q, double a[], double *c[3], double u[], const double v[]) {
902  double h;
903  int i;
904  double qh;
905 
906  qh = q / (avh * avh);
907 
908  for (i = 0; i < n; ++i) {
909  a[i] = y[i] - p * dy[i] * v[i];
910  u[i] = qh * u[i];
911  }
912  /*
913  * calculate c
914  */
915  for (i = 0; i < n - 1; ++i) {
916  h = x[i + 1] - x[i];
917  c[2][i] = (u[i + 1] - u[i]) / (3 * h);
918  c[0][i] = (a[i + 1] - a[i]) / h - (h * c[2][i] + u[i]) * h;
919  c[1][i] = u[i];
920  }
921 }
922 
928 static void sperr1(const double x[], double avh, const double dy[], int n, double *r[3], double p, double var,
929  std::vector<double> *se) {
930  double f, g, h;
931  int i;
932  double f1, g1, h1;
933  /*
934  * Initialize
935  */
936  h = avh / (x[1] - x[0]);
937  (*se)[0] = 1 - p * dy[0] * dy[0] * h * h * r[0][1];
938  r[0][0] = r[1][0] = r[2][0] = 0;
939  /*
940  * Calculate diagonal elements
941  */
942  for (i = 1; i < n - 1; ++i) {
943  f = h;
944  h = avh / (x[i + 1] - x[i]);
945  g = -(f + h);
946  f1 = f * r[0][i - 1] + g * r[1][i - 1] + h * r[2][i - 1];
947  g1 = f * r[1][i - 1] + g * r[0][i] + h * r[1][i];
948  h1 = f * r[2][i - 1] + g * r[1][i] + h * r[0][i + 1];
949  (*se)[i] = 1 - p * dy[i] * dy[i] * (f * f1 + g * g1 + h * h1);
950  }
951  (*se)[n - 1] = 1 - p * dy[n - 1] * dy[n - 1] * h * h * r[0][n - 2];
952  /*
953  * Calculate standard error estimates
954  */
955  for (int i = 0; i < n; ++i) {
956  double const tmp = (*se)[i] * var;
957  (*se)[i] = (tmp >= 0) ? sqrt(tmp) * dy[i] : 0;
958  }
959 }
960 
961 void TautSpline::calculateTautSplineEvenOdd(std::vector<double> const &_tau, std::vector<double> const &_gtau,
962  double const gamma,
963  bool const even // ensure Even symmetry
964 ) {
965  const double *tau = &_tau[0];
966  const double *gtau = &_gtau[0];
967  int const ntau = _tau.size(); // size of tau and gtau, must be >= 2
968  std::vector<double> x, y; // tau and gtau, extended to -ve tau
969 
970  if (tau[0] == 0.0f) {
971  int const np = 2 * ntau - 1;
972  x.resize(np);
973  y.resize(np);
974 
975  x[ntau - 1] = tau[0];
976  y[ntau - 1] = gtau[0];
977  for (int i = 1; i != ntau; ++i) {
978  if (even) {
979  x[ntau - 1 + i] = tau[i];
980  y[ntau - 1 + i] = gtau[i];
981  x[ntau - 1 - i] = -tau[i];
982  y[ntau - 1 - i] = gtau[i];
983  } else {
984  x[ntau - 1 + i] = tau[i];
985  y[ntau - 1 + i] = gtau[i];
986  x[ntau - 1 - i] = -tau[i];
987  y[ntau - 1 - i] = -gtau[i];
988  }
989  }
990  } else {
991  int const np = 2 * ntau;
992  x.resize(np);
993  y.resize(np);
994 
995  for (int i = 0; i != ntau; ++i) {
996  if (even) {
997  x[ntau + i] = tau[i];
998  y[ntau + i] = gtau[i];
999  x[ntau - 1 - i] = -tau[i];
1000  y[ntau - 1 - i] = gtau[i];
1001  } else {
1002  x[ntau + i] = tau[i];
1003  y[ntau + i] = gtau[i];
1004  x[ntau - 1 - i] = -tau[i];
1005  y[ntau - 1 - i] = -gtau[i];
1006  }
1007  }
1008  }
1009 
1010  TautSpline sp(x, y, gamma); // fit a taut spline to x, y
1011  /*
1012  * Now repackage that spline to reflect the original points
1013  */
1014  int ii;
1015  for (ii = sp._knots.size() - 1; ii >= 0; --ii) {
1016  if (sp._knots[ii] < 0.0f) {
1017  break;
1018  }
1019  }
1020  int const i0 = ii + 1;
1021  int const nknot = sp._knots.size() - i0;
1022 
1023  _allocateSpline(nknot);
1024 
1025  for (int i = i0; i != static_cast<int>(sp._knots.size()); ++i) {
1026  _knots[i - i0] = sp._knots[i];
1027  for (int j = 0; j != 4; ++j) {
1028  _coeffs[j][i - i0] = sp._coeffs[j][i];
1029  }
1030  }
1031 }
1032 
1038 static int search_array(double z, double const *x, int n, int i) {
1039  int lo, hi, mid;
1040  double xm;
1041 
1042  if (i < 0 || i >= n) { /* initial guess is useless */
1043  lo = -1;
1044  hi = n;
1045  } else {
1046  unsigned int step = 1; /* how much to step up/down */
1047 
1048  if (z > x[i]) { /* expand search upwards */
1049  if (i == n - 1) { /* off top of array */
1050  return (n - 1);
1051  }
1052 
1053  lo = i;
1054  hi = lo + 1;
1055  while (z >= x[hi]) {
1056  lo = hi;
1057  step += step; /* double step size */
1058  hi = lo + step;
1059  if (hi >= n) { /* reached top of array */
1060  hi = n - 1;
1061  break;
1062  }
1063  }
1064  } else { /* expand it downwards */
1065  if (i == 0) { /* off bottom of array */
1066  return (-1);
1067  }
1068 
1069  hi = i;
1070  lo = i - 1;
1071  while (z < x[lo]) {
1072  hi = lo;
1073  step += step; /* double step size */
1074  lo = hi - step;
1075  if (lo < 0) { /* off bottom of array */
1076  lo = -1;
1077  break;
1078  }
1079  }
1080  }
1081  }
1082  /*
1083  * perform bisection
1084  */
1085  while (hi - lo > 1) {
1086  mid = (lo + hi) / 2;
1087  xm = x[mid];
1088  if (z <= xm) {
1089  hi = mid;
1090  } else {
1091  lo = mid;
1092  }
1093  }
1094 
1095  if (lo == -1) { /* off the bottom */
1096  return (lo);
1097  }
1098  /*
1099  * If there's a discontinuity (many knots at same x-value), choose the
1100  * largest
1101  */
1102  xm = x[lo];
1103  while (lo < n - 1 && x[lo + 1] == xm) lo++;
1104 
1105  return (lo);
1106 }
1107 
1111 static void keep_valid_roots(std::vector<double> &roots, std::vector<double> &newRoots, double x0,
1112  double x1) {
1113  for (unsigned int i = 0; i != newRoots.size(); ++i) {
1114  if (newRoots[i] >= x0 && newRoots[i] < x1) { // keep this root
1115  roots.push_back(newRoots[i]);
1116  }
1117  }
1118 
1119  newRoots.clear();
1120 }
1121 
1125 namespace {
1126 void do_quadratic(double a, double b, double c, std::vector<double> &roots) {
1129  roots.push_back(-c / b);
1130  }
1131  } else {
1132  double const tmp = b * b - 4 * a * c;
1133 
1134  if (tmp >= 0) {
1135  if (b >= 0) {
1136  roots.push_back((-b - sqrt(tmp)) / (2 * a));
1137  } else {
1138  roots.push_back((-b + sqrt(tmp)) / (2 * a));
1139  }
1140  roots.push_back(c / (a * roots[0]));
1141  /*
1142  * sort roots
1143  */
1144  if (roots[0] > roots[1]) {
1145  double const tmp2 = roots[0];
1146  roots[0] = roots[1];
1147  roots[1] = tmp2;
1148  }
1149  }
1150  }
1151 }
1152 
1156 void do_cubic(double a, double b, double c, double d, std::vector<double> &roots) {
1158  do_quadratic(b, c, d, roots);
1159  return;
1160  }
1161  b /= a;
1162  c /= a;
1163  d /= a;
1164 
1165  double const q = (b * b - 3 * c) / 9;
1166  double const r = (2 * b * b * b - 9 * b * c + 27 * d) / 54;
1167  /*
1168  * n.b. note that the test for the number of roots is carried out on the
1169  * same variables as are used in (e.g.) the acos, as it is possible for
1170  * r*r < q*q*q && r > sq*sq*sq due to rounding.
1171  */
1172  double const sq = (q >= 0) ? sqrt(q) : -sqrt(-q);
1173  double const sq3 = sq * sq * sq;
1174  if (::fabs(r) < sq3) { // three real roots
1175  double const theta = ::acos(r / sq3); // sq3 cannot be zero
1176 
1177  roots.push_back(-2 * sq * cos(theta / 3) - b / 3);
1178  roots.push_back(-2 * sq * cos((theta + lsst::geom::TWOPI) / 3) - b / 3);
1179  roots.push_back(-2 * sq * cos((theta - lsst::geom::TWOPI) / 3) - b / 3);
1180  /*
1181  * sort roots
1182  */
1183  if (roots[0] > roots[1]) {
1184  std::swap(roots[0], roots[1]);
1185  }
1186  if (roots[1] > roots[2]) {
1187  std::swap(roots[1], roots[2]);
1188  }
1189  if (roots[0] > roots[1]) {
1190  std::swap(roots[0], roots[1]);
1191  }
1192 
1193  return;
1194  } else if (::fabs(r) == sq3) { /* no more than two real roots */
1195  double const aa = -((r < 0) ? -::pow(-r, 1.0 / 3.0) : ::pow(r, 1.0 / 3.0));
1196 
1197  if (::fabs(aa) < std::numeric_limits<double>::epsilon()) { /* degenerate case; one real root */
1198  roots.push_back(-b / 3);
1199  return;
1200  } else {
1201  roots.push_back(2 * aa - b / 3);
1202  roots.push_back(-aa - b / 3);
1203 
1204  if (roots[0] > roots[1]) {
1205  std::swap(roots[0], roots[1]);
1206  }
1207 
1208  return;
1209  }
1210  } else { /* only one real root */
1211  double tmp = ::sqrt(r * r - (q > 0 ? sq3 * sq3 : -sq3 * sq3));
1212  tmp = r + (r < 0 ? -tmp : tmp);
1213  double const aa = -((tmp < 0) ? -::pow(-tmp, 1.0 / 3.0) : ::pow(tmp, 1.0 / 3.0));
1214  double const bb = (fabs(aa) < std::numeric_limits<double>::epsilon()) ? 0 : q / aa;
1215 
1216  roots.push_back((aa + bb) - b / 3);
1217 #if 0
1218  roots.push_back(-(aa + bb)/2 - b/3); // the real
1219  roots.push_back(::sqrt(3)/2*(aa - bb)); // and imaginary parts of the complex roots
1220 #endif
1221  return;
1222  }
1223 }
1224 } // namespace
1225 
1226 std::vector<double> Spline::roots(double const value, double a, double const b) const {
1227  /*
1228  * Strategy: we know that the interpolant has the form
1229  * val = coef[0][i] +dx*(coef[1][i] + dx*(coef[2][i]/2 + dx*coef[3][i]/6))
1230  * so we can use the usual analytic solution for a cubic. Note that the
1231  * cubic quoted above returns dx, the distance from the previous knot,
1232  * rather than x itself
1233  */
1234  std::vector<double> roots; /* the roots found */
1235  double x0 = a; // lower end of current range
1236  double const x1 = b;
1237  int const nknot = _knots.size();
1238 
1239  int i0 = search_array(x0, &_knots[0], nknot, -1);
1240  int const i1 = search_array(x1, &_knots[0], nknot, i0);
1241  assert(i1 >= i0 && i1 <= nknot - 1);
1242 
1243  std::vector<double> newRoots; // the roots we find in some interval
1244  /*
1245  * Deal with special case that x0 may be off one end or the other of
1246  * the array of knots.
1247  */
1248  if (i0 < 0) { /* off bottom */
1249  i0 = 0;
1250  do_cubic(_coeffs[3][i0] / 6, _coeffs[2][i0] / 2, _coeffs[1][i0], _coeffs[0][i0] - value, newRoots);
1251  //
1252  // Could use
1253  // std::transform(newRoots.begin(), newRoots.end(), newRoots.begin(),
1254  // std::bind(std::plus<double>(), _1, _knots[i0]));
1255  // but let's not
1256  //
1257  for (unsigned int j = 0; j != newRoots.size(); ++j) {
1258  newRoots[j] += _knots[i0];
1259  }
1260  keep_valid_roots(roots, newRoots, x0, _knots[i0]);
1261 
1262  x0 = _knots[i0];
1263  } else if (i0 >= nknot) { /* off top */
1264  i0 = nknot - 1;
1265  assert(i0 >= 0);
1266  do_cubic(_coeffs[3][i0] / 6, _coeffs[2][i0] / 2, _coeffs[1][i0], _coeffs[0][i0] - value, newRoots);
1267 
1268  for (unsigned int j = 0; j != newRoots.size(); ++j) {
1269  newRoots[j] += _knots[i0];
1270  }
1271  keep_valid_roots(roots, newRoots, x0, x1);
1272 
1273  return roots;
1274  }
1275  /*
1276  * OK, now search in main body of spline. Note that i1 may be nknot - 1, and
1277  * in any case the right hand limit of the last segment is at x1, not a knot
1278  */
1279  for (int i = i0; i <= i1; i++) {
1280  do_cubic(_coeffs[3][i] / 6, _coeffs[2][i] / 2, _coeffs[1][i], _coeffs[0][i] - value, newRoots);
1281 
1282  for (unsigned int j = 0; j != newRoots.size(); ++j) {
1283  newRoots[j] += _knots[i];
1284  }
1285  keep_valid_roots(roots, newRoots, ((i == i0) ? x0 : _knots[i]), ((i == i1) ? x1 : _knots[i + 1]));
1286  }
1287 
1288  return roots;
1289 }
1290 } // namespace detail
1291 } // namespace math
1292 } // namespace afw
1293 } // namespace lsst
int const step
double x
table::Key< int > type
Definition: Detector.cc:163
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
double z
Definition: Match.cc:44
int y
Definition: SpanSet.cc:49
table::Key< int > b
table::Key< int > a
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.
Definition: Spline.cc:534
std::vector< double > _knots
Definition: Spline.h:55
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).
Definition: Spline.cc:1226
void derivative(std::vector< double > const &x, std::vector< double > &dydx) const
Find the derivative of a Spline.
Definition: Spline.cc:57
void _allocateSpline(int const nknot)
Allocate the storage a Spline needs.
Definition: Spline.cc:21
void interpolate(std::vector< double > const &x, std::vector< double > &y) const
Interpolate a Spline.
Definition: Spline.cc:29
std::vector< std::vector< double > > _coeffs
Definition: Spline.h:56
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.
Definition: Spline.cc:86
Reports invalid arguments.
Definition: Runtime.h:66
T clear(T... args)
T epsilon(T... args)
T fabs(T... args)
constexpr double TWOPI
Definition: Angle.h:40
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
double cos(Angle const &a)
Definition: Angle.h:103
A base class for image defects.
T pow(T... args)
T push_back(T... args)
T reserve(T... args)
T resize(T... args)
T size(T... args)
T sqrt(T... args)
T swap(T... args)