LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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
14namespace lsst {
15namespace afw {
16namespace math {
17namespace detail {
18
19static int search_array(double z, double const *x, int n, int i);
20
21void 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
112void 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 */
521static 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
524static void sperr1(const double x[], double avh, const double dy[], int n, double *r[3], double p, double var,
526
527static 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
531static 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 */
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 != nullptr) {
673 sperr1(&x[0], avh, &sdf[0], n, r, p, avar, errs);
674 }
675 /*
676 * clean up
677 */
678 if (chisq != nullptr) {
679 *chisq = n * stat[4];
680 }
681}
682
692static 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
785static 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
900static 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
928static void sperr1(const double x[], double avh, const double dy[], int n, double *r[3], double p, double var,
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
961void 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
1038static 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
1111static 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
1125namespace {
1126void do_quadratic(double a, double b, double c, std::vector<double> &roots) {
1127 if (::fabs(a) < std::numeric_limits<double>::epsilon()) {
1128 if (::fabs(b) >= std::numeric_limits<double>::epsilon()) {
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
1156void do_cubic(double a, double b, double c, double d, std::vector<double> &roots) {
1157 if (::fabs(a) < std::numeric_limits<double>::epsilon()) {
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
1226std::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
#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:48
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=nullptr, std::vector< double > *errs=nullptr)
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)
double constexpr TWOPI
Definition Angle.h:41
T push_back(T... args)
T reserve(T... args)
T resize(T... args)
T size(T... args)
T swap(T... args)