Integrator.cxx

Go to the documentation of this file.

// @(#)root/mathmore:$Id: Integrator.cxx 19826 2007-09-19 19:56:11Z rdm $
// Authors: L. Moneta, M. Slawinska 10/2007
 
 /**********************************************************************
  *                                                                    *
  * Copyright (c) 2004 ROOT Foundation,  CERN/PH-SFT                   *
  *                                                                    *
  *                                                                    *
  **********************************************************************/
 
#include "Math/Integrator.h"
 
#include "Math/Error.h"
#include "Math/IFunction.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/VirtualIntegrator.h"
 
#include "Math/AdaptiveIntegratorMultiDim.h"
 
#include "Math/GaussIntegrator.h"
#include "Math/GaussLegendreIntegrator.h"
 
#include "Math/OneDimFunctionAdapter.h"
 
 
#include "RConfigure.h"
 
#include <algorithm>
#include <functional>
#include <ctype.h>   // need to use c version of tolower defined here
 
 
#include <cassert>
 
#ifndef MATH_NO_PLUGIN_MANAGER
 
#include "TROOT.h"
#include "TPluginManager.h"
 
#else // case no plugin manager is available
#ifdef R__HAS_MATHMORE
#include "Math/GSLIntegrator.h"
#include "Math/GSLMCIntegrator.h"
#endif
 
#endif
 
 
 
namespace ROOT {
namespace Math {
 
IntegrationOneDim::Type IntegratorOneDim::GetType(const char *name) {
   if (name == 0) return IntegrationOneDim::kDEFAULT;
   std::string typeName(name);
   std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );
   if (typeName == "GAUSS") return IntegrationOneDim::kGAUSS;
   if (typeName == "GAUSSLEGENDRE") return IntegrationOneDim::kLEGENDRE;
   if (typeName == "ADAPTIVE") return IntegrationOneDim::kADAPTIVE;
   if (typeName == "ADAPTIVESINGULAR") return IntegrationOneDim::kADAPTIVESINGULAR;
   if (typeName == "NONADAPTIVE") return IntegrationOneDim::kNONADAPTIVE;
   if (!typeName.empty()) MATH_WARN_MSG("IntegratorOneDim::GetType","Invalid type name specified - use default integrator" );
   return IntegrationOneDim::kDEFAULT;
}
 
std::string IntegratorOneDim::GetName(IntegrationOneDim::Type type) {
   if (type == IntegrationOneDim::kDEFAULT) type = GetType(IntegratorOneDimOptions::DefaultIntegrator().c_str() );
   if (type == IntegrationOneDim::kGAUSS) return "Gauss";
   if (type == IntegrationOneDim::kLEGENDRE) return "GaussLegendre";
   if (type == IntegrationOneDim::kADAPTIVE) return "Adaptive";
   if (type == IntegrationOneDim::kADAPTIVESINGULAR) return "AdaptiveSingular";
   if (type == IntegrationOneDim::kNONADAPTIVE) return "NonAdaptive";
   MATH_WARN_MSG("IntegratorOneDim::GetType","Invalid type specified " );
   return std::string("undefined");
}
 
 
IntegrationMultiDim::Type IntegratorMultiDim::GetType(const char *name) {
   if (name == 0) return IntegrationMultiDim::kDEFAULT;
   std::string typeName(name);
   std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );
   if (typeName == "ADAPTIVE") return IntegrationMultiDim::kADAPTIVE;
   if (typeName == "VEGAS") return IntegrationMultiDim::kVEGAS;
   if (typeName == "MISER") return IntegrationMultiDim::kMISER;
   if (typeName == "PLAIN") return IntegrationMultiDim::kPLAIN;
   if (!typeName.empty()) MATH_WARN_MSG("IntegratorMultiDim::GetType","Invalid type name specified - use default integrator " );
   return IntegrationMultiDim::kDEFAULT;
}
 
std::string IntegratorMultiDim::GetName(IntegrationMultiDim::Type type) {
   if (type == IntegrationMultiDim::kDEFAULT) type = GetType(IntegratorMultiDimOptions::DefaultIntegrator().c_str() );
   if (type == IntegrationMultiDim::kADAPTIVE) return "ADAPTIVE";
   if (type == IntegrationMultiDim::kVEGAS) return "VEGAS";
   if (type == IntegrationMultiDim::kMISER) return "MISER";
   if (type == IntegrationMultiDim::kPLAIN) return "PLAIN";
   MATH_WARN_MSG("IntegratorMultiDim::GetType","Invalid type specified " );
   return std::string("Undefined");
}
 
void IntegratorOneDim::SetFunction(const IMultiGenFunction &f, unsigned int icoord , const double * x ) {
   // set function from a multi-dim function
   // pass also x in case of multi-dim function express the other dimensions (which are fixed)
   unsigned int ndim = f.NDim();
   assert (icoord < ndim);
   ROOT::Math::OneDimMultiFunctionAdapter<> adapter(f,ndim,icoord);
   // case I pass a vector x which is needed (for example to compute I(y) = Integral( f(x,y) dx) ) need to setCX
   if (x != 0) adapter.SetX(x, x+ ndim);
   SetFunction(adapter,true); // need to copy this object
}
 
 
// methods to create integrators
 
VirtualIntegratorOneDim * IntegratorOneDim::CreateIntegrator(IntegrationOneDim::Type type , double absTol, double relTol, unsigned int size, int rule) {
   // create the concrete class for one-dimensional integration. Use the plug-in manager if needed
 
   if (type == IntegrationOneDim::kDEFAULT) type = IntegratorOneDimOptions::DefaultIntegratorType();
   if (absTol < 0) absTol = IntegratorOneDimOptions::DefaultAbsTolerance();
   if (relTol < 0) relTol = IntegratorOneDimOptions::DefaultRelTolerance();
   if (size <= 0)  size = IntegratorOneDimOptions::DefaultWKSize();
   if (rule <= 0)  rule = IntegratorOneDimOptions::DefaultNPoints();
   //if (ncall  <= 0) ncall  = IntegratorOneDimOptions::DefaultNCalls();
 
 
 
 
#ifndef R__HAS_MATHMORE
   // default type is GAUSS when Mathmore is not built
   if (type == IntegrationOneDim::kADAPTIVE ||
       type == IntegrationOneDim::kADAPTIVESINGULAR ||
       type == IntegrationOneDim::kNONADAPTIVE )
      type = IntegrationOneDim::kGAUSS;
#endif
 
   if (type == IntegrationOneDim::kGAUSS)
      return new GaussIntegrator(relTol);
   if (type == IntegrationOneDim::kLEGENDRE) {
      return new GaussLegendreIntegrator(rule,relTol);
   }
 
   VirtualIntegratorOneDim * ig = 0;
 
#ifdef MATH_NO_PLUGIN_MANAGER    // no PM available
#ifdef R__HAS_MATHMORE
   ig =  new GSLIntegrator(type, absTol, relTol, size);
#else
   MATH_ERROR_MSG("IntegratorOneDim::CreateIntegrator","Integrator type is not available in MathCore");
#endif
 
#else  // case of using Plugin Manager
 
 
   {
      R__LOCKGUARD(gROOTMutex);
      TPluginHandler *h;
      //gDebug = 3;
      if ((h = gROOT->GetPluginManager()->FindHandler("ROOT::Math::VirtualIntegrator", "GSLIntegrator"))) {
         if (h->LoadPlugin() == -1) {
            MATH_WARN_MSG("IntegratorOneDim::CreateIntegrator","Error loading one dimensional GSL integrator - use Gauss integrator");
            return new GaussIntegrator();
         }
 
         // plugin manager requires a string
         std::string typeName = GetName(type);
 
         ig = reinterpret_cast<ROOT::Math::VirtualIntegratorOneDim *>( h->ExecPlugin(5,typeName.c_str(), rule, absTol, relTol, size ) );
         assert(ig != 0);
      }
#ifdef DEBUG
      std::cout << "Loaded Integrator " << typeid(*ig).name() << std::endl;
#endif
   }
#endif
 
   return ig;
}
 
VirtualIntegratorMultiDim * IntegratorMultiDim::CreateIntegrator(IntegrationMultiDim::Type type , double absTol, double relTol, unsigned int ncall) {
   // create concrete class for multidimensional integration
 
#ifndef R__HAS_MATHMORE
   // when Mathmore is not built only possible type is ADAPTIVE. There is no other choice
   type = IntegrationMultiDim::kADAPTIVE;
#endif
 
   if (type == IntegrationMultiDim::kDEFAULT) type = GetType(IntegratorMultiDimOptions::DefaultIntegrator().c_str());
   if (absTol < 0) absTol = IntegratorMultiDimOptions::DefaultAbsTolerance();
   if (relTol < 0) relTol = IntegratorMultiDimOptions::DefaultRelTolerance();
   if (ncall  <= 0) ncall  = IntegratorMultiDimOptions::DefaultNCalls();
   unsigned int size = IntegratorMultiDimOptions::DefaultWKSize();
 
 
   // no need for PM in the adaptive  case using Genz method (class is in MathCore)
   if (type == IntegrationMultiDim::kADAPTIVE)
      return new AdaptiveIntegratorMultiDim(absTol, relTol, ncall, size);
 
   // use now plugin-manager for creating the GSL integrator
 
   VirtualIntegratorMultiDim * ig = 0;
 
#ifdef MATH_NO_PLUGIN_MANAGER  // no PM available
#ifdef R__HAS_MATHMORE
   ig =  new GSLMCIntegrator(type, absTol, relTol, ncall);
#else
   MATH_ERROR_MSG("IntegratorMultiDim::CreateIntegrator","Integrator type is not available in MathCore");
#endif
 
#else  // use ROOT Plugin-Manager to instantiate GSLMCIntegrator
 
   {
      R__LOCKGUARD(gROOTMutex);
      const char * pluginName = "GSLMCIntegrator";
      TPluginHandler *h = nullptr;
      //gDebug = 3;
      if ((h = gROOT->GetPluginManager()->FindHandler("ROOT::Math::VirtualIntegrator", pluginName))) {
         if (h->LoadPlugin() == -1) {
            MATH_WARN_MSG("IntegratorMultiDim::CreateIntegrator","Error loading GSL MC multidim integrator - use adaptive method");
            return new AdaptiveIntegratorMultiDim(absTol, relTol, ncall);
         }
 
         std::string typeName = GetName(type);
 
         ig = reinterpret_cast<ROOT::Math::VirtualIntegratorMultiDim *>( h->ExecPlugin(4,typeName.c_str(), absTol, relTol, ncall ) );
         assert(ig != 0);
 
#ifdef DEBUG
         std::cout << "Loaded Integrator " << typeid(*ig).name() << std::endl;
#endif
      }
   }
#endif
   return ig;
}
 
 
 
 
} // namespace Math
} // namespace ROOT