--- /dev/null
+
+#include "spline.h"
+
+
+extern"C" {
+ void sgram_(double *sg0, double *sg1, double *sg2,
+ double *sg3, double* knot, int* nk);
+ void stxwx_(double *xs, double *ys, double *ws, int *n,
+ double *knot, int *nk, double *xwy, double *hs0, double *hs1, double *hs2, double *hs3);
+ void sslvrg_(double *penalt, double *dofoff, double *xs, double *ys, double *ws, double *ssw, int *n,
+ double *knot, int *nk,
+ double *coef, double *sz, double *lev, double *crit, int *icrit, double *lspar, double *xwy,
+ double *hs0, double *hs1, double *hs2, double *hs3,
+ double *sg0, double *sg1, double *sg2, double *sg3, double *abd,
+ double *p1ip, double *p2ip, int *ld4, int *ldnk, int *ier);
+}
+/***********************************************************/
+// used as reference - http://www.koders.com/fortran/fid8AA63B49CF22F0138E9B3DBDC405696F4A62C1CF.aspx
+// http://www.koders.com/c/fidD995301A8A5549CE0361F4E7FFDFD3CDC4B4E4A3.aspx
+/* A Cubic B-spline Smoothing routine.
+
+ // sbart.f -- translated by f2c (version 20010821).
+ // ------- and f2c-clean,v 1.9 2000/01/13
+ //
+ // According to the GAMFIT sources, this was derived from code by
+ // Finbarr O'Sullivan.
+ //
+
+
+ The algorithm minimises:
+
+ (1/n) * sum ws(i)^2 * (ys(i)-sz(i))^2 + lambda* int ( s"(x) )^2 dx
+
+ lambda is a function of the spar which is assumed to be between 0 and 1
+
+ INPUT
+ -----
+ penalt A penalty > 1 to be used in the gcv criterion
+ dofoff either `df.offset' for GCV or `df' (to be matched).
+ n number of data points
+ ys(n) vector of length n containing the observations
+ ws(n) vector containing the weights given to each data point
+ NB: the code alters the values here.
+ xs(n) vector containing the ordinates of the observations
+ ssw `centered weighted sum of y^2'
+ nk number of b-spline coefficients to be estimated
+ nk <= n+2
+ knot(nk+4) vector of knot points defining the cubic b-spline basis.
+ To obtain full cubic smoothing splines one might
+ have (provided the xs-values are strictly increasing)
+ spar penalised likelihood smoothing parameter
+ ispar indicating if spar is supplied (ispar=1) or to be estimated
+ lspar, uspar lower and upper values for spar search; 0.,1. are good values
+ tol, eps used in Golden Search routine
+ isetup setup indicator [initially 0
+ icrit indicator saying which cross validation score is to be computed
+ 0: none ; 1: GCV ; 2: CV ; 3: 'df matching'
+ ld4 the leading dimension of abd (ie ld4=4)
+ ldnk the leading dimension of p2ip (not referenced)
+
+ OUTPUT
+ ------
+ coef(nk) vector of spline coefficients
+ sz(n) vector of smoothed z-values
+ lev(n) vector of leverages
+ crit either ordinary or generalized CV score
+ spar if ispar != 1
+ lspar == lambda (a function of spar and the design)
+ iter number of iterations needed for spar search (if ispar != 1)
+ ier error indicator
+ ier = 0 ___ everything fine
+ ier = 1 ___ spar too small or too big
+ problem in cholesky decomposition
+
+ Working arrays/matrix
+ xwy X'Wy
+ hs0,hs1,hs2,hs3 the diagonals of the X'WX matrix
+ sg0,sg1,sg2,sg3 the diagonals of the Gram matrix SIGMA
+ abd (ld4,nk) [ X'WX + lambda*SIGMA ] in diagonal form
+ p1ip(ld4,nk) inner products between columns of L inverse
+ p2ip(ldnk,nk) all inner products between columns of L inverse
+ where L'L = [X'WX + lambda*SIGMA] NOT REFERENCED
+ */
+
+int Spline::sbart
+(double *penalt, double *dofoff,
+ double *xs, double *ys, double *ws, double *ssw,
+ int *n, double *knot, int *nk, double *coef,
+ double *sz, double *lev, double *crit, int *icrit,
+ double *spar, int *ispar, int *iter, double *lspar,
+ double *uspar, double *tol, double *eps, int *isetup,
+ int *ld4, int *ldnk, int *ier)
+{
+ try{
+ /* A Cubic B-spline Smoothing routine.
+
+ The algorithm minimises:
+
+ (1/n) * sum ws(i)^2 * (ys(i)-sz(i))^2 + lambda* int ( s(x) )^2 dx
+
+ lambda is a function of the spar which is assumed to be between 0 and 1
+
+ INPUT
+ -----
+ penalt A penalty > 1 to be used in the gcv criterion
+ dofoff either `df.offset' for GCV or `df' (to be matched).
+ n number of data points
+ ys(n) vector of length n containing the observations
+ ws(n) vector containing the weights given to each data point
+ NB: the code alters the values here.
+ xs(n) vector containing the ordinates of the observations
+ ssw `centered weighted sum of y^2
+ nk number of b-spline coefficients to be estimated
+ nk <= n+2
+ knot(nk+4) vector of knot points defining the cubic b-spline basis.
+ To obtain full cubic smoothing splines one might
+ have (provided the xs-values are strictly increasing)
+ spar penalised likelihood smoothing parameter
+ ispar indicating if spar is supplied (ispar=1) or to be estimated
+ lspar, uspar lower and upper values for spar search; 0.,1. are good values
+ tol, eps used in Golden Search routine
+ isetup setup indicator [initially 0
+ icrit indicator saying which cross validation score is to be computed
+ 0: none ; 1: GCV ; 2: CV ; 3: 'df matching'
+ ld4 the leading dimension of abd (ie ld4=4)
+ ldnk the leading dimension of p2ip (not referenced)
+
+ OUTPUT
+ ------
+ coef(nk) vector of spline coefficients
+ sz(n) vector of smoothed z-values
+ lev(n) vector of leverages
+ crit either ordinary or generalized CV score
+ spar if ispar != 1
+ lspar == lambda (a function of spar and the design)
+ iter number of iterations needed for spar search (if ispar != 1)
+ ier error indicator
+ ier = 0 ___ everything fine
+ ier = 1 ___ spar too small or too big
+ problem in cholesky decomposition
+
+ Working arrays/matrix
+ xwy XWy
+ hs0,hs1,hs2,hs3 the diagonals of the XWX matrix
+ sg0,sg1,sg2,sg3 the diagonals of the Gram matrix SIGMA
+ abd (ld4,nk) [ XWX + lambda*SIGMA ] in diagonal form
+ p1ip(ld4,nk) inner products between columns of L inverse
+ p2ip(ldnk,nk) all inner products between columns of L inverse
+ where LL = [XWX + lambda*SIGMA] NOT REFERENCED
+ */
+
+#define CRIT(FX) (*icrit == 3 ? FX - 3. : FX)
+ /* cancellation in (3 + eps) - 3, but still...informative */
+
+#define BIG_f (1e100)
+
+ /* c_Gold is the squared inverse of the golden ratio */
+ static const double c_Gold = 0.381966011250105151795413165634;
+ /* == (3. - sqrt(5.)) / 2. */
+
+ /* Local variables */
+ static double ratio;/* must be static (not needed in R) */
+
+ double a, b, d, e, p, q, r, u, v, w, x;
+ double ax, fu, fv, fw, fx, bx, xm;
+ double t1, t2, tol1, tol2;
+ double* xwy = new double[*nk];
+ double* hs0 = new double[*nk];
+ double* hs1 = new double[*nk];
+ double* hs2 = new double[*nk];
+ double* hs3 = new double[*nk];
+ double* sg0 = new double[*nk];
+ double* sg1 = new double[*nk];
+ double* sg2 = new double[*nk];
+ double* sg3 = new double[*nk];
+ double* abd = new double[*nk*(*ld4)];
+ double* p1ip = new double[*nk*(*ld4)];
+ double* p2ip = new double[*nk];
+ int i, maxit;
+
+ /* unnecessary initializations to keep -Wall happy */
+ d = 0.; fu = 0.; u = 0.;
+ ratio = 1.;
+
+ /* Compute SIGMA, X' W X, X' W z, trace ratio, s0, s1.
+
+ SIGMA -> sg0,sg1,sg2,sg3
+ X' W X -> hs0,hs1,hs2,hs3
+ X' W Z -> xwy
+ */
+
+ /* trevor fixed this 4/19/88
+ * Note: sbart, i.e. stxwx() and sslvrg() {mostly, not always!}, use
+ * the square of the weights; the following rectifies that */
+ for (i = 0; i < *n; ++i)
+ if (ws[i] > 0.)
+ ws[i] = sqrt(ws[i]);
+
+ if (*isetup == 0) {
+ /* SIGMA[i,j] := Int B''(i,t) B''(j,t) dt {B(k,.) = k-th B-spline} */
+ sgram_(sg0, sg1, sg2, sg3, knot, nk);
+ stxwx_(xs, ys, ws, n,
+ knot, nk,
+ xwy,
+ hs0, hs1, hs2, hs3);
+ /* Compute ratio := tr(X' W X) / tr(SIGMA) */
+ t1 = t2 = 0.;
+ for (i = 3 - 1; i < (*nk - 3); ++i) {
+ t1 += hs0[i];
+ t2 += sg0[i];
+ }
+ ratio = t1 / t2;
+ *isetup = 1;
+ }
+ /* Compute estimate */
+
+ if (*ispar == 1) { /* Value of spar supplied */
+ *lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0));
+ sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n,
+ knot, nk,
+ coef, sz, lev, crit, icrit, lspar, xwy,
+ hs0, hs1, hs2, hs3,
+ sg0, sg1, sg2, sg3, abd,
+ p1ip, p2ip, ld4, ldnk, ier);
+ /* got through check 2 */
+ return 0;
+ }
+
+ /* ELSE ---- spar not supplied --> compute it ! ---------------------------
+
+ Use Forsythe Malcom and Moler routine to MINIMIZE criterion
+ f denotes the value of the criterion
+
+ an approximation x to the point where f attains a minimum on
+ the interval (ax,bx) is determined.
+ */
+ ax = *lspar;
+ bx = *uspar;
+
+ /* INPUT
+
+ ax left endpoint of initial interval
+ bx right endpoint of initial interval
+ f function subprogram which evaluates f(x) for any x
+ in the interval (ax,bx)
+ tol desired length of the interval of uncertainty of the final
+ result ( >= 0 )
+
+ OUTPUT
+
+ fmin abcissa approximating the point where f attains a minimum
+ */
+
+ /*
+ The method used is a combination of golden section search and
+ successive parabolic interpolation. convergence is never much slower
+ than that for a fibonacci search. if f has a continuous second
+ derivative which is positive at the minimum (which is not at ax or
+ bx), then convergence is superlinear, and usually of the order of
+ about 1.324....
+ the function f is never evaluated at two points closer together
+ than eps*abs(fmin) + (tol/3), where eps is approximately the square
+ root of the relative machine precision. if f is a unimodal
+ function and the computed values of f are always unimodal when
+ separated by at least eps*abs(x) + (tol/3), then fmin approximates
+ the abcissa of the global minimum of f on the interval ax,bx with
+ an error less than 3*eps*abs(fmin) + tol. if f is not unimodal,
+ then fmin may approximate a local, but perhaps non-global, minimum to
+ the same accuracy.
+ this function subprogram is a slightly modified version of the
+ algol 60 procedure localmin given in richard brent, algorithms for
+ minimization without derivatives, prentice - hall, inc. (1973).
+
+ Double a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w
+ Double fu,fv,fw,fx,x
+ */
+
+ /* eps is approximately the square root of the relative machine
+ precision.
+
+ - eps = 1e0
+ - 10 eps = eps/2e0
+ - tol1 = 1e0 + eps
+ - if (tol1 > 1e0) go to 10
+ - eps = sqrt(eps)
+ R Version <= 1.3.x had
+ eps = .000244 ( = sqrt(5.954 e-8) )
+ -- now eps is passed as argument
+ */
+
+ /* initialization */
+
+ maxit = *iter;
+ *iter = 0;
+ a = ax;
+ b = bx;
+ v = a + c_Gold * (b - a);
+ w = v;
+ x = v;
+ e = 0.;
+ *spar = x;
+ *lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0));
+ sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n,
+ knot, nk,
+ coef, sz, lev, crit, icrit, lspar, xwy,
+ hs0, hs1, hs2, hs3,
+ sg0, sg1, sg2, sg3, abd,
+ p1ip, p2ip, ld4, ldnk, ier);
+ fx = *crit;
+ fv = fx;
+ fw = fx;
+
+ /* main loop
+ --------- */
+ while(*ier == 0) { /* L20: */
+
+ if (m->control_pressed) { return 0; }
+
+ xm = (a + b) * .5;
+ tol1 = *eps * fabs(x) + *tol / 3.;
+ tol2 = tol1 * 2.;
+ ++(*iter);
+
+
+ /* Check the (somewhat peculiar) stopping criterion: note that
+ the RHS is negative as long as the interval [a,b] is not small:*/
+ if (fabs(x - xm) <= tol2 - (b - a) * .5 || *iter > maxit)
+ goto L_End;
+
+
+ /* is golden-section necessary */
+
+ if (fabs(e) <= tol1 ||
+ /* if had Inf then go to golden-section */
+ fx >= BIG_f || fv >= BIG_f || fw >= BIG_f) goto L_GoldenSect;
+
+ /* Fit Parabola */
+
+ r = (x - w) * (fx - fv);
+ q = (x - v) * (fx - fw);
+ p = (x - v) * q - (x - w) * r;
+ q = (q - r) * 2.;
+ if (q > 0.)
+ p = -p;
+ q = fabs(q);
+ r = e;
+ e = d;
+
+ /* is parabola acceptable? Otherwise do golden-section */
+
+ if (fabs(p) >= fabs(.5 * q * r) ||
+ q == 0.)
+ /* above line added by BDR;
+ * [the abs(.) >= abs() = 0 should have branched..]
+ * in FTN: COMMON above ensures q is NOT a register variable */
+
+ goto L_GoldenSect;
+
+ if (p <= q * (a - x) ||
+ p >= q * (b - x)) goto L_GoldenSect;
+
+
+
+ /* Parabolic Interpolation step */
+ d = p / q;
+ u = x + d;
+
+ /* f must not be evaluated too close to ax or bx */
+ if ((u - a < tol2 || b - u < tol2)) {
+ d = abs(tol1) * sgn(xm - x);
+ }
+
+ goto L50;
+ /*------*/
+
+ L_GoldenSect: /* a golden-section step */
+
+ if (x >= xm) e = a - x;
+ else/* x < xm*/ e = b - x;
+ d = c_Gold * e;
+
+
+ L50:
+ u = x + ((fabs(d) >= tol1) ? d : (abs(tol1)*sgn(d)));
+ /* tol1 check : f must not be evaluated too close to x */
+
+ *spar = u;
+ *lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0));
+ sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n,
+ knot, nk,
+ coef, sz, lev, crit, icrit, lspar, xwy,
+ hs0, hs1, hs2, hs3,
+ sg0, sg1, sg2, sg3, abd,
+ p1ip, p2ip, ld4, ldnk, ier);
+ fu = *crit;
+
+ if(isnan(fu)) {
+ fu = 2. * BIG_f;
+ }
+
+ /* update a, b, v, w, and x */
+
+ if (fu <= fx) {
+ if (u >= x) a = x; else b = x;
+
+ v = w; fv = fw;
+ w = x; fw = fx;
+ x = u; fx = fu;
+ }
+ else {
+ if (u < x) a = u; else b = u;
+
+ if (fu <= fw || w == x) { /* L70: */
+ v = w; fv = fw;
+ w = u; fw = fu;
+ } else if (fu <= fv || v == x || v == w) { /* L80: */
+ v = u; fv = fu;
+ }
+ }
+ }/* end main loop -- goto L20; */
+
+ L_End:
+
+ *spar = x;
+ *crit = fx;
+
+ //free memory
+ delete [] xwy;
+ delete [] hs0;
+ delete [] hs1;
+ delete [] hs2;
+ delete [] hs3;
+ delete [] sg0;
+ delete [] sg1;
+ delete [] sg2;
+ delete [] sg3;
+ delete [] abd;
+ delete [] p1ip;
+ delete [] p2ip;
+
+ return 0;
+ /* sbart */
+
+ }catch(exception& e) {
+ m->errorOut(e, "Spline", "sbart");
+ exit(1);
+ }
+}
+
+/***********************************************************/