spline.cpp

   1
   2 #include "spline.h"
   3
   4
   5 extern"C" {
   6         void sgram_(double *sg0, double *sg1, double *sg2,
   7                                 double *sg3, double* knot, int* nk);
   8         void stxwx_(double *xs, double *ys, double *ws, int *n,
   9                                 double *knot, int *nk, double *xwy, double *hs0, double *hs1, double *hs2, double *hs3);
  10         void sslvrg_(double *penalt, double *dofoff, double *xs, double *ys, double *ws, double *ssw, int *n,
  11                                  double *knot, int *nk,
  12                                  double *coef, double *sz, double *lev, double *crit, int *icrit, double *lspar, double *xwy,
  13                                  double *hs0, double *hs1, double *hs2, double *hs3,
  14                                  double *sg0, double *sg1, double *sg2, double *sg3, double *abd,
  15                                  double *p1ip, double *p2ip, int *ld4, int *ldnk, int *ier);
  16 }
  17 /***********************************************************/
  18 // used as reference - http://www.koders.com/fortran/fid8AA63B49CF22F0138E9B3DBDC405696F4A62C1CF.aspx
  19 //  http://www.koders.com/c/fidD995301A8A5549CE0361F4E7FFDFD3CDC4B4E4A3.aspx
  20 /* A Cubic B-spline Smoothing routine.
  21
  22  // sbart.f -- translated by f2c (version 20010821).
  23  // ------- and f2c-clean,v 1.9 2000/01/13
  24  //
  25  // According to the GAMFIT sources, this was derived from code by
  26  // Finbarr O'Sullivan.
  27  //
  28
  29
  30  The algorithm minimises:
  31
  32  (1/n) * sum ws(i)^2 * (ys(i)-sz(i))^2 + lambda* int ( s"(x) )^2 dx
  33
  34  lambda is a function of the spar which is assumed to be between 0 and 1
  35
  36  INPUT
  37  -----
  38  penalt                 A penalty > 1 to be used in the gcv criterion
  39  dofoff                 either `df.offset' for GCV or `df' (to be matched).
  40  n                              number of data points
  41  ys(n)                  vector of length n containing the observations
  42  ws(n)                  vector containing the weights given to each data point
  43  NB: the code alters the values here.
  44  xs(n)                  vector containing the ordinates of the observations
  45  ssw                    `centered weighted sum of y^2'
  46  nk                             number of b-spline coefficients to be estimated
  47  nk <= n+2
  48  knot(nk+4)             vector of knot points defining the cubic b-spline basis.
  49  To obtain full cubic smoothing splines one might
  50  have (provided the xs-values are strictly increasing)
  51  spar                   penalised likelihood smoothing parameter
  52  ispar                  indicating if spar is supplied (ispar=1) or to be estimated
  53  lspar, uspar   lower and upper values for spar search;  0.,1. are good values
  54  tol, eps               used in Golden Search routine
  55  isetup                 setup indicator [initially 0
  56  icrit                  indicator saying which cross validation score is to be computed
  57  0: none ;  1: GCV ;  2: CV ;  3: 'df matching'
  58  ld4                    the leading dimension of abd (ie ld4=4)
  59  ldnk                   the leading dimension of p2ip (not referenced)
  60
  61  OUTPUT
  62  ------
  63  coef(nk)               vector of spline coefficients
  64  sz(n)                  vector of smoothed z-values
  65  lev(n)                 vector of leverages
  66  crit                   either ordinary or generalized CV score
  67  spar                   if ispar != 1
  68  lspar                  == lambda (a function of spar and the design)
  69  iter                   number of iterations needed for spar search (if ispar != 1)
  70  ier                    error indicator
  71  ier = 0 ___  everything fine
  72  ier = 1 ___  spar too small or too big
  73  problem in cholesky decomposition
  74
  75  Working arrays/matrix
  76  xwy                            X'Wy
  77  hs0,hs1,hs2,hs3        the diagonals of the X'WX matrix
  78  sg0,sg1,sg2,sg3        the diagonals of the Gram matrix SIGMA
  79  abd (ld4,nk)           [ X'WX + lambda*SIGMA ] in diagonal form
  80  p1ip(ld4,nk)           inner products between columns of L inverse
  81  p2ip(ldnk,nk)          all inner products between columns of L inverse
  82  where  L'L = [X'WX + lambda*SIGMA]  NOT REFERENCED
  83  */
  84
  85 int Spline::sbart
  86 (double *penalt, double *dofoff,
  87  double *xs, double *ys, double *ws, double *ssw,
  88  int *n, double *knot, int *nk, double *coef,
  89  double *sz, double *lev, double *crit, int *icrit,
  90  double *spar, int *ispar, int *iter, double *lspar,
  91  double *uspar, double *tol, double *eps, int *isetup,
  92  int *ld4, int *ldnk, int *ier)
  93 {
  94         try{
  95         /* A Cubic B-spline Smoothing routine.
  96
  97          The algorithm minimises:
  98
  99          (1/n) * sum ws(i)^2 * (ys(i)-sz(i))^2 + lambda* int ( s(x) )^2 dx
 100
 101          lambda is a function of the spar which is assumed to be between 0 and 1
 102
 103          INPUT
 104          -----
 105          penalt A penalty > 1 to be used in the gcv criterion
 106          dofoff either `df.offset' for GCV or `df' (to be matched).
 107          n              number of data points
 108          ys(n)  vector of length n containing the observations
 109          ws(n)  vector containing the weights given to each data point
 110          NB: the code alters the values here.
 111          xs(n)  vector containing the ordinates of the observations
 112          ssw          `centered weighted sum of y^2
 113          nk             number of b-spline coefficients to be estimated
 114          nk <= n+2
 115          knot(nk+4)     vector of knot points defining the cubic b-spline basis.
 116          To obtain full cubic smoothing splines one might
 117          have (provided the xs-values are strictly increasing)
 118          spar           penalised likelihood smoothing parameter
 119          ispar  indicating if spar is supplied (ispar=1) or to be estimated
 120          lspar, uspar lower and upper values for spar search;  0.,1. are good values
 121          tol, eps       used in Golden Search routine
 122          isetup setup indicator [initially 0
 123          icrit  indicator saying which cross validation score is to be computed
 124          0: none ;  1: GCV ;  2: CV ;  3: 'df matching'
 125          ld4            the leading dimension of abd (ie ld4=4)
 126          ldnk           the leading dimension of p2ip (not referenced)
 127
 128          OUTPUT
 129          ------
 130          coef(nk)       vector of spline coefficients
 131          sz(n)  vector of smoothed z-values
 132          lev(n) vector of leverages
 133          crit           either ordinary or generalized CV score
 134          spar         if ispar != 1
 135          lspar         == lambda (a function of spar and the design)
 136          iter           number of iterations needed for spar search (if ispar != 1)
 137          ier            error indicator
 138          ier = 0 ___  everything fine
 139          ier = 1 ___  spar too small or too big
 140          problem in cholesky decomposition
 141
 142          Working arrays/matrix
 143          xwy                    XWy
 144          hs0,hs1,hs2,hs3        the diagonals of the XWX matrix
 145          sg0,sg1,sg2,sg3        the diagonals of the Gram matrix SIGMA
 146          abd (ld4,nk)           [ XWX + lambda*SIGMA ] in diagonal form
 147          p1ip(ld4,nk)           inner products between columns of L inverse
 148          p2ip(ldnk,nk)  all inner products between columns of L inverse
 149          where  LL = [XWX + lambda*SIGMA]  NOT REFERENCED
 150          */
 151
 152 #define CRIT(FX) (*icrit == 3 ? FX - 3. : FX)
 153                 /* cancellation in (3 + eps) - 3, but still...informative */
 154
 155 #define BIG_f (1e100)
 156
 157                 /* c_Gold is the squared inverse of the golden ratio */
 158                 static const double c_Gold = 0.381966011250105151795413165634;
 159                 /* == (3. - sqrt(5.)) / 2. */
 160
 161                 /* Local variables */
 162                 static double ratio;/* must be static (not needed in R) */
 163
 164                 double a, b, d, e, p, q, r, u, v, w, x;
 165                 double ax, fu, fv, fw, fx, bx, xm;
 166                 double t1, t2, tol1, tol2;
 167                 double* xwy = new double[*nk];
 168                 double* hs0 = new double[*nk];
 169                 double* hs1 = new double[*nk];
 170                 double* hs2 = new double[*nk];
 171                 double* hs3 = new double[*nk];
 172                 double* sg0 = new double[*nk];
 173                 double* sg1 = new double[*nk];
 174                 double* sg2 = new double[*nk];
 175                 double* sg3 = new double[*nk];
 176                 double* abd = new double[*nk*(*ld4)];
 177                 double* p1ip = new double[*nk*(*ld4)];
 178                 double* p2ip = new double[*nk];
 179                 int i, maxit;
 180
 181                 /* unnecessary initializations to keep  -Wall happy */
 182                 d = 0.; fu = 0.; u = 0.;
 183                 ratio = 1.;
 184
 185                 /*  Compute SIGMA, X' W X, X' W z, trace ratio, s0, s1.
 186
 187                  SIGMA  -> sg0,sg1,sg2,sg3
 188                  X' W X -> hs0,hs1,hs2,hs3
 189                  X' W Z -> xwy
 190                  */
 191
 192                 /* trevor fixed this 4/19/88
 193                  * Note: sbart, i.e. stxwx() and sslvrg() {mostly, not always!}, use
 194                  *       the square of the weights; the following rectifies that */
 195                 for (i = 0; i < *n; ++i)
 196                         if (ws[i] > 0.)
 197                                 ws[i] = sqrt(ws[i]);
 198
 199                 if (*isetup == 0) {
 200                         /* SIGMA[i,j] := Int  B''(i,t) B''(j,t) dt  {B(k,.) = k-th B-spline} */
 201                         sgram_(sg0, sg1, sg2, sg3, knot, nk);
 202                         stxwx_(xs, ys, ws, n,
 203                                                         knot, nk,
 204                                                         xwy,
 205                                                         hs0, hs1, hs2, hs3);
 206                         /* Compute ratio :=  tr(X' W X) / tr(SIGMA) */
 207                         t1 = t2 = 0.;
 208                         for (i = 3 - 1; i < (*nk - 3); ++i) {
 209                                 t1 += hs0[i];
 210                                 t2 += sg0[i];
 211                         }
 212                         ratio = t1 / t2;
 213                         *isetup = 1;
 214                 }
 215                 /*     Compute estimate */
 216
 217                 if (*ispar == 1) { /* Value of spar supplied */
 218                         *lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0));
 219                         sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n,
 220                                                          knot, nk,
 221                                                          coef, sz, lev, crit, icrit, lspar, xwy,
 222                                                          hs0, hs1, hs2, hs3,
 223                                                          sg0, sg1, sg2, sg3, abd,
 224                                                          p1ip, p2ip, ld4, ldnk, ier);
 225                         /* got through check 2 */
 226                         return 0;
 227                 }
 228
 229                 /* ELSE ---- spar not supplied --> compute it ! ---------------------------
 230
 231                  Use Forsythe Malcom and Moler routine to MINIMIZE criterion
 232                  f denotes the value of the criterion
 233
 234                  an approximation       x  to the point where   f  attains a minimum  on
 235                  the interval  (ax,bx)  is determined.
 236                  */
 237                 ax = *lspar;
 238                 bx = *uspar;
 239
 240                 /* INPUT
 241
 242                  ax      left endpoint of initial interval
 243                  bx      right endpoint of initial interval
 244                  f       function subprogram which evaluates  f(x)  for any  x
 245          in the interval  (ax,bx)
 246                  tol     desired length of the interval of uncertainty of the final
 247          result ( >= 0 )
 248
 249                  OUTPUT
 250
 251                  fmin    abcissa approximating the point where  f  attains a minimum
 252                  */
 253
 254                 /*
 255                  The method used is a combination of  golden  section  search  and
 256                  successive parabolic interpolation.    convergence is never much slower
 257                  than    that  for  a  fibonacci search.  if  f  has a continuous second
 258                  derivative which is positive at the minimum (which is not  at  ax  or
 259                  bx),    then  convergence  is  superlinear, and usually of the order of
 260                  about  1.324....
 261                  the function  f  is never evaluated at two points closer together
 262                  than    eps*abs(fmin) + (tol/3), where eps is  approximately the square
 263                  root    of  the  relative  machine  precision.   if   f   is a unimodal
 264                  function and the computed values of     f   are  always  unimodal  when
 265                  separated by at least  eps*abs(x) + (tol/3), then  fmin  approximates
 266                  the abcissa of the global minimum of    f  on the interval  ax,bx  with
 267                  an error less than  3*eps*abs(fmin) + tol.  if   f     is not unimodal,
 268                  then fmin may approximate a local, but perhaps non-global, minimum to
 269                  the same accuracy.
 270                  this function subprogram is a slightly modified        version  of  the
 271                  algol  60 procedure    localmin  given in richard brent, algorithms for
 272                  minimization without derivatives, prentice - hall, inc. (1973).
 273
 274                  Double  a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w
 275                  Double  fu,fv,fw,fx,x
 276                  */
 277
 278                 /*  eps is approximately the square root of the relative machine
 279                  precision.
 280
 281                  -       eps = 1e0
 282                  - 10    eps = eps/2e0
 283                  -       tol1 = 1e0 + eps
 284                  -       if (tol1 > 1e0) go to 10
 285                  -       eps = sqrt(eps)
 286                  R Version <= 1.3.x had
 287                  eps = .000244     ( = sqrt(5.954 e-8) )
 288                  -- now eps is passed as argument
 289                  */
 290
 291                 /* initialization */
 292
 293                 maxit = *iter;
 294                 *iter = 0;
 295                 a = ax;
 296                 b = bx;
 297                 v = a + c_Gold * (b - a);
 298                 w = v;
 299                 x = v;
 300                 e = 0.;
 301                 *spar = x;
 302                 *lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0));
 303                 sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n,
 304                                                  knot, nk,
 305                                                  coef, sz, lev, crit, icrit, lspar, xwy,
 306                                 hs0, hs1, hs2, hs3,
 307                                 sg0, sg1, sg2, sg3, abd,
 308                                 p1ip, p2ip, ld4, ldnk, ier);
 309                 fx = *crit;
 310                 fv = fx;
 311                 fw = fx;
 312
 313                 /* main loop
 314                  --------- */
 315                 while(*ier == 0) { /* L20: */
 316
 317                         if (m->control_pressed) { return 0; }
 318
 319                         xm = (a + b) * .5;
 320                         tol1 = *eps * fabs(x) + *tol / 3.;
 321                         tol2 = tol1 * 2.;
 322                         ++(*iter);
 323
 324
 325                         /* Check the (somewhat peculiar) stopping criterion: note that
 326                          the RHS is negative as long as the interval [a,b] is not small:*/
 327                         if (fabs(x - xm) <= tol2 - (b - a) * .5 || *iter > maxit)
 328                                 goto L_End;
 329
 330
 331                         /* is golden-section necessary */
 332
 333                         if (fabs(e) <= tol1 ||
 334                                 /*  if had Inf then go to golden-section */
 335                                 fx >= BIG_f || fv >= BIG_f || fw >= BIG_f) goto L_GoldenSect;
 336
 337                         /* Fit Parabola */
 338
 339                         r = (x - w) * (fx - fv);
 340                         q = (x - v) * (fx - fw);
 341                         p = (x - v) * q - (x - w) * r;
 342                         q = (q - r) * 2.;
 343                         if (q > 0.)
 344                                 p = -p;
 345                         q = fabs(q);
 346                         r = e;
 347                         e = d;
 348
 349                         /* is parabola acceptable?  Otherwise do golden-section */
 350
 351                         if (fabs(p) >= fabs(.5 * q * r) ||
 352                                 q == 0.)
 353                         /* above line added by BDR;
 354                          * [the abs(.) >= abs() = 0 should have branched..]
 355                          * in FTN: COMMON above ensures q is NOT a register variable */
 356
 357                                 goto L_GoldenSect;
 358
 359                         if (p <= q * (a - x) ||
 360                                 p >= q * (b - x))                       goto L_GoldenSect;
 361
 362
 363
 364                         /* Parabolic Interpolation step */
 365                         d = p / q;
 366                         u = x + d;
 367
 368                         /* f must not be evaluated too close to ax or bx */
 369                         if ((u - a < tol2 || b - u < tol2)) {
 370                                 d = abs(tol1) * sgn(xm - x);
 371                         }
 372
 373                         goto L50;
 374                         /*------*/
 375
 376                 L_GoldenSect: /* a golden-section step */
 377
 378                         if (x >= xm)    e = a - x;
 379                         else/* x < xm*/ e = b - x;
 380                         d = c_Gold * e;
 381
 382
 383                 L50:
 384                         u = x + ((fabs(d) >= tol1) ? d : (abs(tol1)*sgn(d)));
 385                         /*  tol1 check : f must not be evaluated too close to x */
 386
 387                         *spar = u;
 388                         *lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0));
 389                         sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n,
 390                                                          knot, nk,
 391                                                          coef, sz, lev, crit, icrit, lspar, xwy,
 392                                         hs0, hs1, hs2, hs3,
 393                                         sg0, sg1, sg2, sg3, abd,
 394                                         p1ip, p2ip, ld4, ldnk, ier);
 395                         fu = *crit;
 396
 397                         if(isnan(fu)) {
 398                                 fu = 2. * BIG_f;
 399                         }
 400
 401                         /*  update  a, b, v, w, and x */
 402
 403                         if (fu <= fx) {
 404                                 if (u >= x) a = x; else b = x;
 405
 406                                 v = w; fv = fw;
 407                                 w = x; fw = fx;
 408                                 x = u; fx = fu;
 409                         }
 410                         else {
 411                                 if (u < x)  a = u; else b = u;
 412
 413                                 if (fu <= fw || w == x) {                       /* L70: */
 414                                         v = w; fv = fw;
 415                                         w = u; fw = fu;
 416                                 } else if (fu <= fv || v == x || v == w) {      /* L80: */
 417                                         v = u; fv = fu;
 418                                 }
 419                         }
 420                 }/* end main loop -- goto L20; */
 421
 422         L_End:
 423
 424                 *spar = x;
 425                 *crit = fx;
 426
 427                 //free memory
 428                 delete [] xwy;
 429                 delete [] hs0;
 430                 delete [] hs1;
 431                 delete [] hs2;
 432                 delete [] hs3;
 433                 delete [] sg0;
 434                 delete [] sg1;
 435                 delete [] sg2;
 436                 delete [] sg3;
 437                 delete [] abd;
 438                 delete [] p1ip;
 439                 delete [] p2ip;
 440
 441                 return 0;
 442                 /* sbart */
 443
 444         }catch(exception& e) {
 445                 m->errorOut(e, "Spline", "sbart");
 446                 exit(1);
 447         }
 448 }
 449
 450 /***********************************************************/