ADMB Documentation  11.5.3254
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2qnm.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::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 = NULL;
00240 
00241       deallocate();
00242       df1b2variable::pool->deallocate();
00243 
00244       for (int i=1;i<df1b2variable::adpool_counter;i++)
00245       {
00246         delete df1b2variable::adpool_vector[i];
00247         df1b2variable::adpool_vector[i]=0;
00248         df1b2variable::nvar_vector[i]=0;
00249       }
00250       df1b2variable::adpool_counter=0;
00251     }
00252     lapprox=new laplace_approximation_calculator(nvar,_unvar,1,nvar+_unvar,
00253       this);
00254     if (lapprox == NULL)
00255     {
00256       cerr << "Error allocating memory for lapprox" << endl;
00257       ad_exit(1);
00258     }
00259     initial_df1b2params::current_phase=initial_params::current_phase;
00260 
00261     initial_df1b2params::save_varsptr();
00262     allocate();
00263     initial_df1b2params::restore_varsptr();
00264 
00265     df1b2_gradlist::set_no_derivatives();
00266     dvector y(1,initial_params::nvarcalc_all());
00267     initial_params::xinit_all(y);
00268     initial_df1b2params::reset_all(y);
00269 
00270     gradient_structure::set_NO_DERIVATIVES();
00271     //vmon_begin();
00272     // see what kind of hessian we are dealing with
00273     if (lapprox->have_users_hesstype==0)
00274     {
00275       if (initial_df1b2params::separable_flag)
00276       {
00277         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-newht",nopt))>-1)
00278         {
00279           lapprox->check_hessian_type2(this);
00280         }
00281         else
00282         {
00283           lapprox->check_hessian_type(this);
00284         }
00285 
00286         // Print Hessian type and related info
00287         switch(lapprox->hesstype)
00288         {
00289         case 1:
00290           cout << "Hessian type 1 " << endl;
00291           break;
00292         case 2:
00293           cout << "\nBlock diagonal Hessian (Block size =" <<
00294             (lapprox->usize)/(lapprox->num_separable_calls) << ")\n" << endl;
00295           break;
00296         case 3:
00297           cout << "\nBanded Hessian (Band width = " << lapprox->bw << ")\n"
00298                << endl;
00299           break;
00300         case 4:
00301           cout << "Hessian type 4 " << endl;
00302           break;
00303         default:
00304           cerr << "This can't happen" << endl;
00305           ad_exit(1);
00306         }
00307       }
00308     }
00309 
00310    /*
00311     cout << "NEED to remove !!!! " << endl;
00312     initial_df1b2params::separable_flag=0;
00313     lapprox->hesstype=1;
00314    */
00315 
00316     // linear mixed effects optimization
00317     if (laplace_approximation_calculator::variance_components_vector)
00318     {
00319 #ifdef DIAG
00320       if (!lapprox)
00321       {
00322         cerr << "this can't happen" << endl;
00323         ad_exit(1);
00324       }
00325 #endif
00326       lapprox->get_hessian_components_banded_lme(this);
00327     }
00328 
00329     if (negdirections==0)
00330     {
00331       while (fmc.ireturn>=0)
00332       {
00333         fmc.fmin(f,x,g);
00334         if (fmc.ireturn>0)
00335         {
00336           if (ifn_trap)
00337           {
00338             if (ifn_trap==fmc.ifn && itn_trap == fmc.itn)
00339             {
00340               cout << "we are here" << endl;
00341             }
00342           }
00343           g=(*lapprox)(x,f,this);
00344           if (bad_step_flag==1)
00345           {
00346             g=1.e+4;
00347             f=2.*fmc.fbest;
00348             bad_step_flag=0;
00349           }
00350 
00351           if (lapprox->init_switch==0)
00352           {
00353             if (f<fmc.fbest)
00354             {
00355               lapprox->ubest=lapprox->uhat;
00356             }
00357           }
00358         }
00359       }
00360     }
00361     else
00362     {
00363       int i;
00364       int nx=negdirections->indexmax();
00365       independent_variables u(1,nx);
00366       dvector x0(1,x.indexmax());
00367       x0=x;
00368       dvector gu(1,nx);
00369       for (i=1;i<=nx;i++) u(i)=1.e-3;
00370       while (fmc.ireturn>=0)
00371       {
00372         fmc.fmin(f,u,gu);
00373         if (fmc.ireturn>0)
00374         {
00375           x=x0;
00376           for (i=1;i<=nx;i++)
00377           {
00378             x+=u(i)*(*negdirections)(i);
00379           }
00380           g=(*lapprox)(x,f,this);
00381           gu=(*negdirections)*g;
00382         }
00383       }
00384       x=x0;
00385       for (i=1;i<=nx;i++)
00386       {
00387         x+=u(i)*(*negdirections)(i);
00388       }
00389       delete negdirections;
00390       negdirections=0;
00391     }
00392     initial_params::set_inactive_only_random_effects();
00393   }
00394 
00395   if (funnel_init_var::py)
00396   {
00397     delete funnel_init_var::py;
00398     funnel_init_var::py=0;
00399   }
00400   gradient_structure::set_NO_DERIVATIVES();
00401   ffbest=fmc.fbest;
00402   g=fmc.gbest(1,fmc.gbest.indexmax());
00403   iexit=fmc.iexit;
00404   ifn=fmc.ifn;
00405   ihflag=fmc.ihflag;
00406   ihang=fmc.ihang;
00407   maxfn_flag=fmc.maxfn_flag;
00408   quit_flag=fmc.quit_flag;
00409   objective_function_value::gmax=fabs(fmc.gmax);
00410 } // end block for quasi newton minimization