Logo ROOT   6.12/07
Reference Guide
GSLIntegrator.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: L. Moneta, A. Zsenei 08/2005
3 
4  /**********************************************************************
5  * *
6  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
7  * *
8  * This library is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License *
10  * as published by the Free Software Foundation; either version 2 *
11  * of the License, or (at your option) any later version. *
12  * *
13  * This library is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
16  * General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this library (see file COPYING); if not, write *
20  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21  * 330, Boston, MA 02111-1307 USA, or contact the author. *
22  * *
23  **********************************************************************/
24 
25 // Implementation file for class GSLIntegrator
26 //
27 // Created by: Lorenzo Moneta at Thu Nov 11 14:22:32 2004
28 //
29 // Last update: Thu Nov 11 14:22:32 2004
30 //
31 
32 #include "Math/GSLIntegrator.h"
33 
34 #include "gsl/gsl_integration.h"
35 
36 #include "Math/IFunction.h"
38 #include "GSLFunctionWrapper.h"
39 
40 // for toupper
41 #include <algorithm>
42 #include <functional>
43 #include <ctype.h> // need to use c version of tolower defined here
44 
45 
46 #include <iostream>
47 
48 
49 
50 namespace ROOT {
51 namespace Math {
52 
53 
54 
55 
56 GSLIntegrator::GSLIntegrator(const Integration::Type type , const Integration::GKRule rule, double absTol, double relTol, size_t size) :
57  fType(type),
58  fRule(rule),
59  fAbsTol(absTol),
60  fRelTol(relTol),
61  fSize(size),
62  fMaxIntervals(size),
63  fResult(0),fError(0),fStatus(-1),fNEval(-1),
64  fFunction(0),
65  fWorkspace(0)
66 {
67  // constructor for all types of integrations
68  // allocate workspace (only if not adaptive algorithm)
69  if (type != Integration::kNONADAPTIVE)
71 
72 
73 }
74 
75 
76 
77 GSLIntegrator::GSLIntegrator(double absTol, double relTol, size_t size) :
80  fAbsTol(absTol),
81  fRelTol(relTol),
82  fSize(size),
83  fMaxIntervals(size),
84  fResult(0),fError(0),fStatus(-1),fNEval(-1),
85  fFunction(0),
86  fWorkspace(0)
87 {
88  // constructor with default type (ADaptiveSingular) , rule is not needed
90 
91 }
92 
93 
94 
95 GSLIntegrator::GSLIntegrator(const Integration::Type type , double absTol, double relTol, size_t size) :
96  fType(type),
98  fAbsTol(absTol),
99  fRelTol(relTol),
100  fSize(size),
101  fMaxIntervals(size),
102  fResult(0),fError(0),fStatus(-1),fNEval(-1),
103  fFunction(0),
104  fWorkspace(0)
105 {
106 
107  // constructor with default rule (gauss31) passing the type
108  // allocate workspace (only if not adaptive algorithm)
109  if (type != Integration::kNONADAPTIVE)
111 
112 }
113 
114  GSLIntegrator::GSLIntegrator(const char * type , int rule, double absTol, double relTol, size_t size) :
116  fAbsTol(absTol),
117  fRelTol(relTol),
118  fSize(size),
119  fMaxIntervals(size),
120  fResult(0),fError(0),fStatus(-1),fNEval(-1),
121  fFunction(0),
122  fWorkspace(0)
123 {
124  //std::cout << type << std::endl;
125 
127  if (type != 0) { // use this dafault
128  std::string typeName(type);
129  std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );
130  if (typeName == "NONADAPTIVE")
132  else if (typeName == "ADAPTIVE")
134  else {
135  if (typeName != "ADAPTIVESINGULAR")
136  MATH_WARN_MSG("GSLIntegrator","Use default type: AdaptiveSingular");
137  }
138  }
139 
140 
141  // constructor with default rule (gauss31) passing the type
142  // allocate workspace (only if not adaptive algorithm)
145 
147 
148 }
149 
150 
152 {
153  // delete workspace and function
154  if (fFunction) delete fFunction;
155  if (fWorkspace) delete fWorkspace;
156 }
157 
160 {
161  // dummy copy ctr
162 }
163 
165 {
166  // dummy operator=
167  if (this == &rhs) return *this; // time saving self-test
168 
169  return *this;
170 }
171 
172 
173 
174 
176  // fill GSLFunctionWrapper with the pointer to the function
177  if (fFunction ==0) fFunction = new GSLFunctionWrapper();
178  fFunction->SetFuncPointer( fp );
179  fFunction->SetParams ( p );
180 }
181 
183  // set function (make a copy of it)
184  if (fFunction ==0) fFunction = new GSLFunctionWrapper();
186 }
187 
188 // evaluation methods
189 
190 double GSLIntegrator::Integral(double a, double b) {
191  // defined integral evaluation
192  // need here look at all types of algorithms
193  // find more elegant solution ? Use template OK, but need to chose algorithm statically , t.b.i.
194 
195  if (!CheckFunction()) return 0;
196 
198  size_t neval = 0; // need to export this ?
200  fNEval = neval;
201  }
202  else if (fType == Integration::kADAPTIVE) {
204  const int npts[6] = {15,21,31,41,51,61};
205  assert(fRule>=1 && fRule <=6);
206  fNEval = (fWorkspace->GetWS()->size)*npts[fRule-1]; // get size of workspace (number of iterations)
207  }
209 
210  // singular integration - look if we know about singular points
211 
212 
214  fNEval = (fWorkspace->GetWS()->size) * 21; //since 21 point rule is used in qags
215  }
216  else {
217 
218  fResult = 0;
219  fError = 0;
220  fStatus = -1;
221  std::cerr << "GSLIntegrator - Error: Unknown integration type" << std::endl;
222  throw std::exception(); //"Unknown integration type");
223  }
224 
225  return fResult;
226 
227 }
228 
229 //=============================
230 double GSLIntegrator::IntegralCauchy(double a, double b, double c) {
231  //eval integral with Cauchy principal value defined at the value c
232  if (!CheckFunction()) return 0;
233 
234  fStatus = gsl_integration_qawc( fFunction->GetFunc(), a, b , c, fAbsTol, fRelTol, fMaxIntervals, fWorkspace->GetWS(), &fResult, &fError);
235  fNEval = (fWorkspace->GetWS()->size) * 15; // 15 point rule is used ?
236 
237  return fResult;
238 
239 }
240 
241 double GSLIntegrator::IntegralCauchy(const IGenFunction & f, double a, double b, double c) {
242  //eval integral with Cauchy principal value defined at the value c
243 
244  if (!CheckFunction()) return 0;
245  SetFunction(f);
246  return IntegralCauchy(a, b, c);
247 
248 }
249 
250 //==============================
251 
252 double GSLIntegrator::Integral( const std::vector<double> & pts) {
253  // integral eval with singularities
254 
255  if (!CheckFunction()) return 0;
256 
257  if (fType == Integration::kADAPTIVESINGULAR && pts.size() >= 2 ) {
258  // remove constness ( should be const in GSL ? )
259  double * p = const_cast<double *>(&pts.front() );
260  fStatus = gsl_integration_qagp( fFunction->GetFunc(), p, pts.size() , fAbsTol, fRelTol, fMaxIntervals, fWorkspace->GetWS(), &fResult, &fError);
261  fNEval = (fWorkspace->GetWS()->size) * 15; // 15 point rule is used ?
262  }
263  else {
264  fResult = 0;
265  fError = 0;
266  fStatus = -1;
267  std::cerr << "GSLIntegrator - Error: Unknown integration type or not enough singular points defined" << std::endl;
268  return 0;
269  }
270  return fResult;
271 }
272 
273 
275  // Eval for indefined integrals: use QAGI method
276  // if method was chosen NO_ADAPTIVE WS does not exist create it
277 
278  if (!CheckFunction()) return 0;
279 
281 
283  fNEval = (fWorkspace->GetWS()->size) * 15; // 15 point rule is used ?
284 
285  return fResult;
286 }
287 
288 
289 
290 double GSLIntegrator::IntegralUp( double a ) {
291  // Integral between [a, + inf]
292  // if method was chosen NO_ADAPTIVE WS does not exist create it
293 
294  if (!CheckFunction()) return 0;
295 
297 
299  fNEval = (fWorkspace->GetWS()->size) * 21; // 21 point rule is used ?
300 
301  return fResult;
302 }
303 
304 
305 
306 double GSLIntegrator::IntegralLow( double b ) {
307  // Integral between [-inf, + b]
308  // if method was chosen NO_ADAPTIVE WS does not exist create it
309 
310  if (!CheckFunction()) return 0;
311 
313 
315  fNEval = (fWorkspace->GetWS()->size) * 21; // 21 point rule is used ?
316 
317  return fResult;
318 }
319 
320 
321 
322 
323 double GSLIntegrator::Integral(const IGenFunction & f, double a, double b) {
324  // use generic function interface
325  SetFunction(f);
326  return Integral(a,b);
327 }
328 
330  // use generic function interface
331  SetFunction(f);
332  return Integral();
333 }
334 
335 double GSLIntegrator::IntegralUp(const IGenFunction & f, double a) {
336  // use generic function interface
337  SetFunction(f);
338  return IntegralUp(a);
339 }
340 
341 double GSLIntegrator::IntegralLow(const IGenFunction & f, double b) {
342  // use generic function interface
343  SetFunction(f);
344  return IntegralLow(b);
345 }
346 
347 double GSLIntegrator::Integral(const IGenFunction & f, const std::vector<double> & pts) {
348  // use generic function interface
349  SetFunction(f);
350  return Integral(pts);
351 }
352 
353 
354 
355 
356 double GSLIntegrator::Integral( GSLFuncPointer f , void * p, double a, double b) {
357  // use c free function pointer
358  SetFunction( f, p);
359  return Integral(a,b);
360 }
361 
363  // use c free function pointer
364  SetFunction( f, p);
365  return Integral();
366 }
367 
368 double GSLIntegrator::IntegralUp( GSLFuncPointer f, void * p, double a ) {
369  // use c free function pointer
370  SetFunction( f, p);
371  return IntegralUp(a);
372 }
373 
374 double GSLIntegrator::IntegralLow( GSLFuncPointer f, void * p, double b ) {
375  // use c free function pointer
376  SetFunction( f, p);
377  return IntegralLow(b);
378 }
379 
380 double GSLIntegrator::Integral( GSLFuncPointer f, void * p, const std::vector<double> & pts ) {
381  // use c free function pointer
382  SetFunction( f, p);
383  return Integral(pts);
384 }
385 
386 
387 
388 double GSLIntegrator::Result() const { return fResult; }
389 
390 double GSLIntegrator::Error() const { return fError; }
391 
392 int GSLIntegrator::Status() const { return fStatus; }
393 
394 
395 // get and setter methods
396 
397 // double GSLIntegrator::getAbsTolerance() const { return fAbsTol; }
398 
399 void GSLIntegrator::SetAbsTolerance(double absTol){ this->fAbsTol = absTol; }
400 
401 // double GSLIntegrator::getRelTolerance() const { return fRelTol; }
402 
403 void GSLIntegrator::SetRelTolerance(double relTol){ this->fRelTol = relTol; }
404 
405 
407 
409  // check if a function has been previously set.
410  if (fFunction && fFunction->IsValid()) return true;
411  fStatus = -1; fResult = 0; fError = 0;
412  std::cerr << "GSLIntegrator - Error : Function has not been specified " << std::endl;
413  return false;
414 }
415 
417 {
418  // set integration options
419  fType = opt.IntegratorType();
424  MATH_WARN_MSG("GSLIntegrator::SetOptions","Invalid rule options - use default ADAPTIVESINGULAR");
426  }
427  SetAbsTolerance( opt.AbsTolerance() );
428  SetRelTolerance( opt.RelTolerance() );
429  fSize = opt.WKSize();
431  if (fType == Integration::kADAPTIVE) {
432  int npts = opt.NPoints();
433  if ( npts >= Integration::kGAUSS15 && npts <= Integration::kGAUSS61)
434  fRule = (Integration::GKRule) npts;
435  else {
436  MATH_WARN_MSG("GSLIntegrator::SetOptions","Invalid rule options - use default GAUSS31");
438  }
439  }
440 }
441 
446  opt.SetWKSize(fSize);
447  opt.SetIntegrator(GetTypeName() );
448 
450  opt.SetNPoints(fRule);
452  opt.SetNPoints( Integration::kGAUSS31 ); // fixed rule for adaptive singular
453  else
454  opt.SetNPoints( 0 ); // not available for the rest
455 
456  return opt;
457 }
458 
459 const char * GSLIntegrator::GetTypeName() const {
460  if (fType == IntegrationOneDim::kADAPTIVE) return "Adaptive";
461  if (fType == IntegrationOneDim::kADAPTIVESINGULAR) return "AdaptiveSingular";
462  if (fType == IntegrationOneDim::kNONADAPTIVE) return "NonAdaptive";
463  return "Undefined";
464 }
465 
466 
467 } // namespace Math
468 } // namespace ROOT
double Result() const
return the Result of the last Integral calculation
GKRule
enumeration specifying the Gauss-KronRod integration rule for ADAPTIVE integration type ...
Interface (abstract) class for 1D numerical integration It must be implemented by the concrate Integr...
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
int gsl_integration_qng(const gsl_function *f, double a, double b, double epsabs, double epsrel, double *result, double *abserr, size_t *neval)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
int gsl_integration_qagiu(gsl_function *f, double a, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
double IntegralUp(const IGenFunction &f, double a)
evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
void SetWKSize(unsigned int size)
set workspace size
void SetIntegrationRule(Integration::GKRule)
set the integration rule (Gauss-Kronrod rule).
void SetNPoints(unsigned int n)
set number of points rule values of 1,2,3,4,5,6 corresponds to 15,21,31,41,51,61 and they are used in...
void SetFunction(const FuncType &f)
fill the GSL C struct from a generic C++ callable object implementing operator()
int Status() const
return the Error Status of the last Integral calculation
#define MATH_WARN_MSG(loc, str)
Definition: Error.h:47
int gsl_integration_qag(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, int key, gsl_integration_workspace *workspace, double *result, double *abserr)
double RelTolerance() const
absolute tolerance
GSLFunctionWrapper * fFunction
virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions &opt)
set the options
int gsl_integration_qags(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
Class for performing numerical integration of a function in one dimension.
Definition: GSLIntegrator.h:90
IntegrationOneDim::Type IntegratorType() const
type of the integrator (return the enumeration type)
unsigned int WKSize() const
size of the workspace
Integration::Type fType
double IntegralLow(const IGenFunction &f, double b)
evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b)
double(* GSLFuncPointer)(double, void *)
Function pointer corresponding to gsl_function signature.
double AbsTolerance() const
non-static methods for retrivieng options
auto * a
Definition: textangle.C:12
double Integral()
evaluate the Integral over the infinite interval (-inf,+inf) using the function previously set with G...
Numerical one dimensional integration options.
void SetFunction(const IGenFunction &f)
method to set the a generic integration function
unsigned int NPoints() const
maximum number of function calls
bool IsValid()
check if function is valid (has been set)
double Error() const
return the estimate of the absolute Error of the last Integral calculation
Type
enumeration specifying the integration types.
PyObject * fType
GSLIntegrator(double absTol=1.E-9, double relTol=1E-6, size_t size=1000)
Default constructor of GSL Integrator for Adaptive Singular integration.
void SetRelTolerance(double relTolerance)
set the desired relative Error
int type
Definition: TGX11.cxx:120
GSLIntegrationWorkspace * fWorkspace
double IntegralCauchy(double a, double b, double c)
evaluate the Cauchy principal value of the integral of a previously defined function f over the defin...
void SetFuncPointer(GSLFuncPointer f)
set in the GSL C struct the pointer to the function evaluation
Namespace for new Math classes and functions.
int gsl_integration_qagil(gsl_function *f, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
void SetAbsTolerance(double absTolerance)
set the desired absolute Error
gsl_integration_workspace * GetWS()
const char * GetTypeName() const
return the name
GSLIntegrator & operator=(const GSLIntegrator &)
Wrapper class to the gsl_function C structure.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
int gsl_integration_qagi(gsl_function *f, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
Integration::GKRule fRule
void SetParams(void *p)
set in the GSL C struct the extra-object pointer
void SetIntegrator(const char *name)
set 1D integrator name
void SetRelTolerance(double tol)
set the relative tolerance
void SetAbsTolerance(double tol)
non-static methods for setting options