ADMB Documentation  11.1x.2738
 All Classes Files Functions Variables Typedefs Friends Defines
f1b2fnl5.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: f1b2fnl5.cpp 2615 2014-11-10 22:35:40Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fnl.h>
00012 #include <adrndeff.h>
00013 #ifndef OPT_LIB
00014   #include <cassert>
00015   #include <climits>
00016 #endif
00017 
00022 void laplace_approximation_calculator::
00023   get_block_diagonal_hessian(df1b2variable& ff)
00024 {
00025   //***********************************************************
00026   //***********************************************************
00027   num_separable_calls++;
00028   set_dependent_variable(ff);
00029   df1b2_gradlist::set_no_derivatives();
00030   df1b2variable::passnumber=1;
00031   df1b2_gradcalc1();
00032 
00033   init_df1b2vector & locy= *funnel_init_var::py;
00034   imatrix& list=*funnel_init_var::plist;
00035   int num_local_re=0;
00036   int num_fixed_effects=0;
00037 
00038   //cout << list << endl;
00039 #ifndef OPT_LIB
00040   assert(funnel_init_var::num_active_parameters <= INT_MAX);
00041 #endif
00042   ivector lre_index(1,(int)funnel_init_var::num_active_parameters);
00043   ivector lfe_index(1,(int)funnel_init_var::num_active_parameters);
00044 
00045   for (int i=1;i<=(int)funnel_init_var::num_active_parameters;i++)
00046   {
00047     if (list(i,1)>xsize)
00048     {
00049       lre_index(++num_local_re)=i;
00050     }
00051     else if (list(i,1)>0)
00052     {
00053       lfe_index(++num_fixed_effects)=i;
00054     }
00055   }
00056 
00057   if (num_local_re > 0)
00058   {
00059     dmatrix& local_Hess=(*block_diagonal_hessian)(num_separable_calls);
00060     dmatrix& local_Dux=(*block_diagonal_Dux)(num_separable_calls);
00061     ivector& local_re_list=(*block_diagonal_re_list)(num_separable_calls);
00062     ivector& local_fe_list=(*block_diagonal_fe_list)(num_separable_calls);
00063     local_re_list.initialize();
00064     local_fe_list.initialize();
00065     local_Hess.initialize();
00066     for (int i=1;i<=num_local_re;i++)
00067     {
00068       local_re_list(i)=list(lre_index(i),1);
00069     }
00070 
00071     for (int i=1;i<=num_fixed_effects;i++)
00072     {
00073       local_fe_list(i)=list(lfe_index(i),1);
00074     }
00075 
00076     for (int i=1;i<=num_local_re;i++)
00077     {
00078       int lrei=lre_index(i);
00079       for (int j=1;j<=num_local_re;j++)
00080       {
00081         int lrej=lre_index(j);
00082         int i2=list(lrei,2);
00083         int j2=list(lrej,2);
00084         local_Hess(i,j)+=locy(i2).u_bar[j2-1];
00085       }
00086     }
00087 
00088     for (int i=1;i<=num_local_re;i++)
00089     {
00090       for (int j=1;j<=num_fixed_effects;j++)
00091       {
00092         int i2=list(lre_index(i),2);
00093         int j2=list(lfe_index(j),2);
00094         local_Dux(i,j)=locy(i2).u_bar[j2-1];
00095       }
00096     }
00097 
00098     have_bounded_random_effects=0;
00099     if (have_bounded_random_effects)
00100     {
00101       for (int i=1;i<=num_local_re;i++)
00102       {
00103         int lrei=lre_index(i);
00104         int i1=list(lrei,1);
00105         for (int j=1;j<=num_local_re;j++)
00106         {
00107           int lrej=lre_index(j);
00108           int j1=list(lrej,1);
00109           local_Hess(i,j)*=scale(i1-xsize)*scale(j1-xsize);
00110         }
00111       }
00112     }
00113   }
00114 
00115   f1b2gradlist->reset();
00116   f1b2gradlist->list.initialize();
00117   f1b2gradlist->list2.initialize();
00118   f1b2gradlist->list3.initialize();
00119   f1b2gradlist->nlist.initialize();
00120   f1b2gradlist->nlist2.initialize();
00121   f1b2gradlist->nlist3.initialize();
00122   funnel_init_var::num_vars=0;
00123   funnel_init_var::num_active_parameters=0;
00124   funnel_init_var::num_inactive_vars=0;
00125 }
00126 
00131 void laplace_approximation_calculator::
00132   do_separable_stuff_laplace_approximation_importance_sampling_adjoint
00133   (df1b2variable& ff)
00134 {
00135   num_separable_calls++;
00136   set_dependent_variable(ff);
00137   //df1b2_gradlist::set_no_derivatives();
00138   df1b2variable::passnumber=1;
00139   df1b2_gradcalc1();
00140 
00141   init_df1b2vector & locy= *funnel_init_var::py;
00142   imatrix& list=*funnel_init_var::plist;
00143 
00144   int us=0; int xs=0;
00145 #ifndef OPT_LIB
00146   assert(funnel_init_var::num_active_parameters <= INT_MAX);
00147 #endif
00148   ivector lre_index(1,(int)funnel_init_var::num_active_parameters);
00149   ivector lfe_index(1,(int)funnel_init_var::num_active_parameters);
00150 
00151   for (int i=1;i<=(int)funnel_init_var::num_active_parameters;i++)
00152   {
00153     if (list(i,1)>xsize)
00154     {
00155       lre_index(++us)=i;
00156     }
00157     else if (list(i,1)>0)
00158     {
00159       lfe_index(++xs)=i;
00160     }
00161   }
00162 
00163   dvector local_xadjoint(1,xs);
00164 
00165   if (us>0)
00166   {
00167     dmatrix local_Hess(1,us,1,us);
00168     dvector local_grad(1,us);
00169     dmatrix local_Dux(1,us,1,xs);
00170     local_Hess.initialize();
00171     dvector local_uadjoint(1,us);
00172     for (int i=1;i<=us;i++)
00173     {
00174       for (int j=1;j<=us;j++)
00175       {
00176         int i2=list(lre_index(i),2);
00177         int j2=list(lre_index(j),2);
00178         local_Hess(i,j)+=locy(i2).u_bar[j2-1];
00179       }
00180     }
00181 
00182     for (int i=1;i<=us;i++)
00183     {
00184       for (int j=1;j<=xs;j++)
00185       {
00186         int i2=list(lre_index(i),2);
00187         int j2=list(lfe_index(j),2);
00188         local_Dux(i,j)=locy(i2).u_bar[j2-1];
00189       }
00190     }
00191 
00192     double f=0.0;
00193     initial_df1b2params::cobjfun+=f;
00194 
00195     for (int i=1;i<=us;i++)
00196     {
00197       for (int j=1;j<=us;j++)
00198       {
00199         int i2=list(lre_index(i),2);
00200         int j2=list(lre_index(j),2);
00201         locy(i2).get_u_bar_tilde()[j2-1]=
00202           (*block_diagonal_vhessianadjoint)(num_separable_calls)(i,j);
00203       }
00204     }
00205 
00206     df1b2variable::passnumber=2;
00207     df1b2_gradcalc1();
00208 
00209     df1b2variable::passnumber=3;
00210     df1b2_gradcalc1();
00211     dvector xtmp(1,xs);
00212     xtmp.initialize();
00213 
00214     local_uadjoint.initialize();
00215     local_xadjoint.initialize();
00216     for (int i=1;i<=xs;i++)
00217     {
00218       int i2=list(lfe_index(i),2);
00219       xtmp(i)+=locy[i2].u_tilde[0];
00220       local_xadjoint(i)+=locy[i2].u_tilde[0];
00221     }
00222     dvector utmp(1,us);
00223     utmp.initialize();
00224     for (int i=1;i<=us;i++)
00225     {
00226       int i2=list(lre_index(i),2);
00227       utmp(i)+=locy[i2].u_tilde[0];
00228       local_uadjoint(i)+=locy[i2].u_tilde[0];
00229     }
00230     if (xs>0)
00231       local_xadjoint -= solve(local_Hess,local_uadjoint)*local_Dux;
00232   }
00233   for (int i=1;i<=xs;i++)
00234   {
00235     int ii=lfe_index(i);
00236     check_local_xadjoint2(list(ii,1))+=local_xadjoint(i);
00237   }
00238   f1b2gradlist->reset();
00239   f1b2gradlist->list.initialize();
00240   f1b2gradlist->list2.initialize();
00241   f1b2gradlist->list3.initialize();
00242   f1b2gradlist->nlist.initialize();
00243   f1b2gradlist->nlist2.initialize();
00244   f1b2gradlist->nlist3.initialize();
00245   funnel_init_var::num_vars=0;
00246   funnel_init_var::num_active_parameters=0;
00247   funnel_init_var::num_inactive_vars=0;
00248 }