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