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