ADMB Documentation  11.1.1916
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lme.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lme.cpp 1781 2014-03-12 20:36:57Z 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 
00020 void laplace_approximation_calculator::get_hessian_components_banded_lme
00021   (function_minimizer * pfmin)
00022 {
00023   int mmin=variance_components_vector->indexmin();
00024   int mmax=variance_components_vector->indexmax();
00025   dmatrix *tmpHess=0;
00026   switch(hesstype)
00027   {
00028   case 3:
00029     cerr << "error -- Has not been impelmented" << endl;
00030     ad_exit(1);
00031     break;
00032   case 4:
00033     if (Hess_components==0)
00034     {
00035       int umin=Hess.indexmin();
00036       int umax=Hess.indexmax();
00037       tmpHess=new dmatrix(umin,umax,umin,umax);
00038       Hess_components=new d3_array(mmin,mmax,umin,umax,umin,umax);
00039       if (tmpHess==0 || Hess_components==0)
00040       {
00041         cerr << "error allocating memory" << endl;
00042         ad_exit(1);
00043       }
00044       df1b2_gradlist::set_no_derivatives();
00045       int nvar=initial_params::nvarcalc_all();
00046       dvector x(1,nvar);
00047       initial_params::xinit_all(x);
00048       initial_df1b2params::reset_all(x);
00049       for (int i=1;i<=nvar;i++) y(i)=x(i);
00050       step=get_newton_raphson_info_banded(pfmin);
00051       *tmpHess=Hess;
00052     }
00053     break;
00054   default:
00055     cerr << "Illegal value for hesstype here" << endl;
00056     ad_exit(1);
00057   }
00058 
00059   dvector df0=exp(-2.0*value(*variance_components_vector));
00060 
00061   dvector vsave(mmin,mmax);
00062   vsave=value(*variance_components_vector);
00063   for(int ic=mmin;ic<=mmax;ic++)
00064   {
00065     (*variance_components_vector)(ic)+=0.2;
00066 
00067     // test newton raphson
00068     switch(hesstype)
00069     {
00070     case 3:
00071       bHess->initialize();
00072       break;
00073     case 4:
00074       Hess.initialize();
00075       break;
00076     default:
00077       cerr << "Illegal value for hesstype here" << endl;
00078       ad_exit(1);
00079     }
00080 
00081     df1b2_gradlist::set_no_derivatives();
00082     int nvar=initial_params::nvarcalc_all();
00083     dvector x(1,nvar);
00084     initial_params::xinit_all(x);
00085     initial_df1b2params::reset_all(x);
00086     for (int i=1;i<=nvar;i++) y(i)=x(i);
00087     step=get_newton_raphson_info_banded(pfmin);
00088 
00089     switch(hesstype)
00090     {
00091     case 3:
00092       cerr << "error -- Has not been impelmented" << endl;
00093       ad_exit(1);
00094       //(*bHess_components)(ic)= *bHess;
00095       break;
00096     case 4:
00097       if (var_flag==1)
00098       {
00099         (*Hess_components)(ic)= (Hess-*tmpHess)/0.2;
00100       }
00101       else
00102       {
00103         double dfp=
00104         exp(-2.0*value(*variance_components_vector)(ic));
00105         (*Hess_components)(ic)= (Hess-*tmpHess)/(dfp-df0(ic));
00106       }
00107       (*variance_components_vector)(ic)-=0.2;
00108       break;
00109     default:
00110       cerr << "Illegal value for hesstype here" << endl;
00111       ad_exit(1);
00112     }
00113   }
00114   *variance_components_vector=vsave;
00115 }
00116 
00121 dvar_matrix laplace_approximation_calculator::get_hessian_from_components_lme
00122   (function_minimizer * pfmin)
00123 {
00124   int mmin=variance_components_vector->indexmin();
00125   int mmax=variance_components_vector->indexmax();
00126 
00127   initial_params::set_inactive_only_random_effects();
00128   independent_variables xx(1,mmax-mmin+1);
00129   initial_params::xinit(xx);    // get the initial values into the
00130   dvar_vector vxx=dvar_vector(xx);
00131   initial_params::reset(vxx);    // get current x values into the model
00132   int umin=Hess.indexmin();
00133   int umax=Hess.indexmax();
00134   dvar_matrix vHess(umin,umax,umin,umax);
00135   vHess.initialize();
00136   switch(hesstype)
00137   {
00138   case 3:
00139     cerr << "error -- Has not been impelmented" << endl;
00140     ad_exit(1);
00141     break;
00142   case 4:
00143     {
00144       for(int ic=mmin;ic<=mmax;ic++)
00145       {
00146         if (var_flag==1)
00147         {
00148           vHess+=
00149             (*variance_components_vector)(ic)*((*Hess_components)(ic));
00150         }
00151         else
00152         {
00153           vHess+= exp(-2.0*(*variance_components_vector)(ic))*
00154             ((*Hess_components)(ic));
00155         }
00156       }
00157     }
00158     break;
00159   default:
00160     cerr << "Illegal value for hesstype here" << endl;
00161     ad_exit(1);
00162   }
00163   return vHess;
00164 }
00165 
00166 
00171 dvector laplace_approximation_calculator::banded_calculations_lme
00172   (const dvector& _x,const double& _f,function_minimizer * pfmin)
00173 {
00174   // for use when there is no separability
00175   ADUNCONST(dvector,x)
00176   ADUNCONST(double,f)
00177   //int i,j;
00178 
00179   initial_params::set_inactive_only_random_effects();
00180   gradient_structure::set_NO_DERIVATIVES();
00181   initial_params::reset(x);    // get current x values into the model
00182   gradient_structure::set_YES_DERIVATIVES();
00183 
00184   initial_params::set_active_only_random_effects();
00185   dvector g=get_gradient_lme(pfmin);
00186 
00187   reset_gradient_stack();
00188   // this is the main loop to do inner optimization
00189   //for (i=1;i<=xsize;i++) { y(i)=x(i); }
00190   //for (i=1;i<=usize;i++) { y(i+xsize)=uhat(i); }
00191 
00192   dvar_matrix vHess=get_hessian_from_components_lme(pfmin);
00193   //cout << setprecision(12) << "norm(vHess) = " << norm(value(vHess)) << endl;
00194 
00195   dvariable ld;
00196   dvariable tmp=0.0;
00197   dvariable sgn;
00198 
00199   dvector step=value(solve(vHess,g,tmp,sgn));
00200   if (value(sgn)<=0)
00201   {
00202     cerr << "sgn sucks" << endl;
00203   }
00204   int mmin=variance_components_vector->indexmin();
00205   int mmax=variance_components_vector->indexmax();
00206   int nv=mmax-mmin+1;
00207   dvector g1(1,nv);
00208   ld=0.5*tmp;
00209   gradcalc(nv,g1);
00210   uhat-=step;
00211 
00212   initial_params::set_active_only_random_effects();
00213   double maxg=max(fabs(get_gradient_lme(uhat,pfmin)));
00214   if (maxg > 1.e-12)
00215   {
00216     cout << "max g = " << maxg << endl;
00217   }
00218 
00219   double f2=0.0;
00220   dvector g2=get_gradient_lme_hp(f2,pfmin);
00221   f=f2+value(ld);
00222   return g1+g2;
00223 }
00224 
00229 dvector laplace_approximation_calculator::get_gradient_lme
00230   (function_minimizer * pfmin)
00231 {
00232   dvector g(1,usize);
00233   dvector ub(1,usize);
00234   independent_variables u(1,usize);
00235   gradcalc(0,g);
00236   initial_params::xinit(u);    // get the initial values into the
00237   uhat=u;
00238 
00239   dvariable pen=0.0;
00240   dvariable vf=0.0;
00241   pen=initial_params::reset(dvar_vector(u));
00242   *objective_function_value::pobjfun=0.0;
00243   pfmin->AD_uf_inner();
00244   vf+=*objective_function_value::pobjfun;
00245 
00246   objective_function_value::fun_without_pen=value(vf);
00247   vf+=pen;
00248   gradcalc(usize, g, vf);
00249   return g;
00250 }
00251 
00256 dvector laplace_approximation_calculator::get_gradient_lme
00257   (const dvector& x,function_minimizer * pfmin)
00258 {
00259   dvector g(1,usize);
00260   dvector ub(1,usize);
00261   independent_variables u(1,usize);
00262   u=x;
00263   gradcalc(0,g);
00264 
00265   dvariable pen=0.0;
00266   dvariable vf=0.0;
00267   pen=initial_params::reset(dvar_vector(u));
00268   *objective_function_value::pobjfun=0.0;
00269   pfmin->AD_uf_inner();
00270   vf+=*objective_function_value::pobjfun;
00271 
00272   objective_function_value::fun_without_pen=value(vf);
00273   vf+=pen;
00274   gradcalc(usize, g, vf);
00275   return g;
00276 }
00277 
00282 dvector laplace_approximation_calculator::get_gradient_lme_hp
00283   (const double& _f,function_minimizer * pfmin)
00284 {
00285   ADUNCONST(double,f)
00286   //double fb=1.e+100;
00287   dvector g(1,xsize);
00288   dvector ub(1,xsize);
00289   independent_variables u(1,xsize);
00290   //gradcalc(xsize,g);
00291   initial_params::restore_start_phase();
00292   initial_params::set_inactive_random_effects();
00293   initial_params::xinit(u);    // get the initial values into the
00294 
00295   dvariable pen=0.0;
00296   dvariable vf=0.0;
00297   pen=initial_params::reset(dvar_vector(u));
00298   *objective_function_value::pobjfun=0.0;
00299   pfmin->AD_uf_inner();
00300   vf+=*objective_function_value::pobjfun;
00301 
00302   objective_function_value::fun_without_pen=value(vf);
00303   vf+=pen;
00304   f=value(vf);
00305   gradcalc(xsize,g);
00306   return g;
00307 }
00308 #endif   //if defined(USE_LAPLACE)