Logo ROOT   6.10/09
Reference Guide
AdaptiveIntegratorMultiDim.cxx
Go to the documentation of this file.
1 // Implementation file for class
2 // AdaptiveIntegratorMultiDim
3 //
4 #include "Math/IFunction.h"
7 #include "Math/Error.h"
8 
9 #include <cmath>
10 #include <algorithm>
11 
12 namespace ROOT {
13 namespace Math {
14 
15 
16 
17 AdaptiveIntegratorMultiDim::AdaptiveIntegratorMultiDim(double absTol, double relTol, unsigned int maxpts, unsigned int size):
18  fDim(0),
19  fMinPts(0),
20  fMaxPts(maxpts),
21  fSize(size),
22  fAbsTol(absTol),
23  fRelTol(relTol),
24  fResult(0),
25  fError(0), fRelError(0),
26  fNEval(0),
27  fStatus(-1),
28  fFun(0)
29 {
30  // constructor - without passing a function
35 }
36 
37 AdaptiveIntegratorMultiDim::AdaptiveIntegratorMultiDim( const IMultiGenFunction &f, double absTol, double relTol, unsigned int maxpts, unsigned int size):
38  fDim(f.NDim()),
39  fMinPts(0),
40  fMaxPts(maxpts),
41  fSize(size),
42  fAbsTol(absTol),
43  fRelTol(relTol),
44  fResult(0),
45  fError(0), fRelError(0),
46  fNEval(0),
47  fStatus(-1),
48  fFun(&f)
49 {
50  // constructur passing a multi-dimensional function interface
51  // constructor - without passing a function
56 }
57 
58 
59 
60 //double AdaptiveIntegratorMultiDim::Result() const { return fIntegrator->Result(); }
61 //double AdaptiveIntegratorMultiDim::Error() const { return fIntegrator->Error(); }
62 
64 {
65  // set the integration function
66  fFun = &f;
67  fDim = f.NDim();
68 }
69 
70 void AdaptiveIntegratorMultiDim::SetRelTolerance(double relTol){ this->fRelTol = relTol; }
71 
72 
74 
75 
76 double AdaptiveIntegratorMultiDim::DoIntegral(const double* xmin, const double * xmax, bool absValue)
77 {
78  // References:
79  //
80  // 1.A.C. Genz and A.A. Malik, Remarks on algorithm 006:
81  // An adaptive algorithm for numerical integration over
82  // an N-dimensional rectangular region, J. Comput. Appl. Math. 6 (1980) 295-302.
83  // 2.A. van Doren and L. de Ridder, An adaptive algorithm for numerical
84  // integration over an n-dimensional cube, J.Comput. Appl. Math. 2 (1976) 207-217.
85 
86  //to be changed later
87  unsigned int n=fDim;
88  bool kFALSE = false;
89  bool kTRUE = true;
90 
91  double epsrel = fRelTol; //specified relative accuracy
92  double epsabs = fAbsTol; //specified relative accuracy
93  //output parameters
94  fStatus = 0; //report status
95  unsigned int nfnevl; //nr of function evaluations
96  double relerr; //an estimation of the relative accuracy of the result
97 
98 
99  double ctr[15], wth[15], wthl[15], z[15];
100 
101  static const double xl2 = 0.358568582800318073;//lambda_2
102  static const double xl4 = 0.948683298050513796;//lambda_4
103  static const double xl5 = 0.688247201611685289;//lambda_5
104  static const double w2 = 980./6561; //weights/2^n
105  static const double w4 = 200./19683;
106  static const double wp2 = 245./486;//error weights/2^n
107  static const double wp4 = 25./729;
108 
109  static const double wn1[14] = { -0.193872885230909911, -0.555606360818980835,
110  -0.876695625666819078, -1.15714067977442459, -1.39694152314179743,
111  -1.59609815576893754, -1.75461057765584494, -1.87247878880251983,
112  -1.94970278920896201, -1.98628257887517146, -1.98221815780114818,
113  -1.93750952598689219, -1.85215668343240347, -1.72615963013768225};
114 
115  static const double wn3[14] = { 0.0518213686937966768, 0.0314992633236803330,
116  0.0111771579535639891,-0.00914494741655235473,-0.0294670527866686986,
117  -0.0497891581567850424,-0.0701112635269013768, -0.0904333688970177241,
118  -0.110755474267134071, -0.131077579637250419, -0.151399685007366752,
119  -0.171721790377483099, -0.192043895747599447, -0.212366001117715794};
120 
121  static const double wn5[14] = { 0.871183254585174982e-01, 0.435591627292587508e-01,
122  0.217795813646293754e-01, 0.108897906823146873e-01, 0.544489534115734364e-02,
123  0.272244767057867193e-02, 0.136122383528933596e-02, 0.680611917644667955e-03,
124  0.340305958822333977e-03, 0.170152979411166995e-03, 0.850764897055834977e-04,
125  0.425382448527917472e-04, 0.212691224263958736e-04, 0.106345612131979372e-04};
126 
127  static const double wpn1[14] = { -1.33196159122085045, -2.29218106995884763,
128  -3.11522633744855959, -3.80109739368998611, -4.34979423868312742,
129  -4.76131687242798352, -5.03566529492455417, -5.17283950617283939,
130  -5.17283950617283939, -5.03566529492455417, -4.76131687242798352,
131  -4.34979423868312742, -3.80109739368998611, -3.11522633744855959};
132 
133  static const double wpn3[14] = { 0.0445816186556927292, -0.0240054869684499309,
134  -0.0925925925925925875, -0.161179698216735251, -0.229766803840877915,
135  -0.298353909465020564, -0.366941015089163228, -0.435528120713305891,
136  -0.504115226337448555, -0.572702331961591218, -0.641289437585733882,
137  -0.709876543209876532, -0.778463648834019195, -0.847050754458161859};
138 
139  double result = 0;
140  double abserr = 0;
141  fStatus = 3;
142  nfnevl = 0;
143  relerr = 0;
144  // does not work for 1D functions
145  if (n < 2 || n > 15) {
146  MATH_WARN_MSGVAL("AdaptiveIntegratorMultiDim::Integral","Wrong function dimension",n);
147  return 0;
148  }
149 
150  double twondm = std::pow(2.0,static_cast<int>(n));
151  //unsigned int minpts = Int_t(twondm)+ 2*n*(n+1)+1;
152 
153  unsigned int ifncls = 0;
154  bool ldv = kFALSE;
155  unsigned int irgnst = 2*n+3;
156  unsigned int irlcls = (unsigned int)(twondm) +2*n*(n+1)+1;//minimal number of nodes in n dim
157  unsigned int isbrgn = irgnst;
158  unsigned int isbrgs = irgnst;
159 
160 
161  unsigned int minpts = fMinPts;
162  unsigned int maxpts = std::max(fMaxPts, irlcls) ;//specified maximal number of function evaluations
163 
164  if (minpts < 1) minpts = irlcls;
165  if (maxpts < minpts) maxpts = 10*minpts;
166 
167  // The original agorithm expected a working space array WK of length IWK
168  // with IWK Length ( >= (2N + 3) * (1 + MAXPTS/(2**N + 2N(N + 1) + 1))/2).
169  // Here, this array is allocated dynamically
170 
171  unsigned int iwk = std::max( fSize, irgnst*(1 +maxpts/irlcls)/2 );
172  double *wk = new double[iwk+10];
173 
174  unsigned int j;
175  for (j=0; j<n; j++) {
176  ctr[j] = (xmax[j] + xmin[j])*0.5;//center of a hypercube
177  wth[j] = (xmax[j] - xmin[j])*0.5;//its width
178  }
179 
180  double rgnvol, sum1, sum2, sum3, sum4, sum5, difmax, f2, f3, dif, aresult;
181  double rgncmp=0, rgnval, rgnerr;
182 
183  unsigned int j1, k, l, m, idvaxn=0, idvax0=0, isbtmp, isbtpp;
184 
185  //InitArgs(z,fParams);
186 
187 L20:
188  rgnvol = twondm;//=2^n
189  for (j=0; j<n; j++) {
190  rgnvol *= wth[j]; //region volume
191  z[j] = ctr[j]; //temporary node
192  }
193  sum1 = (*fFun)((const double*)z);//EvalPar(z,fParams); //evaluate function
194 
195  difmax = 0;
196  sum2 = 0;
197  sum3 = 0;
198 
199  //loop over coordinates
200  for (j=0; j<n; j++) {
201  z[j] = ctr[j] - xl2*wth[j];
202  if (absValue) f2 = std::abs((*fFun)(z));
203  else f2 = (*fFun)(z);
204  z[j] = ctr[j] + xl2*wth[j];
205  if (absValue) f2 += std::abs((*fFun)(z));
206  else f2 += (*fFun)(z);
207  wthl[j] = xl4*wth[j];
208  z[j] = ctr[j] - wthl[j];
209  if (absValue) f3 = std::abs((*fFun)(z));
210  else f3 = (*fFun)(z);
211  z[j] = ctr[j] + wthl[j];
212  if (absValue) f3 += std::abs((*fFun)(z));
213  else f3 += (*fFun)(z);
214  sum2 += f2;//sum func eval with different weights separately
215  sum3 += f3;//for a given region
216  dif = std::abs(7*f2-f3-12*sum1);
217  //storing dimension with biggest error/difference (?)
218  if (dif >= difmax) {
219  difmax=dif;
220  idvaxn=j+1;
221  }
222  z[j] = ctr[j];
223  }
224 
225  sum4 = 0;
226  for (j=1;j<n;j++) {
227  j1 = j-1;
228  for (k=j;k<n;k++) {
229  for (l=0;l<2;l++) {
230  wthl[j1] = -wthl[j1];
231  z[j1] = ctr[j1] + wthl[j1];
232  for (m=0;m<2;m++) {
233  wthl[k] = -wthl[k];
234  z[k] = ctr[k] + wthl[k];
235  if (absValue) sum4 += std::abs((*fFun)(z));
236  else sum4 += (*fFun)(z);
237  }
238  }
239  z[k] = ctr[k];
240  }
241  z[j1] = ctr[j1];
242  }
243 
244  sum5 = 0;
245 
246  for (j=0;j<n;j++) {
247  wthl[j] = -xl5*wth[j];
248  z[j] = ctr[j] + wthl[j];
249  }
250 L90: //sum over end nodes ~gray codes
251  if (absValue) sum5 += std::abs((*fFun)(z));
252  else sum5 += (*fFun)(z);
253  for (j=0;j<n;j++) {
254  wthl[j] = -wthl[j];
255  z[j] = ctr[j] + wthl[j];
256  if (wthl[j] > 0) goto L90;
257  }
258 
259  rgncmp = rgnvol*(wpn1[n-2]*sum1+wp2*sum2+wpn3[n-2]*sum3+wp4*sum4);
260  rgnval = wn1[n-2]*sum1+w2*sum2+wn3[n-2]*sum3+w4*sum4+wn5[n-2]*sum5;
261  rgnval *= rgnvol;
262  // avoid difference of too small numbers
263  //rgnval = 1.0E-30;
264  //rgnerr = TMath::Max( std::abs(rgnval-rgncmp), TMath::Max(std::abs(rgncmp), std::abs(rgnval) )*4.0E-16 );
265  rgnerr = std::abs(rgnval-rgncmp);//compares estim error with expected error
266 
267  result += rgnval;
268  abserr += rgnerr;
269  ifncls += irlcls;
270  aresult = std::abs(result);
271  //if (result > 0 && aresult< 1e-100) {
272  // delete [] wk;
273  // fStatus = 0; //function is probably symmetric ==> integral is null: not an error
274  // return result;
275  //}
276 
277  //if division
278  if (ldv) {
279  L110:
280  isbtmp = 2*isbrgn;
281  if (isbtmp > isbrgs) goto L160;
282  if (isbtmp < isbrgs) {
283  isbtpp = isbtmp + irgnst;
284  if (wk[isbtmp-1] < wk[isbtpp-1]) isbtmp = isbtpp;
285  }
286  if (rgnerr >= wk[isbtmp-1]) goto L160;
287  for (k=0;k<irgnst;k++) {
288  wk[isbrgn-k-1] = wk[isbtmp-k-1];
289  }
290  isbrgn = isbtmp;
291  goto L110;
292  }
293 L140:
294  isbtmp = (isbrgn/(2*irgnst))*irgnst;
295  if (isbtmp >= irgnst && rgnerr > wk[isbtmp-1]) {
296  for (k=0;k<irgnst;k++) {
297  wk[isbrgn-k-1] = wk[isbtmp-k-1];
298  }
299  isbrgn = isbtmp;
300  goto L140;
301  }
302 
303 L160: //to divide or not
304  wk[isbrgn-1] = rgnerr;//storing value & error in last
305  wk[isbrgn-2] = rgnval;//table records
306  wk[isbrgn-3] = double(idvaxn);//coordinate with biggest error
307  for (j=0;j<n;j++) {
308  isbtmp = isbrgn-2*j-4;
309  wk[isbtmp] = ctr[j];
310  wk[isbtmp-1] = wth[j];
311  }
312  if (ldv) {//divison along chosen coordinate
313  ldv = kFALSE;
314  ctr[idvax0-1] += 2*wth[idvax0-1];
315  isbrgs += irgnst;//updating the number of nodes/regions(?)
316  isbrgn = isbrgs;
317  goto L20;
318  }
319  //if no divisions to be made..
320  relerr = abserr;
321  if (aresult != 0) relerr = abserr/aresult;
322 
323 
324  if (relerr < 1e-1 && aresult < 1e-20) fStatus = 0;
325  if (relerr < 1e-3 && aresult < 1e-10) fStatus = 0;
326  if (relerr < 1e-5 && aresult < 1e-5) fStatus = 0;
327  if (isbrgs+irgnst > iwk) fStatus = 2;
328  if (ifncls+2*irlcls > maxpts) {
329  if (sum1==0 && sum2==0 && sum3==0 && sum4==0 && sum5==0){
330  fStatus = 0;
331  result = 0;
332  }
333  else
334  fStatus = 1;
335  }
336  //..and accuracy appropriare
337  // should not use absolute tolerance especially for sharp peaks
338  if ( ( relerr < epsrel || abserr < epsabs ) && ifncls >= minpts) fStatus = 0;
339 
340 #ifdef DEBUG
341  if (ifncls >= minpts) {
342  if (relerr < epsrel ) {
343  printf("relative tol reached for value %20.10g an rel error %20.10g \n",aresult, relerr);
344  fStatus = 0; // We do not use the absolute error.
345  }
346  if (abserr < epsabs ) {
347  printf("Absolute tol reached for value %20.10g and abs error %20.10g \n",aresult, abserr);
348  fStatus = 0; // We do not use the absolute error.
349  }
350  }
351 #endif
352 
353  if (fStatus == 3) {
354  ldv = kTRUE;
355  isbrgn = irgnst;
356  abserr -= wk[isbrgn-1];
357  result -= wk[isbrgn-2];
358  idvax0 = (unsigned int)(wk[isbrgn-3]);
359  for (j=0;j<n;j++) {
360  isbtmp = isbrgn-2*j-4;
361  ctr[j] = wk[isbtmp];
362  wth[j] = wk[isbtmp-1];
363  }
364  if (idvax0 < 1) {
365  // Can happen for overflows / degenerate floats.
366  idvax0 = 1;
367  ::Error("AdaptiveIntegratorMultiDim::DoIntegral()", "Logic error: idvax0 < 1!");
368  }
369  wth[idvax0-1] = 0.5*wth[idvax0-1];
370  ctr[idvax0-1] -= wth[idvax0-1];
371  goto L20;
372  }
373  nfnevl = ifncls; //number of function evaluations performed.
374  fResult = result;
375  fError = abserr;//wk[isbrgn-1];
376  fRelError = relerr;
377  fNEval = nfnevl;
378  delete [] wk;
379 
380  return result; //an approximate value of the integral
381 }
382 
383 
384 
385 double AdaptiveIntegratorMultiDim::Integral(const IMultiGenFunction &f, const double* xmin, const double * xmax)
386 {
387  // calculate integral passing a function object
388  fFun = &f;
389  return Integral(xmin, xmax);
390 
391 }
392 
394  // return the used options
398  opt.SetNCalls(fMaxPts);
399  opt.SetWKSize(fSize);
400  opt.SetIntegrator("ADAPTIVE");
401  return opt;
402 }
403 
405 {
406  // set integration options
408  MATH_ERROR_MSG("AdaptiveIntegratorMultiDim::SetOptions","Invalid options");
409  return;
410  }
411  SetAbsTolerance( opt.AbsTolerance() );
412  SetRelTolerance( opt.RelTolerance() );
413  SetMaxPts( opt.NCalls() );
414  SetSize( opt.WKSize() );
415 }
416 
417 } // namespace Math
418 } // namespace ROOT
419 
420 
421 
double DoIntegral(const double *xmin, const double *xmax, bool absVal=false)
float xmin
Definition: THbookFile.cxx:93
unsigned int NCalls() const
maximum number of function calls
const double absTol
void SetFunction(const IMultiGenFunction &f)
set the integration function (must implement multi-dim function interface: IBaseFunctionMultiDim) ...
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
ROOT::Math::IntegratorMultiDimOptions Options() const
get the option used for the integration
void SetWKSize(unsigned int size)
set workspace size
void SetMaxPts(unsigned int n)
set max points
void SetAbsTolerance(double absTol)
set absolute tolerance
void SetNCalls(unsigned int calls)
set maximum number of function calls
double RelTolerance() const
absolute tolerance
IntegrationMultiDim::Type IntegratorType() const
type of the integrator (return the enumeration type)
void SetRelTolerance(double relTol)
set relative tolerance
double pow(double, double)
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
double Integral(const double *xmin, const double *xmax)
evaluate the integral with the previously given function between xmin[] and xmax[] ...
unsigned int WKSize() const
size of the workspace
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
Numerical multi dimensional integration options.
double AbsTolerance() const
non-static methods for retrivieng options
void SetIntegrator(const char *name)
set multi-dim integrator name
TMarker * m
Definition: textangle.C:8
double Error() const
return integration error
TLine * l
Definition: textangle.C:4
float xmax
Definition: THbookFile.cxx:93
const Bool_t kFALSE
Definition: RtypesCore.h:92
double f(double x)
#define MATH_WARN_MSGVAL(loc, str, x)
Definition: Error.h:67
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
double f2(const double *x)
double result[121]
void SetSize(unsigned int size)
set workspace size
const Bool_t kTRUE
Definition: RtypesCore.h:91
void SetOptions(const ROOT::Math::IntegratorMultiDimOptions &opt)
set the options
const Int_t n
Definition: legend1.C:16
AdaptiveIntegratorMultiDim(double absTol=0.0, double relTol=1E-9, unsigned int maxpts=100000, unsigned int size=0)
Construct given optionally tolerance (absolute and relative), maximum number of function evaluation (...
void SetRelTolerance(double tol)
set the relative tolerance
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
void SetAbsTolerance(double tol)
non-static methods for setting options