rdiff: rDiff/src/locfit/Source/liblfev.c comparison

comparison rDiff/src/locfit/Source/liblfev.c @ 0:0f80a5141704

version 0.3 uploaded

author	vipints
date	Thu, 14 Feb 2013 23:38:36 -0500
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:0f80a5141704
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "mex.h"
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+extern void fitoptions();
+static double hmin, gmin, sig2, pen, vr, tb;
+static lfit *blf;
+static design *bdes;
+int procvbind(des,lf,v)
+design *des;
+lfit *lf;
+int v;
+{ double s0, s1, bi;
+int i, ii, k;
+k = procv_var(des,lf,v);
+wdiag(&lf->lfd, &lf->sp, des,des->wd,&lf->dv,0,1,0);
+s0 = s1 = 0.0;
+for (i=0; i<des->n; i++)
+{ ii = des->ind[i];
+s0+= prwt(&lf->lfd,ii)*des->wd[i]*des->wd[i];
+bi = prwt(&lf->lfd,ii)*fabs(des->wd[i]*ipower(dist(des,ii),deg(&lf->sp)+1));
+s1+= bi*bi;
+}
+vr += s0;
+tb += s1;
+return(k);
+}
+double bcri(h,c,cri)
+double h;
+int c, cri;
+{ double num, den, res[10];
+int (*pv)();
+if (c==DALP)
+blf->sp.nn = h;
+else
+blf->sp.fixh = h;
+if ((cri&63)==BIND)
+{ pv = procvbind;
+vr = tb = 0.0;
+}
+else pv = procvstd;
+if (cri<64) startlf(bdes,blf,pv,0);
+switch(cri&63)
+{ case BGCV:
+ressumm(blf,bdes,res);
+num = -2*blf->lfd.n*res[0];
+den = blf->lfd.n-res[1];
+return(num/(den*den));
+case BCP:
+ressumm(blf,bdes,res);
+return(-2*res[0]/sig2-blf->lfd.n+pen*res[1]);
+case BIND:
+return(vr+pen*pen*tb);
+}
+LERR(("bcri: unknown criterion"));
+return(0.0);
+}
+void bsel2(h0,g0,ifact,c,cri)
+double h0, g0, ifact;
+int c, cri;
+{ int done, inc;
+double h1, g1;
+h1 = h0; g1 = g0;
+done = inc = 0;
+while (!done)
+{ h1 *= 1+ifact;
+g0 = g1;
+g1 = bcri(h1,c,cri);
+if (g1<gmin) { hmin = h1; gmin = g1; }
+if (g1>g0) inc++; else inc = 0;
+switch(cri)
+{ case BIND:
+done = (inc>=4) & (vr<blf->fp.nv);
+break;
+default:
+done = (inc>=4);
+}
+}
+}
+void bsel3(h0,g0,ifact,c,cri)
+double h0, g0, ifact;
+int c, cri;
+{ double h1, g1;
+int i;
+hmin = h0; gmin = g0;
+for (i=-1; i<=1; i++) if (i!=0)
+{ h1 = h0*(1+i*ifact);
+g1 = bcri(h1,c,cri);
+if (g1<gmin) { hmin = h1; gmin = g1; }
+}
+return;
+}
+void bselect(lf,des,c,cri,pn)
+lfit *lf;
+design *des;
+int c, cri;
+double pn;
+{ double h0, g0, ifact;
+int i;
+pen = pn;
+blf = lf;
+bdes = des;
+if (cri==BIND) pen /= factorial(deg(&lf->sp)+1);
+hmin = h0 = (c==DFXH) ? fixh(&lf->sp) : nn(&lf->sp);
+if (h0==0) LERR(("bselect: initial bandwidth is 0"));
+if (lf_error) return;
+sig2 = 1.0;
+gmin = g0 = bcri(h0,c,cri);
+if (cri==BCP)
+{ sig2 = rv(&lf->fp);
+g0 = gmin = bcri(h0,c,cri+64);
+}
+ifact = 0.3;
+bsel2(h0,g0,ifact,c,cri);
+for (i=0; i<5; i++)
+{ ifact = ifact/2;
+bsel3(hmin,gmin,ifact,c,cri);
+}
+if (c==DFXH)
+fixh(&lf->sp) = hmin;
+else
+nn(&lf->sp) = hmin;
+startlf(des,lf,procvstd,0);
+ressumm(lf,des,lf->fp.kap);
+}
+double compsda(x,h,n)
+double *x, h;
+int n;
+/* n/(n-1) * int( fhat''(x)^2 dx ); bandwidth h */
+{ int i, j;
+double ik, sd, z;
+ik = wint(1,NULL,0,WGAUS);
+sd = 0;
+for (i=0; i<n; i++)
+for (j=i; j<n; j++)
+{ z = (x[i]-x[j])/h;
+sd += (2-(i==j))*Wconv4(z,WGAUS)/(ik*ik);
+}
+sd = sd/(n*(n-1)*h*h*h*h*h);
+return(sd);
+}
+double widthsj(x,lambda,n)
+double *x, lambda;
+int n;
+{ double ik, a, b, td, sa, z, c, c1, c2, c3;
+int i, j;
+a = GFACT*0.920*lambda*exp(-log((double)n)/7)/SQRT2;
+b = GFACT*0.912*lambda*exp(-log((double)n)/9)/SQRT2;
+ik = wint(1,NULL,0,WGAUS);
+td = 0;
+for (i=0; i<n; i++)
+for (j=i; j<n; j++)
+{ z = (x[i]-x[j])/b;
+td += (2-(i==j))*Wconv6(z,WGAUS)/(ik*ik);
+}
+td = -td/(n*(n-1));
+j = 2.0;
+c1 = Wconv4(0.0,WGAUS);
+c2 = wint(1,&j,1,WGAUS);
+c3 = Wconv(0.0,WGAUS);  /* (2*c1/(c2*c3))^(1/7)=1.357 */
+sa = compsda(x,a,n);
+c = b*exp(log(c1*ik/(c2*c3*GFACT*GFACT*GFACT*GFACT)*sa/td)/7)*SQRT2;
+return(c);
+}
+void kdecri(x,h,res,c,k,ker,n)
+double *x, h, *res, c;
+int k, ker, n;
+{ int i, j;
+double degfree, dfd, pen, s, r0, r1, d0, d1, ik, wij;
+if (h<=0) WARN(("kdecri, h = %6.4f",h));
+res[0] = res[1] = 0.0;
+ik = wint(1,NULL,0,ker);
+switch(k)
+{ case 1: /* aic */
+pen = 2.0;
+for (i=0; i<n; i++)
+{ r0 = d0 = 0.0;
+for (j=0; j<n; j++)
+{ s = (x[i]-x[j])/h;
+r0 += W(s,ker);
+d0 += s*s*Wd(s,ker);
+}
+d0 = -(d0+r0)/(n*h*h*ik);  /* d0 = d/dh fhat(xi) */
+r0 /= n*h*ik;              /* r0 = fhat(xi) */
+res[0] += -2*log(r0)+pen*W(0.0,ker)/(n*h*ik*r0);
+res[1] += -2*d0/r0-pen*W(0.0,ker)/(n*h*ik*r0)*(d0/r0+1.0/h);
+}
+return;
+case 2: /* ocv */
+for (i=0; i<n; i++)
+{ r0 = 0.0; d0 = 0.0;
+for (j=0; j<n; j++) if (i!=j)
+{ s = (x[i]-x[j])/h;
+r0 += W(s,ker);
+d0 += s*s*Wd(s,ker);
+}
+d0 = -(d0+r0)/((n-1)*h*h);
+r0 = r0/((n-1)*h);
+res[0] -= log(r0);
+res[1] -= d0/r0;
+}
+return;
+case 3: /* lscv */
+r0 = r1 = d0 = d1 = degfree = 0.0;
+for (i=0; i<n; i++)
+{ dfd = 0.0;
+for (j=0; j<n; j++)
+{ s = (x[i]-x[j])/h;
+wij = W(s,ker);
+dfd += wij;
+/*
+*  r0 = \int fhat * fhat = sum_{i,j} W*W( (Xi-Xj)/h ) / n^2 h.
+*  d0 is it's derivative wrt h.
+*
+*  r1 = 1/n sum( f_{-i}(X_i) ).
+*  d1 is  it's derivative wrt h.
+*
+*  degfree = sum_i ( W_0 / sum_j W( (Xi-Xj)/h ) ) is fitted d.f.
+*/
+r0 += Wconv(s,ker);
+d0 += -s*s*Wconv1(s,ker);
+if (i != j)
+{ r1 += wij;
+d1 += -s*s*Wd(s,ker);
+}
+}
+degfree += 1.0/dfd;
+}
+d1 -= r1;
+d0 -= r0;
+res[0] = r0/(n*n*h*ik*ik)   - 2*r1/(n*(n-1)*h*ik);
+res[1] = d0/(n*n*h*h*ik*ik) - 2*d1/(n*(n-1)*h*h*ik);
+res[2] = degfree;
+return;
+case 4: /* bcv */
+r0 = d0 = 0.0;
+for (i=0; i<n; i++)
+for (j=i+1; j<n; j++)
+{ s = (x[i]-x[j])/h;
+r0 += 2*Wconv4(s,ker);
+d0 += 2*s*Wconv5(s,ker);
+}
+d0 = (-d0-r0)/(n*n*h*h*ik*ik);
+r0 = r0/(n*n*h*ik*ik);
+j = 2.0;
+d1 = wint(1,&j,1,ker);
+r1 = Wconv(0.0,ker);
+res[0] = (d1*d1*r0/4+r1/(n*h))/(ik*ik);
+res[1] = (d1*d1*d0/4-r1/(n*h*h))/(ik*ik);
+return;
+case 5: /* sjpi */
+s = c*exp(5*log(h)/7)/SQRT2;
+d0 = compsda(x,s,n);
+res[0] = d0; /* this is S(alpha) in SJ */
+res[1] = exp(log(Wikk(WGAUS,0)/(d0*n))/5)-h;
+return;
+case 6: /* gas-k-k */
+s = exp(log(1.0*n)/10)*h;
+d0 = compsda(x,s,n);
+res[0] = d0;
+res[1] = exp(log(Wikk(WGAUS,0)/(d0*n))/5)-h;
+return;
+}
+LERR(("kdecri: what???"));
+return;
+}
+double esolve(x,j,h0,h1,k,c,ker,n)
+double *x, h0, h1, c;
+int j, k, ker, n;
+{ double h[7], d[7], r[7], res[4], min, minh, fact;
+int i, nc;
+min = 1.0e30; minh = 0.0;
+fact = 1.00001;
+h[6] = h0; kdecri(x,h[6],res,c,j,ker,n);
+r[6] = res[0]; d[6] = res[1];
+if (lf_error) return(0.0);
+nc = 0;
+for (i=0; i<k; i++)
+{ h[5] = h[6]; r[5] = r[6]; d[5] = d[6];
+h[6] = h0*exp((i+1)*log(h1/h0)/k);
+kdecri(x,h[6],res,c,j,ker,n);
+r[6] = res[0]; d[6] = res[1];
+if (lf_error) return(0.0);
+if (d[5]*d[6]<0)
+{ h[2] = h[0] = h[5]; d[2] = d[0] = d[5]; r[2] = r[0] = r[5];
+h[3] = h[1] = h[6]; d[3] = d[1] = d[6]; r[3] = r[1] = r[6];
+while ((h[3]>fact*h[2])|(h[2]>fact*h[3]))
+{ h[4] = h[3]-d[3]*(h[3]-h[2])/(d[3]-d[2]);
+if ((h[4]<h[0]) | (h[4]>h[1])) h[4] = (h[0]+h[1])/2;
+kdecri(x,h[4],res,c,j,ker,n);
+r[4] = res[0]; d[4] = res[1];
+if (lf_error) return(0.0);
+h[2] = h[3]; h[3] = h[4];
+d[2] = d[3]; d[3] = d[4];
+r[2] = r[3]; r[3] = r[4];
+if (d[4]*d[0]>0) { h[0] = h[4]; d[0] = d[4]; r[0] = r[4]; }
+else { h[1] = h[4]; d[1] = d[4]; r[1] = r[4]; }
+}
+if (j>=4) return(h[4]); /* first min for BCV etc */
+if (r[4]<=min) { min = r[4]; minh = h[4]; }
+nc++;
+}
+}
+if (nc==0) minh = (r[5]<r[6]) ? h0 : h1;
+return(minh);
+}
+void kdeselect(band,x,ind,h0,h1,meth,nm,ker,n)
+double h0, h1, *band, *x;
+int *ind, nm, ker, n, *meth;
+{ double scale, c;
+int i, k;
+k = n/4;
+for (i=0; i<n; i++) ind[i] = i;
+scale = kordstat(x,n+1-k,n,ind) - kordstat(x,k,n,ind);
+c = widthsj(x,scale,n);
+for (i=0; i<nm; i++)
+band[i] = esolve(x,meth[i],h0,h1,10,c,ker,n);
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+/*
+*   The function dens_integrate(lf,des,z) is used to integrate a density
+*   estimate (z=1) or the density squared (z=2). This is used to renormalize
+*   the estimate (function dens_renorm) or in the computation of LSCV
+*   (in modlscv.c). The implementation is presently for d=1.
+*
+*   The computation orders the fit points selected by locfit, and
+*   integrates analytically over each interval. For the log-link,
+*   the interpolant used is peicewise quadratic (with one knot in
+*   the middle of each interval); this differs from the cubic interpolant
+*   used elsewhere in Locfit.
+*
+*   TODO: allow for xlim. What can be done simply in >=2 dimensions?
+*/
+#include "lfev.h"
+/*
+* Finds the order of observations in the array x, and
+* stores in integer array ind.
+* At input, lset l=0 and r=length(x)-1.
+* At output, x[ind[0]] <= x[ind[1]] <= ...
+*/
+void lforder(ind,x,l,r)
+int *ind, l, r;
+double *x;
+{ double piv;
+int i, i0, i1;
+piv = (x[ind[l]]+x[ind[r]])/2;
+i0 = l; i1 = r;
+while (i0<=i1)
+{ while ((i0<=i1) && (x[ind[i0]]<=piv)) i0++;
+while ((i0<=i1) && (x[ind[i1]]>piv))  i1--;
+if (i0<i1)
+{ ISWAP(ind[i0],ind[i1]);
+i0++; i1--;
+}
+}
+/* now, x[ind[l..i1]] <= piv < x[ind[i0..r]].
+put the ties in the middle */
+while ((i1>=l) && (x[ind[i1]]==piv)) i1--;
+for (i=l; i<=i1; i++)
+if (x[ind[i]]==piv)
+{ ISWAP(ind[i],ind[i1]);
+while (x[ind[i1]]==piv) i1--;
+}
+if (l<i1) lforder(ind,x,l,i1);
+if (i0<r) lforder(ind,x,i0,r);
+}
+/*
+*  estdiv integrates the density between fit points x0 and x1.
+*  f0, f1 are function values, d0, d1 are derivatives.
+*/
+double estdiv(x0,x1,f0,f1,d0,d1,lin)
+double x0, x1, f0, f1, d0, d1;
+int lin;
+{ double cf[4], I[2], dlt, e0, e1;
+if (x0==x1) return(0.0);
+if (lin==LIDENT)
+{
+/* cf are integrals of hermite polynomials.
+* Then adjust for x1-x0.
+*/
+cf[0] = cf[1] = 0.5;
+cf[2] = 1.0/12.0; cf[3] = -cf[2];
+return( (cf[0]*f0+cf[1]*f1)*(x1-x0)
++ (cf[2]*d0+cf[3]*d1)*(x1-x0)*(x1-x0) );
+}
+/*
+* this is for LLOG
+*/
+dlt = (x1-x0)/2;
+cf[0] = f0;
+cf[1] = d0;
+cf[2] = ( 2*(f1-f0) - dlt*(d1+3*d0) )/(4*dlt*dlt);
+recurint(0.0,dlt,cf,I,0,WRECT);
+e0 = I[0];
+cf[0] = f1;
+cf[1] = -d1;
+cf[2] = ( 2*(f0-f1) + dlt*(d0+3*d1) )/( 4*dlt*dlt );
+recurint(0.0,dlt,cf,I,0,WRECT);
+e1 = I[0];
+return(e0+e1);
+}
+/*
+*   Evaluate the integral of the density estimate to the power z.
+*   This would be severely messed up, if I ever implement parcomp
+*   for densities.
+*/
+double dens_integrate(lf,des,z)
+lfit *lf;
+design *des;
+int z;
+{ int has_deriv, i, i0, i1, nv, *ind;
+double *xev, *fit, *deriv, sum, term;
+double d0, d1, f0, f1;
+fitpt *fp;
+fp = &lf->fp;
+if (fp->d >= 2)
+{ WARN(("dens_integrate requires d=1"));
+return(0.0);
+}
+has_deriv = (deg(&lf->sp) > 0); /* not right? */
+fit = fp->coef;
+if (has_deriv)
+deriv = &fit[fp->nvm];
+xev = evp(fp);
+/*
+* order the vertices
+*/
+nv = fp->nv;
+if (lf->lfd.n<nv) return(0.0);
+ind = des->ind;
+for (i=0; i<nv; i++) ind[i] = i;
+lforder(ind,xev,0,nv-1);
+sum = 0.0;
+/*
+* Estimate the contribution of the boundaries.
+* should really check flim here.
+*/
+i0 = ind[0]; i1 = ind[1];
+f1 = fit[i0];
+d1 = (has_deriv) ? deriv[i0] :
+(fit[i1]-fit[i0])/(xev[i1]-xev[i0]);
+if (d1 <= 0) WARN(("dens_integrate - ouch!"));
+if (z==2)
+{ if (link(&lf->sp)==LLOG)
+{ f1 *= 2; d1 *= 2; }
+else
+{ d1 = 2*d1*f1; f1 = f1*f1; }
+}
+term = (link(&lf->sp)==LIDENT) ? f1*f1/(2*d1) : exp(f1)/d1;
+sum += term;
+i0 = ind[nv-2]; i1 = ind[nv-1];
+f0 = fit[i1];
+d0 = (has_deriv) ? deriv[i1] :
+(fit[i1]-fit[i0])/(xev[i1]-xev[i0]);
+if (d0 >= 0) WARN(("dens_integrate - ouch!"));
+if (z==2)
+{ if (link(&lf->sp)==LLOG)
+{ f0 *= 2; d0 *= 2; }
+else
+{ d0 = 2*d0*f0; f0 = f0*f0; }
+}
+term = (link(&lf->sp)==LIDENT) ? -f0*f0/(2*d0) : exp(f0)/d0;
+sum += term;
+for (i=1; i<nv; i++)
+{ i0 = ind[i-1]; i1 = ind[i];
+f0 = fit[i0]; f1 = fit[i1];
+d0 = (has_deriv) ? deriv[i0] :
+(f1-f0)/(xev[i1]-xev[i0]);
+d1 = (has_deriv) ? deriv[i1] : d0;
+if (z==2)
+{ if (link(&lf->sp)==LLOG)
+{ f0 *= 2; f1 *= 2; d0 *= 2; d1 *= 2; }
+else
+{ d0 *= 2*f0; d1 *= 2*f1; f0 = f0*f0; f1 = f1*f1; }
+}
+term = estdiv(xev[i0],xev[i1],f0,f1,d0,d1,link(&lf->sp));
+sum += term;
+}
+return(sum);
+}
+void dens_renorm(lf,des)
+lfit *lf;
+design *des;
+{ int i;
+double sum;
+sum = dens_integrate(lf,des,1);
+if (sum==0.0) return;
+sum = log(sum);
+for (i=0; i<lf->fp.nv; i++) lf->fp.coef[i] -= sum;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+/*
+*   This file contains functions for constructing and
+*   interpolating the adaptive tree structure. This is
+*   the default evaluation structure used by Locfit.
+*/
+#include "lfev.h"
+/*
+Guess the number of fitting points.
+Needs improving!
+*/
+void atree_guessnv(evs,nvm,ncm,vc,d,alp)
+evstruc *evs;
+double alp;
+int *nvm, *ncm, *vc, d;
+{ double a0, cu, ifl;
+int i, nv, nc;
+*ncm = 1<<30; *nvm = 1<<30;
+*vc = 1 << d;
+if (alp>0)
+{ a0 = (alp > 1) ? 1 : 1/alp;
+if (cut(evs)<0.01)
+{ WARN(("guessnv: cut too small."));
+cut(evs) = 0.01;
+}
+cu = 1;
+for (i=0; i<d; i++) cu *= MIN(1.0,cut(evs));
+nv = (int)((5*a0/cu)**vc);   /* this allows 10*a0/cu splits */
+nc = (int)(10*a0/cu+1);      /* and 10*a0/cu cells */
+if (nv<*nvm) *nvm = nv;
+if (nc<*ncm) *ncm = nc;
+}
+if (*nvm == 1<<30) /* by default, allow 100 splits */
+{ *nvm = 102**vc;
+*ncm = 201;
+}
+/* inflation based on mk */
+ifl = mk(evs)/100.0;
+*nvm = *vc+(int)(ifl**nvm);
+*ncm = (int)(ifl**ncm);
+}
+/*
+Determine whether a cell in the tree needs splitting.
+If so, return the split variable (0..d-1).
+Otherwise, return -1.
+*/
+int atree_split(lf,ce,le,ll,ur)
+lfit *lf;
+int *ce;
+double *le, *ll, *ur;
+{ int d, vc, i, is;
+double h, hmin, score[MXDIM];
+d = lf->fp.d; vc = 1<<d;
+hmin = 0.0;
+for (i=0; i<vc; i++)
+{ h = lf->fp.h[ce[i]];
+if ((h>0) && ((hmin==0)|(h<hmin))) hmin = h;
+}
+is = 0;
+for (i=0; i<d; i++)
+{ le[i] = (ur[i]-ll[i])/lf->lfd.sca[i];
+if ((lf->lfd.sty[i]==STCPAR) || (hmin==0))
+score[i] = 2*(ur[i]-ll[i])/(lf->evs.fl[i+d]-lf->evs.fl[i]);
+else
+score[i] = le[i]/hmin;
+if (score[i]>score[is]) is = i;
+}
+if (cut(&lf->evs)<score[is]) return(is);
+return(-1);
+}
+/*
+recursively grow the tree structure, begining with the parent cell.
+*/
+void atree_grow(des,lf,ce,ct,term,ll,ur)
+design *des;
+lfit *lf;
+int *ce, *ct, *term;
+double *ll, *ur;
+{ int nce[1<<MXDIM];
+int i, i0, i1, d, vc, pv, tk, ns;
+double le[MXDIM], z;
+d = lf->fp.d; vc = 1<<d;
+/* does this cell need splitting?
+If not, wrap up and return.
+*/
+ns = atree_split(lf,ce,le,ll,ur);
+if (ns==-1)
+{ if (ct != NULL) /* reconstructing terminal cells */
+{ for (i=0; i<vc; i++) term[*ct*vc+i] = ce[i];
+(*ct)++;
+}
+return;
+}
+/* split the cell at the midpoint on side ns */
+tk = 1<<ns;
+for (i=0; i<vc; i++)
+{ if ((i&tk)==0) nce[i] = ce[i];
+else
+{ i0 = ce[i];
+i1 = ce[i-tk];
+pv = (lf->lfd.sty[i]!=STCPAR) &&
+(le[ns] < (cut(&lf->evs)*MIN(lf->fp.h[i0],lf->fp.h[i1])));
+nce[i] = newsplit(des,lf,i0,i1,pv);
+if (lf_error) return;
+}
+}
+z = ur[ns]; ur[ns] = (z+ll[ns])/2;
+atree_grow(des,lf,nce,ct,term,ll,ur);
+if (lf_error) return;
+ur[ns] = z;
+for (i=0; i<vc; i++)
+nce[i] = ((i&tk)== 0) ? nce[i+tk] : ce[i];
+z = ll[ns]; ll[ns] = (z+ur[ns])/2;
+atree_grow(des,lf,nce,ct,term,ll,ur);
+if (lf_error) return;
+ll[ns] = z;
+}
+void atree_start(des,lf)
+design *des;
+lfit *lf;
+{ int d, i, j, k, vc, ncm, nvm;
+double ll[MXDIM], ur[MXDIM];
+if (lf_debug>1) mut_printf(" In atree_start\n");
+d = lf->fp.d;
+atree_guessnv(&lf->evs,&nvm,&ncm,&vc,d,nn(&lf->sp));
+if (lf_debug>2) mut_printf(" atree_start: nvm %d ncm %d\n",nvm,ncm);
+trchck(lf,nvm,ncm,vc);
+/* Set the lower left, upper right limits. */
+for (j=0; j<d; j++)
+{ ll[j] = lf->evs.fl[j];
+ur[j] = lf->evs.fl[j+d];
+}
+/* Set the initial cell; fit at the vertices. */
+for (i=0; i<vc; i++)
+{ j = i;
+for (k=0; k<d; ++k)
+{ evptx(&lf->fp,i,k) = (j%2) ? ur[k] : ll[k];
+j >>= 1;
+}
+lf->evs.ce[i] = i;
+PROC_VERTEX(des,lf,i);
+if (lf_error) return;
+lf->evs.s[i] = 0;
+}
+lf->fp.nv = vc;
+/* build the tree */
+atree_grow(des,lf,lf->evs.ce,NULL,NULL,ll,ur);
+lf->evs.nce = 1;
+}
+double atree_int(lf,x,what)
+lfit *lf;
+double *x;
+int what;
+{ double vv[64][64], *ll, *ur, h, xx[MXDIM];
+int lo, tk, ns, nv, nc, d, i, vc;
+int ce[64];
+fitpt *fp;
+evstruc *evs;
+fp = &lf->fp;
+evs= &lf->evs;
+d = fp->d;
+vc = 1<<d;
+for (i=0; i<vc; i++)
+{ setzero(vv[i],vc);
+nc = exvval(fp,vv[i],i,d,what,1);
+ce[i] = evs->ce[i];
+}
+ns = 0;
+while(ns!=-1)
+{ ll = evpt(fp,ce[0]); ur = evpt(fp,ce[vc-1]);
+ns = atree_split(lf,ce,xx,ll,ur);
+if (ns!=-1)
+{ tk = 1<<ns;
+h = ur[ns]-ll[ns];
+lo = (2*(x[ns]-ll[ns])) < h;
+for (i=0; i<vc; i++) if ((tk&i)==0)
+{ nv = findpt(fp,evs,(int)ce[i],(int)ce[i+tk]);
+if (nv==-1) LERR(("Descend tree problem"));
+if (lf_error) return(0.0);
+if (lo)
+{ ce[i+tk] = nv;
+if (evs->s[nv]) exvvalpv(vv[i+tk],vv[i],vv[i+tk],d,ns,h,nc);
+else exvval(fp,vv[i+tk],nv,d,what,1);
+}
+else
+{ ce[i] = nv;
+if (evs->s[nv]) exvvalpv(vv[i],vv[i],vv[i+tk],d,ns,h,nc);
+else exvval(fp,vv[i],nv,d,what,1);
+} }
+} }
+ll = evpt(fp,ce[0]); ur = evpt(fp,ce[vc-1]);
+return(rectcell_interp(x,vv,ll,ur,d,nc));
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+double linear_interp(h,d,f0,f1)
+double h, d, f0, f1;
+{ if (d==0) return(f0);
+return( ( (d-h)*f0 + h*f1 ) / d );
+}
+void hermite2(x,z,phi)
+double x, z, *phi;
+{ double h;
+if (z==0)
+{ phi[0] = 1.0; phi[1] = phi[2] = phi[3] = 0.0;
+return;
+}
+h = x/z;
+if (h<0)
+{ phi[0] = 1; phi[1] = 0;
+phi[2] = h; phi[3] = 0;
+return;
+}
+if (h>1)
+{ phi[0] = 0; phi[1] = 1;
+phi[2] = 0; phi[3] = h-1;
+return;
+}
+phi[1] = h*h*(3-2*h);
+phi[0] = 1-phi[1];
+phi[2] = h*(1-h)*(1-h);
+phi[3] = h*h*(h - 1);
+}
+double cubic_interp(h,f0,f1,d0,d1)
+double h, f0, f1, d0, d1;
+{ double phi[4];
+hermite2(h,1.0,phi);
+return(phi[0]*f0+phi[1]*f1+phi[2]*d0+phi[3]*d1);
+}
+double cubintd(h,f0,f1,d0,d1)
+double h, f0, f1, d0, d1;
+{ double phi[4];
+phi[1] = 6*h*(1-h);
+phi[0] = -phi[1];
+phi[2] = (1-h)*(1-3*h);
+phi[3] = h*(3*h-2);
+return(phi[0]*f0+phi[1]*f1+phi[2]*d0+phi[3]*d1);
+}
+/*
+interpolate over a rectangular cell.
+x = interpolation point.
+vv = array of vertex values.
+ll = lower left corner.
+ur = upper right corner.
+d = dimension.
+nc = no of coefficients.
+*/
+double rectcell_interp(x,vv,ll,ur,d,nc)
+double *x, vv[64][64], *ll, *ur;
+int d, nc;
+{ double phi[4];
+int i, j, k, tk;
+tk = 1<<d;
+for (i=0; i<tk; i++) if (vv[i][0]==NOSLN) return(NOSLN);
+/* no derivatives - use multilinear interpolation */
+if (nc==1)
+{ for (i=d-1; i>=0; i--)
+{ tk = 1<<i;
+for (j=0; j<tk; j++)
+vv[j][0] = linear_interp(x[i]-ll[i],ur[i]-ll[i],vv[j][0],vv[j+tk][0]);
+}
+return(vv[0][0]);
+}
+/* with derivatives -- use cubic */
+if (nc==d+1)
+{ for (i=d-1; i>=0; i--)
+{ hermite2(x[i]-ll[i],ur[i]-ll[i],phi);
+tk = 1<<i;
+phi[2] *= ur[i]-ll[i];
+phi[3] *= ur[i]-ll[i];
+for (j=0; j<tk; j++)
+{ vv[j][0] = phi[0]*vv[j][0] + phi[1]*vv[j+tk][0]
++ phi[2]*vv[j][i+1] + phi[3]*vv[j+tk][i+1];
+for (k=1; k<=i; k++)
+vv[j][k] = phi[0]*vv[j][k] + phi[1]*vv[j+tk][k];
+}
+}
+return(vv[0][0]);
+}
+/* with all coefs -- use multicubic */
+for (i=d-1; i>=0; i--)
+{ hermite2(x[i]-ll[i],ur[i]-ll[i],phi);
+tk = 1<<i;
+phi[2] *= ur[i]-ll[i];
+phi[3] *= ur[i]-ll[i];
+for (j=0; j<tk; j++)
+for (k=0; k<tk; k++)
+vv[j][k] = phi[0]*vv[j][k] + phi[1]*vv[j+tk][k]
++ phi[2]*vv[j][k+tk] + phi[3]*vv[j+tk][k+tk];
+}
+return(vv[0][0]);
+}
+int exvval(fp,vv,nv,d,what,z)
+fitpt *fp;
+double *vv;
+int nv, d, z, what;
+{ int i, k;
+double *values;
+k = (z) ? 1<<d : d+1;
+for (i=1; i<k; i++) vv[i] = 0.0;
+switch(what)
+{ case PCOEF:
+values = fp->coef;
+break;
+case PVARI:
+case PNLX:
+values = fp->nlx;
+break;
+case PT0:
+values = fp->t0;
+break;
+case PBAND:
+vv[0] = fp->h[nv];
+return(1);
+case PDEGR:
+vv[0] = fp->deg[nv];
+return(1);
+case PLIK:
+vv[0] = fp->lik[nv];
+return(1);
+case PRDF:
+vv[0] = fp->lik[2*fp->nvm+nv];
+return(1);
+default:
+LERR(("Invalid what in exvval"));
+return(0);
+}
+vv[0] = values[nv];
+if (!fp->hasd) return(1);
+if (z)
+{ for (i=0; i<d; i++) vv[1<<i] = values[(i+1)*fp->nvm+nv];
+return(1<<d);
+}
+else
+{ for (i=1; i<=d; i++) vv[i] = values[i*fp->nvm+nv];
+return(d+1);
+}
+}
+void exvvalpv(vv,vl,vr,d,k,dl,nc)
+double *vv, *vl, *vr, dl;
+int d, k, nc;
+{ int i, tk, td;
+double f0, f1;
+if (nc==1)
+{ vv[0] = (vl[0]+vr[0])/2;
+return;
+}
+tk = 1<<k;
+td = 1<<d;
+for (i=0; i<td; i++) if ((i&tk)==0)
+{ f0 = (vl[i]+vr[i])/2 + dl*(vl[i+tk]-vr[i+tk])/8;
+f1 = 1.5*(vr[i]-vl[i])/dl - (vl[i+tk]+vr[i+tk])/4;
+vv[i] = f0;
+vv[i+tk] = f1;
+}
+}
+double grid_int(fp,evs,x,what)
+fitpt *fp;
+evstruc *evs;
+double *x;
+int what;
+{ int d, i, j, jj, nc, sk, v[MXDIM], vc, z0, nce[1<<MXDIM], *mg;
+double *ll, *ur, vv[64][64], z;
+d = fp->d;
+ll = evpt(fp,0); ur = evpt(fp,fp->nv-1);
+mg = mg(evs);
+z0 = 0; vc = 1<<d;
+for (j=d-1; j>=0; j--)
+{ v[j] = (int)((mg[j]-1)*(x[j]-ll[j])/(ur[j]-ll[j]));
+if (v[j]<0) v[j]=0;
+if (v[j]>=mg[j]-1) v[j] = mg[j]-2;
+z0 = z0*mg[j]+v[j];
+}
+nce[0] = z0; nce[1] = z0+1; sk = jj = 1;
+for (i=1; i<d; i++)
+{ sk *= mg[i-1];
+jj<<=1;
+for (j=0; j<jj; j++)
+nce[j+jj] = nce[j]+sk;
+}
+for (i=0; i<vc; i++)
+nc = exvval(fp,vv[i],nce[i],d,what,1);
+ll = evpt(fp,nce[0]);
+ur = evpt(fp,nce[vc-1]);
+z = rectcell_interp(x,vv,ll,ur,d,nc);
+return(z);
+}
+double fitp_int(fp,x,what,i)
+fitpt *fp;
+double *x;
+int what, i;
+{ double vv[1+MXDIM];
+exvval(fp,vv,i,fp->d,what,0);
+return(vv[0]);
+}
+double xbar_int(fp,x,what)
+fitpt *fp;
+double *x;
+int what;
+{ int i, nc;
+double vv[1+MXDIM], f;
+nc = exvval(fp,vv,0,fp->d,what,0);
+f = vv[0];
+if (nc>1)
+for (i=0; i<fp->d; i++)
+f += vv[i+1]*(x[i]-evptx(fp,0,i));
+return(f);
+}
+double dointpoint(lf,x,what,ev,j)
+lfit *lf;
+double *x;
+int what, ev, j;
+{ double xf, f, link[LLEN];
+int i, rt;
+fitpt *fp;
+evstruc *evs;
+fp = &lf->fp; evs = &lf->evs;
+for (i=0; i<fp->d; i++) if (lf->lfd.sty[i]==STANGL)
+{ xf = floor(x[i]/(2*PI*lf->lfd.sca[i]));
+x[i] -= xf*2*PI*lf->lfd.sca[i];
+}
+if (what > 64)
+{ rt = what-64;
+what = PCOEF;
+}
+else rt = 0;
+switch(ev)
+{ case EGRID: f = grid_int(fp,evs,x,what); break;
+case EKDTR: f = kdtre_int(fp,evs,x,what); break;
+case ETREE: f = atree_int(lf,x,what); break;
+case EPHULL: f = triang_int(lf,x,what); break;
+case EFITP: f = fitp_int(fp,x,what,j); break;
+case EXBAR: f = xbar_int(fp,x,what); break;
+case ENONE: f = 0; break;
+case ESPHR: f = sphere_int(lf,x,what); break;
+default: LERR(("dointpoint: cannot interpolate structure %d",ev));
+}
+if (((what==PT0)|(what==PNLX)) && (f<0)) f = 0.0;
+f += addparcomp(lf,x,what);
+if (rt>0)
+{
+stdlinks(link,&lf->lfd,&lf->sp,j,f,rsc(fp));
+f = resid(resp(&lf->lfd,j),prwt(&lf->lfd,j),f,fam(&lf->sp),rt,link);
+}
+return(f);
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+/*
+*   Routines for building and interpolating the kd tree.
+*   Initially, this started from the loess code.
+*
+*   Todo: EKDCE isn't working.
+*/
+#include "lfev.h"
+void newcell();
+static int nterm;
+void kdtre_guessnv(evs,nvm,ncm,vc,n,d,alp)
+evstruc *evs;
+double alp;
+int *nvm, *ncm, *vc, n, d;
+{ int k;
+if (ev(evs) == EKDTR)
+{ nterm = (int)(cut(evs)/4 * n * MIN(alp,1.0) );
+k = 2*n/nterm;
+*vc = 1<<d;
+*ncm = 2*k+1;
+*nvm = (k+2)**vc/2;
+return;
+}
+if (ev(evs) == EKDCE)
+{ nterm = (int)(n * alp);
+*vc = 1;
+*nvm = 1+(int)(2*n/nterm);
+*ncm = 2**nvm+1;
+return;
+}
+*nvm = *ncm = *vc = 0;
+return;
+}
+/*
+Split x[pi[l..r]] into two equal sized sets.
+Let m=(l+r)/2.
+At return,
+x[pi[l..m]] < x[pi[m+1..r]].
+Assuming no ties:
+If l+r is odd, the sets have the same size.
+If l+r is even, the low set is larger by 1.
+If there are ties, all ties go in the low set.
+*/
+int ksmall(l, r, m, x, pi)
+int l, r, m, *pi;
+double *x;
+{
+int il, ir, jl, jr;
+double t;
+while(l<r)
+{ t = x[pi[m]];
+/*
+permute the observations so that
+x[pi[l..il]] < t <= x[pi[ir..r]].
+*/
+ir = l; il = r;
+while (ir<il)
+{ while ((ir<=r) && (x[pi[ir]] < t)) ir++;
+while ((il>=l) && (x[pi[il]]>= t)) il--;
+if (ir<il) ISWAP(pi[ir],pi[il]);
+}
+/*
+move  = t to the middle:
+x[pi[l..il]]  < t
+x[pi[jl..jr]] = t
+x[pi[ir..r]] > t
+*/
+jl = ir; jr = r;
+while (ir<jr)
+{ while ((ir<=r)  && (x[pi[ir]]== t)) ir++;
+while ((jr>=jl) && (x[pi[jr]] > t)) jr--;
+if (ir<jr) ISWAP(pi[ir],pi[jr]);
+}
+/*
+we're done if m is in the middle, jl <= m <= jr.
+*/
+if ((jl<=m) & (jr>=m)) return(jr);
+/*
+update l or r.
+*/
+if (m>=ir) l = ir;
+if (m<=il) r = il;
+}
+if (l==r) return(l);
+LERR(("ksmall failure"));
+return(0);
+}
+int terminal(lf,p,pi,fc,d,m,split_val)
+lfit *lf;
+int p, d, fc, *m, *pi;
+double *split_val;
+{ int i, k, lo, hi, split_var;
+double max, min, score, max_score, t;
+/*
+if there are fewer than fc points in the cell, this cell
+is terminal.
+*/
+lo = lf->evs.lo[p]; hi = lf->evs.hi[p];
+if (hi-lo < fc) return(-1);
+/* determine the split variable */
+max_score = 0.0; split_var = 0;
+for (k=0; k<d; k++)
+{ max = min = datum(&lf->lfd, k, pi[lo]);
+for (i=lo+1; i<=hi; i++)
+{ t = datum(&lf->lfd,k,pi[i]);
+if (t<min) min = t;
+if (t>max) max = t;
+}
+score = (max-min) / lf->lfd.sca[k];
+if (score > max_score)
+{ max_score = score;
+split_var = k;
+}
+}
+if (max_score==0) /* all points in the cell are equal */
+return(-1);
+*m = ksmall(lo,hi,(lo+hi)/2, dvari(&lf->lfd,split_var), pi);
+*split_val = datum(&lf->lfd, split_var, pi[*m]);
+if (*m==hi) /* all observations go lo */
+return(-1);
+return(split_var);
+}
+void kdtre_start(des,lf)
+design *des;
+lfit *lf;
+{
+int i, j, vc, d, nc, nv, ncm, nvm, k, m, n, p, *pi;
+double sv;
+d = lf->lfd.d; n = lf->lfd.n; pi = des->ind;
+kdtre_guessnv(&lf->evs,&nvm,&ncm,&vc,n,d,nn(&lf->sp));
+trchck(lf,nvm,ncm,vc);
+nv = 0;
+if (ev(&lf->evs) != EKDCE)
+{ for (i=0; i<vc; i++)
+{ j = i;
+for (k=0; k<d; ++k)
+{ evptx(&lf->fp,i,k) = lf->evs.fl[d*(j%2)+k];
+j >>= 1;
+}
+}
+nv = vc;
+for (j=0; j<vc; j++) lf->evs.ce[j] = j;
+}
+for (i=0; i<n; i++) pi[i] = i;
+p = 0; nc = 1;
+lf->evs.lo[p] = 0; lf->evs.hi[p] = n-1;
+lf->evs.s[p] = -1;
+while (p<nc)
+{ k = terminal(lf,p,pi,nterm,d,&m,&sv);
+if (k>=0)
+{
+if ((ncm<nc+2) | (2*nvm<2*nv+vc))
+{ WARN(("Insufficient space for full tree"));
+lf->evs.nce = nc; lf->fp.nv = nv;
+return;
+}
+/* new lo cell has obsn's lo[p]..m */
+lf->evs.lo[nc] = lf->evs.lo[p];
+lf->evs.hi[nc] = m;
+lf->evs.s[nc] = -1;
+/* new hi cell has obsn's m+1..hi[p] */
+lf->evs.lo[nc+1] = m+1;
+lf->evs.hi[nc+1] = lf->evs.hi[p];
+lf->evs.s[nc+1] = -1;
+/* cell p is split on variable k, value sv */
+lf->evs.s[p] = k;
+lf->evs.sv[p] = sv;
+lf->evs.lo[p] = nc; lf->evs.hi[p] = nc+1;
+nc=nc+2; i = nv;
+/* now compute the new vertices. */
+if (ev(&lf->evs) != EKDCE)
+newcell(&nv,vc,evp(&lf->fp), d, k, sv,
+&lf->evs.ce[p*vc], &lf->evs.ce[(nc-2)*vc], &lf->evs.ce[(nc-1)*vc]);
+}
+else if (ev(&lf->evs)==EKDCE) /* new vertex at cell center */
+{ sv = 0;
+for (i=0; i<d; i++) evptx(&lf->fp,nv,i) = 0;
+for (j=lf->evs.lo[p]; j<=lf->evs.hi[p]; j++)
+{ sv += prwt(&lf->lfd,(int)pi[j]);
+for (i=0; i<d; i++)
+evptx(&lf->fp,nv,i) += datum(&lf->lfd,i,pi[j])*prwt(&lf->lfd,(int)pi[j]);
+}
+for (i=0; i<d; i++) evptx(&lf->fp,nv,i) /= sv;
+lf->lfd.n = lf->evs.hi[p] - lf->evs.lo[p] + 1;
+des->ind = &pi[lf->evs.lo[p]]; /* why? */
+PROC_VERTEX(des,lf,nv);
+lf->lfd.n = n; des->ind = pi;
+nv++;
+}
+p++;
+}
+/* We've built the tree. Now do the fitting. */
+if (ev(&lf->evs)==EKDTR)
+for (i=0; i<nv; i++) PROC_VERTEX(des,lf,i);
+lf->evs.nce = nc; lf->fp.nv = nv;
+return;
+}
+void newcell(nv,vc,xev, d, k, split_val, cpar, clef, crig)
+double *xev, split_val;
+int *cpar, *clef, *crig;
+int *nv, vc, d, k;
+{ int i, ii, j, j2, tk, match;
+tk = 1<<k;
+for (i=0; i<vc; i++)
+{ if ((i&tk) == 0)
+{ for (j=0; j<d; j++) xev[*nv*d+j] = xev[d*cpar[i]+j];
+xev[*nv*d+k] = split_val;
+match = 0; j = vc; /* no matches in first vc points */
+while ((j<*nv) && (!match))
+{ j2 = 0;
+while ((j2<d) && (xev[*nv*d+j2] == xev[j*d+j2])) j2++;
+match = (j2==d);
+if (!match) j++;
+}
+ii = i+tk;
+clef[i] = cpar[i];
+clef[ii]= crig[i] = j;
+crig[ii]= cpar[ii];
+if (!match) (*nv)++;
+}
+}
+return;
+}
+extern void hermite2();
+double blend(fp,evs,s,x,ll,ur,j,nt,t,what)
+fitpt *fp;
+evstruc *evs;
+double s, *x, *ll, *ur;
+int j, nt, *t, what;
+{
+int *ce, k, k1, m, nc, j0, j1;
+double v0, v1, xibar, g0[3], g1[3], gg[4], gp[4], phi[4];
+ce = evs->ce;
+for (k=0; k<4; k++)  /* North South East West */
+{ k1 = (k>1);
+v0 = ll[k1]; v1 = ur[k1];
+j0 = ce[j+2*(k==0)+(k==2)];
+j1 = ce[j+3-2*(k==1)-(k==3)];
+xibar = (k%2==0) ? ur[k<2] : ll[k<2];
+m = nt;
+while ((m>=0) && ((evs->s[t[m]] != (k<=1)) | (evs->sv[t[m]] != xibar))) m--;
+if (m >= 0)
+{ m = (k%2==1) ? evs->lo[t[m]] : evs->hi[t[m]];
+while (evs->s[m] != -1)
+m = (x[evs->s[m]] < evs->sv[m]) ? evs->lo[m] : evs->hi[m];
+if (v0 < evptx(fp,ce[4*m+2*(k==1)+(k==3)],k1))
+{ j0 = ce[4*m+2*(k==1)+(k==3)];
+v0 = evptx(fp,j0,k1);
+}
+if (evptx(fp,ce[4*m+3-2*(k==0)-(k==2)],k1) < v1)
+{ j1 = ce[4*m+3-2*(k==0)-(k==2)];
+v1 = evptx(fp,j1,k1);
+}
+}
+nc = exvval(fp,g0,j0,2,what,0);
+nc = exvval(fp,g1,j1,2,what,0);
+if (nc==1)
+gg[k] = linear_interp((x[(k>1)]-v0),v1-v0,g0[0],g1[0]);
+else
+{ hermite2(x[(k>1)]-v0,v1-v0,phi);
+gg[k] = phi[0]*g0[0]+phi[1]*g1[0]+(phi[2]*g0[1+k1]+phi[3]*g1[1+k1])*(v1-v0);
+gp[k] = phi[0]*g0[2-k1] + phi[1]*g1[2-k1];
+}
+}
+s = -s;
+if (nc==1)
+for (k=0; k<2; k++)
+s += linear_interp(x[k]-ll[k],ur[k]-ll[k],gg[3-2*k],gg[2-2*k]);
+else
+for (k=0; k<2; k++) /* EW NS */
+{ hermite2(x[k]-ll[k],ur[k]-ll[k],phi);
+s += phi[0]*gg[3-2*k] + phi[1]*gg[2-2*k]
++(phi[2]*gp[3-2*k] + phi[3]*gp[2-2*k]) * (ur[k]-ll[k]);
+}
+return(s);
+}
+double kdtre_int(fp,evs,x,what)
+fitpt *fp;
+evstruc *evs;
+double *x;
+int what;
+{
+int *ce, k, vc, t[20], nt, nc, j, d;
+double *ll, *ur, ff, vv[64][64];
+d = fp->d;
+vc = 1<<d;
+if (d > 6) { LERR(("d too large in kdint")); return(NOSLN); }
+/* descend the tree to find the terminal cell */
+nt = 0; t[nt] = 0; k = 0;
+while (evs->s[k] != -1)
+{ nt++;
+if (nt>=20) { LERR(("Too many levels in kdint")); return(NOSLN); }
+k = t[nt] = (x[evs->s[k]] < evs->sv[k]) ? evs->lo[k] : evs->hi[k];
+}
+ce = &evs->ce[k*vc];
+ll = evpt(fp,ce[0]);
+ur = evpt(fp,ce[vc-1]);
+nc = 0;
+for (j=0; j<vc; j++) nc = exvval(fp,vv[j],(int)ce[j],d,what,0);
+ff = rectcell_interp(x,vv,ll,ur,d,nc);
+if (d==2) ff = blend(fp,evs,ff,x,ll,ur,k*vc,nt,t,what);
+return(ff);
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+/*
+* convert eval. structure string to numeric code.
+*/
+#define NETYPE 11
+static char *etype[NETYPE]= { "tree",     "phull", "data", "grid", "kdtree",
+"kdcenter", "cross", "preset", "xbar", "none",
+"sphere" };
+static int   evals[NETYPE]= { ETREE, EPHULL, EDATA, EGRID, EKDTR,
+EKDCE, ECROS,  EPRES, EXBAR, ENONE, ESPHR };
+int lfevstr(char *z)
+{ return(pmatch(z, etype, evals, NETYPE, -1));
+}
+void evstruc_init(evs)
+evstruc *evs;
+{ int i;
+ev(evs) = ETREE;
+mk(evs) = 100;
+cut(evs) = 0.8;
+for (i=0; i<MXDIM; i++)
+{ evs->fl[i] = evs->fl[i+MXDIM] = 0.0;
+evs->mg[i] = 10;
+}
+evs->nce = evs->ncm = 0;
+}
+int evstruc_reqi(nvm,ncm,vc)
+int nvm, ncm, vc;
+{ return(ncm*vc+3*MAX(ncm,nvm));
+}
+/* al=1: allows dynamic allocation.
+* al=0: inhibit. use with caution.
+*/
+void evstruc_alloc(evs,nvm,ncm,vc,al)
+evstruc *evs;
+int nvm, ncm, vc, al;
+{ int rw, *k;
+if (al)
+{ rw = evstruc_reqi(nvm,ncm,vc);
+if (evs->liw<rw)
+{ evs->iwk = (int *)calloc(rw,sizeof(int));
+if ( evs->iwk == NULL ) {
+printf("Problem allocating memory for evs->iwk\n");fflush(stdout);
+}
+evs->liw = rw;
+}
+}
+k = evs->iwk;
+evs->ce = k; k += vc*ncm;
+evs->s  = k; k += MAX(ncm,nvm);
+evs->lo = k; k += MAX(ncm,nvm);
+evs->hi = k; k += MAX(ncm,nvm);
+}
+void guessnv(evs,sp,mdl,n,d,lw,nvc)
+evstruc *evs;
+smpar *sp;
+module *mdl;
+int n, d, *lw, *nvc;
+{ int i, nvm, ncm, vc;
+npar(sp) = calcp(sp,d);
+switch(ev(evs))
+{ case ETREE:
+atree_guessnv(evs,&nvm,&ncm,&vc,d,nn(sp));
+break;
+case EPHULL:
+nvm = ncm = mk(evs)*d;
+vc = d+1;
+break;
+case EDATA:
+case ECROS:
+nvm = n;
+ncm = vc = 0;
+break;
+case EKDTR:
+case EKDCE:
+kdtre_guessnv(evs,&nvm,&ncm,&vc,n,d,nn(sp));
+break;
+case EGRID:
+nvm = 1;
+for (i=0; i<d; i++) nvm *= evs->mg[i];
+ncm = 0;
+vc = 1<<d;
+break;
+case EXBAR:
+case ENONE:
+nvm = 1;
+ncm = vc = 0;
+break;
+case EPRES:
+nvm = evs->mg[0];
+ncm = vc = 0;
+break;
+default:
+LERR(("guessnv: I don't know this evaluation structure."));
+nvm = ncm = vc = 0;
+}
+lw[0] = des_reqd(n,npar(sp));
+lw[1] = lfit_reqd(d,nvm,ncm,mdl);
+lw[2] = evstruc_reqi(nvm,ncm,vc);
+lw[6] = des_reqi(n,npar(sp));
+lw[3] = pc_reqd(d,npar(sp));
+lw[4] = mdl->keepv;
+lw[5] = mdl->keepc;
+if (nvc==NULL) return;
+nvc[0] = nvm;
+nvc[1] = ncm;
+nvc[2] = vc;
+nvc[3] = nvc[4] = 0;
+}
+/*
+* trchck checks the working space on the lfit structure
+* has space for nvm vertices and ncm cells.
+*/
+void lfit_alloc(lf)
+lfit *lf;
+{ lf->fp.lwk = lf->fp.lev = lf->evs.liw = lf->pc.lwk = 0;
+lf->lf_init_id = LF_INIT_ID;
+}
+int lfit_reqd(d,nvm,ncm,mdl)
+int d, nvm, ncm;
+module *mdl;
+{ int z;
+z = mdl->keepv;
+return(nvm*z+ncm);
+}
+void trchck(lf,nvm,ncm,vc)
+lfit *lf;
+int nvm, ncm, vc;
+{ int rw, d, *k;
+double *z;
+if (lf->lf_init_id != LF_INIT_ID) lfit_alloc(lf);
+lf->fp.nvm = nvm; lf->evs.ncm = ncm;
+d = lf->lfd.d;
+if (lf->fp.lev < d*nvm)
+{ lf->fp.xev = (double *)calloc(d*nvm,sizeof(double));
+if ( lf->fp.xev == NULL ) {
+printf("Problem allocating memory for lf->fp.xev\n");fflush(stdout);
+}
+lf->fp.lev = d*nvm;
+}
+rw = lfit_reqd(d,nvm,ncm,&lf->mdl);
+if (lf->fp.lwk < rw)
+{
+lf->fp.coef = (double *)calloc(rw,sizeof(double));
+if ( lf->fp.coef == NULL ) {
+printf("Problem allocating memory for lf->fp.coef\n");fflush(stdout);
+}
+lf->fp.lwk = rw;
+}
+z = lf->fp.wk = lf->fp.coef;
+lf->fp.h = NULL;
+if (!lf->mdl.isset) mut_printf("module not set.\n");
+else
+{ if (lf->mdl.alloc!=NULL) lf->mdl.alloc(lf);
+z += KEEPV(lf) * nvm;
+}
+lf->evs.sv = z; z += ncm;
+evstruc_alloc(&lf->evs,nvm,ncm,vc,1);
+}
+void data_guessnv(nvm,ncm,vc,n)
+int *nvm, *ncm, *vc, n;
+{ *nvm = n;
+*ncm = *vc = 0;
+}
+void dataf(des,lf)
+design *des;
+lfit *lf;
+{
+int d, i, j, ncm, nv, vc;
+d = lf->lfd.d;
+data_guessnv(&nv,&ncm,&vc,lf->lfd.n);
+trchck(lf,nv,ncm,vc);
+for (i=0; i<nv; i++)
+for (j=0; j<d; j++) evptx(&lf->fp,i,j) = datum(&lf->lfd,j,i);
+for (i=0; i<nv; i++)
+{
+PROC_VERTEX(des,lf,i);
+lf->evs.s[i] = 0;
+}
+lf->fp.nv = lf->fp.nvm = nv; lf->evs.nce = 0;
+}
+void xbar_guessnv(nvm,ncm,vc)
+int *nvm, *ncm, *vc;
+{ *nvm = 1;
+*ncm = *vc = 0;
+return;
+}
+void xbarf(des,lf)
+design *des;
+lfit *lf;
+{ int i, d, nvm, ncm, vc;
+d = lf->lfd.d;
+xbar_guessnv(&nvm,&ncm,&vc);
+trchck(lf,1,0,0);
+for (i=0; i<d; i++) evptx(&lf->fp,0,i) = lf->pc.xbar[i];
+PROC_VERTEX(des,lf,0);
+lf->evs.s[0] = 0;
+lf->fp.nv = 1; lf->evs.nce = 0;
+}
+void preset(des,lf)
+design *des;
+lfit *lf;
+{ int i, nv;
+nv = lf->fp.nvm;
+trchck(lf,nv,0,0);
+for (i=0; i<nv; i++)
+{
+PROC_VERTEX(des,lf,i);
+lf->evs.s[i] = 0;
+}
+lf->fp.nv = nv; lf->evs.nce = 0;
+}
+void crossf(des,lf)
+design *des;
+lfit *lf;
+{ int d, i, j, n, nv, ncm, vc;
+double w;
+n = lf->lfd.n; d = lf->lfd.d;
+data_guessnv(&nv,&ncm,&vc,n);
+trchck(lf,nv,ncm,vc);
+if (lf->lfd.w==NULL) LERR(("crossf() needs prior weights"));
+for (i=0; i<n; i++)
+for (j=0; j<d; j++) evptx(&lf->fp,i,j) = datum(&lf->lfd,j,i);
+for (i=0; i<n; i++)
+{ lf->evs.s[i] = 0;
+w = prwt(&lf->lfd,i);
+lf->lfd.w[i] = 0;
+PROC_VERTEX(des,lf,i);
+lf->lfd.w[i] = w;
+}
+lf->fp.nv = n; lf->evs.nce = 0;
+}
+void gridf(des,lf)
+design *des;
+lfit *lf;
+{ int d, i, j, nv, u0, u1, z;
+nv = 1; d = lf->lfd.d;
+for (i=0; i<d; i++)
+{ if (lf->evs.mg[i]==0)
+lf->evs.mg[i] = 2+(int)((lf->evs.fl[i+d]-lf->evs.fl[i])/(lf->lfd.sca[i]*cut(&lf->evs)));
+nv *= lf->evs.mg[i];
+}
+trchck(lf,nv,0,1<<d);
+for (i=0; i<nv; i++)
+{ z = i;
+for (j=0; j<d; j++)
+{ u0 = z%lf->evs.mg[j];
+u1 = lf->evs.mg[j]-1-u0;
+evptx(&lf->fp,i,j) = (lf->evs.mg[j]==1) ? lf->evs.fl[j] :
+(u1*lf->evs.fl[j]+u0*lf->evs.fl[j+d])/(lf->evs.mg[j]-1);
+z = z/lf->evs.mg[j];
+}
+lf->evs.s[i] = 0;
+PROC_VERTEX(des,lf,i);
+}
+lf->fp.nv = nv; lf->evs.nce = 0;
+}
+int findpt(fp,evs,i0,i1)
+fitpt *fp;
+evstruc *evs;
+int i0, i1;
+{ int i;
+if (i0>i1) ISWAP(i0,i1);
+for (i=i1+1; i<fp->nv; i++)
+if ((evs->lo[i]==i0) && (evs->hi[i]==i1)) return(i);
+return(-1);
+}
+/*
+add a new vertex at the midpoint of (x[i0],x[i1]).
+return the vertex number.
+*/
+int newsplit(des,lf,i0,i1,pv)
+design *des;
+lfit *lf;
+int i0, i1, pv;
+{ int i, nv;
+i = findpt(&lf->fp,&lf->evs,i0,i1);
+if (i>=0) return(i);
+if (i0>i1) ISWAP(i0,i1);
+nv = lf->fp.nv;
+/* the point is new. Now check we have space for the new point. */
+if (nv>=lf->fp.nvm)
+{
+LERR(("newsplit: out of vertex space"));
+return(-1);
+}
+/* compute the new point, and evaluate the fit */
+lf->evs.lo[nv] = i0;
+lf->evs.hi[nv] = i1;
+for (i=0; i<lf->fp.d; i++)
+evptx(&lf->fp,nv,i) = (evptx(&lf->fp,i0,i)+evptx(&lf->fp,i1,i))/2;
+if (pv) /* pseudo vertex */
+{ lf->fp.h[nv] = (lf->fp.h[i0]+lf->fp.h[i1])/2;
+lf->evs.s[nv] = 1; /* pseudo-vertex */
+}
+else /* real vertex */
+{
+PROC_VERTEX(des,lf,nv);
+lf->evs.s[nv] = 0;
+}
+lf->fp.nv++;
+return(nv);
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+/*
+*   Functions for constructing the fit and
+*   interpolating on the circle/sphere. d=2 only.
+*/
+#include "lfev.h"
+/*
+* Guess the number of fitting points.
+*/
+void sphere_guessnv(nvm,ncm,vc,mg)
+int *nvm, *ncm, *vc, *mg;
+{ *nvm = mg[1]*(mg[0]+1);
+*ncm = 0;
+*vc = 0;
+}
+void sphere_start(des,lf)
+design *des;
+lfit *lf;
+{ int d, i, j, ct, nv, ncm, vc, *mg;
+double rmin, rmax, *orig, r, th, c, s;
+mg = mg(&lf->evs);
+sphere_guessnv(&nv,&ncm,&vc,mg);
+trchck(lf,nv,0,0);
+d = lf->lfd.d;
+rmin = lf->evs.fl[0];
+rmax = lf->evs.fl[1];
+orig = &lf->evs.fl[2];
+rmin = 0; rmax = 1; orig[0] = orig[1] = 0.0;
+ct = 0;
+for (i=0; i<mg[1]; i++)
+{ th = 2*PI*i/mg[1];
+c = cos(th);
+s = sin(th);
+for (j=0; j<=mg[0]; j++)
+{ r = rmin + (rmax-rmin)*j/mg[0];
+evptx(&lf->fp,ct,0) = orig[0] + r*c;
+evptx(&lf->fp,ct,1) = orig[1] + r*s;
+PROC_VERTEX(des,lf,ct);
+ct++;
+}
+}
+lf->fp.nv = ct;
+lf->evs.nce = 0;
+}
+double sphere_int(lf,x,what)
+lfit *lf;
+double *x;
+int what;
+{ double rmin, rmax, *orig, dx, dy, r, th, th0, th1;
+double v[64][64], c0, c1, s0, s1, r0, r1, d0, d1;
+double ll[2], ur[2], xx[2];
+int i0, j0, i1, j1, *mg, nc, ce[4];
+rmin = lf->evs.fl[0];
+rmax = lf->evs.fl[1];
+orig = &lf->evs.fl[2];
+rmin = 0; rmax = 1; orig[0] = orig[1] = 0.0;
+mg = mg(&lf->evs);
+dx = x[0] - orig[0];
+dy = x[1] - orig[1];
+r = sqrt(dx*dx+dy*dy);
+th = atan2(dy,dx); /* between -pi and pi */
+i0 = (int)floor(mg[1]*th/(2*PI)) % mg[1];
+j0 = (int)(mg[0]*(r-rmin)/(rmax-rmin));
+i1 = (i0+1) % mg[1];
+j1 = j0+1; if (j1>mg[0]) { j0 = mg[0]-1; j1 = mg[0]; }
+ce[0] = i0*(mg[0]+1)+j0;
+ce[1] = i0*(mg[0]+1)+j1;
+ce[2] = i1*(mg[0]+1)+j0;
+ce[3] = i1*(mg[0]+1)+j1;
+nc = exvval(&lf->fp,v[0],ce[0],2,what,1);
+nc = exvval(&lf->fp,v[1],ce[1],2,what,1);
+nc = exvval(&lf->fp,v[2],ce[2],2,what,1);
+nc = exvval(&lf->fp,v[3],ce[3],2,what,1);
+th0 = 2*PI*i0/mg[1]; c0 = cos(th0); s0 = sin(th0);
+th1 = 2*PI*i1/mg[1]; c1 = cos(th1); s1 = sin(th1);
+r0 = rmin + j0*(rmax-rmin)/mg[0];
+r1 = rmin + j1*(rmax-rmin)/mg[0];
+d0 = c0*v[0][1] + s0*v[0][2];
+d1 = r0*(c0*v[0][2]-s0*v[0][1]);
+v[0][1] = d0; v[0][2] = d1;
+d0 = c0*v[1][1] + s0*v[1][2];
+d1 = r1*(c0*v[1][2]-s0*v[1][1]);
+v[1][1] = d0; v[1][2] = d1;
+d0 = c1*v[2][1] + s1*v[2][2];
+d1 = r0*(c1*v[2][2]-s1*v[2][1]);
+v[2][1] = d0; v[2][2] = d1;
+d0 = c1*v[3][1] + s1*v[3][2];
+d1 = r1*(c1*v[3][2]-s1*v[3][1]);
+v[3][1] = d0; v[3][2] = d1;
+xx[0] = r; xx[1] = th;
+ll[0] = r0; ll[1] = th0;
+ur[0] = r1; ur[1] = th1;
+return(rectcell_interp(xx,v,ll,ur,2,nc));
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+void solve(A,b,d) /* this is crude! A organized by column. */
+double *A, *b;
+int d;
+{ int i, j, k;
+double piv;
+for (i=0; i<d; i++)
+{ piv = A[(d+1)*i];
+for (j=i; j<d; j++) A[j*d+i] /= piv;
+b[i] /= piv;
+for (j=0; j<d; j++) if (j != i)
+{ piv = A[i*d+j];
+A[i*d+j] = 0;
+for (k=i+1; k<d; k++)
+A[k*d+j] -= piv*A[k*d+i];
+b[j] -= piv*b[i];
+}
+}
+}
+void triang_guessnv(nvm,ncm,vc,d,mk)
+int *nvm, *ncm, *vc, d, mk;
+{ *nvm = *ncm = mk*d;
+*vc = d+1;
+return;
+}
+int triang_split(lf,ce,le)
+lfit *lf;
+double *le;
+int *ce;
+{ int d, i, j, k, nts, vc;
+double di, dfx[MXDIM];
+nts = 0; d = lf->fp.d; vc = d+1;
+for (i=0; i<d; i++)
+for (j=i+1; j<=d; j++)
+{ for (k=0; k<d; k++)
+dfx[k] = evptx(&lf->fp,ce[i],k)-evptx(&lf->fp,ce[j],k);
+di = rho(dfx,lf->lfd.sca,d,KSPH,NULL);
+le[i*vc+j] = le[j*vc+i] = di/MIN(lf->fp.h[ce[i]],lf->fp.h[ce[j]]);
+nts = nts || le[i*vc+j]>cut(&lf->evs);
+}
+return(nts);
+}
+void resort(pv,xev,dig)
+double *xev;
+int *pv, *dig;
+{ double d0, d1, d2;
+int i;
+d0 = d1 = d2 = 0;
+for (i=0; i<3; i++)
+{ d0 += (xev[3*pv[11]+i]-xev[3*pv[1]+i])*(xev[3*pv[11]+i]-xev[3*pv[1]+i]);
+d1 += (xev[3*pv[ 7]+i]-xev[3*pv[2]+i])*(xev[3*pv[ 7]+i]-xev[3*pv[2]+i]);
+d2 += (xev[3*pv[ 6]+i]-xev[3*pv[3]+i])*(xev[3*pv[ 6]+i]-xev[3*pv[3]+i]);
+}
+if ((d0<=d1) & (d0<=d2))
+{ dig[0] = pv[1]; dig[1] = pv[11];
+dig[2] = pv[2]; dig[3] = pv[7];
+dig[4] = pv[3]; dig[5] = pv[6];
+}
+else if (d1<=d2)
+{ dig[0] = pv[2]; dig[1] = pv[7];
+dig[2] = pv[1]; dig[3] = pv[11];
+dig[4] = pv[3]; dig[5] = pv[6];
+}
+else
+{ dig[0] = pv[3]; dig[1] = pv[6];
+dig[2] = pv[2]; dig[3] = pv[7];
+dig[4] = pv[1]; dig[5] = pv[11];
+}
+}
+void triang_grow(des,lf,ce,ct,term)
+design *des;
+lfit *lf;
+int *ce, *ct, *term;
+{ double le[(1+MXDIM)*(1+MXDIM)], ml;
+int d, i, j, im, jm, vc, pv[(1+MXDIM)*(1+MXDIM)], dig[6];
+int nce[1+MXDIM];
+if (lf_error) return;
+d = lf->fp.d; vc = d+1;
+if (!triang_split(lf,ce,le))
+{ if (ct != NULL)
+{ for (i=0; i<vc; i++) term[*ct*vc+i] = ce[i];
+(*ct)++;
+}
+return;
+}
+if (d>3)
+{ ml = 0;
+for (i=0; i<d; i++)
+for (j=i+1; j<vc; j++)
+if (le[i*vc+j]>ml) { ml = le[i*vc+j]; im = i; jm = j; }
+pv[0] = newsplit(des,lf,(int)ce[im],(int)ce[jm],0);
+for (i=0; i<vc; i++) nce[i] = ce[i];
+nce[im] = pv[0]; triang_grow(des,lf,nce,ct,term); nce[im] = ce[im];
+nce[jm] = pv[0]; triang_grow(des,lf,nce,ct,term);
+return;
+}
+for (i=0; i<d; i++)
+for (j=i+1; j<=d; j++)
+pv[i*vc+j] = pv[j*vc+i]
+= newsplit(des,lf,(int)ce[i],(int)ce[j],le[i*vc+j]<=cut(&lf->evs));
+for (i=0; i<=d; i++) /* corners */
+{ for (j=0; j<=d; j++) nce[j] = (j==i) ? ce[i] : pv[i*vc+j];
+triang_grow(des,lf,nce,ct,term);
+}
+if (d==2) /* center for d=2 */
+{ nce[0] = pv[5]; nce[1] = pv[2]; nce[2] = pv[1];
+triang_grow(des,lf,nce,ct,term);
+}
+if (d==3) /* center for d=3 */
+{ resort(pv,evp(&lf->fp),dig);
+nce[0] = dig[0]; nce[1] = dig[1];
+nce[2] = dig[2]; nce[3] = dig[4]; triang_grow(des,lf,nce,ct,term);
+nce[2] = dig[5]; nce[3] = dig[3]; triang_grow(des,lf,nce,ct,term);
+nce[2] = dig[2]; nce[3] = dig[5]; triang_grow(des,lf,nce,ct,term);
+nce[2] = dig[4]; nce[3] = dig[3]; triang_grow(des,lf,nce,ct,term);
+}
+}
+void triang_descend(tr,xa,ce)
+lfit *tr;
+double *xa;
+int *ce;
+{ double le[(1+MXDIM)*(1+MXDIM)], ml;
+int d, vc, i, j, im, jm, pv[(1+MXDIM)*(1+MXDIM)];
+design *des;
+des = NULL;
+if (!triang_split(tr,ce,le)) return;
+d = tr->fp.d; vc = d+1;
+if (d>3) /* split longest edge */
+{ ml = 0;
+for (i=0; i<d; i++)
+for (j=i+1; j<vc; j++)
+if (le[i*vc+j]>ml) { ml = le[i*vc+j]; im = i; jm = j; }
+pv[0] = newsplit(des,tr,(int)ce[im],(int)ce[jm],0);
+if (xa[im]>xa[jm])
+{ xa[im] -= xa[jm]; xa[jm] *= 2; ce[jm] = pv[0]; }
+else
+{ xa[jm] -= xa[im]; xa[im] *= 2; ce[im] = pv[0]; }
+triang_descend(tr,xa,ce);
+return;
+}
+for (i=0; i<d; i++)
+for (j=i+1; j<=d; j++)
+pv[i*vc+j] = pv[j*vc+i]
+= newsplit(des,tr,(int)ce[i],(int)ce[j],le[i*d+j]<=cut(&tr->evs));
+for (i=0; i<=d; i++) if (xa[i]>=0.5) /* in corner */
+{ for (j=0; j<=d; j++)
+{ if (i!=j) ce[j] = pv[i*vc+j];
+xa[j] = 2*xa[j];
+}
+xa[i] -= 1;
+triang_descend(tr,xa,ce);
+return;
+}
+if (d==1) { LERR(("weights sum to < 1")); }
+if (d==2) /* center */
+{ ce[0] = pv[5]; xa[0] = 1-2*xa[0];
+ce[1] = pv[2]; xa[1] = 1-2*xa[1];
+ce[2] = pv[1]; xa[2] = 1-2*xa[2];
+triang_descend(tr,xa,ce);
+}
+if (d==3) /* center */
+{ double z; int dig[6];
+resort(pv,evp(&tr->fp),dig);
+ce[0] = dig[0]; ce[1] = dig[1];
+xa[0] *= 2; xa[1] *= 2; xa[2] *= 2; xa[3] *= 2;
+if (xa[0]+xa[2]>=1)
+{ if (xa[0]+xa[3]>=1)
+{ ce[2] = dig[2]; ce[3] = dig[4];
+z = xa[0];
+xa[3] += z-1; xa[2] += z-1; xa[0] = xa[1]; xa[1] = 1-z;
+}
+else
+{ ce[2] = dig[2]; ce[3] = dig[5];
+z = xa[3]; xa[3] = xa[1]+xa[2]-1; xa[1] = z;
+z = xa[2]; xa[2] += xa[0]-1; xa[0] = 1-z;
+} }
+else
+{ if (xa[1]+xa[2]>=1)
+{ ce[2] = dig[5]; ce[3] = dig[3];
+xa[1] = 1-xa[1]; xa[2] -= xa[1]; xa[3] -= xa[1];
+}
+else
+{ ce[2] = dig[4]; ce[3] = dig[3];
+z = xa[3]; xa[3] += xa[1]-1; xa[1] = xa[2];
+xa[2] = z+xa[0]-1; xa[0] = 1-z;
+} }
+triang_descend(tr,xa,ce);
+} }
+void covrofdata(lfd,V,mn) /* covar of data; mean in mn */
+lfdata *lfd;
+double *V, *mn;
+{ int d, i, j, k;
+double s;
+s = 0; d = lfd->d;
+for (i=0; i<d*d; i++) V[i] = 0;
+for (i=0; i<lfd->n; i++)
+{ s += prwt(lfd,i);
+for (j=0; j<d; j++)
+for (k=0; k<d; k++)
+V[j*d+k] += prwt(lfd,i)*(datum(lfd,j,i)-mn[j])*(datum(lfd,k,i)-mn[k]);
+}
+for (i=0; i<d*d; i++) V[i] /= s;
+}
+int intri(x,w,xev,xa,d) /* is x in triangle bounded by xd[0..d-1]? */
+double *x, *xev, *xa;
+int *w, d;
+{ int i, j;
+double eps, *r, xd[MXDIM*MXDIM];
+eps = 1.0e-10;
+r = &xev[w[d]*d];
+for (i=0; i<d; i++)
+{ xa[i] = x[i]-r[i];
+for (j=0; j<d; j++) xd[i*d+j] = xev[w[i]*d+j]-r[j];
+}
+solve(xd,xa,d);
+xa[d] = 1.0;
+for (i=0; i<d; i++) xa[d] -= xa[i];
+for (i=0; i<=d; i++) if ((xa[i]<-eps) | (xa[i]>1+eps)) return(0);
+return(1);
+}
+void triang_start(des,lf) /* Triangulation with polyhedral start */
+design *des;
+lfit *lf;
+{
+int i, j, k, n, d, nc, nvm, ncm, vc;
+int *ce, ed[1+MXDIM];
+double V[MXDIM*MXDIM], P[MXDIM*MXDIM], sigma, z[MXDIM], xa[1+MXDIM], *xev;
+xev = evp(&lf->fp);
+d = lf->lfd.d; n = lf->lfd.n;
+lf->fp.nv = nc = 0;
+triang_guessnv(&nvm,&ncm,&vc,d,mk(&lf->evs));
+trchck(lf,nvm,ncm,vc);
+ce = lf->evs.ce;
+for (j=0; j<d; j++) xev[j] = lf->pc.xbar[j];
+lf->fp.nv = 1;
+covrofdata(&lf->lfd,V,xev); /* fix this with scaling */
+eig_dec(V,P,d);
+for (i=0; i<d; i++) /* add vertices +- 2sigma*eigenvect */
+{ sigma = sqrt(V[i*(d+1)]);
+for (j=0; j<d; j++)
+xev[lf->fp.nv*d+j] = xev[j]-2*sigma*P[j*d+i];
+lf->fp.nv++;
+for (j=0; j<d; j++)
+xev[lf->fp.nv*d+j] = xev[j]+2*sigma*P[j*d+i];
+lf->fp.nv++;
+}
+for (i=0; i<n; i++) /* is point i inside? */
+{ ed[0] = 0;
+for (j=0; j<d; j++)
+{ z[j] = 0;
+for (k=0; k<d; k++) z[j] += P[k*d+j]*(datum(&lf->lfd,k,i)-xev[k]);
+ed[j+1] = 2*j+1+(z[j]>0);
+for (k=0; k<d; k++) z[j] = datum(&lf->lfd,j,i);
+}
+k = intri(z,ed,xev,xa,d);
+if (xa[0]<0)
+{ for (j=1; j<=d; j++)
+for (k=0; k<d; k++)
+xev[ed[j]*d+k] = xa[0]*xev[k]+(1-xa[0])*xev[ed[j]*d+k];
+}
+}
+nc = 1<<d; /* create initial cells */
+for (i=0; i<nc; i++)
+{ ce[i*vc] = 0; k = i;
+for (j=0; j<d; j++)
+{ ce[i*vc+j+1] = 2*j+(k%2)+1;
+k>>=1;
+}
+}
+for (i=0; i<lf->fp.nv; i++)
+{ PROC_VERTEX(des,lf,i);
+if (lf_error) return;
+lf->evs.s[i] = 0;
+}
+for (i=0; i<nc; i++)
+triang_grow(des,lf,&ce[i*vc],NULL,NULL);
+lf->evs.nce = nc;
+}
+double triang_cubicint(v,vv,w,d,nc,xxa)
+double *v, *vv, *xxa;
+int *w, d, nc;
+{ double sa, lb, *vert0, *vert1, *vals0, *vals1, deriv0, deriv1;
+int i, j, k;
+if (nc==1) /* linear interpolate */
+{ sa = 0;
+for (i=0; i<=d; i++) sa += xxa[i]*vv[i];
+return(sa);
+}
+sa = 1.0;
+for (j=d; j>0; j--)  /* eliminate v[w[j]] */
+{ lb = xxa[j]/sa;
+for (k=0; k<j; k++) /* Interpolate edge v[w[k]],v[w[j]] */
+{ vert0 = &v[w[k]*d];
+vert1 = &v[w[j]*d];
+vals0 = &vv[k*nc];
+vals1 = &vv[j*nc];
+deriv0 = deriv1 = 0;
+for (i=0; i<d; i++)
+{ deriv0 += (vert1[i]-vert0[i])*vals0[i+1];
+deriv1 += (vert1[i]-vert0[i])*vals1[i+1];
+}
+vals0[0] = cubic_interp(lb,vals0[0],vals1[0],deriv0,deriv1);
+for (i=1; i<=d; i++)
+vals0[i] = (1-lb)*((1-lb)*vals0[i]+lb*vals1[i]);
+}
+sa -= xxa[j];
+if (sa<=0) j = 0;
+}
+return(vals0[0]);
+}
+double triang_clotoch(xev,vv,ce,p,xxa)
+double *xev, *vv, *xxa;
+int *ce, p;
+{ double cfo[3], cfe[3], cg[9], *va, *vb, *vc,
+l0, nm[3], na, nb, nc, *xl, *xr, *xz, d0, d1, lb, dlt, gam;
+int i, w[3], cfl, cfr;
+if (p==1)
+return(xxa[0]*vv[0]+xxa[1]*vv[1]+xxa[2]*vv[2]);
+if (xxa[2]<=MIN(xxa[0],xxa[1]))
+{ va = &xev[2*ce[0]]; vb = &xev[2*ce[1]]; vc = &xev[2*ce[2]];
+w[0] = 0; w[1] = 3; w[2] = 6;
+}
+else
+if (xxa[1]<xxa[0])
+{ w[0] = 0; w[1] = 6; w[2] = 3;
+va = &xev[2*ce[0]]; vb = &xev[2*ce[2]]; vc = &xev[2*ce[1]];
+lb = xxa[1]; xxa[1] = xxa[2]; xxa[2] = lb;
+}
+else
+{ w[0] = 6; w[1] = 3; w[2] = 0;
+va = &xev[2*ce[2]]; vb = &xev[2*ce[1]]; vc = &xev[2*ce[0]];
+lb = xxa[0]; xxa[0] = xxa[2]; xxa[2] = lb;
+}
+/* set cg to values and derivatives on standard triangle */
+for (i=0; i<3; i++)
+{ cg[3*i] = vv[w[i]];
+cg[3*i+1] = ((vb[0]-va[0])*vv[w[i]+1]
++(vb[1]-va[1])*vv[w[i]+2])/2;  /* df/dx */
+cg[3*i+2] = ((2*vc[0]-vb[0]-va[0])*vv[w[i]+1]
++(2*vc[1]-vb[1]-va[1])*vv[w[i]+2])/2.0; /* sqrt{3} df/dy */
+}
+dlt = (vb[0]-va[0])*(vc[1]-va[1])-(vc[0]-va[0])*(vb[1]-va[1]);
+/* Twice area; +ve if abc antic.wise  -ve is abc c.wise */
+cfo[0] = (cg[0]+cg[3]+cg[6])/3;
+cfo[1] = (2*cg[0]-cg[3]-cg[6])/4;
+cfo[2] = (2*cg[3]-cg[0]-cg[6])/4;
+na = -cg[1]+cg[2];  /* perp. deriv, rel. length 2 */
+nb = -cg[4]-cg[5];
+nc = 2*cg[7];
+cfo[1] += (nb-nc)/16;
+cfo[2] += (nc-na)/16;
+na = -cg[1]-cg[2]/3.0;  /* derivatives back to origin */
+nb =  cg[4]-cg[5]/3.0;
+nc =        cg[8]/1.5;
+cfo[0] -= (na+nb+nc)*7/54;
+cfo[1] += 13*(nb+nc-2*na)/144;
+cfo[2] += 13*(na+nc-2*nb)/144;
+for (i=0; i<3; i++)
+{ /* Outward normals by linear interpolation on original triangle.
+Convert to outward normals on standard triangle.
+Actually, computed to opposite corner */
+switch(i)
+{ case 0: xl = vc; xr = vb; xz = va; cfl = w[2]; cfr = w[1];
+break;
+case 1: xl = va; xr = vc; xz = vb; cfl = w[0]; cfr = w[2];
+break;
+case 2: xl = vb; xr = va; xz = vc; cfl = w[1]; cfr = w[0];
+break;
+}
+na = xr[0]-xl[0]; nb = xr[1]-xl[1];
+lb = na*na+nb*nb;
+d0 = 1.5*(vv[cfr]-vv[cfl]) - 0.25*(na*(vv[cfl+1]+vv[cfr+1])
++nb*(vv[cfl+2]+vv[cfr+2]));
+d1 = 0.5*( na*(vv[cfl+2]+vv[cfr+2])-nb*(vv[cfl+1]+vv[cfr+1]) );
+l0 = (xz[0]-xl[0])*na+(xz[1]-xl[1])*nb-lb/2;
+nm[i] = (d1*dlt-l0*d0)/lb;
+}
+cfo[0] -= (nm[0]+nm[1]+nm[2])*4/81;
+cfo[1] += (2*nm[0]-nm[1]-nm[2])/27;
+cfo[2] += (2*nm[1]-nm[0]-nm[2])/27;
+gam = xxa[0]+xxa[1]-2*xxa[2];
+if (gam==0) return(cfo[0]);
+lb = (xxa[0]-xxa[2])/gam;
+d0 = -2*cg[4]; d1 = -2*cg[1];
+cfe[0] = cubic_interp(lb,cg[3],cg[0],d0,d1);
+cfe[1] = cubintd(lb,cg[3],cg[0],d0,d1);
+cfe[2] = -(1-lb)*(1-2*lb)*cg[5] + 4*lb*(1-lb)*nm[2] - lb*(2*lb-1)*cg[2];
+d0 = 2*(lb*cfo[1]+(1-lb)*cfo[2]);
+d1 = (lb-0.5)*cfe[1]+cfe[2]/3.0;
+return(cubic_interp(gam,cfo[0],cfe[0],d0,d1));
+}
+int triang_getvertexvals(fp,evs,vv,i,what)
+fitpt *fp;
+evstruc *evs;
+double *vv;
+int i, what;
+{ double dx, P, le, vl[1+MXDIM], vh[1+MXDIM];
+int d, il, ih, j, nc;
+d = fp->d;
+if (evs->s[i]==0) return(exvval(fp,vv,i,d,what,0));
+il = evs->lo[i]; nc = triang_getvertexvals(fp,evs,vl,il,what);
+ih = evs->hi[i]; nc = triang_getvertexvals(fp,evs,vh,ih,what);
+vv[0] = (vl[0]+vh[0])/2;
+if (nc==1) return(nc);
+P = 1.5*(vh[0]-vl[0]);
+le = 0.0;
+for (j=0; j<d; j++)
+{ dx = evptx(fp,ih,j)-evptx(fp,il,j);
+vv[0] += dx*(vl[j+1]-vh[j+1])/8;
+vv[j+1] = (vl[j+1]+vh[j+1])/2;
+P -= 1.5*dx*vv[j+1];
+le += dx*dx;
+}
+for (j=0; j<d; j++)
+vv[j+1] += P*(evptx(fp,ih,j)-evptx(fp,il,j))/le;
+return(nc);
+}
+double triang_int(lf,x,what)
+lfit *lf;
+double *x;
+int what;
+{
+int d, i, j, k, vc, nc;
+int *ce, nce[1+MXDIM];
+double xa[1+MXDIM], vv[(1+MXDIM)*(1+MXDIM)], lb;
+fitpt *fp;
+evstruc *evs;
+fp = &lf->fp;
+evs= &lf->evs;
+d = fp->d; vc = d+1;
+ce = evs->ce;
+i = 0;
+while ((i<evs->nce) && (!intri(x,&ce[i*vc],evp(fp),xa,d))) i++;
+if (i==evs->nce) return(NOSLN);
+i *= vc;
+for (j=0; j<vc; j++) nce[j] = ce[i+j];
+triang_descend(lf,xa,nce);
+/* order the vertices -- needed for asymmetric interptr */
+do
+{ k=0;
+for (i=0; i<d; i++)
+if (nce[i]>nce[i+1])
+{ j=nce[i]; nce[i]=nce[i+1]; nce[i+1]=j; k=1;
+lb = xa[i]; xa[i] = xa[i+1]; xa[i+1] = lb;
+}
+} while(k);
+nc = 0;
+for (i=0; i<vc; i++)
+nc =  triang_getvertexvals(fp,evs,&vv[i*nc],nce[i],what);
+return((d==2) ? triang_clotoch(evp(fp),vv,nce,nc,xa) :
+triang_cubicint(evp(fp),vv,nce,d,nc,xa));
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+/*
+* Functions for computing residuals and fitted values from
+* the locfit object.
+*
+* fitted(lf,fit,what,cv,ty) computes fitted values from the
+*   fit structure in lf.
+* resid(y,c,w,th,mi,ty) converts fitted values to residuals
+*/
+#include "lfev.h"
+#define NRT 8
+static char *rtype[NRT] = { "deviance", "d2",    "pearson", "raw",
+"ldot",     "lddot", "fit",     "mean" };
+static int   rvals[NRT] = { RDEV, RDEV2, RPEAR, RRAW, RLDOT, RLDDT, RFIT, RMEAN};
+int restyp(z)
+char *z;
+{ int val;
+val = pmatch(z, rtype, rvals, NRT, -1);
+if (val==-1) LERR(("Unknown type = %s",z));
+return(val);
+}
+double resid(y,w,th,fam,ty,res)
+int fam, ty;
+double y, w, th, *res;
+{ double raw;
+fam = fam & 63;
+if ((fam==TGAUS) | (fam==TROBT) | (fam==TCAUC))
+raw = y-res[ZMEAN];
+else
+raw = y-w*res[ZMEAN];
+switch(ty)
+{ case RDEV:
+if (res[ZDLL]>0) return(sqrt(-2*res[ZLIK]));
+else return(-sqrt(-2*res[ZLIK]));
+case RPEAR:
+if (res[ZDDLL]<=0)
+{ if (res[ZDLL]==0) return(0);
+return(NOSLN);
+}
+return(res[ZDLL]/sqrt(res[ZDDLL]));
+case RRAW:  return(raw);
+case RLDOT: return(res[ZDLL]);
+case RDEV2: return(-2*res[ZLIK]);
+case RLDDT: return(res[ZDDLL]);
+case RFIT:  return(th);
+case RMEAN: return(res[ZMEAN]);
+default: LERR(("resid: unknown residual type %d",ty));
+}
+return(0.0);
+}
+double studentize(res,inl,var,ty,link)
+double res, inl, var, *link;
+int ty;
+{ double den;
+inl *= link[ZDDLL];
+var = var*var*link[ZDDLL];
+if (inl>1) inl = 1;
+if (var>inl) var = inl;
+den = 1-2*inl+var;
+if (den<0) return(0.0);
+switch(ty)
+{ case RDEV:
+case RPEAR:
+case RRAW:
+case RLDOT:
+return(res/sqrt(den));
+case RDEV2:
+return(res/den);
+default: return(res);
+}
+}
+void fitted(lf,fit,what,cv,st,ty)
+lfit *lf;
+double *fit;
+int what, cv, st, ty;
+{ int i, j, d, n, evo;
+double xx[MXDIM], th, inl, var, link[LLEN];
+n = lf->lfd.n;
+d = lf->lfd.d;
+evo = ev(&lf->evs);
+cv &= (evo!=ECROS);
+if ((evo==EDATA)|(evo==ECROS)) evo = EFITP;
+setfamily(&lf->sp);
+for (i=0; i<n; i++)
+{ for (j=0; j<d; j++) xx[j] = datum(&lf->lfd,j,i);
+th = dointpoint(lf,xx,what,evo,i);
+if ((what==PT0)|(what==PVARI)) th = th*th;
+if (what==PCOEF)
+{
+th += base(&lf->lfd,i);
+stdlinks(link,&lf->lfd,&lf->sp,i,th,rsc(&lf->fp));
+if ((cv)|(st))
+{ inl = dointpoint(lf,xx,PT0,evo,i);
+inl = inl*inl;
+if (cv)
+{ th -= inl*link[ZDLL];
+stdlinks(link,&lf->lfd,&lf->sp,i,th,rsc(&lf->fp));
+}
+if (st) var = dointpoint(lf,xx,PNLX,evo,i);
+}
+fit[i] = resid(resp(&lf->lfd,i),prwt(&lf->lfd,i),th,fam(&lf->sp),ty,link);
+if (st) fit[i] = studentize(fit[i],inl,var,ty,link);
+} else fit[i] = th;
+if (lf_error) return;
+}
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+extern double robscale;
+/* special version of ressumm to estimate sigma^2, with derivative estimation */
+void ressummd(lf)
+lfit *lf;
+{ int i;
+double s0, s1;
+s0 = s1 = 0.0;
+if ((fam(&lf->sp)&64)==0)
+{ rv(&lf->fp) = 1.0;
+return;
+}
+for (i=0; i<lf->fp.nv; i++)
+{ s0 += lf->fp.lik[2*lf->fp.nvm+i];
+s1 += lf->fp.lik[i];
+}
+if (s0==0.0)
+rv(&lf->fp) = 0.0;
+else
+rv(&lf->fp) = -2*s1/s0;
+}
+/*
+* res[0] = log-likelihood.
+* res[1] = df0 = tr(H)
+* res[2] = df0 = tr(H'H)
+* res[3] = s^2.
+* res[5] = robustscale.
+*/
+void ressumm(lf,des,res)
+lfit *lf;
+design *des;
+double *res;
+{ int i, j, evo, tg;
+double *oy, pw, r1, r2, rdf, t0, t1, u[MXDIM], link[LLEN];
+fitpt *fp;
+res[0] = res[1] = res[2] = 0.0;
+evo = ev(&lf->evs);
+if ((evo==EKDCE) | (evo==EPRES))
+{ res[3] = 1.0;
+return;
+}
+if (lf->dv.nd>0)
+{ ressummd(lf);
+return;
+}
+r1 = r2 = 0.0;
+if ((evo==EDATA) | (evo==ECROS)) evo = EFITP;
+for (i=0; i<lf->lfd.n; i++)
+{ for (j=0; j<lf->lfd.d; j++) u[j] = datum(&lf->lfd,j,i);
+fitv(des,i) = base(&lf->lfd,i)+dointpoint(lf,u,PCOEF,evo,i);
+des->wd[i] = resp(&(lf->lfd),i) - fitv(des,i);
+wght(des,i) = 1.0;
+des->ind[i] = i;
+}
+tg = fam(&lf->sp);
+res[5] = 1.0;
+if ((tg==TROBT+64) | (tg==TCAUC+64)) /* global robust scale */
+{ oy = lf->lfd.y; lf->lfd.y = des->wd;
+des->xev = lf->pc.xbar;
+locfit(&lf->lfd,des,&lf->sp,1,0,0);
+lf->lfd.y = oy;
+res[5] = robscale;
+}
+for (i=0; i<lf->lfd.n; i++)
+{ for (j=0; j<lf->lfd.d; j++) u[j] = datum(&lf->lfd,j,i);
+t0 = dointpoint(lf,u,PT0,evo,i);
+t1 = dointpoint(lf,u,PNLX,evo,i);
+stdlinks(link,&lf->lfd,&lf->sp,i,fitv(des,i),res[5]);
+t1 = t1*t1*link[ZDDLL];
+t0 = t0*t0*link[ZDDLL];
+if (t1>1) t1 = 1;
+if (t0>1) t0 = 1; /* no observation gives >1 deg.free */
+res[0] += link[ZLIK];
+res[1] += t0;
+res[2] += t1;
+pw = prwt(&lf->lfd,i);
+if (pw>0)
+{ r1 += link[ZDLL]*link[ZDLL]/pw;
+r2 += link[ZDDLL]/pw;
+}
+}
+res[3] = 1.0;
+if ((fam(&lf->sp)&64)==64) /* quasi family */
+{ rdf = lf->lfd.n-2*res[1]+res[2];
+if (rdf<1.0)
+{ WARN(("Estimated rdf < 1.0; not estimating variance"));
+}
+else
+res[3] = r1/r2 * lf->lfd.n / rdf;
+}
+/* try to ensure consistency for family="circ"! */
+if (((fam(&lf->sp)&63)==TCIRC) & (lf->lfd.d==1))
+{ int *ind, nv;
+double dlt, th0, th1;
+ind = des->ind;
+nv = fp->nv;
+for (i=0; i<nv; i++) ind[i] = i;
+lforder(ind,evp(fp),0,nv-1);
+for (i=1; i<nv; i++)
+{ dlt = evptx(fp,ind[i],0)-evptx(fp,ind[i-1],0);
+th0 = fp->coef[ind[i]]-dlt*fp->coef[ind[i]+nv]-fp->coef[ind[i-1]];
+th1 = fp->coef[ind[i]]-dlt*fp->coef[ind[i-1]+nv]-fp->coef[ind[i-1]];
+if ((th0>PI)&(th1>PI))
+{ for (j=0; j<i; j++)
+fp->coef[ind[j]] += 2*PI;
+i--;
+}
+if ((th0<(-PI))&(th1<(-PI)))
+{ for (j=0; j<i; j++)
+fp->coef[ind[j]] -= 2*PI;
+i--;
+}
+}
+}
+}
+double rss(lf,des,df)
+lfit *lf;
+design *des;
+double *df;
+{ double ss, res[10];
+ss = 0;
+ressumm(lf,des,res);
+*df = lf->lfd.n - 2*res[1] + res[2];
+return(-2*res[0]);
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+/*
+*
+*   Derivative corrections. The local slopes are not the derivatives
+*   of the local likelihood estimate; the function dercor() computes
+*   the adjustment to get the correct derivatives under the assumption
+*   that h is constant.
+*
+*   By differentiating the local likelihood equations, one obtains
+*
+*     d ^      ^       T      -1   T  d    .       ^
+*    -- a   =  a  -  (X W V X)    X  -- W  l( Y, X a)
+*    dx  0      1                    dx
+*/
+#include "lfev.h"
+extern double robscale;
+void dercor(lfd,sp,des,coef)
+lfdata *lfd;
+smpar *sp;
+design *des;
+double *coef;
+{ double s1, dc[MXDIM], wd, link[LLEN];
+int i, ii, j, m, d, p;
+if (fam(sp)<=THAZ) return;
+if (ker(sp)==WPARM) return;
+d = lfd->d;
+p = des->p; m = des->n;
+if (lf_debug>1) mut_printf("  Correcting derivatives\n");
+fitfun(lfd, sp, des->xev,des->xev,des->f1,NULL);
+jacob_solve(&des->xtwx,des->f1);
+setzero(dc,d);
+/* correction term is e1^T (XTWVX)^{-1} XTW' ldot. */
+for (i=0; i<m; i++)
+{ s1 = innerprod(des->f1,d_xi(des,i),p);
+ii = des->ind[i];
+stdlinks(link,lfd,sp,ii,fitv(des,ii),robscale);
+for (j=0; j<d; j++)
+{ wd = wght(des,ii)*weightd(datum(lfd,j,ii)-des->xev[j],lfd->sca[j],
+d,ker(sp),kt(sp),des->h,lfd->sty[j],dist(des,ii));
+dc[j] += s1*wd*link[ZDLL];
+}
+}
+for (j=0; j<d; j++) coef[j+1] += dc[j];
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+void allocallcf(lf)
+lfit *lf;
+{ lf->fp.coef = VVEC(lf,0);
+lf->fp.h    = VVEC(lf,NPAR(lf));
+}
+int procvallcf(des,lf,v)
+design *des;
+lfit *lf;
+int v;
+{ int i, p, lf_status;
+lf_status = procv_nov(des,lf,v);
+if (lf_error) return(lf_status);
+p = NPAR(lf);
+for (i=0; i<p; i++) VVAL(lf,v,i) = des->cf[i];
+lf->fp.h[v] = des->h;
+return(0);
+}
+void initallcf(lf)
+lfit *lf;
+{ PROCV(lf) = procvallcf;
+ALLOC(lf) = allocallcf;
+PPROC(lf) = NULL;
+KEEPV(lf) = NPAR(lf)+1;
+KEEPC(lf) = 0;
+NOPC(lf)  = 1;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+void pprocgam(lf,des,res)
+lfit *lf;
+design *des;
+double *res;
+{ int i, j, n, evo, op;
+double *fv, *vr, df, t0, t1, u[MXDIM], link[LLEN];
+n = lf->lfd.n;
+fv = res;
+vr = &res[n];
+df = 0.0;
+evo = ev(&lf->evs);
+if (evo==EDATA) evo = EFITP;
+for (i=0; i<n; i++)
+{ for (j=0; j<lf->lfd.d; j++)  u[j] = datum(&lf->lfd,j,i);
+fitv(des,i) = base(&lf->lfd,i)+dointpoint(lf,u,PCOEF,evo,i);
+lf->lfd.y[i] -= fitv(des,i);
+wght(des,i) = 1.0;
+des->ind[i] = i;
+t0 = dointpoint(lf,u,PT0,evo,i);
+t1 = dointpoint(lf,u,PNLX,evo,i);
+stdlinks(link,&lf->lfd,&lf->sp,i,fitv(des,i),1.0);
+t0 = t0*t0*link[ZDDLL];
+t1 = t1*t1*link[ZDDLL];
+if (t0>1) t0 = 1; /* no observation gives >1 deg.free */
+if (t1>1) t1 = 1; /* no observation gives >1 deg.free */
+vr[i] = t1;
+df += t0;
+}
+des->n = n;
+deg(&lf->sp) = 1;
+op = npar(&lf->sp);
+npar(&lf->sp) = des->p = 1+lf->lfd.d;
+des->xev = lf->pc.xbar;
+locfit(&lf->lfd,des,&lf->sp,1,0,0);
+for (i=0; i<n; i++) fv[i] = resp(&lf->lfd,i) - fitv(des,i);
+for (i=0; i<=lf->lfd.d; i++)
+lf->pc.coef[i] += des->cf[i];
+res[2*n] = df-2;
+npar(&lf->sp) = op;
+}
+void initgam(lf)
+lfit *lf;
+{ initstd(lf);
+PPROC(lf) = pprocgam;
+KEEPC(lf) = 2*NOBS(lf)+1;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+int procvhatm(des,lf,v)
+design *des;
+lfit *lf;
+int v;
+{ int k;
+double *l;
+l = &lf->fp.coef[v*lf->lfd.n];
+if ((ker(&lf->sp)!=WPARM) | (!haspc(&lf->pc)))
+{ k = procv_nov(des,lf,v);
+wdiag(&lf->lfd,&lf->sp,des,l,&lf->dv,0,1,1);
+}
+else
+wdiagp(&lf->lfd,&lf->sp,des,l,&lf->pc,&lf->dv,0,1,1);
+lf->fp.h[v] = des->h;
+return(k);
+}
+void allochatm(lf)
+lfit *lf;
+{ lf->fp.coef = VVEC(lf,0);
+lf->fp.h    = VVEC(lf,NOBS(lf));
+}
+void pprochatm(lf,des,res)
+lfit *lf;
+design *des;
+double *res;
+{ transpose(lf->fp.coef,lf->fp.nvm,lf->lfd.n);
+}
+void inithatm(lf)
+lfit *lf;
+{ PROCV(lf) = procvhatm;
+ALLOC(lf) = allochatm;
+PPROC(lf) = pprochatm;
+KEEPV(lf) = NOBS(lf)+1;
+KEEPC(lf) = 1;
+NOPC(lf) = 1; /* could use pc if mi[MKER] == WPARM */
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+static lfit *lf_scb;
+static lfdata *lfd_scb;
+static smpar  *scb_sp;
+static design *des_scb;
+int scbfitter(x,l,reqd)
+double *x, *l;
+int reqd;
+{
+int m;
+des_scb->xev = x;
+if ((ker(scb_sp)!=WPARM) | (!haspc(&lf_scb->pc)))
+{ locfit(lfd_scb,des_scb,&lf_scb->sp,1,1,0);
+m = wdiag(lfd_scb, scb_sp, des_scb,l,&lf_scb->dv,reqd,2,0);
+}
+else
+m = wdiagp(lfd_scb, scb_sp, des_scb,l,&lf_scb->pc,&lf_scb->dv,reqd,2,0);
+return(m);
+}
+int constants(lf,des,res)
+lfit *lf;
+design *des;
+double *res;
+{
+int d, m, nt;
+double *wk;
+evstruc *evs;
+lf_scb = lf;
+des_scb = des;
+lfd_scb = &lf->lfd;
+scb_sp  = &lf->sp;
+evs = &lf->evs;
+d = lfd_scb->d;
+m = lfd_scb->n;
+trchck(lf,0,0,0);
+if (lf_error) return(0);
+if ((ker(scb_sp) != WPARM) && (lf->sp.nn>0))
+WARN(("constants are approximate for varying h"));
+npar(scb_sp) = calcp(scb_sp,lf->lfd.d);
+des_init(des,m,npar(scb_sp));
+set_scales(&lf->lfd);
+set_flim(&lf->lfd,&lf->evs);
+compparcomp(des,&lf->lfd,&lf->sp,&lf->pc,ker(scb_sp)!=WPARM);
+wk = &res[d+1];
+nt = tube_constants(scbfitter,d,m,ev(evs),mg(evs),evs->fl,
+res,wk,(d>3) ? 4 : d+1,0);
+lf->evs.nce = nt;   /* cheat way to return it to S. */
+return(nt);
+}
+void initkappa(lf)
+lfit *lf;
+{ PROCV(lf) = NULL;
+ALLOC(lf) = NULL;
+PPROC(lf) = (void *)constants;
+KEEPV(lf) = 0;
+KEEPC(lf) = NVAR(lf)+1+k0_reqd(NVAR(lf),NOBS(lf),0);
+NOPC(lf) = 0;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+/*  fix df computation for link=IDENT. */
+void pproclscv(lf,des,res)
+lfit *lf;
+design *des;
+double *res;
+{ double df, fh, fh_cv, infl, z0, z1, x[MXDIM];
+int i, n, j, evo;
+z1 = df = 0.0;
+evo = ev(&lf->evs);
+n = lf->lfd.n;
+if ((evo==EDATA) | (evo==ECROS)) evo = EFITP;
+z0 = dens_integrate(lf,des,2);
+for (i=0; i<n; i++)
+{ for (j=0; j<lf->lfd.d; j++) x[j] = datum(&lf->lfd,j,i);
+fh = base(&lf->lfd,i)+dointpoint(lf,x,PCOEF,evo,i);
+if (link(&lf->sp)==LLOG) fh = exp(fh);
+infl = dointpoint(lf,x,PT0,evo,i);
+infl = infl * infl;
+if (infl>1) infl = 1;
+fh_cv = (link(&lf->sp) == LIDENT) ?
+(n*fh - infl) / (n-1.0) : fh*(1-infl)*n/(n-1.0);
+z1 += fh_cv;
+df += infl;
+}
+res[0] = z0-2*z1/n;
+res[1] = df;
+}
+void initlscv(lf)
+lfit *lf;
+{ initstd(lf);
+KEEPC(lf) = 2;
+PPROC(lf) = pproclscv;
+NOPC(lf) = 1;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+static double pen, sig2;
+void goldensec(f,des,tr,eps,xm,ym,meth)
+double (*f)(), eps, *xm, *ym;
+int meth;
+design *des;
+lfit *tr;
+{ double x[4], y[4], xx[11], yy[11];
+int i, im;
+xx[0] = tr->sp.fixh;
+if (xx[0]<=0)
+{ LERR(("regband: initialize h>0"));
+return;
+}
+for (i=0; i<=10; i++)
+{ if (i>0) xx[i] = (1+gold_rat)*xx[i-1];
+yy[i] = f(xx[i],des,tr,meth);
+if ((i==0) || (yy[i]<yy[im])) im = i;
+}
+if (im==0) im = 1;
+if (im==10)im = 9;
+x[0] = xx[im-1]; y[0] = yy[im-1];
+x[1] = xx[im];   y[1] = yy[im];
+x[3] = xx[im+1]; y[3] = yy[im+1];
+x[2] = gold_rat*x[3]+(1-gold_rat)*x[0];
+y[2] = f(x[2],des,tr,meth);
+while (x[3]-x[0]>eps)
+{ if (y[1]<y[2])
+{ x[3] = x[2]; y[3] = y[2];
+x[2] = x[1]; y[2] = y[1];
+x[1] = gold_rat*x[0]+(1-gold_rat)*x[3];
+y[1] = f(x[1],des,tr,meth);
+}
+else
+{ x[0] = x[1]; y[0] = y[1];
+x[1] = x[2]; y[1] = y[2];
+x[2] = gold_rat*x[3]+(1-gold_rat)*x[0];
+y[2] = f(x[2],des,tr,meth);
+}
+}
+im = 0;
+for (i=1; i<4; i++) if (y[i]<y[im]) im = i;
+*xm = x[im]; *ym = y[im];
+}
+double dnk(x,k)
+double x;
+int k;
+{ double f;
+switch(k)
+{ case 0: f = 1; break;
+case 1: f = -x; break;
+case 2: f = x*x-1; break;
+case 3: f = x*(x*x-3); break;
+case 4: f = 3-x*x*(6-x*x); break;
+case 5: f = -x*(15-x*x*(10-x*x)); break;
+case 6: f = -15+x*x*(45-x*x*(15-x*x)); break;
+default: LERR(("dnk: k=%d too large",k)); return(0.0);
+}
+return(f*exp(-x*x/2)/S2PI);
+}
+double locai(h,des,lf)
+double h;
+design *des;
+lfit *lf;
+{ double cp, res[10];
+nn(&lf->sp) = h;
+lf->mdl.procv = procvstd;
+nstartlf(des,lf);
+ressumm(lf,des,res);
+cp = -2*res[0] + pen*res[1];
+return(cp);
+}
+static int fc;
+double loccp(h,des,lf,m) /* m=1: cp    m=2: gcv */
+double h;
+design *des;
+lfit *lf;
+int m;
+{ double cp, res[10];
+int dg, n;
+n = lf->lfd.n;
+nn(&lf->sp) = 0;
+fixh(&lf->sp) = h;
+lf->mdl.procv = procvstd;
+nstartlf(des,lf);
+ressumm(lf,des,res);
+if (m==1)
+{ if (fc) sig2 = res[3]; /* first call - set sig2 */
+cp = -2*res[0]/sig2 - n + 2*res[1];
+}
+else
+cp = -2*n*res[0]/((n-res[1])*(n-res[1]));
+fc = 0;
+return(cp);
+}
+double cp(des,lf,meth)
+design *des;
+lfit *lf;
+int meth;
+{ double hm, ym;
+fc = 1;
+goldensec(loccp,des,lf,0.001,&hm,&ym,meth);
+return(hm);
+}
+double gkk(des,lf)
+design *des;
+lfit *lf;
+{ double h, h5, nf, th;
+int i, j, n, dg0, dg1;
+ev(&lf->evs)= EDATA;
+nn(&lf->sp) = 0;
+n = lf->lfd.n;
+dg0 = deg0(&lf->sp);     /* target degree */
+dg1 = dg0+1+(dg0%2==0);  /* pilot degree */
+nf = exp(log(1.0*n)/10); /* bandwidth inflation factor */
+h = lf->sp.fixh;         /* start bandwidth */
+for (i=0; i<=10; i++)
+{ deg(&lf->sp) = dg1;
+lf->sp.fixh = h*nf;
+lf->mdl.procv = procvstd;
+nstartlf(des,lf);
+th = 0;
+for (j=10; j<n-10; j++)
+th += lf->fp.coef[dg1*n+j]*lf->fp.coef[dg1*n+j];
+th *= n/(n-20.0);
+h5 = sig2 * Wikk(ker(&lf->sp),dg0) / th;
+h = exp(log(h5)/(2*dg1+1));
+if (lf_error) return(0.0);
+/* mut_printf("pilot %8.5f  sel %8.5f\n",lf->sp.fixh,h); */
+}
+return(h);
+}
+double rsw(des,lf)
+design *des;
+lfit *lf;
+{ int i, j, k, nmax, nvm, n, mk, evo, dg0, dg1;
+double rss[6], cp[6], th22, dx, d2, hh;
+nmax = 5;
+evo = ev(&lf->evs); ev(&lf->evs) = EGRID;
+mk = ker(&lf->sp);  ker(&lf->sp) = WRECT;
+dg0 = deg0(&lf->sp);
+dg1 = 1 + dg0 + (dg0%2==0);
+deg(&lf->sp) = 4;
+for (k=nmax; k>0; k--)
+{ lf->evs.mg[0] = k;
+lf->evs.fl[0] = 1.0/(2*k);
+lf->evs.fl[1] = 1-1.0/(2*k);
+nn(&lf->sp) = 0;
+fixh(&lf->sp) = 1.0/(2*k);
+lf->mdl.procv = procvstd;
+nstartlf(des,lf);
+nvm = lf->fp.nvm;
+rss[k] = 0;
+for (i=0; i<k; i++) rss[k] += -2*lf->fp.lik[i];
+}
+n = lf->lfd.n; k = 1;
+for (i=1; i<=nmax; i++)
+{ /* cp[i] = (n-5*nmax)*rss[i]/rss[nmax]-(n-10*i); */
+cp[i] = rss[i]/sig2-(n-10*i);
+if (cp[i]<cp[k]) k = i;
+}
+lf->evs.mg[0] = k;
+lf->evs.fl[0] = 1.0/(2*k);
+lf->evs.fl[1] = 1-1.0/(2*k);
+nn(&lf->sp) = 0;
+fixh(&lf->sp) = 1.0/(2*k);
+lf->mdl.procv = procvstd;
+nstartlf(des,lf);
+ker(&lf->sp) = mk; ev(&lf->evs) = evo;
+nvm = lf->fp.nvm;
+th22 = 0;
+for (i=10; i<n-10; i++)
+{ j = floor(k*datum(&lf->lfd,0,i));
+if (j>=k) j = k-1;
+dx = datum(&lf->lfd,0,i)-evptx(&lf->fp,0,j);
+if (dg1==2)
+d2 = lf->fp.coef[2*nvm+j]+dx*lf->fp.coef[3*nvm+j]+dx*dx*lf->fp.coef[4*nvm+j]/2;
+else d2 = lf->fp.coef[4*nvm+j];
+th22 += d2*d2;
+}
+hh = Wikk(mk,dg0)*sig2/th22*(n-20.0)/n;
+return(exp(log(hh)/(2*dg1+1)));
+}
+#define MAXMETH 10
+void rband(lf,des,hhat)
+lfit *lf;
+design *des;
+double *hhat;
+{ int i, dg, nmeth, meth[MAXMETH];
+double h0, res[10];
+nmeth = lf->dv.nd;
+for (i=0; i<nmeth; i++) meth[i] = lf->dv.deriv[i];
+lf->dv.nd = 0;
+/* first, estimate sigma^2.
+* this is ridiculously bad. lf->sp.fixh = 0.05????
+*/
+/*  dg = deg(&lf->sp); deg(&lf->sp) = 2;
+h0 = lf->sp.fixh;  lf->sp.fixh = 0.05;
+lf->mdl.procv = procvstd;
+nstartlf(des,lf);
+ressumm(lf,des,res);
+deg(&lf->sp) = dg; lf->sp.fixh = h0;
+sig2 = res[3];  */
+for (i=0; i<nmeth; i++)
+{
+switch(meth[i])
+{ case 0: hhat[i] = cp(des,lf,1);
+break;
+case 1: hhat[i] = cp(des,lf,2);
+break;
+case 2: hhat[i] = gkk(des,lf);
+break;
+case 3: hhat[i] = rsw(des,lf);
+break;
+default: hhat[i] = 0;
+mut_printf("Unknown method %d\n",meth[i]);
+}
+if (lf_error) i = nmeth;
+}
+lf->dv.nd = nmeth;
+}
+void initrband(lf)
+lfit *lf;
+{
+initstd(lf);
+KEEPC(lf) = MAXMETH;
+PROCV(lf) = NULL;
+PPROC(lf) = rband;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+static double scb_crit, *x, c[10], kap[5], kaq[5], max_p2;
+static int side, type;
+design *scb_des;
+double covar_par(lf,des,x1,x2)
+lfit *lf;
+design *des;
+double x1, x2;
+{ double *v1, *v2, *wk;
+paramcomp *pc;
+int i, j, p, ispar;
+v1 = des->f1; v2 = des->ss; wk = des->oc;
+ispar = (ker(&lf->sp)==WPARM) && (haspc(&lf->pc));
+p = npar(&lf->sp);
+/*  for parametric models, the covariance is
+*  A(x1)^T (X^T W V X)^{-1} A(x2)
+*  which we can find easily from the parametric component.
+*/
+if (ispar)
+{ pc = &lf->pc;
+fitfun(&lf->lfd, &lf->sp, &x1,pc->xbar,v1,NULL);
+fitfun(&lf->lfd, &lf->sp, &x2,pc->xbar,v2,NULL);
+jacob_hsolve(&lf->pc.xtwx,v1);
+jacob_hsolve(&lf->pc.xtwx,v2);
+}
+/*  for non-parametric models, we must use the cholseky decomposition
+*  of M2 = X^T W^2 V X. Courtesy of lf_vcov caclulations, should have
+*  des->P = M2^{1/2} M1^{-1}.
+*/
+if (!ispar)
+{ fitfun(&lf->lfd, &lf->sp, &x1,des->xev,wk,NULL);
+for (i=0; i<p; i++)
+{ v1[i] = 0;
+for (j=0; j<p; j++) v1[i] += des->P[i*p+j]*wk[j];
+}
+fitfun(&lf->lfd, &lf->sp, &x2,des->xev,wk,NULL);
+for (i=0; i<p; i++)
+{ v2[i] = 0;
+for (j=0; j<p; j++) v2[i] += des->P[i*p+j]*wk[j];
+}
+}
+return(innerprod(v1,v2,p));
+}
+void cumulant(lf,des,sd)
+lfit *lf;
+design *des;
+double sd;
+{ double b2i, b3i, b3j, b4i;
+double ss, si, sj, uii, uij, ujj, k1;
+int ii, i, j, jj;
+for (i=1; i<10; i++) c[i] = 0.0;
+k1 = 0;
+/* ss = sd*sd; */
+ss = covar_par(lf,des,des->xev[0],des->xev[0]);
+/*
+* this isn't valid for nonparametric models. At a minimum,
+* the sums would have to include weights. Still have to work
+* out the right way.
+*/
+for (i=0; i<lf->lfd.n; i++)
+{ ii = des->ind[i];
+b2i = b2(fitv(des,ii),fam(&lf->sp),prwt(&lf->lfd,ii));
+b3i = b3(fitv(des,ii),fam(&lf->sp),prwt(&lf->lfd,ii));
+b4i = b4(fitv(des,ii),fam(&lf->sp),prwt(&lf->lfd,ii));
+si = covar_par(lf,des,des->xev[0],datum(&lf->lfd,0,ii));
+uii= covar_par(lf,des,datum(&lf->lfd,0,ii),datum(&lf->lfd,0,ii));
+if (lf_error) return;
+c[2] += b4i*si*si*uii;
+c[6] += b4i*si*si*si*si;
+c[7] += b3i*si*uii;
+c[8] += b3i*si*si*si;
+/* c[9] += b2i*si*si*si*si;
+c[9] += b2i*b2i*si*si*si*si; */
+k1 += b3i*si*(si*si/ss-uii);
+/* i=j components */
+c[1] += b3i*b3i*si*si*uii*uii;
+c[3] += b3i*b3i*si*si*si*si*uii;
+c[4] += b3i*b3i*si*si*uii*uii;
+for (j=i+1; j<lf->lfd.n; j++)
+{ jj = des->ind[j];
+b3j = b3(fitv(des,ii),fam(&lf->sp),prwt(&lf->lfd,jj));
+sj = covar_par(lf,des,des->xev[0],datum(&lf->lfd,0,jj));
+uij= covar_par(lf,des,datum(&lf->lfd,0,ii),datum(&lf->lfd,0,jj));
+ujj= covar_par(lf,des,datum(&lf->lfd,0,jj),datum(&lf->lfd,0,jj));
+c[1] += 2*b3i*b3j*si*sj*uij*uij;
+c[3] += 2*b3i*b3j*si*si*sj*sj*uij;
+c[4] += b3i*b3j*uij*(si*si*ujj+sj*sj*uii);
+if (lf_error) return;
+}
+}
+c[5] = c[1];
+c[7] = c[7]*c[8];
+c[8] = c[8]*c[8];
+c[1] /= ss; c[2] /= ss; c[3] /= ss*ss; c[4] /= ss;
+c[5] /= ss; c[6] /= ss*ss; c[7] /= ss*ss;
+c[8] /= ss*ss*ss; c[9] /= ss*ss;
+/* constants used in p(x,z) computation */
+kap[1] = k1/(2*sqrt(ss));
+kap[2] = 1 + 0.5*(c[1]-c[2]+c[4]-c[7]) - 3*c[3] + c[6] + 1.75*c[8];
+kap[4] = -9*c[3] + 3*c[6] + 6*c[8] + 3*c[9];
+/* constants used in q(x,u) computation */
+kaq[2] = c[3] - 1.5*c[8] - c[5] - c[4] + 0.5*c[7] + c[6] - c[2];
+kaq[4] = -3*c[3] - 6*c[4] - 6*c[5] + 3*c[6] + 3*c[7] - 3*c[8] + 3*c[9];
+}
+/* q2(u) := u+q2(x,u) in paper */
+double q2(u)
+double u;
+{ return(u-u*(36.0*kaq[2] + 3*kaq[4]*(u*u-3) + c[8]*((u*u-10)*u*u+15))/72.0);
+}
+/*  p2(u) := p2(x,u) in paper */
+double p2(u)
+double u;
+{ return( -u*( 36*(kap[2]-1+kap[1]*kap[1])
++ 3*(kap[4]+4*kap[1]*sqrt(kap[3]))*(u*u-3)
++ c[8]*((u*u-10)*u*u+15) ) / 72 );
+}
+extern int likereg();
+double gldn_like(a)
+double a;
+{ int err;
+scb_des->fix[0] = 1;
+scb_des->cf[0] = a;
+max_nr(likereg, scb_des->cf, scb_des->oc, scb_des->res, scb_des->f1,
+&scb_des->xtwx, scb_des->p, lf_maxit, 1.0e-6, &err);
+scb_des->fix[0] = 0;
+return(scb_des->llk);
+}
+/* v1/v2 is correct for deg=0 only */
+void get_gldn(fp,des,lo,hi,v)
+fitpt *fp;
+design *des;
+double *lo, *hi;
+int v;
+{ double v1, v2, c, tlk;
+int err;
+v1 = fp->nlx[v];
+v2 = fp->t0[v];
+c = scb_crit * v1 / v2;
+tlk = des->llk - c*c/2;
+mut_printf("v %8.5f %8.5f  c %8.5f  tlk %8.5f   llk %8.5f\n",v1,v2,c,tlk,des->llk);
+/* want: { a : l(a) >= l(a-hat) - c*c/2 } */
+lo[v] = fp->coef[v] - scb_crit*v1;
+hi[v] = fp->coef[v] + scb_crit*v1;
+err = 0;
+mut_printf("lo %2d\n",v);
+lo[v] = solve_secant(gldn_like,tlk,lo[v],fp->coef[v],1e-8,BDF_EXPLEFT,&err);
+if (err>0) mut_printf("solve_secant error\n");
+mut_printf("hi %2d\n",v);
+hi[v] = solve_secant(gldn_like,tlk,fp->coef[v],hi[v],1e-8,BDF_EXPRIGHT,&err);
+if (err>0) mut_printf("solve_secant error\n");
+}
+int procvscb2(des,lf,v)
+design *des;
+lfit *lf;
+int v;
+{ double thhat, sd, *lo, *hi, u;
+int err, st, tmp;
+x = des->xev = evpt(&lf->fp,v);
+tmp = haspc(&lf->pc);
+/* if ((ker(&lf->sp)==WPARM) && (haspc(&lf->pc)))
+{ lf->coef[v] = thhat = addparcomp(lf,des->xev,PCOEF);
+lf->nlx[v] = lf->t0[v] = sd = addparcomp(lf,des->xev,PNLX);
+}
+else */
+{ haspc(&lf->pc) = 0;
+st = procvstd(des,lf,v);
+thhat = lf->fp.coef[v];
+sd = lf->fp.nlx[v];
+}
+if ((type==GLM2) | (type==GLM3) | (type==GLM4))
+{ if (ker(&lf->sp) != WPARM)
+WARN(("nonparametric fit; correction is invalid"));
+cumulant(lf,des,sd);
+}
+haspc(&lf->pc) = tmp;
+lo = lf->fp.t0;
+hi = &lo[lf->fp.nvm];
+switch(type)
+{
+case GLM1:
+return(st);
+case GLM2: /* centered scr */
+lo[v] = kap[1];
+hi[v] = sqrt(kap[2]);
+return(st);
+case GLM3: /* corrected 2 */
+lo[v] = solve_secant(q2,scb_crit,0.0,2*scb_crit,0.000001,BDF_NONE,&err);
+return(st);
+case GLM4: /* corrected 2' */
+u = fabs(p2(scb_crit));
+max_p2 = MAX(max_p2,u);
+return(st);
+case GLDN:
+get_gldn(&lf->fp,des,lo,hi,v);
+return(st);
+}
+LERR(("procvscb2: invalid type"));
+return(st);
+}
+void scb(lf,des,res)
+lfit *lf;
+design *des;
+double *res;
+{ double k1, k2, *lo, *hi, sig, thhat, nlx, rss[10];
+int i, nterms;
+scb_des= des;
+npar(&lf->sp) = calcp(&lf->sp,lf->lfd.d);
+des_init(des,lf->lfd.n,npar(&lf->sp));
+type = geth(&lf->fp);
+if (type >= 80) /* simultaneous */
+{
+nterms = constants(lf,des,res);
+scb_crit = critval(0.05,res,nterms,lf->lfd.d,TWO_SIDED,0.0,GAUSS);
+type -= 10;
+}
+else /* pointwise */
+{ res[0] = 1;
+scb_crit = critval(0.05,res,1,lf->lfd.d,TWO_SIDED,0.0,GAUSS);
+}
+max_p2 = 0.0;
+lf->mdl.procv = procvscb2;
+nstartlf(des,lf);
+if ((fam(&lf->sp)&64)==64)
+{ i = haspc(&lf->pc); haspc(&lf->pc) = 0;
+ressumm(lf,des,rss);
+haspc(&lf->pc) = i;
+sig = sqrt(rss[3]);
+}
+else sig = 1.0;
+lo = lf->fp.t0;
+hi = &lo[lf->fp.nvm];
+for (i=0; i<lf->fp.nv; i++)
+{ thhat = lf->fp.coef[i];
+nlx = lf->fp.nlx[i];
+switch(type)
+{
+case GLM1:  /* basic scb */
+lo[i] = thhat - scb_crit * sig * nlx;
+hi[i] = thhat + scb_crit * sig * nlx;
+break;
+case GLM2:
+k1 = lo[i];
+k2 = hi[i];
+lo[i] = thhat - k1*nlx - scb_crit*nlx*k2;
+hi[i] = thhat - k1*nlx + scb_crit*nlx*k2;
+break;
+case GLM3:
+k1 = lo[i];
+lo[i] = thhat - k1*nlx;
+hi[i] = thhat + k1*nlx;
+case GLM4:  /* corrected 2' */
+lo[i] = thhat - (scb_crit-max_p2)*lf->fp.nlx[i];
+hi[i] = thhat + (scb_crit-max_p2)*lf->fp.nlx[i];
+break;
+case GLDN:
+break;
+}
+}
+}
+void initscb(lf)
+lfit *lf;
+{ initstd(lf);
+PROCV(lf) = NULL;
+KEEPC(lf) = NVAR(lf)+1+k0_reqd(NVAR(lf),NOBS(lf),0);
+PPROC(lf) = scb;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+int procvsimple(des,lf,v)
+design *des;
+lfit *lf;
+int v;
+{ int lf_status;
+lf_status = procv_nov(des,lf,v);
+VVAL(lf,v,0) = des->cf[cfn(des,0)];
+return(lf_status);
+}
+void allocsimple(lf)
+lfit *lf;
+{ lf->fp.coef = VVEC(lf,0);
+lf->fp.h = NULL;
+}
+void initsimple(lf)
+lfit *lf;
+{
+PROCV(lf) = procvsimple;
+ALLOC(lf) = allocsimple;
+PPROC(lf) = NULL;
+KEEPV(lf) = 1;
+KEEPC(lf) = 1;
+NOPC(lf)  = 1;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+/* 3d+8 variables to keep:
+* d+1 coef+derivs.
+* d+1 sd's + derivs.
+* d+1 infl + derivs.
+*   3 likelihood and d.f's.
+*   1 bandwidth h
+*   1 degree.
+*/
+void allocstd(lf)
+lfit *lf;
+{ int d;
+d = NVAR(lf);
+lf->fp.coef = VVEC(lf,0);
+lf->fp.nlx  = VVEC(lf,d+1);
+lf->fp.t0   = VVEC(lf,2*d+2);
+lf->fp.lik  = VVEC(lf,3*d+3);
+lf->fp.h    = VVEC(lf,3*d+6);
+lf->fp.deg  = VVEC(lf,3*d+7);
+}
+int procvstd(des,lf,v)
+design *des;
+lfit *lf;
+int v;
+{ int d, p, nvm, i, k;
+double t0[1+MXDIM], vari[1+MXDIM];
+k = procv_var(des,lf,v);
+if (lf_error) return(k);
+d = lf->lfd.d;
+p = npar(&lf->sp);
+nvm = lf->fp.nvm;
+if (k != LF_OK) lf_status_msg(k);
+lf->fp.lik[v] = des->llk;
+lf->fp.lik[nvm+v] = des->tr2;
+lf->fp.lik[2*nvm+v] = des->tr0 - des->tr2;
+for (i=0; i<des->ncoef; i++)
+vari[i] = des->V[p*cfn(des,0) + cfn(des,i)];
+vari[0] = sqrt(vari[0]);
+if (vari[0]>0) for (i=1; i<des->ncoef; i++) vari[i] /= vari[0];
+t0[0] = sqrt(des->f1[0]);
+if (t0[0]>0) for (i=1; i<des->ncoef; i++) t0[i] = des->f1[i]/t0[0];
+if (dc(&lf->fp)) dercor(&lf->lfd,&lf->sp,des,des->cf);
+subparcomp(des,lf,des->cf);
+for (i=0; i<des->ncoef; i++)
+lf->fp.coef[i*lf->fp.nvm+v] =  des->cf[cfn(des,i)];
+subparcomp2(des,lf,vari,t0);
+for (i=0; i<des->ncoef; i++)
+{ lf->fp.nlx[i*nvm+v] = vari[i];
+lf->fp.t0[i*nvm+v]  = t0[i];
+}
+lf->fp.deg[v] = deg(&lf->sp);
+return(k);
+}
+void pprocstd(lf,des,res)
+lfit *lf;
+design *des;
+double *res;
+{
+ressumm(lf,des,res);
+}
+void initstd(lf)
+lfit *lf;
+{ PROCV(lf) = procvstd;
+ALLOC(lf) = allocstd;
+PPROC(lf) = pprocstd;
+KEEPV(lf) = 3*NVAR(lf) + 8;
+KEEPC(lf) = 6;
+NOPC(lf)  = 0;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+extern void procstd(), allocstd();
+extern double robscale;
+double vocri(lk,t0,t2,pen)
+double lk, t0, t2, pen;
+{ if (pen==0) return(-2*t0*lk/((t0-t2)*(t0-t2)));
+return((-2*lk+pen*t2)/t0);
+}
+double intvo(des,lf,c0,c1,a,p,t0,t20,t21)
+design *des;
+lfit *lf;
+double *c0, *c1, a, t0, t20, t21;
+int p;
+{ double th, lk, link[LLEN];
+int i, ii;
+lk = 0;
+for (i=0; i<des->n; i++)
+{ ii = des->ind[i];
+th = (1-a)*innerprod(c0,d_xi(des,i),p) + a*innerprod(c1,d_xi(des,i),p);
+stdlinks(link,&lf->lfd,&lf->sp,ii,th,robscale);
+lk += wght(des,ii)*link[ZLIK];
+}
+des->llk = lk;
+return(vocri(des->llk,t0,(1-a)*t20+a*t21,pen(&lf->sp)));
+}
+int procvvord(des,lf,v)
+design *des;
+lfit *lf;
+int v;
+{ double tr[6], gcv, g0, ap, coef[4][10], t2[4], th, md;
+int i, j, k, d1, i0, p1, ip;
+des->xev = evpt(&lf->fp,v);
+ap = pen(&lf->sp);
+if ((ap==0) & ((fam(&lf->sp)&63)!=TGAUS)) ap = 2.0;
+d1 = deg(&lf->sp); p1 = npar(&lf->sp);
+for (i=0; i<p1; i++) coef[0][i] = coef[1][i] = coef[2][i] = coef[3][i] = 0.0;
+i0 = 0; g0 = 0;
+ip = 1;
+for (i=deg0(&lf->sp); i<=d1; i++)
+{ deg(&lf->sp) = i;
+des->p = npar(&lf->sp) = calcp(&lf->sp,lf->lfd.d);
+k = locfit(&lf->lfd,des,&lf->sp,0, i==deg0(&lf->sp),0);
+local_df(&lf->lfd,&lf->sp,des,tr);
+gcv = vocri(des->llk,tr[0],tr[2],ap);
+if ((i==deg0(&lf->sp)) || (gcv<g0)) { i0 = i; g0 = gcv; md = i; }
+for (j=0; j<des->p; j++) coef[i][j] = des->cf[j];
+t2[i] = tr[2];
+#ifdef RESEARCH
+if ((ip) && (i>deg0(&lf->sp)))
+{ for (j=1; j<10; j++)
+{ gcv = intvo(des,lf,coef[i-1],coef[i],j/10.0,des->p,tr[0],t2[i-1],t2[i]);
+if (gcv<g0) { g0 = gcv; md = i-1+j/10.0; }
+}
+}
+#endif
+}
+lf->fp.h[v] = des->h;
+if (lf->fp.h[v]<=0) WARN(("zero bandwidth in procvvord"));
+if (i0<d1) /* recompute the best fit */
+{ deg(&lf->sp) = i0;
+des->p = npar(&lf->sp) = calcp(&lf->sp,lf->lfd.d);
+k = locfit(&lf->lfd,des,&lf->sp,0,0,0);
+for (i=npar(&lf->sp); i<p1; i++) des->cf[i] = 0.0;
+i0 = md; if (i0==d1) i0--;
+th = md-i0;
+for (i=0; i<p1; i++) des->cf[i] = (1-th)*coef[i0][i]+th*coef[i0+1][i];
+deg(&lf->sp) = d1; npar(&lf->sp) = p1;
+}
+for (i=0; i<p1; i++) lf->fp.coef[i*lf->fp.nvm+v] = des->cf[i];
+lf->fp.deg[v] = md;
+return(k);
+}
+void initvord(lf)
+lfit *lf;
+{ initstd(lf);
+PROCV(lf) = procvvord;
+ALLOC(lf) = allocstd;
+PPROC(lf) = NULL;
+KEEPC(lf) = 0;
+NOPC(lf)  = 1;
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+/*
+* functions for computing and subtracting, adding the
+* parametric component
+*/
+#include "lfev.h"
+int noparcomp(sp)
+smpar *sp;
+{ int tg;
+if (ubas(sp)) return(1);
+tg = fam(sp) & 63;
+if (tg<=THAZ) return(1);
+if (tg==TROBT) return(1);
+if (tg==TCAUC) return(1);
+if (tg==TQUANT) return(1);
+return(0);
+}
+int pc_reqd(d,p)
+int d, p;
+{ return(d + 2*p + jac_reqd(p));
+}
+void pcchk(pc,d,p,lc)
+paramcomp *pc;
+int d, p, lc;
+{ int rw;
+double *z;
+rw = pc_reqd(d,p);
+if (pc->lwk < rw)
+{ pc->wk = (double *)calloc(rw,sizeof(double));
+if ( pc->wk == NULL ) {
+printf("Problem allocating memory for pc->wk\n");fflush(stdout);
+}
+pc->lwk= rw;
+}
+z = pc->wk;
+pc->xbar = z; z += d;
+pc->coef = z; z += p;
+pc->f    = z; z += p;
+z = jac_alloc(&pc->xtwx,p,z);
+pc->xtwx.p = p;
+}
+void compparcomp(des,lfd,sp,pc,nopc)
+design *des;
+lfdata *lfd;
+smpar *sp;
+paramcomp *pc;
+int nopc;
+{ int i, j, k, p;
+double wt, sw;
+if (lf_debug>1) mut_printf(" compparcomp:\n");
+p = des->p;
+pcchk(pc,lfd->d,p,1);
+for (i=0; i<lfd->d; i++) pc->xbar[i] = 0.0;
+sw = 0.0;
+for (i=0; i<lfd->n; i++)
+{
+wt = prwt(lfd,i);
+sw += wt;
+for (j=0; j<lfd->d; j++)
+pc->xbar[j] += datum(lfd,j,i)*wt;
+des->ind[i] = i;
+wght(des,i) = 1.0;
+}
+for (i=0; i<lfd->d; i++) pc->xbar[i] /= sw;
+if ((nopc) || noparcomp(sp))
+{ haspc(pc) = 0;
+return;
+}
+haspc(pc) = 1;
+des->xev = pc->xbar;
+k = locfit(lfd,des,sp,0,0,0);
+if (k != LF_OK) lf_status_msg(k);
+if (lf_error) return;
+switch(k)
+{ case LF_NOPT: return;
+case LF_INFA: return;
+case LF_NCON: return;
+case LF_OOB: return;
+case LF_NSLN: return;
+case LF_PF:
+WARN(("compparcomp: perfect fit"));
+case LF_OK:
+case LF_DONE:
+for (i=0; i<p; i++)
+{ pc->coef[i] = des->cf[i];
+pc->xtwx.dg[i] = des->xtwx.dg[i];
+pc->xtwx.wk[i] = des->xtwx.wk[i];
+}
+for (i=0; i<p*p; i++)
+{ pc->xtwx.Z[i] = des->xtwx.Z[i];
+pc->xtwx.Q[i] = des->xtwx.Q[i];
+}
+pc->xtwx.sm = des->xtwx.sm;
+pc->xtwx.st = des->xtwx.st;
+return;
+default:
+LERR(("compparcomp: locfit unknown return status %d",k));
+return;
+}
+}
+void subparcomp(des,lf,coef)
+design *des;
+lfit *lf;
+double *coef;
+{ int i, nd;
+deriv *dv;
+paramcomp *pc;
+pc = &lf->pc;
+if (!haspc(pc)) return;
+dv = &lf->dv; nd = dv->nd;
+fitfun(&lf->lfd, &lf->sp, des->xev,pc->xbar,des->f1,dv);
+coef[0] -= innerprod(pc->coef,des->f1,pc->xtwx.p);
+if (des->ncoef == 1) return;
+dv->nd = nd+1;
+for (i=0; i<lf->lfd.d; i++)
+{ dv->deriv[nd] = i;
+fitfun(&lf->lfd, &lf->sp, des->xev,pc->xbar,des->f1,dv);
+coef[i+1] -= innerprod(pc->coef,des->f1,pc->xtwx.p);
+}
+dv->nd = nd;
+}
+void subparcomp2(des,lf,vr,il)
+design *des;
+lfit *lf;
+double *vr, *il;
+{ double t0, t1;
+int i, nd;
+deriv *dv;
+paramcomp *pc;
+pc = &lf->pc;
+if (!haspc(pc)) return;
+dv = &lf->dv; nd = dv->nd;
+fitfun(&lf->lfd, &lf->sp, des->xev,pc->xbar,des->f1,dv);
+for (i=0; i<npar(&lf->sp); i++) pc->f[i] = des->f1[i];
+jacob_solve(&pc->xtwx,des->f1);
+t0 = sqrt(innerprod(pc->f,des->f1,pc->xtwx.p));
+vr[0] -= t0;
+il[0] -= t0;
+if ((t0==0) | (des->ncoef==1)) return;
+dv->nd = nd+1;
+for (i=0; i<lf->lfd.d; i++)
+{ dv->deriv[nd] = i;
+fitfun(&lf->lfd, &lf->sp, des->xev,pc->xbar,pc->f,dv);
+t1 = innerprod(pc->f,des->f1,pc->xtwx.p)/t0;
+vr[i+1] -= t1;
+il[i+1] -= t1;
+}
+dv->nd = nd;
+}
+double addparcomp(lf,x,c)
+lfit *lf;
+double *x;
+int c;
+{ double y;
+paramcomp *pc;
+pc = &lf->pc;
+if (!haspc(pc)) return(0.0);
+fitfun(&lf->lfd, &lf->sp, x,pc->xbar,pc->f,&lf->dv);
+if (c==PCOEF) return(innerprod(pc->coef,pc->f,pc->xtwx.p));
+if ((c==PNLX)|(c==PT0)|(c==PVARI))
+{ y = sqrt(jacob_qf(&pc->xtwx,pc->f));
+return(y);
+}
+return(0.0);
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+/*
+preplot():  interpolates the fit to a new set of points.
+lf  -- the fit structure.
+x   -- the points to predict at.
+f   -- vector to return the predictions.
+se  -- vector to return std errors (NULL if not req'd)
+band-- char for conf band type. ('n'=none, 'g'=global etc.)
+n   -- no of predictions (or vector of margin lengths for grid)
+where -- where to predict:
+1 = points in the array x.
+2 = grid defined by margins in x.
+3 = data points from lf (ignore x).
+4 = fit points from lf (ignore x).
+what -- what to predict.
+(PCOEF etc; see lfcons.h file)
+*/
+#define NWH 8
+static char *whtyp[NWH] = { "coef", "nlx", "infl", "band",
+"degr", "like", "rdf", "vari" };
+static int   whval[NWH] = { PCOEF, PNLX, PT0, PBAND, PDEGR, PLIK, PRDF, PVARI };
+int ppwhat(z)
+char *z;
+{ int val;
+val = pmatch(z, whtyp, whval, NWH, -1);
+if (val==-1) LERR(("Unknown what = %s",z));
+return(val);
+}
+static char cb;
+double *sef, *fit, sigmahat;
+void predptall(lf,x,what,ev,i)
+lfit *lf;
+double *x;
+int what, ev, i;
+{ double lik, rdf;
+fit[i] = dointpoint(lf,x,what,ev,i);
+if (cb=='n') return;
+sef[i] = dointpoint(lf,x,PNLX,ev,i);
+if (cb=='g')
+{ sef[i] *= sigmahat;
+return;
+}
+if (cb=='l')
+{ lik = dointpoint(lf,x,PLIK,ev,i);
+rdf = dointpoint(lf,x,PRDF,ev,i);
+sef[i] *= sqrt(-2*lik/rdf);
+return;
+}
+if (cb=='p')
+{ sef[i] = sigmahat*sqrt(1+sef[i]*sef[i]);
+return;
+}
+}
+void predptdir(lf,des,x,what,i)
+lfit *lf;
+design *des;
+double *x;
+int what, i;
+{ int needcv;
+des->xev = x;
+needcv = (what==PVARI) | (what==PNLX) | (what==PT0) | (what==PRDF);
+locfit(&lf->lfd,des,&lf->sp,0,1,needcv);
+switch(what)
+{ case PCOEF: fit[i] = des->cf[0]; break;
+case PVARI: fit[i] = des->V[0]; break;
+case PNLX:  fit[i] = sqrt(des->V[0]); break;
+case PT0:   fit[i] = des->f1[0]; break;
+case PBAND: fit[i] = des->h; break;
+case PDEGR: fit[i] = deg(&lf->sp); break;
+case PLIK:  fit[i] = des->llk; break;
+case PRDF:  fit[i] = des->tr0 - des->tr2; break;
+default: LERR(("unknown what in predptdir"));
+}
+}
+void prepvector(lf,des,x,n,what,dir) /* interpolate a vector */
+lfit *lf;
+design *des;
+double **x;
+int n, what, dir;
+{ int i, j;
+double xx[MXDIM];
+for (i=0; i<n; i++)
+{ for (j=0; j<lf->fp.d; j++) xx[j] = x[j][i];
+if (dir)
+predptdir(lf,des,xx,what,i);
+else
+predptall(lf,xx,what,ev(&lf->evs),i);
+if (lf_error) return;
+}
+}
+void prepfitp(lf,what)
+lfit *lf;
+int what;
+{ int  i;
+for (i=0; i<lf->fp.nv; i++)
+{ predptall(lf,evpt(&lf->fp,i),what,EFITP,i);
+if (lf_error) return;
+}
+}
+void prepgrid(lf,des,x,mg,n,what,dir) /* interpolate a grid given margins */
+lfit *lf;
+design *des;
+double **x;
+int *mg, dir, n, what;
+{ int i, ii, j, d;
+double xv[MXDIM];
+d = lf->fp.d;
+for (i=0; i<n; i++)
+{ ii = i;
+for (j=0; j<d; j++)
+{ xv[j] = x[j][ii%mg[j]];
+ii /= mg[j];
+}
+if (dir)
+predptdir(lf,des,xv,what,i);
+else
+predptall(lf,xv,what,ev(&lf->evs),i);
+if (lf_error) return;
+}
+}
+void preplot(lf,x,f,se,band,mg,where,what,dir)
+lfit *lf;
+double **x, *f, *se;
+int *mg, where, what, dir;
+char band;
+{ int d, i, n;
+double *xx[MXDIM];
+design ppdes;
+d = lf->fp.d;
+fit = f;
+sef = se;
+cb = band;
+if (cb!='n') sigmahat = sqrt(rv(&lf->fp));
+if (dir) des_init(&ppdes,lf->lfd.n,npar(&lf->sp));
+switch(where)
+{ case 1: /* vector */
+n = mg[0];
+prepvector(lf,&ppdes,x,n,what,dir);
+break;
+case 2: /* grid */
+n = 1;
+for (i=0; i<d; i++) n *= mg[i];
+prepgrid(lf,&ppdes,x,mg,n,what,dir);
+break;
+case 3: /* data */
+n = lf->lfd.n;
+if ((ev(&lf->evs)==EDATA) | (ev(&lf->evs)==ECROS))
+{ prepfitp(lf,what);
+dir = 0;
+}
+else
+{ for (i=0; i<d; i++) xx[i] = dvari(&lf->lfd,i);
+prepvector(lf,&ppdes,xx,n,what,dir);
+}
+break;
+case 4: /* fit points */
+n = lf->fp.nv; dir = 0;
+prepfitp(lf,what);
+break;
+default:
+LERR(("unknown where in preplot"));
+}
+if ((!dir) && ((what==PT0)|(what==PVARI)))
+for (i=0; i<n; i++) f[i] = f[i]*f[i];
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+#include "lfev.h"
+int procv_nov(des,lf,v)
+design *des;
+lfit *lf;
+int v;
+{ int lf_status;
+if (lf_debug>1) mut_printf(" procveraw: %d\n",v);
+des->xev = evpt(&lf->fp,v);
+if (acri(&lf->sp)==ANONE)
+lf_status = locfit(&lf->lfd,des,&lf->sp,0,1,0);
+else
+lf_status = alocfit(&lf->lfd,&lf->sp,&lf->dv,des,0);
+if (lf->fp.h != NULL) lf->fp.h[v] = des->h;
+return(lf_status);
+}
+int procv_var(des,lf,v)
+design *des;
+lfit *lf;
+int v;
+{ int i, lf_status;
+if (lf_debug>1) mut_printf(" procvraw: %d\n",v);
+des->xev = evpt(&lf->fp,v);
+if (acri(&lf->sp)==ANONE)
+lf_status = locfit(&lf->lfd,des,&lf->sp,0,1,1);
+else
+lf_status = alocfit(&lf->lfd,&lf->sp,&lf->dv,des,1);
+if (lf->fp.h != NULL) lf->fp.h[v] = des->h;
+return(lf_status);
+}
+/*
+* Copyright 1996-2006 Catherine Loader.
+*/
+/*
+* startmodule(lf,des,mod,dir) -- the standard entry point.
+*   des and lf are pointers to the design and fit structures.
+*   mod - module name. Set to NULL if the module is already
+*         initialized.
+*   dir - for dynamic modules, the directory.
+*
+* initmodule(mdl,mod,dir,lf)
+*   direct call for module initialization.
+*
+*/
+#include "lfev.h"
+#ifdef WINDOWS
+#define DIRSEP '\\'
+#define PATHSEP ';'
+#else
+#define DIRSEP '/'
+#define PATHSEP ':'
+#endif
+#ifdef ALLOW_MODULES
+#ifdef WINDOWS
+#include <windows.h>
+#define DLEXT "dll"
+#define DLOPEN(x) LoadLibrary(x)
+#define DLSYM GetProcAddress
+#else
+#include <dlfcn.h>
+#define DLEXT "so"
+#define DLOPEN(x) dlopen(x,RTLD_LAZY)
+#define DLSYM dlsym
+#endif
+#endif
+static double fpkap[6];
+void fitpt_init(fp)
+fitpt *fp;
+{
+dc(fp) = 0;
+geth(fp) = GSTD;
+fp->nv = fp->nvm = 0;
+if (fp->kap==NULL) fp->kap = fpkap;
+}
+void lfit_init(lf)
+lfit *lf;
+{
+lfdata_init(&lf->lfd);
+evstruc_init(&lf->evs);
+smpar_init(&lf->sp,&lf->lfd);
+deriv_init(&lf->dv);
+fitpt_init(&lf->fp);
+lf->mdl.np = 0;
+}
+void fitdefault(lf)
+lfit *lf;
+{ WARN(("fitdefault deprecated -- use lfit_init()"));
+lfit_init(lf);
+}
+void set_flim(lfd,evs)
+lfdata *lfd;
+evstruc *evs;
+{ int i, j, d, n;
+double z, mx, mn, *bx;
+if (ev(evs)==ESPHR) return;
+d = lfd->d; n = lfd->n;
+bx = evs->fl;
+for (i=0; i<d; i++)
+if (bx[i]==bx[i+d])
+{ if (lfd->sty[i]==STANGL)
+{ bx[i] = 0.0; bx[i+d] = 2*PI*lfd->sca[i];
+}
+else
+{ mx = mn = datum(lfd,i,0);
+for (j=1; j<n; j++)
+{ mx = MAX(mx,datum(lfd,i,j));
+mn = MIN(mn,datum(lfd,i,j));
+}
+if (lfd->xl[i]<lfd->xl[i+d]) /* user set xlim; maybe use them. */
+{ z = mx-mn;
+if (mn-0.2*z < lfd->xl[i]) mn = lfd->xl[i];
+if (mx+0.2*z > lfd->xl[i+d]) mx = lfd->xl[i+d];
+}
+bx[i] = mn;
+bx[i+d] = mx;
+}
+}
+}
+double vvari(v,n)
+double *v;
+int n;
+{ int i;
+double xb, s2;
+xb = s2 = 0.0;
+for (i=0; i<n; i++) xb += v[i];
+xb /= n;
+for (i=0; i<n; i++) s2 += SQR(v[i]-xb);
+return(s2/(n-1));
+}
+void set_scales(lfd)
+lfdata *lfd;
+{ int i;
+for (i=0; i<lfd->d; i++)
+if (lfd->sca[i]<=0) /* set automatic scales */
+{ if (lfd->sty[i]==STANGL)
+lfd->sca[i] = 1.0;
+else lfd->sca[i] = sqrt(vvari(lfd->x[i],lfd->n));
+}
+}
+void nstartlf(des,lf)
+design *des;
+lfit *lf;
+{ int i, d, n;
+if (lf_debug>0) mut_printf("nstartlf\n");
+n = lf->lfd.n;
+d = lf->lfd.d;
+npar(&lf->sp) = calcp(&lf->sp,d);
+des_init(des,n,npar(&lf->sp));
+set_scales(&lf->lfd);
+set_flim(&lf->lfd,&lf->evs);
+compparcomp(des,&lf->lfd,&lf->sp,&lf->pc,lf->mdl.nopc);
+if (lf_error) return;
+makecfn(&lf->sp,des,&lf->dv,lf->lfd.d);
+lf->lfd.ord = 0;
+if ((d==1) && (lf->lfd.sty[0]!=STANGL))
+{ i = 1;
+while ((i<n) && (datum(&lf->lfd,0,i)>=datum(&lf->lfd,0,i-1))) i++;
+lf->lfd.ord = (i==n);
+}
+for (i=0; i<npar(&lf->sp); i++) des->fix[i] = 0;
+lf->fp.d = lf->lfd.d;
+lf->fp.hasd = (des->ncoef==(1+lf->fp.d));
+lf->fp.nv = lf->evs.nce = 0;
+if (lf_debug>1) mut_printf("call eval structure %d\n",ev(&lf->evs));
+switch(ev(&lf->evs))
+{ case EPHULL: triang_start(des,lf); break;
+case EDATA:  dataf(des,lf); break;
+case ECROS:  crossf(des,lf); break;
+case EGRID:  gridf(des,lf); break;
+case ETREE:  atree_start(des,lf); break;
+case EKDCE:  kt(&lf->sp) = KCE;
+case EKDTR:  kdtre_start(des,lf); break;
+case EPRES:  preset(des,lf); break;
+case EXBAR:  xbarf(des,lf); break;
+case ENONE:  return;
+case ESPHR:  sphere_start(des,lf); break;
+case ESPEC:  lf->evs.espec(des,lf); break;
+default: LERR(("startlf: Invalid evaluation structure %d",ev(&lf->evs)));
+}
+/* renormalize for family=density */
+if ((de_renorm) && (fam(&lf->sp)==TDEN)) dens_renorm(lf,des);
+}
+/*
+* getnextdir() gets the next dir from a string dirpath="dir1:dir2:dir3:..."
+*   (;-separated on windows).
+*   The directory is returned through dirnext, and the function returns
+*   a pointer to the next string.
+*   typical usage is recursive, dirpath = getnextdir(dirpath,dirnext).
+*   with the above example, this sets dirnext="dir1" and dirpath="dir2:dir3:...".
+* if the input dirpath has no :, then it is copied to dirnext, and return is "".
+* if input dirpath is "", dirnext is set to "", and null pointer returned.
+*/
+char *getnextdir(dirpath,dirnext)
+char *dirpath, *dirnext;
+{ char *z;
+if (strlen(dirpath)==0)
+{ sprintf(dirnext,"");
+return(NULL);
+}
+z = strchr(dirpath,PATHSEP);
+if (z==NULL)
+{ sprintf(dirnext,"%s%c",dirpath,DIRSEP);
+return(&dirpath[strlen(dirnext)]);
+}
+*z = '\0';
+sprintf(dirnext,"%s%c",dirpath,DIRSEP);
+return(++z);
+}
+int initmodule(mdl, mod, dir, lf)
+module *mdl;
+lfit *lf;
+char *mod, *dir;
+{ int n, d, p;
+#ifdef ALLOW_MODULES
+#ifdef WINDOWS
+HINSTANCE res;
+typedef void (CALLBACK* DLLFN)();
+DLLFN init;
+#else
+void *res;
+void (*init)();
+#endif
+char distname[500];
+#endif
+n = lf->lfd.n;
+d = lf->lfd.d;
+p = npar(&lf->sp) = calcp(&lf->sp,d);
+mdl->isset = 1;
+PPROC(lf) = NULL;
+if (strcmp(mod,"std")==0)    { initstd(lf); return(1); }
+if (strcmp(mod,"simple")==0) { initsimple(lf); return(1); }
+if (strcmp(mod,"allcf")==0)  { initallcf(lf); return(1); }
+if (strcmp(mod,"hatm")==0)   { inithatm(lf); return(1); }
+if (strcmp(mod,"kappa")==0)  { initkappa(lf); return(1); }
+if (strcmp(mod,"lscv")==0)   { initlscv(lf); return(1); }
+if (strcmp(mod,"gamf")==0)   { initgam(lf); return(1); }
+if (strcmp(mod,"gamp")==0)   { initgam(lf); return(1); }
+if (strcmp(mod,"rband")==0)  { initrband(lf); return(1); }
+if (strcmp(mod,"scb")==0)    { initscb(lf); return(1); }
+if (strcmp(mod,"vord")==0)   { initvord(lf); return(1); }
+#ifdef ALLOW_MODULES
+while (dir != NULL)
+{
+dir = getnextdir(dir,distname);
+sprintf(&distname[strlen(distname)],"mod%s.%s",mod,DLEXT);
+res = DLOPEN(distname);
+if (res==NULL)
+{
+#ifdef WINDOWS
+mut_printf("LoadLibrary failed: %s, %d\n",distname,GetLastError());
+#else
+mut_printf("dlopen failed: %s\n",dlerror());
+#endif
+}
+else
+{
+#ifdef WINDOWS
+mut_printf("LoadLibrary success: %s\n",distname);
+#else
+mut_printf("dlopen success: %s\n",distname);
+#endif
+sprintf(distname,"init%s",mod);
+init = (void *)DLSYM(res,distname);
+if (init==NULL)
+{ mut_printf("I can't find %s() function.\n",distname);
+mdl->isset = 0;
+return(0);
+}
+init(lf);
+return(1);
+}
+}
+#endif
+mdl->isset = 0;
+return(0);
+}
+/*
+* startmodule is the entry point to launch the fit.
+* if mod is provided, will first initialize the module.
+* if mod==NULL, assumes the module has been initialized separately.
+*/
+void startmodule(lf,des,mod,dir)
+lfit *lf;
+design *des;
+char *mod, *dir;
+{ int z;
+if (mod != NULL)
+{ z = initmodule(&lf->mdl,mod,dir,lf);
+if (!z) return;
+}
+lf->fp.nv = lf->evs.nce = 0;
+if (lf_error) return;
+if (PROCV(lf) != NULL) nstartlf(des,lf);
+if (lf_error) return;
+if (PPROC(lf) != NULL) PPROC(lf)(lf,des,lf->fp.kap);
+}
+/* for compatability, more or less. */
+void startlf(des,lf,vfun,nopc)
+design *des;
+lfit *lf;
+int (*vfun)(), nopc;
+{ int z;
+z = initmodule(&lf->mdl,"std",NULL,lf);
+if (!z) return;
+lf->mdl.procv = vfun;
+lf->mdl.nopc = nopc;
+nstartlf(des,lf);
+}

Mercurial > repos > vipints > rdiff

comparison rDiff/src/locfit/Source/liblfev.c @ 0:0f80a5141704