Revision 1109 trunk/src/df1b2-separable/df1b2im3f.cpp

df1b2im3f.cpp (revision 1109)
2 2
 * $Id$
3 3
 *
4 4
 * Author: David Fournier
5
 * Copyright (c) 2008-2012 Regents of the University of California 
5
 * Copyright (c) 2008-2012 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
     {
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff