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