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