/* Copyright (C) 1986-2000 Simon N. Wood  snw@st-and.ac.uk

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License   
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
(www.gnu.org/copyleft/gpl.html)

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
USA.*/


/* various spline routines */


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "matrix.h"
#include "spline.h"
#define SQR(a) ((a)*(a))
#define CUBE(a) ((a)*(a)*(a))
#define QUAD(a) ((a)*(a)*(a)*(a))


/* The next 4 functions are basis functions for the 1st derivative
   representation of a cubic spline */

void ErrorMessage(char *msg,int fatal);

double b0(x0,x1,x) double x0,x1,x;

/* multiplies function value at x0 */

{ double res,h,xx1;
  h=x1-x0;xx1=x-x1;
  res=2.0*(x-x0+0.5*h)*xx1*xx1/(h*h*h);
  return(res);
}

double b1(x0,x1,x) double x0,x1,x;

/* multiplies function value at x1 */

{ double res,h,xx0;
  h=x1-x0;xx0=x-x0;
  res= -2.0*(x-x1-0.5*h)*xx0*xx0/(h*h*h);
  return(res);
}

double d0(x0,x1,x) double x0,x1,x;

/* multiplies gradient at x0 */

{ double res,h,xx1;
  h=x1-x0;xx1=x-x1;
  res=(x-x0)*xx1*xx1/(h*h);
  return(res);
}

double d1(x0,x1,x) double x0,x1,x;

/* multiplies gradient at x1 */

{ double res,h,xx0;
  h=x1-x0;xx0=x-x0;
  res=xx0*xx0*(x-x1)/(h*h);
  return(res);
}


/* The next 4 functions are derivatives of the spline basis functions used
   above. */

double db0(x0,x1,x) double x0,x1,x;

{ double res,h,xx1;
  h=x1-x0;xx1=x-x1;
  res=2.0*(2.0*(x-x0+0.5*h)*xx1+xx1*xx1)/(h*h*h);
  return(res);
}
double db1(x0,x1,x) double x0,x1,x;

{ double res,h,xx0;
  h=x1-x0;xx0=x-x0;
  res= -2.0*(2.0*(x-x1-0.5*h)*xx0+xx0*xx0)/(h*h*h);
  return(res);
}

double dd0(x0,x1,x) double x0,x1,x;

{ double res,h,xx1;
  h=x1-x0;xx1=x-x1;
  res=(xx1*xx1+2.0*(x-x0)*xx1)/(h*h);
  return(res);
}

double dd1(x0,x1,x) double x0,x1,x;

{ double res,h,xx0;
  h=x1-x0;xx0=x-x0;
  res=(xx0*xx0+2.0*xx0*(x-x1))/(h*h);
  return(res);
}

/* the next 4 functions are the integrals of the basis functions given above
   from a to b */

double intb0(x0,x1,a,b) double x0,x1,a,b;

{ double h,r;
  h=x1-x0;
  r=2.0*((QUAD(b)-QUAD(a))/4.0
	     +(CUBE(b)-CUBE(a))*(h/2.0-2.0*x1-x0)/3.0
	     +(SQR(b)-SQR(a))*(SQR(x1)+2.0*x1*x0-h*x1)/2.0
	     +(b-a)*(h*SQR(x1)/2.0-SQR(x1)*x0))/CUBE(h);
  return(r);
}



double intb1(x0,x1,a,b) double x0,x1,a,b;

{ double h,r;
  h=x1-x0;
  r= -2*((QUAD(b)-QUAD(a))/4.0-
	      (CUBE(b)-CUBE(a))*(h/2.0+2.0*x0+x1)/3.0
		+(SQR(b)-SQR(a))*(SQR(x0)+2.0*x1*x0+h*x0)/2.0
		-(b-a)*(h*SQR(x0)/2.0+SQR(x0)*x1))/CUBE(h);
  return(r);
}



double intd0(x0,x1,a,b) double x0,x1,a,b;

{ double h,r;
  h=x1-x0;
  r=((QUAD(b)-QUAD(a))/4.0-(CUBE(b)-CUBE(a))*(2.0*x1+x0)/3.0+
      (SQR(b)-SQR(a))*(SQR(x1)+2.0*x0*x1)/2.0-(b-a)*x0*SQR(x1))/SQR(h);
  return(r);
}


double intd1(x0,x1,a,b) double x0,x1,a,b;

{ double h,r;
  h=x1-x0;
  r=((QUAD(b)-QUAD(a))/4.0-(CUBE(b)-CUBE(a))*(2.0*x0+x1)/3.0+
      (SQR(b)-SQR(a))*(SQR(x0)+2.0*x0*x1)/2.0-(b-a)*x1*SQR(x0))/SQR(h);
  return(r);
}


matrix getD(h,nak) matrix h;int nak;

/* the matrix mapping the value of the spline to the gradients at the knots.
   nak is true for 'not-a-knot' end conditions at the early end, otherwise
   'natural' end conditions are used. If there are only 2 knots then the spline
   is taken as a straight line if only 1 a constant. */

{ long i,j,n;
  matrix T,D,Res;
  n=h.r+1;
  T=initmat(n,n);D=initmat(n,n);Res=initmat(n,n);
  for (i=0;i<n;i++) for (j=0;j<n;j++)
  { T.M[i][j]=0.0;D.M[i][j]=0.0;}
  if (n==1L)
  { Res.M[0][0]=0.0;
  } else
  if (n==2L)
  { Res.M[0][0]=Res.M[1][0]=-1.0/h.V[0];
    Res.M[0][1]=Res.M[1][1]=1.0/h.V[0];
  } else
  { for (i=0;i<n;i++) T.M[i][i]=2.0;
    for (i=1;i<n-1;i++)
    { T.M[i][i-1]=h.V[i]/(h.V[i]+h.V[i-1]);
      T.M[i][i+1]=1.0-T.M[i][i-1];
      D.M[i][i-1]= -3.0*T.M[i][i-1]/h.V[i-1];
      D.M[i][i+1]=3.0*T.M[i][i+1]/h.V[i];
      D.M[i][i]= -(D.M[i][i+1]+D.M[i][i-1]);
    }
    if (!nak)
    { T.M[0][1]=1.0;D.M[0][0]= -3.0/h.V[0];D.M[0][1]= -D.M[0][0];}
    else
    { T.M[0][1]=2.0*(h.V[0]+h.V[1])/h.V[1];
      D.M[0][0]= -2.0*(3.0*h.V[0]+2.0*h.V[1])/
		(h.V[0]*(h.V[0]+h.V[1]));
      D.M[0][2]=2.0*h.V[0]*h.V[0]/
      (h.V[1]*h.V[1]*(h.V[0]+h.V[1]));
      D.M[0][1]= -D.M[0][0]-D.M[0][2];
    }
    T.M[n-1][n-2]=1.0;D.M[n-1][n-2]= -3.0/h.V[n-2];
    D.M[n-1][n-1]= -D.M[n-1][n-2];
    invert(&T);
    matmult(Res,T,D,0,0);
  }
  freemat(T);freemat(D);
  return(Res);
}

double bsbf(int m,double *k,double x,int i,int n)

/* recursive calculation of b-spline basis functions -
   see, for example, Lancaster and Salkauskas "Curve and Surface Fitting" 1986
   or Chui "Multivariate Splines" 1988 or deBoor 1978 */

{ double z0,z1;
  int j1,j0;
  if (m==0) /* then recursion has reached the Oth order spline */
  { if (i==n-2&&(x>=k[i]&&x<=k[i+1])) return(1.0);
    if (i>n-2||i<0) return(0.0);
    if (x<k[i]||x>=k[i+1]) return(0.0);
    else return(1.0);
  } else
  { j1=i+m;if (j1>n-1) j1=n-1;j0=i;if (j0<0) j0=0;
    if (j0>=j1) z0=0.0;else z0=(x-k[j0])/(k[j1]-k[j0]);
    j1=i+m+1;if (j1>n-1) j1=n-1;j0=i+1;if (j0<0) j0=0;
    if (j0>=j1) z1=0.0;else z1=(k[j1]-x)/(k[j1]-k[j0]);
    return(z0*bsbf(m-1,k,x,i,n)+z1*bsbf(m-1,k,x,i+1,n));
  }
}

