Revision 751
trunk/src/df1b2separable/f1b2fnl3.cpp (revision 751)  

101  101 
} 
102  102  
103  103 
/** 
104 
* Description not yet available. 

105 
* \param 

104 
* Calculates the Laplace approximation for a single separable function 

105 
* in the "block diagonal", i.e. each u(i) occurs only in a single separable function 

106 
* \param ff value of separable function (???) 

107 
* This function will be called multiple times (once for each separable function). 

108 
* Notation: x = fixed effects (parameters) and u = random effects. 

106  109 
*/ 
107  110 
void laplace_approximation_calculator:: 
108  111 
do_separable_stuff_laplace_approximation_block_diagonal(df1b2variable& ff) 
109  112 
{ 
110 
set_dependent_variable(ff); 

113 
set_dependent_variable(ff); // Initializs "dot_bar" for reverse mode AD


111  114 
df1b2_gradlist::set_no_derivatives(); 
112  115 
df1b2variable::passnumber=1; 
113 
df1b2_gradcalc1(); 

116 
df1b2_gradcalc1(); // Forward mode AD follow by a series of reverse sweeps


114  117 

115 
init_df1b2vector & locy= *funnel_init_var::py; 

116 
imatrix& list=*funnel_init_var::plist; 

118 
init_df1b2vector & locy= *funnel_init_var::py; // Independent variables for separable function


119 
imatrix& list=*funnel_init_var::plist; // Index into "locy"


117  120  
118 
int i; int j; int us=0; int xs=0; 

121 
int i; int j; int us=0; int xs=0; // us = #u's and xs = #x's


119  122 
ivector lre_index(1,funnel_init_var::num_active_parameters); 
120  123 
ivector lfe_index(1,funnel_init_var::num_active_parameters); 
121  124  
125 
// count to find us and xs, and find indexes of fixed and random effects 

122  126 
for (i=1;i<=funnel_init_var::num_active_parameters;i++) 
123  127 
{ 
124 
if (list(i,1)>xsize) 

128 
if (list(i,1)>xsize) // x's are stored first in the joint vector


125  129 
{ 
126  130 
lre_index(++us)=i; 
127  131 
} 
...  ...  
131  135 
} 
132  136 
} 
133  137 

134 
dvector local_xadjoint(1,xs); 

138 
dvector local_xadjoint(1,xs); // First order derivative of ff wrt x


135  139 
for (j=1;j<=xs;j++) 
136  140 
{ 
137  141 
int j2=list(lfe_index(j),2); 
138 
local_xadjoint(j)=ff.u_dot[j21]; 

142 
local_xadjoint(j)=ff.u_dot[j21]; // u_dot is the result of forward AD


139  143 
} 
140  144 

141  145 
if (us>0) 
142  146 
{ 
147 
// Find Hessian matrix needed for Laplace approximation 

143  148 
dmatrix local_Hess(1,us,1,us); 
144  149 
dvector local_grad(1,us); 
145  150 
dmatrix local_Dux(1,us,1,xs); 
...  ...  
154  159 
local_Hess(i,j)+=locy(i2).u_bar[j21]; 
155  160 
} 
156  161 
} 
162 


163 
// First order derivative of separable function wrt u 

157  164 
for (i=1;i<=us;i++) 
158  165 
{ 
159  166 
int i2=list(lre_index(i),2); 
160  167 
local_uadjoint(i)= ff.u_dot[i21]; 
161  168 
} 
162  169 

170 
// Mixed derivatives wrt x and u needed in the sensitivity of u_hat wrt x 

163  171 
for (i=1;i<=us;i++) 
164  172 
{ 
165  173 
for (j=1;j<=xs;j++) 
...  ...  
170  178 
} 
171  179 
} 
172  180 

181 
// Enter calculations for the derivative of log(det(Hessian)) 

173  182 

174  183 
//if (initial_df1b2params::separable_calculation_type==3) 
175  184 
{ 
176  185 

177  186 
//int nvar=us*us; 
178 
double f; 

179 
dmatrix Hessadjoint=get_gradient_for_hessian_calcs(local_Hess,f); 

180 
initial_df1b2params::cobjfun+=f; 

187 
double f; // 0.5*log(det(local_Hess))


188 
dmatrix Hessadjoint=get_gradient_for_hessian_calcs(local_Hess,f);


189 
initial_df1b2params::cobjfun+=f; // Adds 0.5*log(det(local_Hess))


181  190 

182  191 
for (i=1;i<=us;i++) 
183  192 
{ 
...  ...  
211  220 
local_uadjoint(i)+=locy[i2].u_tilde[0]; 
212  221 
} 
Also available in: Unified diff