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