int bsmap(double *b,double *k,double x,int n)

/* given a knot sequence k_0 ..... k_{n-1} this routine returns the 4 non-zero
   cubic b splines basis functions at x - these are returned in b. The routine
   also returns the index of the last basis function in the series of 4.....
   so s(x) = b[0]*p_{i-3} + b[1]*p_{i-2} + b[2]*p_{i-1} + b[3]*p_i
   where p is the parameter vector, which has 3 more coefficients than there are
   knots. This routine is not suitable for extrapolation beyond [k_0,k_{n-1}].
*/
{ static ik=0L;
  int i;
  /* locate the knot interval */
  if (x<k[0]||x>k[n-1]) ErrorMessage("You are trying to extrapolate with bsmap()",1);
  while (k[ik]>x) ik--; while (ik<n-1&&k[ik+1]<=x) ik++;
  for (i=0;i<4;i++) b[i]=bsbf(3,k,x,ik-3+i,n);
  return(ik);
}


void rkmap(matrix k,matrix c, double x)

/* uses spline basis:

   s(x) = p_0+p_1 x + \sum_{i=2}^{n-1} p_i [(x-k_{i-2})^+]^3

   where n= number of knots+2 - the total dimension of the paramter vector p.
   Returns coefficient vector c, such that s(x)=c'p.....
   ( Appears to be what U uses in bs(), but due to excellent nature of
     documentation, I'm not certain of this.)
*/

{ int i;
  double z;
  c.V[0]=1.0;c.V[1]=x;
  for (i=2;i<c.r;i++)
  { z=x-k.V[i-2];z=fabs(z);c.V[i]=z*z*z;}
}



void imap(im,t,a,b,kill)
matrix im,t;double a,b;int kill;

/* im maps the function values of the spline at the t.M[i]s to the integral
   of the spline from a to b. setting kill=1 clears this routines internal
   storage after use. Set kill =1 BEFORE a change in the vector t, which
   should contain the knot positions. im must be initialized before calling
   ( im =initmat(t.r,1L) }*/

{ static matrix D;static char first=1;
  matrix h;
  long i,j,ia,ib;
  if (first)
  { first=0;h=initmat(t.r-1,1L);
    for (i=0L;i<t.r-1;i++) h.V[i]=t.V[i+1]-t.V[i];
    D=getD(h,0); /* time trajectories always have natural end conditions */
    freemat(h);
  }
  ia=0L;while((a>t.M[ia+1][0])&&(ia<t.r-2)) ia++;
  ib=0L;while((b>t.M[ib+1][0])&&(ib<t.r-2)) ib++;
  if (ia==ib)
  { for (i=0;i<im.r;i++)
    im.V[i]=D.M[ia][i]*intd0(t.V[ia],t.V[ia+1],a,b)+
	    D.M[ia+1][i]*intd1(t.V[ia],t.V[ia+1],a,b);
    im.V[ia]+=intb0(t.V[ia],t.V[ia+1],a,b);
    im.V[ia+1]+=intb1(t.V[ia],t.V[ia+1],a,b);
  } else
  { for (i=0;i<im.r;i++)
    im.V[i]=D.M[ia][i]*intd0(t.V[ia],t.V[ia+1],a,t.V[ia+1])+
	    D.M[ia+1][i]*intd1(t.V[ia],t.V[ia+1],a,t.V[ia+1]);
    im.V[ia]+=intb0(t.V[ia],t.V[ia+1],a,t.V[ia+1]);
    im.V[ia+1]+=intb1(t.V[ia],t.V[ia+1],a,t.V[ia+1]);
    for (j=ia+1;j<ib;j++)
    { for (i=0;i<im.r;i++)
      im.V[i]+=D.M[j][i]*intd0(t.V[j],t.V[j+1],t.V[j],t.V[j+1])+
      D.M[j+1][i]*intd1(t.V[j],t.V[j+1],t.V[j],t.V[j+1]);
      im.V[j]+=intb0(t.V[j],t.V[j+1],t.V[j],t.V[j+1]);
      im.V[j+1]+=intb1(t.V[j],t.V[j+1],t.V[j],t.V[j+1]);
    }
    for (i=0;i<im.r;i++)
    im.V[i]+=D.M[ib][i]*intd0(t.V[ib],t.V[ib+1],t.V[ib],b)+
	      D.M[ib+1][i]*intd1(t.V[ib],t.V[ib+1],t.V[ib],b);
    im.V[ib]+=intb0(t.V[ib],t.V[ib+1],t.V[ib],b);
    im.V[ib+1]+=intb1(t.V[ib],t.V[ib+1],t.V[ib],b);
  }
  if (kill)
  { freemat(D);first=1;}
}

void tmap(tm,tgm,t,time,kill)
matrix tm,tgm,t;
double time;
int kill; /* to release static matrix allocation set to 1 otherwise 0 and
	     prepare for a new sequence of knot positions in t*/

/* tm maps values of a function at the t values contained in vector t to
   the value of a spline through those points at 'time' ;tgm does the same
   for the gradient of the spline */

{ static matrix D;static char first=1;
  matrix h;
  long i,k;
  if (first)
  { first=0;h=initmat(t.r-1,1L);
    for (i=0L;i<t.r-1;i++) h.V[i]=t.V[i+1]-t.V[i];
    D=getD(h,0); /* time trajectories always have natural end conditions */
    freemat(h);
  }
  if (t.r==1L)
  { tm.V[0]=1.0;tgm.V[0]=0.0;}
  else
  { i=0L;while((time>t.V[i+1])&&(i<t.r-2)) i++;
    for (k=0;k<t.r;k++)
    tm.V[k]=D.M[i][k]*d0(t.V[i],t.V[i+1],time)+
	    D.M[i+1][k]*d1(t.V[i],t.V[i+1],time);
    tm.V[i]+=b0(t.V[i],t.V[i+1],time);
    tm.V[i+1]+=b1(t.V[i],t.V[i+1],time);
    for (k=0;k<t.r;k++)
    tgm.V[k]=D.M[i][k]*dd0(t.V[i],t.V[i+1],time)+
	     D.M[i+1][k]*dd1(t.V[i],t.V[i+1],time);
    tgm.V[i]+=db0(t.V[i],t.V[i+1],time);
    tgm.V[i+1]+=db1(t.V[i],t.V[i+1],time);
  }
  if (kill)
  { first=1;
    freemat(D);
  }
}

void poscon(matrix *A,matrix *b,matrix x,long nop,int perknot)

/* Gets a set of matrix constraints ensuring non-negativity of a spline
   with knots at the elements of x at each of a set of nop evenly spaced
   points. If perknot=1 then nop is constraints per section of cubic,
   otherwise it's the total number of evenly spaced constraints */

