ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2qnm.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2qnm.cpp 2791 2014-12-11 22:24:16Z 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 static int no_stuff=0;
00018 
00023 void function_minimizer::quasi_newton_block(int nvar,int _crit,
00024   independent_variables& x,const dvector& _g,const double& _f)
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 <<
00035         "Usage -fntrap option needs two non-negative integers  -- ignored.\n";
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   fmm fmc(negdirections ? negdirections->indexmax() : nvar);
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     int _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==0)
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 
00284         // Print Hessian type and related info
00285         switch(lapprox->hesstype)
00286         {
00287         case 1:
00288           cout << "Hessian type 1 " << endl;
00289           break;
00290         case 2:
00291           cout << "\nBlock diagonal Hessian (Block size =" <<
00292             (lapprox->usize)/(lapprox->num_separable_calls) << ")\n" << endl;
00293           break;
00294         case 3:
00295           cout << "\nBanded Hessian (Band width = " << lapprox->bw << ")\n"
00296                << endl;
00297           break;
00298         case 4:
00299           cout << "Hessian type 4 " << endl;
00300           break;
00301         default:
00302           cerr << "This can't happen" << endl;
00303           ad_exit(1);
00304         }
00305       }
00306     }
00307 
00308    /*
00309     cout << "NEED to remove !!!! " << endl;
00310     initial_df1b2params::separable_flag=0;
00311     lapprox->hesstype=1;
00312    */
00313 
00314     // linear mixed effects optimization
00315     if (laplace_approximation_calculator::variance_components_vector)
00316     {
00317       if (!lapprox)
00318       {
00319         cerr << "this can't happen" << endl;
00320         ad_exit(1);
00321       }
00322       lapprox->get_hessian_components_banded_lme(this);
00323     }
00324 
00325     if (negdirections==0)
00326     {
00327       while (fmc.ireturn>=0)
00328       {
00329         fmc.fmin(f,x,g);
00330         if (fmc.ireturn>0)
00331         {
00332           if (ifn_trap)
00333           {
00334             if (ifn_trap==fmc.ifn && itn_trap == fmc.itn)
00335             {
00336               cout << "we are here" << endl;
00337             }
00338           }
00339           g=(*lapprox)(x,f,this);
00340           if (bad_step_flag==1)
00341           {
00342             g=1.e+4;
00343             f=2.*fmc.fbest;
00344             bad_step_flag=0;
00345           }
00346 
00347           if (lapprox->init_switch==0)
00348           {
00349             if (f<fmc.fbest)
00350             {
00351               lapprox->ubest=lapprox->uhat;
00352             }
00353           }
00354         }
00355       }
00356     }
00357     else
00358     {
00359       int i;
00360       int nx=negdirections->indexmax();
00361       independent_variables u(1,nx);
00362       dvector x0(1,x.indexmax());
00363       x0=x;
00364       dvector gu(1,nx);
00365       for (i=1;i<=nx;i++) u(i)=1.e-3;
00366       while (fmc.ireturn>=0)
00367       {
00368         fmc.fmin(f,u,gu);
00369         if (fmc.ireturn>0)
00370         {
00371           x=x0;
00372           for (i=1;i<=nx;i++)
00373           {
00374             x+=u(i)*(*negdirections)(i);
00375           }
00376           g=(*lapprox)(x,f,this);
00377           gu=(*negdirections)*g;
00378         }
00379       }
00380       x=x0;
00381       for (i=1;i<=nx;i++)
00382       {
00383         x+=u(i)*(*negdirections)(i);
00384       }
00385       delete negdirections;
00386       negdirections=0;
00387     }
00388     initial_params::set_inactive_only_random_effects();
00389   }
00390 
00391   if (funnel_init_var::py)
00392   {
00393     delete funnel_init_var::py;
00394     funnel_init_var::py=0;
00395   }
00396   gradient_structure::set_NO_DERIVATIVES();
00397   ffbest=fmc.fbest;
00398   g=fmc.gbest(1,fmc.gbest.indexmax());
00399   iexit=fmc.iexit;
00400   ifn=fmc.ifn;
00401   ihflag=fmc.ihflag;
00402   ihang=fmc.ihang;
00403   maxfn_flag=fmc.maxfn_flag;
00404   quit_flag=fmc.quit_flag;
00405   objective_function_value::gmax=fabs(fmc.gmax);
00406 } // end block for quasi newton minimization