ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
dfqromb.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dfqromb.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  */
00012 #include <admodel.h>
00013 //#define EPS 1.0e-4
00014 #define JMAX 50
00015 //#define JMAXP JMAX+1
00016 //#define K 5
00017 
00018 class model_parameters;
00019 
00020 dvariable trapzd(dvariable (model_parameters::*func)(const dvariable&),
00021   double a,double b,int n);
00022 dvariable trapzd(dvariable (model_parameters::*func)(const dvariable&),
00023   double a, const dvariable& b,int n);
00024 dvariable trapzd(dvariable (model_parameters::*func)(const dvariable&),
00025   const dvariable& a, const dvariable& b, int n);
00026 dvariable trapzd(dvariable (model_parameters::*func)(const dvariable&),
00027   const dvariable& a, double b, int n);
00028 
00029 void polint(const dvector& xa, const dvar_vector& ya,int n,double x,
00030   const dvariable& y, const dvariable& dy);
00031 
00039 dvariable function_minimizer::adromb(
00040   dvariable (model_parameters::*func)(const dvariable&),
00041   double a, double b, int ns)
00042 {
00043   const double base = 4;
00044   int MAXN = min(JMAX, ns);
00045   dvar_vector s(1,MAXN+1);
00046 
00047   for(int j=1; j<=MAXN+1; j++)
00048   {
00049     s[j] = trapzd(func,a,b,j);
00050   }
00051 
00052   for(int iter=1; iter<=MAXN+1; iter++)
00053   {
00054     for(int j=1; j<=MAXN+1-iter; j++)
00055     {
00056       s[j] = (pow(base,iter)*s[j+1]-s[j])/(pow(base,iter)-1);
00057     }
00058   }
00059 
00060   return s[1];
00061 }
00062 
00070 dvariable function_minimizer::adromb(
00071   dvariable (model_parameters::*func)(const dvariable&),
00072   const dvariable& a, double b, int ns)
00073 {
00074   const double base = 4;
00075   int MAXN = min(JMAX, ns);
00076   dvar_vector s(1,MAXN+1);
00077 
00078   for(int j=1; j<=MAXN+1; j++)
00079   {
00080     s[j] = trapzd(func,a,b,j);
00081   }
00082 
00083   for(int iter=1; iter<=MAXN+1; iter++)
00084   {
00085     for(int j=1; j<=MAXN+1-iter; j++)
00086     {
00087       s[j] = (pow(base,iter)*s[j+1]-s[j])/(pow(base,iter)-1);
00088     }
00089   }
00090 
00091   return s[1];
00092 }
00093 
00101 dvariable function_minimizer::adromb(
00102   dvariable (model_parameters::*func)(const dvariable&),
00103   double a, const dvariable& b, int ns)
00104 {
00105   const double base = 4;
00106   int MAXN = min(JMAX, ns);
00107   dvar_vector s(1,MAXN+1);
00108 
00109   for(int j=1; j<=MAXN+1; j++)
00110   {
00111     s[j] = trapzd(func,a,b,j);
00112   }
00113 
00114   for(int iter=1; iter<=MAXN+1; iter++)
00115   {
00116     for(int j=1; j<=MAXN+1-iter; j++)
00117     {
00118       s[j] = (pow(base,iter)*s[j+1]-s[j])/(pow(base,iter)-1);
00119     }
00120   }
00121 
00122   return s[1];
00123 }
00124 
00132 dvariable function_minimizer::adromb(
00133   dvariable (model_parameters::*func)(const dvariable&),
00134   const dvariable& a, const dvariable& b, int ns)
00135 {
00136   const double base = 4;
00137   int MAXN = min(JMAX, ns);
00138   dvar_vector s(1,MAXN+1);
00139 
00140   for(int j=1; j<=MAXN+1; j++)
00141   {
00142     s[j] = trapzd(func,a,b,j);
00143   }
00144 
00145   for(int iter=1; iter<=MAXN+1; iter++)
00146   {
00147     for(int j=1; j<=MAXN+1-iter; j++)
00148     {
00149       s[j] = (pow(base,iter)*s[j+1]-s[j])/(pow(base,iter)-1);
00150     }
00151   }
00152 
00153   return s[1];
00154 }
00155 
00163 dvariable function_minimizer::trapzd(
00164   dvariable (model_parameters::*func)(const dvariable&),
00165   double a, double b, int n)
00166 {
00167   double x,num_interval,hn;
00168   dvariable sum;
00169   static dvariable s;
00170   static int interval;
00171   int j;
00172   model_parameters * ptr= (model_parameters *) mycast();
00173   if (n == 1) {
00174     interval=1;
00175     return (s=0.5*(b-a)*((ptr->*func)(a)+(ptr->*func)(b)));
00176   } else {
00177     num_interval=interval;
00178     hn=(b-a)/num_interval;
00179     x=a+0.5*hn;
00180     for (sum=0.0,j=1;j<=interval;j++,x+=hn) sum += (ptr->*func)(x);
00181     interval *= 2;
00182     s=0.5*(s+(b-a)*sum/num_interval);
00183     return s;
00184   }
00185 }
00186 
00194 dvariable function_minimizer::trapzd(
00195   dvariable (model_parameters::*func)(const dvariable&),
00196   const dvariable& a, double b, int n)
00197 {
00198   double num_interval;
00199   dvariable x,sum,hn;
00200   static dvariable s;
00201   static int interval;
00202   int j;
00203   model_parameters * ptr= (model_parameters *) mycast();
00204   if (n == 1) {
00205     interval=1;
00206     return (s=0.5*(b-a)*((ptr->*func)(a)+(ptr->*func)(b)));
00207   } else {
00208     num_interval=interval;
00209     hn=(b-a)/num_interval;
00210     x=a+0.5*hn;
00211     for (sum=0.0,j=1;j<=interval;j++,x+=hn) sum += (ptr->*func)(x);
00212     interval *= 2;
00213     s=0.5*(s+(b-a)*sum/num_interval);
00214     return s;
00215   }
00216 }
00217 
00224 dvariable function_minimizer::trapzd(
00225   dvariable (model_parameters::*func)(const dvariable&),
00226   double a, const dvariable& b, int n)
00227 {
00228   double num_interval;
00229   dvariable sum,hn,x;
00230   static dvariable s;
00231 
00232   static int interval;
00233   int j;
00234 
00235   model_parameters * ptr= (model_parameters *) mycast();
00236   if (n == 1) {
00237     interval=1;
00238     return (s=0.5*(b-a)*((ptr->*func)(a)+(ptr->*func)(b)));
00239   } else {
00240     num_interval=interval;
00241     hn=(b-a)/num_interval;
00242     x=a+0.5*hn;
00243     for (sum=0.0,j=1;j<=interval;j++,x+=hn) sum += (ptr->*func)(x);
00244     interval *= 2;
00245     s=0.5*(s+(b-a)*sum/num_interval);
00246     return s;
00247   }
00248 }
00249 
00256 dvariable function_minimizer::trapzd(
00257   dvariable (model_parameters::*func)(const dvariable&),
00258   const dvariable& a, const dvariable& b, int n)
00259 {
00260   double num_interval;
00261   dvariable sum,hn,x;
00262   static dvariable s;
00263   static int interval;
00264   int j;
00265 
00266   model_parameters * ptr= (model_parameters *) mycast();
00267   if (n == 1) {
00268     interval=1;
00269     return (s=0.5*(b-a)*((ptr->*func)(a)+(ptr->*func)(b)));
00270   } else {
00271     num_interval=interval;
00272     hn=(b-a)/num_interval;
00273     x=a+0.5*hn;
00274     for (sum=0.0,j=1;j<=interval;j++,x+=hn) sum += (ptr->*func)(x);
00275     interval *= 2;
00276     s=0.5*(s+(b-a)*sum/num_interval);
00277     return s;
00278   }
00279 }
00280 //#undef EPS
00281 #undef JMAX
00282 //#undef JMAXP
00283 //#undef K
00284 
00285  //Not used elsewhere
00286 void polint(const dvector& xa, const dvar_vector& ya,int n,double x,
00287   const dvariable& _y, const dvariable& _dy)
00288 {
00289   dvariable& y=(dvariable&) _y;
00290   dvariable& dy=(dvariable&) _dy;
00291   int i,m,ns=1;
00292   double dif,dift,ho,hp;
00293   dvariable den,w;
00294 
00295   dif=fabs(x-xa[1]);
00296   dvar_vector c(1,n);
00297   dvar_vector d(1,n);
00298   for (i=1;i<=n;i++) {
00299     if ( (dift=fabs(x-xa[i])) < dif) {
00300       ns=i;
00301       dif=dift;
00302     }
00303     c[i]=ya[i];
00304     d[i]=ya[i];
00305   }
00306   y=ya[ns--];
00307   for (m=1;m<n;m++) {
00308     for (i=1;i<=n-m;i++) {
00309       ho=xa[i]-x;
00310       hp=xa[i+m]-x;
00311       w=c[i+1]-d[i];
00312       if ( (den=ho-hp) == 0.0)
00313       {
00314         cerr << "Error in routine POLINT" << endl;
00315         exit(1);
00316       }
00317       den=w/den;
00318       d[i]=hp*den;
00319       c[i]=ho*den;
00320     }
00321     y += (dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
00322   }
00323 }