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

df1b2im2.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
......
31 31
  gradient_structure::set_YES_DERIVATIVES();
32 32
  int nvar=x.size()+u0.size()+((bw+1)*(2*u0.size()-bw))/2;
33 33
  independent_variables y(1,nvar);
34
  
34

  
35 35
  // need to set random effects active together with whatever
36 36
  // init parameters should be active in this phase
37
  initial_params::set_inactive_only_random_effects(); 
38
  initial_params::set_active_random_effects(); 
39
  /*int onvar=*/initial_params::nvarcalc(); 
37
  initial_params::set_inactive_only_random_effects();
38
  initial_params::set_active_random_effects();
39
  /*int onvar=*/initial_params::nvarcalc();
40 40
  initial_params::xinit(y);    // get the initial values into the
41 41
  y(1,xs)=x;
42 42

  
......
51 51
    y(ii++)=bHess(i,j);
52 52
  }
53 53

  
54
  dvar_vector vy=dvar_vector(y); 
54
  dvar_vector vy=dvar_vector(y);
55 55
  banded_symmetric_dvar_matrix vHess(1,us,bw);
56 56
  initial_params::reset(vy);    // get the initial values into the
57 57
  ii=xs+us+1;
......
65 65
   int nsamp=pmin->lapprox->num_importance_samples;
66 66
   dvar_vector sample_value(1,nsamp);
67 67
   sample_value.initialize();
68
   
68

  
69 69
   int ierr = 0;
70 70
   banded_lower_triangular_dvar_matrix ch=choleski_decomp(vHess,ierr);
71 71
   if (ierr)
......
81 81
     for (int is=1;is<=nsamp;is++)
82 82
     {
83 83
       dvar_vector tau=solve_trans(ch,pmin->lapprox->epsilon(is));
84
      
84

  
85 85
       vy(xs+1,xs+us).shift(1)+=tau;
86 86
       initial_params::reset(vy);    // get the values into the model
87 87
       vy(xs+1,xs+us).shift(1)-=tau;
88
  
88

  
89 89
       *objective_function_value::pobjfun=0.0;
90 90
       pmin->AD_uf_outer();
91
  
91

  
92 92
       sample_value(is)=*objective_function_value::pobjfun
93 93
         -0.5*norm2(pmin->lapprox->epsilon(is));
94 94
     }
......
100 100
     for (is=1;is<=nsamp;is++)
101 101
     {
102 102
       dvar_vector tau=solve_trans(ch,pmin->lapprox->epsilon(is));
103
      
103

  
104 104
       vy(xs+1,xs+us).shift(1)+=tau;
105 105
       initial_params::reset(vy);    // get the values into the model
106 106
       vy(xs+1,xs+us).shift(1)-=tau;
107
  
107

  
108 108
       *objective_function_value::pobjfun=0.0;
109 109
       pmin->AD_uf_outer();
110
  
110

  
111 111
       sample_value(is)=*objective_function_value::pobjfun;
112 112
       normal_weight(is)=0.5*norm2(pmin->lapprox->epsilon(is));
113 113
     }
......
126 126
     cout << "The normalized sample value normal weight pairs are " << endl;
127 127
     for (is=1;is<=nsamp;is++)
128 128
     {
129
       cout << normal_weight(is) << "  " 
129
       cout << normal_weight(is) << "  "
130 130
            << sample_value(is)-normal_weight(is) << endl;
131 131
     }
132 132
     ad_exit(1);
133 133
   }
134 134

  
135 135
   dvariable min_vf=min(sample_value);
136
   vf=min_vf-log(mean(exp(min_vf-sample_value))); 
137
   vf-=us*0.91893853320467241; 
138
   
136
   vf=min_vf-log(mean(exp(min_vf-sample_value)));
137
   vf-=us*0.91893853320467241;
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff