EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GaussianRegGradHessian.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GaussianRegGradHessian.cpp
2 #include "Pincushion.h"
3 
4 #include <Eigen/LU>
5 
6 #include <cmath>
7 #include <memory>
8 
9 using namespace std;
10 using namespace Eigen;
11 using namespace SeamStress;
12 
13 namespace FitNewton
14 {
15  GaussianRegGradHessian::GaussianRegGradHessian(FunctionGradHessian* func_instance, double sig, unsigned long int numthreads) : FunctionGradHessian(func_instance->nPars(), 0), func(func_instance), data_has_errors(false), errors_are_weights(false), sigma(sig)
16  {
17  Seamstress::init_vector(numthreads, vss);
18 
19  vssp = new vector<Seamstress*>();
20  for(unsigned long int i=0;i<vss.size();i++){vssp->push_back(&(vss[i]));}
21 
23 
24  nthreads = numthreads;
25 
26  pthread_mutexattr_init(&mattr);
27  pthread_mutexattr_settype(&mattr, PTHREAD_MUTEX_NORMAL);
28  pthread_mutex_init(&mutex, &mattr);
29 
30  vector<vector<double> >* tempvec = NULL;
31  thread_points.assign(nthreads, tempvec);
32  vector<double>* tempvec2 = NULL;
33  thread_data_errors.assign(nthreads, tempvec2);
34  thread_data.assign(nthreads, tempvec2);
35 
37  }
38 
39 
41  {
42  for(unsigned int i=0;i<nthreads;i++)
43  {
44  vss[i].stop();
45  }
46  delete pins;
47  delete vssp;
48 
49  for(unsigned long int i=0;i<nthreads;i++)
50  {
51  if(thread_points[i] != NULL){delete thread_points[i];}
52  if(thread_data_errors[i] != NULL){delete thread_data_errors[i];}
53  if(thread_data[i] != NULL){delete thread_data[i];}
54  }
55  }
56 
57 
58  void GaussianRegGradHessian::setPoints(const vector<vector<double> >& POINTS)
59  {
60  if(POINTS.size() == 0){return;}
61  points = &POINTS;
62  niter = points->size();
64  else{thread_tot=nthreads;}
66  }
67 
68 
70  {
71  unsigned long int w = (*((unsigned long int *)arg));
72  unsigned long int part = niter/thread_tot;
73  unsigned long int rem = niter - part*thread_tot;
74  unsigned long int start = 0;
75  unsigned long int end = 0;
76  if(w<rem)
77  {
78  start = (part+1)*w;
79  end = start + part;
80  }
81  else
82  {
83  start = (part+1)*rem;
84  start += (part*(w-rem));
85  end = start + part - 1;
86  }
87  if(thread_points[w] != NULL){delete thread_points[w];}
88  thread_points[w] = new vector<vector<double> >();
89  vector<double> tempvec;
90  tempvec.assign((*points)[0].size(), 0.);
91  for(unsigned long int i=(start);i<=(end);i++)
92  {
93  thread_points[w]->push_back(tempvec);
94  for(unsigned long int j=0;j<thread_points[w]->back().size();j++)
95  {
96  thread_points[w]->back()[j] = (*points)[i][j];
97  }
98  }
99  }
100 
101 
102  void GaussianRegGradHessian::setErrors(const vector<double>& ERROR)
103  {
104  if(ERROR.size() == 0){return;}
106  niter = data_errors->size();
108  else{thread_tot=nthreads;}
110  }
111 
112 
114  {
115  unsigned long int w = (*((unsigned long int *)arg));
116  unsigned long int part = niter/thread_tot;
117  unsigned long int rem = niter - part*thread_tot;
118  unsigned long int start = 0;
119  unsigned long int end = 0;
120  if(w<rem)
121  {
122  start = (part+1)*w;
123  end = start + part;
124  }
125  else
126  {
127  start = (part+1)*rem;
128  start += (part*(w-rem));
129  end = start + part - 1;
130  }
131 
132  if(thread_data_errors[w] != NULL){delete thread_data_errors[w];}
133  thread_data_errors[w] = new vector<double>();
134  for(unsigned long int i=(start);i<=(end);i++)
135  {
136  thread_data_errors[w]->push_back((*data_errors)[i]);
137  }
138  }
139 
140 
141  void GaussianRegGradHessian::setData(const vector<double>& DATA)
142  {
143  if(DATA.size() == 0){return;}
144  data = &DATA;
145  niter = data->size();
147  else{thread_tot=nthreads;}
149  }
150 
151 
153  {
154  unsigned long int w = (*((unsigned long int *)arg));
155  unsigned long int part = niter/thread_tot;
156  unsigned long int rem = niter - part*thread_tot;
157  unsigned long int start = 0;
158  unsigned long int end = 0;
159  if(w<rem)
160  {
161  start = (part+1)*w;
162  end = start + part;
163  }
164  else
165  {
166  start = (part+1)*rem;
167  start += (part*(w-rem));
168  end = start + part - 1;
169  }
170 
171  if(thread_data[w] != NULL){delete thread_data[w];}
172  thread_data[w] = new vector<double>();
173  for(unsigned long int i=(start);i<=(end);i++)
174  {
175  thread_data[w]->push_back((*data)[i]);
176  }
177  }
178 
179 
181  {
183  clone->setData(*data);
184  clone->setPoints(*points);
185  if(data_has_errors==true){clone->setErrors(*data_errors);}
186  return clone;
187  }
188 
189 
190  void GaussianRegGradHessian::computeCovariance(const double& val, const MatrixXd& hessian)
191  {
192  covariance = hessian.inverse();
193 // hessian.computeInverse(&covariance);
194  if(data_has_errors==false){covariance *= (val/( (double)(points->size() - npars) ) );}
195  else if(errors_are_weights==true)
196  {
197  double scale = 0.;
198  for(unsigned int i=0;i<data_errors->size();i++)
199  {
200  double temp = 1./((*data_errors)[i]);temp*=temp;
201  scale += temp;
202  }
203  covariance *= val/(scale - (double)npars);
204  }
205  }
206 
207 
209  {
210  cov = covariance;
211  }
212 
213 
214  bool GaussianRegGradHessian::calcValGradHessian(const VectorXd& x, double& val, VectorXd& grad, MatrixXd& hessian)
215  {
216  current_eval = &x;
217 
218  val=0.;
219  grad = VectorXd::Zero(npars);
220  hessian = MatrixXd::Zero(npars, npars);
221 
222  val_output = &val;
223  grad_output = &grad;
224  hessian_output = &hessian;
225  covariance = MatrixXd::Zero(npars, npars);
226  if(points == NULL){return true;}
227  niter = points->size();
229  else{thread_tot=nthreads;}
231  covariance = hessian.inverse();
232 // hessian.computeInverse(&covariance);
233  if(data_has_errors==false){covariance *= (val/( (double)(points->size() - npars) ) );}
234  else if(errors_are_weights==true)
235  {
236  double scale = 0.;
237  for(unsigned int i=0;i<data_errors->size();i++)
238  {
239  double temp = 1./((*data_errors)[i]);temp*=temp;
240  scale += temp;
241  }
242 // scale/=( (double)(points->size() - npars) );
243 // covariance *= scale;
244  covariance *= val/(scale - (double)npars);
245  }
246  return true;
247  }
248 
249 
251  {
252  unsigned long int w = (*((unsigned long int *)arg));
253  //initialize the val, grad, and hessian to be added to those of the main thread
254  double temp_val=0.;
255  VectorXd temp_grad = VectorXd::Zero(npars);
256  MatrixXd temp_hessian = MatrixXd::Zero(npars, npars);
257 
258  //initialize the val, grad, and hessian of the function with fixed parameters at the given points vector
259  double temp_val2=0.;
260  VectorXd temp_grad2 = VectorXd::Zero(npars);
261  MatrixXd temp_hessian2 = MatrixXd::Zero(npars, npars);
262 
263  VectorXd cur_grad = VectorXd::Zero(npars);
264 
265  // f = -exp(-(t-v)^2/(2*s^2))
266  // df/dx = -f*(1/s^2)*(t-v)*dt/dx
267  // d^2f/dxdy = (-1/s^2)*[ (df/dy)*(t-v)*(dt/dx) + f*(dt/dy)*(dt/dx) + f*(t-v)*(d^2t/dxdy) ]
268 
269  FunctionGradHessian* func_instance = func->Clone();
270 
271  double inv_sig2 = 1./(sigma*sigma);
272 
273 
274  for(unsigned long int i=0;i<thread_points[w]->size();i++)
275  {
276  for(unsigned int j=0;j<(*(thread_points[w]))[i].size();j++){func_instance->setFixedPar(j, (*(thread_points[w]))[i][j]);}
277  func_instance->calcValGradHessian((*current_eval), temp_val2, temp_grad2, temp_hessian2);
278  double inv_variance = 1.;
279  if(data_has_errors==true){inv_variance = 1./(*(thread_data_errors[w]))[i];inv_variance*=inv_variance;}
280 
281  double temp1 = (temp_val2 - (*(thread_data[w]))[i]);
282  double temp2 = -exp(-temp1*temp1*0.5*inv_sig2);
283 
284  temp_val += temp2*inv_variance;
285 
286  for(unsigned int j=0;j<npars;j++)
287  {
288  cur_grad(j) = -temp2*inv_sig2*temp1*temp_grad2(j)*inv_variance;
289  temp_grad(j) += cur_grad(j);
290  }
291 
292  for(unsigned int k=0;k<npars;k++)
293  {
294  for(unsigned int j=0;j<npars;j++)
295  {
296  temp_hessian(j, k) += -inv_sig2*( cur_grad(k)*temp1*temp_grad2(j) + temp2*temp_grad2(k)*temp_grad2(j) + temp2*temp1*temp_hessian2(j, k) )*inv_variance;
297  }
298  }
299  }
300 
301 
302  pthread_mutex_lock(&mutex);
303  (*val_output) += temp_val;
304  (*grad_output) += temp_grad;
305  (*hessian_output) += temp_hessian;
306  pthread_mutex_unlock(&mutex);
307 
308 
309 
310  delete func_instance;
311  }
312 }
313 
314 
315 
316 
317 
318