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