ADMB Documentation  11.1.2274
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lme.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lme.cpp 1935 2014-04-26 02:02:58Z 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 (var_flag==1)
00097       {
00098         (*Hess_components)(ic)= (Hess-*tmpHess)/0.2;
00099       }
00100       else
00101       {
00102         double dfp=
00103         exp(-2.0*value(*variance_components_vector)(ic));
00104         (*Hess_components)(ic)= (Hess-*tmpHess)/(dfp-df0(ic));
00105       }
00106       (*variance_components_vector)(ic)-=0.2;
00107       break;
00108     default:
00109       cerr << "Illegal value for hesstype here" << endl;
00110       ad_exit(1);
00111     }
00112   }
00113   *variance_components_vector=vsave;
00114 }
00115 
00120 dvar_matrix laplace_approximation_calculator::get_hessian_from_components_lme
00121   (function_minimizer * pfmin)
00122 {
00123   int mmin=variance_components_vector->indexmin();
00124   int mmax=variance_components_vector->indexmax();
00125 
00126   initial_params::set_inactive_only_random_effects();
00127   independent_variables xx(1,mmax-mmin+1);
00128   initial_params::xinit(xx);    // get the initial values into the
00129   dvar_vector vxx=dvar_vector(xx);
00130   initial_params::reset(vxx);    // get current x values into the model
00131   int umin=Hess.indexmin();
00132   int umax=Hess.indexmax();
00133   dvar_matrix vHess(umin,umax,umin,umax);
00134   vHess.initialize();
00135   switch(hesstype)
00136   {
00137   case 3:
00138     cerr << "error -- Has not been impelmented" << endl;
00139     ad_exit(1);
00140     break;
00141   case 4:
00142     {
00143       for(int ic=mmin;ic<=mmax;ic++)
00144       {
00145         if (var_flag==1)
00146         {
00147           vHess+=
00148             (*variance_components_vector)(ic)*((*Hess_components)(ic));
00149         }
00150         else
00151         {
00152           vHess+= exp(-2.0*(*variance_components_vector)(ic))*
00153             ((*Hess_components)(ic));
00154         }
00155       }
00156     }
00157     break;
00158   default:
00159     cerr << "Illegal value for hesstype here" << endl;
00160     ad_exit(1);
00161   }
00162   return vHess;
00163 }
00164 
00165 
00170 dvector laplace_approximation_calculator::banded_calculations_lme
00171   (const dvector& _x,const double& _f,function_minimizer * pfmin)
00172 {
00173   // for use when there is no separability
00174   ADUNCONST(dvector,x)
00175   ADUNCONST(double,f)
00176   //int i,j;
00177 
00178   initial_params::set_inactive_only_random_effects();
00179   gradient_structure::set_NO_DERIVATIVES();
00180   initial_params::reset(x);    // get current x values into the model
00181   gradient_structure::set_YES_DERIVATIVES();
00182 
00183   initial_params::set_active_only_random_effects();
00184   dvector g=get_gradient_lme(pfmin);
00185 
00186   reset_gradient_stack();
00187   // this is the main loop to do inner optimization
00188   //for (i=1;i<=xsize;i++) { y(i)=x(i); }
00189   //for (i=1;i<=usize;i++) { y(i+xsize)=uhat(i); }
00190 
00191   dvar_matrix vHess=get_hessian_from_components_lme(pfmin);
00192   //cout << setprecision(12) << "norm(vHess) = " << norm(value(vHess)) << endl;
00193 
00194   dvariable ld;
00195   dvariable tmp=0.0;
00196   dvariable sgn;
00197 
00198   dvector step=value(solve(vHess,g,tmp,sgn));
00199   if (value(sgn)<=0)
00200   {
00201     cerr << "sgn sucks" << endl;
00202   }
00203   int mmin=variance_components_vector->indexmin();
00204   int mmax=variance_components_vector->indexmax();
00205   int nv=mmax-mmin+1;
00206   dvector g1(1,nv);
00207   ld=0.5*tmp;
00208   gradcalc(nv,g1);
00209   uhat-=step;
00210 
00211   initial_params::set_active_only_random_effects();
00212   double maxg=max(fabs(get_gradient_lme(uhat,pfmin)));
00213   if (maxg > 1.e-12)
00214   {
00215     cout << "max g = " << maxg << endl;
00216   }
00217 
00218   double f2=0.0;
00219   dvector g2=get_gradient_lme_hp(f2,pfmin);
00220   f=f2+value(ld);
00221   return g1+g2;
00222 }
00223 
00228 dvector laplace_approximation_calculator::get_gradient_lme
00229   (function_minimizer * pfmin)
00230 {
00231   dvector g(1,usize);
00232   dvector ub(1,usize);
00233   independent_variables u(1,usize);
00234   gradcalc(0,g);
00235   initial_params::xinit(u);    // get the initial values into the
00236   uhat=u;
00237 
00238   dvariable pen=0.0;
00239   dvariable vf=0.0;
00240   pen=initial_params::reset(dvar_vector(u));
00241   *objective_function_value::pobjfun=0.0;
00242   pfmin->AD_uf_inner();
00243   vf+=*objective_function_value::pobjfun;
00244 
00245   objective_function_value::fun_without_pen=value(vf);
00246   vf+=pen;
00247   gradcalc(usize, g, vf);
00248   return g;
00249 }
00250 
00255 dvector laplace_approximation_calculator::get_gradient_lme
00256   (const dvector& x,function_minimizer * pfmin)
00257 {
00258   dvector g(1,usize);
00259   dvector ub(1,usize);
00260   independent_variables u(1,usize);
00261   u=x;
00262   gradcalc(0,g);
00263 
00264   dvariable pen=0.0;
00265   dvariable vf=0.0;
00266   pen=initial_params::reset(dvar_vector(u));
00267   *objective_function_value::pobjfun=0.0;
00268   pfmin->AD_uf_inner();
00269   vf+=*objective_function_value::pobjfun;
00270 
00271   objective_function_value::fun_without_pen=value(vf);
00272   vf+=pen;
00273   gradcalc(usize, g, vf);
00274   return g;
00275 }
00276 
00281 dvector laplace_approximation_calculator::get_gradient_lme_hp
00282   (const double& _f,function_minimizer * pfmin)
00283 {
00284   ADUNCONST(double,f)
00285   //double fb=1.e+100;
00286   dvector g(1,xsize);
00287   dvector ub(1,xsize);
00288   independent_variables u(1,xsize);
00289   //gradcalc(xsize,g);
00290   initial_params::restore_start_phase();
00291   initial_params::set_inactive_random_effects();
00292   initial_params::xinit(u);    // get the initial values into the
00293 
00294   dvariable pen=0.0;
00295   dvariable vf=0.0;
00296   pen=initial_params::reset(dvar_vector(u));
00297   *objective_function_value::pobjfun=0.0;
00298   pfmin->AD_uf_inner();
00299   vf+=*objective_function_value::pobjfun;
00300 
00301   objective_function_value::fun_without_pen=value(vf);
00302   vf+=pen;
00303   f=value(vf);
00304   gradcalc(xsize,g);
00305   return g;
00306 }