ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
randeff.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: randeff.cpp 2581 2014-11-07 03:21:24Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 #if defined(max)
00010   #undef max
00011 #endif
00012 
00013 #ifdef ISZERO
00014   #undef ISZERO
00015 #endif
00016 #define ISZERO(d) ((d)==0.0)
00017 
00018 typedef int integer;
00019 typedef long int logical;
00020 
00021 #ifdef __cplusplus
00022 extern "C" {
00023 #endif
00024 
00025 int xlbfgs_(integer *n, integer *m, dvar_vector & x, dvariable & f,
00026   dvar_vector & g, logical *diagco, dvar_vector & diag, integer *iprint,
00027   double*  eps, double*  xtol, dvar_vector & w, integer *iflag, integer* iter);
00028 
00029 #ifdef __cplusplus
00030 }
00031 #endif
00032 
00033 dvariable function_minimizer::random_effects_maximization(const dvar_vector& _x)
00034 {
00035   integer m=5;
00036   double crit=0.0;
00037   //double crit=1.e-3;
00038   int maxfn=400;
00039   int maxiter=50;
00040   int iprint=-10;
00041 
00042   dvar_vector& x = (dvar_vector&)(_x);
00043 
00044   integer nvar=x.indexmax()-x.indexmin()+1; // get the number of active
00045   if (m<=0)
00046   {
00047     cerr << "illegal value for number of steps in limited memory"
00048       " quasi-newton method -- set equal to 10" << endl;
00049   }
00050   integer noprintx=0;
00051   int on;
00052   int maxfn_option=0;
00053   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-maxfn"))>-1)
00054   {
00055     maxfn_option=1;
00056   }
00057   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-nox"))>-1)
00058   {
00059     noprintx=1;
00060   }
00061   double _crit=0;
00062   integer itn=0;
00063   int ifn=0;
00064   int nopt = 0;
00065   // set the convergence criterion by command line
00066   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-crit",nopt))>-1)
00067   {
00068     if (!nopt)
00069     {
00070       cerr << "Usage -crit option needs number  -- ignored" << endl;
00071     }
00072     else
00073     {
00074       char * end;
00075       _crit=strtod(ad_comm::argv[on+1],&end);
00076       if (_crit<=0)
00077       {
00078         cerr << "Usage -crit option needs positive number  -- ignored" << endl;
00079         _crit=0.0;
00080       }
00081     }
00082   }
00083   gradient_structure::set_YES_DERIVATIVES();
00084   // set convergence criterion for this phase
00085   dvector convergence_criteria;
00086   if (!(!convergence_criteria))
00087   {
00088     int ind=min(convergence_criteria.indexmax(),
00089     initial_params::current_phase);
00090     crit=convergence_criteria(ind);
00091   }
00092   if (!ISZERO(_crit))
00093   {
00094     crit = _crit;
00095   }
00096   dvector maximum_function_evaluations;
00097   if (!(!maximum_function_evaluations) && !maxfn_option)
00098   {
00099     int ind=min(maximum_function_evaluations.indexmax(),
00100       initial_params::current_phase);
00101     maxfn= (int) maximum_function_evaluations(ind);
00102   }
00103   dvariable vf=0.0;
00104 
00105   double xtol;
00106   dvariable f;
00107   dvar_vector diag(1,nvar);
00108   //int j, n, icall;
00109   int icall;
00110   integer iflag;
00111   dvariable fbest=1.e+100;
00112   dvar_vector g(1,nvar);
00113   dvar_vector xbest(1,nvar);
00114   dvar_vector gbest(1,nvar);
00115   g.initialize();
00116   //double t1, t2;
00117   long int diagco = 0;
00118   integer iprintx[2];
00119   //double epsx;
00120   //m = 35;
00121   dvar_vector w(1,nvar+2*m+2*nvar*m);
00122   iprintx[0] = iprint;
00123   iprintx[1] = 0;
00124   crit = 1e-5;
00125   xtol = 1e-16;
00126   icall = 0;
00127   iflag = 0;
00128 
00129 L20:
00130   f = 0.;
00131   //vf=initial_params::reset(dvar_vector(x));
00132   f=user_randeff(x);
00133   ifn++;
00134   if (f<fbest)
00135   {
00136     fbest=f;
00137     xbest=x;
00138     gbest=g;
00139   }
00140 
00141   //gradcalc(nvar,g);
00142   g=user_dfrandeff(x);
00143 
00144 #if defined(USE_DDOUBLE)
00145 #undef double
00146   if(fmod(double(itn),double(iprint)) == 0)
00147 #define double dd_real
00148 #else
00149   if(fmod(double(itn),double(iprint)) == 0)
00150 #endif
00151   {
00152     if (iprint>0)
00153     {
00154       if (!itn)
00155       {
00156         if (ad_printf) (*ad_printf)("\nInitial statistics: ");
00157       }
00158       else
00159       {
00160         if (ad_printf) (*ad_printf)("\nIntermediate statistics: ");
00161       }
00162 
00163       if (ad_printf)
00164         (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00165         nvar, itn, ifn);
00166 
00167       if (!itn)
00168       {
00169         double xf=value(f);
00170         double xg=max(value(g));
00171         if (ad_printf)
00172           (*ad_printf)(
00173             "Function value %12.4le; maximum gradient component mag %12.4le\n",
00174             xf, xg);
00175       }
00176       else
00177       {
00178         double xf=value(fbest);
00179         double xg=max(value(gbest));
00180         if (ad_printf)
00181           (*ad_printf)(
00182             "Function value %12.4le; maximum gradient component mag %12.4le\n",
00183             xf, xg);
00184       }
00185       if (!noprintx)
00186       {
00187         if (!itn)
00188           fmmdisp(value(x), value(g), nvar, 0,noprintx);
00189         else
00190           fmmdisp(value(xbest), value(gbest), nvar, 0,noprintx);
00191       }
00192     }
00193   }
00194   crit=1.e-15;
00195   xlbfgs_(&nvar, &m, x , f, g, &diagco, diag,
00196     iprintx, &crit, &xtol, w, &iflag,&itn);
00197 
00198   if (iflag <= 0)
00199   {
00200     goto L50;
00201   }
00202   if (iflag > 1)
00203   {
00204     goto L50;
00205   }
00206   ++icall;
00207   if (icall > maxfn)
00208   if (itn > maxiter)
00209   {
00210     cout << "Exceeded maxiter" << endl;
00211     goto L50;
00212   }
00213   goto L20;
00214 L50:
00215   if (f<fbest)
00216   {
00217     fbest=f;
00218     xbest=x;
00219     gbest=g;
00220   }
00221   if (iprint>0)
00222   {
00223     double xf=value(f);
00224     double xg=max(value(g));
00225     if (ad_printf)
00226     {
00227       (*ad_printf)("\nfinal statistics: ");
00228       (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00229         nvar, itn, ifn);
00230       (*ad_printf)(
00231         "Function value %12.4le; maximum gradient component mag %12.4le\n",
00232         xf, xg);
00233     }
00234     fmmdisp(value(x),value(g), nvar, 0,noprintx);
00235   }
00236   //gradient_structure::set_NO_DERIVATIVES();
00237   //objective_function_value::gmax=fabs(max(value(gbest)));
00238   //user_d2frandeff(x);
00239   int sgn=1;
00240   dvar_matrix hess=symmetrize(user_d2frandeff(x));
00241   dvariable ld=ln_det(hess,sgn);
00242   if (sgn==-1)
00243   {
00244     dmatrix v=eigenvectors(value(hess));
00245     dvector e=eigenvalues(value(hess));
00246     cout << endl << e << endl;
00247     for (int i=e.indexmin();i<=e.indexmax();i++)
00248     {
00249       if (e(i)<0.0)
00250         cout << i << " " << v(i) << endl;
00251     }
00252   }
00253   x=xbest;
00254   return fbest+0.5*ln_det(user_d2frandeff(xbest),sgn);
00255 }