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