{ double dt;
  long i,j,totn,k;
  matrix xm,xmg;
  if (perknot) totn=(x.r-1)*nop+1;else totn=nop;
  (*A)=initmat(totn,x.r);
  xm=initmat(x.r,1L);
  xmg=initmat(x.r,1L);
  tmap(xm,xmg,x,x.V[0],1); /* resetting tmap routine */
  if (perknot)
  { for (i=0;i<x.r-1;i++)
    { dt=(x.V[i+1]-x.V[i])/nop;
      for (k=0;k<nop;k++)
      { tmap(xm,xmg,x,x.V[i]+k*dt,0);
        for (j=0;j<x.r;j++) A->M[i*nop+k][j]=xm.V[j];
      }
    }
    tmap(xm,xmg,x,x.V[x.r-1],0);
    for (j=0;j<x.r;j++) A->M[A->r-1][j]=xm.V[j];
  }
  else
  { dt=(x.V[x.r-1]-x.V[0])/(nop-1);
    for (i=0;i<nop;i++)
    { tmap(xm,xmg,x,x.V[0]+i*dt,0);
      for (j=0;j<x.r;j++) A->M[i][j]=xm.V[j];
    }
  }
  tmap(xm,xmg,x,x.V[0],1); /* resetting tmap routine */
  (*b)=initmat(totn,1L);
  freemat(xm);freemat(xmg);
}




void getHBH(HBH,h,nak,rescale) matrix *HBH,h;int nak,rescale;

/* Generates the wiggliness measure matrix for vector h; nak=0 for natural
   end conditions or nak=1 to use the not a knot condition at the lower end;
   set rescale=1 to produce a measure rescaled for the unit interval, set to
   zero otherwise */

{ long n,i,j;
  matrix C,B,BI,H,hn;
  double interval=0.0;
  n=h.r;
  if (rescale)
  { for (i=0;i<h.r;i++) interval+=h.V[i];
    hn=initmat(h.r,1L);
    for (i=0;i<h.r;i++) hn.V[i]=h.V[i]/interval;
  } else hn=h;
  (*HBH)=initmat(n+1,n+1);
  if (!nak)
  { C=initmat(n-1,n+1);
    B=initmat(n-1,n-1);
    H=initmat(n-1,n+1);
    for (i=0;i<n-1;i++)
    { for (j=0;j<n-1;j++)
      { B.M[i][j]=0.0;
	     H.M[i][j]=0.0;
      }
      H.M[i][n-1]=0.0;
      H.M[i][n]=0.0;
    }
    for (i=0;i<n-1;i++)
    { B.M[i][i]=(hn.V[i]+hn.V[i+1])/3.0;
      H.M[i][i]=1.0/hn.V[i];
      H.M[i][i+1]= -1.0/hn.V[i]-1.0/hn.V[i+1];
      H.M[i][i+2]=1.0/hn.V[i+1];
    }
    for (i=0;i<n-2;i++)
    { B.M[i][i+1]=hn.V[i+1]/6.0;
      B.M[i+1][i]=hn.V[i+1]/6.0;
    }
    invert(&B);
    matmult(C,B,H,0,0);
    matmult((*HBH),H,C,1,0);
    freemat(C);freemat(B);freemat(H);
  } else
  { H=initmat(n,n+1);
    BI=initmat(n,n);B=initmat(n,n);
    for (i=0;i<H.r;i++) for (j=0;j<H.c;j++) H.M[i][j]=0.0;
    for (i=1;i<n;i++)
    { H.M[i][i-1]=1.0/hn.V[i-1];H.M[i][i]= -1.0/hn.V[i-1]-1.0/hn.V[i];
      H.M[i][i+1]=1.0/hn.V[i];
    }
    for (i=0;i<n;i++) for (j=0;j<n;j++)
    { BI.M[i][j]=0.0;B.M[i][j]=0.0;}
    for (i=1;i<n;i++)
    { B.M[i][i-1]=hn.V[i-1]/6.0;B.M[i][i]=(hn.V[i-1]+hn.V[i])/3.0;
      if (i<(n-1))
      { B.M[i][i+1]=hn.V[i]/6.0;
	BI.M[i][i+1]=B.M[i][i+1];
      }
      for (j=0;j<2;j++) BI.M[i][j+i-1]=B.M[i][j+i-1];
    }
    B.M[0][0]= -hn.V[1];B.M[0][1]=hn.V[0]+hn.V[1];
    B.M[0][2]= -hn.V[0];
    BI.M[0][0]=hn.V[0]/3.0;BI.M[0][1]=hn.V[0]/6.0;
    C=initmat(n,n);
    invert(&B);
    matmult(C,BI,B,0,0);
    matmult(BI,B,C,1,0);
    freemat(B);freemat(C);
    C=initmat(n,n+1);
    matmult(C,BI,H,0,0);
    matmult((*HBH),H,C,1,0);
    freemat(C);freemat(BI);freemat(H);
  }
  if (rescale) freemat(hn);
}


void getSmooth(S,x,rescale) matrix *S,x;int rescale;

/* gets a natural wiggliness measure for a spline with knots at the elements
   of vector x. Set rescale not zero to pretend that the domain is the unit
   interval. */

{ matrix h;
  long i;
  h=initmat(x.r-1L,1L);
  for (i=0;i<x.r-1;i++) h.V[i]=x.V[i+1]-x.V[i];
  getHBH(S,h,0,rescale);
  freemat(h);
}

void createA(D,A,h,up,pos) matrix D,*A,h;int up,pos;

/* gets matrix A for constraints of the form Ay>=b ensuring monotonic

   change  of the cubic spline interpolating (x ,y ) where h = x   -x
					       i  i         i   i+1   i
   D is the matrix mapping y to gradients at x, which must be supplied to the

   routine along with h. Up should be non-zero for increase, zero for

   decrease. pos should be non-zero for a non-negative spline.  */

{ long i,j,n;
  double m;
  n=h.r;
  if (up) m= -1.0; else m=1.0;
  if (pos) (*A)=initmat(4*n+1,n+1); else (*A)=initmat(4*n,n+1);
  for (i=0;i<n;i++)
  { for (j=0;j<n+1;j++)
    { if (j==i)
      { A->M[i][j]=(D.M[i][j]+3.0/h.V[i])*m;   /**not certain of d.M update**/
	A->M[i+n][j]=(D.M[i+1][j]+3.0/h.V[i])*m;
	A->M[i+2*n][j]=m;
	A->M[i+3*n][j]= -D.M[i][j]*m;
      } else
      if (j==(i+1))
      { A->M[i][j]=(D.M[i][j]-3.0/h.V[i])*m;
	A->M[i+n][j]=(D.M[i+1][j]-3.0/h.V[i])*m;
	A->M[i+2*n][j]= -m;
	A->M[i+3*n][j]= -D.M[i][j]*m;
      } else
      { A->M[i][j]=D.M[i][j]*m;
	A->M[i+n][j]=D.M[i+1][j]*m;
	A->M[i+2*n][j]=0.0;
	A->M[i+3*n][j]= -D.M[i][j]*m;
      }
    }
  }
  if (pos)
  { for (j=0;j<n+1;j++) A->M[4*n][j]=0.0;
    if (up) A->M[4*n][0]=1.0; else A->M[4*n][n]=1.0;
  }
}


double spline(x,p,g,z)
matrix x,p,g;
double z;

/* evaluates the spline with knots at the positions in x, values at the knots
   in p and gradients at the knots in g at position z */

{ static int i=0; /* x.M[i]<=z<x.M[i+1]*/
  double res,*xV,*gV,*pV;
  xV=x.V;pV=p.V;gV=g.V; /* 2/2/97 */
  while ((xV[i]>z)&&(i>0)) i--;
  while ((xV[i+1]<=z)&&(i<x.r-2)) i++;
  res=b0(xV[i],xV[i+1],z)*pV[i]+b1(xV[i],xV[i+1],z)*pV[i+1]+
      d0(xV[i],xV[i+1],z)*gV[i]+d1(xV[i],xV[i+1],z)*gV[i+1];
  return(res);
}

double gspline(x,p,g,z)
matrix x,p,g;
double z;

/* evaluates the spline with knots at the positions in x, values at the knots
   in p and gradients at the knots in g at position z */

