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