|
0
|
1 /*
|
|
|
2 * Copyright 1996-2006 Catherine Loader.
|
|
|
3 */
|
|
|
4
|
|
|
5 #include "mex.h"
|
|
|
6 /*
|
|
|
7 * Copyright 1996-2006 Catherine Loader.
|
|
|
8 */
|
|
|
9 #include "lfev.h"
|
|
|
10
|
|
|
11 extern void fitoptions();
|
|
|
12
|
|
|
13 static double hmin, gmin, sig2, pen, vr, tb;
|
|
|
14 static lfit *blf;
|
|
|
15 static design *bdes;
|
|
|
16
|
|
|
17 int procvbind(des,lf,v)
|
|
|
18 design *des;
|
|
|
19 lfit *lf;
|
|
|
20 int v;
|
|
|
21 { double s0, s1, bi;
|
|
|
22 int i, ii, k;
|
|
|
23 k = procv_var(des,lf,v);
|
|
|
24 wdiag(&lf->lfd, &lf->sp, des,des->wd,&lf->dv,0,1,0);
|
|
|
25 s0 = s1 = 0.0;
|
|
|
26 for (i=0; i<des->n; i++)
|
|
|
27 { ii = des->ind[i];
|
|
|
28 s0+= prwt(&lf->lfd,ii)*des->wd[i]*des->wd[i];
|
|
|
29 bi = prwt(&lf->lfd,ii)*fabs(des->wd[i]*ipower(dist(des,ii),deg(&lf->sp)+1));
|
|
|
30 s1+= bi*bi;
|
|
|
31 }
|
|
|
32 vr += s0;
|
|
|
33 tb += s1;
|
|
|
34 return(k);
|
|
|
35 }
|
|
|
36
|
|
|
37 double bcri(h,c,cri)
|
|
|
38 double h;
|
|
|
39 int c, cri;
|
|
|
40 { double num, den, res[10];
|
|
|
41 int (*pv)();
|
|
|
42 if (c==DALP)
|
|
|
43 blf->sp.nn = h;
|
|
|
44 else
|
|
|
45 blf->sp.fixh = h;
|
|
|
46 if ((cri&63)==BIND)
|
|
|
47 { pv = procvbind;
|
|
|
48 vr = tb = 0.0;
|
|
|
49 }
|
|
|
50 else pv = procvstd;
|
|
|
51 if (cri<64) startlf(bdes,blf,pv,0);
|
|
|
52 switch(cri&63)
|
|
|
53 { case BGCV:
|
|
|
54 ressumm(blf,bdes,res);
|
|
|
55 num = -2*blf->lfd.n*res[0];
|
|
|
56 den = blf->lfd.n-res[1];
|
|
|
57 return(num/(den*den));
|
|
|
58 case BCP:
|
|
|
59 ressumm(blf,bdes,res);
|
|
|
60 return(-2*res[0]/sig2-blf->lfd.n+pen*res[1]);
|
|
|
61 case BIND:
|
|
|
62 return(vr+pen*pen*tb);
|
|
|
63 }
|
|
|
64 LERR(("bcri: unknown criterion"));
|
|
|
65 return(0.0);
|
|
|
66 }
|
|
|
67
|
|
|
68 void bsel2(h0,g0,ifact,c,cri)
|
|
|
69 double h0, g0, ifact;
|
|
|
70 int c, cri;
|
|
|
71 { int done, inc;
|
|
|
72 double h1, g1;
|
|
|
73 h1 = h0; g1 = g0;
|
|
|
74 done = inc = 0;
|
|
|
75 while (!done)
|
|
|
76 { h1 *= 1+ifact;
|
|
|
77 g0 = g1;
|
|
|
78 g1 = bcri(h1,c,cri);
|
|
|
79 if (g1<gmin) { hmin = h1; gmin = g1; }
|
|
|
80 if (g1>g0) inc++; else inc = 0;
|
|
|
81 switch(cri)
|
|
|
82 { case BIND:
|
|
|
83 done = (inc>=4) & (vr<blf->fp.nv);
|
|
|
84 break;
|
|
|
85 default:
|
|
|
86 done = (inc>=4);
|
|
|
87 }
|
|
|
88 }
|
|
|
89 }
|
|
|
90
|
|
|
91 void bsel3(h0,g0,ifact,c,cri)
|
|
|
92 double h0, g0, ifact;
|
|
|
93 int c, cri;
|
|
|
94 { double h1, g1;
|
|
|
95 int i;
|
|
|
96 hmin = h0; gmin = g0;
|
|
|
97 for (i=-1; i<=1; i++) if (i!=0)
|
|
|
98 { h1 = h0*(1+i*ifact);
|
|
|
99 g1 = bcri(h1,c,cri);
|
|
|
100 if (g1<gmin) { hmin = h1; gmin = g1; }
|
|
|
101 }
|
|
|
102 return;
|
|
|
103 }
|
|
|
104
|
|
|
105 void bselect(lf,des,c,cri,pn)
|
|
|
106 lfit *lf;
|
|
|
107 design *des;
|
|
|
108 int c, cri;
|
|
|
109 double pn;
|
|
|
110 { double h0, g0, ifact;
|
|
|
111 int i;
|
|
|
112 pen = pn;
|
|
|
113 blf = lf;
|
|
|
114 bdes = des;
|
|
|
115 if (cri==BIND) pen /= factorial(deg(&lf->sp)+1);
|
|
|
116 hmin = h0 = (c==DFXH) ? fixh(&lf->sp) : nn(&lf->sp);
|
|
|
117 if (h0==0) LERR(("bselect: initial bandwidth is 0"));
|
|
|
118 if (lf_error) return;
|
|
|
119 sig2 = 1.0;
|
|
|
120
|
|
|
121 gmin = g0 = bcri(h0,c,cri);
|
|
|
122 if (cri==BCP)
|
|
|
123 { sig2 = rv(&lf->fp);
|
|
|
124 g0 = gmin = bcri(h0,c,cri+64);
|
|
|
125 }
|
|
|
126
|
|
|
127 ifact = 0.3;
|
|
|
128 bsel2(h0,g0,ifact,c,cri);
|
|
|
129
|
|
|
130 for (i=0; i<5; i++)
|
|
|
131 { ifact = ifact/2;
|
|
|
132 bsel3(hmin,gmin,ifact,c,cri);
|
|
|
133 }
|
|
|
134 if (c==DFXH)
|
|
|
135 fixh(&lf->sp) = hmin;
|
|
|
136 else
|
|
|
137 nn(&lf->sp) = hmin;
|
|
|
138 startlf(des,lf,procvstd,0);
|
|
|
139 ressumm(lf,des,lf->fp.kap);
|
|
|
140 }
|
|
|
141
|
|
|
142 double compsda(x,h,n)
|
|
|
143 double *x, h;
|
|
|
144 int n;
|
|
|
145 /* n/(n-1) * int( fhat''(x)^2 dx ); bandwidth h */
|
|
|
146 { int i, j;
|
|
|
147 double ik, sd, z;
|
|
|
148 ik = wint(1,NULL,0,WGAUS);
|
|
|
149 sd = 0;
|
|
|
150
|
|
|
151 for (i=0; i<n; i++)
|
|
|
152 for (j=i; j<n; j++)
|
|
|
153 { z = (x[i]-x[j])/h;
|
|
|
154 sd += (2-(i==j))*Wconv4(z,WGAUS)/(ik*ik);
|
|
|
155 }
|
|
|
156 sd = sd/(n*(n-1)*h*h*h*h*h);
|
|
|
157 return(sd);
|
|
|
158 }
|
|
|
159
|
|
|
160 double widthsj(x,lambda,n)
|
|
|
161 double *x, lambda;
|
|
|
162 int n;
|
|
|
163 { double ik, a, b, td, sa, z, c, c1, c2, c3;
|
|
|
164 int i, j;
|
|
|
165 a = GFACT*0.920*lambda*exp(-log((double)n)/7)/SQRT2;
|
|
|
166 b = GFACT*0.912*lambda*exp(-log((double)n)/9)/SQRT2;
|
|
|
167 ik = wint(1,NULL,0,WGAUS);
|
|
|
168
|
|
|
169 td = 0;
|
|
|
170 for (i=0; i<n; i++)
|
|
|
171 for (j=i; j<n; j++)
|
|
|
172 { z = (x[i]-x[j])/b;
|
|
|
173 td += (2-(i==j))*Wconv6(z,WGAUS)/(ik*ik);
|
|
|
174 }
|
|
|
175
|
|
|
176 td = -td/(n*(n-1));
|
|
|
177 j = 2.0;
|
|
|
178 c1 = Wconv4(0.0,WGAUS);
|
|
|
179 c2 = wint(1,&j,1,WGAUS);
|
|
|
180 c3 = Wconv(0.0,WGAUS); /* (2*c1/(c2*c3))^(1/7)=1.357 */
|
|
|
181 sa = compsda(x,a,n);
|
|
|
182 c = b*exp(log(c1*ik/(c2*c3*GFACT*GFACT*GFACT*GFACT)*sa/td)/7)*SQRT2;
|
|
|
183 return(c);
|
|
|
184 }
|
|
|
185
|
|
|
186 void kdecri(x,h,res,c,k,ker,n)
|
|
|
187 double *x, h, *res, c;
|
|
|
188 int k, ker, n;
|
|
|
189 { int i, j;
|
|
|
190 double degfree, dfd, pen, s, r0, r1, d0, d1, ik, wij;
|
|
|
191
|
|
|
192 if (h<=0) WARN(("kdecri, h = %6.4f",h));
|
|
|
193
|
|
|
194 res[0] = res[1] = 0.0;
|
|
|
195 ik = wint(1,NULL,0,ker);
|
|
|
196 switch(k)
|
|
|
197 { case 1: /* aic */
|
|
|
198 pen = 2.0;
|
|
|
199 for (i=0; i<n; i++)
|
|
|
200 { r0 = d0 = 0.0;
|
|
|
201 for (j=0; j<n; j++)
|
|
|
202 { s = (x[i]-x[j])/h;
|
|
|
203 r0 += W(s,ker);
|
|
|
204 d0 += s*s*Wd(s,ker);
|
|
|
205 }
|
|
|
206 d0 = -(d0+r0)/(n*h*h*ik); /* d0 = d/dh fhat(xi) */
|
|
|
207 r0 /= n*h*ik; /* r0 = fhat(xi) */
|
|
|
208 res[0] += -2*log(r0)+pen*W(0.0,ker)/(n*h*ik*r0);
|
|
|
209 res[1] += -2*d0/r0-pen*W(0.0,ker)/(n*h*ik*r0)*(d0/r0+1.0/h);
|
|
|
210 }
|
|
|
211 return;
|
|
|
212 case 2: /* ocv */
|
|
|
213 for (i=0; i<n; i++)
|
|
|
214 { r0 = 0.0; d0 = 0.0;
|
|
|
215 for (j=0; j<n; j++) if (i!=j)
|
|
|
216 { s = (x[i]-x[j])/h;
|
|
|
217 r0 += W(s,ker);
|
|
|
218 d0 += s*s*Wd(s,ker);
|
|
|
219 }
|
|
|
220 d0 = -(d0+r0)/((n-1)*h*h);
|
|
|
221 r0 = r0/((n-1)*h);
|
|
|
222 res[0] -= log(r0);
|
|
|
223 res[1] -= d0/r0;
|
|
|
224 }
|
|
|
225 return;
|
|
|
226 case 3: /* lscv */
|
|
|
227 r0 = r1 = d0 = d1 = degfree = 0.0;
|
|
|
228 for (i=0; i<n; i++)
|
|
|
229 { dfd = 0.0;
|
|
|
230 for (j=0; j<n; j++)
|
|
|
231 { s = (x[i]-x[j])/h;
|
|
|
232 wij = W(s,ker);
|
|
|
233 dfd += wij;
|
|
|
234 /*
|
|
|
235 * r0 = \int fhat * fhat = sum_{i,j} W*W( (Xi-Xj)/h ) / n^2 h.
|
|
|
236 * d0 is it's derivative wrt h.
|
|
|
237 *
|
|
|
238 * r1 = 1/n sum( f_{-i}(X_i) ).
|
|
|
239 * d1 is it's derivative wrt h.
|
|
|
240 *
|
|
|
241 * degfree = sum_i ( W_0 / sum_j W( (Xi-Xj)/h ) ) is fitted d.f.
|
|
|
242 */
|
|
|
243 r0 += Wconv(s,ker);
|
|
|
244 d0 += -s*s*Wconv1(s,ker);
|
|
|
245 if (i != j)
|
|
|
246 { r1 += wij;
|
|
|
247 d1 += -s*s*Wd(s,ker);
|
|
|
248 }
|
|
|
249 }
|
|
|
250 degfree += 1.0/dfd;
|
|
|
251 }
|
|
|
252 d1 -= r1;
|
|
|
253 d0 -= r0;
|
|
|
254 res[0] = r0/(n*n*h*ik*ik) - 2*r1/(n*(n-1)*h*ik);
|
|
|
255 res[1] = d0/(n*n*h*h*ik*ik) - 2*d1/(n*(n-1)*h*h*ik);
|
|
|
256 res[2] = degfree;
|
|
|
257 return;
|
|
|
258 case 4: /* bcv */
|
|
|
259 r0 = d0 = 0.0;
|
|
|
260 for (i=0; i<n; i++)
|
|
|
261 for (j=i+1; j<n; j++)
|
|
|
262 { s = (x[i]-x[j])/h;
|
|
|
263 r0 += 2*Wconv4(s,ker);
|
|
|
264 d0 += 2*s*Wconv5(s,ker);
|
|
|
265 }
|
|
|
266 d0 = (-d0-r0)/(n*n*h*h*ik*ik);
|
|
|
267 r0 = r0/(n*n*h*ik*ik);
|
|
|
268 j = 2.0;
|
|
|
269 d1 = wint(1,&j,1,ker);
|
|
|
270 r1 = Wconv(0.0,ker);
|
|
|
271 res[0] = (d1*d1*r0/4+r1/(n*h))/(ik*ik);
|
|
|
272 res[1] = (d1*d1*d0/4-r1/(n*h*h))/(ik*ik);
|
|
|
273 return;
|
|
|
274 case 5: /* sjpi */
|
|
|
275 s = c*exp(5*log(h)/7)/SQRT2;
|
|
|
276 d0 = compsda(x,s,n);
|
|
|
277 res[0] = d0; /* this is S(alpha) in SJ */
|
|
|
278 res[1] = exp(log(Wikk(WGAUS,0)/(d0*n))/5)-h;
|
|
|
279 return;
|
|
|
280 case 6: /* gas-k-k */
|
|
|
281 s = exp(log(1.0*n)/10)*h;
|
|
|
282 d0 = compsda(x,s,n);
|
|
|
283 res[0] = d0;
|
|
|
284 res[1] = exp(log(Wikk(WGAUS,0)/(d0*n))/5)-h;
|
|
|
285 return;
|
|
|
286 }
|
|
|
287 LERR(("kdecri: what???"));
|
|
|
288 return;
|
|
|
289 }
|
|
|
290
|
|
|
291 double esolve(x,j,h0,h1,k,c,ker,n)
|
|
|
292 double *x, h0, h1, c;
|
|
|
293 int j, k, ker, n;
|
|
|
294 { double h[7], d[7], r[7], res[4], min, minh, fact;
|
|
|
295 int i, nc;
|
|
|
296 min = 1.0e30; minh = 0.0;
|
|
|
297 fact = 1.00001;
|
|
|
298 h[6] = h0; kdecri(x,h[6],res,c,j,ker,n);
|
|
|
299 r[6] = res[0]; d[6] = res[1];
|
|
|
300 if (lf_error) return(0.0);
|
|
|
301 nc = 0;
|
|
|
302 for (i=0; i<k; i++)
|
|
|
303 { h[5] = h[6]; r[5] = r[6]; d[5] = d[6];
|
|
|
304 h[6] = h0*exp((i+1)*log(h1/h0)/k);
|
|
|
305 kdecri(x,h[6],res,c,j,ker,n);
|
|
|
306 r[6] = res[0]; d[6] = res[1];
|
|
|
307 if (lf_error) return(0.0);
|
|
|
308 if (d[5]*d[6]<0)
|
|
|
309 { h[2] = h[0] = h[5]; d[2] = d[0] = d[5]; r[2] = r[0] = r[5];
|
|
|
310 h[3] = h[1] = h[6]; d[3] = d[1] = d[6]; r[3] = r[1] = r[6];
|
|
|
311 while ((h[3]>fact*h[2])|(h[2]>fact*h[3]))
|
|
|
312 { h[4] = h[3]-d[3]*(h[3]-h[2])/(d[3]-d[2]);
|
|
|
313 if ((h[4]<h[0]) | (h[4]>h[1])) h[4] = (h[0]+h[1])/2;
|
|
|
314 kdecri(x,h[4],res,c,j,ker,n);
|
|
|
315 r[4] = res[0]; d[4] = res[1];
|
|
|
316 if (lf_error) return(0.0);
|
|
|
317 h[2] = h[3]; h[3] = h[4];
|
|
|
318 d[2] = d[3]; d[3] = d[4];
|
|
|
319 r[2] = r[3]; r[3] = r[4];
|
|
|
320 if (d[4]*d[0]>0) { h[0] = h[4]; d[0] = d[4]; r[0] = r[4]; }
|
|
|
321 else { h[1] = h[4]; d[1] = d[4]; r[1] = r[4]; }
|
|
|
322 }
|
|
|
323 if (j>=4) return(h[4]); /* first min for BCV etc */
|
|
|
324 if (r[4]<=min) { min = r[4]; minh = h[4]; }
|
|
|
325 nc++;
|
|
|
326 }
|
|
|
327 }
|
|
|
328 if (nc==0) minh = (r[5]<r[6]) ? h0 : h1;
|
|
|
329 return(minh);
|
|
|
330 }
|
|
|
331
|
|
|
332 void kdeselect(band,x,ind,h0,h1,meth,nm,ker,n)
|
|
|
333 double h0, h1, *band, *x;
|
|
|
334 int *ind, nm, ker, n, *meth;
|
|
|
335 { double scale, c;
|
|
|
336 int i, k;
|
|
|
337 k = n/4;
|
|
|
338 for (i=0; i<n; i++) ind[i] = i;
|
|
|
339 scale = kordstat(x,n+1-k,n,ind) - kordstat(x,k,n,ind);
|
|
|
340 c = widthsj(x,scale,n);
|
|
|
341 for (i=0; i<nm; i++)
|
|
|
342 band[i] = esolve(x,meth[i],h0,h1,10,c,ker,n);
|
|
|
343 }
|
|
|
344 /*
|
|
|
345 * Copyright 1996-2006 Catherine Loader.
|
|
|
346 */
|
|
|
347 /*
|
|
|
348 * The function dens_integrate(lf,des,z) is used to integrate a density
|
|
|
349 * estimate (z=1) or the density squared (z=2). This is used to renormalize
|
|
|
350 * the estimate (function dens_renorm) or in the computation of LSCV
|
|
|
351 * (in modlscv.c). The implementation is presently for d=1.
|
|
|
352 *
|
|
|
353 * The computation orders the fit points selected by locfit, and
|
|
|
354 * integrates analytically over each interval. For the log-link,
|
|
|
355 * the interpolant used is peicewise quadratic (with one knot in
|
|
|
356 * the middle of each interval); this differs from the cubic interpolant
|
|
|
357 * used elsewhere in Locfit.
|
|
|
358 *
|
|
|
359 * TODO: allow for xlim. What can be done simply in >=2 dimensions?
|
|
|
360 */
|
|
|
361
|
|
|
362 #include "lfev.h"
|
|
|
363
|
|
|
364 /*
|
|
|
365 * Finds the order of observations in the array x, and
|
|
|
366 * stores in integer array ind.
|
|
|
367 * At input, lset l=0 and r=length(x)-1.
|
|
|
368 * At output, x[ind[0]] <= x[ind[1]] <= ...
|
|
|
369 */
|
|
|
370 void lforder(ind,x,l,r)
|
|
|
371 int *ind, l, r;
|
|
|
372 double *x;
|
|
|
373 { double piv;
|
|
|
374 int i, i0, i1;
|
|
|
375 piv = (x[ind[l]]+x[ind[r]])/2;
|
|
|
376 i0 = l; i1 = r;
|
|
|
377 while (i0<=i1)
|
|
|
378 { while ((i0<=i1) && (x[ind[i0]]<=piv)) i0++;
|
|
|
379 while ((i0<=i1) && (x[ind[i1]]>piv)) i1--;
|
|
|
380 if (i0<i1)
|
|
|
381 { ISWAP(ind[i0],ind[i1]);
|
|
|
382 i0++; i1--;
|
|
|
383 }
|
|
|
384 }
|
|
|
385 /* now, x[ind[l..i1]] <= piv < x[ind[i0..r]].
|
|
|
386 put the ties in the middle */
|
|
|
387 while ((i1>=l) && (x[ind[i1]]==piv)) i1--;
|
|
|
388 for (i=l; i<=i1; i++)
|
|
|
389 if (x[ind[i]]==piv)
|
|
|
390 { ISWAP(ind[i],ind[i1]);
|
|
|
391 while (x[ind[i1]]==piv) i1--;
|
|
|
392 }
|
|
|
393
|
|
|
394 if (l<i1) lforder(ind,x,l,i1);
|
|
|
395 if (i0<r) lforder(ind,x,i0,r);
|
|
|
396 }
|
|
|
397
|
|
|
398 /*
|
|
|
399 * estdiv integrates the density between fit points x0 and x1.
|
|
|
400 * f0, f1 are function values, d0, d1 are derivatives.
|
|
|
401 */
|
|
|
402 double estdiv(x0,x1,f0,f1,d0,d1,lin)
|
|
|
403 double x0, x1, f0, f1, d0, d1;
|
|
|
404 int lin;
|
|
|
405 { double cf[4], I[2], dlt, e0, e1;
|
|
|
406
|
|
|
407 if (x0==x1) return(0.0);
|
|
|
408
|
|
|
409 if (lin==LIDENT)
|
|
|
410 {
|
|
|
411 /* cf are integrals of hermite polynomials.
|
|
|
412 * Then adjust for x1-x0.
|
|
|
413 */
|
|
|
414 cf[0] = cf[1] = 0.5;
|
|
|
415 cf[2] = 1.0/12.0; cf[3] = -cf[2];
|
|
|
416 return( (cf[0]*f0+cf[1]*f1)*(x1-x0)
|
|
|
417 + (cf[2]*d0+cf[3]*d1)*(x1-x0)*(x1-x0) );
|
|
|
418 }
|
|
|
419
|
|
|
420 /*
|
|
|
421 * this is for LLOG
|
|
|
422 */
|
|
|
423
|
|
|
424 dlt = (x1-x0)/2;
|
|
|
425 cf[0] = f0;
|
|
|
426 cf[1] = d0;
|
|
|
427 cf[2] = ( 2*(f1-f0) - dlt*(d1+3*d0) )/(4*dlt*dlt);
|
|
|
428 recurint(0.0,dlt,cf,I,0,WRECT);
|
|
|
429 e0 = I[0];
|
|
|
430
|
|
|
431 cf[0] = f1;
|
|
|
432 cf[1] = -d1;
|
|
|
433 cf[2] = ( 2*(f0-f1) + dlt*(d0+3*d1) )/( 4*dlt*dlt );
|
|
|
434 recurint(0.0,dlt,cf,I,0,WRECT);
|
|
|
435 e1 = I[0];
|
|
|
436
|
|
|
437 return(e0+e1);
|
|
|
438 }
|
|
|
439
|
|
|
440 /*
|
|
|
441 * Evaluate the integral of the density estimate to the power z.
|
|
|
442 * This would be severely messed up, if I ever implement parcomp
|
|
|
443 * for densities.
|
|
|
444 */
|
|
|
445 double dens_integrate(lf,des,z)
|
|
|
446 lfit *lf;
|
|
|
447 design *des;
|
|
|
448 int z;
|
|
|
449 { int has_deriv, i, i0, i1, nv, *ind;
|
|
|
450 double *xev, *fit, *deriv, sum, term;
|
|
|
451 double d0, d1, f0, f1;
|
|
|
452 fitpt *fp;
|
|
|
453
|
|
|
454 fp = &lf->fp;
|
|
|
455
|
|
|
456 if (fp->d >= 2)
|
|
|
457 { WARN(("dens_integrate requires d=1"));
|
|
|
458 return(0.0);
|
|
|
459 }
|
|
|
460
|
|
|
461 has_deriv = (deg(&lf->sp) > 0); /* not right? */
|
|
|
462 fit = fp->coef;
|
|
|
463 if (has_deriv)
|
|
|
464 deriv = &fit[fp->nvm];
|
|
|
465 xev = evp(fp);
|
|
|
466
|
|
|
467 /*
|
|
|
468 * order the vertices
|
|
|
469 */
|
|
|
470 nv = fp->nv;
|
|
|
471 if (lf->lfd.n<nv) return(0.0);
|
|
|
472 ind = des->ind;
|
|
|
473 for (i=0; i<nv; i++) ind[i] = i;
|
|
|
474 lforder(ind,xev,0,nv-1);
|
|
|
475 sum = 0.0;
|
|
|
476
|
|
|
477 /*
|
|
|
478 * Estimate the contribution of the boundaries.
|
|
|
479 * should really check flim here.
|
|
|
480 */
|
|
|
481 i0 = ind[0]; i1 = ind[1];
|
|
|
482 f1 = fit[i0];
|
|
|
483 d1 = (has_deriv) ? deriv[i0] :
|
|
|
484 (fit[i1]-fit[i0])/(xev[i1]-xev[i0]);
|
|
|
485 if (d1 <= 0) WARN(("dens_integrate - ouch!"));
|
|
|
486 if (z==2)
|
|
|
487 { if (link(&lf->sp)==LLOG)
|
|
|
488 { f1 *= 2; d1 *= 2; }
|
|
|
489 else
|
|
|
490 { d1 = 2*d1*f1; f1 = f1*f1; }
|
|
|
491 }
|
|
|
492 term = (link(&lf->sp)==LIDENT) ? f1*f1/(2*d1) : exp(f1)/d1;
|
|
|
493 sum += term;
|
|
|
494
|
|
|
495 i0 = ind[nv-2]; i1 = ind[nv-1];
|
|
|
496 f0 = fit[i1];
|
|
|
497 d0 = (has_deriv) ? deriv[i1] :
|
|
|
498 (fit[i1]-fit[i0])/(xev[i1]-xev[i0]);
|
|
|
499 if (d0 >= 0) WARN(("dens_integrate - ouch!"));
|
|
|
500 if (z==2)
|
|
|
501 { if (link(&lf->sp)==LLOG)
|
|
|
502 { f0 *= 2; d0 *= 2; }
|
|
|
503 else
|
|
|
504 { d0 = 2*d0*f0; f0 = f0*f0; }
|
|
|
505 }
|
|
|
506 term = (link(&lf->sp)==LIDENT) ? -f0*f0/(2*d0) : exp(f0)/d0;
|
|
|
507 sum += term;
|
|
|
508
|
|
|
509 for (i=1; i<nv; i++)
|
|
|
510 { i0 = ind[i-1]; i1 = ind[i];
|
|
|
511 f0 = fit[i0]; f1 = fit[i1];
|
|
|
512 d0 = (has_deriv) ? deriv[i0] :
|
|
|
513 (f1-f0)/(xev[i1]-xev[i0]);
|
|
|
514 d1 = (has_deriv) ? deriv[i1] : d0;
|
|
|
515 if (z==2)
|
|
|
516 { if (link(&lf->sp)==LLOG)
|
|
|
517 { f0 *= 2; f1 *= 2; d0 *= 2; d1 *= 2; }
|
|
|
518 else
|
|
|
519 { d0 *= 2*f0; d1 *= 2*f1; f0 = f0*f0; f1 = f1*f1; }
|
|
|
520 }
|
|
|
521 term = estdiv(xev[i0],xev[i1],f0,f1,d0,d1,link(&lf->sp));
|
|
|
522 sum += term;
|
|
|
523 }
|
|
|
524
|
|
|
525 return(sum);
|
|
|
526 }
|
|
|
527
|
|
|
528 void dens_renorm(lf,des)
|
|
|
529 lfit *lf;
|
|
|
530 design *des;
|
|
|
531 { int i;
|
|
|
532 double sum;
|
|
|
533 sum = dens_integrate(lf,des,1);
|
|
|
534 if (sum==0.0) return;
|
|
|
535 sum = log(sum);
|
|
|
536 for (i=0; i<lf->fp.nv; i++) lf->fp.coef[i] -= sum;
|
|
|
537 }
|
|
|
538 /*
|
|
|
539 * Copyright 1996-2006 Catherine Loader.
|
|
|
540 */
|
|
|
541 /*
|
|
|
542 * This file contains functions for constructing and
|
|
|
543 * interpolating the adaptive tree structure. This is
|
|
|
544 * the default evaluation structure used by Locfit.
|
|
|
545 */
|
|
|
546
|
|
|
547 #include "lfev.h"
|
|
|
548
|
|
|
549 /*
|
|
|
550 Guess the number of fitting points.
|
|
|
551 Needs improving!
|
|
|
552 */
|
|
|
553 void atree_guessnv(evs,nvm,ncm,vc,d,alp)
|
|
|
554 evstruc *evs;
|
|
|
555 double alp;
|
|
|
556 int *nvm, *ncm, *vc, d;
|
|
|
557 { double a0, cu, ifl;
|
|
|
558 int i, nv, nc;
|
|
|
559
|
|
|
560 *ncm = 1<<30; *nvm = 1<<30;
|
|
|
561 *vc = 1 << d;
|
|
|
562
|
|
|
563 if (alp>0)
|
|
|
564 { a0 = (alp > 1) ? 1 : 1/alp;
|
|
|
565 if (cut(evs)<0.01)
|
|
|
566 { WARN(("guessnv: cut too small."));
|
|
|
567 cut(evs) = 0.01;
|
|
|
568 }
|
|
|
569 cu = 1;
|
|
|
570 for (i=0; i<d; i++) cu *= MIN(1.0,cut(evs));
|
|
|
571 nv = (int)((5*a0/cu)**vc); /* this allows 10*a0/cu splits */
|
|
|
572 nc = (int)(10*a0/cu+1); /* and 10*a0/cu cells */
|
|
|
573 if (nv<*nvm) *nvm = nv;
|
|
|
574 if (nc<*ncm) *ncm = nc;
|
|
|
575 }
|
|
|
576
|
|
|
577 if (*nvm == 1<<30) /* by default, allow 100 splits */
|
|
|
578 { *nvm = 102**vc;
|
|
|
579 *ncm = 201;
|
|
|
580 }
|
|
|
581
|
|
|
582 /* inflation based on mk */
|
|
|
583 ifl = mk(evs)/100.0;
|
|
|
584 *nvm = *vc+(int)(ifl**nvm);
|
|
|
585 *ncm = (int)(ifl**ncm);
|
|
|
586
|
|
|
587 }
|
|
|
588
|
|
|
589 /*
|
|
|
590 Determine whether a cell in the tree needs splitting.
|
|
|
591 If so, return the split variable (0..d-1).
|
|
|
592 Otherwise, return -1.
|
|
|
593 */
|
|
|
594 int atree_split(lf,ce,le,ll,ur)
|
|
|
595 lfit *lf;
|
|
|
596 int *ce;
|
|
|
597 double *le, *ll, *ur;
|
|
|
598 { int d, vc, i, is;
|
|
|
599 double h, hmin, score[MXDIM];
|
|
|
600 d = lf->fp.d; vc = 1<<d;
|
|
|
601
|
|
|
602 hmin = 0.0;
|
|
|
603 for (i=0; i<vc; i++)
|
|
|
604 { h = lf->fp.h[ce[i]];
|
|
|
605 if ((h>0) && ((hmin==0)|(h<hmin))) hmin = h;
|
|
|
606 }
|
|
|
607
|
|
|
608 is = 0;
|
|
|
609 for (i=0; i<d; i++)
|
|
|
610 { le[i] = (ur[i]-ll[i])/lf->lfd.sca[i];
|
|
|
611 if ((lf->lfd.sty[i]==STCPAR) || (hmin==0))
|
|
|
612 score[i] = 2*(ur[i]-ll[i])/(lf->evs.fl[i+d]-lf->evs.fl[i]);
|
|
|
613 else
|
|
|
614 score[i] = le[i]/hmin;
|
|
|
615 if (score[i]>score[is]) is = i;
|
|
|
616 }
|
|
|
617 if (cut(&lf->evs)<score[is]) return(is);
|
|
|
618 return(-1);
|
|
|
619 }
|
|
|
620
|
|
|
621 /*
|
|
|
622 recursively grow the tree structure, begining with the parent cell.
|
|
|
623 */
|
|
|
624 void atree_grow(des,lf,ce,ct,term,ll,ur)
|
|
|
625 design *des;
|
|
|
626 lfit *lf;
|
|
|
627 int *ce, *ct, *term;
|
|
|
628 double *ll, *ur;
|
|
|
629 { int nce[1<<MXDIM];
|
|
|
630 int i, i0, i1, d, vc, pv, tk, ns;
|
|
|
631 double le[MXDIM], z;
|
|
|
632 d = lf->fp.d; vc = 1<<d;
|
|
|
633
|
|
|
634 /* does this cell need splitting?
|
|
|
635 If not, wrap up and return.
|
|
|
636 */
|
|
|
637 ns = atree_split(lf,ce,le,ll,ur);
|
|
|
638 if (ns==-1)
|
|
|
639 { if (ct != NULL) /* reconstructing terminal cells */
|
|
|
640 { for (i=0; i<vc; i++) term[*ct*vc+i] = ce[i];
|
|
|
641 (*ct)++;
|
|
|
642 }
|
|
|
643 return;
|
|
|
644 }
|
|
|
645
|
|
|
646 /* split the cell at the midpoint on side ns */
|
|
|
647 tk = 1<<ns;
|
|
|
648 for (i=0; i<vc; i++)
|
|
|
649 { if ((i&tk)==0) nce[i] = ce[i];
|
|
|
650 else
|
|
|
651 { i0 = ce[i];
|
|
|
652 i1 = ce[i-tk];
|
|
|
653 pv = (lf->lfd.sty[i]!=STCPAR) &&
|
|
|
654 (le[ns] < (cut(&lf->evs)*MIN(lf->fp.h[i0],lf->fp.h[i1])));
|
|
|
655 nce[i] = newsplit(des,lf,i0,i1,pv);
|
|
|
656 if (lf_error) return;
|
|
|
657 }
|
|
|
658 }
|
|
|
659 z = ur[ns]; ur[ns] = (z+ll[ns])/2;
|
|
|
660 atree_grow(des,lf,nce,ct,term,ll,ur);
|
|
|
661 if (lf_error) return;
|
|
|
662 ur[ns] = z;
|
|
|
663 for (i=0; i<vc; i++)
|
|
|
664 nce[i] = ((i&tk)== 0) ? nce[i+tk] : ce[i];
|
|
|
665 z = ll[ns]; ll[ns] = (z+ur[ns])/2;
|
|
|
666 atree_grow(des,lf,nce,ct,term,ll,ur);
|
|
|
667 if (lf_error) return;
|
|
|
668 ll[ns] = z;
|
|
|
669 }
|
|
|
670
|
|
|
671 void atree_start(des,lf)
|
|
|
672 design *des;
|
|
|
673 lfit *lf;
|
|
|
674 { int d, i, j, k, vc, ncm, nvm;
|
|
|
675 double ll[MXDIM], ur[MXDIM];
|
|
|
676
|
|
|
677 if (lf_debug>1) mut_printf(" In atree_start\n");
|
|
|
678 d = lf->fp.d;
|
|
|
679 atree_guessnv(&lf->evs,&nvm,&ncm,&vc,d,nn(&lf->sp));
|
|
|
680 if (lf_debug>2) mut_printf(" atree_start: nvm %d ncm %d\n",nvm,ncm);
|
|
|
681 trchck(lf,nvm,ncm,vc);
|
|
|
682
|
|
|
683 /* Set the lower left, upper right limits. */
|
|
|
684 for (j=0; j<d; j++)
|
|
|
685 { ll[j] = lf->evs.fl[j];
|
|
|
686 ur[j] = lf->evs.fl[j+d];
|
|
|
687 }
|
|
|
688
|
|
|
689 /* Set the initial cell; fit at the vertices. */
|
|
|
690 for (i=0; i<vc; i++)
|
|
|
691 { j = i;
|
|
|
692 for (k=0; k<d; ++k)
|
|
|
693 { evptx(&lf->fp,i,k) = (j%2) ? ur[k] : ll[k];
|
|
|
694 j >>= 1;
|
|
|
695 }
|
|
|
696 lf->evs.ce[i] = i;
|
|
|
697 PROC_VERTEX(des,lf,i);
|
|
|
698 if (lf_error) return;
|
|
|
699 lf->evs.s[i] = 0;
|
|
|
700 }
|
|
|
701 lf->fp.nv = vc;
|
|
|
702
|
|
|
703 /* build the tree */
|
|
|
704 atree_grow(des,lf,lf->evs.ce,NULL,NULL,ll,ur);
|
|
|
705 lf->evs.nce = 1;
|
|
|
706 }
|
|
|
707
|
|
|
708 double atree_int(lf,x,what)
|
|
|
709 lfit *lf;
|
|
|
710 double *x;
|
|
|
711 int what;
|
|
|
712 { double vv[64][64], *ll, *ur, h, xx[MXDIM];
|
|
|
713 int lo, tk, ns, nv, nc, d, i, vc;
|
|
|
714 int ce[64];
|
|
|
715 fitpt *fp;
|
|
|
716 evstruc *evs;
|
|
|
717
|
|
|
718 fp = &lf->fp;
|
|
|
719 evs= &lf->evs;
|
|
|
720 d = fp->d;
|
|
|
721 vc = 1<<d;
|
|
|
722
|
|
|
723 for (i=0; i<vc; i++)
|
|
|
724 { setzero(vv[i],vc);
|
|
|
725 nc = exvval(fp,vv[i],i,d,what,1);
|
|
|
726 ce[i] = evs->ce[i];
|
|
|
727 }
|
|
|
728 ns = 0;
|
|
|
729 while(ns!=-1)
|
|
|
730 { ll = evpt(fp,ce[0]); ur = evpt(fp,ce[vc-1]);
|
|
|
731 ns = atree_split(lf,ce,xx,ll,ur);
|
|
|
732 if (ns!=-1)
|
|
|
733 { tk = 1<<ns;
|
|
|
734 h = ur[ns]-ll[ns];
|
|
|
735 lo = (2*(x[ns]-ll[ns])) < h;
|
|
|
736 for (i=0; i<vc; i++) if ((tk&i)==0)
|
|
|
737 { nv = findpt(fp,evs,(int)ce[i],(int)ce[i+tk]);
|
|
|
738 if (nv==-1) LERR(("Descend tree problem"));
|
|
|
739 if (lf_error) return(0.0);
|
|
|
740 if (lo)
|
|
|
741 { ce[i+tk] = nv;
|
|
|
742 if (evs->s[nv]) exvvalpv(vv[i+tk],vv[i],vv[i+tk],d,ns,h,nc);
|
|
|
743 else exvval(fp,vv[i+tk],nv,d,what,1);
|
|
|
744 }
|
|
|
745 else
|
|
|
746 { ce[i] = nv;
|
|
|
747 if (evs->s[nv]) exvvalpv(vv[i],vv[i],vv[i+tk],d,ns,h,nc);
|
|
|
748 else exvval(fp,vv[i],nv,d,what,1);
|
|
|
749 } }
|
|
|
750 } }
|
|
|
751 ll = evpt(fp,ce[0]); ur = evpt(fp,ce[vc-1]);
|
|
|
752 return(rectcell_interp(x,vv,ll,ur,d,nc));
|
|
|
753 }
|
|
|
754 /*
|
|
|
755 * Copyright 1996-2006 Catherine Loader.
|
|
|
756 */
|
|
|
757 #include "lfev.h"
|
|
|
758
|
|
|
759 double linear_interp(h,d,f0,f1)
|
|
|
760 double h, d, f0, f1;
|
|
|
761 { if (d==0) return(f0);
|
|
|
762 return( ( (d-h)*f0 + h*f1 ) / d );
|
|
|
763 }
|
|
|
764
|
|
|
765 void hermite2(x,z,phi)
|
|
|
766 double x, z, *phi;
|
|
|
767 { double h;
|
|
|
768 if (z==0)
|
|
|
769 { phi[0] = 1.0; phi[1] = phi[2] = phi[3] = 0.0;
|
|
|
770 return;
|
|
|
771 }
|
|
|
772 h = x/z;
|
|
|
773 if (h<0)
|
|
|
774 { phi[0] = 1; phi[1] = 0;
|
|
|
775 phi[2] = h; phi[3] = 0;
|
|
|
776 return;
|
|
|
777 }
|
|
|
778 if (h>1)
|
|
|
779 { phi[0] = 0; phi[1] = 1;
|
|
|
780 phi[2] = 0; phi[3] = h-1;
|
|
|
781 return;
|
|
|
782 }
|
|
|
783 phi[1] = h*h*(3-2*h);
|
|
|
784 phi[0] = 1-phi[1];
|
|
|
785 phi[2] = h*(1-h)*(1-h);
|
|
|
786 phi[3] = h*h*(h - 1);
|
|
|
787 }
|
|
|
788
|
|
|
789 double cubic_interp(h,f0,f1,d0,d1)
|
|
|
790 double h, f0, f1, d0, d1;
|
|
|
791 { double phi[4];
|
|
|
792 hermite2(h,1.0,phi);
|
|
|
793 return(phi[0]*f0+phi[1]*f1+phi[2]*d0+phi[3]*d1);
|
|
|
794 }
|
|
|
795
|
|
|
796 double cubintd(h,f0,f1,d0,d1)
|
|
|
797 double h, f0, f1, d0, d1;
|
|
|
798 { double phi[4];
|
|
|
799 phi[1] = 6*h*(1-h);
|
|
|
800 phi[0] = -phi[1];
|
|
|
801 phi[2] = (1-h)*(1-3*h);
|
|
|
802 phi[3] = h*(3*h-2);
|
|
|
803 return(phi[0]*f0+phi[1]*f1+phi[2]*d0+phi[3]*d1);
|
|
|
804 }
|
|
|
805
|
|
|
806 /*
|
|
|
807 interpolate over a rectangular cell.
|
|
|
808 x = interpolation point.
|
|
|
809 vv = array of vertex values.
|
|
|
810 ll = lower left corner.
|
|
|
811 ur = upper right corner.
|
|
|
812 d = dimension.
|
|
|
813 nc = no of coefficients.
|
|
|
814 */
|
|
|
815 double rectcell_interp(x,vv,ll,ur,d,nc)
|
|
|
816 double *x, vv[64][64], *ll, *ur;
|
|
|
817 int d, nc;
|
|
|
818 { double phi[4];
|
|
|
819 int i, j, k, tk;
|
|
|
820
|
|
|
821 tk = 1<<d;
|
|
|
822 for (i=0; i<tk; i++) if (vv[i][0]==NOSLN) return(NOSLN);
|
|
|
823
|
|
|
824 /* no derivatives - use multilinear interpolation */
|
|
|
825 if (nc==1)
|
|
|
826 { for (i=d-1; i>=0; i--)
|
|
|
827 { tk = 1<<i;
|
|
|
828 for (j=0; j<tk; j++)
|
|
|
829 vv[j][0] = linear_interp(x[i]-ll[i],ur[i]-ll[i],vv[j][0],vv[j+tk][0]);
|
|
|
830 }
|
|
|
831 return(vv[0][0]);
|
|
|
832 }
|
|
|
833
|
|
|
834 /* with derivatives -- use cubic */
|
|
|
835 if (nc==d+1)
|
|
|
836 { for (i=d-1; i>=0; i--)
|
|
|
837 { hermite2(x[i]-ll[i],ur[i]-ll[i],phi);
|
|
|
838 tk = 1<<i;
|
|
|
839 phi[2] *= ur[i]-ll[i];
|
|
|
840 phi[3] *= ur[i]-ll[i];
|
|
|
841 for (j=0; j<tk; j++)
|
|
|
842 { vv[j][0] = phi[0]*vv[j][0] + phi[1]*vv[j+tk][0]
|
|
|
843 + phi[2]*vv[j][i+1] + phi[3]*vv[j+tk][i+1];
|
|
|
844 for (k=1; k<=i; k++)
|
|
|
845 vv[j][k] = phi[0]*vv[j][k] + phi[1]*vv[j+tk][k];
|
|
|
846 }
|
|
|
847 }
|
|
|
848 return(vv[0][0]);
|
|
|
849 }
|
|
|
850
|
|
|
851 /* with all coefs -- use multicubic */
|
|
|
852 for (i=d-1; i>=0; i--)
|
|
|
853 { hermite2(x[i]-ll[i],ur[i]-ll[i],phi);
|
|
|
854 tk = 1<<i;
|
|
|
855 phi[2] *= ur[i]-ll[i];
|
|
|
856 phi[3] *= ur[i]-ll[i];
|
|
|
857 for (j=0; j<tk; j++)
|
|
|
858 for (k=0; k<tk; k++)
|
|
|
859 vv[j][k] = phi[0]*vv[j][k] + phi[1]*vv[j+tk][k]
|
|
|
860 + phi[2]*vv[j][k+tk] + phi[3]*vv[j+tk][k+tk];
|
|
|
861 }
|
|
|
862 return(vv[0][0]);
|
|
|
863 }
|
|
|
864
|
|
|
865 int exvval(fp,vv,nv,d,what,z)
|
|
|
866 fitpt *fp;
|
|
|
867 double *vv;
|
|
|
868 int nv, d, z, what;
|
|
|
869 { int i, k;
|
|
|
870 double *values;
|
|
|
871
|
|
|
872 k = (z) ? 1<<d : d+1;
|
|
|
873 for (i=1; i<k; i++) vv[i] = 0.0;
|
|
|
874 switch(what)
|
|
|
875 { case PCOEF:
|
|
|
876 values = fp->coef;
|
|
|
877 break;
|
|
|
878 case PVARI:
|
|
|
879 case PNLX:
|
|
|
880 values = fp->nlx;
|
|
|
881 break;
|
|
|
882 case PT0:
|
|
|
883 values = fp->t0;
|
|
|
884 break;
|
|
|
885 case PBAND:
|
|
|
886 vv[0] = fp->h[nv];
|
|
|
887 return(1);
|
|
|
888 case PDEGR:
|
|
|
889 vv[0] = fp->deg[nv];
|
|
|
890 return(1);
|
|
|
891 case PLIK:
|
|
|
892 vv[0] = fp->lik[nv];
|
|
|
893 return(1);
|
|
|
894 case PRDF:
|
|
|
895 vv[0] = fp->lik[2*fp->nvm+nv];
|
|
|
896 return(1);
|
|
|
897 default:
|
|
|
898 LERR(("Invalid what in exvval"));
|
|
|
899 return(0);
|
|
|
900 }
|
|
|
901 vv[0] = values[nv];
|
|
|
902 if (!fp->hasd) return(1);
|
|
|
903 if (z)
|
|
|
904 { for (i=0; i<d; i++) vv[1<<i] = values[(i+1)*fp->nvm+nv];
|
|
|
905 return(1<<d);
|
|
|
906 }
|
|
|
907 else
|
|
|
908 { for (i=1; i<=d; i++) vv[i] = values[i*fp->nvm+nv];
|
|
|
909 return(d+1);
|
|
|
910 }
|
|
|
911 }
|
|
|
912
|
|
|
913 void exvvalpv(vv,vl,vr,d,k,dl,nc)
|
|
|
914 double *vv, *vl, *vr, dl;
|
|
|
915 int d, k, nc;
|
|
|
916 { int i, tk, td;
|
|
|
917 double f0, f1;
|
|
|
918 if (nc==1)
|
|
|
919 { vv[0] = (vl[0]+vr[0])/2;
|
|
|
920 return;
|
|
|
921 }
|
|
|
922 tk = 1<<k;
|
|
|
923 td = 1<<d;
|
|
|
924 for (i=0; i<td; i++) if ((i&tk)==0)
|
|
|
925 { f0 = (vl[i]+vr[i])/2 + dl*(vl[i+tk]-vr[i+tk])/8;
|
|
|
926 f1 = 1.5*(vr[i]-vl[i])/dl - (vl[i+tk]+vr[i+tk])/4;
|
|
|
927 vv[i] = f0;
|
|
|
928 vv[i+tk] = f1;
|
|
|
929 }
|
|
|
930 }
|
|
|
931
|
|
|
932 double grid_int(fp,evs,x,what)
|
|
|
933 fitpt *fp;
|
|
|
934 evstruc *evs;
|
|
|
935 double *x;
|
|
|
936 int what;
|
|
|
937 { int d, i, j, jj, nc, sk, v[MXDIM], vc, z0, nce[1<<MXDIM], *mg;
|
|
|
938 double *ll, *ur, vv[64][64], z;
|
|
|
939
|
|
|
940 d = fp->d;
|
|
|
941 ll = evpt(fp,0); ur = evpt(fp,fp->nv-1);
|
|
|
942 mg = mg(evs);
|
|
|
943
|
|
|
944 z0 = 0; vc = 1<<d;
|
|
|
945 for (j=d-1; j>=0; j--)
|
|
|
946 { v[j] = (int)((mg[j]-1)*(x[j]-ll[j])/(ur[j]-ll[j]));
|
|
|
947 if (v[j]<0) v[j]=0;
|
|
|
948 if (v[j]>=mg[j]-1) v[j] = mg[j]-2;
|
|
|
949 z0 = z0*mg[j]+v[j];
|
|
|
950 }
|
|
|
951 nce[0] = z0; nce[1] = z0+1; sk = jj = 1;
|
|
|
952 for (i=1; i<d; i++)
|
|
|
953 { sk *= mg[i-1];
|
|
|
954 jj<<=1;
|
|
|
955 for (j=0; j<jj; j++)
|
|
|
956 nce[j+jj] = nce[j]+sk;
|
|
|
957 }
|
|
|
958 for (i=0; i<vc; i++)
|
|
|
959 nc = exvval(fp,vv[i],nce[i],d,what,1);
|
|
|
960 ll = evpt(fp,nce[0]);
|
|
|
961 ur = evpt(fp,nce[vc-1]);
|
|
|
962 z = rectcell_interp(x,vv,ll,ur,d,nc);
|
|
|
963 return(z);
|
|
|
964 }
|
|
|
965
|
|
|
966 double fitp_int(fp,x,what,i)
|
|
|
967 fitpt *fp;
|
|
|
968 double *x;
|
|
|
969 int what, i;
|
|
|
970 { double vv[1+MXDIM];
|
|
|
971 exvval(fp,vv,i,fp->d,what,0);
|
|
|
972 return(vv[0]);
|
|
|
973 }
|
|
|
974
|
|
|
975 double xbar_int(fp,x,what)
|
|
|
976 fitpt *fp;
|
|
|
977 double *x;
|
|
|
978 int what;
|
|
|
979 { int i, nc;
|
|
|
980 double vv[1+MXDIM], f;
|
|
|
981 nc = exvval(fp,vv,0,fp->d,what,0);
|
|
|
982 f = vv[0];
|
|
|
983 if (nc>1)
|
|
|
984 for (i=0; i<fp->d; i++)
|
|
|
985 f += vv[i+1]*(x[i]-evptx(fp,0,i));
|
|
|
986 return(f);
|
|
|
987 }
|
|
|
988
|
|
|
989 double dointpoint(lf,x,what,ev,j)
|
|
|
990 lfit *lf;
|
|
|
991 double *x;
|
|
|
992 int what, ev, j;
|
|
|
993 { double xf, f, link[LLEN];
|
|
|
994 int i, rt;
|
|
|
995 fitpt *fp;
|
|
|
996 evstruc *evs;
|
|
|
997
|
|
|
998 fp = &lf->fp; evs = &lf->evs;
|
|
|
999 for (i=0; i<fp->d; i++) if (lf->lfd.sty[i]==STANGL)
|
|
|
1000 { xf = floor(x[i]/(2*PI*lf->lfd.sca[i]));
|
|
|
1001 x[i] -= xf*2*PI*lf->lfd.sca[i];
|
|
|
1002 }
|
|
|
1003 if (what > 64)
|
|
|
1004 { rt = what-64;
|
|
|
1005 what = PCOEF;
|
|
|
1006 }
|
|
|
1007 else rt = 0;
|
|
|
1008
|
|
|
1009 switch(ev)
|
|
|
1010 { case EGRID: f = grid_int(fp,evs,x,what); break;
|
|
|
1011 case EKDTR: f = kdtre_int(fp,evs,x,what); break;
|
|
|
1012 case ETREE: f = atree_int(lf,x,what); break;
|
|
|
1013 case EPHULL: f = triang_int(lf,x,what); break;
|
|
|
1014 case EFITP: f = fitp_int(fp,x,what,j); break;
|
|
|
1015 case EXBAR: f = xbar_int(fp,x,what); break;
|
|
|
1016 case ENONE: f = 0; break;
|
|
|
1017 case ESPHR: f = sphere_int(lf,x,what); break;
|
|
|
1018 default: LERR(("dointpoint: cannot interpolate structure %d",ev));
|
|
|
1019 }
|
|
|
1020 if (((what==PT0)|(what==PNLX)) && (f<0)) f = 0.0;
|
|
|
1021 f += addparcomp(lf,x,what);
|
|
|
1022
|
|
|
1023 if (rt>0)
|
|
|
1024 {
|
|
|
1025 stdlinks(link,&lf->lfd,&lf->sp,j,f,rsc(fp));
|
|
|
1026 f = resid(resp(&lf->lfd,j),prwt(&lf->lfd,j),f,fam(&lf->sp),rt,link);
|
|
|
1027 }
|
|
|
1028
|
|
|
1029 return(f);
|
|
|
1030 }
|
|
|
1031 /*
|
|
|
1032 * Copyright 1996-2006 Catherine Loader.
|
|
|
1033 */
|
|
|
1034 /*
|
|
|
1035 * Routines for building and interpolating the kd tree.
|
|
|
1036 * Initially, this started from the loess code.
|
|
|
1037 *
|
|
|
1038 * Todo: EKDCE isn't working.
|
|
|
1039 */
|
|
|
1040
|
|
|
1041 #include "lfev.h"
|
|
|
1042
|
|
|
1043 void newcell();
|
|
|
1044 static int nterm;
|
|
|
1045
|
|
|
1046 void kdtre_guessnv(evs,nvm,ncm,vc,n,d,alp)
|
|
|
1047 evstruc *evs;
|
|
|
1048 double alp;
|
|
|
1049 int *nvm, *ncm, *vc, n, d;
|
|
|
1050 { int k;
|
|
|
1051 if (ev(evs) == EKDTR)
|
|
|
1052 { nterm = (int)(cut(evs)/4 * n * MIN(alp,1.0) );
|
|
|
1053 k = 2*n/nterm;
|
|
|
1054 *vc = 1<<d;
|
|
|
1055 *ncm = 2*k+1;
|
|
|
1056 *nvm = (k+2)**vc/2;
|
|
|
1057 return;
|
|
|
1058 }
|
|
|
1059 if (ev(evs) == EKDCE)
|
|
|
1060 { nterm = (int)(n * alp);
|
|
|
1061 *vc = 1;
|
|
|
1062 *nvm = 1+(int)(2*n/nterm);
|
|
|
1063 *ncm = 2**nvm+1;
|
|
|
1064 return;
|
|
|
1065 }
|
|
|
1066 *nvm = *ncm = *vc = 0;
|
|
|
1067 return;
|
|
|
1068 }
|
|
|
1069
|
|
|
1070 /*
|
|
|
1071 Split x[pi[l..r]] into two equal sized sets.
|
|
|
1072
|
|
|
1073 Let m=(l+r)/2.
|
|
|
1074 At return,
|
|
|
1075 x[pi[l..m]] < x[pi[m+1..r]].
|
|
|
1076 Assuming no ties:
|
|
|
1077 If l+r is odd, the sets have the same size.
|
|
|
1078 If l+r is even, the low set is larger by 1.
|
|
|
1079 If there are ties, all ties go in the low set.
|
|
|
1080 */
|
|
|
1081 int ksmall(l, r, m, x, pi)
|
|
|
1082 int l, r, m, *pi;
|
|
|
1083 double *x;
|
|
|
1084 {
|
|
|
1085 int il, ir, jl, jr;
|
|
|
1086 double t;
|
|
|
1087
|
|
|
1088
|
|
|
1089 while(l<r)
|
|
|
1090 { t = x[pi[m]];
|
|
|
1091
|
|
|
1092 /*
|
|
|
1093 permute the observations so that
|
|
|
1094 x[pi[l..il]] < t <= x[pi[ir..r]].
|
|
|
1095 */
|
|
|
1096 ir = l; il = r;
|
|
|
1097 while (ir<il)
|
|
|
1098 { while ((ir<=r) && (x[pi[ir]] < t)) ir++;
|
|
|
1099 while ((il>=l) && (x[pi[il]]>= t)) il--;
|
|
|
1100 if (ir<il) ISWAP(pi[ir],pi[il]);
|
|
|
1101 }
|
|
|
1102
|
|
|
1103 /*
|
|
|
1104 move = t to the middle:
|
|
|
1105 x[pi[l..il]] < t
|
|
|
1106 x[pi[jl..jr]] = t
|
|
|
1107 x[pi[ir..r]] > t
|
|
|
1108 */
|
|
|
1109 jl = ir; jr = r;
|
|
|
1110 while (ir<jr)
|
|
|
1111 { while ((ir<=r) && (x[pi[ir]]== t)) ir++;
|
|
|
1112 while ((jr>=jl) && (x[pi[jr]] > t)) jr--;
|
|
|
1113 if (ir<jr) ISWAP(pi[ir],pi[jr]);
|
|
|
1114 }
|
|
|
1115
|
|
|
1116 /*
|
|
|
1117 we're done if m is in the middle, jl <= m <= jr.
|
|
|
1118 */
|
|
|
1119 if ((jl<=m) & (jr>=m)) return(jr);
|
|
|
1120
|
|
|
1121 /*
|
|
|
1122 update l or r.
|
|
|
1123 */
|
|
|
1124 if (m>=ir) l = ir;
|
|
|
1125 if (m<=il) r = il;
|
|
|
1126 }
|
|
|
1127 if (l==r) return(l);
|
|
|
1128 LERR(("ksmall failure"));
|
|
|
1129 return(0);
|
|
|
1130 }
|
|
|
1131
|
|
|
1132 int terminal(lf,p,pi,fc,d,m,split_val)
|
|
|
1133 lfit *lf;
|
|
|
1134 int p, d, fc, *m, *pi;
|
|
|
1135 double *split_val;
|
|
|
1136 { int i, k, lo, hi, split_var;
|
|
|
1137 double max, min, score, max_score, t;
|
|
|
1138
|
|
|
1139 /*
|
|
|
1140 if there are fewer than fc points in the cell, this cell
|
|
|
1141 is terminal.
|
|
|
1142 */
|
|
|
1143 lo = lf->evs.lo[p]; hi = lf->evs.hi[p];
|
|
|
1144 if (hi-lo < fc) return(-1);
|
|
|
1145
|
|
|
1146 /* determine the split variable */
|
|
|
1147 max_score = 0.0; split_var = 0;
|
|
|
1148 for (k=0; k<d; k++)
|
|
|
1149 { max = min = datum(&lf->lfd, k, pi[lo]);
|
|
|
1150 for (i=lo+1; i<=hi; i++)
|
|
|
1151 { t = datum(&lf->lfd,k,pi[i]);
|
|
|
1152 if (t<min) min = t;
|
|
|
1153 if (t>max) max = t;
|
|
|
1154 }
|
|
|
1155 score = (max-min) / lf->lfd.sca[k];
|
|
|
1156 if (score > max_score)
|
|
|
1157 { max_score = score;
|
|
|
1158 split_var = k;
|
|
|
1159 }
|
|
|
1160 }
|
|
|
1161 if (max_score==0) /* all points in the cell are equal */
|
|
|
1162 return(-1);
|
|
|
1163
|
|
|
1164 *m = ksmall(lo,hi,(lo+hi)/2, dvari(&lf->lfd,split_var), pi);
|
|
|
1165 *split_val = datum(&lf->lfd, split_var, pi[*m]);
|
|
|
1166
|
|
|
1167 if (*m==hi) /* all observations go lo */
|
|
|
1168 return(-1);
|
|
|
1169 return(split_var);
|
|
|
1170 }
|
|
|
1171
|
|
|
1172 void kdtre_start(des,lf)
|
|
|
1173 design *des;
|
|
|
1174 lfit *lf;
|
|
|
1175 {
|
|
|
1176 int i, j, vc, d, nc, nv, ncm, nvm, k, m, n, p, *pi;
|
|
|
1177 double sv;
|
|
|
1178 d = lf->lfd.d; n = lf->lfd.n; pi = des->ind;
|
|
|
1179 kdtre_guessnv(&lf->evs,&nvm,&ncm,&vc,n,d,nn(&lf->sp));
|
|
|
1180 trchck(lf,nvm,ncm,vc);
|
|
|
1181
|
|
|
1182 nv = 0;
|
|
|
1183 if (ev(&lf->evs) != EKDCE)
|
|
|
1184 { for (i=0; i<vc; i++)
|
|
|
1185 { j = i;
|
|
|
1186 for (k=0; k<d; ++k)
|
|
|
1187 { evptx(&lf->fp,i,k) = lf->evs.fl[d*(j%2)+k];
|
|
|
1188 j >>= 1;
|
|
|
1189 }
|
|
|
1190 }
|
|
|
1191 nv = vc;
|
|
|
1192 for (j=0; j<vc; j++) lf->evs.ce[j] = j;
|
|
|
1193 }
|
|
|
1194
|
|
|
1195 for (i=0; i<n; i++) pi[i] = i;
|
|
|
1196 p = 0; nc = 1;
|
|
|
1197 lf->evs.lo[p] = 0; lf->evs.hi[p] = n-1;
|
|
|
1198 lf->evs.s[p] = -1;
|
|
|
1199 while (p<nc)
|
|
|
1200 { k = terminal(lf,p,pi,nterm,d,&m,&sv);
|
|
|
1201 if (k>=0)
|
|
|
1202 {
|
|
|
1203 if ((ncm<nc+2) | (2*nvm<2*nv+vc))
|
|
|
1204 { WARN(("Insufficient space for full tree"));
|
|
|
1205 lf->evs.nce = nc; lf->fp.nv = nv;
|
|
|
1206 return;
|
|
|
1207 }
|
|
|
1208
|
|
|
1209 /* new lo cell has obsn's lo[p]..m */
|
|
|
1210 lf->evs.lo[nc] = lf->evs.lo[p];
|
|
|
1211 lf->evs.hi[nc] = m;
|
|
|
1212 lf->evs.s[nc] = -1;
|
|
|
1213
|
|
|
1214 /* new hi cell has obsn's m+1..hi[p] */
|
|
|
1215 lf->evs.lo[nc+1] = m+1;
|
|
|
1216 lf->evs.hi[nc+1] = lf->evs.hi[p];
|
|
|
1217 lf->evs.s[nc+1] = -1;
|
|
|
1218
|
|
|
1219 /* cell p is split on variable k, value sv */
|
|
|
1220 lf->evs.s[p] = k;
|
|
|
1221 lf->evs.sv[p] = sv;
|
|
|
1222 lf->evs.lo[p] = nc; lf->evs.hi[p] = nc+1;
|
|
|
1223
|
|
|
1224 nc=nc+2; i = nv;
|
|
|
1225
|
|
|
1226 /* now compute the new vertices. */
|
|
|
1227 if (ev(&lf->evs) != EKDCE)
|
|
|
1228 newcell(&nv,vc,evp(&lf->fp), d, k, sv,
|
|
|
1229 &lf->evs.ce[p*vc], &lf->evs.ce[(nc-2)*vc], &lf->evs.ce[(nc-1)*vc]);
|
|
|
1230
|
|
|
1231 }
|
|
|
1232 else if (ev(&lf->evs)==EKDCE) /* new vertex at cell center */
|
|
|
1233 { sv = 0;
|
|
|
1234 for (i=0; i<d; i++) evptx(&lf->fp,nv,i) = 0;
|
|
|
1235 for (j=lf->evs.lo[p]; j<=lf->evs.hi[p]; j++)
|
|
|
1236 { sv += prwt(&lf->lfd,(int)pi[j]);
|
|
|
1237 for (i=0; i<d; i++)
|
|
|
1238 evptx(&lf->fp,nv,i) += datum(&lf->lfd,i,pi[j])*prwt(&lf->lfd,(int)pi[j]);
|
|
|
1239 }
|
|
|
1240 for (i=0; i<d; i++) evptx(&lf->fp,nv,i) /= sv;
|
|
|
1241 lf->lfd.n = lf->evs.hi[p] - lf->evs.lo[p] + 1;
|
|
|
1242 des->ind = &pi[lf->evs.lo[p]]; /* why? */
|
|
|
1243 PROC_VERTEX(des,lf,nv);
|
|
|
1244 lf->lfd.n = n; des->ind = pi;
|
|
|
1245 nv++;
|
|
|
1246 }
|
|
|
1247 p++;
|
|
|
1248 }
|
|
|
1249
|
|
|
1250 /* We've built the tree. Now do the fitting. */
|
|
|
1251 if (ev(&lf->evs)==EKDTR)
|
|
|
1252 for (i=0; i<nv; i++) PROC_VERTEX(des,lf,i);
|
|
|
1253
|
|
|
1254 lf->evs.nce = nc; lf->fp.nv = nv;
|
|
|
1255 return;
|
|
|
1256 }
|
|
|
1257
|
|
|
1258 void newcell(nv,vc,xev, d, k, split_val, cpar, clef, crig)
|
|
|
1259 double *xev, split_val;
|
|
|
1260 int *cpar, *clef, *crig;
|
|
|
1261 int *nv, vc, d, k;
|
|
|
1262 { int i, ii, j, j2, tk, match;
|
|
|
1263 tk = 1<<k;
|
|
|
1264 for (i=0; i<vc; i++)
|
|
|
1265 { if ((i&tk) == 0)
|
|
|
1266 { for (j=0; j<d; j++) xev[*nv*d+j] = xev[d*cpar[i]+j];
|
|
|
1267 xev[*nv*d+k] = split_val;
|
|
|
1268 match = 0; j = vc; /* no matches in first vc points */
|
|
|
1269 while ((j<*nv) && (!match))
|
|
|
1270 { j2 = 0;
|
|
|
1271 while ((j2<d) && (xev[*nv*d+j2] == xev[j*d+j2])) j2++;
|
|
|
1272 match = (j2==d);
|
|
|
1273 if (!match) j++;
|
|
|
1274 }
|
|
|
1275 ii = i+tk;
|
|
|
1276 clef[i] = cpar[i];
|
|
|
1277 clef[ii]= crig[i] = j;
|
|
|
1278 crig[ii]= cpar[ii];
|
|
|
1279 if (!match) (*nv)++;
|
|
|
1280 }
|
|
|
1281 }
|
|
|
1282 return;
|
|
|
1283 }
|
|
|
1284
|
|
|
1285 extern void hermite2();
|
|
|
1286
|
|
|
1287 double blend(fp,evs,s,x,ll,ur,j,nt,t,what)
|
|
|
1288 fitpt *fp;
|
|
|
1289 evstruc *evs;
|
|
|
1290 double s, *x, *ll, *ur;
|
|
|
1291 int j, nt, *t, what;
|
|
|
1292 {
|
|
|
1293 int *ce, k, k1, m, nc, j0, j1;
|
|
|
1294 double v0, v1, xibar, g0[3], g1[3], gg[4], gp[4], phi[4];
|
|
|
1295 ce = evs->ce;
|
|
|
1296 for (k=0; k<4; k++) /* North South East West */
|
|
|
1297 { k1 = (k>1);
|
|
|
1298 v0 = ll[k1]; v1 = ur[k1];
|
|
|
1299 j0 = ce[j+2*(k==0)+(k==2)];
|
|
|
1300 j1 = ce[j+3-2*(k==1)-(k==3)];
|
|
|
1301 xibar = (k%2==0) ? ur[k<2] : ll[k<2];
|
|
|
1302 m = nt;
|
|
|
1303 while ((m>=0) && ((evs->s[t[m]] != (k<=1)) | (evs->sv[t[m]] != xibar))) m--;
|
|
|
1304 if (m >= 0)
|
|
|
1305 { m = (k%2==1) ? evs->lo[t[m]] : evs->hi[t[m]];
|
|
|
1306 while (evs->s[m] != -1)
|
|
|
1307 m = (x[evs->s[m]] < evs->sv[m]) ? evs->lo[m] : evs->hi[m];
|
|
|
1308 if (v0 < evptx(fp,ce[4*m+2*(k==1)+(k==3)],k1))
|
|
|
1309 { j0 = ce[4*m+2*(k==1)+(k==3)];
|
|
|
1310 v0 = evptx(fp,j0,k1);
|
|
|
1311 }
|
|
|
1312 if (evptx(fp,ce[4*m+3-2*(k==0)-(k==2)],k1) < v1)
|
|
|
1313 { j1 = ce[4*m+3-2*(k==0)-(k==2)];
|
|
|
1314 v1 = evptx(fp,j1,k1);
|
|
|
1315 }
|
|
|
1316 }
|
|
|
1317 nc = exvval(fp,g0,j0,2,what,0);
|
|
|
1318 nc = exvval(fp,g1,j1,2,what,0);
|
|
|
1319 if (nc==1)
|
|
|
1320 gg[k] = linear_interp((x[(k>1)]-v0),v1-v0,g0[0],g1[0]);
|
|
|
1321 else
|
|
|
1322 { hermite2(x[(k>1)]-v0,v1-v0,phi);
|
|
|
1323 gg[k] = phi[0]*g0[0]+phi[1]*g1[0]+(phi[2]*g0[1+k1]+phi[3]*g1[1+k1])*(v1-v0);
|
|
|
1324 gp[k] = phi[0]*g0[2-k1] + phi[1]*g1[2-k1];
|
|
|
1325 }
|
|
|
1326 }
|
|
|
1327 s = -s;
|
|
|
1328 if (nc==1)
|
|
|
1329 for (k=0; k<2; k++)
|
|
|
1330 s += linear_interp(x[k]-ll[k],ur[k]-ll[k],gg[3-2*k],gg[2-2*k]);
|
|
|
1331 else
|
|
|
1332 for (k=0; k<2; k++) /* EW NS */
|
|
|
1333 { hermite2(x[k]-ll[k],ur[k]-ll[k],phi);
|
|
|
1334 s += phi[0]*gg[3-2*k] + phi[1]*gg[2-2*k]
|
|
|
1335 +(phi[2]*gp[3-2*k] + phi[3]*gp[2-2*k]) * (ur[k]-ll[k]);
|
|
|
1336 }
|
|
|
1337 return(s);
|
|
|
1338 }
|
|
|
1339
|
|
|
1340 double kdtre_int(fp,evs,x,what)
|
|
|
1341 fitpt *fp;
|
|
|
1342 evstruc *evs;
|
|
|
1343 double *x;
|
|
|
1344 int what;
|
|
|
1345 {
|
|
|
1346 int *ce, k, vc, t[20], nt, nc, j, d;
|
|
|
1347 double *ll, *ur, ff, vv[64][64];
|
|
|
1348 d = fp->d;
|
|
|
1349 vc = 1<<d;
|
|
|
1350 if (d > 6) { LERR(("d too large in kdint")); return(NOSLN); }
|
|
|
1351
|
|
|
1352 /* descend the tree to find the terminal cell */
|
|
|
1353 nt = 0; t[nt] = 0; k = 0;
|
|
|
1354 while (evs->s[k] != -1)
|
|
|
1355 { nt++;
|
|
|
1356 if (nt>=20) { LERR(("Too many levels in kdint")); return(NOSLN); }
|
|
|
1357 k = t[nt] = (x[evs->s[k]] < evs->sv[k]) ? evs->lo[k] : evs->hi[k];
|
|
|
1358 }
|
|
|
1359
|
|
|
1360 ce = &evs->ce[k*vc];
|
|
|
1361 ll = evpt(fp,ce[0]);
|
|
|
1362 ur = evpt(fp,ce[vc-1]);
|
|
|
1363 nc = 0;
|
|
|
1364 for (j=0; j<vc; j++) nc = exvval(fp,vv[j],(int)ce[j],d,what,0);
|
|
|
1365 ff = rectcell_interp(x,vv,ll,ur,d,nc);
|
|
|
1366
|
|
|
1367 if (d==2) ff = blend(fp,evs,ff,x,ll,ur,k*vc,nt,t,what);
|
|
|
1368 return(ff);
|
|
|
1369 }
|
|
|
1370 /*
|
|
|
1371 * Copyright 1996-2006 Catherine Loader.
|
|
|
1372 */
|
|
|
1373 #include "lfev.h"
|
|
|
1374
|
|
|
1375 /*
|
|
|
1376 * convert eval. structure string to numeric code.
|
|
|
1377 */
|
|
|
1378 #define NETYPE 11
|
|
|
1379 static char *etype[NETYPE]= { "tree", "phull", "data", "grid", "kdtree",
|
|
|
1380 "kdcenter", "cross", "preset", "xbar", "none",
|
|
|
1381 "sphere" };
|
|
|
1382 static int evals[NETYPE]= { ETREE, EPHULL, EDATA, EGRID, EKDTR,
|
|
|
1383 EKDCE, ECROS, EPRES, EXBAR, ENONE, ESPHR };
|
|
|
1384 int lfevstr(char *z)
|
|
|
1385 { return(pmatch(z, etype, evals, NETYPE, -1));
|
|
|
1386 }
|
|
|
1387
|
|
|
1388 void evstruc_init(evs)
|
|
|
1389 evstruc *evs;
|
|
|
1390 { int i;
|
|
|
1391 ev(evs) = ETREE;
|
|
|
1392 mk(evs) = 100;
|
|
|
1393 cut(evs) = 0.8;
|
|
|
1394 for (i=0; i<MXDIM; i++)
|
|
|
1395 { evs->fl[i] = evs->fl[i+MXDIM] = 0.0;
|
|
|
1396 evs->mg[i] = 10;
|
|
|
1397 }
|
|
|
1398 evs->nce = evs->ncm = 0;
|
|
|
1399 }
|
|
|
1400
|
|
|
1401 int evstruc_reqi(nvm,ncm,vc)
|
|
|
1402 int nvm, ncm, vc;
|
|
|
1403 { return(ncm*vc+3*MAX(ncm,nvm));
|
|
|
1404 }
|
|
|
1405
|
|
|
1406 /* al=1: allows dynamic allocation.
|
|
|
1407 * al=0: inhibit. use with caution.
|
|
|
1408 */
|
|
|
1409 void evstruc_alloc(evs,nvm,ncm,vc,al)
|
|
|
1410 evstruc *evs;
|
|
|
1411 int nvm, ncm, vc, al;
|
|
|
1412 { int rw, *k;
|
|
|
1413
|
|
|
1414 if (al)
|
|
|
1415 { rw = evstruc_reqi(nvm,ncm,vc);
|
|
|
1416 if (evs->liw<rw)
|
|
|
1417 { evs->iwk = (int *)calloc(rw,sizeof(int));
|
|
|
1418 if ( evs->iwk == NULL ) {
|
|
|
1419 printf("Problem allocating memory for evs->iwk\n");fflush(stdout);
|
|
|
1420 }
|
|
|
1421 evs->liw = rw;
|
|
|
1422 }
|
|
|
1423 }
|
|
|
1424 k = evs->iwk;
|
|
|
1425 evs->ce = k; k += vc*ncm;
|
|
|
1426 evs->s = k; k += MAX(ncm,nvm);
|
|
|
1427 evs->lo = k; k += MAX(ncm,nvm);
|
|
|
1428 evs->hi = k; k += MAX(ncm,nvm);
|
|
|
1429 }
|
|
|
1430
|
|
|
1431 void guessnv(evs,sp,mdl,n,d,lw,nvc)
|
|
|
1432 evstruc *evs;
|
|
|
1433 smpar *sp;
|
|
|
1434 module *mdl;
|
|
|
1435 int n, d, *lw, *nvc;
|
|
|
1436 { int i, nvm, ncm, vc;
|
|
|
1437
|
|
|
1438 npar(sp) = calcp(sp,d);
|
|
|
1439 switch(ev(evs))
|
|
|
1440 { case ETREE:
|
|
|
1441 atree_guessnv(evs,&nvm,&ncm,&vc,d,nn(sp));
|
|
|
1442 break;
|
|
|
1443 case EPHULL:
|
|
|
1444 nvm = ncm = mk(evs)*d;
|
|
|
1445 vc = d+1;
|
|
|
1446 break;
|
|
|
1447 case EDATA:
|
|
|
1448 case ECROS:
|
|
|
1449 nvm = n;
|
|
|
1450 ncm = vc = 0;
|
|
|
1451 break;
|
|
|
1452 case EKDTR:
|
|
|
1453 case EKDCE:
|
|
|
1454 kdtre_guessnv(evs,&nvm,&ncm,&vc,n,d,nn(sp));
|
|
|
1455 break;
|
|
|
1456 case EGRID:
|
|
|
1457 nvm = 1;
|
|
|
1458 for (i=0; i<d; i++) nvm *= evs->mg[i];
|
|
|
1459 ncm = 0;
|
|
|
1460 vc = 1<<d;
|
|
|
1461 break;
|
|
|
1462 case EXBAR:
|
|
|
1463 case ENONE:
|
|
|
1464 nvm = 1;
|
|
|
1465 ncm = vc = 0;
|
|
|
1466 break;
|
|
|
1467 case EPRES:
|
|
|
1468 nvm = evs->mg[0];
|
|
|
1469 ncm = vc = 0;
|
|
|
1470 break;
|
|
|
1471 default:
|
|
|
1472 LERR(("guessnv: I don't know this evaluation structure."));
|
|
|
1473 nvm = ncm = vc = 0;
|
|
|
1474 }
|
|
|
1475
|
|
|
1476 lw[0] = des_reqd(n,npar(sp));
|
|
|
1477 lw[1] = lfit_reqd(d,nvm,ncm,mdl);
|
|
|
1478 lw[2] = evstruc_reqi(nvm,ncm,vc);
|
|
|
1479 lw[6] = des_reqi(n,npar(sp));
|
|
|
1480 lw[3] = pc_reqd(d,npar(sp));
|
|
|
1481 lw[4] = mdl->keepv;
|
|
|
1482 lw[5] = mdl->keepc;
|
|
|
1483
|
|
|
1484 if (nvc==NULL) return;
|
|
|
1485 nvc[0] = nvm;
|
|
|
1486 nvc[1] = ncm;
|
|
|
1487 nvc[2] = vc;
|
|
|
1488 nvc[3] = nvc[4] = 0;
|
|
|
1489 }
|
|
|
1490
|
|
|
1491 /*
|
|
|
1492 * trchck checks the working space on the lfit structure
|
|
|
1493 * has space for nvm vertices and ncm cells.
|
|
|
1494 */
|
|
|
1495 void lfit_alloc(lf)
|
|
|
1496 lfit *lf;
|
|
|
1497 { lf->fp.lwk = lf->fp.lev = lf->evs.liw = lf->pc.lwk = 0;
|
|
|
1498 lf->lf_init_id = LF_INIT_ID;
|
|
|
1499 }
|
|
|
1500 int lfit_reqd(d,nvm,ncm,mdl)
|
|
|
1501 int d, nvm, ncm;
|
|
|
1502 module *mdl;
|
|
|
1503 { int z;
|
|
|
1504 z = mdl->keepv;
|
|
|
1505 return(nvm*z+ncm);
|
|
|
1506 }
|
|
|
1507
|
|
|
1508 void trchck(lf,nvm,ncm,vc)
|
|
|
1509 lfit *lf;
|
|
|
1510 int nvm, ncm, vc;
|
|
|
1511 { int rw, d, *k;
|
|
|
1512 double *z;
|
|
|
1513
|
|
|
1514 if (lf->lf_init_id != LF_INIT_ID) lfit_alloc(lf);
|
|
|
1515
|
|
|
1516 lf->fp.nvm = nvm; lf->evs.ncm = ncm;
|
|
|
1517 d = lf->lfd.d;
|
|
|
1518
|
|
|
1519 if (lf->fp.lev < d*nvm)
|
|
|
1520 { lf->fp.xev = (double *)calloc(d*nvm,sizeof(double));
|
|
|
1521 if ( lf->fp.xev == NULL ) {
|
|
|
1522 printf("Problem allocating memory for lf->fp.xev\n");fflush(stdout);
|
|
|
1523 }
|
|
|
1524 lf->fp.lev = d*nvm;
|
|
|
1525 }
|
|
|
1526
|
|
|
1527 rw = lfit_reqd(d,nvm,ncm,&lf->mdl);
|
|
|
1528 if (lf->fp.lwk < rw)
|
|
|
1529 {
|
|
|
1530 lf->fp.coef = (double *)calloc(rw,sizeof(double));
|
|
|
1531 if ( lf->fp.coef == NULL ) {
|
|
|
1532 printf("Problem allocating memory for lf->fp.coef\n");fflush(stdout);
|
|
|
1533 }
|
|
|
1534 lf->fp.lwk = rw;
|
|
|
1535 }
|
|
|
1536 z = lf->fp.wk = lf->fp.coef;
|
|
|
1537
|
|
|
1538 lf->fp.h = NULL;
|
|
|
1539 if (!lf->mdl.isset) mut_printf("module not set.\n");
|
|
|
1540 else
|
|
|
1541 { if (lf->mdl.alloc!=NULL) lf->mdl.alloc(lf);
|
|
|
1542 z += KEEPV(lf) * nvm;
|
|
|
1543 }
|
|
|
1544 lf->evs.sv = z; z += ncm;
|
|
|
1545
|
|
|
1546 evstruc_alloc(&lf->evs,nvm,ncm,vc,1);
|
|
|
1547 }
|
|
|
1548
|
|
|
1549 void data_guessnv(nvm,ncm,vc,n)
|
|
|
1550 int *nvm, *ncm, *vc, n;
|
|
|
1551 { *nvm = n;
|
|
|
1552 *ncm = *vc = 0;
|
|
|
1553 }
|
|
|
1554
|
|
|
1555 void dataf(des,lf)
|
|
|
1556 design *des;
|
|
|
1557 lfit *lf;
|
|
|
1558 {
|
|
|
1559 int d, i, j, ncm, nv, vc;
|
|
|
1560
|
|
|
1561 d = lf->lfd.d;
|
|
|
1562 data_guessnv(&nv,&ncm,&vc,lf->lfd.n);
|
|
|
1563 trchck(lf,nv,ncm,vc);
|
|
|
1564
|
|
|
1565 for (i=0; i<nv; i++)
|
|
|
1566 for (j=0; j<d; j++) evptx(&lf->fp,i,j) = datum(&lf->lfd,j,i);
|
|
|
1567 for (i=0; i<nv; i++)
|
|
|
1568 {
|
|
|
1569 PROC_VERTEX(des,lf,i);
|
|
|
1570 lf->evs.s[i] = 0;
|
|
|
1571 }
|
|
|
1572 lf->fp.nv = lf->fp.nvm = nv; lf->evs.nce = 0;
|
|
|
1573 }
|
|
|
1574
|
|
|
1575 void xbar_guessnv(nvm,ncm,vc)
|
|
|
1576 int *nvm, *ncm, *vc;
|
|
|
1577 { *nvm = 1;
|
|
|
1578 *ncm = *vc = 0;
|
|
|
1579 return;
|
|
|
1580 }
|
|
|
1581
|
|
|
1582 void xbarf(des,lf)
|
|
|
1583 design *des;
|
|
|
1584 lfit *lf;
|
|
|
1585 { int i, d, nvm, ncm, vc;
|
|
|
1586 d = lf->lfd.d;
|
|
|
1587 xbar_guessnv(&nvm,&ncm,&vc);
|
|
|
1588 trchck(lf,1,0,0);
|
|
|
1589 for (i=0; i<d; i++) evptx(&lf->fp,0,i) = lf->pc.xbar[i];
|
|
|
1590 PROC_VERTEX(des,lf,0);
|
|
|
1591 lf->evs.s[0] = 0;
|
|
|
1592 lf->fp.nv = 1; lf->evs.nce = 0;
|
|
|
1593 }
|
|
|
1594
|
|
|
1595 void preset(des,lf)
|
|
|
1596 design *des;
|
|
|
1597 lfit *lf;
|
|
|
1598 { int i, nv;
|
|
|
1599
|
|
|
1600 nv = lf->fp.nvm;
|
|
|
1601 trchck(lf,nv,0,0);
|
|
|
1602 for (i=0; i<nv; i++)
|
|
|
1603 {
|
|
|
1604 PROC_VERTEX(des,lf,i);
|
|
|
1605 lf->evs.s[i] = 0;
|
|
|
1606 }
|
|
|
1607 lf->fp.nv = nv; lf->evs.nce = 0;
|
|
|
1608 }
|
|
|
1609
|
|
|
1610 void crossf(des,lf)
|
|
|
1611 design *des;
|
|
|
1612 lfit *lf;
|
|
|
1613 { int d, i, j, n, nv, ncm, vc;
|
|
|
1614 double w;
|
|
|
1615
|
|
|
1616 n = lf->lfd.n; d = lf->lfd.d;
|
|
|
1617 data_guessnv(&nv,&ncm,&vc,n);
|
|
|
1618 trchck(lf,nv,ncm,vc);
|
|
|
1619
|
|
|
1620 if (lf->lfd.w==NULL) LERR(("crossf() needs prior weights"));
|
|
|
1621 for (i=0; i<n; i++)
|
|
|
1622 for (j=0; j<d; j++) evptx(&lf->fp,i,j) = datum(&lf->lfd,j,i);
|
|
|
1623 for (i=0; i<n; i++)
|
|
|
1624 { lf->evs.s[i] = 0;
|
|
|
1625 w = prwt(&lf->lfd,i);
|
|
|
1626 lf->lfd.w[i] = 0;
|
|
|
1627 PROC_VERTEX(des,lf,i);
|
|
|
1628 lf->lfd.w[i] = w;
|
|
|
1629 }
|
|
|
1630 lf->fp.nv = n; lf->evs.nce = 0;
|
|
|
1631 }
|
|
|
1632
|
|
|
1633 void gridf(des,lf)
|
|
|
1634 design *des;
|
|
|
1635 lfit *lf;
|
|
|
1636 { int d, i, j, nv, u0, u1, z;
|
|
|
1637 nv = 1; d = lf->lfd.d;
|
|
|
1638 for (i=0; i<d; i++)
|
|
|
1639 { if (lf->evs.mg[i]==0)
|
|
|
1640 lf->evs.mg[i] = 2+(int)((lf->evs.fl[i+d]-lf->evs.fl[i])/(lf->lfd.sca[i]*cut(&lf->evs)));
|
|
|
1641 nv *= lf->evs.mg[i];
|
|
|
1642 }
|
|
|
1643 trchck(lf,nv,0,1<<d);
|
|
|
1644 for (i=0; i<nv; i++)
|
|
|
1645 { z = i;
|
|
|
1646 for (j=0; j<d; j++)
|
|
|
1647 { u0 = z%lf->evs.mg[j];
|
|
|
1648 u1 = lf->evs.mg[j]-1-u0;
|
|
|
1649 evptx(&lf->fp,i,j) = (lf->evs.mg[j]==1) ? lf->evs.fl[j] :
|
|
|
1650 (u1*lf->evs.fl[j]+u0*lf->evs.fl[j+d])/(lf->evs.mg[j]-1);
|
|
|
1651 z = z/lf->evs.mg[j];
|
|
|
1652 }
|
|
|
1653 lf->evs.s[i] = 0;
|
|
|
1654 PROC_VERTEX(des,lf,i);
|
|
|
1655 }
|
|
|
1656 lf->fp.nv = nv; lf->evs.nce = 0;
|
|
|
1657 }
|
|
|
1658
|
|
|
1659 int findpt(fp,evs,i0,i1)
|
|
|
1660 fitpt *fp;
|
|
|
1661 evstruc *evs;
|
|
|
1662 int i0, i1;
|
|
|
1663 { int i;
|
|
|
1664 if (i0>i1) ISWAP(i0,i1);
|
|
|
1665 for (i=i1+1; i<fp->nv; i++)
|
|
|
1666 if ((evs->lo[i]==i0) && (evs->hi[i]==i1)) return(i);
|
|
|
1667 return(-1);
|
|
|
1668 }
|
|
|
1669
|
|
|
1670 /*
|
|
|
1671 add a new vertex at the midpoint of (x[i0],x[i1]).
|
|
|
1672 return the vertex number.
|
|
|
1673 */
|
|
|
1674 int newsplit(des,lf,i0,i1,pv)
|
|
|
1675 design *des;
|
|
|
1676 lfit *lf;
|
|
|
1677 int i0, i1, pv;
|
|
|
1678 { int i, nv;
|
|
|
1679
|
|
|
1680 i = findpt(&lf->fp,&lf->evs,i0,i1);
|
|
|
1681 if (i>=0) return(i);
|
|
|
1682
|
|
|
1683 if (i0>i1) ISWAP(i0,i1);
|
|
|
1684 nv = lf->fp.nv;
|
|
|
1685
|
|
|
1686 /* the point is new. Now check we have space for the new point. */
|
|
|
1687 if (nv>=lf->fp.nvm)
|
|
|
1688 {
|
|
|
1689 LERR(("newsplit: out of vertex space"));
|
|
|
1690 return(-1);
|
|
|
1691 }
|
|
|
1692
|
|
|
1693 /* compute the new point, and evaluate the fit */
|
|
|
1694 lf->evs.lo[nv] = i0;
|
|
|
1695 lf->evs.hi[nv] = i1;
|
|
|
1696 for (i=0; i<lf->fp.d; i++)
|
|
|
1697 evptx(&lf->fp,nv,i) = (evptx(&lf->fp,i0,i)+evptx(&lf->fp,i1,i))/2;
|
|
|
1698 if (pv) /* pseudo vertex */
|
|
|
1699 { lf->fp.h[nv] = (lf->fp.h[i0]+lf->fp.h[i1])/2;
|
|
|
1700 lf->evs.s[nv] = 1; /* pseudo-vertex */
|
|
|
1701 }
|
|
|
1702 else /* real vertex */
|
|
|
1703 {
|
|
|
1704 PROC_VERTEX(des,lf,nv);
|
|
|
1705 lf->evs.s[nv] = 0;
|
|
|
1706 }
|
|
|
1707 lf->fp.nv++;
|
|
|
1708
|
|
|
1709 return(nv);
|
|
|
1710 }
|
|
|
1711 /*
|
|
|
1712 * Copyright 1996-2006 Catherine Loader.
|
|
|
1713 */
|
|
|
1714 /*
|
|
|
1715 * Functions for constructing the fit and
|
|
|
1716 * interpolating on the circle/sphere. d=2 only.
|
|
|
1717 */
|
|
|
1718
|
|
|
1719 #include "lfev.h"
|
|
|
1720
|
|
|
1721 /*
|
|
|
1722 * Guess the number of fitting points.
|
|
|
1723 */
|
|
|
1724 void sphere_guessnv(nvm,ncm,vc,mg)
|
|
|
1725 int *nvm, *ncm, *vc, *mg;
|
|
|
1726 { *nvm = mg[1]*(mg[0]+1);
|
|
|
1727 *ncm = 0;
|
|
|
1728 *vc = 0;
|
|
|
1729 }
|
|
|
1730
|
|
|
1731 void sphere_start(des,lf)
|
|
|
1732 design *des;
|
|
|
1733 lfit *lf;
|
|
|
1734 { int d, i, j, ct, nv, ncm, vc, *mg;
|
|
|
1735 double rmin, rmax, *orig, r, th, c, s;
|
|
|
1736
|
|
|
1737 mg = mg(&lf->evs);
|
|
|
1738 sphere_guessnv(&nv,&ncm,&vc,mg);
|
|
|
1739 trchck(lf,nv,0,0);
|
|
|
1740 d = lf->lfd.d;
|
|
|
1741
|
|
|
1742 rmin = lf->evs.fl[0];
|
|
|
1743 rmax = lf->evs.fl[1];
|
|
|
1744 orig = &lf->evs.fl[2];
|
|
|
1745 rmin = 0; rmax = 1; orig[0] = orig[1] = 0.0;
|
|
|
1746
|
|
|
1747 ct = 0;
|
|
|
1748 for (i=0; i<mg[1]; i++)
|
|
|
1749 { th = 2*PI*i/mg[1];
|
|
|
1750 c = cos(th);
|
|
|
1751 s = sin(th);
|
|
|
1752 for (j=0; j<=mg[0]; j++)
|
|
|
1753 { r = rmin + (rmax-rmin)*j/mg[0];
|
|
|
1754 evptx(&lf->fp,ct,0) = orig[0] + r*c;
|
|
|
1755 evptx(&lf->fp,ct,1) = orig[1] + r*s;
|
|
|
1756 PROC_VERTEX(des,lf,ct);
|
|
|
1757 ct++;
|
|
|
1758 }
|
|
|
1759 }
|
|
|
1760 lf->fp.nv = ct;
|
|
|
1761 lf->evs.nce = 0;
|
|
|
1762 }
|
|
|
1763
|
|
|
1764 double sphere_int(lf,x,what)
|
|
|
1765 lfit *lf;
|
|
|
1766 double *x;
|
|
|
1767 int what;
|
|
|
1768 { double rmin, rmax, *orig, dx, dy, r, th, th0, th1;
|
|
|
1769 double v[64][64], c0, c1, s0, s1, r0, r1, d0, d1;
|
|
|
1770 double ll[2], ur[2], xx[2];
|
|
|
1771 int i0, j0, i1, j1, *mg, nc, ce[4];
|
|
|
1772
|
|
|
1773 rmin = lf->evs.fl[0];
|
|
|
1774 rmax = lf->evs.fl[1];
|
|
|
1775 orig = &lf->evs.fl[2];
|
|
|
1776 rmin = 0; rmax = 1; orig[0] = orig[1] = 0.0;
|
|
|
1777 mg = mg(&lf->evs);
|
|
|
1778
|
|
|
1779 dx = x[0] - orig[0];
|
|
|
1780 dy = x[1] - orig[1];
|
|
|
1781 r = sqrt(dx*dx+dy*dy);
|
|
|
1782 th = atan2(dy,dx); /* between -pi and pi */
|
|
|
1783
|
|
|
1784 i0 = (int)floor(mg[1]*th/(2*PI)) % mg[1];
|
|
|
1785 j0 = (int)(mg[0]*(r-rmin)/(rmax-rmin));
|
|
|
1786
|
|
|
1787 i1 = (i0+1) % mg[1];
|
|
|
1788 j1 = j0+1; if (j1>mg[0]) { j0 = mg[0]-1; j1 = mg[0]; }
|
|
|
1789
|
|
|
1790 ce[0] = i0*(mg[0]+1)+j0;
|
|
|
1791 ce[1] = i0*(mg[0]+1)+j1;
|
|
|
1792 ce[2] = i1*(mg[0]+1)+j0;
|
|
|
1793 ce[3] = i1*(mg[0]+1)+j1;
|
|
|
1794 nc = exvval(&lf->fp,v[0],ce[0],2,what,1);
|
|
|
1795 nc = exvval(&lf->fp,v[1],ce[1],2,what,1);
|
|
|
1796 nc = exvval(&lf->fp,v[2],ce[2],2,what,1);
|
|
|
1797 nc = exvval(&lf->fp,v[3],ce[3],2,what,1);
|
|
|
1798
|
|
|
1799 th0 = 2*PI*i0/mg[1]; c0 = cos(th0); s0 = sin(th0);
|
|
|
1800 th1 = 2*PI*i1/mg[1]; c1 = cos(th1); s1 = sin(th1);
|
|
|
1801 r0 = rmin + j0*(rmax-rmin)/mg[0];
|
|
|
1802 r1 = rmin + j1*(rmax-rmin)/mg[0];
|
|
|
1803
|
|
|
1804 d0 = c0*v[0][1] + s0*v[0][2];
|
|
|
1805 d1 = r0*(c0*v[0][2]-s0*v[0][1]);
|
|
|
1806 v[0][1] = d0; v[0][2] = d1;
|
|
|
1807
|
|
|
1808 d0 = c0*v[1][1] + s0*v[1][2];
|
|
|
1809 d1 = r1*(c0*v[1][2]-s0*v[1][1]);
|
|
|
1810 v[1][1] = d0; v[1][2] = d1;
|
|
|
1811
|
|
|
1812 d0 = c1*v[2][1] + s1*v[2][2];
|
|
|
1813 d1 = r0*(c1*v[2][2]-s1*v[2][1]);
|
|
|
1814 v[2][1] = d0; v[2][2] = d1;
|
|
|
1815
|
|
|
1816 d0 = c1*v[3][1] + s1*v[3][2];
|
|
|
1817 d1 = r1*(c1*v[3][2]-s1*v[3][1]);
|
|
|
1818 v[3][1] = d0; v[3][2] = d1;
|
|
|
1819
|
|
|
1820 xx[0] = r; xx[1] = th;
|
|
|
1821 ll[0] = r0; ll[1] = th0;
|
|
|
1822 ur[0] = r1; ur[1] = th1;
|
|
|
1823 return(rectcell_interp(xx,v,ll,ur,2,nc));
|
|
|
1824 }
|
|
|
1825 /*
|
|
|
1826 * Copyright 1996-2006 Catherine Loader.
|
|
|
1827 */
|
|
|
1828 #include "lfev.h"
|
|
|
1829
|
|
|
1830 void solve(A,b,d) /* this is crude! A organized by column. */
|
|
|
1831 double *A, *b;
|
|
|
1832 int d;
|
|
|
1833 { int i, j, k;
|
|
|
1834 double piv;
|
|
|
1835 for (i=0; i<d; i++)
|
|
|
1836 { piv = A[(d+1)*i];
|
|
|
1837 for (j=i; j<d; j++) A[j*d+i] /= piv;
|
|
|
1838 b[i] /= piv;
|
|
|
1839 for (j=0; j<d; j++) if (j != i)
|
|
|
1840 { piv = A[i*d+j];
|
|
|
1841 A[i*d+j] = 0;
|
|
|
1842 for (k=i+1; k<d; k++)
|
|
|
1843 A[k*d+j] -= piv*A[k*d+i];
|
|
|
1844 b[j] -= piv*b[i];
|
|
|
1845 }
|
|
|
1846 }
|
|
|
1847 }
|
|
|
1848
|
|
|
1849 void triang_guessnv(nvm,ncm,vc,d,mk)
|
|
|
1850 int *nvm, *ncm, *vc, d, mk;
|
|
|
1851 { *nvm = *ncm = mk*d;
|
|
|
1852 *vc = d+1;
|
|
|
1853 return;
|
|
|
1854 }
|
|
|
1855
|
|
|
1856 int triang_split(lf,ce,le)
|
|
|
1857 lfit *lf;
|
|
|
1858 double *le;
|
|
|
1859 int *ce;
|
|
|
1860 { int d, i, j, k, nts, vc;
|
|
|
1861 double di, dfx[MXDIM];
|
|
|
1862 nts = 0; d = lf->fp.d; vc = d+1;
|
|
|
1863 for (i=0; i<d; i++)
|
|
|
1864 for (j=i+1; j<=d; j++)
|
|
|
1865 { for (k=0; k<d; k++)
|
|
|
1866 dfx[k] = evptx(&lf->fp,ce[i],k)-evptx(&lf->fp,ce[j],k);
|
|
|
1867 di = rho(dfx,lf->lfd.sca,d,KSPH,NULL);
|
|
|
1868 le[i*vc+j] = le[j*vc+i] = di/MIN(lf->fp.h[ce[i]],lf->fp.h[ce[j]]);
|
|
|
1869 nts = nts || le[i*vc+j]>cut(&lf->evs);
|
|
|
1870 }
|
|
|
1871 return(nts);
|
|
|
1872 }
|
|
|
1873
|
|
|
1874 void resort(pv,xev,dig)
|
|
|
1875 double *xev;
|
|
|
1876 int *pv, *dig;
|
|
|
1877 { double d0, d1, d2;
|
|
|
1878 int i;
|
|
|
1879 d0 = d1 = d2 = 0;
|
|
|
1880 for (i=0; i<3; i++)
|
|
|
1881 { d0 += (xev[3*pv[11]+i]-xev[3*pv[1]+i])*(xev[3*pv[11]+i]-xev[3*pv[1]+i]);
|
|
|
1882 d1 += (xev[3*pv[ 7]+i]-xev[3*pv[2]+i])*(xev[3*pv[ 7]+i]-xev[3*pv[2]+i]);
|
|
|
1883 d2 += (xev[3*pv[ 6]+i]-xev[3*pv[3]+i])*(xev[3*pv[ 6]+i]-xev[3*pv[3]+i]);
|
|
|
1884 }
|
|
|
1885 if ((d0<=d1) & (d0<=d2))
|
|
|
1886 { dig[0] = pv[1]; dig[1] = pv[11];
|
|
|
1887 dig[2] = pv[2]; dig[3] = pv[7];
|
|
|
1888 dig[4] = pv[3]; dig[5] = pv[6];
|
|
|
1889 }
|
|
|
1890 else if (d1<=d2)
|
|
|
1891 { dig[0] = pv[2]; dig[1] = pv[7];
|
|
|
1892 dig[2] = pv[1]; dig[3] = pv[11];
|
|
|
1893 dig[4] = pv[3]; dig[5] = pv[6];
|
|
|
1894 }
|
|
|
1895 else
|
|
|
1896 { dig[0] = pv[3]; dig[1] = pv[6];
|
|
|
1897 dig[2] = pv[2]; dig[3] = pv[7];
|
|
|
1898 dig[4] = pv[1]; dig[5] = pv[11];
|
|
|
1899 }
|
|
|
1900 }
|
|
|
1901
|
|
|
1902 void triang_grow(des,lf,ce,ct,term)
|
|
|
1903 design *des;
|
|
|
1904 lfit *lf;
|
|
|
1905 int *ce, *ct, *term;
|
|
|
1906 { double le[(1+MXDIM)*(1+MXDIM)], ml;
|
|
|
1907 int d, i, j, im, jm, vc, pv[(1+MXDIM)*(1+MXDIM)], dig[6];
|
|
|
1908 int nce[1+MXDIM];
|
|
|
1909 if (lf_error) return;
|
|
|
1910 d = lf->fp.d; vc = d+1;
|
|
|
1911 if (!triang_split(lf,ce,le))
|
|
|
1912 { if (ct != NULL)
|
|
|
1913 { for (i=0; i<vc; i++) term[*ct*vc+i] = ce[i];
|
|
|
1914 (*ct)++;
|
|
|
1915 }
|
|
|
1916 return;
|
|
|
1917 }
|
|
|
1918 if (d>3)
|
|
|
1919 { ml = 0;
|
|
|
1920 for (i=0; i<d; i++)
|
|
|
1921 for (j=i+1; j<vc; j++)
|
|
|
1922 if (le[i*vc+j]>ml) { ml = le[i*vc+j]; im = i; jm = j; }
|
|
|
1923 pv[0] = newsplit(des,lf,(int)ce[im],(int)ce[jm],0);
|
|
|
1924 for (i=0; i<vc; i++) nce[i] = ce[i];
|
|
|
1925 nce[im] = pv[0]; triang_grow(des,lf,nce,ct,term); nce[im] = ce[im];
|
|
|
1926 nce[jm] = pv[0]; triang_grow(des,lf,nce,ct,term);
|
|
|
1927 return;
|
|
|
1928 }
|
|
|
1929
|
|
|
1930 for (i=0; i<d; i++)
|
|
|
1931 for (j=i+1; j<=d; j++)
|
|
|
1932 pv[i*vc+j] = pv[j*vc+i]
|
|
|
1933 = newsplit(des,lf,(int)ce[i],(int)ce[j],le[i*vc+j]<=cut(&lf->evs));
|
|
|
1934 for (i=0; i<=d; i++) /* corners */
|
|
|
1935 { for (j=0; j<=d; j++) nce[j] = (j==i) ? ce[i] : pv[i*vc+j];
|
|
|
1936 triang_grow(des,lf,nce,ct,term);
|
|
|
1937 }
|
|
|
1938
|
|
|
1939 if (d==2) /* center for d=2 */
|
|
|
1940 { nce[0] = pv[5]; nce[1] = pv[2]; nce[2] = pv[1];
|
|
|
1941 triang_grow(des,lf,nce,ct,term);
|
|
|
1942 }
|
|
|
1943 if (d==3) /* center for d=3 */
|
|
|
1944 { resort(pv,evp(&lf->fp),dig);
|
|
|
1945 nce[0] = dig[0]; nce[1] = dig[1];
|
|
|
1946 nce[2] = dig[2]; nce[3] = dig[4]; triang_grow(des,lf,nce,ct,term);
|
|
|
1947 nce[2] = dig[5]; nce[3] = dig[3]; triang_grow(des,lf,nce,ct,term);
|
|
|
1948 nce[2] = dig[2]; nce[3] = dig[5]; triang_grow(des,lf,nce,ct,term);
|
|
|
1949 nce[2] = dig[4]; nce[3] = dig[3]; triang_grow(des,lf,nce,ct,term);
|
|
|
1950 }
|
|
|
1951 }
|
|
|
1952
|
|
|
1953 void triang_descend(tr,xa,ce)
|
|
|
1954 lfit *tr;
|
|
|
1955 double *xa;
|
|
|
1956 int *ce;
|
|
|
1957 { double le[(1+MXDIM)*(1+MXDIM)], ml;
|
|
|
1958 int d, vc, i, j, im, jm, pv[(1+MXDIM)*(1+MXDIM)];
|
|
|
1959 design *des;
|
|
|
1960 des = NULL;
|
|
|
1961 if (!triang_split(tr,ce,le)) return;
|
|
|
1962 d = tr->fp.d; vc = d+1;
|
|
|
1963
|
|
|
1964 if (d>3) /* split longest edge */
|
|
|
1965 { ml = 0;
|
|
|
1966 for (i=0; i<d; i++)
|
|
|
1967 for (j=i+1; j<vc; j++)
|
|
|
1968 if (le[i*vc+j]>ml) { ml = le[i*vc+j]; im = i; jm = j; }
|
|
|
1969 pv[0] = newsplit(des,tr,(int)ce[im],(int)ce[jm],0);
|
|
|
1970 if (xa[im]>xa[jm])
|
|
|
1971 { xa[im] -= xa[jm]; xa[jm] *= 2; ce[jm] = pv[0]; }
|
|
|
1972 else
|
|
|
1973 { xa[jm] -= xa[im]; xa[im] *= 2; ce[im] = pv[0]; }
|
|
|
1974 triang_descend(tr,xa,ce);
|
|
|
1975 return;
|
|
|
1976 }
|
|
|
1977
|
|
|
1978 for (i=0; i<d; i++)
|
|
|
1979 for (j=i+1; j<=d; j++)
|
|
|
1980 pv[i*vc+j] = pv[j*vc+i]
|
|
|
1981 = newsplit(des,tr,(int)ce[i],(int)ce[j],le[i*d+j]<=cut(&tr->evs));
|
|
|
1982 for (i=0; i<=d; i++) if (xa[i]>=0.5) /* in corner */
|
|
|
1983 { for (j=0; j<=d; j++)
|
|
|
1984 { if (i!=j) ce[j] = pv[i*vc+j];
|
|
|
1985 xa[j] = 2*xa[j];
|
|
|
1986 }
|
|
|
1987 xa[i] -= 1;
|
|
|
1988 triang_descend(tr,xa,ce);
|
|
|
1989 return;
|
|
|
1990 }
|
|
|
1991 if (d==1) { LERR(("weights sum to < 1")); }
|
|
|
1992 if (d==2) /* center */
|
|
|
1993 { ce[0] = pv[5]; xa[0] = 1-2*xa[0];
|
|
|
1994 ce[1] = pv[2]; xa[1] = 1-2*xa[1];
|
|
|
1995 ce[2] = pv[1]; xa[2] = 1-2*xa[2];
|
|
|
1996 triang_descend(tr,xa,ce);
|
|
|
1997 }
|
|
|
1998 if (d==3) /* center */
|
|
|
1999 { double z; int dig[6];
|
|
|
2000 resort(pv,evp(&tr->fp),dig);
|
|
|
2001 ce[0] = dig[0]; ce[1] = dig[1];
|
|
|
2002 xa[0] *= 2; xa[1] *= 2; xa[2] *= 2; xa[3] *= 2;
|
|
|
2003 if (xa[0]+xa[2]>=1)
|
|
|
2004 { if (xa[0]+xa[3]>=1)
|
|
|
2005 { ce[2] = dig[2]; ce[3] = dig[4];
|
|
|
2006 z = xa[0];
|
|
|
2007 xa[3] += z-1; xa[2] += z-1; xa[0] = xa[1]; xa[1] = 1-z;
|
|
|
2008 }
|
|
|
2009 else
|
|
|
2010 { ce[2] = dig[2]; ce[3] = dig[5];
|
|
|
2011 z = xa[3]; xa[3] = xa[1]+xa[2]-1; xa[1] = z;
|
|
|
2012 z = xa[2]; xa[2] += xa[0]-1; xa[0] = 1-z;
|
|
|
2013 } }
|
|
|
2014 else
|
|
|
2015 { if (xa[1]+xa[2]>=1)
|
|
|
2016 { ce[2] = dig[5]; ce[3] = dig[3];
|
|
|
2017 xa[1] = 1-xa[1]; xa[2] -= xa[1]; xa[3] -= xa[1];
|
|
|
2018 }
|
|
|
2019 else
|
|
|
2020 { ce[2] = dig[4]; ce[3] = dig[3];
|
|
|
2021 z = xa[3]; xa[3] += xa[1]-1; xa[1] = xa[2];
|
|
|
2022 xa[2] = z+xa[0]-1; xa[0] = 1-z;
|
|
|
2023 } }
|
|
|
2024 triang_descend(tr,xa,ce);
|
|
|
2025 } }
|
|
|
2026
|
|
|
2027 void covrofdata(lfd,V,mn) /* covar of data; mean in mn */
|
|
|
2028 lfdata *lfd;
|
|
|
2029 double *V, *mn;
|
|
|
2030 { int d, i, j, k;
|
|
|
2031 double s;
|
|
|
2032 s = 0; d = lfd->d;
|
|
|
2033 for (i=0; i<d*d; i++) V[i] = 0;
|
|
|
2034 for (i=0; i<lfd->n; i++)
|
|
|
2035 { s += prwt(lfd,i);
|
|
|
2036 for (j=0; j<d; j++)
|
|
|
2037 for (k=0; k<d; k++)
|
|
|
2038 V[j*d+k] += prwt(lfd,i)*(datum(lfd,j,i)-mn[j])*(datum(lfd,k,i)-mn[k]);
|
|
|
2039 }
|
|
|
2040 for (i=0; i<d*d; i++) V[i] /= s;
|
|
|
2041 }
|
|
|
2042
|
|
|
2043 int intri(x,w,xev,xa,d) /* is x in triangle bounded by xd[0..d-1]? */
|
|
|
2044 double *x, *xev, *xa;
|
|
|
2045 int *w, d;
|
|
|
2046 { int i, j;
|
|
|
2047 double eps, *r, xd[MXDIM*MXDIM];
|
|
|
2048 eps = 1.0e-10;
|
|
|
2049 r = &xev[w[d]*d];
|
|
|
2050 for (i=0; i<d; i++)
|
|
|
2051 { xa[i] = x[i]-r[i];
|
|
|
2052 for (j=0; j<d; j++) xd[i*d+j] = xev[w[i]*d+j]-r[j];
|
|
|
2053 }
|
|
|
2054 solve(xd,xa,d);
|
|
|
2055 xa[d] = 1.0;
|
|
|
2056 for (i=0; i<d; i++) xa[d] -= xa[i];
|
|
|
2057 for (i=0; i<=d; i++) if ((xa[i]<-eps) | (xa[i]>1+eps)) return(0);
|
|
|
2058 return(1);
|
|
|
2059 }
|
|
|
2060
|
|
|
2061 void triang_start(des,lf) /* Triangulation with polyhedral start */
|
|
|
2062 design *des;
|
|
|
2063 lfit *lf;
|
|
|
2064 {
|
|
|
2065 int i, j, k, n, d, nc, nvm, ncm, vc;
|
|
|
2066 int *ce, ed[1+MXDIM];
|
|
|
2067 double V[MXDIM*MXDIM], P[MXDIM*MXDIM], sigma, z[MXDIM], xa[1+MXDIM], *xev;
|
|
|
2068 xev = evp(&lf->fp);
|
|
|
2069 d = lf->lfd.d; n = lf->lfd.n;
|
|
|
2070 lf->fp.nv = nc = 0;
|
|
|
2071
|
|
|
2072 triang_guessnv(&nvm,&ncm,&vc,d,mk(&lf->evs));
|
|
|
2073 trchck(lf,nvm,ncm,vc);
|
|
|
2074
|
|
|
2075 ce = lf->evs.ce;
|
|
|
2076 for (j=0; j<d; j++) xev[j] = lf->pc.xbar[j];
|
|
|
2077 lf->fp.nv = 1;
|
|
|
2078 covrofdata(&lf->lfd,V,xev); /* fix this with scaling */
|
|
|
2079 eig_dec(V,P,d);
|
|
|
2080
|
|
|
2081 for (i=0; i<d; i++) /* add vertices +- 2sigma*eigenvect */
|
|
|
2082 { sigma = sqrt(V[i*(d+1)]);
|
|
|
2083 for (j=0; j<d; j++)
|
|
|
2084 xev[lf->fp.nv*d+j] = xev[j]-2*sigma*P[j*d+i];
|
|
|
2085 lf->fp.nv++;
|
|
|
2086 for (j=0; j<d; j++)
|
|
|
2087 xev[lf->fp.nv*d+j] = xev[j]+2*sigma*P[j*d+i];
|
|
|
2088 lf->fp.nv++;
|
|
|
2089 }
|
|
|
2090
|
|
|
2091 for (i=0; i<n; i++) /* is point i inside? */
|
|
|
2092 { ed[0] = 0;
|
|
|
2093 for (j=0; j<d; j++)
|
|
|
2094 { z[j] = 0;
|
|
|
2095 for (k=0; k<d; k++) z[j] += P[k*d+j]*(datum(&lf->lfd,k,i)-xev[k]);
|
|
|
2096 ed[j+1] = 2*j+1+(z[j]>0);
|
|
|
2097 for (k=0; k<d; k++) z[j] = datum(&lf->lfd,j,i);
|
|
|
2098 }
|
|
|
2099 k = intri(z,ed,xev,xa,d);
|
|
|
2100 if (xa[0]<0)
|
|
|
2101 { for (j=1; j<=d; j++)
|
|
|
2102 for (k=0; k<d; k++)
|
|
|
2103 xev[ed[j]*d+k] = xa[0]*xev[k]+(1-xa[0])*xev[ed[j]*d+k];
|
|
|
2104 }
|
|
|
2105 }
|
|
|
2106
|
|
|
2107 nc = 1<<d; /* create initial cells */
|
|
|
2108 for (i=0; i<nc; i++)
|
|
|
2109 { ce[i*vc] = 0; k = i;
|
|
|
2110 for (j=0; j<d; j++)
|
|
|
2111 { ce[i*vc+j+1] = 2*j+(k%2)+1;
|
|
|
2112 k>>=1;
|
|
|
2113 }
|
|
|
2114 }
|
|
|
2115
|
|
|
2116 for (i=0; i<lf->fp.nv; i++)
|
|
|
2117 { PROC_VERTEX(des,lf,i);
|
|
|
2118 if (lf_error) return;
|
|
|
2119 lf->evs.s[i] = 0;
|
|
|
2120 }
|
|
|
2121 for (i=0; i<nc; i++)
|
|
|
2122 triang_grow(des,lf,&ce[i*vc],NULL,NULL);
|
|
|
2123 lf->evs.nce = nc;
|
|
|
2124 }
|
|
|
2125
|
|
|
2126 double triang_cubicint(v,vv,w,d,nc,xxa)
|
|
|
2127 double *v, *vv, *xxa;
|
|
|
2128 int *w, d, nc;
|
|
|
2129 { double sa, lb, *vert0, *vert1, *vals0, *vals1, deriv0, deriv1;
|
|
|
2130 int i, j, k;
|
|
|
2131 if (nc==1) /* linear interpolate */
|
|
|
2132 { sa = 0;
|
|
|
2133 for (i=0; i<=d; i++) sa += xxa[i]*vv[i];
|
|
|
2134 return(sa);
|
|
|
2135 }
|
|
|
2136 sa = 1.0;
|
|
|
2137 for (j=d; j>0; j--) /* eliminate v[w[j]] */
|
|
|
2138 { lb = xxa[j]/sa;
|
|
|
2139 for (k=0; k<j; k++) /* Interpolate edge v[w[k]],v[w[j]] */
|
|
|
2140 { vert0 = &v[w[k]*d];
|
|
|
2141 vert1 = &v[w[j]*d];
|
|
|
2142 vals0 = &vv[k*nc];
|
|
|
2143 vals1 = &vv[j*nc];
|
|
|
2144 deriv0 = deriv1 = 0;
|
|
|
2145 for (i=0; i<d; i++)
|
|
|
2146 { deriv0 += (vert1[i]-vert0[i])*vals0[i+1];
|
|
|
2147 deriv1 += (vert1[i]-vert0[i])*vals1[i+1];
|
|
|
2148 }
|
|
|
2149 vals0[0] = cubic_interp(lb,vals0[0],vals1[0],deriv0,deriv1);
|
|
|
2150 for (i=1; i<=d; i++)
|
|
|
2151 vals0[i] = (1-lb)*((1-lb)*vals0[i]+lb*vals1[i]);
|
|
|
2152 }
|
|
|
2153 sa -= xxa[j];
|
|
|
2154 if (sa<=0) j = 0;
|
|
|
2155 }
|
|
|
2156 return(vals0[0]);
|
|
|
2157 }
|
|
|
2158
|
|
|
2159 double triang_clotoch(xev,vv,ce,p,xxa)
|
|
|
2160 double *xev, *vv, *xxa;
|
|
|
2161 int *ce, p;
|
|
|
2162 { double cfo[3], cfe[3], cg[9], *va, *vb, *vc,
|
|
|
2163 l0, nm[3], na, nb, nc, *xl, *xr, *xz, d0, d1, lb, dlt, gam;
|
|
|
2164 int i, w[3], cfl, cfr;
|
|
|
2165 if (p==1)
|
|
|
2166 return(xxa[0]*vv[0]+xxa[1]*vv[1]+xxa[2]*vv[2]);
|
|
|
2167 if (xxa[2]<=MIN(xxa[0],xxa[1]))
|
|
|
2168 { va = &xev[2*ce[0]]; vb = &xev[2*ce[1]]; vc = &xev[2*ce[2]];
|
|
|
2169 w[0] = 0; w[1] = 3; w[2] = 6;
|
|
|
2170 }
|
|
|
2171 else
|
|
|
2172 if (xxa[1]<xxa[0])
|
|
|
2173 { w[0] = 0; w[1] = 6; w[2] = 3;
|
|
|
2174 va = &xev[2*ce[0]]; vb = &xev[2*ce[2]]; vc = &xev[2*ce[1]];
|
|
|
2175 lb = xxa[1]; xxa[1] = xxa[2]; xxa[2] = lb;
|
|
|
2176 }
|
|
|
2177 else
|
|
|
2178 { w[0] = 6; w[1] = 3; w[2] = 0;
|
|
|
2179 va = &xev[2*ce[2]]; vb = &xev[2*ce[1]]; vc = &xev[2*ce[0]];
|
|
|
2180 lb = xxa[0]; xxa[0] = xxa[2]; xxa[2] = lb;
|
|
|
2181 }
|
|
|
2182
|
|
|
2183 /* set cg to values and derivatives on standard triangle */
|
|
|
2184 for (i=0; i<3; i++)
|
|
|
2185 { cg[3*i] = vv[w[i]];
|
|
|
2186 cg[3*i+1] = ((vb[0]-va[0])*vv[w[i]+1]
|
|
|
2187 +(vb[1]-va[1])*vv[w[i]+2])/2; /* df/dx */
|
|
|
2188 cg[3*i+2] = ((2*vc[0]-vb[0]-va[0])*vv[w[i]+1]
|
|
|
2189 +(2*vc[1]-vb[1]-va[1])*vv[w[i]+2])/2.0; /* sqrt{3} df/dy */
|
|
|
2190 }
|
|
|
2191 dlt = (vb[0]-va[0])*(vc[1]-va[1])-(vc[0]-va[0])*(vb[1]-va[1]);
|
|
|
2192 /* Twice area; +ve if abc antic.wise -ve is abc c.wise */
|
|
|
2193 cfo[0] = (cg[0]+cg[3]+cg[6])/3;
|
|
|
2194 cfo[1] = (2*cg[0]-cg[3]-cg[6])/4;
|
|
|
2195 cfo[2] = (2*cg[3]-cg[0]-cg[6])/4;
|
|
|
2196 na = -cg[1]+cg[2]; /* perp. deriv, rel. length 2 */
|
|
|
2197 nb = -cg[4]-cg[5];
|
|
|
2198 nc = 2*cg[7];
|
|
|
2199 cfo[1] += (nb-nc)/16;
|
|
|
2200 cfo[2] += (nc-na)/16;
|
|
|
2201 na = -cg[1]-cg[2]/3.0; /* derivatives back to origin */
|
|
|
2202 nb = cg[4]-cg[5]/3.0;
|
|
|
2203 nc = cg[8]/1.5;
|
|
|
2204 cfo[0] -= (na+nb+nc)*7/54;
|
|
|
2205 cfo[1] += 13*(nb+nc-2*na)/144;
|
|
|
2206 cfo[2] += 13*(na+nc-2*nb)/144;
|
|
|
2207 for (i=0; i<3; i++)
|
|
|
2208 { /* Outward normals by linear interpolation on original triangle.
|
|
|
2209 Convert to outward normals on standard triangle.
|
|
|
2210 Actually, computed to opposite corner */
|
|
|
2211 switch(i)
|
|
|
2212 { case 0: xl = vc; xr = vb; xz = va; cfl = w[2]; cfr = w[1];
|
|
|
2213 break;
|
|
|
2214 case 1: xl = va; xr = vc; xz = vb; cfl = w[0]; cfr = w[2];
|
|
|
2215 break;
|
|
|
2216 case 2: xl = vb; xr = va; xz = vc; cfl = w[1]; cfr = w[0];
|
|
|
2217 break;
|
|
|
2218 }
|
|
|
2219 na = xr[0]-xl[0]; nb = xr[1]-xl[1];
|
|
|
2220 lb = na*na+nb*nb;
|
|
|
2221 d0 = 1.5*(vv[cfr]-vv[cfl]) - 0.25*(na*(vv[cfl+1]+vv[cfr+1])
|
|
|
2222 +nb*(vv[cfl+2]+vv[cfr+2]));
|
|
|
2223 d1 = 0.5*( na*(vv[cfl+2]+vv[cfr+2])-nb*(vv[cfl+1]+vv[cfr+1]) );
|
|
|
2224 l0 = (xz[0]-xl[0])*na+(xz[1]-xl[1])*nb-lb/2;
|
|
|
2225 nm[i] = (d1*dlt-l0*d0)/lb;
|
|
|
2226 }
|
|
|
2227 cfo[0] -= (nm[0]+nm[1]+nm[2])*4/81;
|
|
|
2228 cfo[1] += (2*nm[0]-nm[1]-nm[2])/27;
|
|
|
2229 cfo[2] += (2*nm[1]-nm[0]-nm[2])/27;
|
|
|
2230
|
|
|
2231 gam = xxa[0]+xxa[1]-2*xxa[2];
|
|
|
2232 if (gam==0) return(cfo[0]);
|
|
|
2233 lb = (xxa[0]-xxa[2])/gam;
|
|
|
2234 d0 = -2*cg[4]; d1 = -2*cg[1];
|
|
|
2235 cfe[0] = cubic_interp(lb,cg[3],cg[0],d0,d1);
|
|
|
2236 cfe[1] = cubintd(lb,cg[3],cg[0],d0,d1);
|
|
|
2237 cfe[2] = -(1-lb)*(1-2*lb)*cg[5] + 4*lb*(1-lb)*nm[2] - lb*(2*lb-1)*cg[2];
|
|
|
2238 d0 = 2*(lb*cfo[1]+(1-lb)*cfo[2]);
|
|
|
2239 d1 = (lb-0.5)*cfe[1]+cfe[2]/3.0;
|
|
|
2240 return(cubic_interp(gam,cfo[0],cfe[0],d0,d1));
|
|
|
2241 }
|
|
|
2242
|
|
|
2243 int triang_getvertexvals(fp,evs,vv,i,what)
|
|
|
2244 fitpt *fp;
|
|
|
2245 evstruc *evs;
|
|
|
2246 double *vv;
|
|
|
2247 int i, what;
|
|
|
2248 { double dx, P, le, vl[1+MXDIM], vh[1+MXDIM];
|
|
|
2249 int d, il, ih, j, nc;
|
|
|
2250 d = fp->d;
|
|
|
2251 if (evs->s[i]==0) return(exvval(fp,vv,i,d,what,0));
|
|
|
2252
|
|
|
2253 il = evs->lo[i]; nc = triang_getvertexvals(fp,evs,vl,il,what);
|
|
|
2254 ih = evs->hi[i]; nc = triang_getvertexvals(fp,evs,vh,ih,what);
|
|
|
2255 vv[0] = (vl[0]+vh[0])/2;
|
|
|
2256 if (nc==1) return(nc);
|
|
|
2257 P = 1.5*(vh[0]-vl[0]);
|
|
|
2258 le = 0.0;
|
|
|
2259 for (j=0; j<d; j++)
|
|
|
2260 { dx = evptx(fp,ih,j)-evptx(fp,il,j);
|
|
|
2261 vv[0] += dx*(vl[j+1]-vh[j+1])/8;
|
|
|
2262 vv[j+1] = (vl[j+1]+vh[j+1])/2;
|
|
|
2263 P -= 1.5*dx*vv[j+1];
|
|
|
2264 le += dx*dx;
|
|
|
2265 }
|
|
|
2266 for (j=0; j<d; j++)
|
|
|
2267 vv[j+1] += P*(evptx(fp,ih,j)-evptx(fp,il,j))/le;
|
|
|
2268 return(nc);
|
|
|
2269 }
|
|
|
2270
|
|
|
2271 double triang_int(lf,x,what)
|
|
|
2272 lfit *lf;
|
|
|
2273 double *x;
|
|
|
2274 int what;
|
|
|
2275 {
|
|
|
2276 int d, i, j, k, vc, nc;
|
|
|
2277 int *ce, nce[1+MXDIM];
|
|
|
2278 double xa[1+MXDIM], vv[(1+MXDIM)*(1+MXDIM)], lb;
|
|
|
2279 fitpt *fp;
|
|
|
2280 evstruc *evs;
|
|
|
2281 fp = &lf->fp;
|
|
|
2282 evs= &lf->evs;
|
|
|
2283
|
|
|
2284 d = fp->d; vc = d+1;
|
|
|
2285 ce = evs->ce;
|
|
|
2286 i = 0;
|
|
|
2287 while ((i<evs->nce) && (!intri(x,&ce[i*vc],evp(fp),xa,d))) i++;
|
|
|
2288 if (i==evs->nce) return(NOSLN);
|
|
|
2289 i *= vc;
|
|
|
2290 for (j=0; j<vc; j++) nce[j] = ce[i+j];
|
|
|
2291 triang_descend(lf,xa,nce);
|
|
|
2292
|
|
|
2293 /* order the vertices -- needed for asymmetric interptr */
|
|
|
2294 do
|
|
|
2295 { k=0;
|
|
|
2296 for (i=0; i<d; i++)
|
|
|
2297 if (nce[i]>nce[i+1])
|
|
|
2298 { j=nce[i]; nce[i]=nce[i+1]; nce[i+1]=j; k=1;
|
|
|
2299 lb = xa[i]; xa[i] = xa[i+1]; xa[i+1] = lb;
|
|
|
2300 }
|
|
|
2301 } while(k);
|
|
|
2302 nc = 0;
|
|
|
2303 for (i=0; i<vc; i++)
|
|
|
2304 nc = triang_getvertexvals(fp,evs,&vv[i*nc],nce[i],what);
|
|
|
2305 return((d==2) ? triang_clotoch(evp(fp),vv,nce,nc,xa) :
|
|
|
2306 triang_cubicint(evp(fp),vv,nce,d,nc,xa));
|
|
|
2307 }
|
|
|
2308 /*
|
|
|
2309 * Copyright 1996-2006 Catherine Loader.
|
|
|
2310 */
|
|
|
2311 /*
|
|
|
2312 * Functions for computing residuals and fitted values from
|
|
|
2313 * the locfit object.
|
|
|
2314 *
|
|
|
2315 * fitted(lf,fit,what,cv,ty) computes fitted values from the
|
|
|
2316 * fit structure in lf.
|
|
|
2317 * resid(y,c,w,th,mi,ty) converts fitted values to residuals
|
|
|
2318 */
|
|
|
2319
|
|
|
2320 #include "lfev.h"
|
|
|
2321
|
|
|
2322 #define NRT 8
|
|
|
2323 static char *rtype[NRT] = { "deviance", "d2", "pearson", "raw",
|
|
|
2324 "ldot", "lddot", "fit", "mean" };
|
|
|
2325 static int rvals[NRT] = { RDEV, RDEV2, RPEAR, RRAW, RLDOT, RLDDT, RFIT, RMEAN};
|
|
|
2326 int restyp(z)
|
|
|
2327 char *z;
|
|
|
2328 { int val;
|
|
|
2329
|
|
|
2330 val = pmatch(z, rtype, rvals, NRT, -1);
|
|
|
2331 if (val==-1) LERR(("Unknown type = %s",z));
|
|
|
2332 return(val);
|
|
|
2333 }
|
|
|
2334
|
|
|
2335 double resid(y,w,th,fam,ty,res)
|
|
|
2336 int fam, ty;
|
|
|
2337 double y, w, th, *res;
|
|
|
2338 { double raw;
|
|
|
2339
|
|
|
2340 fam = fam & 63;
|
|
|
2341 if ((fam==TGAUS) | (fam==TROBT) | (fam==TCAUC))
|
|
|
2342 raw = y-res[ZMEAN];
|
|
|
2343 else
|
|
|
2344 raw = y-w*res[ZMEAN];
|
|
|
2345 switch(ty)
|
|
|
2346 { case RDEV:
|
|
|
2347 if (res[ZDLL]>0) return(sqrt(-2*res[ZLIK]));
|
|
|
2348 else return(-sqrt(-2*res[ZLIK]));
|
|
|
2349 case RPEAR:
|
|
|
2350 if (res[ZDDLL]<=0)
|
|
|
2351 { if (res[ZDLL]==0) return(0);
|
|
|
2352 return(NOSLN);
|
|
|
2353 }
|
|
|
2354 return(res[ZDLL]/sqrt(res[ZDDLL]));
|
|
|
2355 case RRAW: return(raw);
|
|
|
2356 case RLDOT: return(res[ZDLL]);
|
|
|
2357 case RDEV2: return(-2*res[ZLIK]);
|
|
|
2358 case RLDDT: return(res[ZDDLL]);
|
|
|
2359 case RFIT: return(th);
|
|
|
2360 case RMEAN: return(res[ZMEAN]);
|
|
|
2361 default: LERR(("resid: unknown residual type %d",ty));
|
|
|
2362 }
|
|
|
2363 return(0.0);
|
|
|
2364 }
|
|
|
2365
|
|
|
2366 double studentize(res,inl,var,ty,link)
|
|
|
2367 double res, inl, var, *link;
|
|
|
2368 int ty;
|
|
|
2369 { double den;
|
|
|
2370 inl *= link[ZDDLL];
|
|
|
2371 var = var*var*link[ZDDLL];
|
|
|
2372 if (inl>1) inl = 1;
|
|
|
2373 if (var>inl) var = inl;
|
|
|
2374 den = 1-2*inl+var;
|
|
|
2375 if (den<0) return(0.0);
|
|
|
2376 switch(ty)
|
|
|
2377 { case RDEV:
|
|
|
2378 case RPEAR:
|
|
|
2379 case RRAW:
|
|
|
2380 case RLDOT:
|
|
|
2381 return(res/sqrt(den));
|
|
|
2382 case RDEV2:
|
|
|
2383 return(res/den);
|
|
|
2384 default: return(res);
|
|
|
2385 }
|
|
|
2386 }
|
|
|
2387
|
|
|
2388 void fitted(lf,fit,what,cv,st,ty)
|
|
|
2389 lfit *lf;
|
|
|
2390 double *fit;
|
|
|
2391 int what, cv, st, ty;
|
|
|
2392 { int i, j, d, n, evo;
|
|
|
2393 double xx[MXDIM], th, inl, var, link[LLEN];
|
|
|
2394 n = lf->lfd.n;
|
|
|
2395 d = lf->lfd.d;
|
|
|
2396 evo = ev(&lf->evs);
|
|
|
2397 cv &= (evo!=ECROS);
|
|
|
2398 if ((evo==EDATA)|(evo==ECROS)) evo = EFITP;
|
|
|
2399 setfamily(&lf->sp);
|
|
|
2400
|
|
|
2401 for (i=0; i<n; i++)
|
|
|
2402 { for (j=0; j<d; j++) xx[j] = datum(&lf->lfd,j,i);
|
|
|
2403 th = dointpoint(lf,xx,what,evo,i);
|
|
|
2404 if ((what==PT0)|(what==PVARI)) th = th*th;
|
|
|
2405 if (what==PCOEF)
|
|
|
2406 {
|
|
|
2407 th += base(&lf->lfd,i);
|
|
|
2408 stdlinks(link,&lf->lfd,&lf->sp,i,th,rsc(&lf->fp));
|
|
|
2409 if ((cv)|(st))
|
|
|
2410 { inl = dointpoint(lf,xx,PT0,evo,i);
|
|
|
2411 inl = inl*inl;
|
|
|
2412 if (cv)
|
|
|
2413 { th -= inl*link[ZDLL];
|
|
|
2414 stdlinks(link,&lf->lfd,&lf->sp,i,th,rsc(&lf->fp));
|
|
|
2415 }
|
|
|
2416 if (st) var = dointpoint(lf,xx,PNLX,evo,i);
|
|
|
2417 }
|
|
|
2418 fit[i] = resid(resp(&lf->lfd,i),prwt(&lf->lfd,i),th,fam(&lf->sp),ty,link);
|
|
|
2419 if (st) fit[i] = studentize(fit[i],inl,var,ty,link);
|
|
|
2420 } else fit[i] = th;
|
|
|
2421 if (lf_error) return;
|
|
|
2422 }
|
|
|
2423 }
|
|
|
2424 /*
|
|
|
2425 * Copyright 1996-2006 Catherine Loader.
|
|
|
2426 */
|
|
|
2427 #include "lfev.h"
|
|
|
2428
|
|
|
2429 extern double robscale;
|
|
|
2430
|
|
|
2431 /* special version of ressumm to estimate sigma^2, with derivative estimation */
|
|
|
2432 void ressummd(lf)
|
|
|
2433 lfit *lf;
|
|
|
2434 { int i;
|
|
|
2435 double s0, s1;
|
|
|
2436 s0 = s1 = 0.0;
|
|
|
2437 if ((fam(&lf->sp)&64)==0)
|
|
|
2438 { rv(&lf->fp) = 1.0;
|
|
|
2439 return;
|
|
|
2440 }
|
|
|
2441 for (i=0; i<lf->fp.nv; i++)
|
|
|
2442 { s0 += lf->fp.lik[2*lf->fp.nvm+i];
|
|
|
2443 s1 += lf->fp.lik[i];
|
|
|
2444 }
|
|
|
2445 if (s0==0.0)
|
|
|
2446 rv(&lf->fp) = 0.0;
|
|
|
2447 else
|
|
|
2448 rv(&lf->fp) = -2*s1/s0;
|
|
|
2449 }
|
|
|
2450
|
|
|
2451 /*
|
|
|
2452 * res[0] = log-likelihood.
|
|
|
2453 * res[1] = df0 = tr(H)
|
|
|
2454 * res[2] = df0 = tr(H'H)
|
|
|
2455 * res[3] = s^2.
|
|
|
2456 * res[5] = robustscale.
|
|
|
2457 */
|
|
|
2458 void ressumm(lf,des,res)
|
|
|
2459 lfit *lf;
|
|
|
2460 design *des;
|
|
|
2461 double *res;
|
|
|
2462 { int i, j, evo, tg;
|
|
|
2463 double *oy, pw, r1, r2, rdf, t0, t1, u[MXDIM], link[LLEN];
|
|
|
2464 fitpt *fp;
|
|
|
2465
|
|
|
2466 res[0] = res[1] = res[2] = 0.0;
|
|
|
2467
|
|
|
2468 evo = ev(&lf->evs);
|
|
|
2469 if ((evo==EKDCE) | (evo==EPRES))
|
|
|
2470 { res[3] = 1.0;
|
|
|
2471 return;
|
|
|
2472 }
|
|
|
2473 if (lf->dv.nd>0)
|
|
|
2474 { ressummd(lf);
|
|
|
2475 return;
|
|
|
2476 }
|
|
|
2477
|
|
|
2478 r1 = r2 = 0.0;
|
|
|
2479 if ((evo==EDATA) | (evo==ECROS)) evo = EFITP;
|
|
|
2480
|
|
|
2481 for (i=0; i<lf->lfd.n; i++)
|
|
|
2482 { for (j=0; j<lf->lfd.d; j++) u[j] = datum(&lf->lfd,j,i);
|
|
|
2483 fitv(des,i) = base(&lf->lfd,i)+dointpoint(lf,u,PCOEF,evo,i);
|
|
|
2484 des->wd[i] = resp(&(lf->lfd),i) - fitv(des,i);
|
|
|
2485 wght(des,i) = 1.0;
|
|
|
2486 des->ind[i] = i;
|
|
|
2487 }
|
|
|
2488
|
|
|
2489 tg = fam(&lf->sp);
|
|
|
2490 res[5] = 1.0;
|
|
|
2491 if ((tg==TROBT+64) | (tg==TCAUC+64)) /* global robust scale */
|
|
|
2492 { oy = lf->lfd.y; lf->lfd.y = des->wd;
|
|
|
2493 des->xev = lf->pc.xbar;
|
|
|
2494 locfit(&lf->lfd,des,&lf->sp,1,0,0);
|
|
|
2495 lf->lfd.y = oy;
|
|
|
2496 res[5] = robscale;
|
|
|
2497 }
|
|
|
2498
|
|
|
2499 for (i=0; i<lf->lfd.n; i++)
|
|
|
2500 { for (j=0; j<lf->lfd.d; j++) u[j] = datum(&lf->lfd,j,i);
|
|
|
2501 t0 = dointpoint(lf,u,PT0,evo,i);
|
|
|
2502 t1 = dointpoint(lf,u,PNLX,evo,i);
|
|
|
2503 stdlinks(link,&lf->lfd,&lf->sp,i,fitv(des,i),res[5]);
|
|
|
2504 t1 = t1*t1*link[ZDDLL];
|
|
|
2505 t0 = t0*t0*link[ZDDLL];
|
|
|
2506 if (t1>1) t1 = 1;
|
|
|
2507 if (t0>1) t0 = 1; /* no observation gives >1 deg.free */
|
|
|
2508 res[0] += link[ZLIK];
|
|
|
2509 res[1] += t0;
|
|
|
2510 res[2] += t1;
|
|
|
2511 pw = prwt(&lf->lfd,i);
|
|
|
2512 if (pw>0)
|
|
|
2513 { r1 += link[ZDLL]*link[ZDLL]/pw;
|
|
|
2514 r2 += link[ZDDLL]/pw;
|
|
|
2515 }
|
|
|
2516 }
|
|
|
2517
|
|
|
2518 res[3] = 1.0;
|
|
|
2519 if ((fam(&lf->sp)&64)==64) /* quasi family */
|
|
|
2520 { rdf = lf->lfd.n-2*res[1]+res[2];
|
|
|
2521 if (rdf<1.0)
|
|
|
2522 { WARN(("Estimated rdf < 1.0; not estimating variance"));
|
|
|
2523 }
|
|
|
2524 else
|
|
|
2525 res[3] = r1/r2 * lf->lfd.n / rdf;
|
|
|
2526 }
|
|
|
2527
|
|
|
2528 /* try to ensure consistency for family="circ"! */
|
|
|
2529 if (((fam(&lf->sp)&63)==TCIRC) & (lf->lfd.d==1))
|
|
|
2530 { int *ind, nv;
|
|
|
2531 double dlt, th0, th1;
|
|
|
2532 ind = des->ind;
|
|
|
2533 nv = fp->nv;
|
|
|
2534 for (i=0; i<nv; i++) ind[i] = i;
|
|
|
2535 lforder(ind,evp(fp),0,nv-1);
|
|
|
2536 for (i=1; i<nv; i++)
|
|
|
2537 { dlt = evptx(fp,ind[i],0)-evptx(fp,ind[i-1],0);
|
|
|
2538 th0 = fp->coef[ind[i]]-dlt*fp->coef[ind[i]+nv]-fp->coef[ind[i-1]];
|
|
|
2539 th1 = fp->coef[ind[i]]-dlt*fp->coef[ind[i-1]+nv]-fp->coef[ind[i-1]];
|
|
|
2540 if ((th0>PI)&(th1>PI))
|
|
|
2541 { for (j=0; j<i; j++)
|
|
|
2542 fp->coef[ind[j]] += 2*PI;
|
|
|
2543 i--;
|
|
|
2544 }
|
|
|
2545 if ((th0<(-PI))&(th1<(-PI)))
|
|
|
2546 { for (j=0; j<i; j++)
|
|
|
2547 fp->coef[ind[j]] -= 2*PI;
|
|
|
2548 i--;
|
|
|
2549 }
|
|
|
2550 }
|
|
|
2551 }
|
|
|
2552
|
|
|
2553 }
|
|
|
2554
|
|
|
2555 double rss(lf,des,df)
|
|
|
2556 lfit *lf;
|
|
|
2557 design *des;
|
|
|
2558 double *df;
|
|
|
2559 { double ss, res[10];
|
|
|
2560 ss = 0;
|
|
|
2561 ressumm(lf,des,res);
|
|
|
2562 *df = lf->lfd.n - 2*res[1] + res[2];
|
|
|
2563 return(-2*res[0]);
|
|
|
2564 }
|
|
|
2565 /*
|
|
|
2566 * Copyright 1996-2006 Catherine Loader.
|
|
|
2567 */
|
|
|
2568 /*
|
|
|
2569 *
|
|
|
2570 * Derivative corrections. The local slopes are not the derivatives
|
|
|
2571 * of the local likelihood estimate; the function dercor() computes
|
|
|
2572 * the adjustment to get the correct derivatives under the assumption
|
|
|
2573 * that h is constant.
|
|
|
2574 *
|
|
|
2575 * By differentiating the local likelihood equations, one obtains
|
|
|
2576 *
|
|
|
2577 * d ^ ^ T -1 T d . ^
|
|
|
2578 * -- a = a - (X W V X) X -- W l( Y, X a)
|
|
|
2579 * dx 0 1 dx
|
|
|
2580 */
|
|
|
2581
|
|
|
2582 #include "lfev.h"
|
|
|
2583 extern double robscale;
|
|
|
2584
|
|
|
2585 void dercor(lfd,sp,des,coef)
|
|
|
2586 lfdata *lfd;
|
|
|
2587 smpar *sp;
|
|
|
2588 design *des;
|
|
|
2589 double *coef;
|
|
|
2590 { double s1, dc[MXDIM], wd, link[LLEN];
|
|
|
2591 int i, ii, j, m, d, p;
|
|
|
2592 if (fam(sp)<=THAZ) return;
|
|
|
2593 if (ker(sp)==WPARM) return;
|
|
|
2594
|
|
|
2595 d = lfd->d;
|
|
|
2596 p = des->p; m = des->n;
|
|
|
2597
|
|
|
2598 if (lf_debug>1) mut_printf(" Correcting derivatives\n");
|
|
|
2599 fitfun(lfd, sp, des->xev,des->xev,des->f1,NULL);
|
|
|
2600 jacob_solve(&des->xtwx,des->f1);
|
|
|
2601 setzero(dc,d);
|
|
|
2602
|
|
|
2603 /* correction term is e1^T (XTWVX)^{-1} XTW' ldot. */
|
|
|
2604 for (i=0; i<m; i++)
|
|
|
2605 { s1 = innerprod(des->f1,d_xi(des,i),p);
|
|
|
2606 ii = des->ind[i];
|
|
|
2607 stdlinks(link,lfd,sp,ii,fitv(des,ii),robscale);
|
|
|
2608 for (j=0; j<d; j++)
|
|
|
2609 { wd = wght(des,ii)*weightd(datum(lfd,j,ii)-des->xev[j],lfd->sca[j],
|
|
|
2610 d,ker(sp),kt(sp),des->h,lfd->sty[j],dist(des,ii));
|
|
|
2611 dc[j] += s1*wd*link[ZDLL];
|
|
|
2612 }
|
|
|
2613
|
|
|
2614 }
|
|
|
2615 for (j=0; j<d; j++) coef[j+1] += dc[j];
|
|
|
2616 }
|
|
|
2617 /*
|
|
|
2618 * Copyright 1996-2006 Catherine Loader.
|
|
|
2619 */
|
|
|
2620 #include "lfev.h"
|
|
|
2621
|
|
|
2622 void allocallcf(lf)
|
|
|
2623 lfit *lf;
|
|
|
2624 { lf->fp.coef = VVEC(lf,0);
|
|
|
2625 lf->fp.h = VVEC(lf,NPAR(lf));
|
|
|
2626 }
|
|
|
2627
|
|
|
2628 int procvallcf(des,lf,v)
|
|
|
2629 design *des;
|
|
|
2630 lfit *lf;
|
|
|
2631 int v;
|
|
|
2632 { int i, p, lf_status;
|
|
|
2633
|
|
|
2634 lf_status = procv_nov(des,lf,v);
|
|
|
2635 if (lf_error) return(lf_status);
|
|
|
2636
|
|
|
2637 p = NPAR(lf);
|
|
|
2638 for (i=0; i<p; i++) VVAL(lf,v,i) = des->cf[i];
|
|
|
2639 lf->fp.h[v] = des->h;
|
|
|
2640
|
|
|
2641 return(0);
|
|
|
2642 }
|
|
|
2643
|
|
|
2644 void initallcf(lf)
|
|
|
2645 lfit *lf;
|
|
|
2646 { PROCV(lf) = procvallcf;
|
|
|
2647 ALLOC(lf) = allocallcf;
|
|
|
2648 PPROC(lf) = NULL;
|
|
|
2649 KEEPV(lf) = NPAR(lf)+1;
|
|
|
2650 KEEPC(lf) = 0;
|
|
|
2651 NOPC(lf) = 1;
|
|
|
2652 }
|
|
|
2653 /*
|
|
|
2654 * Copyright 1996-2006 Catherine Loader.
|
|
|
2655 */
|
|
|
2656 #include "lfev.h"
|
|
|
2657
|
|
|
2658 void pprocgam(lf,des,res)
|
|
|
2659 lfit *lf;
|
|
|
2660 design *des;
|
|
|
2661 double *res;
|
|
|
2662 { int i, j, n, evo, op;
|
|
|
2663 double *fv, *vr, df, t0, t1, u[MXDIM], link[LLEN];
|
|
|
2664
|
|
|
2665 n = lf->lfd.n;
|
|
|
2666 fv = res;
|
|
|
2667 vr = &res[n];
|
|
|
2668 df = 0.0;
|
|
|
2669
|
|
|
2670 evo = ev(&lf->evs);
|
|
|
2671 if (evo==EDATA) evo = EFITP;
|
|
|
2672
|
|
|
2673 for (i=0; i<n; i++)
|
|
|
2674 { for (j=0; j<lf->lfd.d; j++) u[j] = datum(&lf->lfd,j,i);
|
|
|
2675 fitv(des,i) = base(&lf->lfd,i)+dointpoint(lf,u,PCOEF,evo,i);
|
|
|
2676 lf->lfd.y[i] -= fitv(des,i);
|
|
|
2677 wght(des,i) = 1.0;
|
|
|
2678 des->ind[i] = i;
|
|
|
2679
|
|
|
2680 t0 = dointpoint(lf,u,PT0,evo,i);
|
|
|
2681 t1 = dointpoint(lf,u,PNLX,evo,i);
|
|
|
2682 stdlinks(link,&lf->lfd,&lf->sp,i,fitv(des,i),1.0);
|
|
|
2683 t0 = t0*t0*link[ZDDLL];
|
|
|
2684 t1 = t1*t1*link[ZDDLL];
|
|
|
2685 if (t0>1) t0 = 1; /* no observation gives >1 deg.free */
|
|
|
2686 if (t1>1) t1 = 1; /* no observation gives >1 deg.free */
|
|
|
2687 vr[i] = t1;
|
|
|
2688 df += t0;
|
|
|
2689 }
|
|
|
2690
|
|
|
2691 des->n = n;
|
|
|
2692 deg(&lf->sp) = 1;
|
|
|
2693 op = npar(&lf->sp);
|
|
|
2694 npar(&lf->sp) = des->p = 1+lf->lfd.d;
|
|
|
2695 des->xev = lf->pc.xbar;
|
|
|
2696 locfit(&lf->lfd,des,&lf->sp,1,0,0);
|
|
|
2697
|
|
|
2698 for (i=0; i<n; i++) fv[i] = resp(&lf->lfd,i) - fitv(des,i);
|
|
|
2699 for (i=0; i<=lf->lfd.d; i++)
|
|
|
2700 lf->pc.coef[i] += des->cf[i];
|
|
|
2701 res[2*n] = df-2;
|
|
|
2702 npar(&lf->sp) = op;
|
|
|
2703 }
|
|
|
2704
|
|
|
2705 void initgam(lf)
|
|
|
2706 lfit *lf;
|
|
|
2707 { initstd(lf);
|
|
|
2708 PPROC(lf) = pprocgam;
|
|
|
2709 KEEPC(lf) = 2*NOBS(lf)+1;
|
|
|
2710 }
|
|
|
2711 /*
|
|
|
2712 * Copyright 1996-2006 Catherine Loader.
|
|
|
2713 */
|
|
|
2714 #include "lfev.h"
|
|
|
2715
|
|
|
2716 int procvhatm(des,lf,v)
|
|
|
2717 design *des;
|
|
|
2718 lfit *lf;
|
|
|
2719 int v;
|
|
|
2720 { int k;
|
|
|
2721 double *l;
|
|
|
2722 l = &lf->fp.coef[v*lf->lfd.n];
|
|
|
2723 if ((ker(&lf->sp)!=WPARM) | (!haspc(&lf->pc)))
|
|
|
2724 { k = procv_nov(des,lf,v);
|
|
|
2725 wdiag(&lf->lfd,&lf->sp,des,l,&lf->dv,0,1,1);
|
|
|
2726 }
|
|
|
2727 else
|
|
|
2728 wdiagp(&lf->lfd,&lf->sp,des,l,&lf->pc,&lf->dv,0,1,1);
|
|
|
2729 lf->fp.h[v] = des->h;
|
|
|
2730 return(k);
|
|
|
2731 }
|
|
|
2732
|
|
|
2733 void allochatm(lf)
|
|
|
2734 lfit *lf;
|
|
|
2735 { lf->fp.coef = VVEC(lf,0);
|
|
|
2736 lf->fp.h = VVEC(lf,NOBS(lf));
|
|
|
2737 }
|
|
|
2738
|
|
|
2739 void pprochatm(lf,des,res)
|
|
|
2740 lfit *lf;
|
|
|
2741 design *des;
|
|
|
2742 double *res;
|
|
|
2743 { transpose(lf->fp.coef,lf->fp.nvm,lf->lfd.n);
|
|
|
2744 }
|
|
|
2745
|
|
|
2746 void inithatm(lf)
|
|
|
2747 lfit *lf;
|
|
|
2748 { PROCV(lf) = procvhatm;
|
|
|
2749 ALLOC(lf) = allochatm;
|
|
|
2750 PPROC(lf) = pprochatm;
|
|
|
2751 KEEPV(lf) = NOBS(lf)+1;
|
|
|
2752 KEEPC(lf) = 1;
|
|
|
2753 NOPC(lf) = 1; /* could use pc if mi[MKER] == WPARM */
|
|
|
2754 }
|
|
|
2755 /*
|
|
|
2756 * Copyright 1996-2006 Catherine Loader.
|
|
|
2757 */
|
|
|
2758 #include "lfev.h"
|
|
|
2759
|
|
|
2760 static lfit *lf_scb;
|
|
|
2761 static lfdata *lfd_scb;
|
|
|
2762 static smpar *scb_sp;
|
|
|
2763 static design *des_scb;
|
|
|
2764
|
|
|
2765 int scbfitter(x,l,reqd)
|
|
|
2766 double *x, *l;
|
|
|
2767 int reqd;
|
|
|
2768 {
|
|
|
2769 int m;
|
|
|
2770 des_scb->xev = x;
|
|
|
2771 if ((ker(scb_sp)!=WPARM) | (!haspc(&lf_scb->pc)))
|
|
|
2772 { locfit(lfd_scb,des_scb,&lf_scb->sp,1,1,0);
|
|
|
2773 m = wdiag(lfd_scb, scb_sp, des_scb,l,&lf_scb->dv,reqd,2,0);
|
|
|
2774 }
|
|
|
2775 else
|
|
|
2776 m = wdiagp(lfd_scb, scb_sp, des_scb,l,&lf_scb->pc,&lf_scb->dv,reqd,2,0);
|
|
|
2777 return(m);
|
|
|
2778 }
|
|
|
2779
|
|
|
2780 int constants(lf,des,res)
|
|
|
2781 lfit *lf;
|
|
|
2782 design *des;
|
|
|
2783 double *res;
|
|
|
2784 {
|
|
|
2785 int d, m, nt;
|
|
|
2786 double *wk;
|
|
|
2787 evstruc *evs;
|
|
|
2788
|
|
|
2789 lf_scb = lf;
|
|
|
2790 des_scb = des;
|
|
|
2791 lfd_scb = &lf->lfd;
|
|
|
2792 scb_sp = &lf->sp;
|
|
|
2793
|
|
|
2794 evs = &lf->evs;
|
|
|
2795 d = lfd_scb->d;
|
|
|
2796 m = lfd_scb->n;
|
|
|
2797 trchck(lf,0,0,0);
|
|
|
2798
|
|
|
2799 if (lf_error) return(0);
|
|
|
2800 if ((ker(scb_sp) != WPARM) && (lf->sp.nn>0))
|
|
|
2801 WARN(("constants are approximate for varying h"));
|
|
|
2802 npar(scb_sp) = calcp(scb_sp,lf->lfd.d);
|
|
|
2803 des_init(des,m,npar(scb_sp));
|
|
|
2804 set_scales(&lf->lfd);
|
|
|
2805 set_flim(&lf->lfd,&lf->evs);
|
|
|
2806 compparcomp(des,&lf->lfd,&lf->sp,&lf->pc,ker(scb_sp)!=WPARM);
|
|
|
2807
|
|
|
2808 wk = &res[d+1];
|
|
|
2809 nt = tube_constants(scbfitter,d,m,ev(evs),mg(evs),evs->fl,
|
|
|
2810 res,wk,(d>3) ? 4 : d+1,0);
|
|
|
2811 lf->evs.nce = nt; /* cheat way to return it to S. */
|
|
|
2812 return(nt);
|
|
|
2813 }
|
|
|
2814
|
|
|
2815 void initkappa(lf)
|
|
|
2816 lfit *lf;
|
|
|
2817 { PROCV(lf) = NULL;
|
|
|
2818 ALLOC(lf) = NULL;
|
|
|
2819 PPROC(lf) = (void *)constants;
|
|
|
2820 KEEPV(lf) = 0;
|
|
|
2821 KEEPC(lf) = NVAR(lf)+1+k0_reqd(NVAR(lf),NOBS(lf),0);
|
|
|
2822 NOPC(lf) = 0;
|
|
|
2823 }
|
|
|
2824 /*
|
|
|
2825 * Copyright 1996-2006 Catherine Loader.
|
|
|
2826 */
|
|
|
2827 #include "lfev.h"
|
|
|
2828
|
|
|
2829 /* fix df computation for link=IDENT. */
|
|
|
2830 void pproclscv(lf,des,res)
|
|
|
2831 lfit *lf;
|
|
|
2832 design *des;
|
|
|
2833 double *res;
|
|
|
2834 { double df, fh, fh_cv, infl, z0, z1, x[MXDIM];
|
|
|
2835 int i, n, j, evo;
|
|
|
2836 z1 = df = 0.0;
|
|
|
2837 evo = ev(&lf->evs);
|
|
|
2838 n = lf->lfd.n;
|
|
|
2839 if ((evo==EDATA) | (evo==ECROS)) evo = EFITP;
|
|
|
2840
|
|
|
2841 z0 = dens_integrate(lf,des,2);
|
|
|
2842
|
|
|
2843 for (i=0; i<n; i++)
|
|
|
2844 { for (j=0; j<lf->lfd.d; j++) x[j] = datum(&lf->lfd,j,i);
|
|
|
2845 fh = base(&lf->lfd,i)+dointpoint(lf,x,PCOEF,evo,i);
|
|
|
2846 if (link(&lf->sp)==LLOG) fh = exp(fh);
|
|
|
2847 infl = dointpoint(lf,x,PT0,evo,i);
|
|
|
2848 infl = infl * infl;
|
|
|
2849 if (infl>1) infl = 1;
|
|
|
2850 fh_cv = (link(&lf->sp) == LIDENT) ?
|
|
|
2851 (n*fh - infl) / (n-1.0) : fh*(1-infl)*n/(n-1.0);
|
|
|
2852 z1 += fh_cv;
|
|
|
2853 df += infl;
|
|
|
2854 }
|
|
|
2855
|
|
|
2856 res[0] = z0-2*z1/n;
|
|
|
2857 res[1] = df;
|
|
|
2858 }
|
|
|
2859
|
|
|
2860 void initlscv(lf)
|
|
|
2861 lfit *lf;
|
|
|
2862 { initstd(lf);
|
|
|
2863 KEEPC(lf) = 2;
|
|
|
2864 PPROC(lf) = pproclscv;
|
|
|
2865 NOPC(lf) = 1;
|
|
|
2866 }
|
|
|
2867 /*
|
|
|
2868 * Copyright 1996-2006 Catherine Loader.
|
|
|
2869 */
|
|
|
2870 #include "lfev.h"
|
|
|
2871
|
|
|
2872 static double pen, sig2;
|
|
|
2873
|
|
|
2874 void goldensec(f,des,tr,eps,xm,ym,meth)
|
|
|
2875 double (*f)(), eps, *xm, *ym;
|
|
|
2876 int meth;
|
|
|
2877 design *des;
|
|
|
2878 lfit *tr;
|
|
|
2879 { double x[4], y[4], xx[11], yy[11];
|
|
|
2880 int i, im;
|
|
|
2881 xx[0] = tr->sp.fixh;
|
|
|
2882 if (xx[0]<=0)
|
|
|
2883 { LERR(("regband: initialize h>0"));
|
|
|
2884 return;
|
|
|
2885 }
|
|
|
2886 for (i=0; i<=10; i++)
|
|
|
2887 { if (i>0) xx[i] = (1+gold_rat)*xx[i-1];
|
|
|
2888 yy[i] = f(xx[i],des,tr,meth);
|
|
|
2889 if ((i==0) || (yy[i]<yy[im])) im = i;
|
|
|
2890 }
|
|
|
2891 if (im==0) im = 1;
|
|
|
2892 if (im==10)im = 9;
|
|
|
2893 x[0] = xx[im-1]; y[0] = yy[im-1];
|
|
|
2894 x[1] = xx[im]; y[1] = yy[im];
|
|
|
2895 x[3] = xx[im+1]; y[3] = yy[im+1];
|
|
|
2896 x[2] = gold_rat*x[3]+(1-gold_rat)*x[0];
|
|
|
2897 y[2] = f(x[2],des,tr,meth);
|
|
|
2898 while (x[3]-x[0]>eps)
|
|
|
2899 { if (y[1]<y[2])
|
|
|
2900 { x[3] = x[2]; y[3] = y[2];
|
|
|
2901 x[2] = x[1]; y[2] = y[1];
|
|
|
2902 x[1] = gold_rat*x[0]+(1-gold_rat)*x[3];
|
|
|
2903 y[1] = f(x[1],des,tr,meth);
|
|
|
2904 }
|
|
|
2905 else
|
|
|
2906 { x[0] = x[1]; y[0] = y[1];
|
|
|
2907 x[1] = x[2]; y[1] = y[2];
|
|
|
2908 x[2] = gold_rat*x[3]+(1-gold_rat)*x[0];
|
|
|
2909 y[2] = f(x[2],des,tr,meth);
|
|
|
2910 }
|
|
|
2911 }
|
|
|
2912 im = 0;
|
|
|
2913 for (i=1; i<4; i++) if (y[i]<y[im]) im = i;
|
|
|
2914 *xm = x[im]; *ym = y[im];
|
|
|
2915 }
|
|
|
2916
|
|
|
2917 double dnk(x,k)
|
|
|
2918 double x;
|
|
|
2919 int k;
|
|
|
2920 { double f;
|
|
|
2921 switch(k)
|
|
|
2922 { case 0: f = 1; break;
|
|
|
2923 case 1: f = -x; break;
|
|
|
2924 case 2: f = x*x-1; break;
|
|
|
2925 case 3: f = x*(x*x-3); break;
|
|
|
2926 case 4: f = 3-x*x*(6-x*x); break;
|
|
|
2927 case 5: f = -x*(15-x*x*(10-x*x)); break;
|
|
|
2928 case 6: f = -15+x*x*(45-x*x*(15-x*x)); break;
|
|
|
2929 default: LERR(("dnk: k=%d too large",k)); return(0.0);
|
|
|
2930 }
|
|
|
2931 return(f*exp(-x*x/2)/S2PI);
|
|
|
2932 }
|
|
|
2933
|
|
|
2934 double locai(h,des,lf)
|
|
|
2935 double h;
|
|
|
2936 design *des;
|
|
|
2937 lfit *lf;
|
|
|
2938 { double cp, res[10];
|
|
|
2939 nn(&lf->sp) = h;
|
|
|
2940 lf->mdl.procv = procvstd;
|
|
|
2941 nstartlf(des,lf);
|
|
|
2942 ressumm(lf,des,res);
|
|
|
2943 cp = -2*res[0] + pen*res[1];
|
|
|
2944 return(cp);
|
|
|
2945 }
|
|
|
2946
|
|
|
2947 static int fc;
|
|
|
2948
|
|
|
2949 double loccp(h,des,lf,m) /* m=1: cp m=2: gcv */
|
|
|
2950 double h;
|
|
|
2951 design *des;
|
|
|
2952 lfit *lf;
|
|
|
2953 int m;
|
|
|
2954 { double cp, res[10];
|
|
|
2955 int dg, n;
|
|
|
2956
|
|
|
2957 n = lf->lfd.n;
|
|
|
2958 nn(&lf->sp) = 0;
|
|
|
2959 fixh(&lf->sp) = h;
|
|
|
2960 lf->mdl.procv = procvstd;
|
|
|
2961 nstartlf(des,lf);
|
|
|
2962 ressumm(lf,des,res);
|
|
|
2963 if (m==1)
|
|
|
2964 { if (fc) sig2 = res[3]; /* first call - set sig2 */
|
|
|
2965 cp = -2*res[0]/sig2 - n + 2*res[1];
|
|
|
2966 }
|
|
|
2967 else
|
|
|
2968 cp = -2*n*res[0]/((n-res[1])*(n-res[1]));
|
|
|
2969 fc = 0;
|
|
|
2970 return(cp);
|
|
|
2971 }
|
|
|
2972
|
|
|
2973 double cp(des,lf,meth)
|
|
|
2974 design *des;
|
|
|
2975 lfit *lf;
|
|
|
2976 int meth;
|
|
|
2977 { double hm, ym;
|
|
|
2978 fc = 1;
|
|
|
2979 goldensec(loccp,des,lf,0.001,&hm,&ym,meth);
|
|
|
2980 return(hm);
|
|
|
2981 }
|
|
|
2982
|
|
|
2983 double gkk(des,lf)
|
|
|
2984 design *des;
|
|
|
2985 lfit *lf;
|
|
|
2986 { double h, h5, nf, th;
|
|
|
2987 int i, j, n, dg0, dg1;
|
|
|
2988 ev(&lf->evs)= EDATA;
|
|
|
2989 nn(&lf->sp) = 0;
|
|
|
2990 n = lf->lfd.n;
|
|
|
2991 dg0 = deg0(&lf->sp); /* target degree */
|
|
|
2992 dg1 = dg0+1+(dg0%2==0); /* pilot degree */
|
|
|
2993 nf = exp(log(1.0*n)/10); /* bandwidth inflation factor */
|
|
|
2994 h = lf->sp.fixh; /* start bandwidth */
|
|
|
2995 for (i=0; i<=10; i++)
|
|
|
2996 { deg(&lf->sp) = dg1;
|
|
|
2997 lf->sp.fixh = h*nf;
|
|
|
2998 lf->mdl.procv = procvstd;
|
|
|
2999 nstartlf(des,lf);
|
|
|
3000 th = 0;
|
|
|
3001 for (j=10; j<n-10; j++)
|
|
|
3002 th += lf->fp.coef[dg1*n+j]*lf->fp.coef[dg1*n+j];
|
|
|
3003 th *= n/(n-20.0);
|
|
|
3004 h5 = sig2 * Wikk(ker(&lf->sp),dg0) / th;
|
|
|
3005 h = exp(log(h5)/(2*dg1+1));
|
|
|
3006 if (lf_error) return(0.0);
|
|
|
3007 /* mut_printf("pilot %8.5f sel %8.5f\n",lf->sp.fixh,h); */
|
|
|
3008 }
|
|
|
3009 return(h);
|
|
|
3010 }
|
|
|
3011
|
|
|
3012 double rsw(des,lf)
|
|
|
3013 design *des;
|
|
|
3014 lfit *lf;
|
|
|
3015 { int i, j, k, nmax, nvm, n, mk, evo, dg0, dg1;
|
|
|
3016 double rss[6], cp[6], th22, dx, d2, hh;
|
|
|
3017 nmax = 5;
|
|
|
3018 evo = ev(&lf->evs); ev(&lf->evs) = EGRID;
|
|
|
3019 mk = ker(&lf->sp); ker(&lf->sp) = WRECT;
|
|
|
3020 dg0 = deg0(&lf->sp);
|
|
|
3021 dg1 = 1 + dg0 + (dg0%2==0);
|
|
|
3022 deg(&lf->sp) = 4;
|
|
|
3023 for (k=nmax; k>0; k--)
|
|
|
3024 { lf->evs.mg[0] = k;
|
|
|
3025 lf->evs.fl[0] = 1.0/(2*k);
|
|
|
3026 lf->evs.fl[1] = 1-1.0/(2*k);
|
|
|
3027 nn(&lf->sp) = 0;
|
|
|
3028 fixh(&lf->sp) = 1.0/(2*k);
|
|
|
3029 lf->mdl.procv = procvstd;
|
|
|
3030 nstartlf(des,lf);
|
|
|
3031 nvm = lf->fp.nvm;
|
|
|
3032 rss[k] = 0;
|
|
|
3033 for (i=0; i<k; i++) rss[k] += -2*lf->fp.lik[i];
|
|
|
3034 }
|
|
|
3035 n = lf->lfd.n; k = 1;
|
|
|
3036 for (i=1; i<=nmax; i++)
|
|
|
3037 { /* cp[i] = (n-5*nmax)*rss[i]/rss[nmax]-(n-10*i); */
|
|
|
3038 cp[i] = rss[i]/sig2-(n-10*i);
|
|
|
3039 if (cp[i]<cp[k]) k = i;
|
|
|
3040 }
|
|
|
3041 lf->evs.mg[0] = k;
|
|
|
3042 lf->evs.fl[0] = 1.0/(2*k);
|
|
|
3043 lf->evs.fl[1] = 1-1.0/(2*k);
|
|
|
3044 nn(&lf->sp) = 0;
|
|
|
3045 fixh(&lf->sp) = 1.0/(2*k);
|
|
|
3046 lf->mdl.procv = procvstd;
|
|
|
3047 nstartlf(des,lf);
|
|
|
3048 ker(&lf->sp) = mk; ev(&lf->evs) = evo;
|
|
|
3049 nvm = lf->fp.nvm;
|
|
|
3050 th22 = 0;
|
|
|
3051 for (i=10; i<n-10; i++)
|
|
|
3052 { j = floor(k*datum(&lf->lfd,0,i));
|
|
|
3053 if (j>=k) j = k-1;
|
|
|
3054 dx = datum(&lf->lfd,0,i)-evptx(&lf->fp,0,j);
|
|
|
3055 if (dg1==2)
|
|
|
3056 d2 = lf->fp.coef[2*nvm+j]+dx*lf->fp.coef[3*nvm+j]+dx*dx*lf->fp.coef[4*nvm+j]/2;
|
|
|
3057 else d2 = lf->fp.coef[4*nvm+j];
|
|
|
3058 th22 += d2*d2;
|
|
|
3059 }
|
|
|
3060 hh = Wikk(mk,dg0)*sig2/th22*(n-20.0)/n;
|
|
|
3061 return(exp(log(hh)/(2*dg1+1)));
|
|
|
3062 }
|
|
|
3063
|
|
|
3064 #define MAXMETH 10
|
|
|
3065
|
|
|
3066 void rband(lf,des,hhat)
|
|
|
3067 lfit *lf;
|
|
|
3068 design *des;
|
|
|
3069 double *hhat;
|
|
|
3070 { int i, dg, nmeth, meth[MAXMETH];
|
|
|
3071 double h0, res[10];
|
|
|
3072
|
|
|
3073 nmeth = lf->dv.nd;
|
|
|
3074 for (i=0; i<nmeth; i++) meth[i] = lf->dv.deriv[i];
|
|
|
3075 lf->dv.nd = 0;
|
|
|
3076
|
|
|
3077 /* first, estimate sigma^2.
|
|
|
3078 * this is ridiculously bad. lf->sp.fixh = 0.05????
|
|
|
3079 */
|
|
|
3080 /* dg = deg(&lf->sp); deg(&lf->sp) = 2;
|
|
|
3081 h0 = lf->sp.fixh; lf->sp.fixh = 0.05;
|
|
|
3082 lf->mdl.procv = procvstd;
|
|
|
3083 nstartlf(des,lf);
|
|
|
3084 ressumm(lf,des,res);
|
|
|
3085 deg(&lf->sp) = dg; lf->sp.fixh = h0;
|
|
|
3086 sig2 = res[3]; */
|
|
|
3087
|
|
|
3088 for (i=0; i<nmeth; i++)
|
|
|
3089 {
|
|
|
3090 switch(meth[i])
|
|
|
3091 { case 0: hhat[i] = cp(des,lf,1);
|
|
|
3092 break;
|
|
|
3093 case 1: hhat[i] = cp(des,lf,2);
|
|
|
3094 break;
|
|
|
3095 case 2: hhat[i] = gkk(des,lf);
|
|
|
3096 break;
|
|
|
3097 case 3: hhat[i] = rsw(des,lf);
|
|
|
3098 break;
|
|
|
3099 default: hhat[i] = 0;
|
|
|
3100 mut_printf("Unknown method %d\n",meth[i]);
|
|
|
3101 }
|
|
|
3102 if (lf_error) i = nmeth;
|
|
|
3103 }
|
|
|
3104 lf->dv.nd = nmeth;
|
|
|
3105 }
|
|
|
3106
|
|
|
3107 void initrband(lf)
|
|
|
3108 lfit *lf;
|
|
|
3109 {
|
|
|
3110 initstd(lf);
|
|
|
3111 KEEPC(lf) = MAXMETH;
|
|
|
3112 PROCV(lf) = NULL;
|
|
|
3113 PPROC(lf) = rband;
|
|
|
3114 }
|
|
|
3115 /*
|
|
|
3116 * Copyright 1996-2006 Catherine Loader.
|
|
|
3117 */
|
|
|
3118 #include "lfev.h"
|
|
|
3119 static double scb_crit, *x, c[10], kap[5], kaq[5], max_p2;
|
|
|
3120 static int side, type;
|
|
|
3121 design *scb_des;
|
|
|
3122
|
|
|
3123 double covar_par(lf,des,x1,x2)
|
|
|
3124 lfit *lf;
|
|
|
3125 design *des;
|
|
|
3126 double x1, x2;
|
|
|
3127 { double *v1, *v2, *wk;
|
|
|
3128 paramcomp *pc;
|
|
|
3129 int i, j, p, ispar;
|
|
|
3130
|
|
|
3131 v1 = des->f1; v2 = des->ss; wk = des->oc;
|
|
|
3132 ispar = (ker(&lf->sp)==WPARM) && (haspc(&lf->pc));
|
|
|
3133 p = npar(&lf->sp);
|
|
|
3134
|
|
|
3135 /* for parametric models, the covariance is
|
|
|
3136 * A(x1)^T (X^T W V X)^{-1} A(x2)
|
|
|
3137 * which we can find easily from the parametric component.
|
|
|
3138 */
|
|
|
3139 if (ispar)
|
|
|
3140 { pc = &lf->pc;
|
|
|
3141 fitfun(&lf->lfd, &lf->sp, &x1,pc->xbar,v1,NULL);
|
|
|
3142 fitfun(&lf->lfd, &lf->sp, &x2,pc->xbar,v2,NULL);
|
|
|
3143 jacob_hsolve(&lf->pc.xtwx,v1);
|
|
|
3144 jacob_hsolve(&lf->pc.xtwx,v2);
|
|
|
3145 }
|
|
|
3146
|
|
|
3147 /* for non-parametric models, we must use the cholseky decomposition
|
|
|
3148 * of M2 = X^T W^2 V X. Courtesy of lf_vcov caclulations, should have
|
|
|
3149 * des->P = M2^{1/2} M1^{-1}.
|
|
|
3150 */
|
|
|
3151 if (!ispar)
|
|
|
3152 { fitfun(&lf->lfd, &lf->sp, &x1,des->xev,wk,NULL);
|
|
|
3153 for (i=0; i<p; i++)
|
|
|
3154 { v1[i] = 0;
|
|
|
3155 for (j=0; j<p; j++) v1[i] += des->P[i*p+j]*wk[j];
|
|
|
3156 }
|
|
|
3157 fitfun(&lf->lfd, &lf->sp, &x2,des->xev,wk,NULL);
|
|
|
3158 for (i=0; i<p; i++)
|
|
|
3159 { v2[i] = 0;
|
|
|
3160 for (j=0; j<p; j++) v2[i] += des->P[i*p+j]*wk[j];
|
|
|
3161 }
|
|
|
3162 }
|
|
|
3163
|
|
|
3164 return(innerprod(v1,v2,p));
|
|
|
3165 }
|
|
|
3166
|
|
|
3167 void cumulant(lf,des,sd)
|
|
|
3168 lfit *lf;
|
|
|
3169 design *des;
|
|
|
3170 double sd;
|
|
|
3171 { double b2i, b3i, b3j, b4i;
|
|
|
3172 double ss, si, sj, uii, uij, ujj, k1;
|
|
|
3173 int ii, i, j, jj;
|
|
|
3174 for (i=1; i<10; i++) c[i] = 0.0;
|
|
|
3175 k1 = 0;
|
|
|
3176
|
|
|
3177 /* ss = sd*sd; */
|
|
|
3178 ss = covar_par(lf,des,des->xev[0],des->xev[0]);
|
|
|
3179
|
|
|
3180 /*
|
|
|
3181 * this isn't valid for nonparametric models. At a minimum,
|
|
|
3182 * the sums would have to include weights. Still have to work
|
|
|
3183 * out the right way.
|
|
|
3184 */
|
|
|
3185 for (i=0; i<lf->lfd.n; i++)
|
|
|
3186 { ii = des->ind[i];
|
|
|
3187 b2i = b2(fitv(des,ii),fam(&lf->sp),prwt(&lf->lfd,ii));
|
|
|
3188 b3i = b3(fitv(des,ii),fam(&lf->sp),prwt(&lf->lfd,ii));
|
|
|
3189 b4i = b4(fitv(des,ii),fam(&lf->sp),prwt(&lf->lfd,ii));
|
|
|
3190 si = covar_par(lf,des,des->xev[0],datum(&lf->lfd,0,ii));
|
|
|
3191 uii= covar_par(lf,des,datum(&lf->lfd,0,ii),datum(&lf->lfd,0,ii));
|
|
|
3192 if (lf_error) return;
|
|
|
3193
|
|
|
3194 c[2] += b4i*si*si*uii;
|
|
|
3195 c[6] += b4i*si*si*si*si;
|
|
|
3196 c[7] += b3i*si*uii;
|
|
|
3197 c[8] += b3i*si*si*si;
|
|
|
3198 /* c[9] += b2i*si*si*si*si;
|
|
|
3199 c[9] += b2i*b2i*si*si*si*si; */
|
|
|
3200 k1 += b3i*si*(si*si/ss-uii);
|
|
|
3201
|
|
|
3202 /* i=j components */
|
|
|
3203 c[1] += b3i*b3i*si*si*uii*uii;
|
|
|
3204 c[3] += b3i*b3i*si*si*si*si*uii;
|
|
|
3205 c[4] += b3i*b3i*si*si*uii*uii;
|
|
|
3206
|
|
|
3207 for (j=i+1; j<lf->lfd.n; j++)
|
|
|
3208 { jj = des->ind[j];
|
|
|
3209 b3j = b3(fitv(des,ii),fam(&lf->sp),prwt(&lf->lfd,jj));
|
|
|
3210 sj = covar_par(lf,des,des->xev[0],datum(&lf->lfd,0,jj));
|
|
|
3211 uij= covar_par(lf,des,datum(&lf->lfd,0,ii),datum(&lf->lfd,0,jj));
|
|
|
3212 ujj= covar_par(lf,des,datum(&lf->lfd,0,jj),datum(&lf->lfd,0,jj));
|
|
|
3213
|
|
|
3214 c[1] += 2*b3i*b3j*si*sj*uij*uij;
|
|
|
3215 c[3] += 2*b3i*b3j*si*si*sj*sj*uij;
|
|
|
3216 c[4] += b3i*b3j*uij*(si*si*ujj+sj*sj*uii);
|
|
|
3217 if (lf_error) return;
|
|
|
3218 }
|
|
|
3219 }
|
|
|
3220 c[5] = c[1];
|
|
|
3221 c[7] = c[7]*c[8];
|
|
|
3222 c[8] = c[8]*c[8];
|
|
|
3223
|
|
|
3224 c[1] /= ss; c[2] /= ss; c[3] /= ss*ss; c[4] /= ss;
|
|
|
3225 c[5] /= ss; c[6] /= ss*ss; c[7] /= ss*ss;
|
|
|
3226 c[8] /= ss*ss*ss; c[9] /= ss*ss;
|
|
|
3227
|
|
|
3228 /* constants used in p(x,z) computation */
|
|
|
3229 kap[1] = k1/(2*sqrt(ss));
|
|
|
3230 kap[2] = 1 + 0.5*(c[1]-c[2]+c[4]-c[7]) - 3*c[3] + c[6] + 1.75*c[8];
|
|
|
3231 kap[4] = -9*c[3] + 3*c[6] + 6*c[8] + 3*c[9];
|
|
|
3232
|
|
|
3233 /* constants used in q(x,u) computation */
|
|
|
3234 kaq[2] = c[3] - 1.5*c[8] - c[5] - c[4] + 0.5*c[7] + c[6] - c[2];
|
|
|
3235 kaq[4] = -3*c[3] - 6*c[4] - 6*c[5] + 3*c[6] + 3*c[7] - 3*c[8] + 3*c[9];
|
|
|
3236 }
|
|
|
3237
|
|
|
3238 /* q2(u) := u+q2(x,u) in paper */
|
|
|
3239 double q2(u)
|
|
|
3240 double u;
|
|
|
3241 { return(u-u*(36.0*kaq[2] + 3*kaq[4]*(u*u-3) + c[8]*((u*u-10)*u*u+15))/72.0);
|
|
|
3242 }
|
|
|
3243
|
|
|
3244 /* p2(u) := p2(x,u) in paper */
|
|
|
3245 double p2(u)
|
|
|
3246 double u;
|
|
|
3247 { return( -u*( 36*(kap[2]-1+kap[1]*kap[1])
|
|
|
3248 + 3*(kap[4]+4*kap[1]*sqrt(kap[3]))*(u*u-3)
|
|
|
3249 + c[8]*((u*u-10)*u*u+15) ) / 72 );
|
|
|
3250 }
|
|
|
3251
|
|
|
3252 extern int likereg();
|
|
|
3253 double gldn_like(a)
|
|
|
3254 double a;
|
|
|
3255 { int err;
|
|
|
3256
|
|
|
3257 scb_des->fix[0] = 1;
|
|
|
3258 scb_des->cf[0] = a;
|
|
|
3259 max_nr(likereg, scb_des->cf, scb_des->oc, scb_des->res, scb_des->f1,
|
|
|
3260 &scb_des->xtwx, scb_des->p, lf_maxit, 1.0e-6, &err);
|
|
|
3261 scb_des->fix[0] = 0;
|
|
|
3262
|
|
|
3263 return(scb_des->llk);
|
|
|
3264 }
|
|
|
3265
|
|
|
3266 /* v1/v2 is correct for deg=0 only */
|
|
|
3267 void get_gldn(fp,des,lo,hi,v)
|
|
|
3268 fitpt *fp;
|
|
|
3269 design *des;
|
|
|
3270 double *lo, *hi;
|
|
|
3271 int v;
|
|
|
3272 { double v1, v2, c, tlk;
|
|
|
3273 int err;
|
|
|
3274
|
|
|
3275 v1 = fp->nlx[v];
|
|
|
3276 v2 = fp->t0[v];
|
|
|
3277 c = scb_crit * v1 / v2;
|
|
|
3278 tlk = des->llk - c*c/2;
|
|
|
3279 mut_printf("v %8.5f %8.5f c %8.5f tlk %8.5f llk %8.5f\n",v1,v2,c,tlk,des->llk);
|
|
|
3280
|
|
|
3281 /* want: { a : l(a) >= l(a-hat) - c*c/2 } */
|
|
|
3282 lo[v] = fp->coef[v] - scb_crit*v1;
|
|
|
3283 hi[v] = fp->coef[v] + scb_crit*v1;
|
|
|
3284
|
|
|
3285 err = 0;
|
|
|
3286
|
|
|
3287 mut_printf("lo %2d\n",v);
|
|
|
3288 lo[v] = solve_secant(gldn_like,tlk,lo[v],fp->coef[v],1e-8,BDF_EXPLEFT,&err);
|
|
|
3289 if (err>0) mut_printf("solve_secant error\n");
|
|
|
3290 mut_printf("hi %2d\n",v);
|
|
|
3291 hi[v] = solve_secant(gldn_like,tlk,fp->coef[v],hi[v],1e-8,BDF_EXPRIGHT,&err);
|
|
|
3292 if (err>0) mut_printf("solve_secant error\n");
|
|
|
3293 }
|
|
|
3294
|
|
|
3295 int procvscb2(des,lf,v)
|
|
|
3296 design *des;
|
|
|
3297 lfit *lf;
|
|
|
3298 int v;
|
|
|
3299 { double thhat, sd, *lo, *hi, u;
|
|
|
3300 int err, st, tmp;
|
|
|
3301 x = des->xev = evpt(&lf->fp,v);
|
|
|
3302 tmp = haspc(&lf->pc);
|
|
|
3303 /* if ((ker(&lf->sp)==WPARM) && (haspc(&lf->pc)))
|
|
|
3304 { lf->coef[v] = thhat = addparcomp(lf,des->xev,PCOEF);
|
|
|
3305 lf->nlx[v] = lf->t0[v] = sd = addparcomp(lf,des->xev,PNLX);
|
|
|
3306 }
|
|
|
3307 else */
|
|
|
3308 { haspc(&lf->pc) = 0;
|
|
|
3309 st = procvstd(des,lf,v);
|
|
|
3310 thhat = lf->fp.coef[v];
|
|
|
3311 sd = lf->fp.nlx[v];
|
|
|
3312 }
|
|
|
3313 if ((type==GLM2) | (type==GLM3) | (type==GLM4))
|
|
|
3314 { if (ker(&lf->sp) != WPARM)
|
|
|
3315 WARN(("nonparametric fit; correction is invalid"));
|
|
|
3316 cumulant(lf,des,sd);
|
|
|
3317 }
|
|
|
3318 haspc(&lf->pc) = tmp;
|
|
|
3319 lo = lf->fp.t0;
|
|
|
3320 hi = &lo[lf->fp.nvm];
|
|
|
3321 switch(type)
|
|
|
3322 {
|
|
|
3323 case GLM1:
|
|
|
3324 return(st);
|
|
|
3325 case GLM2: /* centered scr */
|
|
|
3326 lo[v] = kap[1];
|
|
|
3327 hi[v] = sqrt(kap[2]);
|
|
|
3328 return(st);
|
|
|
3329 case GLM3: /* corrected 2 */
|
|
|
3330 lo[v] = solve_secant(q2,scb_crit,0.0,2*scb_crit,0.000001,BDF_NONE,&err);
|
|
|
3331 return(st);
|
|
|
3332 case GLM4: /* corrected 2' */
|
|
|
3333 u = fabs(p2(scb_crit));
|
|
|
3334 max_p2 = MAX(max_p2,u);
|
|
|
3335 return(st);
|
|
|
3336 case GLDN:
|
|
|
3337 get_gldn(&lf->fp,des,lo,hi,v);
|
|
|
3338 return(st);
|
|
|
3339 }
|
|
|
3340 LERR(("procvscb2: invalid type"));
|
|
|
3341 return(st);
|
|
|
3342 }
|
|
|
3343
|
|
|
3344 void scb(lf,des,res)
|
|
|
3345 lfit *lf;
|
|
|
3346 design *des;
|
|
|
3347 double *res;
|
|
|
3348 { double k1, k2, *lo, *hi, sig, thhat, nlx, rss[10];
|
|
|
3349 int i, nterms;
|
|
|
3350
|
|
|
3351 scb_des= des;
|
|
|
3352
|
|
|
3353 npar(&lf->sp) = calcp(&lf->sp,lf->lfd.d);
|
|
|
3354 des_init(des,lf->lfd.n,npar(&lf->sp));
|
|
|
3355
|
|
|
3356 type = geth(&lf->fp);
|
|
|
3357
|
|
|
3358 if (type >= 80) /* simultaneous */
|
|
|
3359 {
|
|
|
3360 nterms = constants(lf,des,res);
|
|
|
3361 scb_crit = critval(0.05,res,nterms,lf->lfd.d,TWO_SIDED,0.0,GAUSS);
|
|
|
3362 type -= 10;
|
|
|
3363 }
|
|
|
3364 else /* pointwise */
|
|
|
3365 { res[0] = 1;
|
|
|
3366 scb_crit = critval(0.05,res,1,lf->lfd.d,TWO_SIDED,0.0,GAUSS);
|
|
|
3367 }
|
|
|
3368
|
|
|
3369 max_p2 = 0.0;
|
|
|
3370 lf->mdl.procv = procvscb2;
|
|
|
3371 nstartlf(des,lf);
|
|
|
3372
|
|
|
3373 if ((fam(&lf->sp)&64)==64)
|
|
|
3374 { i = haspc(&lf->pc); haspc(&lf->pc) = 0;
|
|
|
3375 ressumm(lf,des,rss);
|
|
|
3376 haspc(&lf->pc) = i;
|
|
|
3377 sig = sqrt(rss[3]);
|
|
|
3378 }
|
|
|
3379 else sig = 1.0;
|
|
|
3380
|
|
|
3381 lo = lf->fp.t0;
|
|
|
3382 hi = &lo[lf->fp.nvm];
|
|
|
3383 for (i=0; i<lf->fp.nv; i++)
|
|
|
3384 { thhat = lf->fp.coef[i];
|
|
|
3385 nlx = lf->fp.nlx[i];
|
|
|
3386 switch(type)
|
|
|
3387 {
|
|
|
3388 case GLM1: /* basic scb */
|
|
|
3389 lo[i] = thhat - scb_crit * sig * nlx;
|
|
|
3390 hi[i] = thhat + scb_crit * sig * nlx;
|
|
|
3391 break;
|
|
|
3392 case GLM2:
|
|
|
3393 k1 = lo[i];
|
|
|
3394 k2 = hi[i];
|
|
|
3395 lo[i] = thhat - k1*nlx - scb_crit*nlx*k2;
|
|
|
3396 hi[i] = thhat - k1*nlx + scb_crit*nlx*k2;
|
|
|
3397 break;
|
|
|
3398 case GLM3:
|
|
|
3399 k1 = lo[i];
|
|
|
3400 lo[i] = thhat - k1*nlx;
|
|
|
3401 hi[i] = thhat + k1*nlx;
|
|
|
3402 case GLM4: /* corrected 2' */
|
|
|
3403 lo[i] = thhat - (scb_crit-max_p2)*lf->fp.nlx[i];
|
|
|
3404 hi[i] = thhat + (scb_crit-max_p2)*lf->fp.nlx[i];
|
|
|
3405 break;
|
|
|
3406 case GLDN:
|
|
|
3407 break;
|
|
|
3408 }
|
|
|
3409 }
|
|
|
3410 }
|
|
|
3411
|
|
|
3412 void initscb(lf)
|
|
|
3413 lfit *lf;
|
|
|
3414 { initstd(lf);
|
|
|
3415 PROCV(lf) = NULL;
|
|
|
3416 KEEPC(lf) = NVAR(lf)+1+k0_reqd(NVAR(lf),NOBS(lf),0);
|
|
|
3417 PPROC(lf) = scb;
|
|
|
3418 }
|
|
|
3419 /*
|
|
|
3420 * Copyright 1996-2006 Catherine Loader.
|
|
|
3421 */
|
|
|
3422 #include "lfev.h"
|
|
|
3423
|
|
|
3424 int procvsimple(des,lf,v)
|
|
|
3425 design *des;
|
|
|
3426 lfit *lf;
|
|
|
3427 int v;
|
|
|
3428 { int lf_status;
|
|
|
3429 lf_status = procv_nov(des,lf,v);
|
|
|
3430 VVAL(lf,v,0) = des->cf[cfn(des,0)];
|
|
|
3431 return(lf_status);
|
|
|
3432 }
|
|
|
3433
|
|
|
3434 void allocsimple(lf)
|
|
|
3435 lfit *lf;
|
|
|
3436 { lf->fp.coef = VVEC(lf,0);
|
|
|
3437 lf->fp.h = NULL;
|
|
|
3438 }
|
|
|
3439
|
|
|
3440 void initsimple(lf)
|
|
|
3441 lfit *lf;
|
|
|
3442 {
|
|
|
3443 PROCV(lf) = procvsimple;
|
|
|
3444 ALLOC(lf) = allocsimple;
|
|
|
3445 PPROC(lf) = NULL;
|
|
|
3446 KEEPV(lf) = 1;
|
|
|
3447 KEEPC(lf) = 1;
|
|
|
3448 NOPC(lf) = 1;
|
|
|
3449 }
|
|
|
3450 /*
|
|
|
3451 * Copyright 1996-2006 Catherine Loader.
|
|
|
3452 */
|
|
|
3453 #include "lfev.h"
|
|
|
3454
|
|
|
3455 /* 3d+8 variables to keep:
|
|
|
3456 * d+1 coef+derivs.
|
|
|
3457 * d+1 sd's + derivs.
|
|
|
3458 * d+1 infl + derivs.
|
|
|
3459 * 3 likelihood and d.f's.
|
|
|
3460 * 1 bandwidth h
|
|
|
3461 * 1 degree.
|
|
|
3462 */
|
|
|
3463
|
|
|
3464 void allocstd(lf)
|
|
|
3465 lfit *lf;
|
|
|
3466 { int d;
|
|
|
3467 d = NVAR(lf);
|
|
|
3468 lf->fp.coef = VVEC(lf,0);
|
|
|
3469 lf->fp.nlx = VVEC(lf,d+1);
|
|
|
3470 lf->fp.t0 = VVEC(lf,2*d+2);
|
|
|
3471 lf->fp.lik = VVEC(lf,3*d+3);
|
|
|
3472 lf->fp.h = VVEC(lf,3*d+6);
|
|
|
3473 lf->fp.deg = VVEC(lf,3*d+7);
|
|
|
3474 }
|
|
|
3475
|
|
|
3476 int procvstd(des,lf,v)
|
|
|
3477 design *des;
|
|
|
3478 lfit *lf;
|
|
|
3479 int v;
|
|
|
3480 { int d, p, nvm, i, k;
|
|
|
3481 double t0[1+MXDIM], vari[1+MXDIM];
|
|
|
3482 k = procv_var(des,lf,v);
|
|
|
3483 if (lf_error) return(k);
|
|
|
3484
|
|
|
3485 d = lf->lfd.d;
|
|
|
3486 p = npar(&lf->sp);
|
|
|
3487 nvm = lf->fp.nvm;
|
|
|
3488
|
|
|
3489 if (k != LF_OK) lf_status_msg(k);
|
|
|
3490
|
|
|
3491 lf->fp.lik[v] = des->llk;
|
|
|
3492 lf->fp.lik[nvm+v] = des->tr2;
|
|
|
3493 lf->fp.lik[2*nvm+v] = des->tr0 - des->tr2;
|
|
|
3494
|
|
|
3495 for (i=0; i<des->ncoef; i++)
|
|
|
3496 vari[i] = des->V[p*cfn(des,0) + cfn(des,i)];
|
|
|
3497 vari[0] = sqrt(vari[0]);
|
|
|
3498 if (vari[0]>0) for (i=1; i<des->ncoef; i++) vari[i] /= vari[0];
|
|
|
3499
|
|
|
3500 t0[0] = sqrt(des->f1[0]);
|
|
|
3501 if (t0[0]>0) for (i=1; i<des->ncoef; i++) t0[i] = des->f1[i]/t0[0];
|
|
|
3502
|
|
|
3503 if (dc(&lf->fp)) dercor(&lf->lfd,&lf->sp,des,des->cf);
|
|
|
3504 subparcomp(des,lf,des->cf);
|
|
|
3505 for (i=0; i<des->ncoef; i++)
|
|
|
3506 lf->fp.coef[i*lf->fp.nvm+v] = des->cf[cfn(des,i)];
|
|
|
3507
|
|
|
3508 subparcomp2(des,lf,vari,t0);
|
|
|
3509 for (i=0; i<des->ncoef; i++)
|
|
|
3510 { lf->fp.nlx[i*nvm+v] = vari[i];
|
|
|
3511 lf->fp.t0[i*nvm+v] = t0[i];
|
|
|
3512 }
|
|
|
3513
|
|
|
3514 lf->fp.deg[v] = deg(&lf->sp);
|
|
|
3515 return(k);
|
|
|
3516 }
|
|
|
3517
|
|
|
3518 void pprocstd(lf,des,res)
|
|
|
3519 lfit *lf;
|
|
|
3520 design *des;
|
|
|
3521 double *res;
|
|
|
3522 {
|
|
|
3523 ressumm(lf,des,res);
|
|
|
3524 }
|
|
|
3525
|
|
|
3526 void initstd(lf)
|
|
|
3527 lfit *lf;
|
|
|
3528 { PROCV(lf) = procvstd;
|
|
|
3529 ALLOC(lf) = allocstd;
|
|
|
3530 PPROC(lf) = pprocstd;
|
|
|
3531 KEEPV(lf) = 3*NVAR(lf) + 8;
|
|
|
3532 KEEPC(lf) = 6;
|
|
|
3533 NOPC(lf) = 0;
|
|
|
3534 }
|
|
|
3535 /*
|
|
|
3536 * Copyright 1996-2006 Catherine Loader.
|
|
|
3537 */
|
|
|
3538 #include "lfev.h"
|
|
|
3539
|
|
|
3540 extern void procstd(), allocstd();
|
|
|
3541 extern double robscale;
|
|
|
3542
|
|
|
3543 double vocri(lk,t0,t2,pen)
|
|
|
3544 double lk, t0, t2, pen;
|
|
|
3545 { if (pen==0) return(-2*t0*lk/((t0-t2)*(t0-t2)));
|
|
|
3546 return((-2*lk+pen*t2)/t0);
|
|
|
3547 }
|
|
|
3548
|
|
|
3549 double intvo(des,lf,c0,c1,a,p,t0,t20,t21)
|
|
|
3550 design *des;
|
|
|
3551 lfit *lf;
|
|
|
3552 double *c0, *c1, a, t0, t20, t21;
|
|
|
3553 int p;
|
|
|
3554 { double th, lk, link[LLEN];
|
|
|
3555 int i, ii;
|
|
|
3556 lk = 0;
|
|
|
3557 for (i=0; i<des->n; i++)
|
|
|
3558 { ii = des->ind[i];
|
|
|
3559 th = (1-a)*innerprod(c0,d_xi(des,i),p) + a*innerprod(c1,d_xi(des,i),p);
|
|
|
3560 stdlinks(link,&lf->lfd,&lf->sp,ii,th,robscale);
|
|
|
3561 lk += wght(des,ii)*link[ZLIK];
|
|
|
3562 }
|
|
|
3563 des->llk = lk;
|
|
|
3564 return(vocri(des->llk,t0,(1-a)*t20+a*t21,pen(&lf->sp)));
|
|
|
3565 }
|
|
|
3566
|
|
|
3567 int procvvord(des,lf,v)
|
|
|
3568 design *des;
|
|
|
3569 lfit *lf;
|
|
|
3570 int v;
|
|
|
3571 { double tr[6], gcv, g0, ap, coef[4][10], t2[4], th, md;
|
|
|
3572 int i, j, k, d1, i0, p1, ip;
|
|
|
3573 des->xev = evpt(&lf->fp,v);
|
|
|
3574
|
|
|
3575 ap = pen(&lf->sp);
|
|
|
3576 if ((ap==0) & ((fam(&lf->sp)&63)!=TGAUS)) ap = 2.0;
|
|
|
3577 d1 = deg(&lf->sp); p1 = npar(&lf->sp);
|
|
|
3578 for (i=0; i<p1; i++) coef[0][i] = coef[1][i] = coef[2][i] = coef[3][i] = 0.0;
|
|
|
3579 i0 = 0; g0 = 0;
|
|
|
3580 ip = 1;
|
|
|
3581
|
|
|
3582 for (i=deg0(&lf->sp); i<=d1; i++)
|
|
|
3583 { deg(&lf->sp) = i;
|
|
|
3584 des->p = npar(&lf->sp) = calcp(&lf->sp,lf->lfd.d);
|
|
|
3585 k = locfit(&lf->lfd,des,&lf->sp,0, i==deg0(&lf->sp),0);
|
|
|
3586
|
|
|
3587 local_df(&lf->lfd,&lf->sp,des,tr);
|
|
|
3588 gcv = vocri(des->llk,tr[0],tr[2],ap);
|
|
|
3589 if ((i==deg0(&lf->sp)) || (gcv<g0)) { i0 = i; g0 = gcv; md = i; }
|
|
|
3590
|
|
|
3591 for (j=0; j<des->p; j++) coef[i][j] = des->cf[j];
|
|
|
3592 t2[i] = tr[2];
|
|
|
3593
|
|
|
3594 #ifdef RESEARCH
|
|
|
3595 if ((ip) && (i>deg0(&lf->sp)))
|
|
|
3596 { for (j=1; j<10; j++)
|
|
|
3597 { gcv = intvo(des,lf,coef[i-1],coef[i],j/10.0,des->p,tr[0],t2[i-1],t2[i]);
|
|
|
3598 if (gcv<g0) { g0 = gcv; md = i-1+j/10.0; }
|
|
|
3599 }
|
|
|
3600 }
|
|
|
3601 #endif
|
|
|
3602 }
|
|
|
3603 lf->fp.h[v] = des->h;
|
|
|
3604 if (lf->fp.h[v]<=0) WARN(("zero bandwidth in procvvord"));
|
|
|
3605
|
|
|
3606 if (i0<d1) /* recompute the best fit */
|
|
|
3607 { deg(&lf->sp) = i0;
|
|
|
3608 des->p = npar(&lf->sp) = calcp(&lf->sp,lf->lfd.d);
|
|
|
3609 k = locfit(&lf->lfd,des,&lf->sp,0,0,0);
|
|
|
3610 for (i=npar(&lf->sp); i<p1; i++) des->cf[i] = 0.0;
|
|
|
3611 i0 = md; if (i0==d1) i0--;
|
|
|
3612 th = md-i0;
|
|
|
3613 for (i=0; i<p1; i++) des->cf[i] = (1-th)*coef[i0][i]+th*coef[i0+1][i];
|
|
|
3614 deg(&lf->sp) = d1; npar(&lf->sp) = p1;
|
|
|
3615 }
|
|
|
3616
|
|
|
3617 for (i=0; i<p1; i++) lf->fp.coef[i*lf->fp.nvm+v] = des->cf[i];
|
|
|
3618 lf->fp.deg[v] = md;
|
|
|
3619 return(k);
|
|
|
3620 }
|
|
|
3621
|
|
|
3622 void initvord(lf)
|
|
|
3623 lfit *lf;
|
|
|
3624 { initstd(lf);
|
|
|
3625 PROCV(lf) = procvvord;
|
|
|
3626 ALLOC(lf) = allocstd;
|
|
|
3627 PPROC(lf) = NULL;
|
|
|
3628 KEEPC(lf) = 0;
|
|
|
3629 NOPC(lf) = 1;
|
|
|
3630 }
|
|
|
3631 /*
|
|
|
3632 * Copyright 1996-2006 Catherine Loader.
|
|
|
3633 */
|
|
|
3634 /*
|
|
|
3635 * functions for computing and subtracting, adding the
|
|
|
3636 * parametric component
|
|
|
3637 */
|
|
|
3638
|
|
|
3639 #include "lfev.h"
|
|
|
3640
|
|
|
3641 int noparcomp(sp)
|
|
|
3642 smpar *sp;
|
|
|
3643 { int tg;
|
|
|
3644 if (ubas(sp)) return(1);
|
|
|
3645 tg = fam(sp) & 63;
|
|
|
3646 if (tg<=THAZ) return(1);
|
|
|
3647 if (tg==TROBT) return(1);
|
|
|
3648 if (tg==TCAUC) return(1);
|
|
|
3649 if (tg==TQUANT) return(1);
|
|
|
3650 return(0);
|
|
|
3651 }
|
|
|
3652
|
|
|
3653 int pc_reqd(d,p)
|
|
|
3654 int d, p;
|
|
|
3655 { return(d + 2*p + jac_reqd(p));
|
|
|
3656 }
|
|
|
3657
|
|
|
3658 void pcchk(pc,d,p,lc)
|
|
|
3659 paramcomp *pc;
|
|
|
3660 int d, p, lc;
|
|
|
3661 { int rw;
|
|
|
3662 double *z;
|
|
|
3663
|
|
|
3664 rw = pc_reqd(d,p);
|
|
|
3665 if (pc->lwk < rw)
|
|
|
3666 { pc->wk = (double *)calloc(rw,sizeof(double));
|
|
|
3667 if ( pc->wk == NULL ) {
|
|
|
3668 printf("Problem allocating memory for pc->wk\n");fflush(stdout);
|
|
|
3669 }
|
|
|
3670 pc->lwk= rw;
|
|
|
3671 }
|
|
|
3672 z = pc->wk;
|
|
|
3673
|
|
|
3674 pc->xbar = z; z += d;
|
|
|
3675 pc->coef = z; z += p;
|
|
|
3676 pc->f = z; z += p;
|
|
|
3677
|
|
|
3678 z = jac_alloc(&pc->xtwx,p,z);
|
|
|
3679 pc->xtwx.p = p;
|
|
|
3680 }
|
|
|
3681
|
|
|
3682 void compparcomp(des,lfd,sp,pc,nopc)
|
|
|
3683 design *des;
|
|
|
3684 lfdata *lfd;
|
|
|
3685 smpar *sp;
|
|
|
3686 paramcomp *pc;
|
|
|
3687 int nopc;
|
|
|
3688 { int i, j, k, p;
|
|
|
3689 double wt, sw;
|
|
|
3690
|
|
|
3691 if (lf_debug>1) mut_printf(" compparcomp:\n");
|
|
|
3692 p = des->p;
|
|
|
3693 pcchk(pc,lfd->d,p,1);
|
|
|
3694
|
|
|
3695 for (i=0; i<lfd->d; i++) pc->xbar[i] = 0.0;
|
|
|
3696 sw = 0.0;
|
|
|
3697 for (i=0; i<lfd->n; i++)
|
|
|
3698 {
|
|
|
3699 wt = prwt(lfd,i);
|
|
|
3700 sw += wt;
|
|
|
3701 for (j=0; j<lfd->d; j++)
|
|
|
3702 pc->xbar[j] += datum(lfd,j,i)*wt;
|
|
|
3703 des->ind[i] = i;
|
|
|
3704 wght(des,i) = 1.0;
|
|
|
3705 }
|
|
|
3706 for (i=0; i<lfd->d; i++) pc->xbar[i] /= sw;
|
|
|
3707 if ((nopc) || noparcomp(sp))
|
|
|
3708 { haspc(pc) = 0;
|
|
|
3709 return;
|
|
|
3710 }
|
|
|
3711 haspc(pc) = 1;
|
|
|
3712 des->xev = pc->xbar;
|
|
|
3713 k = locfit(lfd,des,sp,0,0,0);
|
|
|
3714 if (k != LF_OK) lf_status_msg(k);
|
|
|
3715 if (lf_error) return;
|
|
|
3716 switch(k)
|
|
|
3717 { case LF_NOPT: return;
|
|
|
3718 case LF_INFA: return;
|
|
|
3719 case LF_NCON: return;
|
|
|
3720 case LF_OOB: return;
|
|
|
3721 case LF_NSLN: return;
|
|
|
3722 case LF_PF:
|
|
|
3723 WARN(("compparcomp: perfect fit"));
|
|
|
3724 case LF_OK:
|
|
|
3725 case LF_DONE:
|
|
|
3726 for (i=0; i<p; i++)
|
|
|
3727 { pc->coef[i] = des->cf[i];
|
|
|
3728 pc->xtwx.dg[i] = des->xtwx.dg[i];
|
|
|
3729 pc->xtwx.wk[i] = des->xtwx.wk[i];
|
|
|
3730 }
|
|
|
3731 for (i=0; i<p*p; i++)
|
|
|
3732 { pc->xtwx.Z[i] = des->xtwx.Z[i];
|
|
|
3733 pc->xtwx.Q[i] = des->xtwx.Q[i];
|
|
|
3734 }
|
|
|
3735 pc->xtwx.sm = des->xtwx.sm;
|
|
|
3736 pc->xtwx.st = des->xtwx.st;
|
|
|
3737 return;
|
|
|
3738 default:
|
|
|
3739 LERR(("compparcomp: locfit unknown return status %d",k));
|
|
|
3740 return;
|
|
|
3741 }
|
|
|
3742 }
|
|
|
3743
|
|
|
3744 void subparcomp(des,lf,coef)
|
|
|
3745 design *des;
|
|
|
3746 lfit *lf;
|
|
|
3747 double *coef;
|
|
|
3748 { int i, nd;
|
|
|
3749 deriv *dv;
|
|
|
3750 paramcomp *pc;
|
|
|
3751
|
|
|
3752 pc = &lf->pc;
|
|
|
3753 if (!haspc(pc)) return;
|
|
|
3754
|
|
|
3755 dv = &lf->dv; nd = dv->nd;
|
|
|
3756 fitfun(&lf->lfd, &lf->sp, des->xev,pc->xbar,des->f1,dv);
|
|
|
3757 coef[0] -= innerprod(pc->coef,des->f1,pc->xtwx.p);
|
|
|
3758 if (des->ncoef == 1) return;
|
|
|
3759
|
|
|
3760 dv->nd = nd+1;
|
|
|
3761 for (i=0; i<lf->lfd.d; i++)
|
|
|
3762 { dv->deriv[nd] = i;
|
|
|
3763 fitfun(&lf->lfd, &lf->sp, des->xev,pc->xbar,des->f1,dv);
|
|
|
3764 coef[i+1] -= innerprod(pc->coef,des->f1,pc->xtwx.p);
|
|
|
3765 }
|
|
|
3766 dv->nd = nd;
|
|
|
3767 }
|
|
|
3768
|
|
|
3769 void subparcomp2(des,lf,vr,il)
|
|
|
3770 design *des;
|
|
|
3771 lfit *lf;
|
|
|
3772 double *vr, *il;
|
|
|
3773 { double t0, t1;
|
|
|
3774 int i, nd;
|
|
|
3775 deriv *dv;
|
|
|
3776 paramcomp *pc;
|
|
|
3777
|
|
|
3778 pc = &lf->pc;
|
|
|
3779 if (!haspc(pc)) return;
|
|
|
3780
|
|
|
3781 dv = &lf->dv; nd = dv->nd;
|
|
|
3782
|
|
|
3783 fitfun(&lf->lfd, &lf->sp, des->xev,pc->xbar,des->f1,dv);
|
|
|
3784 for (i=0; i<npar(&lf->sp); i++) pc->f[i] = des->f1[i];
|
|
|
3785 jacob_solve(&pc->xtwx,des->f1);
|
|
|
3786 t0 = sqrt(innerprod(pc->f,des->f1,pc->xtwx.p));
|
|
|
3787 vr[0] -= t0;
|
|
|
3788 il[0] -= t0;
|
|
|
3789 if ((t0==0) | (des->ncoef==1)) return;
|
|
|
3790
|
|
|
3791 dv->nd = nd+1;
|
|
|
3792 for (i=0; i<lf->lfd.d; i++)
|
|
|
3793 { dv->deriv[nd] = i;
|
|
|
3794 fitfun(&lf->lfd, &lf->sp, des->xev,pc->xbar,pc->f,dv);
|
|
|
3795 t1 = innerprod(pc->f,des->f1,pc->xtwx.p)/t0;
|
|
|
3796 vr[i+1] -= t1;
|
|
|
3797 il[i+1] -= t1;
|
|
|
3798 }
|
|
|
3799 dv->nd = nd;
|
|
|
3800 }
|
|
|
3801
|
|
|
3802 double addparcomp(lf,x,c)
|
|
|
3803 lfit *lf;
|
|
|
3804 double *x;
|
|
|
3805 int c;
|
|
|
3806 { double y;
|
|
|
3807 paramcomp *pc;
|
|
|
3808
|
|
|
3809 pc = &lf->pc;
|
|
|
3810 if (!haspc(pc)) return(0.0);
|
|
|
3811 fitfun(&lf->lfd, &lf->sp, x,pc->xbar,pc->f,&lf->dv);
|
|
|
3812 if (c==PCOEF) return(innerprod(pc->coef,pc->f,pc->xtwx.p));
|
|
|
3813 if ((c==PNLX)|(c==PT0)|(c==PVARI))
|
|
|
3814 { y = sqrt(jacob_qf(&pc->xtwx,pc->f));
|
|
|
3815 return(y);
|
|
|
3816 }
|
|
|
3817 return(0.0);
|
|
|
3818 }
|
|
|
3819 /*
|
|
|
3820 * Copyright 1996-2006 Catherine Loader.
|
|
|
3821 */
|
|
|
3822 #include "lfev.h"
|
|
|
3823
|
|
|
3824 /*
|
|
|
3825 preplot(): interpolates the fit to a new set of points.
|
|
|
3826 lf -- the fit structure.
|
|
|
3827 x -- the points to predict at.
|
|
|
3828 f -- vector to return the predictions.
|
|
|
3829 se -- vector to return std errors (NULL if not req'd)
|
|
|
3830 band-- char for conf band type. ('n'=none, 'g'=global etc.)
|
|
|
3831 n -- no of predictions (or vector of margin lengths for grid)
|
|
|
3832 where -- where to predict:
|
|
|
3833 1 = points in the array x.
|
|
|
3834 2 = grid defined by margins in x.
|
|
|
3835 3 = data points from lf (ignore x).
|
|
|
3836 4 = fit points from lf (ignore x).
|
|
|
3837 what -- what to predict.
|
|
|
3838 (PCOEF etc; see lfcons.h file)
|
|
|
3839
|
|
|
3840 */
|
|
|
3841
|
|
|
3842 #define NWH 8
|
|
|
3843 static char *whtyp[NWH] = { "coef", "nlx", "infl", "band",
|
|
|
3844 "degr", "like", "rdf", "vari" };
|
|
|
3845 static int whval[NWH] = { PCOEF, PNLX, PT0, PBAND, PDEGR, PLIK, PRDF, PVARI };
|
|
|
3846 int ppwhat(z)
|
|
|
3847 char *z;
|
|
|
3848 { int val;
|
|
|
3849
|
|
|
3850 val = pmatch(z, whtyp, whval, NWH, -1);
|
|
|
3851 if (val==-1) LERR(("Unknown what = %s",z));
|
|
|
3852 return(val);
|
|
|
3853 }
|
|
|
3854
|
|
|
3855 static char cb;
|
|
|
3856 double *sef, *fit, sigmahat;
|
|
|
3857
|
|
|
3858 void predptall(lf,x,what,ev,i)
|
|
|
3859 lfit *lf;
|
|
|
3860 double *x;
|
|
|
3861 int what, ev, i;
|
|
|
3862 { double lik, rdf;
|
|
|
3863 fit[i] = dointpoint(lf,x,what,ev,i);
|
|
|
3864 if (cb=='n') return;
|
|
|
3865 sef[i] = dointpoint(lf,x,PNLX,ev,i);
|
|
|
3866 if (cb=='g')
|
|
|
3867 { sef[i] *= sigmahat;
|
|
|
3868 return;
|
|
|
3869 }
|
|
|
3870 if (cb=='l')
|
|
|
3871 { lik = dointpoint(lf,x,PLIK,ev,i);
|
|
|
3872 rdf = dointpoint(lf,x,PRDF,ev,i);
|
|
|
3873 sef[i] *= sqrt(-2*lik/rdf);
|
|
|
3874 return;
|
|
|
3875 }
|
|
|
3876 if (cb=='p')
|
|
|
3877 { sef[i] = sigmahat*sqrt(1+sef[i]*sef[i]);
|
|
|
3878 return;
|
|
|
3879 }
|
|
|
3880 }
|
|
|
3881
|
|
|
3882 void predptdir(lf,des,x,what,i)
|
|
|
3883 lfit *lf;
|
|
|
3884 design *des;
|
|
|
3885 double *x;
|
|
|
3886 int what, i;
|
|
|
3887 { int needcv;
|
|
|
3888 des->xev = x;
|
|
|
3889 needcv = (what==PVARI) | (what==PNLX) | (what==PT0) | (what==PRDF);
|
|
|
3890 locfit(&lf->lfd,des,&lf->sp,0,1,needcv);
|
|
|
3891 switch(what)
|
|
|
3892 { case PCOEF: fit[i] = des->cf[0]; break;
|
|
|
3893 case PVARI: fit[i] = des->V[0]; break;
|
|
|
3894 case PNLX: fit[i] = sqrt(des->V[0]); break;
|
|
|
3895 case PT0: fit[i] = des->f1[0]; break;
|
|
|
3896 case PBAND: fit[i] = des->h; break;
|
|
|
3897 case PDEGR: fit[i] = deg(&lf->sp); break;
|
|
|
3898 case PLIK: fit[i] = des->llk; break;
|
|
|
3899 case PRDF: fit[i] = des->tr0 - des->tr2; break;
|
|
|
3900 default: LERR(("unknown what in predptdir"));
|
|
|
3901 }
|
|
|
3902 }
|
|
|
3903
|
|
|
3904 void prepvector(lf,des,x,n,what,dir) /* interpolate a vector */
|
|
|
3905 lfit *lf;
|
|
|
3906 design *des;
|
|
|
3907 double **x;
|
|
|
3908 int n, what, dir;
|
|
|
3909 { int i, j;
|
|
|
3910 double xx[MXDIM];
|
|
|
3911 for (i=0; i<n; i++)
|
|
|
3912 { for (j=0; j<lf->fp.d; j++) xx[j] = x[j][i];
|
|
|
3913 if (dir)
|
|
|
3914 predptdir(lf,des,xx,what,i);
|
|
|
3915 else
|
|
|
3916 predptall(lf,xx,what,ev(&lf->evs),i);
|
|
|
3917 if (lf_error) return;
|
|
|
3918 }
|
|
|
3919 }
|
|
|
3920
|
|
|
3921 void prepfitp(lf,what)
|
|
|
3922 lfit *lf;
|
|
|
3923 int what;
|
|
|
3924 { int i;
|
|
|
3925 for (i=0; i<lf->fp.nv; i++)
|
|
|
3926 { predptall(lf,evpt(&lf->fp,i),what,EFITP,i);
|
|
|
3927 if (lf_error) return;
|
|
|
3928 }
|
|
|
3929 }
|
|
|
3930
|
|
|
3931 void prepgrid(lf,des,x,mg,n,what,dir) /* interpolate a grid given margins */
|
|
|
3932 lfit *lf;
|
|
|
3933 design *des;
|
|
|
3934 double **x;
|
|
|
3935 int *mg, dir, n, what;
|
|
|
3936 { int i, ii, j, d;
|
|
|
3937 double xv[MXDIM];
|
|
|
3938 d = lf->fp.d;
|
|
|
3939 for (i=0; i<n; i++)
|
|
|
3940 { ii = i;
|
|
|
3941 for (j=0; j<d; j++)
|
|
|
3942 { xv[j] = x[j][ii%mg[j]];
|
|
|
3943 ii /= mg[j];
|
|
|
3944 }
|
|
|
3945 if (dir)
|
|
|
3946 predptdir(lf,des,xv,what,i);
|
|
|
3947 else
|
|
|
3948 predptall(lf,xv,what,ev(&lf->evs),i);
|
|
|
3949 if (lf_error) return;
|
|
|
3950 }
|
|
|
3951 }
|
|
|
3952
|
|
|
3953 void preplot(lf,x,f,se,band,mg,where,what,dir)
|
|
|
3954 lfit *lf;
|
|
|
3955 double **x, *f, *se;
|
|
|
3956 int *mg, where, what, dir;
|
|
|
3957 char band;
|
|
|
3958 { int d, i, n;
|
|
|
3959 double *xx[MXDIM];
|
|
|
3960 design ppdes;
|
|
|
3961 d = lf->fp.d;
|
|
|
3962 fit = f;
|
|
|
3963 sef = se;
|
|
|
3964 cb = band;
|
|
|
3965 if (cb!='n') sigmahat = sqrt(rv(&lf->fp));
|
|
|
3966 if (dir) des_init(&ppdes,lf->lfd.n,npar(&lf->sp));
|
|
|
3967
|
|
|
3968 switch(where)
|
|
|
3969 { case 1: /* vector */
|
|
|
3970 n = mg[0];
|
|
|
3971 prepvector(lf,&ppdes,x,n,what,dir);
|
|
|
3972 break;
|
|
|
3973 case 2: /* grid */
|
|
|
3974 n = 1;
|
|
|
3975 for (i=0; i<d; i++) n *= mg[i];
|
|
|
3976 prepgrid(lf,&ppdes,x,mg,n,what,dir);
|
|
|
3977 break;
|
|
|
3978 case 3: /* data */
|
|
|
3979 n = lf->lfd.n;
|
|
|
3980 if ((ev(&lf->evs)==EDATA) | (ev(&lf->evs)==ECROS))
|
|
|
3981 { prepfitp(lf,what);
|
|
|
3982 dir = 0;
|
|
|
3983 }
|
|
|
3984 else
|
|
|
3985 { for (i=0; i<d; i++) xx[i] = dvari(&lf->lfd,i);
|
|
|
3986 prepvector(lf,&ppdes,xx,n,what,dir);
|
|
|
3987 }
|
|
|
3988 break;
|
|
|
3989 case 4: /* fit points */
|
|
|
3990 n = lf->fp.nv; dir = 0;
|
|
|
3991 prepfitp(lf,what);
|
|
|
3992 break;
|
|
|
3993 default:
|
|
|
3994 LERR(("unknown where in preplot"));
|
|
|
3995 }
|
|
|
3996
|
|
|
3997 if ((!dir) && ((what==PT0)|(what==PVARI)))
|
|
|
3998 for (i=0; i<n; i++) f[i] = f[i]*f[i];
|
|
|
3999 }
|
|
|
4000 /*
|
|
|
4001 * Copyright 1996-2006 Catherine Loader.
|
|
|
4002 */
|
|
|
4003 #include "lfev.h"
|
|
|
4004
|
|
|
4005 int procv_nov(des,lf,v)
|
|
|
4006 design *des;
|
|
|
4007 lfit *lf;
|
|
|
4008 int v;
|
|
|
4009 { int lf_status;
|
|
|
4010
|
|
|
4011 if (lf_debug>1) mut_printf(" procveraw: %d\n",v);
|
|
|
4012 des->xev = evpt(&lf->fp,v);
|
|
|
4013
|
|
|
4014 if (acri(&lf->sp)==ANONE)
|
|
|
4015 lf_status = locfit(&lf->lfd,des,&lf->sp,0,1,0);
|
|
|
4016 else
|
|
|
4017 lf_status = alocfit(&lf->lfd,&lf->sp,&lf->dv,des,0);
|
|
|
4018 if (lf->fp.h != NULL) lf->fp.h[v] = des->h;
|
|
|
4019
|
|
|
4020 return(lf_status);
|
|
|
4021 }
|
|
|
4022
|
|
|
4023 int procv_var(des,lf,v)
|
|
|
4024 design *des;
|
|
|
4025 lfit *lf;
|
|
|
4026 int v;
|
|
|
4027 { int i, lf_status;
|
|
|
4028
|
|
|
4029 if (lf_debug>1) mut_printf(" procvraw: %d\n",v);
|
|
|
4030 des->xev = evpt(&lf->fp,v);
|
|
|
4031
|
|
|
4032 if (acri(&lf->sp)==ANONE)
|
|
|
4033 lf_status = locfit(&lf->lfd,des,&lf->sp,0,1,1);
|
|
|
4034 else
|
|
|
4035 lf_status = alocfit(&lf->lfd,&lf->sp,&lf->dv,des,1);
|
|
|
4036 if (lf->fp.h != NULL) lf->fp.h[v] = des->h;
|
|
|
4037
|
|
|
4038 return(lf_status);
|
|
|
4039 }
|
|
|
4040 /*
|
|
|
4041 * Copyright 1996-2006 Catherine Loader.
|
|
|
4042 */
|
|
|
4043 /*
|
|
|
4044 * startmodule(lf,des,mod,dir) -- the standard entry point.
|
|
|
4045 * des and lf are pointers to the design and fit structures.
|
|
|
4046 * mod - module name. Set to NULL if the module is already
|
|
|
4047 * initialized.
|
|
|
4048 * dir - for dynamic modules, the directory.
|
|
|
4049 *
|
|
|
4050 * initmodule(mdl,mod,dir,lf)
|
|
|
4051 * direct call for module initialization.
|
|
|
4052 *
|
|
|
4053 */
|
|
|
4054
|
|
|
4055 #include "lfev.h"
|
|
|
4056
|
|
|
4057 #ifdef WINDOWS
|
|
|
4058
|
|
|
4059 #define DIRSEP '\\'
|
|
|
4060 #define PATHSEP ';'
|
|
|
4061
|
|
|
4062 #else
|
|
|
4063
|
|
|
4064 #define DIRSEP '/'
|
|
|
4065 #define PATHSEP ':'
|
|
|
4066
|
|
|
4067 #endif
|
|
|
4068
|
|
|
4069
|
|
|
4070 #ifdef ALLOW_MODULES
|
|
|
4071
|
|
|
4072 #ifdef WINDOWS
|
|
|
4073 #include <windows.h>
|
|
|
4074 #define DLEXT "dll"
|
|
|
4075 #define DLOPEN(x) LoadLibrary(x)
|
|
|
4076 #define DLSYM GetProcAddress
|
|
|
4077
|
|
|
4078 #else
|
|
|
4079
|
|
|
4080 #include <dlfcn.h>
|
|
|
4081 #define DLEXT "so"
|
|
|
4082 #define DLOPEN(x) dlopen(x,RTLD_LAZY)
|
|
|
4083 #define DLSYM dlsym
|
|
|
4084 #endif
|
|
|
4085
|
|
|
4086 #endif
|
|
|
4087
|
|
|
4088 static double fpkap[6];
|
|
|
4089 void fitpt_init(fp)
|
|
|
4090 fitpt *fp;
|
|
|
4091 {
|
|
|
4092 dc(fp) = 0;
|
|
|
4093 geth(fp) = GSTD;
|
|
|
4094 fp->nv = fp->nvm = 0;
|
|
|
4095 if (fp->kap==NULL) fp->kap = fpkap;
|
|
|
4096 }
|
|
|
4097
|
|
|
4098 void lfit_init(lf)
|
|
|
4099 lfit *lf;
|
|
|
4100 {
|
|
|
4101 lfdata_init(&lf->lfd);
|
|
|
4102 evstruc_init(&lf->evs);
|
|
|
4103 smpar_init(&lf->sp,&lf->lfd);
|
|
|
4104 deriv_init(&lf->dv);
|
|
|
4105 fitpt_init(&lf->fp);
|
|
|
4106 lf->mdl.np = 0;
|
|
|
4107 }
|
|
|
4108
|
|
|
4109 void fitdefault(lf)
|
|
|
4110 lfit *lf;
|
|
|
4111 { WARN(("fitdefault deprecated -- use lfit_init()"));
|
|
|
4112 lfit_init(lf);
|
|
|
4113 }
|
|
|
4114
|
|
|
4115 void set_flim(lfd,evs)
|
|
|
4116 lfdata *lfd;
|
|
|
4117 evstruc *evs;
|
|
|
4118 { int i, j, d, n;
|
|
|
4119 double z, mx, mn, *bx;
|
|
|
4120
|
|
|
4121 if (ev(evs)==ESPHR) return;
|
|
|
4122 d = lfd->d; n = lfd->n;
|
|
|
4123 bx = evs->fl;
|
|
|
4124 for (i=0; i<d; i++)
|
|
|
4125 if (bx[i]==bx[i+d])
|
|
|
4126 { if (lfd->sty[i]==STANGL)
|
|
|
4127 { bx[i] = 0.0; bx[i+d] = 2*PI*lfd->sca[i];
|
|
|
4128 }
|
|
|
4129 else
|
|
|
4130 { mx = mn = datum(lfd,i,0);
|
|
|
4131 for (j=1; j<n; j++)
|
|
|
4132 { mx = MAX(mx,datum(lfd,i,j));
|
|
|
4133 mn = MIN(mn,datum(lfd,i,j));
|
|
|
4134 }
|
|
|
4135 if (lfd->xl[i]<lfd->xl[i+d]) /* user set xlim; maybe use them. */
|
|
|
4136 { z = mx-mn;
|
|
|
4137 if (mn-0.2*z < lfd->xl[i]) mn = lfd->xl[i];
|
|
|
4138 if (mx+0.2*z > lfd->xl[i+d]) mx = lfd->xl[i+d];
|
|
|
4139 }
|
|
|
4140 bx[i] = mn;
|
|
|
4141 bx[i+d] = mx;
|
|
|
4142 }
|
|
|
4143 }
|
|
|
4144 }
|
|
|
4145
|
|
|
4146 double vvari(v,n)
|
|
|
4147 double *v;
|
|
|
4148 int n;
|
|
|
4149 { int i;
|
|
|
4150 double xb, s2;
|
|
|
4151 xb = s2 = 0.0;
|
|
|
4152 for (i=0; i<n; i++) xb += v[i];
|
|
|
4153 xb /= n;
|
|
|
4154 for (i=0; i<n; i++) s2 += SQR(v[i]-xb);
|
|
|
4155 return(s2/(n-1));
|
|
|
4156 }
|
|
|
4157
|
|
|
4158 void set_scales(lfd)
|
|
|
4159 lfdata *lfd;
|
|
|
4160 { int i;
|
|
|
4161 for (i=0; i<lfd->d; i++)
|
|
|
4162 if (lfd->sca[i]<=0) /* set automatic scales */
|
|
|
4163 { if (lfd->sty[i]==STANGL)
|
|
|
4164 lfd->sca[i] = 1.0;
|
|
|
4165 else lfd->sca[i] = sqrt(vvari(lfd->x[i],lfd->n));
|
|
|
4166 }
|
|
|
4167 }
|
|
|
4168
|
|
|
4169 void nstartlf(des,lf)
|
|
|
4170 design *des;
|
|
|
4171 lfit *lf;
|
|
|
4172 { int i, d, n;
|
|
|
4173
|
|
|
4174 if (lf_debug>0) mut_printf("nstartlf\n");
|
|
|
4175 n = lf->lfd.n;
|
|
|
4176 d = lf->lfd.d;
|
|
|
4177 npar(&lf->sp) = calcp(&lf->sp,d);
|
|
|
4178
|
|
|
4179 des_init(des,n,npar(&lf->sp));
|
|
|
4180 set_scales(&lf->lfd);
|
|
|
4181 set_flim(&lf->lfd,&lf->evs);
|
|
|
4182 compparcomp(des,&lf->lfd,&lf->sp,&lf->pc,lf->mdl.nopc);
|
|
|
4183 if (lf_error) return;
|
|
|
4184 makecfn(&lf->sp,des,&lf->dv,lf->lfd.d);
|
|
|
4185
|
|
|
4186 lf->lfd.ord = 0;
|
|
|
4187 if ((d==1) && (lf->lfd.sty[0]!=STANGL))
|
|
|
4188 { i = 1;
|
|
|
4189 while ((i<n) && (datum(&lf->lfd,0,i)>=datum(&lf->lfd,0,i-1))) i++;
|
|
|
4190 lf->lfd.ord = (i==n);
|
|
|
4191 }
|
|
|
4192 for (i=0; i<npar(&lf->sp); i++) des->fix[i] = 0;
|
|
|
4193
|
|
|
4194 lf->fp.d = lf->lfd.d;
|
|
|
4195 lf->fp.hasd = (des->ncoef==(1+lf->fp.d));
|
|
|
4196 lf->fp.nv = lf->evs.nce = 0;
|
|
|
4197
|
|
|
4198 if (lf_debug>1) mut_printf("call eval structure %d\n",ev(&lf->evs));
|
|
|
4199 switch(ev(&lf->evs))
|
|
|
4200 { case EPHULL: triang_start(des,lf); break;
|
|
|
4201 case EDATA: dataf(des,lf); break;
|
|
|
4202 case ECROS: crossf(des,lf); break;
|
|
|
4203 case EGRID: gridf(des,lf); break;
|
|
|
4204 case ETREE: atree_start(des,lf); break;
|
|
|
4205 case EKDCE: kt(&lf->sp) = KCE;
|
|
|
4206 case EKDTR: kdtre_start(des,lf); break;
|
|
|
4207 case EPRES: preset(des,lf); break;
|
|
|
4208 case EXBAR: xbarf(des,lf); break;
|
|
|
4209 case ENONE: return;
|
|
|
4210 case ESPHR: sphere_start(des,lf); break;
|
|
|
4211 case ESPEC: lf->evs.espec(des,lf); break;
|
|
|
4212 default: LERR(("startlf: Invalid evaluation structure %d",ev(&lf->evs)));
|
|
|
4213 }
|
|
|
4214
|
|
|
4215 /* renormalize for family=density */
|
|
|
4216 if ((de_renorm) && (fam(&lf->sp)==TDEN)) dens_renorm(lf,des);
|
|
|
4217 }
|
|
|
4218
|
|
|
4219 /*
|
|
|
4220 * getnextdir() gets the next dir from a string dirpath="dir1:dir2:dir3:..."
|
|
|
4221 * (;-separated on windows).
|
|
|
4222 * The directory is returned through dirnext, and the function returns
|
|
|
4223 * a pointer to the next string.
|
|
|
4224 * typical usage is recursive, dirpath = getnextdir(dirpath,dirnext).
|
|
|
4225 * with the above example, this sets dirnext="dir1" and dirpath="dir2:dir3:...".
|
|
|
4226 * if the input dirpath has no :, then it is copied to dirnext, and return is "".
|
|
|
4227 * if input dirpath is "", dirnext is set to "", and null pointer returned.
|
|
|
4228 */
|
|
|
4229 char *getnextdir(dirpath,dirnext)
|
|
|
4230 char *dirpath, *dirnext;
|
|
|
4231 { char *z;
|
|
|
4232 if (strlen(dirpath)==0)
|
|
|
4233 { sprintf(dirnext,"");
|
|
|
4234 return(NULL);
|
|
|
4235 }
|
|
|
4236
|
|
|
4237 z = strchr(dirpath,PATHSEP);
|
|
|
4238 if (z==NULL)
|
|
|
4239 { sprintf(dirnext,"%s%c",dirpath,DIRSEP);
|
|
|
4240 return(&dirpath[strlen(dirnext)]);
|
|
|
4241 }
|
|
|
4242
|
|
|
4243 *z = '\0';
|
|
|
4244 sprintf(dirnext,"%s%c",dirpath,DIRSEP);
|
|
|
4245 return(++z);
|
|
|
4246 }
|
|
|
4247
|
|
|
4248 int initmodule(mdl, mod, dir, lf)
|
|
|
4249 module *mdl;
|
|
|
4250 lfit *lf;
|
|
|
4251 char *mod, *dir;
|
|
|
4252 { int n, d, p;
|
|
|
4253 #ifdef ALLOW_MODULES
|
|
|
4254 #ifdef WINDOWS
|
|
|
4255 HINSTANCE res;
|
|
|
4256 typedef void (CALLBACK* DLLFN)();
|
|
|
4257 DLLFN init;
|
|
|
4258 #else
|
|
|
4259 void *res;
|
|
|
4260 void (*init)();
|
|
|
4261 #endif
|
|
|
4262 char distname[500];
|
|
|
4263 #endif
|
|
|
4264
|
|
|
4265 n = lf->lfd.n;
|
|
|
4266 d = lf->lfd.d;
|
|
|
4267 p = npar(&lf->sp) = calcp(&lf->sp,d);
|
|
|
4268
|
|
|
4269 mdl->isset = 1;
|
|
|
4270 PPROC(lf) = NULL;
|
|
|
4271 if (strcmp(mod,"std")==0) { initstd(lf); return(1); }
|
|
|
4272 if (strcmp(mod,"simple")==0) { initsimple(lf); return(1); }
|
|
|
4273 if (strcmp(mod,"allcf")==0) { initallcf(lf); return(1); }
|
|
|
4274 if (strcmp(mod,"hatm")==0) { inithatm(lf); return(1); }
|
|
|
4275 if (strcmp(mod,"kappa")==0) { initkappa(lf); return(1); }
|
|
|
4276 if (strcmp(mod,"lscv")==0) { initlscv(lf); return(1); }
|
|
|
4277 if (strcmp(mod,"gamf")==0) { initgam(lf); return(1); }
|
|
|
4278 if (strcmp(mod,"gamp")==0) { initgam(lf); return(1); }
|
|
|
4279 if (strcmp(mod,"rband")==0) { initrband(lf); return(1); }
|
|
|
4280 if (strcmp(mod,"scb")==0) { initscb(lf); return(1); }
|
|
|
4281 if (strcmp(mod,"vord")==0) { initvord(lf); return(1); }
|
|
|
4282
|
|
|
4283 #ifdef ALLOW_MODULES
|
|
|
4284 while (dir != NULL)
|
|
|
4285 {
|
|
|
4286 dir = getnextdir(dir,distname);
|
|
|
4287 sprintf(&distname[strlen(distname)],"mod%s.%s",mod,DLEXT);
|
|
|
4288 res = DLOPEN(distname);
|
|
|
4289 if (res==NULL)
|
|
|
4290 {
|
|
|
4291 #ifdef WINDOWS
|
|
|
4292 mut_printf("LoadLibrary failed: %s, %d\n",distname,GetLastError());
|
|
|
4293 #else
|
|
|
4294 mut_printf("dlopen failed: %s\n",dlerror());
|
|
|
4295 #endif
|
|
|
4296 }
|
|
|
4297 else
|
|
|
4298 {
|
|
|
4299 #ifdef WINDOWS
|
|
|
4300 mut_printf("LoadLibrary success: %s\n",distname);
|
|
|
4301 #else
|
|
|
4302 mut_printf("dlopen success: %s\n",distname);
|
|
|
4303 #endif
|
|
|
4304 sprintf(distname,"init%s",mod);
|
|
|
4305 init = (void *)DLSYM(res,distname);
|
|
|
4306 if (init==NULL)
|
|
|
4307 { mut_printf("I can't find %s() function.\n",distname);
|
|
|
4308 mdl->isset = 0;
|
|
|
4309 return(0);
|
|
|
4310 }
|
|
|
4311 init(lf);
|
|
|
4312 return(1);
|
|
|
4313 }
|
|
|
4314 }
|
|
|
4315 #endif
|
|
|
4316
|
|
|
4317 mdl->isset = 0;
|
|
|
4318 return(0);
|
|
|
4319 }
|
|
|
4320
|
|
|
4321 /*
|
|
|
4322 * startmodule is the entry point to launch the fit.
|
|
|
4323 * if mod is provided, will first initialize the module.
|
|
|
4324 * if mod==NULL, assumes the module has been initialized separately.
|
|
|
4325 */
|
|
|
4326 void startmodule(lf,des,mod,dir)
|
|
|
4327 lfit *lf;
|
|
|
4328 design *des;
|
|
|
4329 char *mod, *dir;
|
|
|
4330 { int z;
|
|
|
4331
|
|
|
4332 if (mod != NULL)
|
|
|
4333 { z = initmodule(&lf->mdl,mod,dir,lf);
|
|
|
4334 if (!z) return;
|
|
|
4335 }
|
|
|
4336
|
|
|
4337 lf->fp.nv = lf->evs.nce = 0;
|
|
|
4338 if (lf_error) return;
|
|
|
4339 if (PROCV(lf) != NULL) nstartlf(des,lf);
|
|
|
4340 if (lf_error) return;
|
|
|
4341 if (PPROC(lf) != NULL) PPROC(lf)(lf,des,lf->fp.kap);
|
|
|
4342 }
|
|
|
4343
|
|
|
4344 /* for compatability, more or less. */
|
|
|
4345 void startlf(des,lf,vfun,nopc)
|
|
|
4346 design *des;
|
|
|
4347 lfit *lf;
|
|
|
4348 int (*vfun)(), nopc;
|
|
|
4349 { int z;
|
|
|
4350 z = initmodule(&lf->mdl,"std",NULL,lf);
|
|
|
4351 if (!z) return;
|
|
|
4352 lf->mdl.procv = vfun;
|
|
|
4353 lf->mdl.nopc = nopc;
|
|
|
4354 nstartlf(des,lf);
|
|
|
4355 }
|