{ static int i=0; /* x.M[i]<=z<x.M[i+1]*/
  double res,*xV,*gV,*pV;
  xV=x.V;pV=p.V;gV=g.V; /* 2/2/97 */
  while ((xV[i]>z)&&(i>0)) i--;
  while ((xV[i+1]<=z)&&(i<x.r-2)) i++;
  res=db0(xV[i],xV[i+1],z)*pV[i]+db1(xV[i],xV[i+1],z)*pV[i+1]+
      dd0(xV[i],xV[i+1],z)*gV[i]+dd1(xV[i],xV[i+1],z)*gV[i+1];
  return(res);
}

void hyman_filter(x,a,b) matrix x,a,b;

/* Filters spline function to yield co-monotonicity in accordance with
   Hyman (1983) SIAM J. Sci. Stat. Comput. 4(4):645-654, x is knot position
   a is value at knot b is gradient at knot. See also Dougherty, Edelman and
   Hyman 1989 Mathematics of Computation 52: 471-494. Co-monotonicity logic
   corrected 28/1/99: sig calculation incorrect.*/

{ double *S,*aV,*xV,*bV,sig;
  int i;
  aV=a.V;bV=b.V;xV=x.V;   /* 2/2/97 */
  S=(double *)calloc(x.r+1,sizeof(double));
  if ((a.r!=x.r)||(b.r!=x.r))
  { ErrorMessage("Incompatible matrices in hyman_filter",1);}
  for (i=0;i<a.r-1;i++) S[i+1]=(aV[i+1]-aV[i])/(xV[i+1]-xV[i]);
  S[0]=S[1];S[x.r]=S[x.r-1];
  for (i=0;i<a.r;i++)
  { if ((S[i]*S[i+1]>0.0)) sig=S[i+1];
    else sig=bV[i];
    if (sig>=0.0)
    bV[i]=min(max(0.0,bV[i]),3.0*min(fabs(S[i]),fabs(S[i+1])));
    else
    bV[i]=max(min(0.0,bV[i]),-3.0*min(fabs(S[i]),fabs(S[i+1])));
  }
#ifdef SUNC
  cfree(S);
#else
  free(S);
#endif
}


void grad_limit(x,y,g,lim) matrix x,y,g;double lim;

/* This routine was designed for use with the 'mort.c' program and attempts
   to make sure that the gradient of a piecewise polynomial is always less
   than or equal to the value given by 'lim'. Algorithm is quite cunning and
   should always work. 1/97 */

{ long i;
  /* first subtract x*lim from the data and lim from the gradients -
     then Hyman filter the results for co-monotonicity, before adding the
     corrections back on again */
  for (i=0;i<x.r;i++)
  { g.V[i] -= lim; /* subtract max rate of increase */
    y.V[i] -= lim*(x.V[i]-x.V[0]); /* subtract from values */
  }
  /* data is now monotonically decreasing - assuming it is correct */
  for (i=0;i<x.r-1;i++)
  { if (y.V[i+1]>y.V[i])
    { ErrorMessage("Faulty maturation at age data",1);}
  }
  /* now apply Hyman filter */
  hyman_filter(x,y,g);
  /* ... and add on line previously subtracted */
  for (i=0;i<x.r;i++)
  { g.V[i] += lim; /* add max rate of increase */
    y.V[i] += lim*(x.V[i]-x.V[0]); /* add to values */
  }
}

/*void splcoeffs(x,a,b,c,d) matrix x,a,b,c,d;*/

/* returns coefficients for
				2          3
   s(x)=a  + b (x-x ) + c (x-x ) + d (x-x ) .
	 i    i    i     i    i     i    i

   given the vectors x and a. Assumes that a spline with 1 knot is a constant
   and that 2 knots implies a straight line */
/*
{ matrix l,h,al,mu,z;
  long i,n;
  n=x.r;
  if (n==1L) { b.V[0]=c.V[0]=d.V[0]=0.0;return ;}
  if (n==2L)
  { b.V[1]=b.V[0]=(a.V[1]-a.V[0])/(x.V[1]-x.V[0]);c.V[0]=d.V[0]=0.0;return ;}
  h=initmat(n-1L,1L);al=initmat(n,1L);
  z=initmat(n,1L);mu=initmat(n,1L);
  l=initmat(n,1L);
  for (i=0;i<n-1;i++) h.V[i]=x.V[i+1]-x.V[i];
  for (i=1;i<n-1;i++)
  al.V[i]=3.0*(a.V[i+1]*h.V[i-1]-a.V[i]*(x.V[i+1]-x.V[i-1])+
  a.V[i-1]*h.V[i])/(h.V[i-1]*h.V[i]);
  l.V[0]=1;mu.V[0]=0.0;z.V[0]=0.0;
  for (i=1;i<n-1;i++)
  { l.V[i]=2.0*(x.V[i+1]-x.V[i-1])-h.V[i-1]*mu.V[i-1];
    mu.V[i]=h.V[i]/l.V[i];
    z.V[i]=(al.V[i]-h.V[i-1]*z.V[i-1])/l.V[i];
  }
  l.V[n-1]=1.0;z.V[n-1]=0.0;c.V[n-1]=0.0;
  for (i=n-2;i>=0;i--)
  { c.V[i]=z.V[i]-mu.V[i]*c.V[i+1];
    b.V[i]=(a.V[i+1]-a.V[i])/h.V[i]-
            h.V[i]*(c.V[i+1]+2.0*c.V[i])/3.0;
    d.V[i]=(c.V[i+1]-c.V[i])/(3.0*h.V[i]);
  }
  freemat(z);freemat(mu);freemat(al);freemat(h);freemat(l);
} */
void splcoeffs(x,a,b,c,d) matrix x,a,b,c,d;

/* returns coefficients for
				2          3
   s(x)=a  + b (x-x ) + c (x-x ) + d (x-x ) .
	 i    i    i     i    i     i    i

   given the vectors x and a. Assumes that a spline with 1 knot is a constant
   and that 2 knots implies a straight line. This version allocates no extra
   memory on the heap. */

{ long i,n;
  double h0,h1,li,*xV,*aV,*bV,*cV,*dV;
  n=x.r;
  xV=x.V;aV=a.V;bV=b.V;cV=c.V;dV=d.V;
  if (n==1L) { bV[0]=cV[0]=dV[0]=0.0;return ;}
  if (n==2L)
  { bV[1]=bV[0]=(aV[1]-aV[0])/(xV[1]-xV[0]);cV[0]=dV[0]=0.0;return ;}
  h1=xV[1]-xV[0];
  for (i=1;i<n-1;i++)
  { h0=h1;h1=xV[i+1]-xV[i];
    cV[i]=3.0*(aV[i+1]*h0-aV[i]*(xV[i+1]-xV[i-1])+
    aV[i-1]*h1)/(h0*h1);
  }
  dV[0]=0.0;bV[0]=0.0;h1=xV[1]-xV[0];
  for (i=1;i<n-1;i++)
  { h0=h1;h1=xV[i+1]-xV[i];
    li=2.0*(xV[i+1]-xV[i-1])-h0*dV[i-1];
    dV[i]=h1/li;
    bV[i]=(cV[i]-h0*bV[i-1])/li;
  }
  bV[n-1]=0.0;cV[n-1]=0.0;
  for (i=n-2;i>=0;i--)
  { h1=xV[i+1]-xV[i];
    cV[i]=bV[i]-dV[i]*cV[i+1];
    bV[i]=(aV[i+1]-aV[i])/h1-
            h1*(cV[i+1]+2.0*cV[i])/3.0;
    dV[i]=(cV[i+1]-cV[i])/(3.0*h1);
  }
}




double spl(no,x,y,x0,x1,mode,resetx,resety)
int no,mode,resetx,resety;matrix x,y;double x0,x1;

