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