ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lmn2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lmn2.cpp 2792 2014-12-12 19:47:41Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <sstream>
00012 using std::istringstream;
00013 
00014 #  include <admodel.h>
00015 #  include <df1b2fun.h>
00016 #  include <adrndeff.h>
00017 //#include <vmon.h>
00018 static int no_stuff=0;
00019 //static void xxxy(void) {}
00020 
00025 void function_minimizer::limited_memory_quasi_newton_block(int nvar,int _crit,
00026   independent_variables& x,const dvector& _g,const double& _f,int nsteps)
00027 {
00028   int ifn_trap=0;
00029   int itn_trap=0;
00030   int on=0;
00031   int nopt=0;
00032   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-fntrap",nopt))>-1)
00033   {
00034     if (nopt !=2)
00035     {
00036       cerr << "Usage -fntrap option needs two non-negative integers  -- ignored"
00037            << endl;
00038     }
00039     else
00040     {
00041       ifn_trap=atoi(ad_comm::argv[on+1]);
00042       itn_trap=atoi(ad_comm::argv[on+2]);
00043     }
00044   }
00045 
00046   double & f= (double&)_f;
00047   dvector & g= (dvector&)_g;
00048   // *********************************************************
00049   // block for quasi-newton minimization
00050   //int itnold=0;
00051   fmmt1 fmc(negdirections ? negdirections->indexmax() : nvar, nsteps);
00052   int on1;
00053   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,"-nox"))>-1)
00054   {
00055     fmc.noprintx=1;
00056   }
00057   fmc.maxfn= maxfn;
00058   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,"-dd",nopt))>-1)
00059   {
00060     if (!nopt)
00061     {
00062       cerr << "Usage -iprint option needs integer  -- ignored" << endl;
00063       //fmc.iprint=iprint;
00064     }
00065     else
00066     {
00067       int jj=atoi(ad_comm::argv[on1+1]);
00068       fmc.dcheck_flag=jj;
00069     }
00070   }
00071   nopt=0;
00072   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,"-iprint",nopt))>-1)
00073   {
00074     if (!nopt)
00075     {
00076       cerr << "Usage -iprint option needs integer  -- ignored" << endl;
00077       fmc.iprint=iprint;
00078     }
00079     else
00080     {
00081       int jj=atoi(ad_comm::argv[on1+1]);
00082       fmc.iprint=jj;
00083     }
00084   }
00085   else
00086   {
00087     fmc.iprint= iprint;
00088   }
00089   fmc.crit = crit;
00090   fmc.imax = imax;
00091   fmc.dfn= dfn;
00092 
00093   double _dfn=1.e-2;
00094   if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-dfn",nopt))>-1)
00095   {
00096     if (!nopt)
00097     {
00098       cerr << "Usage -dfn option needs number  -- ignored" << endl;
00099     }
00100     else
00101     {
00102       istringstream ist(ad_comm::argv[on+1]);
00103       ist >> _dfn;
00104 
00105       if (_dfn<0)
00106       {
00107         cerr << "Usage -dfn option needs positive number  -- ignored" << endl;
00108         _dfn=0.0;
00109       }
00110     }
00111   }
00112   if (_dfn>=0)
00113   {
00114     fmc.dfn=_dfn;
00115   }
00116 
00117   fmc.scroll_flag= scroll_flag;
00118   fmc.min_improve=min_improve;
00119   gradient_structure::set_YES_DERIVATIVES();
00120   // set convergence criterion for this phase
00121   if (_crit)
00122   {
00123     fmc.crit = _crit;
00124   }
00125   if (!(!convergence_criteria))
00126   {
00127     int ind=min(convergence_criteria.indexmax(),
00128       initial_params::current_phase);
00129     fmc.crit=convergence_criteria(ind);
00130   }
00131   if (!(!maximum_function_evaluations))
00132   {
00133     int ind=min(maximum_function_evaluations.indexmax(),
00134       initial_params::current_phase);
00135     fmc.maxfn= (int) maximum_function_evaluations(ind);
00136   }
00137   int unvar=1;
00138   if (random_effects_flag)
00139   {
00140     initial_params::set_active_only_random_effects();
00141     //cout << nvar << endl;
00142     unvar=initial_params::nvarcalc(); // get the number of active
00143     initial_params::restore_start_phase();
00144     initial_params::set_inactive_random_effects();
00145     int nvar1=initial_params::nvarcalc(); // get the number of active
00146     if (nvar1 != nvar)
00147     {
00148       cerr << "failed sanity check in "
00149        "void function_minimizer::quasi_newton_block" << endl;
00150       ad_exit(1);
00151     }
00152   }
00153 
00154   if (!random_effects_flag || !unvar)
00155   {
00156     dvariable xf=initial_params::reset(dvar_vector(x));
00157     reset_gradient_stack();
00158     gradcalc(0,g);
00159     if (negdirections==0)
00160     {
00161       while (fmc.ireturn>=0)
00162       {
00163         fmc.fmin(f,x,g);
00164         if (fmc.ireturn>0)
00165         {
00166           dvariable vf=0.0;
00167           vf=initial_params::reset(dvar_vector(x));
00168           *objective_function_value::pobjfun=0.0;
00169           pre_userfunction();
00170           if ( no_stuff ==0 && quadratic_prior::get_num_quadratic_prior()>0)
00171           {
00172             quadratic_prior::get_M_calculations();
00173           }
00174           vf+=*objective_function_value::pobjfun;
00175           f=value(vf);
00176           gradcalc(nvar,g);
00177         }
00178       }
00179     }
00180     else
00181     {
00182       int i;
00183       int nx=negdirections->indexmax();
00184       independent_variables u(1,nx);
00185       dvector gu(1,nx);
00186       for (i=1;i<=nx;i++) u(i)=1.e-3;
00187       dvector x0(1,x.indexmax());
00188       x0=x;
00189       while (fmc.ireturn>=0)
00190       {
00191         fmc.fmin(f,u,gu);
00192         if (fmc.ireturn>0)
00193         {
00194           dvariable vf=0.0;
00195           x=x0;
00196           for (i=1;i<=nx;i++)
00197           {
00198             x+=u(i)*(*negdirections)(i);
00199           }
00200           vf=initial_params::reset(dvar_vector(x));
00201           *objective_function_value::pobjfun=0.0;
00202           pre_userfunction();
00203           if ( no_stuff ==0 && quadratic_prior::get_num_quadratic_prior()>0)
00204           {
00205             quadratic_prior::get_M_calculations();
00206           }
00207           vf+=*objective_function_value::pobjfun;
00208           f=value(vf);
00209           gradcalc(nvar,g);
00210           gu=(*negdirections)*g;
00211         }
00212       }
00213       x=x0;
00214       for (i=1;i<=nx;i++)
00215       {
00216         x+=u(i)*(*negdirections)(i);
00217       }
00218       delete negdirections;
00219       negdirections=0;
00220     }
00221   }
00222   else
00223   {
00224     // calculate the number of random effects unvar
00225     // this turns on random effects variables and turns off
00226     // everything else
00227     //cout << nvar << endl;
00228     initial_params::set_active_only_random_effects();
00229     //cout << nvar << endl;
00230     unvar=initial_params::nvarcalc(); // get the number of active
00231     //df1b2_gradlist::set_no_derivatives();
00232 
00233     if (funnel_init_var::py)
00234     {
00235       delete funnel_init_var::py;
00236       funnel_init_var::py=0;
00237     }
00238     if (lapprox)
00239     {
00240       delete lapprox;
00241       lapprox=0;
00242       df1b2variable::pool->deallocate();
00243 
00244       for (int i=1;i<df1b2variable::adpool_counter;i++)
00245       {
00246         delete df1b2variable::adpool_vector[i];
00247         df1b2variable::adpool_vector[i]=0;
00248         df1b2variable::nvar_vector[i]=0;
00249       }
00250       df1b2variable::adpool_counter=0;
00251     }
00252     lapprox=new laplace_approximation_calculator(nvar,unvar,1,nvar+unvar,
00253       this);
00254     if (lapprox==0)
00255     {
00256       cerr << "Error allocating memory for lapprox" << endl;
00257       ad_exit(1);
00258     }
00259     initial_df1b2params::current_phase=initial_params::current_phase;
00260 
00261     initial_df1b2params::save_varsptr();
00262     allocate();
00263     initial_df1b2params::restore_varsptr();
00264 
00265     df1b2_gradlist::set_no_derivatives();
00266     dvector y(1,initial_params::nvarcalc_all());
00267     initial_params::xinit_all(y);
00268     initial_df1b2params::reset_all(y);
00269 
00270     gradient_structure::set_NO_DERIVATIVES();
00271     //vmon_begin();
00272     // see what kind of hessian we are dealing with
00273     if (lapprox->have_users_hesstype==0)
00274     {
00275       if (initial_df1b2params::separable_flag)
00276       {
00277         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-newht",nopt))>-1)
00278         {
00279           lapprox->check_hessian_type2(this);
00280         }
00281         else
00282         {
00283           lapprox->check_hessian_type(this);
00284         }
00285         cout << "Hessian type = " << lapprox->hesstype << endl;
00286       }
00287     }
00288 
00289    /*
00290     cout << "NEED to remove !!!! " << endl;
00291     initial_df1b2params::separable_flag=0;
00292     lapprox->hesstype=1;
00293    */
00294 
00295     // linear mixed effects optimization
00296     if (laplace_approximation_calculator::variance_components_vector)
00297     {
00298       if (!lapprox)
00299       {
00300         cerr << "this can't happen" << endl;
00301         ad_exit(1);
00302       }
00303       lapprox->get_hessian_components_banded_lme(this);
00304     }
00305 
00306     if (negdirections==0)
00307     {
00308       while (fmc.ireturn>=0)
00309       {
00310         fmc.fmin(f,x,g);
00311         if (fmc.ireturn>0)
00312         {
00313           if (ifn_trap)
00314           {
00315             if (ifn_trap==fmc.ifn && itn_trap == fmc.itn)
00316             {
00317               cout << "we are here" << endl;
00318             }
00319           }
00320           g=(*lapprox)(x,f,this);
00321          /*   don't have this in official yet
00322           if (bad_step_flag==1)
00323           {
00324             g=1.e+4;
00325             f=100.*fmc.fbest;
00326           }
00327          */
00328           if (lapprox->init_switch==0)
00329           {
00330             if (f<fmc.fbest)
00331             {
00332               lapprox->ubest=lapprox->uhat;
00333             }
00334           }
00335         }
00336       }
00337     }
00338     else
00339     {
00340       int i;
00341       int nx=negdirections->indexmax();
00342       independent_variables u(1,nx);
00343       dvector x0(1,x.indexmax());
00344       x0=x;
00345       dvector gu(1,nx);
00346       for (i=1;i<=nx;i++) u(i)=1.e-3;
00347       while (fmc.ireturn>=0)
00348       {
00349         fmc.fmin(f,u,gu);
00350         if (fmc.ireturn>0)
00351         {
00352           x=x0;
00353           for (i=1;i<=nx;i++)
00354           {
00355             x+=u(i)*(*negdirections)(i);
00356           }
00357           g=(*lapprox)(x,f,this);
00358           gu=(*negdirections)*g;
00359         }
00360       }
00361       x=x0;
00362       for (i=1;i<=nx;i++)
00363       {
00364         x+=u(i)*(*negdirections)(i);
00365       }
00366       delete negdirections;
00367       negdirections=0;
00368     }
00369     initial_params::set_inactive_only_random_effects();
00370   }
00371 
00372   if (funnel_init_var::py)
00373   {
00374     delete funnel_init_var::py;
00375     funnel_init_var::py=0;
00376   }
00377   gradient_structure::set_NO_DERIVATIVES();
00378   ffbest=fmc.fbest;
00379   g=fmc.gbest(1,fmc.gbest.indexmax());
00380   iexit=fmc.iexit;
00381   ifn=fmc.ifn;
00382   ihflag=fmc.ihflag;
00383   ihang=fmc.ihang;
00384   maxfn_flag=fmc.maxfn_flag;
00385   quit_flag=fmc.quit_flag;
00386   objective_function_value::gmax=fabs(fmc.gmax);
00387 } // end block for quasi newton minimization