Logo ROOT   6.19/01
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 
73 void AdaptiveIntegratorMultiDim::SetAbsTolerance(double absTol){ this->fAbsTol = absTol; }
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 absolute 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  MATH_ERROR_MSG("AdaptiveIntegratorMultiDim::DoIntegral()", "Logic error: idvax0 < 1!");
368  //::Error("AdaptiveIntegratorMultiDim::DoIntegral()", "Logic error: idvax0 < 1!");
369  }
370  wth[idvax0-1] = 0.5*wth[idvax0-1];
371  ctr[idvax0-1] -= wth[idvax0-1];
372  goto L20;
373  }
374  nfnevl = ifncls; //number of function evaluations performed.
375  fResult = result;
376  fError = abserr;//wk[isbrgn-1];
377  fRelError = relerr;
378  fNEval = nfnevl;
379  delete [] wk;
380 
381  return result; //an approximate value of the integral
382 }
383 
384 
385 
386 double AdaptiveIntegratorMultiDim::Integral(const IMultiGenFunction &f, const double* xmin, const double * xmax)
387 {
388  // calculate integral passing a function object
389  fFun = &f;
390  return Integral(xmin, xmax);
391 
392 }
393 
395  // return the used options
399  opt.SetNCalls(fMaxPts);
400  opt.SetWKSize(fSize);
401  opt.SetIntegrator("ADAPTIVE");
402  return opt;
403 }
404 
406 {
407  // set integration options
409  MATH_ERROR_MSG("AdaptiveIntegratorMultiDim::SetOptions","Invalid options");
410  return;
411  }
412  SetAbsTolerance( opt.AbsTolerance() );
413  SetRelTolerance( opt.RelTolerance() );
414  SetMaxPts( opt.NCalls() );
415  SetSize( opt.WKSize() );
416 }
417 
418 } // namespace Math
419 } // namespace ROOT
420 
421 
422 
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
auto * m
Definition: textangle.C:8
void SetFunction(const IMultiGenFunction &f)
set the integration function (must implement multi-dim function interface: IBaseFunctionMultiDim) ...
VSD Structures.
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
#define MATH_WARN_MSGVAL(loc, txt, x)
Definition: Error.h:104
void SetAbsTolerance(double absTol)
set absolute tolerance
void SetNCalls(unsigned int calls)
set maximum number of function calls
#define f(i)
Definition: RSha256.hxx:104
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:82
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
float xmax
Definition: THbookFile.cxx:93
const Bool_t kFALSE
Definition: RtypesCore.h:88
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
auto * l
Definition: textangle.C:4
void SetSize(unsigned int size)
set workspace size
const Bool_t kTRUE
Definition: RtypesCore.h:87
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
void SetAbsTolerance(double tol)
non-static methods for setting options