ADMB Documentation  11.1.2490
 All Classes Files Functions Variables Typedefs Friends Defines
f1b2fnl3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: f1b2fnl3.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   do_separable_stuff_x_u_block_diagonal(df1b2variable& ff)
00020 {
00021    // we need g_xu g_uu f_x and f_u
00022   set_dependent_variable(ff);
00023   df1b2_gradlist::set_no_derivatives();
00024   df1b2variable::passnumber=1;
00025   df1b2_gradcalc1();
00026 
00027   init_df1b2vector & locy= *funnel_init_var::py;
00028   imatrix& list=*funnel_init_var::plist;
00029 
00030   int i; int j; int us=0; int xs=0;
00031   ivector lre_index(1,funnel_init_var::num_active_parameters);
00032   ivector lfe_index(1,funnel_init_var::num_active_parameters);
00033 
00034   for (i=1;i<=funnel_init_var::num_active_parameters;i++)
00035   {
00036     if (list(i,1)>xsize)
00037     {
00038       lre_index(++us)=i;
00039     }
00040     else if (list(i,1)>0)
00041     {
00042       lfe_index(++xs)=i;
00043     }
00044   }
00045 
00046   dvector local_xadjoint(1,xs);
00047   dvector local_uadjoint(1,us);
00048   for (i=1;i<=xs;i++)
00049   {
00050     int ii=lfe_index(i);
00051     local_xadjoint(i)=(*grad_x_u)(list(ii,1));
00052   }
00053   for (i=1;i<=us;i++)
00054   {
00055     int ii=lre_index(i);
00056     local_uadjoint(i)=(*grad_x_u)(list(ii,1));
00057   }
00058   dvector tmp;
00059   if (us>0)
00060   {
00061     dmatrix local_Hess(1,us,1,us);
00062     dvector local_grad(1,us);
00063     dmatrix local_Dux(1,us,1,xs);
00064     local_Hess.initialize();
00065     for (i=1;i<=us;i++)
00066     {
00067       for (j=1;j<=us;j++)
00068       {
00069         int i2=list(lre_index(i),2);
00070         int j2=list(lre_index(j),2);
00071         local_Hess(i,j)+=locy(i2).u_bar[j2-1];
00072       }
00073     }
00074     for (i=1;i<=us;i++)
00075     {
00076       for (j=1;j<=xs;j++)
00077       {
00078         int i2=list(lre_index(i),2);
00079         int j2=list(lfe_index(j),2);
00080         local_Dux(i,j)=locy(i2).u_bar[j2-1];
00081       }
00082     }
00083     tmp=solve(local_Hess,local_uadjoint)*local_Dux;
00084   }
00085 
00086   for (i=1;i<=xs;i++)
00087   {
00088     int ii=lfe_index(i);
00089     (*grad_x)(list(ii,1))+=tmp(i);
00090   }
00091   f1b2gradlist->reset();
00092   f1b2gradlist->list.initialize();
00093   f1b2gradlist->list2.initialize();
00094   f1b2gradlist->list3.initialize();
00095   f1b2gradlist->nlist.initialize();
00096   f1b2gradlist->nlist2.initialize();
00097   f1b2gradlist->nlist3.initialize();
00098   funnel_init_var::num_vars=0;
00099   funnel_init_var::num_active_parameters=0;
00100   funnel_init_var::num_inactive_vars=0;
00101 }
00102 
00111 void laplace_approximation_calculator::
00112   do_separable_stuff_laplace_approximation_block_diagonal(df1b2variable& ff)
00113 {
00114   set_dependent_variable(ff);// Initializs "dot_bar" for reverse mode AD
00115   df1b2_gradlist::set_no_derivatives();
00116   df1b2variable::passnumber=1;
00117   df1b2_gradcalc1();// Forward mode AD follow by a series of reverse sweeps
00118 
00119   // Independent variables for separable function
00120   init_df1b2vector& locy= *funnel_init_var::py;
00121   imatrix& list=*funnel_init_var::plist; // Index into "locy"
00122 
00123   int i; int j; int us=0; int xs=0; // us = #u's and xs = #x's
00124   ivector lre_index(1,funnel_init_var::num_active_parameters);
00125   ivector lfe_index(1,funnel_init_var::num_active_parameters);
00126 
00127   // count to find us and xs, and find indexes of fixed and random effects
00128   for (i=1;i<=funnel_init_var::num_active_parameters;i++)
00129   {
00130     if (list(i,1)>xsize) // x's are stored first in the joint vector
00131     {
00132       lre_index(++us)=i;
00133     }
00134     else if (list(i,1)>0)
00135     {
00136       lfe_index(++xs)=i;
00137     }
00138   }
00139 
00140   dvector local_xadjoint(1,xs);  // First order derivative of ff wrt x
00141   for (j=1;j<=xs;j++)
00142   {
00143     int j2=list(lfe_index(j),2);
00144     local_xadjoint(j)=ff.u_dot[j2-1];  // u_dot is the result of forward AD
00145   }
00146 
00147   if (us>0)
00148   {
00149     // Find Hessian matrix needed for Laplace approximation
00150     dmatrix local_Hess(1,us,1,us);
00151     dvector local_grad(1,us);
00152     dmatrix local_Dux(1,us,1,xs);
00153     local_Hess.initialize();
00154     dvector local_uadjoint(1,us);
00155     for (i=1;i<=us;i++)
00156     {
00157       for (j=1;j<=us;j++)
00158       {
00159         int i2=list(lre_index(i),2);
00160         int j2=list(lre_index(j),2);
00161         local_Hess(i,j)+=locy(i2).u_bar[j2-1];
00162       }
00163     }
00164 
00165     // First order derivative of separable function wrt u
00166     for (i=1;i<=us;i++)
00167     {
00168       int i2=list(lre_index(i),2);
00169       local_uadjoint(i)= ff.u_dot[i2-1];
00170     }
00171 
00172     // Mixed derivatives wrt x and u needed in the sensitivity of u_hat wrt x
00173     for (i=1;i<=us;i++)
00174     {
00175       for (j=1;j<=xs;j++)
00176       {
00177         int i2=list(lre_index(i),2);
00178         int j2=list(lfe_index(j),2);
00179         local_Dux(i,j)=locy(i2).u_bar[j2-1];
00180       }
00181     }
00182 
00183     // Enter calculations for the derivative of log(det(Hessian))
00184 
00185     //if (initial_df1b2params::separable_calculation_type==3)
00186     {
00187     //int nvar=us*us;
00188     double f;// 0.5*log(det(local_Hess))
00189     dmatrix Hessadjoint=get_gradient_for_hessian_calcs(local_Hess,f);
00190     initial_df1b2params::cobjfun+=f;  // Adds 0.5*log(det(local_Hess))
00191 
00192     for (i=1;i<=us;i++)
00193     {
00194       for (j=1;j<=us;j++)
00195       {
00196         int i2=list(lre_index(i),2);
00197         int j2=list(lre_index(j),2);
00198         locy(i2).get_u_bar_tilde()[j2-1]=Hessadjoint(i,j);
00199       }
00200     }
00201 
00202      df1b2variable::passnumber=2;
00203      df1b2_gradcalc1();
00204 
00205      df1b2variable::passnumber=3;
00206      df1b2_gradcalc1();
00207       dvector xtmp(1,xs);
00208       xtmp.initialize();
00209       for (i=1;i<=xs;i++)
00210       {
00211         int i2=list(lfe_index(i),2);
00212         xtmp(i)+=locy[i2].u_tilde[0];
00213         local_xadjoint(i)+=locy[i2].u_tilde[0];
00214       }
00215       dvector utmp(1,us);
00216       utmp.initialize();
00217       for (i=1;i<=us;i++)
00218       {
00219         int i2=list(lre_index(i),2);
00220         utmp(i)+=locy[i2].u_tilde[0];
00221         local_uadjoint(i)+=locy[i2].u_tilde[0];
00222       }
00223       if (xs>0)
00224         local_xadjoint -= local_uadjoint*inv(local_Hess)*local_Dux;
00225     }
00226   }
00227   for (i=1;i<=xs;i++)
00228   {
00229     int ii=lfe_index(i);
00230     // Ads contribution to "global" gradient vector
00231     xadjoint(list(ii,1))+=local_xadjoint(i);
00232   }
00233   f1b2gradlist->reset();
00234   f1b2gradlist->list.initialize();
00235   f1b2gradlist->list2.initialize();
00236   f1b2gradlist->list3.initialize();
00237   f1b2gradlist->nlist.initialize();
00238   f1b2gradlist->nlist2.initialize();
00239   f1b2gradlist->nlist3.initialize();
00240   funnel_init_var::num_vars=0;
00241   funnel_init_var::num_active_parameters=0;
00242   funnel_init_var::num_inactive_vars=0;
00243 }
00244 
00249 dmatrix laplace_approximation_calculator::get_gradient_for_hessian_calcs
00250   (const dmatrix& local_Hess,double & f)
00251 {
00252   int us=local_Hess.indexmax();
00253   int nvar=us*us;
00254   independent_variables cy(1,nvar);
00255   cy.initialize();
00256   int ii=1; int i,j;
00257   for (i=1;i<=us;i++)
00258     for (j=1;j<=us;j++)
00259       cy(ii++)=local_Hess(i,j);
00260 
00261   dvar_vector vy=dvar_vector(cy);
00262   dvar_matrix vHess(1,us,1,us);
00263 
00264   ii=1;
00265   for (i=1;i<=us;i++)
00266     for (j=1;j<=us;j++)
00267       vHess(i,j)=vy(ii++);
00268 
00269   dvariable vf=0.0;
00270   int sgn=0;
00271   if (pmin->multinomial_weights==0)
00272   {
00273     vf+=0.5*ln_det(vHess,sgn);
00274   }
00275   else
00276   {
00277     dvector & w= *(pmin->multinomial_weights);
00278     double w_i=w[separable_calls_counter];
00279     double d=vHess.indexmax()-vHess.indexmin()+1;
00280     vf+=w_i*(0.5*ln_det(vHess,sgn)-0.5*d*log(w_i));
00281     vf-=w_i*d*.91893853320467241;
00282   }
00283   f=value(vf);
00284   dvector g(1,nvar);
00285   gradcalc(nvar,g);
00286   dmatrix hessadjoint(1,us,1,us);
00287   ii=1;
00288   for (i=1;i<=us;i++)
00289     for (j=1;j<=us;j++)
00290       hessadjoint(i,j)=g(ii++);
00291 
00292   return hessadjoint;
00293 }