1 /*
2 * Copyright (c) 1996-2001 Lucent Technologies.
3 * See README file for details.
4 *
5 *
6 *
7 startlf(des,lf,vfun,nopc) -- starting point for locfit.
8 des and lf are pointers to the design and fit structures.
9 vfun is the vertex processing function.
10 nopc=1 inhibits computation of parametric component.
11 fitdefault(lf,n,d) -- fit default parameters.
12 lf is pointer to fit; n and d are sample size and dimension.
13 deschk() -- assignment function for the design structure.
14 preproc() -- fit preprocessing (limits, scales, paramcomp etc.)
15 bbox() -- compute bounding box.
16
17 fitoptions()
18 clocfit() -- start point for CLocfit - interpret cmd line etc.
19 */
20
21 #include "local.h"
22
23 extern INT cvi;
24
fitdefault(lf,n,d)25 void fitdefault(lf,n,d)
26 lfit *lf;
27 INT n, d;
28 { INT i, *mi;
29 double *dp;
30
31 dp = lf->dp;
32 mi = lf->mi;
33
34 mi[MTG] = TNUL;
35 mi[MTG] = (lf->y==NULL) ? TDEN : (64+TGAUS);
36 mi[MLINK] = LDEFAU;
37 mi[MACRI] = ANONE;
38 mi[MDEG] = mi[MDEG0] = 2;
39 mi[MEV] = (ident==1) ? EDATA : ETREE;
40 mi[MKT] = KSPH; mi[MKER] = WTCUB;
41 mi[MIT] = IDEFA; mi[MDC] = mi[MREN] = 0;
42 mi[MK] = 100; mi[MMINT] = 20;
43 mi[MMXIT] = 20;
44 mi[MN] = n; mi[MDIM] = d;
45 mi[MDEB] = 0;
46 mi[MUBAS] = 0;
47
48
49 dp[DALP] = 0.7; dp[DFXH] = dp[DADP] = 0.0;
50 dp[DCUT] = 0.8;
51
52 if (d<=0)
53 ERROR(("must set MDIM before calling fitdefault"));
54 for (i=0; i<d; i++)
55 { lf->sca[i] = 1.0;
56 lf->xl[i] = lf->xl[i+d] = 0.0;
57 lf->fl[i] = lf->fl[i+d] = 0.0;
58 }
59 }
60
des_reqd(n,p)61 int des_reqd(n,p)
62 INT n, p;
63 {
64 return (n*(p+5)+2*p*p+4*p + jac_reqd(p));
65 }
des_reqi(INT n)66 int des_reqi(INT n) { return(n); }
67
deschk(des,n,p)68 void deschk(des,n,p)
69 design *des;
70 INT n, p;
71 {
72 double *z;
73 des->dw = checkvarlen(des->dw,des_reqd(n,p),"_deswork",VDOUBLE);
74 z = vdptr(des->dw);
75 des->X = z; z += n*p;
76 setzero(des->X, n*p);
77
78 des->w = z; z += n;
79 setzero(des->w, n);
80
81 des->res=z; z += n;
82 setzero(des->res, n);
83
84 des->di =z; z += n;
85 setzero(des->di, n);
86
87 des->th =z; z += n;
88 setzero(des->th, n);
89
90 des->wd =z; z += n;
91 setzero(des->wd, n);
92
93 des->V =z; z += p*p;
94 setzero(des->V, p*p);
95
96 des->P =z; z += p*p;
97 setzero(des->P, p*p);
98
99 des->f1 =z; z += p;
100 setzero(des->f1, p);
101
102 des->ss =z; z += p;
103 setzero(des->ss, p);
104
105 des->oc =z; z += p;
106 setzero(des->oc, p);
107
108 des->cf =z; z += p;
109 setzero(des->cf, p);
110
111 z = jac_alloc(&des->xtwx,p,z);
112
113 des->index = checkvarlen(des->index,des_reqi(n),"_desidx",VINT);
114 des->ind = (INT *)vdptr(des->index);
115 des->n = n;
116 des->p = p;
117 des->xtwx.p = p;
118 }
119
bbox(lf,bx)120 void bbox(lf,bx)
121 lfit *lf;
122 double *bx;
123 { INT i, j, d, n;
124 double z, mx, mn;
125 d = lf->mi[MDIM]; n = lf->mi[MN];
126 for (i=0; i<d; i++)
127 if (bx[i]==bx[i+d])
128 { if (lf->sty[i]==STANGL)
129 { bx[i] = 0.0; bx[i+d] = 2*PI*lf->sca[i];
130 }
131 else
132 { mx = mn = datum(lf,i,0);
133 for (j=1; j<n; j++)
134 { mx = MAX(mx,datum(lf,i,j));
135 mn = MIN(mn,datum(lf,i,j));
136 }
137 if (lf->xl[i]<lf->xl[i+d]) /* user set xlim; maybe use them. */
138 { z = mx-mn;
139 if (mn-0.2*z < lf->xl[i]) mn = lf->xl[i];
140 if (mx+0.2*z > lf->xl[i+d]) mx = lf->xl[i+d];
141 }
142 bx[i] = mn;
143 bx[i+d] = mx;
144 }
145 }
146 }
147
preproc(des,lf,nopc)148 void preproc(des,lf,nopc)
149 design *des;
150 lfit *lf;
151 INT nopc;
152 { INT d, i, j, n;
153 double xb;
154 d = lf->mi[MDIM]; n = lf->mi[MN];
155 lf->mi[MLINK] = defaultlink(lf->mi[MLINK],lf->mi[MTG]);
156 if (!validlinks(lf->mi[MLINK],lf->mi[MTG]))
157 { ERROR(("Invalid family/link combination"));
158 return;
159 }
160 compparcomp(des,lf,nopc);
161 if (lf->w==NULL)
162 lf->dp[DSWT] = lf->mi[MN];
163 else
164 { lf->dp[DSWT] = 0;
165 for (i=0; i<lf->mi[MN]; i++) lf->dp[DSWT] += prwt(lf,i);
166 }
167 for (i=0; i<d; i++)
168 if (lf->sca[i]<=0) /* set automatic scales */
169 { if (lf->sty[i]==STANGL) lf->sca[i] = 1.0;
170 else
171 { xb = lf->sca[i] = 0.0;
172 for (j=0; j<n; j++) xb += datum(lf,i,j);
173 xb /= n;
174 for (j=0; j<n; j++) lf->sca[i] += SQR(datum(lf,i,j)-xb);
175 lf->sca[i] = sqrt(lf->sca[i]/(n-1));
176 }
177 }
178 bbox(lf,lf->fl);
179 }
180
181 #ifdef CVERSION
182 extern void do_scbsim();
183 #endif
184
startlf(des,lf,vfun,nopc)185 void startlf(des,lf,vfun,nopc)
186 design *des;
187 lfit *lf;
188 INT (*vfun)(), nopc;
189 {
190 INT i, *mi;
191 des->vfun = vfun;
192 mi = lf->mi;
193 mi[MP] = calcp(mi,mi[MDEG]);
194 des->pref = 0;
195 cvi = -1; /* inhibit cross validation */
196 deschk(des,mi[MN],mi[MP]);
197 if (mi[MDEB]>0) printf("preprocess\n");
198 preproc(des,lf,nopc);
199 if (mi[MDEB]>0) printf("preprocess ok\n");
200 if (lf_error) return;
201 lf->ord = 0;
202 makecfn(des,lf);
203 if ((mi[MDIM]==1) && (lf->sty[0]!=STANGL))
204 { i = 1;
205 while ((i<mi[MN]) && (datum(lf,0,i)>=datum(lf,0,i-1))) i++;
206 lf->ord = (i==mi[MN]);
207 }
208
209 if (mi[MDEB]>0) printf("call eval structure\n");
210 switch(mi[MEV])
211 { case EPHULL: triang_start(des,lf); break;
212 case EDATA: dataf(des,lf); break;
213 case ECROS: crossf(des,lf); break;
214 case EGRID: gridf(des,lf); break;
215 case ETREE: atree_start(des,lf); break;
216 case EKDCE: mi[MKT] = KCE;
217 case EKDTR: kdtre_start(des,lf); break;
218 case EPRES: preset(des,lf); break;
219 case EXBAR: xbarf(des,lf); break;
220 case ENONE: lf->nv = lf->nce = 0;
221 return;
222 #ifdef CVERSION
223 case 100: do_scbsim(des,lf); break;
224 #endif
225 default: ERROR(("startlf: Invalid evaluation structure"));
226 }
227
228 /* renormalize for family=density */
229 if ((mi[MREN]) && (mi[MTG]==TDEN)) dens_renorm(lf,des);
230 }
231
232 #ifdef CVERSION
233 extern lfit lf;
234 extern design des;
235 extern plots pl[];
236 int curwin;
237 vari *vb;
238
nofit()239 INT nofit()
240 { if (lf.mi==NULL) return(1);
241 return(lf.mi[MEV]==ENULL);
242 }
243
endfit()244 void endfit()
245 { INT i;
246 for (i=0; i<MAXWIN; i++)
247 if (pl[i].track != NULL)
248 { curwin = i;
249 cmdint(pl[i].track);
250 }
251 }
252
drl(key,dv,mi)253 INT drl(key,dv,mi)
254 char *key;
255 INT *dv, *mi;
256 { INT i, nd;
257 nd = readilist(dv,key,0,mi[MDEG],0);
258 for (i=0; i<nd; i++)
259 { if ((dv[i]<1) | (dv[i]>mi[MDIM]))
260 ERROR(("drl: Invalid derivatives %s",key));
261 dv[i]--;
262 }
263 return(nd);
264 }
265
fitoptions(lf,vc,re)266 void fitoptions(lf,vc,re)
267 lfit *lf;
268 vari *vc;
269 INT re;
270 {
271 INT d = 0, n, i, i0, i1, *mi;
272 char kc, *key;
273 vari *v;
274
275 re &= (!nofit());
276 i0 = getarg(vc,"formula",1);
277 if ((!re) && (i0==0)) { ERROR(("no formula")); return; }
278 i1 = getarg(vc,"data",1);
279 if (i1>0) doreaddata(argval(vc,i1),(INT)0);
280 if (re)
281 recondat(0,&lf->mi[MN]);
282 else
283 { lf->base = lf->y = lf->c = lf->w = NULL;
284 lf->nd = 0;
285 strcpy(lf->yname,"_NuLl");
286 strcpy(lf->wname,"_NuLl");
287 strcpy(lf->bname,"_NuLl");
288 strcpy(lf->cname,"_NuLl");
289 }
290 if (i0>0) /* interpret formula */
291 { key = argval(vc,i0);
292 n = -1;
293 i0 = i1 = 0; d = 0;
294 while ((i0<strlen(key)) && (key[i0]!='~')) i0++;
295 if (key[i0] != '~') { ERROR(("invalid formula %s",key)); return; }
296 if (i0>0)
297 { key[i0] = '\0';
298 lf->y = vdptr(findvar(key,1,&n));
299 strcpy(lf->yname,key);
300 key[i0] = '~';
301 }
302 i1 = i0 = i0+1;
303 while (i1<strlen(key))
304 { while ((i1<strlen(key)) && (key[i1]!='+')) i1++;
305 kc = key[i1]; key[i1] = '\0';
306 lf->sty[d] = KPROD;
307 if (stm(&key[i0],"left(",5))
308 { lf->sty[d] = STLEFT;
309 i0 = i0+5; key[i1-1] = '\0';
310 }
311 else if (stm(&key[i0],"right(",6))
312 { lf->sty[d] = STRIGH;
313 i0 = i0+6; key[i1-1] = '\0';
314 }
315 else if (stm(&key[i0],"ang(",4))
316 { lf->sty[d] = STANGL;
317 i0 = i0+4; key[i1-1] = '\0';
318 }
319 else if (stm(&key[i0],"cpar(",5))
320 { lf->sty[d] = STCPAR;
321 i0 = i0+5; key[i1-1] = '\0';
322 }
323 dvari(lf,d) = vdptr(findvar(&key[i0],1,&n));
324 strcpy(lf->xname[d],&key[i0]);
325 if (lf->sty[d]!=KPROD) key[i1-1] = ')';
326 d++; key[i1] = kc;
327 i0 = i1 = i1+1;
328 }
329 fitdefault(lf,n,d);
330 }
331 mi = lf->mi;
332
333 i = getarg(vc,"weights",1);
334 if (i>0)
335 { lf->w = vdptr(findvar(argval(vc,i),1,&mi[MN]));
336 strcpy(lf->wname,argval(vc,i));
337 }
338 i = getarg(vc,"cens",1);
339 if (i>0)
340 { lf->c = vdptr(findvar(argval(vc,i),1,&mi[MN]));
341 strcpy(lf->cname,argval(vc,i));
342 }
343 i = getarg(vc,"base",1);
344 if (i>0)
345 { lf->base = vdptr(findvar(argval(vc,i),1,&mi[MN]));
346 strcpy(lf->bname,argval(vc,i));
347 }
348
349 i = getarg(vc,"scale",1);
350 if (i>0)
351 { if (argvalis(vc,i,"T"))
352 for (i=0; i<d; i++) lf->sca[i] = 0;
353 else if (argvalis(vc,i,"F"))
354 for (i=0; i<d; i++) lf->sca[i] = 1;
355 else
356 arvect(argval(vc,i),lf->sca,d,0);
357 }
358
359 i = getarg(vc,"vb",0);
360 if (i>0)
361 { lf->dp[DALP] = -1;
362 vb = arbuild(argval(vc,i),0,strlen(argval(vc,i))-1,NULL,0,1);
363 setvarname(vb,"_varband");
364 }
365 else
366 { i = getarg(vc,"alpha",1);
367 if (i>0) arvect(argval(vc,i),&lf->dp[DALP],3,1);
368 }
369
370 i = getarg(vc,"deg",1);
371 if (i>0)
372 { i = readilist(&mi[MDEG0],argval(vc,i),1,2,0);
373 if (i==1) mi[MDEG] = mi[MDEG0];
374 }
375
376 i = getarg(vc,"family",1);if (i>0) setstrval(mi,MTG,argval(vc,i));
377 i = getarg(vc,"link",1); if (i>0) setstrval(mi,MLINK,argval(vc,i));
378 i = getarg(vc,"ev",1);
379 if (i>0)
380 { v = findvar(argval(vc,i),0,NULL);
381 if (v!=NULL)
382 { mi[MEV] = EPRES;
383 lf->xxev= v;
384 lf->nvm = v->n;
385 }
386 else
387 setstrval(mi,MEV,argval(vc,i));
388 }
389 i = getarg(vc,"acri",1); if (i>0) setstrval(mi,MACRI,argval(vc,i));
390
391 i = getarg(vc,"mg",1);
392 if (i>0) readilist(lf->mg,argval(vc,i),1,MXDIM,1);
393
394 i = getarg(vc,"kt",1); if (i>0) setstrval(mi,MKT, argval(vc,i));
395 i = getarg(vc,"kern",1); if (i>0) setstrval(mi,MKER,argval(vc,i));
396 i = getarg(vc,"itype",1);if (i>0) setstrval(mi,MIT, argval(vc,i));
397
398 i = getarg(vc,"cut",1);
399 if (i>0) lf->dp[DCUT] = darith(argval(vc,i));
400
401 i = getarg(vc,"flim",1);
402 if (i>0) arvect(argval(vc,i),lf->fl,2*d,2);
403
404 i = getarg(vc,"xlim",1);
405 if (i>0) arvect(argval(vc,i),lf->xl,2*d,2);
406
407 i = getarg(vc,"deriv",0);
408 if (i>0) lf->nd = drl(argval(vc,i),lf->deriv,lf->mi);
409 i = getarg(vc,"dc",1); if (i>0) mi[MDC] = getlogic(vc,i);
410 i = getarg(vc,"maxk",1); if (i>0) readilist(&mi[MK],argval(vc,i),1,1,0);
411 i = getarg(vc,"mint",1); if (i>0) readilist(&mi[MMINT],argval(vc,i),1,1,0);
412 i = getarg(vc,"maxit",1); if (i>0) readilist(&mi[MMXIT],argval(vc,i),1,1,0);
413 i = getarg(vc,"renorm",1);if (i>0) mi[MREN] = getlogic(vc,i);
414 i = getarg(vc,"debug",1); if (i>0) readilist(&mi[MDEB],argval(vc,i),1,1,0);
415 }
416
clocfit(v,re)417 void clocfit(v,re)
418 INT re;
419 vari *v;
420 {
421 lf.ord = 0;
422 lf.kap[0] = lf.kap[1] = lf.kap[2] = 0.0; lf.nk = 0;
423 fitoptions(&lf,v,re);
424 if (lf_error)
425 { if (lf.mi!=NULL) lf.mi[MEV] = ENULL;
426 return;
427 }
428
429
430 lf.nv = 0;
431 if (lf.mi[MDEG0]==lf.mi[MDEG])
432 { startlf(&des,&lf,procv,0);
433 if (!lf_error) ressumm(&lf,&des);
434 }
435 else
436 startlf(&des,&lf,procvvord,0);
437 if (lf_error)
438 { if (!re) lf.mi[MEV] = ENULL;
439 return;
440 }
441
442 //printf("Evaluation structure %d, %d points.\n",lf.mi[MEV],lf.nv);
443 if (argarg(v,0) != NULL) dosavefit(&lf,argarg(v,0),"wb",(INT)0);
444 endfit();
445 }
446
447 #endif
448