Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
Integrator.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id: Integrator.cxx 19826 2007-09-19 19:56:11Z rdm $
2 // Authors: L. Moneta, M. Slawinska 10/2007
3 
4  /**********************************************************************
5  * *
6  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
7  * *
8  * *
9  **********************************************************************/
10 
11 #include "Math/Integrator.h"
12 
13 #include "Math/Error.h"
14 #include "Math/IFunction.h"
15 #include "Math/IntegratorMultiDim.h"
16 #include "Math/VirtualIntegrator.h"
17 
18 #include "Math/AdaptiveIntegratorMultiDim.h"
19 
20 #include "Math/GaussIntegrator.h"
21 #include "Math/GaussLegendreIntegrator.h"
22 
23 #include "Math/OneDimFunctionAdapter.h"
24 
25 
26 #include "RConfigure.h"
27 
28 #include <algorithm>
29 #include <functional>
30 #include <ctype.h> // need to use c version of tolower defined here
31 
32 
33 #include <cassert>
34 
35 #ifndef MATH_NO_PLUGIN_MANAGER
36 
37 #include "TROOT.h"
38 #include "TPluginManager.h"
39 
40 #else // case no plugin manager is available
41 #ifdef R__HAS_MATHMORE
42 #include "Math/GSLIntegrator.h"
43 #include "Math/GSLMCIntegrator.h"
44 #endif
45 
46 #endif
47 
48 
49 
50 namespace ROOT {
51 namespace Math {
52 
53 IntegrationOneDim::Type IntegratorOneDim::GetType(const char *name) {
54  if (name == 0) return IntegrationOneDim::kDEFAULT;
55  std::string typeName(name);
56  std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );
57  if (typeName == "GAUSS") return IntegrationOneDim::kGAUSS;
58  if (typeName == "GAUSSLEGENDRE") return IntegrationOneDim::kLEGENDRE;
59  if (typeName == "ADAPTIVE") return IntegrationOneDim::kADAPTIVE;
60  if (typeName == "ADAPTIVESINGULAR") return IntegrationOneDim::kADAPTIVESINGULAR;
61  if (typeName == "NONADAPTIVE") return IntegrationOneDim::kNONADAPTIVE;
62  if (!typeName.empty()) MATH_WARN_MSG("IntegratorOneDim::GetType","Invalid type name specified - use default integrator" );
63  return IntegrationOneDim::kDEFAULT;
64 }
65 
66 std::string IntegratorOneDim::GetName(IntegrationOneDim::Type type) {
67  if (type == IntegrationOneDim::kDEFAULT) type = GetType(IntegratorOneDimOptions::DefaultIntegrator().c_str() );
68  if (type == IntegrationOneDim::kGAUSS) return "Gauss";
69  if (type == IntegrationOneDim::kLEGENDRE) return "GaussLegendre";
70  if (type == IntegrationOneDim::kADAPTIVE) return "Adaptive";
71  if (type == IntegrationOneDim::kADAPTIVESINGULAR) return "AdaptiveSingular";
72  if (type == IntegrationOneDim::kNONADAPTIVE) return "NonAdaptive";
73  MATH_WARN_MSG("IntegratorOneDim::GetType","Invalid type specified " );
74  return std::string("undefined");
75 }
76 
77 
78 IntegrationMultiDim::Type IntegratorMultiDim::GetType(const char *name) {
79  if (name == 0) return IntegrationMultiDim::kDEFAULT;
80  std::string typeName(name);
81  std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );
82  if (typeName == "ADAPTIVE") return IntegrationMultiDim::kADAPTIVE;
83  if (typeName == "VEGAS") return IntegrationMultiDim::kVEGAS;
84  if (typeName == "MISER") return IntegrationMultiDim::kMISER;
85  if (typeName == "PLAIN") return IntegrationMultiDim::kPLAIN;
86  if (!typeName.empty()) MATH_WARN_MSG("IntegratorMultiDim::GetType","Invalid type name specified - use default integrator " );
87  return IntegrationMultiDim::kDEFAULT;
88 }
89 
90 std::string IntegratorMultiDim::GetName(IntegrationMultiDim::Type type) {
91  if (type == IntegrationMultiDim::kDEFAULT) type = GetType(IntegratorMultiDimOptions::DefaultIntegrator().c_str() );
92  if (type == IntegrationMultiDim::kADAPTIVE) return "ADAPTIVE";
93  if (type == IntegrationMultiDim::kVEGAS) return "VEGAS";
94  if (type == IntegrationMultiDim::kMISER) return "MISER";
95  if (type == IntegrationMultiDim::kPLAIN) return "PLAIN";
96  MATH_WARN_MSG("IntegratorMultiDim::GetType","Invalid type specified " );
97  return std::string("Undefined");
98 }
99 
100 void IntegratorOneDim::SetFunction(const IMultiGenFunction &f, unsigned int icoord , const double * x ) {
101  // set function from a multi-dim function
102  // pass also x in case of multi-dim function express the other dimensions (which are fixed)
103  unsigned int ndim = f.NDim();
104  assert (icoord < ndim);
105  ROOT::Math::OneDimMultiFunctionAdapter<> adapter(f,ndim,icoord);
106  // case I pass a vector x which is needed (for example to compute I(y) = Integral( f(x,y) dx) ) need to setCX
107  if (x != 0) adapter.SetX(x, x+ ndim);
108  SetFunction(adapter,true); // need to copy this object
109 }
110 
111 
112 // methods to create integrators
113 
114 VirtualIntegratorOneDim * IntegratorOneDim::CreateIntegrator(IntegrationOneDim::Type type , double absTol, double relTol, unsigned int size, int rule) {
115  // create the concrete class for one-dimensional integration. Use the plug-in manager if needed
116 
117  if (type == IntegrationOneDim::kDEFAULT) type = IntegratorOneDimOptions::DefaultIntegratorType();
118  if (absTol < 0) absTol = IntegratorOneDimOptions::DefaultAbsTolerance();
119  if (relTol < 0) relTol = IntegratorOneDimOptions::DefaultRelTolerance();
120  if (size <= 0) size = IntegratorOneDimOptions::DefaultWKSize();
121  if (rule <= 0) rule = IntegratorOneDimOptions::DefaultNPoints();
122  //if (ncall <= 0) ncall = IntegratorOneDimOptions::DefaultNCalls();
123 
124 
125 
126 
127 #ifndef R__HAS_MATHMORE
128  // default type is GAUSS when Mathmore is not built
129  if (type == IntegrationOneDim::kADAPTIVE ||
130  type == IntegrationOneDim::kADAPTIVESINGULAR ||
131  type == IntegrationOneDim::kNONADAPTIVE )
132  type = IntegrationOneDim::kGAUSS;
133 #endif
134 
135  if (type == IntegrationOneDim::kGAUSS)
136  return new GaussIntegrator(relTol);
137  if (type == IntegrationOneDim::kLEGENDRE) {
138  return new GaussLegendreIntegrator(rule,relTol);
139  }
140 
141  VirtualIntegratorOneDim * ig = 0;
142 
143 #ifdef MATH_NO_PLUGIN_MANAGER // no PM available
144 #ifdef R__HAS_MATHMORE
145  ig = new GSLIntegrator(type, absTol, relTol, size);
146 #else
147  MATH_ERROR_MSG("IntegratorOneDim::CreateIntegrator","Integrator type is not available in MathCore");
148 #endif
149 
150 #else // case of using Plugin Manager
151 
152 
153  {
154  R__LOCKGUARD(gROOTMutex);
155  TPluginHandler *h;
156  //gDebug = 3;
157  if ((h = gROOT->GetPluginManager()->FindHandler("ROOT::Math::VirtualIntegrator", "GSLIntegrator"))) {
158  if (h->LoadPlugin() == -1) {
159  MATH_WARN_MSG("IntegratorOneDim::CreateIntegrator","Error loading one dimensional GSL integrator - use Gauss integrator");
160  return new GaussIntegrator();
161  }
162 
163  // plugin manager requires a string
164  std::string typeName = GetName(type);
165 
166  ig = reinterpret_cast<ROOT::Math::VirtualIntegratorOneDim *>( h->ExecPlugin(5,typeName.c_str(), rule, absTol, relTol, size ) );
167  assert(ig != 0);
168  }
169 #ifdef DEBUG
170  std::cout << "Loaded Integrator " << typeid(*ig).name() << std::endl;
171 #endif
172  }
173 #endif
174 
175  return ig;
176 }
177 
178 VirtualIntegratorMultiDim * IntegratorMultiDim::CreateIntegrator(IntegrationMultiDim::Type type , double absTol, double relTol, unsigned int ncall) {
179  // create concrete class for multidimensional integration
180 
181 #ifndef R__HAS_MATHMORE
182  // when Mathmore is not built only possible type is ADAPTIVE. There is no other choice
183  type = IntegrationMultiDim::kADAPTIVE;
184 #endif
185 
186  if (type == IntegrationMultiDim::kDEFAULT) type = GetType(IntegratorMultiDimOptions::DefaultIntegrator().c_str());
187  if (absTol < 0) absTol = IntegratorMultiDimOptions::DefaultAbsTolerance();
188  if (relTol < 0) relTol = IntegratorMultiDimOptions::DefaultRelTolerance();
189  if (ncall <= 0) ncall = IntegratorMultiDimOptions::DefaultNCalls();
190  unsigned int size = IntegratorMultiDimOptions::DefaultWKSize();
191 
192 
193  // no need for PM in the adaptive case using Genz method (class is in MathCore)
194  if (type == IntegrationMultiDim::kADAPTIVE)
195  return new AdaptiveIntegratorMultiDim(absTol, relTol, ncall, size);
196 
197  // use now plugin-manager for creating the GSL integrator
198 
199  VirtualIntegratorMultiDim * ig = 0;
200 
201 #ifdef MATH_NO_PLUGIN_MANAGER // no PM available
202 #ifdef R__HAS_MATHMORE
203  ig = new GSLMCIntegrator(type, absTol, relTol, ncall);
204 #else
205  MATH_ERROR_MSG("IntegratorMultiDim::CreateIntegrator","Integrator type is not available in MathCore");
206 #endif
207 
208 #else // use ROOT Plugin-Manager to instantiate GSLMCIntegrator
209 
210  {
211  R__LOCKGUARD(gROOTMutex);
212  const char * pluginName = "GSLMCIntegrator";
213  TPluginHandler *h = nullptr;
214  //gDebug = 3;
215  if ((h = gROOT->GetPluginManager()->FindHandler("ROOT::Math::VirtualIntegrator", pluginName))) {
216  if (h->LoadPlugin() == -1) {
217  MATH_WARN_MSG("IntegratorMultiDim::CreateIntegrator","Error loading GSL MC multidim integrator - use adaptive method");
218  return new AdaptiveIntegratorMultiDim(absTol, relTol, ncall);
219  }
220 
221  std::string typeName = GetName(type);
222 
223  ig = reinterpret_cast<ROOT::Math::VirtualIntegratorMultiDim *>( h->ExecPlugin(4,typeName.c_str(), absTol, relTol, ncall ) );
224  assert(ig != 0);
225 
226 #ifdef DEBUG
227  std::cout << "Loaded Integrator " << typeid(*ig).name() << std::endl;
228 #endif
229  }
230  }
231 #endif
232  return ig;
233 }
234 
235 
236 
237 
238 } // namespace Math
239 } // namespace ROOT