ADMB Documentation  11.1.1920
 All Classes Files Functions Variables Typedefs Friends Defines
lmnewton.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: lmnewton.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 #ifdef __cplusplus
00010 extern "C" {
00011 #endif
00012 
00013 int lbfgs_(long int *, long int *, double *, double *, double *, long int *,
00014   double *, long int *, double *, double *, double *, long int *,long int *,
00015   long int *);
00016 
00017 #ifdef __cplusplus
00018 }
00019 #endif
00020 
00021 void function_minimizer::limited_memory_quasi_newton(
00022   const independent_variables& _x, int m)
00023 {
00024   independent_variables& x = (independent_variables&) _x;
00025   if (m<=0)
00026   {
00027     cerr << "illegal value for number of steps in limited memory"
00028       " quasi-newton method -- set equal to 10" << endl;
00029   }
00030   int noprintx=0;
00031   int on;
00032   int maxfn_option=0;
00033   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-maxfn"))>-1)
00034   {
00035     maxfn_option=1;
00036   }
00037   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-nox"))>-1)
00038   {
00039     noprintx=1;
00040   }
00041   int nvar=initial_params::nvarcalc(); // get the number of active
00042   double _crit=0;
00043   int itn=0;
00044   int ifn=0;
00045   int nopt = 0;
00046   // set the convergence criterion by command line
00047   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-crit",nopt))>-1)
00048   {
00049     if (!nopt)
00050     {
00051       cerr << "Usage -crit option needs number  -- ignored" << endl;
00052     }
00053     else
00054     {
00055       char * end;
00056       _crit=strtod(ad_comm::argv[on+1],&end);
00057       if (_crit<=0)
00058       {
00059         cerr << "Usage -crit option needs positive number  -- ignored" << endl;
00060         _crit=0.0;
00061       }
00062     }
00063   }
00064 
00065   gradient_structure::set_YES_DERIVATIVES();
00066   // set convergence criterion for this phase
00067   if (!(!convergence_criteria))
00068   {
00069     int ind=min(convergence_criteria.indexmax(),
00070     initial_params::current_phase);
00071     crit=convergence_criteria(ind);
00072   }
00073   if (_crit)
00074   {
00075     crit = _crit;
00076   }
00077   if (!(!maximum_function_evaluations) && !maxfn_option)
00078   {
00079     int ind=min(maximum_function_evaluations.indexmax(),
00080       initial_params::current_phase);
00081     maxfn= (int) maximum_function_evaluations(ind);
00082   }
00083   dvariable vf=0.0;
00084 
00085   double xtol, f;
00086   dvector diag(1,nvar);
00087   //int j, n, iflag, icall;
00088   int iflag, icall;
00089   double fbest=1.e+100;
00090   dvector g(1,nvar);
00091   g.initialize();
00092   dvector xbest(1,nvar);
00093   dvector gbest(1,nvar);
00094   //double t1, t2;
00095   long int diagco=0;
00096   int iprintx[2];
00097   //double epsx;
00098   //m = 35;
00099   dvector w(1,nvar+2*m+2*nvar*m);
00100   nopt=0;
00101   int on1;
00102   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,"-iprint",nopt))>-1)
00103   {
00104     if (!nopt)
00105     {
00106       cerr << "Usage -iprint option needs integer  -- ignored" << endl;
00107     }
00108     else
00109     {
00110       int jj=atoi(ad_comm::argv[on1+1]);
00111       iprint=jj;
00112     }
00113   }
00114   iprintx[0] = iprint;
00115   iprintx[1] = 0;
00116   crit = 1e-5;
00117   xtol = 1e-16;
00118   icall = 0;
00119   iflag = 0;
00120   long int linfo=0;
00121 
00122 L20:
00123   f = 0.;
00124   vf=initial_params::reset(dvar_vector(x));
00125   *objective_function_value::pobjfun=0.0;
00126   userfunction();
00127   vf+=*objective_function_value::pobjfun;
00128   f=value(vf);
00129   ifn++;
00130   if (f<fbest)
00131   {
00132     fbest=f;
00133     xbest=x;
00134     gbest=g;
00135   }
00136 
00137   gradcalc(nvar,g);
00138   if(fmod(double(itn),double(iprint)) == 0)
00139   {
00140     if (iprint>0)
00141     {
00142       if (!itn)
00143       {
00144         if (ad_printf) (*ad_printf)("\nInitial statistics: ");
00145       }
00146       else
00147       {
00148         if (ad_printf) (*ad_printf)("\nIntermediate statistics: ");
00149       }
00150 
00151       if (ad_printf) (*ad_printf)(
00152         "%d variables; iteration %ld; function evaluation %ld\n",
00153         nvar, itn, ifn);
00154 
00155       if (!itn)
00156       {
00157         if (ad_printf) (*ad_printf)(
00158           "Function value %12.4le; maximum gradient component mag %12.4le\n",
00159           f, max(fabs(g)));
00160       }
00161       else
00162       {
00163         if (ad_printf)
00164           (*ad_printf)(
00165             "Function value %12.4le; maximum gradient component mag %12.4le\n",
00166             fbest, max(gbest)
00167           );
00168       }
00169       if (!noprintx)
00170       {
00171         if (!itn)
00172           fmmdisp(x, g, nvar, 0,noprintx);
00173         else
00174           fmmdisp(xbest, gbest, nvar, 0,noprintx);
00175       }
00176     }
00177   }
00178 
00179   long int lnvar=nvar;
00180   long int litn=itn;
00181   long int liprintx= *iprintx;
00182   long int liflag=iflag;
00183   long int lm=m;
00184   lbfgs_(&lnvar, &lm, &(x[1]) , &f, &(g[1]), &diagco, &(diag[1]),
00185     &liprintx, &crit, &xtol, &(w[1]), &liflag,&litn,&linfo);
00186   itn=int(litn);
00187   iflag=int(liflag);
00188 
00189   if (iflag <= 0)
00190   {
00191     goto L50;
00192   }
00193   ++icall;
00194   if (icall > maxfn)
00195   {
00196     goto L50;
00197   }
00198   goto L20;
00199 L50:
00200   if (iprint>0)
00201   {
00202     if (ad_printf) (*ad_printf)("\nfinal statistics: ");
00203 
00204     if (ad_printf) (*ad_printf)(
00205       "%d variables; iteration %ld; function evaluation %ld\n",
00206       nvar, itn, ifn);
00207     if (ad_printf) (*ad_printf)(
00208       "Function value %12.4le; maximum gradient component mag %12.4le\n",
00209       f, max(g));
00210     fmmdisp(x, g, nvar, 0,noprintx);
00211   }
00212   gradient_structure::set_NO_DERIVATIVES();
00213   x=xbest;
00214   objective_function_value::gmax=fabs(max(gbest));
00215 }
00216 
00217 void function_minimizer::limited_memory_quasi_newton
00218   (double& f, const independent_variables& _x, int m, int noprintx,
00219   int maxfn, double crit)
00220 {
00221   independent_variables& x = (independent_variables&) _x;
00222   if (m<=0)
00223   {
00224     cerr << "illegal value for number of steps in limited memory"
00225       " quasi-newton method -- set equal to 10" << endl;
00226   }
00227   int on;
00228 
00229   int nvar=initial_params::nvarcalc(); // get the number of active
00230   double _crit=0;
00231   int itn=0;
00232   int ifn=0;
00233   int nopt = 0;
00234   // set the convergence criterion by command line
00235   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-crit",nopt))>-1)
00236   {
00237     if (!nopt)
00238     {
00239       cerr << "Usage -crit option needs number  -- ignored" << endl;
00240     }
00241     else
00242     {
00243       char* end;
00244       _crit=strtod(ad_comm::argv[on+1],&end);
00245       if (_crit<=0)
00246       {
00247         cerr << "Usage -crit option needs positive number  -- ignored" << endl;
00248       }
00249     }
00250   }
00251   gradient_structure::set_YES_DERIVATIVES();
00252   dvariable vf=0.0;
00253 
00254   double xtol;
00255   dvector diag(1,nvar);
00256   //int j, n, iflag, icall;
00257   int iflag, icall;
00258   double fbest=1.e+100;
00259   dvector g(1,nvar);
00260   g.initialize();
00261   dvector xbest(1,nvar);
00262   dvector gbest(1,nvar);
00263   //double t1, t2;
00264   long int diagco=0;
00265   int iprintx[2];
00266   //double epsx;
00267   //m = 35;
00268   dvector w(1,nvar+2*m+2*nvar*m);
00269   iprintx[0] = iprint;
00270   iprintx[1] = 0;
00271   crit = 1e-5;
00272   xtol = 1e-16;
00273   icall = 0;
00274   iflag = 0;
00275   long int linfo=0;
00276 
00277 L20:
00278   f = 0.;
00279   vf=initial_params::reset(dvar_vector(x));
00280   *objective_function_value::pobjfun=0.0;
00281   userfunction();
00282   vf+=*objective_function_value::pobjfun;
00283   f=value(vf);
00284   ifn++;
00285   if (f<fbest)
00286   {
00287     fbest=f;
00288     xbest=x;
00289     gbest=g;
00290   }
00291 
00292   gradcalc(nvar,g);
00293   if(fmod(double(itn),double(iprint)) == 0)
00294   {
00295     if (iprint>0)
00296     {
00297       if (!itn)
00298       {
00299         if (ad_printf) (*ad_printf)("\nInitial statistics: ");
00300       }
00301       else
00302       {
00303         if (ad_printf) (*ad_printf)("\nIntermediate statistics: ");
00304       }
00305 
00306       if (ad_printf) (*ad_printf)(
00307         "%d variables; iteration %ld; function evaluation %ld\n",
00308         nvar, itn, ifn);
00309 
00310       if (!itn)
00311       {
00312         if (ad_printf) (*ad_printf)(
00313           "Function value %12.4le; maximum gradient component mag %12.4le\n",
00314           f, max(g));
00315       }
00316       else
00317       {
00318         if (ad_printf) (*ad_printf)(
00319           "Function value %12.4le; maximum gradient component mag %12.4le\n",
00320           fbest, max(gbest));
00321       }
00322       if (!noprintx)
00323       {
00324         if (!itn)
00325           fmmdisp(x, g, nvar, 0,noprintx);
00326         else
00327           fmmdisp(xbest, gbest, nvar, 0,noprintx);
00328       }
00329     }
00330   }
00331   long int lnvar=nvar;
00332   long int litn=itn;
00333   long int liprintx= *iprintx;
00334   long int liflag=iflag;
00335   long int lm=m;
00336   lbfgs_(&lnvar, &lm, &(x[1]) , &f, &(g[1]), &diagco, &(diag[1]),
00337     &liprintx, &crit, &xtol, &(w[1]), &liflag,&litn,&linfo);
00338   itn=int(litn);
00339   iflag=int(liflag);
00340 
00341   if (iflag <= 0)
00342   {
00343     goto L50;
00344   }
00345   ++icall;
00346   if (icall > maxfn)
00347   {
00348     goto L50;
00349   }
00350   goto L20;
00351 L50:
00352   if (iprint>0)
00353   {
00354     if (ad_printf)
00355     {
00356       (*ad_printf)("\nfinal statistics: ");
00357       (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00358         nvar, itn, ifn);
00359       (*ad_printf)(
00360         "Function value %12.4le; maximum gradient component mag %12.4le\n",
00361         f, max(g)
00362       );
00363     }
00364     fmmdisp(x, g, nvar, 0,noprintx);
00365   }
00366   //gradient_structure::set_NO_DERIVATIVES();
00367   x=xbest;
00368   f=fbest;
00369   objective_function_value::gmax=fabs(max(gbest));
00370 }