/* this is a routine that stores a number of spline functions for evaluation
   no is the index number for the splines coefficients.
   x and y are the x-y values interpolated, the x & y values need only be
   supplied when resetx/y is non-zero, they are stored between times;
   mode determines what sort of evaluation.
    0   for integration from x0 to x1.
    1   for evaluation at x0
    2   for the derivative at x0.
   resetx/y is set to 1 to reset the coefficients for spline 'no'.
   To clear space resetx=1 resety=1 and x.r= 0.
*/

{ static matrix a[125],b[125],c[125],d[125],xn[125];
  static int on[125];
  static long ipos[125],jpos[125];
  long i;
  double res,xp;
  if ((no>124)||(no<0)) { ErrorMessage("Silly index passed to spl.",1);}
  if ((resetx)||(resety))
  { if (((resetx)&&(y.r!=x.r))||((!resetx)&&(y.r!=xn[no].r)))
    { ErrorMessage("x and y vectors passed to spl() are not the same length.",1);
    }
    if (on[no])
    { if (resetx) freemat(xn[no]);
      freemat(b[no]);freemat(c[no]);freemat(d[no]);
      if (resety) freemat(a[no]);
    }
    if (x.r==0 && resetx==1) { on[no]=0; return(0.0);}
    if (resetx) xn[no]=initmat(x.r,1L);b[no]=initmat(y.r,1L);
    c[no]=initmat(b[no].r,1L);d[no]=initmat(b[no].r,1L);
    if (resety) a[no]=initmat(y.r,1L);
    if (resetx) for (i=0;i<x.r;i++) xn[no].V[i]=x.V[i];
    if (resety) for (i=0;i<y.r;i++) a[no].V[i]=y.V[i];
    splcoeffs(xn[no],a[no],b[no],c[no],d[no]);ipos[no]=0L;jpos[no]=0L;
    on[no]=1;
  } else
  { if (a[no].r!=xn[no].r)
    { ErrorMessage("y vector wrong length in spl().",1);}
  }
  if (xn[no].r==1L)
  { if (!mode) return(a[no].V[0]*(x1-x0));else
    if (mode==1L) return(a[no].V[0]);else
    if (mode==2L) return(0.0);
  } else
  if (xn[no].r==2L)
  { if (!mode) return((a[no].V[0]+b[no].V[0]*((x1+x0)*0.5-xn[no].V[0]))*(x1-x0));else
    if (mode==1L) return(a[no].V[0]+b[no].V[0]*(x0-xn[no].V[0]));else
    if (mode==2L) return(b[no].V[0]);
  }
/*  if ((x0<xn[no].V[0])||(x0>xn[no].V[xn[no].r-1]))
  ErrorMessage("x0 out of bound in spl().");
  if ((!mode)&&((x1<xn[no].V[0])||(x1>xn[no].V[xn[no].r-1])))
  ErrorMessage("x1 out of bound in spl().");*/
  x0=max(x0,xn[no].V[0]);x0=min(x0,xn[no].V[xn[no].r-1]);
  x1=max(x1,xn[no].V[0]);x1=min(x1,xn[no].V[xn[no].r-1]);
  while ((xn[no].V[ipos[no]]>x0)&&(ipos[no]>0)) ipos[no]--;
  while ((xn[no].V[ipos[no]+1]<x0)&&(ipos[no]<xn[no].r-2)) ipos[no]++;
  if (!mode)
  { while ((xn[no].V[jpos[no]]>x1)&&(jpos[no]>0)) jpos[no]--;
    while ((xn[no].V[jpos[no]+1]<x1)&&(jpos[no]<xn[no].r-2)) jpos[no]++;
    if (jpos[no]==ipos[no])
    { res=a[no].V[ipos[no]]*intb0(xn[no].V[ipos[no]],xn[no].V[ipos[no]+1],x0,x1)+
	   a[no].V[ipos[no]+1]*intb1(xn[no].V[ipos[no]],xn[no].V[ipos[no]+1],x0,x1)+
	   b[no].V[ipos[no]]*intd0(xn[no].V[ipos[no]],xn[no].V[ipos[no]+1],x0,x1)+
	   b[no].V[ipos[no]+1]*intd1(xn[no].V[ipos[no]],xn[no].V[ipos[no]+1],x0,x1);
    } else
    { res=a[no].V[ipos[no]]*intb0(xn[no].V[ipos[no]],xn[no].V[ipos[no]+1],x0,xn[no].V[ipos[no]+1])+
	   a[no].V[ipos[no]+1]*intb1(xn[no].V[ipos[no]],xn[no].V[ipos[no]+1],x0,xn[no].V[ipos[no]+1])+
	   b[no].V[ipos[no]]*intd0(xn[no].V[ipos[no]],xn[no].V[ipos[no]+1],x0,xn[no].V[ipos[no]+1])+
	   b[no].V[ipos[no]+1]*intd1(xn[no].V[ipos[no]],xn[no].V[ipos[no]+1],x0,xn[no].V[ipos[no]+1]);
      for (i=ipos[no]+1;i<jpos[no];i++)
      { res+=a[no].V[i]*intb0(xn[no].V[i],xn[no].V[i+1],xn[no].V[i],xn[no].V[i+1])+
	     a[no].V[i+1]*intb1(xn[no].V[i],xn[no].V[i+1],xn[no].V[i],xn[no].V[i+1])+
	     b[no].V[i]*intd0(xn[no].V[i],xn[no].V[i+1],xn[no].V[i],xn[no].V[i+1])+
	     b[no].V[i+1]*intd1(xn[no].V[i],xn[no].V[i+1],xn[no].V[i],xn[no].V[i+1]);
      }
      i=jpos[no];
      res+=a[no].V[i]*intb0(xn[no].V[i],xn[no].V[i+1],xn[no].V[i],x1)+
	   a[no].V[i+1]*intb1(xn[no].V[i],xn[no].V[i+1],xn[no].V[i],x1)+
	   b[no].V[i]*intd0(xn[no].V[i],xn[no].V[i+1],xn[no].V[i],x1)+
	   b[no].V[i+1]*intd1(xn[no].V[i],xn[no].V[i+1],xn[no].V[i],x1);
    }
  } else if (mode==1)
  { i=ipos[no];
    xp=x0-xn[no].V[i];
    res=a[no].V[i]+(b[no].V[i]+(c[no].V[i]+d[no].V[i]*xp)*xp)*xp;
  } else if (mode==2)
  { i=ipos[no];
    xp=x0-xn[no].V[i];
    res=b[no].V[i]+(c[no].V[i]*2.0+d[no].V[i]*3.0*xp)*xp;
  }
  return(res);
}




matrix quadcon(n,type,rescale) long n; int rescale,type;

/* produces spline like constraint matrices for measuring smoothness of
   variation of n equally spaced datapoints */

{ matrix S;
  long i;
  double mult;
  if (rescale) mult=(double)n*(double)n;else mult=1.0;
  S=initmat(n,n);
  if (type==1)
  { for (i=0;i<n;i++) S.M[0][i]=mult;
  } else
  if (type==2)
  { for (i=0;i<n-1;i++)
    { S.M[i][i]=2.0*mult;S.M[i][i+1]=S.M[i+1][i]= -mult;}
    S.M[0][0]=S.M[n-1][n-1]=mult;
  } else
  if (type==3)
  { for (i=0;i<n-2;i++)
    { S.M[i+1][i+1]=6.0*mult;
      S.M[i+2][i+1]=S.M[i+1][i+2]= -4.0*mult;
      S.M[i][i+2]=S.M[i+2][i]=mult;
    }
    S.M[0][1]=S.M[1][0]=S.M[n-1][n-2]=S.M[n-2][n-1]= -2.0*mult;
    S.M[0][0]=S.M[n-1][n-1]=mult;S.M[1][1]=S.M[n-2][n-2]=5.0*mult;
  }
  return(S);
}

void GetInf(x,W,A,lam) matrix x,W,A;double lam;


