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
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
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
ROOT::Math::IntegratorMultiDimOptions Options() const override
get the option used for the integration
unsigned int fMinPts
minimum number of function evaluation requested
double DoIntegral(const double *xmin, const double *xmax, bool absVal=false)
double Integral(const double *xmin, const double *xmax) override
evaluate the integral with the previously given function between xmin[] and xmax[]
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 (...
unsigned int fMaxPts
maximum number of function evaluation requested
void SetSize(unsigned int size)
set workspace size
void SetRelTolerance(double relTol) override
set relative tolerance
void SetFunction(const IMultiGenFunction &f) override
set the integration function (must implement multi-dim function interface: IBaseFunctionMultiDim)
void SetMaxPts(unsigned int n)
set max points
void SetAbsTolerance(double absTol) override
set absolute tolerance
unsigned int fSize
max size of working array (explode with dimension)
void SetOptions(const ROOT::Math::IntegratorMultiDimOptions &opt) override
set the options
int fStatus
status of algorithm (error if not zero)
unsigned int fDim
dimensionality of integrand
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
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1739
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1753
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.
TMarker m
Definition: textangle.C:8
TLine l
Definition: textangle.C:4