ADMB Documentation  11.1.1920
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp2.cpp 1919 2014-04-22 22:02:01Z 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 void get_second_ders(int xs,int us,const init_df1b2vector y,dmatrix& Hess,
00016   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00017   laplace_approximation_calculator* lap);
00018 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
00019   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00020   const dmatrix& _Hessadjoint,function_minimizer * pmin);
00021 
00022 static void xxx(ivector re_list,ivector fe_list){}
00023 
00028 dvector laplace_approximation_calculator::block_diagonal_calculations
00029   (const dvector& _x,const double& _f,function_minimizer * pfmin)
00030 {
00031   // for use when there is no separability
00032   ADUNCONST(dvector,x)
00033   ADUNCONST(double,f)
00034   //int i,j;
00035   int i;
00036 
00037   initial_params::set_inactive_only_random_effects();
00038   gradient_structure::set_NO_DERIVATIVES();
00039   initial_params::reset(x);    // get current x values into the model
00040   gradient_structure::set_YES_DERIVATIVES();
00041 
00042   initial_params::set_active_only_random_effects();
00043   //int lmn_flag=0;
00044   if (!inner_lmnflag)
00045   {
00046     if (!ADqd_flag)
00047       uhat=get_uhat_quasi_newton_block_diagonal(x,pfmin);
00048     else
00049       uhat=get_uhat_quasi_newton_qd(x,pfmin);
00050   }
00051   else
00052   {
00053     uhat=get_uhat_lm_newton(x,pfmin);
00054   }
00055   if (!allocated(scale))
00056   {
00057     scale.allocate(1,uhat.indexmax());
00058   }
00059   else
00060   {
00061     if (scale.indexmax() != uhat.indexmax())
00062     {
00063       scale.deallocate();
00064       scale.allocate(1,uhat.indexmax());
00065     }
00066   }
00067 
00068   if (!allocated(curv))
00069   {
00070     curv.allocate(1,uhat.indexmax());
00071   }
00072   else
00073   {
00074     if (curv.indexmax() != uhat.indexmax())
00075     {
00076       curv.deallocate();
00077       curv.allocate(1,uhat.indexmax());
00078     }
00079   }
00080 
00081   if (sparse_hessian_flag==0)
00082   {
00083     for (i=1;i<=xsize;i++)
00084     {
00085       y(i)=x(i);
00086     }
00087     for (int i=1;i<=usize;i++)
00088     {
00089       y(i+xsize)=uhat(i);
00090     }
00091   }
00092   else
00093   {
00094     for (i=1;i<=xsize;i++)
00095     {
00096       value(y(i))=x(i);
00097     }
00098     for (int i=1;i<=usize;i++)
00099     {
00100       value(y(i+xsize))=uhat(i);
00101     }
00102   }
00103 
00104   for(int ii=1;ii<=num_nr_iters;ii++)
00105   {
00106     {
00107       // test newton raphson
00108       //Hess.initialize();
00109       /*int check=*/initial_params::stddev_scale(scale,uhat);
00110       /*check=*/initial_params::stddev_curvscale(curv,uhat);
00111       max_separable_g=0.0;
00112       pmin->inner_opt_flag=1;
00113       step=get_newton_raphson_info_block_diagonal(pfmin);
00114       cout << "max separable g " << max_separable_g << endl;
00115       cout << "Newton raphson " << ii << endl;
00116       uhat+=step;
00117 
00118       evaluate_function(uhat,pfmin);
00119       pmin->inner_opt_flag=0;
00120     }
00121 
00122     if (sparse_hessian_flag==0)
00123     {
00124       for (int i=1;i<=usize;i++)
00125       {
00126         y(i+xsize)=uhat(i);
00127       }
00128     }
00129     else
00130     {
00131       for (int i=1;i<=usize;i++)
00132       {
00133         value(y(i+xsize))=uhat(i);
00134       }
00135     }
00136   }
00137 
00138   cout << initial_df1b2params::cobjfun << endl;
00139   xadjoint.initialize();
00140   uadjoint.initialize();
00141   block_diagonal_flag=2;
00142   used_flags.initialize();
00143   funnel_init_var::lapprox=this;
00144   if (use_gauss_hermite>0)
00145   {
00146     df1b2variable pen=0.0;
00147     initial_df1b2params::reset(y,pen);
00148     initial_params::straight_through_flag=0;
00149     block_diagonal_flag=6;
00150     num_separable_calls=0;
00151     // get the block diagonal hessians to use in importance sampling
00152     pfmin->user_function();
00153     //cout << (*pfmin->lapprox->block_diagonal_hessian) << endl;
00154     block_diagonal_flag=2;
00155     initial_params::straight_through_flag=0;
00156 
00157     // do importance sampling and get ders bakc to Hessian adjoint
00158     // new stuff for more than one random effect in each separable call
00159     //  Apr 17 07
00160     if (multi_random_effects==0)
00161     {
00162       f=do_gauss_hermite_block_diagonal(x,uhat,Hess,xadjoint,
00163         uadjoint,Hessadjoint,pfmin);
00164     }
00165     else
00166     {
00167       f=do_gauss_hermite_block_diagonal_multi(x,uhat,Hess,xadjoint,
00168         uadjoint,Hessadjoint,pfmin);
00169     }
00170     int xmax=xadjoint.indexmax();
00171     dvector x_con(1,xmax);
00172     x_con.initialize();
00173     dvector dscale(1,nvar);   // need to get scale from somewhere
00174     dscale=1.0;
00175     /*int check=*/initial_params::stddev_scale(dscale,x);
00176     dvector sscale=dscale(1,xsize);
00177     // *******************************************************
00178     // *******************************************************
00179     // *******************************************************
00180     // derivatives from hessian adjoint back
00181     {
00182       x_con.initialize();
00183 
00184       for (int i=1;i<=num_separable_calls;i++)
00185       {
00186         ivector& re_list=(*block_diagonal_re_list)(i);
00187         ivector& fe_list=(*block_diagonal_fe_list)(i);
00188         dmatrix& Dux=(*block_diagonal_Dux)(i);
00189         dmatrix& H=(*block_diagonal_hessian)(i);
00190         xxx(re_list,fe_list);
00191         int mmax=re_list.indexmax();
00192         dvector tmp(1,mmax);
00193 
00194         int j;
00195         for (j=1;j<=re_list.indexmax();j++)
00196         {
00197           tmp(j)=uadjoint(re_list(j)-xmax);
00198         }
00199 
00200         if (allocated(fe_list))
00201         {
00202           if (allocated(H))
00203           {
00204             dvector tmp1=solve(H,tmp);
00205             dvector xtmp=tmp1*Dux;
00206             for (j=1;j<=fe_list.indexmax();j++)
00207             {
00208               x_con(fe_list(j))+=xtmp(j);
00209             }
00210           }
00211         }
00212       }
00213       if (initial_df1b2params::separable_flag)
00214       {
00215         x_con=elem_prod(x_con,sscale);
00216       }
00217     }
00218     xadjoint-=x_con;
00219     // *******************************************************
00220     // *******************************************************
00221     // *******************************************************
00222 
00223     block_diagonal_flag=3;
00224     //pfmin->lapprox->xadjoint.initialize();
00225     //pfmin->lapprox->uadjoint.initialize();
00226     pfmin->lapprox->num_separable_calls=0;
00227     pfmin->lapprox->check_local_xadjoint.initialize();
00228     pfmin->lapprox->check_local_xadjoint2.initialize();
00229     pfmin->lapprox->check_local_uadjoint.initialize();
00230     pfmin->lapprox->check_local_uadjoint2.initialize();
00231     //df1b2_gradlist::set_yes_derivatives();
00232     //initial_df1b2params::reset(y,pen);
00233     pfmin->user_function();
00234     dvector lcx=elem_prod(check_local_xadjoint2,sscale);
00235     xadjoint+=lcx;
00236     //df1b2_gradlist::set_no_derivatives();
00237     funnel_init_var::lapprox=0;
00238     block_diagonal_flag=0;
00239     initial_params::set_inactive_only_random_effects();
00240   }
00241   else if (num_importance_samples>0)
00242   {
00243     df1b2variable pen=0.0;
00244     initial_df1b2params::reset(y,pen);
00245     initial_params::straight_through_flag=0;
00246     block_diagonal_flag=6;
00247     num_separable_calls=0;
00248     // get the block diagonal hessians to use in importance sampling
00249     pfmin->user_function();
00250     //cout << (*pfmin->lapprox->block_diagonal_hessian) << endl;
00251     block_diagonal_flag=2;
00252     initial_params::straight_through_flag=0;
00253     // do importance sampling and get ders bakc to Hessian adjoint
00254     if (isfunnel_flag==0)
00255     {
00256       if (antiflag==0)
00257       {
00258         f=calculate_importance_sample_block_diagonal_option2(x,uhat,Hess,
00259           xadjoint,uadjoint,Hessadjoint,pfmin);
00260       }
00261       else
00262       {
00263         f=calculate_importance_sample_block_diagonal_option_antithetical(x,
00264           uhat,Hess,xadjoint,uadjoint,Hessadjoint,pfmin);
00265       }
00266     }
00267     else
00268     {
00269       f=calculate_importance_sample_block_diagonal_funnel(x,uhat,Hess,xadjoint,
00270         uadjoint,Hessadjoint,pfmin);
00271     }
00272 
00273     int xmax=xadjoint.indexmax();
00274     dvector x_con(1,xmax);
00275     x_con.initialize();
00276     dvector dscale(1,nvar);   // need to get scale from somewhere
00277     dscale=1.0;
00278     /*int check=*/initial_params::stddev_scale(dscale,x);
00279     dvector sscale=dscale(1,xsize);
00280     // *******************************************************
00281     // *******************************************************
00282     // *******************************************************
00283     // derivatives from hessian adjoint back
00284     {
00285       x_con.initialize();
00286 
00287       for (int i=1;i<=num_separable_calls;i++)
00288       {
00289         dmatrix& H=(*block_diagonal_hessian)(i);
00290         if (allocated(H))
00291         {
00292           ivector& re_list=(*block_diagonal_re_list)(i);
00293           ivector& fe_list=(*block_diagonal_fe_list)(i);
00294           dmatrix& Dux=(*block_diagonal_Dux)(i);
00295           xxx(re_list,fe_list);
00296           int mmax=re_list.indexmax();
00297           dvector tmp(1,mmax);
00298 
00299           int j;
00300           for (j=1;j<=re_list.indexmax();j++)
00301           {
00302             tmp(j)=uadjoint(re_list(j)-xmax);
00303           }
00304 
00305           if (allocated(fe_list))
00306           {
00307             dvector tmp1=solve(H,tmp);
00308             dvector xtmp=tmp1*Dux;
00309             for (j=1;j<=fe_list.indexmax();j++)
00310             {
00311               x_con(fe_list(j))+=xtmp(j);
00312             }
00313           }
00314         }
00315       }
00316       if (initial_df1b2params::separable_flag)
00317       {
00318         //for (i=1;i<=usize;i++)
00319         //{
00320         //  Dux(i)=elem_prod(Dux(i),sscale);
00321         //}
00322         x_con=elem_prod(x_con,sscale);
00323       }
00324     }
00325     xadjoint-=x_con;
00326     // *******************************************************
00327     // *******************************************************
00328     // *******************************************************
00329 
00330     block_diagonal_flag=3;
00331     //pfmin->lapprox->xadjoint.initialize();
00332     //pfmin->lapprox->uadjoint.initialize();
00333     pfmin->lapprox->num_separable_calls=0;
00334     pfmin->lapprox->check_local_xadjoint.initialize();
00335     pfmin->lapprox->check_local_xadjoint2.initialize();
00336     pfmin->lapprox->check_local_uadjoint.initialize();
00337     pfmin->lapprox->check_local_uadjoint2.initialize();
00338     //df1b2_gradlist::set_yes_derivatives();
00339     //initial_df1b2params::reset(y,pen);
00340     pfmin->user_function();
00341     dvector lcx=elem_prod(check_local_xadjoint2,sscale);
00342     xadjoint+=lcx;
00343     //df1b2_gradlist::set_no_derivatives();
00344     funnel_init_var::lapprox=0;
00345     block_diagonal_flag=0;
00346     initial_params::set_inactive_only_random_effects();
00347   }
00348   else
00349   {
00350     if (function_minimizer::first_hessian_flag)
00351     {
00352       // need hessin of random effects for stddeV report
00353       df1b2variable pen=0.0;
00354       initial_df1b2params::reset(y,pen);
00355       initial_params::straight_through_flag=0;
00356       block_diagonal_flag=6;
00357       allocate_block_diagonal_stuff();
00358       num_separable_calls=0;
00359       // get the block diagonal hessians to use in importance sampling
00360       pfmin->user_function();
00361       //cout << (*pfmin->lapprox->block_diagonal_hessian) << endl;
00362       block_diagonal_flag=2;
00363     }
00364     initial_params::straight_through_flag=1;
00365     df1b2variable pen=0.0;
00366     if (saddlepointflag==2)
00367     {
00368       pmin->inner_opt_flag=0;
00369       f=get_fx_fu(pfmin);
00370     }
00371     initial_df1b2params::reset(y,pen);
00372     pmin->inner_opt_flag=1;
00373     pfmin->pre_user_function();
00374     pmin->inner_opt_flag=0;
00375     initial_params::straight_through_flag=0;
00376     if (saddlepointflag!=2)
00377     {
00378       f=initial_df1b2params::cobjfun;
00379     }
00380     else
00381     {
00382       xadjoint=(*grad_x_u)(1,xsize)-(*grad_x);
00383     }
00384 
00385     if (saddlepointflag!=2 && pmin->multinomial_weights==0)
00386     {
00387       f-=usize*.91893853320467241;
00388     }
00389 
00390     funnel_init_var::lapprox=0;
00391     block_diagonal_flag=0;
00392     dvector scale1(1,xsize);   // need to get scale from somewhere
00393     initial_params::set_inactive_only_random_effects();
00394     /*int check=*/initial_params::stddev_scale(scale1,x);
00395     for (i=1;i<=xadjoint.indexmax();i++)
00396       xadjoint(i)*=scale1(i);
00397   }
00398   //cout << initial_df1b2params::cobjfun << endl;
00399   //f=initial_df1b2params::cobjfun;
00400   return xadjoint;
00401 }
00402 
00407 dvector laplace_approximation_calculator::get_newton_raphson_info_block_diagonal
00408   (function_minimizer * pfmin)
00409 {
00410   int nv=initial_df1b2params::set_index();
00411   if (allocated(used_flags))
00412   {
00413     if (used_flags.indexmax() != nv)
00414     {
00415       used_flags.safe_deallocate();
00416     }
00417   }
00418   if (!allocated(used_flags))
00419   {
00420     used_flags.safe_allocate(1,nv);
00421   }
00422 
00423   //for (ip=1;ip<=num_der_blocks;ip++)
00424   {
00425     used_flags.initialize();
00426 
00427     // do we need to reallocate memory for df1b2variables?
00428     int ip = 0;
00429     check_for_need_to_reallocate(ip);
00430 
00431     df1b2_gradlist::set_no_derivatives();
00432     //cout << re_objective_function_value::pobjfun << endl;
00433     //cout << re_objective_function_value::pobjfun->ptr << endl;
00434     (*re_objective_function_value::pobjfun)=0;
00435     df1b2variable pen=0.0;
00436     df1b2variable zz=0.0;
00437 
00438     initial_df1b2params::reset(y,pen);
00439 
00440     // call function to do block diagonal newton-raphson
00441     // the step vector from the newton-raphson is in the vector step
00442     df1b2_gradlist::set_yes_derivatives();
00443 
00444     funnel_init_var::lapprox=this;
00445     //cout << funnel_init_var::lapprox << endl;
00446     block_diagonal_flag=1;
00447     pfmin->pre_user_function();
00448     //pfmin->user_function();
00449     funnel_init_var::lapprox=0;
00450     block_diagonal_flag=0;
00451   }
00452   return step;
00453 }
00454 
00455 #endif //#if defined(USE_LAPLACE)