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