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