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