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

version 0.3 uploaded

0: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);
}