Cubic spline data smoother.
Algorithm 642 collected algorithms from ACM. Algorithm appeared in Acm-Trans. Math. Software, vol.12, no. 2, Jun., 1986, p. 150.
Translated from fortran by a combination of f2c and RHL.
536 {
537 float var = 1;
538 int const n =
x.size();
539 double const ratio = 2.0;
540 double const tau = 1.618033989;
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
546
548
550 double *c[3];
554
556
557 double *r[3];
558 r[0] = &scratch[0] + 1;
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
568
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) {
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 {
582
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;
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;
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
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
663
664 spcof1(&
x[0], avh, &f[0], &sdf[0], n, p, q,
y, c, u, v);
665
666 stat[2] /= avdf * avdf;
667 stat[3] /= avdf * avdf;
668 stat[4] /= avdf * avdf;
669
670
671
672 if (errs != nullptr) {
673 sperr1(&
x[0], avh, &sdf[0], n, r, p, avar, errs);
674 }
675
676
677
678 if (chisq != nullptr) {
679 *chisq = n * stat[4];
680 }
681}
void _allocateSpline(int const nknot)
Allocate the storage a Spline needs.
std::vector< std::vector< double > > _coeffs