]> git.donarmstrong.com Git - mothur.git/blob - spline.cpp
added sequence name to error string in fastq.info. Changed np_shannon to npshannon.
[mothur.git] / 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 /***********************************************************/