Revision 1109 trunk/src/df1b2separable/df1b2im3f.cpp
df1b2im3f.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) 
...  ...  
101  101 
} 
102  102  
103  103 
dvariable vf=0.0; 
104 


104  
105  105 
int nsamp=pmin>lapprox>num_importance_samples; 
106  106 
dvar_vector sample_value(1,nsamp); 
107  107 
sample_value.initialize(); 
...  ...  
137  137 
pmin>lapprox>epsilon(is)(offset+1,offset+lus).shift(1); 
138  138 
offset+=lus; 
139  139 
} 
140 


140  
141  141 
// have to reorder the terms to match the block diagonal hessian 
142  142 
imatrix & ls=*(pmin>lapprox>block_diagonal_re_list); 
143  143 
int mmin=ls.indexmin(); 
144  144 
int mmax=ls.indexmax(); 
145 


145  
146  146 
int ii=1; 
147  147 
int i; 
148  148 
for (i=mmin;i<=mmax;i++) 
...  ...  
170  170 
vy(ls(i,j))=tau(ii++); 
171  171 
} 
172  172 
} 
173 


173  
174  174 
*objective_function_value::pobjfun=0.0; 
175  175 
pmin>AD_uf_outer(); 
176 


176  
177  177 
if (pmin>lapprox>use_outliers==0) 
178  178 
{ 
179  179 
sample_value(icount)=*objective_function_value::pobjfun 
...  ...  
182  182 
else 
183  183 
{ 
184  184 
dvector& e=pmin>lapprox>epsilon(is); 
185 


185  
186  186 
sample_value(icount)=*objective_function_value::pobjfun 
187 
+sum(log(.95*exp(0.5*square(e))+.0166666667*exp(square(e)/18.0)))


187 
+sum(log(.95*exp(0.5*square(e))+.0166666667*exp(square(e)/18.0))) 

188  188 
.91893853320467274177; 
189  189 
} 
190  190 
} 
191 


191  
192  192 
if (icount>0) 
193  193 
{ 
Also available in: Unified diff