ADMB Documentation  11.1.2402
 All Classes Files Functions Variables Typedefs Friends Defines
vcubicspline.cpp
Go to the documentation of this file.
00001 #include "statsLib.h"
00002 
00014 vcubic_spline_function_array::vcubic_spline_function_array( const dvector & x)
00015  {
00016    int mmin=indexmin();
00017    int mmax=indexmax();
00018 
00019    int n=mmax-mmin+1;
00020    ptr = new pvcubic_spline_function[n];
00021    ptr-=mmin;
00022    for (int i=mmin;i<=mmax;i++)
00023    {
00024      // ptr[i]= new  vcubic_spline_function(x); // not sure how to call this...should return matrix ...
00025    }
00026  }
00027 
00028 vcubic_spline_function_array::vcubic_spline_function_array(int mmin,int mmax,
00029  const dvector & x, const dvar_matrix& y)
00030  {
00031    index_min=mmin;
00032    index_max=mmax;
00033    int n=mmax-mmin+1;
00034    ptr = new pvcubic_spline_function[n];
00035    ptr-=mmin;
00036    for (int i=mmin;i<=mmax;i++)
00037    {
00038      ptr[i]= new  vcubic_spline_function(x,y(i));
00039    }
00040  }
00041 
00042 vcubic_spline_function & vcubic_spline_function_array::operator () (int i)
00043 {
00044   if (i<indexmin() || i> indexmax())
00045    {
00046      cerr << "index out of range in function"
00047           " vcubic_spline_function & operator () (int i)"
00048       << endl;
00049      ad_exit(1);
00050    }
00051    return *(ptr[i]);
00052  }
00053 
00054 
00055 dvar_matrix vcubic_spline_function_array::operator( ) (const dvector & v)
00056  {
00057    int mmin=indexmin();
00058    int mmax=indexmax();
00059    dvar_matrix tmp(mmin,mmax,v.indexmin(),v.indexmax());
00060    for (int i=mmin;i<=mmax;i++)
00061    {
00062      tmp(i)=(*this)(i)(v);
00063    }
00064    return tmp;
00065  }
00066 
00067 vcubic_spline_function_array::~vcubic_spline_function_array()
00068 {
00069    int mmin=indexmin();
00070    int mmax=indexmax();
00071 
00072    for (int i=mmin;i<=mmax;i++)
00073    {
00074      delete ptr[i];
00075    }
00076    ptr+=mmin;
00077    delete ptr;
00078    ptr=0;
00079  }
00080 
00081 
00082 
00083 
00093 dvar_vector cubic_spline(const dvar_vector& spline_nodes, const dvector& ip)
00094 {
00095   RETURN_ARRAYS_INCREMENT();                                                              
00096   int nodes=size_count(spline_nodes);
00097   dvector ia(1,nodes);
00098   ia.fill_seqadd(0,1./(nodes-1));
00099   dvector fa = (ip-min(ip))/(max(ip)-min(ip));
00100   vcubic_spline_function ffa(ia,spline_nodes);
00101   RETURN_ARRAYS_DECREMENT();
00102   return(ffa(fa));
00103 }
00104 
00105 
00106 
00107 void bicubic_spline(const dvector& x, const dvector& y, dvar_matrix& knots, dvar_matrix& S)
00108 {
00109   /*
00110   Author:  Steven Martell
00111   Date: July 29, 2010
00112   Comments:  Based on code from Numerical Recipies.
00113 
00114   This function returns matrix S which is the interpolated values of knots
00115   over knots[1..m][1..n] grid.
00116 
00117   first call splie2 to get second-derivatives at knot points
00118   void splie2(const dvector& _x1a,const dvector& _x2a,const dmatrix& _ya,dvar_matrix& _y2a)
00119 
00120   then run the splin2 to get the spline points
00121   dvariable splin2(const dvector& _x1a,const dvector* _x2a, const dmatrix _ya,
00122     dvar_matrix& _y2a, const double& x1,const double& x2)
00123   */
00124 
00125   RETURN_ARRAYS_INCREMENT();
00126   int i,j;
00127   int m=knots.rowmax();
00128   int n=knots.colmax();
00129 
00130   int mm=S.rowmax()-S.rowmin()+1;
00131   int nn=S.colmax()-S.colmin()+1;
00132 
00133   dvar_matrix shift_S(1,mm,1,nn);
00134 
00135   dvector im(1,mm); im.fill_seqadd(0,1./(mm-1.));
00136   dvector in(1,nn); in.fill_seqadd(0,1./(nn-1.));
00137   dvar_matrix y2(1,m,1,n);  //matrix of second-derivatives
00138   y2=splie2(x,y,knots);
00139 
00140   for(i=1;i<=mm;i++){
00141     for(j=1;j<=nn;j++){
00142       shift_S(i,j)=splin2(x,y,knots,y2,in(j),im(i));
00143     }
00144   }
00145 
00146   int ii,jj;
00147   ii=0;
00148   for(i=S.rowmin();i<=S.rowmax();i++)
00149   {
00150     ii++; jj=0;
00151     for(j=S.colmin();j<=S.colmax();j++)
00152     {
00153       jj++;
00154       S(i,j)=shift_S(ii,jj);
00155     }
00156   }
00157 
00158   //cout<<shift_S<<endl;
00159   RETURN_ARRAYS_DECREMENT();
00160   //cout<<"Bicubic"<<endl;
00161 }
00162 
00163 
00164 
00165   // dvar_vector spline(const dvector &_x,const dvar_vector&_y,double yp1,double ypn)
00166   // {
00167   //  RETURN_ARRAYS_INCREMENT();
00168   //  dvector& x=(dvector&) _x;
00169   //  dvar_vector& y=(dvar_vector&) _y;
00170   //  int orig_min=x.indexmin();
00171   //  x.shift(1);
00172   //  y.shift(1);
00173   //  // need to check that x is monotone increasing;
00174   //  if  ( x.indexmax() != y.indexmax() )
00175   //  {
00176   //    cerr << " Incompatible bounds in input to spline, fix it" << endl;
00177   //  }
00178   //  int n=x.indexmax();
00179   //  dvar_vector y2(1,n);
00180   //  int i,k;
00181   //  dvariable  p,qn,sig,un;
00182   //  dvar_vector u(1,n-1);
00183   //  if (yp1 > 0.99e30)
00184   //   {
00185   //     y2[1]=u[1]=0.0;
00186   //   }
00187   //   else
00188   //   {
00189   //     y2[1] = -0.5;
00190   //     u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
00191   //   }
00192   //   for (i=2;i<=n-1;i++)
00193   //   {
00194   //     sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
00195   //     p=sig*y2[i-1]+2.0;
00196   //     y2[i]=(sig-1.0)/p;
00197   //     u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
00198   //     u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
00199   //   }
00200   //   if (ypn > 0.99e30)
00201   //   {
00202   //     qn=un=0.0;
00203   //   }
00204   //   else
00205   //   {
00206   //     qn=0.5;
00207   //     un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
00208   //   }
00209   //   y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
00210   //   for (k=n-1;k>=1;k--)
00211   //   {
00212   //     y2[k]=y2[k]*y2[k+1]+u[k];
00213   //   }
00214   //   x.shift(orig_min);
00215   //   y.shift(orig_min);
00216   //   y2.shift(orig_min);
00217   //   RETURN_ARRAYS_DECREMENT();
00218   //   return y2;
00219   // }
00220   //
00221   //
00222   // dvariable splint(const dvector& _xa,const dvar_vector& _ya, const dvar_vector& _y2a,double x)
00223   // {
00224   //   RETURN_ARRAYS_INCREMENT();
00225   //   dvector& xa=(dvector&) _xa;
00226   //   dvar_vector& ya=(dvar_vector&) _ya;
00227   //   dvar_vector& y2a=(dvar_vector&) _y2a;
00228   //   int orig_min=xa.indexmin();
00229   //   xa.shift(1);
00230   //   ya.shift(1);
00231   //   y2a.shift(1);
00232   //   dvariable y;
00233   //   int n = xa.indexmax();
00234   //   int klo,khi,k;
00235   //   dvariable h,b,a;
00236   //
00237   //   klo=1;
00238   //   khi=n;
00239   //   while (khi-klo > 1)
00240   //   {
00241   //     k=(khi+klo) >> 1;
00242   //     if (xa[k] > x) khi=k;
00243   //     else klo=k;
00244   //   }
00245   //   h=xa[khi]-xa[klo];
00246   //   a=(xa[khi]-x)/h;
00247   //   b=(x-xa[klo])/h;
00248   //   y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
00249   //   xa.shift(orig_min);
00250   //   ya.shift(orig_min);
00251   //   y2a.shift(orig_min);
00252   //   RETURN_ARRAYS_DECREMENT();
00253   //   return y;
00254   // }
00255   //
00256 
00257 dvar_matrix splie2(const dvector& _x1a,const dvector& _x2a,const dvar_matrix& _ya)//,dvar_matrix& _y2a)
00258 {
00259 /*  NR code:
00260   void splie2(float x1a[], float x2a[], float **ya, int m, int n, float **y2a)
00261   Given an m by n tabulated function ya[1..m][1..n], and tabulated independent variables
00262   x2a[1..n], this routine constructs one-dimensional natural cubic splines of the rows of ya
00263   and returns the second-derivatives in the array y2a[1..m][1..n]. (The array x1a[1..m] is
00264   included in the argument list merely for consistency with routine splin2.)
00265   {
00266   void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]);
00267   int j;
00268   for (j=1;j<=m;j++)
00269   spline(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]); Values 1x1030 signal a nat-
00270   } */
00271   RETURN_ARRAYS_INCREMENT();
00272   dvector& x1a=(dvector&) _x1a;
00273   dvector& x2a=(dvector&) _x2a;
00274   dvar_matrix& ya=(dvar_matrix&) _ya;
00275   //dvar_matrix& y2a=(dvar_matrix&) _y2a;
00276   int m=ya.rowmax();
00277   int n=ya.colmax();
00278   dvar_matrix y2a(1,m,1,n);
00279   int j;
00280   for(j=1;j<=m;j++)
00281     y2a(j)=spline(x1a,ya(j),1.0e30,1.e30);
00282   //function should return second-derivatives in y2a[1..m][1..n]
00283   RETURN_ARRAYS_DECREMENT();
00284   return y2a;
00285 }
00286 
00287 
00288 
00289 dvariable splin2(const dvector& _x1a,const dvector& _x2a, const dvar_matrix _ya,
00290   dvar_matrix& _y2a, const double& x1,const double& x2)//,dvariable& y)
00291 {
00292   /*
00293   Original NR code:
00294   void splin2(float x1a[], float x2a[], float **ya, float **y2a, int m, int n,
00295   float x1, float x2, float *y)
00296   Given x1a, x2a, ya, m, n as described in splie2 and y2a as produced by that routine; and
00297   given a desired interpolating point x1,x2; this routine returns an interpolated function value y
00298   by bicubic spline interpolation.
00299   {
00300   void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]);
00301   void splint(float xa[], float ya[], float y2a[], int n, float x, float *y);
00302   int j;
00303   float *ytmp,*yytmp;
00304   ytmp=vector(1,m);
00305   yytmp=vector(1,m); Perform m evaluations of the row splines constructed by
00306   splie2, using the one-dimensional spline evaluator
00307   splint.
00308   for (j=1;j<=m;j++)
00309   splint(x2a,ya[j],y2a[j],n,x2,&yytmp[j]);
00310   spline(x1a,yytmp,m,1.0e30,1.0e30,ytmp); Construct the one-dimensional colsplint(
00311   x1a,yytmp,ytmp,m,x1,y); umn spline and evaluate it.
00312   free_vector(yytmp,1,m);
00313   free_vector(ytmp,1,m);
00314   }
00315   */
00316   RETURN_ARRAYS_INCREMENT();
00317   dvector& x1a=(dvector&) _x1a;
00318   dvector& x2a=(dvector&) _x2a;
00319   dvar_matrix& ya=(dvar_matrix&) _ya;
00320   dvar_matrix& y2a=(dvar_matrix&) _y2a;
00321   int j;
00322   int m=ya.rowmax();
00323   int n=ya.colmax();
00324   dvariable y;
00325   dvar_vector ytmp(1,m);
00326   dvar_vector yytmp(1,m);
00327   for (j=1;j<=m;j++)
00328     yytmp[j]=splint(x1a,ya[j],y2a[j],x2);
00329 
00330   ytmp=spline(x2a,yytmp,1.0e30,1.0e30);
00331   y=splint(x2a,yytmp,ytmp,x1);
00332 
00333   RETURN_ARRAYS_DECREMENT();
00334   return(y);
00335 }
00336