/* Obtains the influence matrix for a 1-d smoothing spline to data at
   x, given smoothing parameter lam */

{ matrix S;
  long i,j;
  getSmooth(&S,x,0);
  for (i=0;i<A.r;i++)
  { for (j=0;j<A.c;j++)
    A.M[i][j]=lam*S.M[i][j];
    A.M[i][i]+=W.V[i];
  }
  choleski(A,S,1,1);
  for (i=0;i<A.r;i++) for (j=0;j<A.c;j++) A.M[i][j]*=W.V[i];
  freemat(S);
}


/***************************************************************************/
/** Routines to fit cubic smoothing spline to data with error variance    **/
/** of the data unknown, by the method of Generalised Cross Validation    **/
/** introduced by Craven and Wahba (1979). The Spline is parametrised in  **/
/** the manner of Reinsch (1967) and calculated using the method of       **/
/** orthogonal transformations (Givens rotations) developed by Hutchinson **/
/** and deHoog (1987).                                                    **/
/***************************************************************************/


void ss_setupstore(n,ub,lb,U,V,a,b,c,d,TrA,Wy,h,ci)
long n;
double ***ub,***lb,**a,**b,**c,**d,**TrA,**Wy,**h,**ci;
Ru **U;Rv **V;

{ long i;
  *ub=(double **)calloc((size_t)n+1,sizeof(double *));
  *lb=(double **)calloc((size_t)n+1,sizeof(double *));
  for (i=0;i<=n;i++)
  { (*lb)[i]=(double *)calloc((size_t)3,sizeof(double));
    (*ub)[i]=(double *)calloc((size_t)4,sizeof(double));
  }
  *U=(Ru *)calloc((size_t)n+1,sizeof(Ru));
  *V=(Rv *)calloc((size_t)n+1,sizeof(Rv));
  *a=(double *)calloc((size_t)n+1,sizeof(double));
  *b=(double *)calloc((size_t)n+1,sizeof(double));
  *c=(double *)calloc((size_t)n+1,sizeof(double));
  *d=(double *)calloc((size_t)n+1,sizeof(double));
  *TrA=(double *)calloc((size_t)n+1,sizeof(double));
  *Wy=(double *)calloc((size_t)2*n+2,sizeof(double));
  *h=(double *)calloc((size_t)n+1,sizeof(double));
  *ci=(double *)calloc((size_t)n+1,sizeof(double));
  if (!(*ci))
  { ErrorMessage("Out of memory in ss_setupstore. Exiting ....",1);}
}

void ss_freestore(n,ub,lb,U,V,TrA,Wy,h)
long n;
double ***ub,***lb,**TrA,**Wy,**h;
Ru **U;Rv **V;

{ long i;
  for (i=0;i<=n;i++)
  { free((*lb)[i]);free((*ub)[i]);
  }
  free(*ub);free(*lb);free(*U);free(*V);
  free(*TrA);free(*Wy);free(*h);
}



/****************************************************************************/
/**  QTz performs the Givens rotation specified by i,j,c and s on the      **/
/**  vector Wy which originally holds vector z (3.6) of H&dH and finally   **/
/**  holds the residual vector.                                            **/
/****************************************************************************/

void QTz(i,j,c,s,Wy)long i,j;double c,s,*Wy;

{ double temp;
  temp=Wy[i]*c+Wy[j]*s;
  Wy[j]=Wy[j]*c-Wy[i]*s;
  Wy[i]=temp;
}

/***************************************************************************/
/** givens calculates the parameters s and c for the Givens rotation that **/
/** will annihilate b.                                                    **/
/***************************************************************************/

void givens(a,b,c,s)double a,b,*c,*s;

/* Calculates c and s to annihilate b */

{ double t;
  if (a==0.0) { *c=1.0;*s=0.0;} else
  if (fabs(a)<=fabs(b))
  { t=a/b;
    *s=1/sqrt(1+t*t);
    *c= (*s)*t;
  } else
  { t=b/a;
    *c=1/sqrt(1+t*t);
    *s=(*c)*t;
  }
}


/****************************************************************************/
/** GCV calculates the generalised cross validation function derived by    **/
/** Craven and Wahba as described by (2.14) of H&dH (1985). The method of  **/
/** calculation is the algorithm described by H&dH (1987) as an extension  **/
/** of Eldens (1984) least squares algorithm.                              **/
/****************************************************************************/

double GCV(p,ub,lb,U,V,TrA,Wy,h,ci,w,x,y,sig2,n,method)
double p,**ub,**lb,*TrA,*Wy,*h,*ci,*w,*x,*y,sig2;long n;int method;
Ru *U;Rv *V;
/* Note that this routine modifies ub and lb */


