Logo ROOT  
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
12namespace ROOT {
13namespace Math {
14
15
16
17AdaptiveIntegratorMultiDim::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
37AdaptiveIntegratorMultiDim::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
70void AdaptiveIntegratorMultiDim::SetRelTolerance(double relTol){ this->fRelTol = relTol; }
71
72
73void AdaptiveIntegratorMultiDim::SetAbsTolerance(double absTol){ this->fAbsTol = absTol; }
74
75
76double 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
187L20:
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 }
250L90: //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 }
293L140:
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
303L160: //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
386double 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 }
414 SetMaxPts( opt.NCalls() );
415 SetSize( opt.WKSize() );
416}
417
418} // namespace Math
419} // namespace ROOT
420
421
422
double
Definition: Converters.cxx:939
size_t fSize
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:83
#define MATH_WARN_MSGVAL(loc, txt, x)
Definition: Error.h:105
#define f(i)
Definition: RSha256.hxx:104
#define e(i)
Definition: RSha256.hxx:103
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
void SetOptions(const ROOT::Math::IntegratorMultiDimOptions &opt)
set the options
void SetAbsTolerance(double absTol)
set absolute tolerance
double DoIntegral(const double *xmin, const double *xmax, bool absVal=false)
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 SetSize(unsigned int size)
set workspace size
void SetMaxPts(unsigned int n)
set max points
void SetFunction(const IMultiGenFunction &f)
set the integration function (must implement multi-dim function interface: IBaseFunctionMultiDim)
ROOT::Math::IntegratorMultiDimOptions Options() const
get the option used for the integration
void SetRelTolerance(double relTol)
set relative tolerance
double Integral(const double *xmin, const double *xmax)
evaluate the integral with the previously given function between xmin[] and xmax[]
void SetAbsTolerance(double tol)
non-static methods for setting options
double RelTolerance() const
absolute tolerance
void SetRelTolerance(double tol)
set the relative tolerance
unsigned int WKSize() const
size of the workspace
double AbsTolerance() const
non-static methods for retrieving options
void SetWKSize(unsigned int size)
set workspace size
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
Numerical multi dimensional integration options.
void SetIntegrator(const char *name)
set multi-dim integrator name
void SetNCalls(unsigned int calls)
set maximum number of function calls
unsigned int NCalls() const
maximum number of function calls
IntegrationMultiDim::Type IntegratorType() const
type of the integrator (return the enumeration type)
@ kADAPTIVE
adaptive multi-dimensional integration
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4