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