root/Math/Integrator.h

0001 // @(#)root/mathmore:$Id$
0002 // Authors: L. Moneta, M. Slawinska 10/2007
0003
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2007 ROOT Foundation,  CERN/PH-SFT                   *
0007   *                                                                    *
0008   *                                                                    *
0009   **********************************************************************/
0010
0011 // Header file for class Integrator
0012 //
0013 //
0014 #ifndef ROOT_Math_Integrator
0015 #define ROOT_Math_Integrator
0016
0017 #include "Math/AllIntegrationTypes.h"
0018
0019 #include "Math/IntegratorOptions.h"
0020
0021 #include "Math/IFunction.h"
0022
0023 #include "Math/VirtualIntegrator.h"
0024
0025 #include <memory>
0026 #include <vector>
0027 #include <string>
0028
0029
0030 /**
0031 @defgroup NumAlgo Numerical Algorithms
0032
0033 Numerical Algorithm classes from the \ref MathCore and \ref MathMore libraries.
0034
0035 @ingroup MathCore
0036 @ingroup MathMore
0037
0038 */
0039
0040
0041 /**
0042
0043 @defgroup Integration Numerical Integration
0044
0045 Classes for numerical integration of functions.
0046 These classes provide algorithms for integration of one-dimensional functions, with several adaptive and non-adaptive methods
0047 and for integration of multi-dimensional function using an adaptive method or MonteCarlo Integration (GSLMCIntegrator).
0048 The basic classes ROOT::Math::IntegratorOneDim provides a common interface for the one-dimensional methods while the class
0049 ROOT::Math::IntegratorMultiDim provides the interface for the multi-dimensional ones.
0050 The methods can be configured (e.g  setting the default method with its default parameters) using the ROOT::Math::IntegratorOneDimOptions and
0051 ROOT::Math::IntegratorMultiDimOptions classes.
0052
0053 @ingroup  NumAlgo
0054
0055 */
0056
0057
0058
0059 namespace ROOT {
0060 namespace Math {
0061
0062
0063
0064
0065 //____________________________________________________________________________________________
0066 /**
0067
0068 User Class for performing numerical integration of a function in one dimension.
0069 It uses the plug-in manager to load advanced numerical integration algorithms from GSL, which reimplements the
0070 algorithms used in the QUADPACK, a numerical integration package written in Fortran.
0071
0072 Various types of adaptive and non-adaptive integration are supported. These include
0073 integration over infinite and semi-infinite ranges and singular integrals.
0074
0075 The integration type is selected using the Integration::type enumeration
0076 in the class constructor.
0077 The default type is adaptive integration with singularity
0078 (ADAPTIVESINGULAR or QAGS in the QUADPACK convention) applying a Gauss-Kronrod 21-point integration rule.
0079 In the case of ADAPTIVE type, the integration rule can also be specified via the
0080 Integration::GKRule. The default rule is 31 points.
0081
0082 In the case of integration over infinite and semi-infinite ranges, the type used is always
0083 ADAPTIVESINGULAR applying a transformation from the original interval into (0,1).
0084
0085 The ADAPTIVESINGULAR type is the most sophisticated type. When performances are
0086 important, it is then recommended to use the NONADAPTIVE type in case of smooth functions or
0087  ADAPTIVE with a lower Gauss-Kronrod rule.
0088
0089 For detailed description on GSL integration algorithms see the
0090 <A HREF="http://www.gnu.org/software/gsl/manual/gsl-ref_16.html#SEC248">GSL Manual</A>.
0091
0092
0093 @ingroup Integration
0094
0095 */
0096
0097
0098 class IntegratorOneDim {
0099
0100 public:
0101
0102    typedef IntegrationOneDim::Type Type;   // for the enumerations defining the types
0103
0104     // constructors
0105
0106
0107     /**
0108         Constructor of one dimensional Integrator, default type is adaptive
0109
0110        @param type   integration type (adaptive, non-adaptive, etc..)
0111        @param absTol desired absolute Error
0112        @param relTol desired relative Error
0113        @param size maximum number of sub-intervals
0114        @param rule  Gauss-Kronrod integration rule (only for GSL kADAPTIVE type)
0115
0116        Possible type values are : kGAUSS (simple Gauss method), kADAPTIVE (from GSL), kADAPTIVESINGULAR (from GSL), kNONADAPTIVE (from GSL)
0117        Possible rule values are kGAUS15 (rule = 1), kGAUS21( rule = 2), kGAUS31(rule =3), kGAUS41 (rule=4), kGAUS51 (rule =5), kGAUS61(rule =6)
0118        lower rules are indicated for singular functions while higher for smooth functions to get better accuracies
0119
0120        NOTE: When the default values are passed, the values used are taken from the default defined in ROOT::Math::IntegratorOneDimOptions
0121     */
0122     explicit
0123     IntegratorOneDim(IntegrationOneDim::Type type = IntegrationOneDim::kDEFAULT, double absTol = -1, double relTol = -1, unsigned int size = 0, unsigned int rule = 0) :
0124        fIntegrator(nullptr), fFunc(nullptr)
0125    {
0126       fIntegrator = CreateIntegrator(type, absTol, relTol, size, rule);
0127    }
0128
0129    /**
0130        Constructor of one dimensional Integrator passing a function interface
0131
0132        @param f      integration function (1D interface). It is copied inside
0133        @param type   integration type (adaptive, non-adaptive, etc..)
0134        @param absTol desired absolute tolerance. The algorithm will stop when either the absolute OR the relative tolerance are satisfied.
0135        @param relTol desired relative tolerance
0136        @param size maximum number of sub-intervals
0137        @param rule Gauss-Kronrod integration rule (only for GSL ADAPTIVE type)
0138
0139        NOTE: When no values are passed, the values used are taken from the default defined in ROOT::Math::IntegratorOneDimOptions
0140     */
0141    explicit
0142    IntegratorOneDim(const IGenFunction &f, IntegrationOneDim::Type type = IntegrationOneDim::kDEFAULT, double absTol = -1, double relTol = -1, unsigned int size = 0, int rule = 0) :
0143       fIntegrator(nullptr), fFunc(nullptr)
0144    {
0145       fIntegrator = CreateIntegrator(type, absTol, relTol, size, rule);
0146       SetFunction(f,true);
0147    }
0148
0149    /**
0150         Template Constructor of one dimensional Integrator passing a generic function object
0151
0152        @param f      integration function (any C++ callable object implementing operator()(double x)
0153        @param type   integration type (adaptive, non-adaptive, etc..)
0154        @param absTol desired absolute tolerance. The algorithm will stop when either the absolute OR the relative tolerance are satisfied.
0155        @param relTol desired relative tolerance
0156        @param size maximum number of sub-intervals
0157        @param rule Gauss-Kronrod integration rule (only for GSL ADAPTIVE type)
0158
0159        NOTE: When no values are passed, the values used are taken from the default defined in ROOT::Math::IntegratorOneDimOptions
0160
0161     */
0162
0163    template<class Function>
0164    explicit
0165    IntegratorOneDim(Function & f, IntegrationOneDim::Type type = IntegrationOneDim::kDEFAULT, double absTol = -1, double relTol = -1, unsigned int size = 0, int rule = 0) :
0166       fIntegrator(nullptr), fFunc(nullptr)
0167    {
0168       fIntegrator = CreateIntegrator(type, absTol, relTol, size, rule);
0169       SetFunction(f);
0170    }
0171
0172    /// destructor (will delete contained pointers)
0173    virtual ~IntegratorOneDim() {
0174       if (fIntegrator) delete fIntegrator;
0175       if (fFunc) delete fFunc;
0176    }
0177
0178    // disable copy constructor and assignment operator
0179
0180 private:
0181    IntegratorOneDim(const IntegratorOneDim &) : fIntegrator(nullptr), fFunc(nullptr) {}
0182    IntegratorOneDim & operator=(const IntegratorOneDim &) { return *this; }
0183
0184 public:
0185
0186
0187    // template methods for generic functors
0188
0189    /**
0190       method to set the a generic integration function
0191       @param f integration function. The function type must implement the assignment operator, <em>  double  operator() (  double  x ) </em>
0192
0193    */
0194
0195
0196    template<class Function>
0197    inline void SetFunction(Function & f);
0198
0199    /**
0200        set one dimensional function for 1D integration
0201     */
0202    void SetFunction  (const IGenFunction &f, bool copy = false) {
0203       if (!fIntegrator) return;
0204       if (copy) {
0205          if (fFunc) delete fFunc;
0206          fFunc = f.Clone();
0207          fIntegrator->SetFunction(*fFunc);
0208          return;
0209       }
0210       fIntegrator->SetFunction(f);
0211    }
0212
0213
0214    /**
0215       Set integration function from a multi-dim function type.
0216       Can be used in case of having 1D function implementing the generic interface
0217        @param f      integration function
0218        @param icoord index of coordinate on which the integration is performed
0219        @param x  array of the passed variables values. In case of dim=1 a 0 can be passed
0220    */
0221    void SetFunction(const IMultiGenFunction &f, unsigned int icoord , const double * x );
0222
0223     // integration methods using a function
0224
0225     /**
0226        evaluate the Integral of a function f over the defined interval (a,b)
0227        @param f integration function. The function type must be a C++ callable object implementing operator()(double x)
0228        @param a lower value of the integration interval
0229        @param b upper value of the integration interval
0230     */
0231    template<class Function>
0232    double Integral(Function & f, double a, double b);
0233
0234
0235     /**
0236        evaluate the Integral of a function f over the defined interval (a,b)
0237        @param f integration function. The function type must implement the mathlib::IGenFunction interface
0238        @param a lower value of the integration interval
0239        @param b upper value of the integration interval
0240     */
0241    double Integral(const IGenFunction & f, double a, double b) {
0242      SetFunction(f,false);
0243       return Integral(a,b);
0244    }
0245
0246
0247 //    /**
0248 //      evaluate the Integral of a function f over the infinite interval (-inf,+inf)
0249 //       @param f integration function. The function type must be a C++ callable object implementing operator()(double x)
0250 //    */
0251 //    template<class Function>
0252 //    double Integral(const Function & f);
0253
0254    /**
0255       evaluate the Integral of a function f over the infinite interval (-inf,+inf)
0256       @param f integration function. The function type must implement the mathlib::IGenFunction interface
0257    */
0258    double Integral(const IGenFunction & f) {
0259       SetFunction(f,false);
0260       return Integral();
0261    }
0262
0263
0264 //    /**
0265 //      evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
0266 //      @param f integration function. The function type must be a C++ callable object implementing operator()(double x)
0267 //      @param a lower value of the integration interval
0268 //    */
0269 //    template<class Function>
0270 //    double IntegralUp(Function & f, double a);
0271
0272    /**
0273       evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
0274       @param f integration function. The function type must implement the mathlib::IGenFunction interface
0275       @param a lower value of the integration interval
0276
0277    */
0278    double IntegralUp(const IGenFunction & f, double a ) {
0279      SetFunction(f,false);
0280       return IntegralUp(a);
0281    }
0282
0283 //    /**
0284 //      evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b)
0285 //      @param f integration function. The function type must be a C++ callable object implementing operator()(double x)
0286 //      @param b upper value of the integration interval
0287 //    */
0288 //    template<class Function>
0289 //    double IntegralLow(Function & f, double b);
0290
0291    /**
0292       evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b)
0293       @param f integration function. The function type must implement the mathlib::IGenFunction interface
0294       @param b upper value of the integration interval
0295    */
0296    double IntegralLow(const IGenFunction & f, double b ) {
0297      SetFunction(f,false);
0298       return IntegralLow(b);
0299    }
0300
0301    /**
0302       evaluate the Integral of a function f with known singular points over the defined Integral (a,b)
0303       @param f integration function. The function type must be a C++ callable object implementing operator()(double x)
0304       @param pts vector containing both the function singular points and the lower/upper edges of the interval. The vector must have as first element the lower edge of the integration Integral ( \a a) and last element the upper value.
0305
0306    */
0307    template<class Function>
0308    double Integral(Function & f, const std::vector<double> & pts );
0309
0310    /**
0311       evaluate the Integral of a function f with known singular points over the defined Integral (a,b)
0312       @param f integration function. The function type must implement the mathlib::IGenFunction interface
0313       @param pts vector containing both the function singular points and the lower/upper edges of the interval. The vector must have as first element the lower edge of the integration Integral ( \a a) and last element the upper value.
0314
0315    */
0316    double Integral(const IGenFunction & f, const std::vector<double> & pts ) {
0317      SetFunction(f,false);
0318       return Integral(pts);
0319    }
0320
0321    /**
0322       evaluate the Cauchy principal value of the integral of  a function f over the defined interval (a,b) with a singularity at c
0323       @param f integration function. The function type must be a C++ callable object implementing operator()(double x)
0324       @param a lower value of the integration interval
0325       @param b upper value of the integration interval
0326       @param c position of singularity
0327
0328    */
0329    template<class Function>
0330    double IntegralCauchy(Function & f, double a, double b, double c);
0331
0332    /**
0333       evaluate the Cauchy principal value of the integral of  a function f over the defined interval (a,b) with a singularity at c
0334        @param f integration function. The function type must implement the mathlib::IGenFunction interface
0335        @param a lower value of the integration interval
0336        @param b upper value of the integration interval
0337        @param c position of singularity
0338
0339    */
0340    double IntegralCauchy(const IGenFunction & f, double a, double b, double c) {
0341      SetFunction(f,false);
0342       return IntegralCauchy(a,b,c);
0343    }
0344
0345
0346
0347    // integration method using cached function
0348
0349    /**
0350       evaluate the Integral over the defined interval (a,b) using the function previously set with Integrator::SetFunction method
0351       @param a lower value of the integration interval
0352       @param b upper value of the integration interval
0353    */
0354
0355    double Integral(double a, double b) {
0356       return !fIntegrator ? 0 : fIntegrator->Integral(a,b);
0357    }
0358
0359
0360    /**
0361       evaluate the Integral over the infinite interval (-inf,+inf) using the function previously set with Integrator::SetFunction method.
0362    */
0363
0364    double Integral( ) {
0365       return !fIntegrator ? 0 : fIntegrator->Integral();
0366    }
0367
0368    /**
0369       evaluate the Integral of a function f over the semi-infinite interval (a,+inf) using the function previously set with Integrator::SetFunction method.
0370       @param a lower value of the integration interval
0371    */
0372    double IntegralUp(double a ) {
0373       return !fIntegrator ? 0 : fIntegrator->IntegralUp(a);
0374    }
0375
0376    /**
0377       evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b) using the function previously set with Integrator::SetFunction method.
0378       @param b upper value of the integration interval
0379    */
0380    double IntegralLow( double b ) {
0381       return !fIntegrator ? 0 : fIntegrator->IntegralLow(b);
0382    }
0383    /**
0384        define operator() for IntegralLow
0385     */
0386    double operator() (double x) {
0387       return IntegralLow(x);
0388    }
0389
0390
0391    /**
0392       evaluate the Integral over the defined interval (a,b) using the function previously set with Integrator::SetFunction method. The function has known singular points.
0393       @param pts vector containing both the function singular points and the lower/upper edges of the interval. The vector must have as first element the lower edge of the integration Integral ( \a a) and last element the upper value.
0394
0395    */
0396    double Integral( const std::vector<double> & pts) {
0397       return !fIntegrator ? 0 : fIntegrator->Integral(pts);
0398    }
0399
0400    /**
0401       evaluate the Cauchy principal value of the integral of  a function f over the defined interval (a,b) with a singularity at c
0402
0403    */
0404    double IntegralCauchy(double a, double b, double c) {
0405       return !fIntegrator ? 0 : fIntegrator->IntegralCauchy(a,b,c);
0406    }
0407
0408    /**
0409       return  the Result of the last Integral calculation
0410    */
0411    double Result() const { return !fIntegrator ? 0 : fIntegrator->Result(); }
0412
0413    /**
0414       return the estimate of the absolute Error of the last Integral calculation
0415    */
0416    double Error() const { return !fIntegrator ? 0 : fIntegrator->Error(); }
0417
0418    /**
0419       return the Error Status of the last Integral calculation
0420    */
0421    int Status() const { return !fIntegrator ? -1 : fIntegrator->Status(); }
0422
0423    /**
0424       return number of function evaluations in calculating the integral
0425       (if integrator do not implement this function returns -1)
0426    */
0427    int NEval() const { return !fIntegrator ? -1 : fIntegrator->NEval(); }
0428
0429
0430    // setter for control Parameters  (getters are not needed so far )
0431
0432    /**
0433       set the desired relative Error
0434    */
0435    void SetRelTolerance(double relTolerance) { if (fIntegrator) fIntegrator->SetRelTolerance(relTolerance); }
0436
0437
0438    /**
0439       set the desired absolute Error
0440    */
0441    void SetAbsTolerance(double absTolerance) { if (fIntegrator) fIntegrator->SetAbsTolerance(absTolerance); }
0442
0443    /**
0444       return a pointer to integrator object
0445    */
0446    VirtualIntegratorOneDim * GetIntegrator() { return fIntegrator; }
0447
0448    /**
0449        set the options
0450    */
0451    void SetOptions(const ROOT::Math::IntegratorOneDimOptions & opt) { if (fIntegrator) fIntegrator->SetOptions(opt); }
0452
0453    /**
0454        retrieve the options
0455    */
0456    ROOT::Math::IntegratorOneDimOptions Options() const { return (fIntegrator) ? fIntegrator->Options() : IntegratorOneDimOptions(); }
0457
0458    /// return name of integrator
0459    std::string Name() const { return (fIntegrator) ? Options().Integrator() : std::string(""); }
0460
0461    /// static function to get the enumeration from a string
0462    static IntegrationOneDim::Type GetType(const char * name);
0463
0464    /// static function to get a string from the enumeration
0465    static std::string GetName(IntegrationOneDim::Type);
0466
0467
0468 protected:
0469
0470    VirtualIntegratorOneDim * CreateIntegrator(IntegrationOneDim::Type type , double absTol, double relTol, unsigned int size, int rule);
0471
0472 private:
0473
0474    VirtualIntegratorOneDim * fIntegrator;   ///< pointer to integrator interface class
0475    IGenFunction            * fFunc;         ///< pointer to owned function
0476
0477 };
0478
0479
0480    typedef IntegratorOneDim Integrator;
0481
0482
0483 } // namespace Math
0484 } // namespace ROOT
0485
0486
0487
0488
0489 #include "Math/WrappedFunction.h"
0490
0491 template<class Function>
0492 void ROOT::Math::IntegratorOneDim::SetFunction( Function & f) {
0493   ::ROOT::Math::WrappedFunction<Function &> wf(f);
0494   // need to copy the wrapper function, the instance created here will be deleted after SetFunction()
0495   SetFunction(wf, true);
0496 }
0497
0498 template<class Function>
0499 double ROOT::Math::IntegratorOneDim::Integral(Function & f, double a, double b) {
0500    ::ROOT::Math::WrappedFunction< Function &> wf(f);
0501    SetFunction(wf,false); // no copy is needed in this case
0502    return Integral(a,b);
0503 }
0504
0505 // remove because can create ambiguities
0506
0507 // template<class Function>
0508 // double ROOT::Math::IntegratorOneDim::Integral(const  Function & f) {
0509 //    ROOT::Math::WrappedFunction<const Function &> wf(f);
0510 //    SetFunction(wf,false); // no copy is needed in this case
0511 //    return Integral();
0512 // }
0513
0514 // template<class Function>
0515 // double ROOT::Math::IntegratorOneDim::IntegralLow(Function  & f, double x) {
0516 //    ROOT::Math::WrappedFunction< Function &> wf(f);
0517 //    SetFunction(wf,false); // no copy is needed in this case
0518 //    return IntegralLow(x);
0519 // }
0520
0521 // template<class Function>
0522 // double ROOT::Math::IntegratorOneDim::IntegralUp(Function & f, double x) {
0523 //    ROOT::Math::WrappedFunction<Function &> wf(f);
0524 //    SetFunction(wf,false); // no copy is needed in this case
0525 //    return IntegralUp(x);
0526 // }
0527
0528 template<class Function>
0529 double ROOT::Math::IntegratorOneDim::Integral(Function & f, const std::vector<double> & pts) {
0530    ::ROOT::Math::WrappedFunction<Function &> wf(f);
0531    SetFunction(wf,false); // no copy is needed in this case
0532    return Integral(pts);
0533 }
0534
0535 template<class Function>
0536 double ROOT::Math::IntegratorOneDim::IntegralCauchy(Function & f, double a, double b, double c) {
0537    ::ROOT::Math::WrappedFunction<Function & > wf(f);
0538    SetFunction(wf,false); // no copy is needed in this case
0539    return IntegralCauchy(a,b,c);
0540 }
0541
0542
0543
0544
0545
0546 #endif /* ROOT_Math_Integrator */