ADMB Documentation  11.1.2371
 All Classes Files Functions Variables Typedefs Friends Defines
amoeba.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: amoeba.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB foundation
00006  */
00011 #include <admodel.h>
00012 #include <math.h>
00013 #define NRANSI
00014 //#include "nrutil.h"
00015 #define NMAX 5000
00016 
00020 #define GET_PSUM \
00021           for (j=1;j<=ndim;j++) {\
00022           for (sum=0.0,i=1;i<=mpts;i++) sum += p[i][j];\
00023           psum[j]=sum;}
00024 
00028 #define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
00029 
00030 
00044 void function_minimizer::adamoeba(const dmatrix& _p, const dvector& _y,
00045   int ndim, double ftol,int nfunk)
00046 {
00047   dmatrix& p=(dmatrix&) _p;
00048   dvector& y=(dvector&) _y;
00049   int i,ihi,ilo,inhi,j,mpts=ndim+1;
00050   double rtol,sum,swap,ysave,ytry;
00051 
00052   dvector psum(1,ndim);
00053   nfunk=0;
00054   GET_PSUM
00055   for (;;) {
00056     ilo=1;
00057     ihi = y[1]>y[2] ? (inhi=2,1) : (inhi=1,2);
00058     for (i=1;i<=mpts;i++) {
00059       if (y[i] <= y[ilo]) ilo=i;
00060       if (y[i] > y[ihi]) {
00061         inhi=ihi;
00062         ihi=i;
00063       } else if (y[i] > y[inhi] && i != ihi) inhi=i;
00064     }
00065     rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]));
00066     if (rtol < ftol) {
00067       SWAP(y[1],y[ilo])
00068       for (i=1;i<=ndim;i++) SWAP(p[1][i],p[ilo][i])
00069       break;
00070     }
00071     if (nfunk >= NMAX)
00072     {
00073       cerr << "NMAX exceeded" << endl;
00074     }
00075     nfunk += 2;
00076     ytry=amxxx(p,y,psum,ndim,ihi,-1.0);
00077     if (ytry <= y[ilo])
00078       ytry=amxxx(p,y,psum,ndim,ihi,2.0);
00079     else if (ytry >= y[inhi]) {
00080       ysave=y[ihi];
00081       ytry=amxxx(p,y,psum,ndim,ihi,0.5);
00082       if (ytry >= ysave) {
00083         for (i=1;i<=mpts;i++) {
00084           if (i != ilo) {
00085             for (j=1;j<=ndim;j++)
00086               p[i][j]=psum[j]=0.5*(p[i][j]+p[ilo][j]);
00087 
00088             dvariable vf=0.0;
00089             vf=initial_params::reset(dvar_vector(psum));
00090             *objective_function_value::pobjfun=0.0;
00091             userfunction();
00092             vf+=*objective_function_value::pobjfun;
00093 
00094             y[i]=value(vf);
00095           }
00096         }
00097         nfunk += ndim;
00098         GET_PSUM
00099       }
00100     } else --(nfunk);
00101   }
00102 }
00103 #undef SWAP
00104 #undef GET_PSUM
00105 #undef NMAX
00106 #undef NRANSI
00107 #define NRANSI
00108 
00113 double function_minimizer::amxxx(const dmatrix& _p, const dvector& _y,
00114   const dvector& _psum, int ndim, int ihi, double fac)
00115 {
00116   dmatrix& p=(dmatrix&) _p;
00117   dvector& y=(dvector&) _y;
00118   dvector& psum=(dvector&) _psum;
00119   int j;
00120   double fac1,fac2,ytry;
00121 
00122   dvector ptry(1,ndim);
00123   fac1=(1.0-fac)/ndim;
00124   fac2=fac1-fac;
00125   for (j=1;j<=ndim;j++) ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
00126 
00127   dvariable vf=0.0;
00128   vf=initial_params::reset(dvar_vector(ptry));
00129   *objective_function_value::pobjfun=0.0;
00130   userfunction();
00131   vf+=*objective_function_value::pobjfun;
00132   ytry=value(vf);
00133 
00134   if (ytry < y[ihi]) {
00135     y[ihi]=ytry;
00136     for (j=1;j<=ndim;j++) {
00137       psum[j] += ptry[j]-p[ihi][j];
00138       p[ihi][j]=ptry[j];
00139       }
00140   }
00141   return ytry;
00142 }
00143 #undef NRANSI
00144 
00156 void function_minimizer::neldmead(int n, dvector& _start, dvector& _xmin,
00157   double *ynewlo, double reqmin, double delta,int *icount, int *numres,
00158   int *ifault)
00159 //
00160 //  Purpose:
00161 //
00162 //    NELMIN minimizes a function using the Nelder-Mead algorithm.
00163 //
00164 //  Discussion:
00165 //
00166 //    This routine seeks the minimum value of a user-specified function.
00167 //
00168 //    Simplex function minimisation procedure due to Nelder+Mead(1965),
00169 //    as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with
00170 //    subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976,
00171 //    25, 97) and Hill(1978, 27, 380-2)
00172 //
00173 //    This routine does not include a termination test using the
00174 //    fitting of a quadratic surface.
00175 //
00176 //  Modified:
00177 //
00178 //    October 2010 by Derek Seiple
00179 //
00180 //  Author:
00181 //
00182 //    FORTRAN77 version by R ONeill
00183 //    C++ version by John Burkardt
00184 //
00185 //  Reference:
00186 //
00187 //    John Nelder, Roger Mead,
00188 //    A simplex method for function minimization,
00189 //    Computer Journal,
00190 //    Volume 7, 1965, pages 308-313.
00191 //
00192 //    R ONeill,
00193 //    Algorithm AS 47:
00194 //    Function Minimization Using a Simplex Procedure,
00195 //    Applied Statistics,
00196 //    Volume 20, Number 3, 1971, pages 338-345.
00197 //
00198 //  Parameters:
00199 //
00200 //    Input, int N, the number of variables.
00201 //
00202 //    Input/output, double START[N].  On input, a starting point
00203 //    for the iteration.  On output, this data may have been overwritten.
00204 //
00205 //    Output, double XMIN[N], the coordinates of the point which
00206 //    is estimated to minimize the function.
00207 //
00208 //    Output, double YNEWLO, the minimum value of the function.
00209 //
00210 //    Input, double REQMIN, the terminating limit for the variance
00211 //    of function values.
00212 //
00213 //    Input, int KONVGE, the convergence check is carried out
00214 //    every KONVGE iterations.
00215 //
00216 //    Input, int KCOUNT, the maximum number of function
00217 //    evaluations.
00218 //
00219 //    Output, int *ICOUNT, the number of function evaluations
00220 //    used.
00221 //
00222 //    Output, int *NUMRES, the number of restarts.
00223 //
00224 //    Output, int *IFAULT, error indicator.
00225 //    0, no errors detected.
00226 //    1, REQMIN, N, or KONVGE has an illegal value.
00227 //    2, iteration terminated because KCOUNT was exceeded without convergence.
00228 //
00229 {
00230   dvector& start=(dvector&) _start;
00231   dvector& xmin=(dvector&) _xmin;
00232   int slb = start.indexmin();
00233   int xlb = xmin.indexmin();
00234   start.shift(0);
00235   xmin.shift(0);
00236 
00237   int i;
00238   dvector step(0,n-1);
00239   for(i=0;i<n;i++)
00240   {
00241     step(i)=delta;
00242   }
00243 
00244   int konvge = 10;
00245   int kcount = 5000;
00246   double ccoeff = 0.5;
00247   double del;
00248   double dn;
00249   double dnn;
00250   double ecoeff = 2.0;
00251   double eps = 0.001;
00252 
00253   int ihi;
00254   int ilo;
00255   int j;
00256   int jcount;
00257   int l;
00258   int nn;
00259   double rcoeff = 1.0;
00260   double rq;
00261   double x;
00262   double y2star;
00263   double ylo;
00264   double ystar;
00265   double z;
00266 
00267   // Check the input parameters.
00268   if ( reqmin <= 0.0 )
00269   {
00270     *ifault = 1;
00271     return;
00272   }
00273 
00274   if ( n < 1 )
00275   {
00276     *ifault = 1;
00277     return;
00278   }
00279 
00280   if ( konvge < 1 )
00281   {
00282     *ifault = 1;
00283     return;
00284   }
00285 
00286   dvector p(0,n*(n+1)-1);
00287   dvector pstar(0,n-1);
00288   dvector p2star(0,n-1);
00289   dvector pbar(0,n-1);
00290   dvector y(0,n);
00291 
00292   *icount = 0;
00293   *numres = 0;
00294 
00295   jcount = konvge;
00296   dn = ( double ) ( n );
00297   nn = n + 1;
00298   dnn = ( double ) ( nn );
00299   del = 1.0;
00300   rq = reqmin * dn;
00301 
00302   // Initial or restarted loop.
00303   dvariable vf=0.0;
00304 
00305   for ( ; ; )
00306   {
00307     for ( i = 0; i < n; i++ )
00308     {
00309       p[i+n*n] = start[i];
00310     }
00311     start.shift(1);
00312     vf=0.0;
00313     vf=initial_params::reset(dvar_vector(start));
00314     *objective_function_value::pobjfun=0.0;
00315     userfunction();
00316     vf+=*objective_function_value::pobjfun;
00317     y[n] = value(vf);
00318     start.shift(0);
00319 
00320     *icount = *icount + 1;
00321 
00322     for ( j = 0; j < n; j++ )
00323     {
00324       x = start[j];
00325       start[j] = start[j] + step[j] * del;
00326       for ( i = 0; i < n; i++ )
00327       {
00328         p[i+j*n] = start[i];
00329       }
00330       start.shift(1);
00331       vf=0.0;
00332       vf=initial_params::reset(dvar_vector(start));
00333       *objective_function_value::pobjfun=0.0;
00334       userfunction();
00335       vf+=*objective_function_value::pobjfun;
00336       y[j]=value(vf);
00337       start.shift(0);
00338 
00339       *icount = *icount + 1;
00340       start[j] = x;
00341     }
00342     // The simplex construction is complete.
00343 
00344     // Find highest and lowest Y values.  YNEWLO = Y(IHI) indicates
00345     // the vertex of the simplex to be replaced.
00346     ylo = y[0];
00347     ilo = 0;
00348 
00349     for ( i = 1; i < nn; i++ )
00350     {
00351       if ( y[i] < ylo )
00352       {
00353         ylo = y[i];
00354         ilo = i;
00355       }
00356     }
00357 
00358     // Inner loop.
00359     for ( ; ; )
00360     {
00361       if ( kcount <= *icount )
00362       {
00363         break;
00364       }
00365       *ynewlo = y[0];
00366       ihi = 0;
00367 
00368       for ( i = 1; i < nn; i++ )
00369       {
00370         if ( *ynewlo < y[i] )
00371         {
00372           *ynewlo = y[i];
00373           ihi = i;
00374         }
00375       }
00376 
00377       // Calculate PBAR, the centroid of the simplex vertices
00378       // excepting the vertex with Y value YNEWLO.
00379       for ( i = 0; i < n; i++ )
00380       {
00381         z = 0.0;
00382         for ( j = 0; j < nn; j++ )
00383         {
00384           z = z + p[i+j*n];
00385         }
00386         z = z - p[i+ihi*n];
00387         pbar[i] = z / dn;
00388       }
00389 
00390       // Reflection through the centroid.
00391       for ( i = 0; i < n; i++ )
00392       {
00393         pstar[i] = pbar[i] + rcoeff * ( pbar[i] - p[i+ihi*n] );
00394       }
00395       pstar.shift(1);
00396       vf=0.0;
00397       vf=initial_params::reset(dvar_vector(pstar));
00398       *objective_function_value::pobjfun=0.0;
00399       userfunction();
00400       vf+=*objective_function_value::pobjfun;
00401       ystar = value(vf);
00402       pstar.shift(0);
00403 
00404       *icount = *icount + 1;
00405 
00406       // Successful reflection, so extension.
00407       if ( ystar < ylo )
00408       {
00409         for ( i = 0; i < n; i++ )
00410         {
00411           p2star[i] = pbar[i] + ecoeff * ( pstar[i] - pbar[i] );
00412         }
00413         p2star.shift(1);
00414         vf=0.0;
00415         vf=initial_params::reset(dvar_vector(p2star));
00416         *objective_function_value::pobjfun=0.0;
00417         userfunction();
00418         vf+=*objective_function_value::pobjfun;
00419         y2star = value(vf);
00420         p2star.shift(0);
00421 
00422         *icount = *icount + 1;
00423 
00424         // Check extension.
00425         if ( ystar < y2star )
00426         {
00427           for ( i = 0; i < n; i++ )
00428           {
00429             p[i+ihi*n] = pstar[i];
00430           }
00431           y[ihi] = ystar;
00432         }
00433 
00434         // Retain extension or contraction.
00435         else
00436         {
00437           for ( i = 0; i < n; i++ )
00438           {
00439             p[i+ihi*n] = p2star[i];
00440           }
00441           y[ihi] = y2star;
00442         }
00443       }
00444 
00445       // No extension.
00446       else
00447       {
00448         l = 0;
00449         for ( i = 0; i < nn; i++ )
00450         {
00451           if ( ystar < y[i] )
00452           {
00453             l = l + 1;
00454           }
00455         }
00456 
00457         if ( 1 < l )
00458         {
00459           for ( i = 0; i < n; i++ )
00460           {
00461             p[i+ihi*n] = pstar[i];
00462           }
00463           y[ihi] = ystar;
00464         }
00465 
00466         // Contraction on the Y(IHI) side of the centroid.
00467         else if ( l == 0 )
00468         {
00469           for ( i = 0; i < n; i++ )
00470           {
00471             p2star[i] = pbar[i] + ccoeff * ( p[i+ihi*n] - pbar[i] );
00472           }
00473           p2star.shift(1);
00474           vf=0.0;
00475           vf=initial_params::reset(dvar_vector(p2star));
00476           *objective_function_value::pobjfun=0.0;
00477           userfunction();
00478           vf+=*objective_function_value::pobjfun;
00479           y2star = value(vf);
00480           p2star.shift(0);
00481 
00482           *icount = *icount + 1;
00483 
00484           // Contract the whole simplex.
00485           if ( y[ihi] < y2star )
00486           {
00487             for ( j = 0; j < nn; j++ )
00488             {
00489               for ( i = 0; i < n; i++ )
00490               {
00491                 p[i+j*n] = ( p[i+j*n] + p[i+ilo*n] ) * 0.5;
00492                 xmin[i] = p[i+j*n];
00493               }
00494               xmin.shift(1);
00495               vf=0.0;
00496               vf=initial_params::reset(dvar_vector(xmin));
00497               *objective_function_value::pobjfun=0.0;
00498               userfunction();
00499               vf+=*objective_function_value::pobjfun;
00500               y[j] = value(vf);
00501               xmin.shift(0);
00502 
00503               *icount = *icount + 1;
00504             }
00505             ylo = y[0];
00506             ilo = 0;
00507 
00508             for ( i = 1; i < nn; i++ )
00509             {
00510               if ( y[i] < ylo )
00511               {
00512                 ylo = y[i];
00513                 ilo = i;
00514               }
00515             }
00516             continue;
00517           }
00518 
00519           //  Retain contraction.
00520           else
00521           {
00522             for ( i = 0; i < n; i++ )
00523             {
00524               p[i+ihi*n] = p2star[i];
00525             }
00526             y[ihi] = y2star;
00527           }
00528         }
00529 
00530         // Contraction on the reflection side of the centroid.
00531         else if ( l == 1 )
00532         {
00533           for ( i = 0; i < n; i++ )
00534           {
00535             p2star[i] = pbar[i] + ccoeff * ( pstar[i] - pbar[i] );
00536           }
00537           p2star.shift(1);
00538           vf=0.0;
00539           vf=initial_params::reset(dvar_vector(p2star));
00540           *objective_function_value::pobjfun=0.0;
00541           userfunction();
00542           vf+=*objective_function_value::pobjfun;
00543           y2star = value(vf);
00544           p2star.shift(0);
00545 
00546           *icount = *icount + 1;
00547 
00548           // Retain reflection?
00549           if ( y2star <= ystar )
00550           {
00551             for ( i = 0; i < n; i++ )
00552             {
00553               p[i+ihi*n] = p2star[i];
00554             }
00555             y[ihi] = y2star;
00556           }
00557           else
00558           {
00559             for ( i = 0; i < n; i++ )
00560             {
00561               p[i+ihi*n] = pstar[i];
00562             }
00563             y[ihi] = ystar;
00564           }
00565         }
00566       }
00567 
00568       // Check if YLO improved.
00569       if ( y[ihi] < ylo )
00570       {
00571         ylo = y[ihi];
00572         ilo = ihi;
00573       }
00574       jcount = jcount - 1;
00575 
00576       if ( 0 < jcount )
00577       {
00578         continue;
00579       }
00580 
00581       // Check to see if minimum reached.
00582       if ( *icount <= kcount )
00583       {
00584         jcount = konvge;
00585 
00586         z = 0.0;
00587         for ( i = 0; i < nn; i++ )
00588         {
00589           z = z + y[i];
00590         }
00591         x = z / dnn;
00592 
00593         z = 0.0;
00594         for ( i = 0; i < nn; i++ )
00595         {
00596           z = z + pow ( y[i] - x, 2 );
00597         }
00598 
00599         if ( z <= rq )
00600         {
00601           break;
00602         }
00603       }
00604     }
00605 
00606     // Factorial tests to check that YNEWLO is a local minimum.
00607     for ( i = 0; i < n; i++ )
00608     {
00609       xmin[i] = p[i+ilo*n];
00610     }
00611     *ynewlo = y[ilo];
00612 
00613     if ( kcount < *icount )
00614     {
00615       *ifault = 2;
00616       break;
00617     }
00618 
00619     *ifault = 0;
00620 
00621     for ( i = 0; i < n; i++ )
00622     {
00623       del = step[i] * eps;
00624       xmin[i] = xmin[i] + del;
00625       xmin.shift(1);
00626       vf=0.0;
00627       vf=initial_params::reset(dvar_vector(xmin));
00628       *objective_function_value::pobjfun=0.0;
00629       userfunction();
00630       vf+=*objective_function_value::pobjfun;
00631       z = value(vf);
00632       xmin.shift(0);
00633 
00634       *icount = *icount + 1;
00635       if ( z < *ynewlo )
00636       {
00637         *ifault = 2;
00638         break;
00639       }
00640       xmin[i] = xmin[i] - del - del;
00641       xmin.shift(1);
00642       vf=0.0;
00643       vf=initial_params::reset(dvar_vector(xmin));
00644       *objective_function_value::pobjfun=0.0;
00645       userfunction();
00646       vf+=*objective_function_value::pobjfun;
00647       z = value(vf);
00648       xmin.shift(0);
00649 
00650       *icount = *icount + 1;
00651       if ( z < *ynewlo )
00652       {
00653         *ifault = 2;
00654         break;
00655       }
00656       xmin[i] = xmin[i] + del;
00657     }
00658 
00659     if ( *ifault == 0 )
00660     {
00661       break;
00662     }
00663 
00664     // Restart the procedure.
00665     for ( i = 0; i < n; i++ )
00666     {
00667       start[i] = xmin[i];
00668     }
00669     del = eps;
00670     *numres = *numres + 1;
00671   }
00672 
00673   start.shift(slb);
00674   xmin.shift(xlb);
00675 
00676   return;
00677 }