Revision 1109 trunk/src/df1b2separable/df1b2ghmult.cpp
df1b2ghmult.cpp (revision 1109)  

2  2 
* $Id$ 
3  3 
* 
4  4 
* Author: David Fournier 
5 
* Copyright (c) 20082012 Regents of the University of California


5 
* Copyright (c) 20082012 Regents of the University of California 

6  6 
*/ 
7  7 
/** 
8  8 
* \file 
...  ...  
33  33 
int hroom = sum(square(lrea)); 
34  34 
int nvar=x.size()+u0.size()+hroom; 
35  35 
independent_variables y(1,nvar); 
36 


36  
37  37 
// need to set random effects active together with whatever 
38  38 
// init parameters should be active in this phase 
39 
initial_params::set_inactive_only_random_effects();


40 
initial_params::set_active_random_effects();


41 
/*int onvar=*/initial_params::nvarcalc();


39 
initial_params::set_inactive_only_random_effects(); 

40 
initial_params::set_active_random_effects(); 

41 
/*int onvar=*/initial_params::nvarcalc(); 

42  42 
initial_params::xinit(y); // get the initial values into the 
43  43 
// do we need this next line? 
44  44 
y(1,xs)=x; 
...  ...  
53  53 
quadratic_prior::get_cHessian_contribution(Hess,vxs); 
54  54 
} 
55  55 
// Here need hooks for sparse matrix structures 
56 


56  
57  57 
dvar3_array & block_diagonal_vhessian= 
58  58 
*pmin>lapprox>block_diagonal_vhessian; 
59  59 
block_diagonal_vhessian.initialize(); 
...  ...  
74  74 
dvector g(1,nvar); 
75  75 
gradcalc(0,g); 
76  76 
gradient_structure::set_YES_DERIVATIVES(); 
77 
dvar_vector vy=dvar_vector(y);


77 
dvar_vector vy=dvar_vector(y); 

78  78 
//initial_params::stddev_vscale(d,vy); 
79  79 
ii=xs+us+1; 
80  80 
if (initial_df1b2params::have_bounded_random_effects) 
...  ...  
117  117 
delete pmin>lapprox>gh>mi; 
118  118 
pmin>lapprox>gh>mi=0; 
119  119 
} 
120 


120  
121  121 
pmin>lapprox>gh>mi=new multi_index(1,nsamp, 
122  122 
pmin>lapprox>multi_random_effects); 
123  123  
...  ...  
144  144 
xv(iu)= xx(mi()(iu)); 
145  145 
} 
146  146 
tau(offset+1,offset+lus).shift(1)=block_diagonal_ch(ic)*xv; 
147 


147  
148  148 
offset+=lus; 
149  149 
} 
150  150 
} 
151 


151  
152  152 
// have to reorder the terms to match the block diagonal hessian 
153  153 
imatrix & ls=*(pmin>lapprox>block_diagonal_re_list); 
154  154 
int mmin=ls.indexmin(); 
155  155 
int mmax=ls.indexmax(); 
156 


156  
157  157 
int ii=1; 
158  158 
int i; 
159  159 
for (i=mmin;i<=mmax;i++) 
...  ...  
219  219 
} 
220  220  
221  221 
vf+=ld; 
222 
//vf+=us*0.91893853320467241;


222 
//vf+=us*0.91893853320467241; 

223  223  
224  224 
double f=value(vf); 
225  225 
gradcalc(nvar,g); 
...  ...  
229  229 
vy(xs+1,xs+us).shift(1)=u0; 
230  230 
initial_params::reset(vy); // get the values into the model 
231  231 
gradient_structure::set_YES_DERIVATIVES(); 
232 


232  
233  233 
pmin>lapprox>in_gauss_hermite_phase=0; 
234 


234  
235  235 
ii=1; 
236  236 
for (i=1;i<=xs;i++) 
237  237 
xadjoint(i)=g(ii++); 
Also available in: Unified diff