{ long i,k,j;
  double Lt,temp,c,s,residual,trace,upper,L[4][4],X[4];
  ss_setup(ub,lb,x,h,w,n);
  for (i=1;i<=n;i++) Wy[i]=y[i]/w[i];
  for (i=n+1;i<2*n+2;i++) Wy[i]=0.0;   /*NOTE trial fix 2*n+2 from 2*n-2!!*/
  p=1/sqrt(p);
/****************************************************************************/
/** multiplies the upper band of the matrix by the new smoothing parameter **/
/****************************************************************************/
  for (i=1;i<=n-2;i++) for (j=1;j<=3;j++) ub[i][j]=p*ub[i][j];
/****************************************************************************/
/** Performs the QR decomposition of the two banded matrix consisting of   **/
/** ub and lb storing the Givens rotations which make up Q in the record   **/
/** arrays V and U where the notation and indeces are defined by (3.13) -  **/
/** (3.15) of H&dH (1987).                                                 **/
/****************************************************************************/
  for (i=1;i<=n-3;i++)
  { givens(ub[i+1][1],lb[i][2],&c,&s);
    QTz(i+1,n+i,c,s,Wy);
    temp=c*lb[i][1]-s*ub[i][2];
    ub[i+1][1]=s*lb[i][2]+c*ub[i+1][1];
    ub[i][2]=s*lb[i][1]+c*ub[i][2];
    lb[i][1]=temp;
    U[i].Rni[1][0]=-s;U[i].Rni[1][1]=c;

    givens(ub[i][1],lb[i][1],&c,&s);
    QTz(i,n+i,c,s,Wy);
    ub[i][1]=c*ub[i][1]+s*lb[i][1];
    U[i].Rni[0][0]=-s;U[i].Rni[0][1]=c;

    givens(ub[i][1],ub[i][2],&c,&s);
    QTz(i,i+1,c,s,Wy);
    ub[i][1]=c*ub[i][1]+s*ub[i][2];
    upper=s*ub[i+1][1];
    ub[i+1][1]=c*ub[i+1][1];
    V[i].Ri[1][0]=-s;V[i].Ri[1][1]=c;

    givens(ub[i][1],ub[i][3],&c,&s);
    QTz(i,i+2,c,s,Wy);
    ub[i+1][2]=c*ub[i+1][2]-s*upper;
    if (i!=(n-3)) ub[i+2][1]=c*ub[i+2][1];
    V[i].Ri[2][0]=-s;V[i].Ri[2][1]=c;
    /* check that this should not be in if */
  }
  i=n-2;
  givens(ub[i][1],lb[i][1],&c,&s);
  QTz(i,n+i,c,s,Wy);
  ub[i][1]=c*ub[i][1]+s*lb[i][1];
  U[i].Rni[0][0]=-s;U[i].Rni[0][1]=c;

  givens(ub[i][1],ub[i][2],&c,&s);
  QTz(i,i+1,c,s,Wy);
  ub[i][1]=c*ub[i][1]+s*ub[i][2];
  V[i].Ri[1][0]=-s;V[i].Ri[1][1]=c;

  givens(ub[i][1],ub[i][3],&c,&s);
  QTz(i,i+2,c,s,Wy);
  V[i].Ri[2][0]=-s;V[i].Ri[2][1]=c;

  /** Calculates Weighted Residual **/

  for (i=(n-1);i<(2*n+2);i++) Wy[i]=0.0;
  for (i=(n-2);i>=1;i--)
  { QTz(i,i+2,V[i].Ri[2][1],V[i].Ri[2][0],Wy);
    QTz(i,i+1,V[i].Ri[1][1],V[i].Ri[1][0],Wy);
    QTz(i,n+i,U[i].Rni[0][1],U[i].Rni[0][0],Wy);
    if (i!=(n-2)) QTz(i+1,n+i,U[i].Rni[1][1],U[i].Rni[1][0],Wy);
  }
  residual=0.0;
  for (i=1;i<=n;i++) residual+=Wy[i]*Wy[i];

  /** Calculate the Trace of the influence matrix **/

  L[3][1]=V[n-2].Ri[2][1];L[3][3]=-V[n-2].Ri[2][0];
  L[3][2]=-L[3][1]*V[n-2].Ri[1][0];L[3][1]=L[3][1]*V[n-2].Ri[1][1];

  X[3]=-L[3][1]*U[n-2].Rni[0][0];
  L[3][1]=L[3][1]*U[n-2].Rni[0][1];

  TrA[n]=L[3][3]*L[3][3];
  L[3][3]=L[3][2];L[3][2]=L[3][1];

  L[2][1]=V[n-3].Ri[2][1];L[2][3]=-V[n-3].Ri[2][0];
  L[3][1]=L[3][3]*V[n-3].Ri[2][0];
  L[3][3]=L[3][3]*V[n-3].Ri[2][1];
  L[2][2]=-L[2][1]*V[n-3].Ri[1][0];
  L[2][1]=L[2][1]*V[n-3].Ri[1][1];
  Lt=L[3][1]*V[n-3].Ri[1][1]+L[3][2]*V[n-3].Ri[1][0];
  L[3][2]=L[3][2]*V[n-3].Ri[1][1]-L[3][1]*V[n-3].Ri[1][0];
  L[3][1]=Lt;

  X[2]=-L[2][1]*U[n-3].Rni[0][0];L[2][1]=L[2][1]*U[n-3].Rni[0][1];
  X[3]=-L[3][1]*U[n-3].Rni[0][0];L[3][1]=L[3][1]*U[n-3].Rni[0][1];
  L[2][2]=L[2][2]*U[n-3].Rni[1][1]+X[2]*U[n-3].Rni[1][0];
  L[3][2]=L[3][2]*U[n-3].Rni[1][1]+X[3]*U[n-3].Rni[1][0];

  TrA[n-1]=(L[3][3]*L[3][3]+L[2][3]*L[2][3]);
  givens(L[2][1],L[3][1],&c,&s);        /** The succesive rotation  **/
  L[2][1]=L[2][1]*c+L[3][1]*s;
  Lt=L[2][2]*c+L[3][2]*s;
  L[3][2]=L[3][2]*c-L[2][2]*s;L[2][2]=Lt;
  L[3][3]=L[3][2];L[2][3]=L[2][2];L[2][2]=L[2][1];

  for (i=(n-4);i>=1;i--)
  { L[1][3]=-V[i].Ri[2][0];L[1][1]=V[i].Ri[2][1];
    L[2][1]=L[2][3]*V[i].Ri[2][0];L[2][3]=L[2][3]*V[i].Ri[2][1];
    L[3][1]=L[3][3]*V[i].Ri[2][0];L[3][3]=L[3][3]*V[i].Ri[2][1];
    givens(L[1][1],L[3][1],&c,&s);s=-s; /** Rotation to remove upper
		  element BEFORE it propagates **/
    L[1][1]=L[1][1]*c-L[3][1]*s;
    L[1][2]=-L[1][1]*V[i].Ri[1][0];
    L[1][1]=L[1][1]*V[i].Ri[1][1];
    Lt=L[2][1]*V[i].Ri[1][1]+L[2][2]*V[i].Ri[1][0];
    L[2][2]=L[2][2]*V[i].Ri[1][1]-L[2][1]*V[i].Ri[1][0];
    L[2][1]=Lt;

    for (k=1;k<=2;k++)
    { X[k]=-L[k][1]*U[i].Rni[0][0];
      L[k][1]=L[k][1]*U[i].Rni[0][1];
      L[k][2]=L[k][2]*U[i].Rni[1][1]+X[k]*U[i].Rni[1][0];
    }
    givens(L[1][1],L[2][1],&c,&s);      /** Second rotation removing
	     upper element **/
    L[1][1]=L[1][1]*c+L[2][1]*s;
    Lt=L[1][2]*c+L[2][2]*s;
    L[2][2]=L[2][2]*c-L[1][2]*s;L[1][2]=Lt;
    TrA[i+2]=L[3][3]*L[3][3]+L[2][3]*L[2][3]+L[1][3]*L[1][3];
    if (i!=1)
    { L[3][3]=L[2][2];L[2][3]=L[1][2];L[2][2]=L[1][1];
    }
  }
  TrA[2]=L[2][2]*L[2][2]+L[1][2]*L[1][2];
  TrA[1]=L[1][1]*L[1][1];
  trace=0.0;
  for (i=1;i<=n;i++) trace+=TrA[i];
  for (i=1;i<=n;i++)   ci[i]=1.96*sqrt(residual/trace)*sqrt(1.0-TrA[i]);
  if (method==0)
  { residual=0.0;
    for (i=1;i<=n;i++) residual+=Wy[i]*Wy[i]/(TrA[i]*TrA[i]);
    return(residual/n);               /* ordinary cross validation */
  } else
  if (method==1)
  return(residual+2*sig2*(n-trace));  /* expected loss function */
  else
  return(n*residual/(trace*trace));   /* generalized cross validation */
}

/****************************************************************************/
/** The following routine  performs a global search for the minimiser of   **/
/** GCV, OCV or optimal loss score starting at low p                       **/
/** i.e. very smooth this is necessary because at very high p with very    **/
/** noisy data the gcv function seems to tend to zero. Once the approximate**/
/** minimum is found a bisection search is performed, for the best p.      **/
/** Failure of the algorithm to converge should be countered in the first  **/
/** instance by lowering the starting value of p.                          **/
/****************************************************************************/

double ss_minimiser(ub,lb,U,V,TrA,Wy,h,ci,w,x,y,sig2,n,method)
double **ub,**lb,*TrA,*Wy,*h,*ci,*w,*x,*y,sig2;
long n;int method;Ru *U;Rv *V;

{ double p0,p1,p2,cv1,cv2=0.0,max=0.0,pt,p1t,ft,f1t,tau;
  long i;
  char down=0;
  p1=1e-7;
  for (i=0;i<82L;i++)
  { cv1=GCV(p1,ub,lb,U,V,TrA,Wy,h,ci,w,x,y,sig2,n,method);
    if (cv1>max) max=cv1;
    if (cv1<max*0.999) down=1;
    if ((!i)||(cv1<cv2))
    { p2=p1;cv2=cv1;}
    if ((down)&&(cv1>cv2*1.01)&&(p1>20*p2)) i=82L;
    p1*=1.4;
  }
  /* golden section search to polish minimisation */
  p0=p2/1.4;p1=p2*1.4;
  tau=2.0/(1.0+sqrt(5.0));
  pt=p0+(p1-p0)*tau;
  ft=GCV(pt,ub,lb,U,V,TrA,Wy,h,ci,w,x,y,sig2,n,method);
  p1t=p0+(p1-p0)*(1.0-tau);
  f1t=GCV(p1t,ub,lb,U,V,TrA,Wy,h,ci,w,x,y,sig2,n,method);
  while ((pt-p1t)>1e-5*fabs(pt+p1t))
  { if (ft<f1t)
    { p0=p1t;p1t=pt;f1t=ft;pt=p0+(p1-p0)*tau;
      ft=GCV(pt,ub,lb,U,V,TrA,Wy,h,ci,w,x,y,sig2,n,method);
    } else
    { p1=pt;pt=p1t;ft=f1t;p1t=p0+(p1-p0)*(1.0-tau);
      f1t=GCV(p1t,ub,lb,U,V,TrA,Wy,h,ci,w,x,y,sig2,n,method);
    }
  }
  p2=pt;
  GCV(p2,ub,lb,U,V,TrA,Wy,h,ci,w,x,y,sig2,n,method);
  return(p2);
 
}


