Logo ROOT   6.08/07
Reference Guide
FitUtilParallel.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Author: L. Moneta Tue Nov 28 10:52:47 2006
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7  * *
8  * *
9  **********************************************************************/
10 
11 // Implementation file for class FitUtil
12 
13 #ifdef ROOT_FIT_PARALLEL
14 
15 #include "Fit/FitUtilParallel.h"
16 
17 #include "Fit/BinData.h"
18 #include "Fit/UnBinData.h"
19 #include "Fit/FitUtil.h"
20 
21 #include "Math/IParamFunction.h"
22 
23 int ncalls = 0;
24 
25 #include <limits>
26 #include <cmath>
27 #include <numeric>
28 
29 //#define DEBUG
30 #ifdef DEBUG
31 #include <iostream>
32 #endif
33 
34 #ifdef USE_PTHREAD
35 
36 #include <pthread.h>
37 
38 
39 #define NUMBER_OF_THREADS 2
40 
41 #else
42 #include <omp.h>
43 // maximum number of threads for the array
44 #define MAX_THREAD 8
45 
46 #endif
47 
48 namespace ROOT {
49 
50  namespace Fit {
51 
52  namespace FitUtilParallel {
53 
54 #ifdef USE_PTHREAD
55 
56 class ThreadData {
57 
58 public:
59 
60  ThreadData() :
61  fBegin(0),fEnd(0) ,
62  fData(0),
63  fFunc(0)
64  {}
65 
66  ThreadData(unsigned int i1, unsigned int i2, const BinData & data, IModelFunction &func) :
67  fBegin(i1), fEnd(i2),
68  fData(&data),
69  fFunc(&func)
70  {}
71 
72  void Set(unsigned int nrej, double sum) {
73  fNRej = nrej;
74  fSum = sum;
75  }
76 
77  const BinData & Data() const { return *fData; }
78 
79  IModelFunction & Func() { return *fFunc; }
80 
81  double Sum() const { return fSum; }
82 
83  unsigned int NRej() const { return fNRej; }
84 
85  unsigned int Begin() const { return fBegin; }
86  unsigned int End() const { return fEnd; }
87 
88 private:
89 
90  const unsigned int fBegin;
91  const unsigned int fEnd;
92  const BinData * fData;
93  IModelFunction * fFunc;
94  double fSum;
95  unsigned int fNRej;
96 };
97 
98 
99 // function used by the threads
100 void *EvaluateResidual(void * ptr) {
101 
102  ThreadData * t = (ThreadData *) ptr;
103 
104  unsigned int istart = t->Begin();
105  unsigned int iend = t->End();
106  double chi2 = 0;
107  unsigned int nRejected = 0;
108  const int nthreads = NUMBER_OF_THREADS;
109 
110  const BinData & data = t->Data();
111  IModelFunction & func = t->Func();
112  for (unsigned int i = istart; i < iend; i+=nthreads) {
113  const double * x = data.Coords(i);
114  double y = data.Value(i);
115  double invError = data.InvError(i);
116  double fval = 0;
117  fval = func ( x );
118 
119 // #ifdef DEBUG
120 // std::cout << x[0] << " " << y << " " << 1./invError << " params : ";
121 // for (int ipar = 0; ipar < func.NPar(); ++ipar)
122 // std::cout << p[ipar] << "\t";
123 // std::cout << "\tfval = " << fval << std::endl;
124 // #endif
125 
126  // avoid singularity in the function (infinity and nan ) in the chi2 sum
127  // eventually add possibility of excluding some points (like singularity)
128  if (fval > - std::numeric_limits<double>::max() && fval < std::numeric_limits<double>::max() ) {
129  // calculat chi2 point
130  double tmp = ( y -fval )* invError;
131  chi2 += tmp*tmp;
132  }
133  else
134  nRejected++;
135 
136  }
137 
138 #ifdef DEBUG
139  std::cout << "end loop " << istart << " " << iend << " chi2 = " << chi2 << " nrej " << nRejected << std::endl;
140 #endif
141  t->Set(nRejected,chi2);
142  return 0;
143 }
144 
145 double EvaluateChi2(IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints) {
146  // evaluate the chi2 given a function reference , the data and returns the value and also in nPoints
147  // the actual number of used points
148 
149  const int nthreads = NUMBER_OF_THREADS;
150 
151  unsigned int n = data.Size();
152 
153 #ifdef DEBUG
154  std::cout << "\n\nFit data size = " << n << std::endl;
155 #endif
156 
157  func.SetParameters(p);
158 
159  // start the threads
160  pthread_t thread[nthreads];
161  ThreadData * td[nthreads];
162  unsigned int istart = 0;
163  for (int ithread = 0; ithread < nthreads; ++ithread) {
164 // int n_th = n/nthreads;
165 // if (ithread == 0 ) n_th += n%nthreads;
166 // int iend = istart + n_th;
167  int iend = n;
168  istart = ithread;
169  td[ithread] = new ThreadData(istart,iend,data,func);
170  pthread_create(&thread[ithread], NULL, EvaluateResidual, td[ithread]);
171  //istart = iend;
172  }
173 
174  for (int ithread = 0; ithread < nthreads; ++ithread)
175  pthread_join(thread[ithread], NULL);
176 
177  // sum finally the results of the various threads
178 
179  double chi2 = 0;
180  int nRejected = 0;
181  for (int ithread = 0; ithread < nthreads; ++ithread) {
182  nRejected += td[ithread]->NRej();
183  chi2 += td[ithread]->Sum();
184  delete td[ithread];
185  }
186 
187 
188 #ifdef DEBUG
189  std::cout << "chi2 = " << chi2 << " n = " << nRejected << std::endl;
190 #endif
191 
192  nPoints = n - nRejected;
193  return chi2;
194 
195 }
196 
197 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
198 // use open MP
199 #else
200 
201 double EvaluateChi2(IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints) {
202  return FitUtil::EvaluateChi2(func,data,p,nPoints);
203 }
204 
205 
206 // use openMP (log-likelihood calculation)
207 
208 inline double EvalLogF(double fval) {
209  // evaluate the log with a protections against negative argument to the log
210  // smooth linear extrapolation below function values smaller than epsilon
211  // (better than a simple cut-off)
212  const static double epsilon = 2.*std::numeric_limits<double>::min();
213  if(fval<= epsilon)
214  return fval/epsilon + std::log(epsilon) - 1;
215  else
216  return std::log(fval);
217 }
218 
219 
220 double EvaluateLogL(IModelFunction & func, const UnBinData & data, const double * p, unsigned int &nPoints) {
221  // evaluate the LogLikelihood
222 
223  unsigned int n = data.Size();
224 
225 #ifdef DEBUG
226  std::cout << "\n\nFit data size = " << n << std::endl;
227  //std::cout << "func pointer is " << typeid(func).name() << std::endl;
228  std::cout << "\tpar = [ " << func.NPar() << " ] = ";
229  for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
230  std::cout << p[ipar] << ", ";
231  std::cout <<"---------------------------\n";
232 #endif
233 
234  double logl = 0;
235  int nRejected = 0;
236 // func.SetParameters(p);
237 
238  //std::vector<double> sum( MAX_THREAD );
239 
240 
241 #pragma omp parallel
242 #pragma omp for reduction (+:logl,nRejected)
243 //#pragma omp for reduction (+:logl,nRejected) schedule (static, 10)
244 //#pragma omp reduction (+:nRejected)
245 
246  for (unsigned int i = 0; i < n; ++ i) {
247 
248  //int ith = omp_get_thread_num();
249  //func.SetParameters(p);
250 
251  const double * x = data.Coords(i);
252  double fval = func ( x, p ); // this is thread safe passing the params
253 
254  if (fval < 0) {
255  nRejected++; // reject points with negative pdf (cannot exist)
256  }
257  else
258  //sum[ith] += EvalLogF( fval);
259  logl += EvalLogF( fval);
260 
261 #ifdef DEBUG
262 #pragma omp critical
263  { std::cout << " ==== i = " << i << " thread " << omp_get_thread_num()
264  << "fval = " << fval << " logl = " << logl << std::endl;}
265 // std::cout << "x [ " << data.PointSize() << " ] = ";
266 // for (unsigned int j = 0; j < data.PointSize(); ++j)
267 // std::cout << x[j] << "\t";
268 #endif
269  }
270 
271  // reset the number of fitting data points
272  if (nRejected != 0) nPoints = n - nRejected;
273 
274 #ifdef DEBUG
275  ncalls++;
276  int pr = std::cout.precision(15);
277  std::cout << "ncalls " << ncalls << " Logl = " << logl << " np = " << nPoints << std::endl;
278  std::cout.precision(pr);
279  assert(ncalls<3);
280 #endif
281 
282  // logl = std::accumulate(sum.begin(), sum.end(),0. );
283 
284  double result = - logl;
285  return result;
286 }
287 
288 #endif
289 
290 
291  } // end namespace FitUtilParallel
292 
293  } // end namespace Fit
294 
295 } // end namespace ROOT
296 
297 
298 
299 
300 #endif
void Begin(Int_t type)
static long int sum(long int i)
Definition: Factory.cxx:1786
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
ROOT::Math::IParamMultiFunction Func
Double_t x[n]
Definition: legend1.C:17
std::vector< std::vector< double > > Data
double EvaluateLogL(const IModelFunction &func, const UnBinData &data, const double *x, int iWeight, bool extended, unsigned int &nPoints)
evaluate the LogL given a model function and the data at the point x.
Definition: FitUtil.cxx:893
REAL epsilon
Definition: triangle.c:617
ROOT::Math::IParamMultiFunction IModelFunction
Definition: FitUtil.h:40
double EvaluateChi2(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints)
Chi2 Functions.
Definition: FitUtil.cxx:338
Double_t Sum(const double *x, const double *p)
TFitResultPtr Fit(FitObject *h1, TF1 *f1, Foption_t &option, const ROOT::Math::MinimizerOptions &moption, const char *goption, ROOT::Fit::DataRange &range)
Definition: HFitImpl.cxx:134
Double_t y[n]
Definition: legend1.C:17
double func(double *x, double *p)
Definition: stressTF1.cxx:213
#define NULL
Definition: Rtypes.h:82
void End()
double result[121]
const Int_t n
Definition: legend1.C:16
double log(double)