ADMB Documentation  11.1.1913
 All Classes Files Functions Variables Typedefs Friends Defines
f1b2trst.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: f1b2trst.cpp 1682 2014-02-25 23:43:06Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #if defined(USE_LAPLACE)
00012 #  include <admodel.h>
00013 #  include <df1b2fun.h>
00014 #  include <adrndeff.h>
00015 //#include <vmon.h>
00016 
00021 void function_minimizer::trust_region_update(int nvar,int _crit,
00022   independent_variables& x,const dvector& _g,const double& _f)
00023 {
00024   double & f= (double&)_f;
00025   dvector & g= (dvector&)_g;
00026   fmm fmc(nvar);
00027   if (random_effects_flag)
00028   {
00029     initial_params::set_active_only_random_effects();
00030     //int unvar=initial_params::nvarcalc(); // get the number of active
00031     initial_params::restore_start_phase();
00032     initial_params::set_inactive_random_effects();
00033     int nvar1=initial_params::nvarcalc(); // get the number of active
00034     if (nvar1 != nvar)
00035     {
00036       cerr << "failed sanity check in "
00037        "void function_minimizer::quasi_newton_block" << endl;
00038       ad_exit(1);
00039     }
00040   }
00041 
00042   uistream uis("admodel.hes");
00043   if (!uis)
00044   {
00045     cerr << "Error trying to open file admodel.hes" << endl;
00046     ad_exit(1);
00047   }
00048   int hnvar = 0;
00049   uis >> hnvar;
00050   dmatrix Hess(1,hnvar,1,hnvar);
00051   uis >> Hess;
00052   if (!uis)
00053   {
00054     cerr << "Error trying to read Hessian from admodel.hes" << endl;
00055     ad_exit(1);
00056   }
00057 
00058   dmatrix tester(1,hnvar,1,hnvar);
00059 
00060   tester=Hess;
00061 
00062   dvector e=sort(eigenvalues(Hess));
00063 
00064   double lambda=-e(1)+100.;
00065 
00066   for (int i=1;i<=hnvar;i++)
00067   {
00068     tester(i,i)+=lambda;
00069   }
00070   dvector step =  x-solve(tester,g);
00071 
00072   {
00073     // calculate the number of random effects unvar
00074     // this turns on random effects variables and turns off
00075     // everything else
00076     //cout << nvar << endl;
00077     initial_params::set_active_only_random_effects();
00078     //cout << nvar << endl;
00079     int unvar=initial_params::nvarcalc(); // get the number of active
00080     //df1b2_gradlist::set_no_derivatives();
00081 
00082     if (lapprox)
00083     {
00084       delete lapprox;
00085       lapprox=0;
00086       df1b2variable::pool->deallocate();
00087 
00088       for (int i=0;i<df1b2variable::adpool_counter;i++)
00089       {
00090         delete df1b2variable::adpool_vector[i];
00091         df1b2variable::adpool_vector[i]=0;
00092         df1b2variable::nvar_vector[i]=0;
00093         df1b2variable::adpool_counter=0;
00094       }
00095     }
00096     lapprox=new laplace_approximation_calculator(nvar,unvar,1,nvar+unvar,
00097       this);
00098     initial_df1b2params::current_phase=initial_params::current_phase;
00099 
00100     initial_df1b2params::save_varsptr();
00101     allocate();
00102     initial_df1b2params::restore_varsptr();
00103 
00104     df1b2_gradlist::set_no_derivatives();
00105     int nvar=initial_params::nvarcalc_all();
00106     dvector y(1,nvar);
00107     initial_params::xinit_all(y);
00108     initial_df1b2params::reset_all(y);
00109 
00110     g=(*lapprox)(step,f,this);
00111   }
00112 } // end block for quasi newton minimization
00113 
00114 #endif