/****************************************************************************/
/** The spline coefficients as described in Reinsch 1967 are calculated    **/
/** from the residual vector as described in H&dH (1987) (3.8), (3.10) and **/
/** Reinsch.                                                               **/
/****************************************************************************/

void ss_coeffs(lb,a,b,c,d,h,y,w,Wy,n)
double **lb,*a,*b,*c,*d,*h,*y,*w,*Wy;long n;

{ double *GTA,*z;
  long i;
  GTA=(double *)calloc((size_t)n+1,sizeof(double));
  z=(double *)calloc((size_t)n+1,sizeof(double));
  if (!z)
  { ErrorMessage("Out of memory in ss_coeffs. Exiting .... ",1);}
  for (i=1;i<=n;i++) a[i]=y[i]-w[i]*Wy[i];
  for (i=1;i<=n-2;i++)
  GTA[i]=a[i]/h[i]-a[i+1]*(1/h[i]+1/h[i+1])+a[i+2]/h[i+1];
  z[1]=GTA[1]/lb[1][1];
  for (i=2;i<=n-2;i++) z[i]=(GTA[i]-lb[i-1][2]*z[i-1])/lb[i][1];
  c[n-1]=z[n-2]/lb[n-2][1];c[n]=0.0;c[1]=0.0;
  for (i=n-3;i>=1;i--)
  c[i+1]=(z[i]-lb[i][2]*c[i+2])/lb[i][1];
  for (i=1;i<=n-1;i++)
  { d[i]=(c[i+1]-c[i])/(3*h[i]);
    b[i]=(a[i+1]-a[i])/h[i]-c[i]*h[i]-d[i]*h[i]*h[i];
  }
  free(GTA);free(z);
}

/****************************************************************************/
/**  The upper and lower bands ,ub and lb of the matrix C of H&dH (3.2) are**/
/** set up the lower band coming from the Choleski decomposition of H (3.1)**/
/****************************************************************************/

void ss_setup(ub,lb,x,h,w,n)
double **ub,**lb, *x,*h,*w;long n;


{ double *hh,*hh1;
  long i;
  hh=(double *)calloc((size_t)n+1,sizeof(double));
  hh1=(double *)calloc((size_t)n+1,sizeof(double));
  if (!hh1) { ErrorMessage("Out of memory in ss_setup. Exiting .... ",1);}
  for (i=1;i<=n-1;i++) h[i]=x[i+1]-x[i];

  for (i=1;i<=n-2;i++) hh[i]=2.0*(h[i]+h[i+1])/3.0;
  for (i=1;i<=n-3;i++) hh1[i]=h[i+1]/3.0;
  lb[1][1]=sqrt(hh[1]);lb[1][2]=hh1[1]/lb[1][1];
  for (i=2;i<=n-3;i++)
  { lb[i][1]= sqrt(hh[i]-lb[i-1][2]*lb[i-1][2]);
    lb[i][2]= hh1[i]/lb[i][1];
  }
  lb[n-2][1]=sqrt(hh[n-2]-lb[n-3][2]*lb[n-3][2]);
  for (i=1;i<=n-2;i++)
  { ub[i][1]=w[i]/h[i];
    ub[i][2]=-w[i+1]*(1/h[i]+1/h[i+1]);
    ub[i][3]=w[i+2]/h[i+1];
  }
  free(hh);free(hh1);
}

/****************************************************************************/
/** The spline function.                                                   **/
/****************************************************************************/

double polyspl(xx,a,b,c,d,x,i) double xx,*a,*b,*c,*d,*x;long i;

{ double xxi;
  xxi=xx-x[i];
  return(a[i]+b[i]*xxi+c[i]*xxi*xxi+d[i]*xxi*xxi*xxi);
}


double fitss_array(lam,n,x,y,w,a,b,c,d,ci,method)
long n; double lam,*x,*y,*w,**a,**b,**c,**d,**ci;int method;

/* fits cubic smoothing spline to data in x,y using (relative)
   variance estimates in w. Note that this was translated from pascal
   code - and all arrays start from 1 as a result.
   a,b,c,d and ci are initialised in this routine.
*/

{ double **ub,**lb;
  Ru *U;
  Rv *V;
  double *TrA,*Wy,*h,sig2=0.0;
  long i;
  for (i=1;i<=n;i++) sig2+=w[i];sig2/=n;
  for (i=1;i<=n;i++) w[i]=sqrt(w[i]/sig2);
  /* algorithm uses strange convention that if the model is
     Yi=f(xi)+ei where ei is an error term then E(ei*ei)=w[i]*w[i]*var
     where var is the error variance */
  ss_setupstore(n,&ub,&lb,&U,&V,a,b,c,d,&TrA,&Wy,&h,ci);
  ss_setup(ub,lb,x,h,w,n);
  if (method==3) { GCV(lam,ub,lb,U,V,TrA,Wy,h,*ci,w,x,y,sig2,n,1);}
  else lam=ss_minimiser(ub,lb,U,V,TrA,Wy,h,*ci,w,x,y,sig2,n,method);
  ss_setup(ub,lb,x,h,w,n); /* needed to get back original lb for .... */
  ss_coeffs(lb,*a,*b,*c,*d,h,y,w,Wy,n);
  for (i=1;i<=n;i++) w[i]=w[i]*w[i]*sig2;   /* restoring w to original values */
  ss_freestore(n,&ub,&lb,&U,&V,&TrA,&Wy,&h);
  return(lam);
}


double ss_fit(x,y,w,a,b,c,d,method,justa) matrix x,y,w,a,b,c,d;int method,justa;
/* Routine for calculating cubic smoothing splines by Hutchinson & deHoog's
   method. This is really an interface routine to fitss_array, which works
   with arrays rather than matrices. Set justa=1 if only the fitted
   values are required, otherwise set justa=0.
   For method: 0 - OCV
	       1 - Optimal loss function (w must contain absolute variances)
	       2 - GCV
          3 - use s.p. in a.V[0]
   All matrices must be initialised prior to calling (except b,c,d if justa=1)
   returns smoothing parameter.
*/

{ double *ar,*br,*cr,*dr,*ci,lam;
  long i;
  x.V--;y.V--;w.V--; /* array version requires arrays starting at 1 */
  lam=a.V[0];   // only used if method==3
  lam =fitss_array(lam,x.r,x.V,y.V,w.V,&ar,&br,&cr,&dr,&ci,method);
  x.V++;y.V++;w.V++; /* resetting addresses to correct values */
  if (justa)
  for (i=0;i<x.r;i++) a.V[i]=ar[i+1];
  else
  for (i=0;i<x.r;i++)
  { a.V[i]=ar[i+1];
    b.V[i]=br[i+1];
    c.V[i]=cr[i+1];d.V[i]=dr[i+1];
  }
  free(ar);free(br);free(cr);free(dr);free(ci);
  return(lam);
}

/************************************************************************************/
/* NOTES..................                                                          */
/************************************************************************************/

/* 
9/5/2000 spl() modified so that freeing the memory by a call with resetx==1, resety==1 
         and x.r==0, results in the spline being recorded as uninitialised (which it 
         needs to be to avoid attempts to free random chunks of memory later!)
*/