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