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