## Revision 751

trunk/src/df1b2-separable/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[j2-1];
```
142
```    local_xadjoint(j)=ff.u_dot[j2-1];  // 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[j2-1];
```
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[i2-1];
```
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
```      }
```
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff