Logo ROOT   6.13/01
Reference Guide
GSLMCIntegrator.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Author: Magdalena Slawinska 08/2007
3 
4  /**********************************************************************
5  * *
6  * Copyright (c) 2007 ROOT Foundation, CERN/PH-SFT *
7  * *
8  * This library is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License *
10  * as published by the Free Software Foundation; either version 2 *
11  * of the License, or (at your option) any later version. *
12  * *
13  * This library is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
16  * General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this library (see file COPYING); if not, write *
20  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21  * 330, Boston, MA 02111-1307 USA, or contact the author. *
22  * *
23  **********************************************************************/
24 //
25 // implementation file for class GSLMCIntegrator
26 // Author: Magdalena Slawinska
27 //
28 //
29 
30 #include "Math/IFunction.h"
31 #include "Math/Error.h"
32 
34 
35 #include "Math/GSLMCIntegrator.h"
36 #include "Math/GSLRndmEngines.h"
38 #include "GSLRngWrapper.h"
39 
40 #include <algorithm>
41 #include <functional>
42 #include <ctype.h> // need to use c version of tolower defined here
43 
44 
45 #include "gsl/gsl_monte_vegas.h"
46 #include "gsl/gsl_monte_miser.h"
47 #include "gsl/gsl_monte_plain.h"
48 
49 
50 
51 namespace ROOT {
52 namespace Math {
53 
54 
55 
56 // constructors
57 
58 // GSLMCIntegrator::GSLMCIntegrator():
59 // fResult(0),fError(0),fStatus(-1),
60 // fWorkspace(0),
61 // fFunction(0)
62 // {
63 // // constructor of GSL MCIntegrator.Vegas MC is set as default integration type
64 // //set random number generator
65 // fRng = new GSLRngWrapper();
66 // fRng->Allocate();
67 // // use the default options
68 // ROOT::Math::IntegratorMultiDimOptions opts; // this create the default options
69 // SetOptions(opts);
70 // }
71 
72 
73 GSLMCIntegrator::GSLMCIntegrator(MCIntegration::Type type, double absTol, double relTol, unsigned int calls):
74  fType(type),
75  fDim(0),
76  fCalls((calls > 0) ? calls : IntegratorMultiDimOptions::DefaultNCalls()),
77  fAbsTol((absTol >= 0) ? absTol : IntegratorMultiDimOptions::DefaultAbsTolerance() ),
78  fRelTol((relTol >= 0) ? relTol : IntegratorMultiDimOptions::DefaultRelTolerance() ),
79  fResult(0),fError(0),fStatus(-1),
80  fExtGen(false),
81  fWorkspace(0),
82  fFunction(0)
83 {
84  // constructor of GSL MCIntegrator using enumeration as type
85  SetType(type);
86  //set random number generator
87  fRng = new GSLRngWrapper();
88  fRng->Allocate();
89  // use the default options for the needed extra parameters
90  // use the default options for the needed extra parameters
93  if (opts != 0) SetParameters( VegasParameters(*opts) );
94  }
95  else if (fType == MCIntegration::kMISER) {
97  if (opts != 0) SetParameters( MiserParameters(*opts) );
98  }
99 
100 }
101 
102 GSLMCIntegrator::GSLMCIntegrator(const char * type, double absTol, double relTol, unsigned int calls):
103  fDim(0),
104  fCalls(calls),
105  fAbsTol(absTol),
106  fRelTol(relTol),
107  fResult(0),fError(0),fStatus(-1),
108  fExtGen(false),
109  fWorkspace(0),
110  fFunction(0)
111 {
112  // constructor of GSL MCIntegrator. Vegas MC is set as default integration type if type == 0
113  SetTypeName(type);
114 
115  //set random number generator
116  fRng = new GSLRngWrapper();
117  fRng->Allocate();
118  // use the default options for the needed extra parameters
119  if (fType == MCIntegration::kVEGAS) {
121  if (opts != 0) SetParameters( VegasParameters(*opts) );
122  }
123  else if (fType == MCIntegration::kMISER) {
125  if (opts != 0) SetParameters( MiserParameters(*opts) );
126  }
127 
128 }
129 
130 
131 
133 {
134  // delete workspace
135  if (fWorkspace) delete fWorkspace;
136  if (fRng != 0 && !fExtGen) delete fRng;
137  if (fFunction != 0) delete fFunction;
138  fRng = 0;
139 
140 }
141 
142 
143 // disable copy ctrs
144 
145 
148 {}
149 
151 
152 
153 
154 
155 
157 {
158  // method to set the a generic integration function
161  fDim = f.NDim();
162 }
163 
164 void GSLMCIntegrator::SetFunction( GSLMonteFuncPointer f, unsigned int dim, void * p )
165 {
166  // method to set the a generic integration function
169  fFunction->SetParams ( p );
170  fDim = dim;
171 }
172 
173 
174 
175 double GSLMCIntegrator::Integral(const double* a, const double* b)
176 {
177  // evaluate the Integral of a over the defined interval (a[],b[])
178  assert(fRng != 0);
179  gsl_rng* fr = fRng->Rng();
180  assert(fr != 0);
181  if (!CheckFunction()) return 0;
182 
183  // initialize by creating the right WS
184  // (if dimension and type are different than previous calculation)
185  DoInitialize();
186 
188  {
190  assert(ws != 0);
191  fStatus = gsl_monte_vegas_integrate( fFunction->GetFunc(), (double *) a, (double*) b , fDim, fCalls, fr, ws->GetWS(), &fResult, &fError);
192  }
193  else if (fType == MCIntegration::kMISER)
194  {
196  assert(ws != 0);
197  fStatus = gsl_monte_miser_integrate( fFunction->GetFunc(), (double *) a, (double *) b , fDim, fCalls, fr, ws->GetWS(), &fResult, &fError);
198  }
199  else if (fType == MCIntegration::kPLAIN)
200  {
202  assert(ws != 0);
203  fStatus = gsl_monte_plain_integrate( fFunction->GetFunc(), (double *) a, (double *) b , fDim, fCalls, fr, ws->GetWS(), &fResult, &fError);
204  }
205  /**/
206  else
207  {
208 
209  fResult = 0;
210  fError = 0;
211  fStatus = -1;
212  std::cerr << "GSLIntegrator - Error: Unknown integration type" << std::endl;
213  throw std::exception();
214  }
215 
216  return fResult;
217 
218 }
219 
220 
221 double GSLMCIntegrator::Integral(const GSLMonteFuncPointer & f, unsigned int dim, double* a, double* b, void * p )
222 {
223  // evaluate the Integral for a function f over the defined interval (a[],b[])
224  SetFunction(f,dim,p);
225  return Integral(a,b);
226 }
227 
228 
229 /* to be added later
230  double GSLMCIntegrator::Integral(GSLMonteFuncPointer f, void * p, double* a, double* b)
231  {
232 
233  }
234 
235 */
236 //MCIntegration::Type GSLMCIntegrator::MCType() const {return fType;}
237 
238 /**
239  return the Result of the last Integral calculation
240 */
241 double GSLMCIntegrator::Result() const { return fResult; }
242 
243 /**
244  return the estimate of the absolute Error of the last Integral calculation
245 */
246 double GSLMCIntegrator::Error() const { return fError; }
247 
248 /**
249  return the Error Status of the last Integral calculation
250 */
251 int GSLMCIntegrator::Status() const { return fStatus; }
252 
253 
254 // setter for control Parameters (getters are not needed so far )
255 
256 /**
257  set the desired relative Error
258 */
259 void GSLMCIntegrator::SetRelTolerance(double relTol){ this->fRelTol = relTol; }
260 
261 /**
262  set the desired absolute Error
263 */
264 void GSLMCIntegrator::SetAbsTolerance(double absTol){ this->fAbsTol = absTol; }
265 
267  // delete previous exist generator
268  if (fRng && !fExtGen) delete fRng;
269  fRng = r.Engine();
270  fExtGen = true;
271 }
272 
274 {
275  // create workspace according to the type
276  fType=type;
277  if (fWorkspace != 0) {
278  if (type == fWorkspace->Type() ) return;
279  delete fWorkspace; // delete because is a different type
280  fWorkspace = 0;
281  }
282  //create Workspace according to type
283  if (type == MCIntegration::kPLAIN) {
285  }
286  else if (type == MCIntegration::kMISER) {
288  }
289  else {
290  // default: use VEGAS
291  if (type != MCIntegration::kVEGAS) {
292  MATH_WARN_MSG("GSLMCIntegration","Invalid integration type : use Vegas as default");
294  }
296  }
297 }
298 
299 void GSLMCIntegrator::SetTypeName(const char * type)
300 {
301  // set the integration type using a string
302  std::string typeName = (type!=0) ? type : "VEGAS";
303  if (type == 0) MATH_INFO_MSG("GSLMCIntegration::SetTypeName","use default Vegas integrator method");
304  std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );
305 
306  MCIntegration::Type integType = MCIntegration::kVEGAS; // default
307 
308  if (typeName == "PLAIN") {
309  integType = MCIntegration::kPLAIN;
310  }
311  else if (typeName == "MISER") {
312  integType = MCIntegration::kMISER;
313  }
314  else if (typeName != "VEGAS") {
315  MATH_WARN_MSG("GSLMCIntegration::SetTypeName","Invalid integration type : use Vegas as default");
316  }
317 
318  // create the fWorkspace object
319  if (integType != fType) SetType(integType);
320 }
321 
322 
324 {
325  // set integration mode for VEGAS method
327  {
329  assert(ws != 0);
330  if(mode == MCIntegration::kIMPORTANCE) ws->GetWS()->mode = GSL_VEGAS_MODE_IMPORTANCE;
331  else if(mode == MCIntegration::kSTRATIFIED) ws->GetWS()->mode = GSL_VEGAS_MODE_STRATIFIED;
332  else if(mode == MCIntegration::kIMPORTANCE_ONLY) ws->GetWS()->mode = GSL_VEGAS_MODE_IMPORTANCE_ONLY;
333  }
334 
335  else std::cerr << "Mode not matching integration type";
336 }
337 
339 {
340  // set integration options
341  SetTypeName(opt.Integrator().c_str() );
342  SetAbsTolerance( opt.AbsTolerance() );
343  SetRelTolerance( opt.RelTolerance() );
344  fCalls = opt.NCalls();
345 
346  //std::cout << fType << " " << MCIntegration::kVEGAS << std::endl;
347 
348  // specific options
349  ROOT::Math::IOptions * extraOpt = opt.ExtraOptions();
350  if (extraOpt) {
351  if (fType == MCIntegration::kVEGAS ) {
352  VegasParameters p(*extraOpt);
353  SetParameters(p);
354  }
355  else if (fType == MCIntegration::kMISER ) {
356  MiserParameters p(fDim); // if possible pass dimension
357  p = (*extraOpt);
358  SetParameters(p);
359  }
360  else {
361  MATH_WARN_MSG("GSLMCIntegrator::SetOptions","Invalid options set for the chosen integration type - ignore them");
362  }
363  }
364 }
365 
366 
368 {
369  // set method parameters
371  {
373  assert(ws != 0);
374  ws->SetParameters(p);
375  }
376  else
377  MATH_ERROR_MSG("GSLIntegrator::SetParameters"," Parameters not matching integration type");
378 }
379 
381 {
382  // set method parameters
384  {
386  assert(ws != 0);
387  ws->SetParameters(p);
388  }
389  else
390  MATH_ERROR_MSG("GSLIntegrator::SetParameters"," Parameters not matching integration type");
391 }
392 
393 
395 {
396  // initialize by setting integration type
397 
398  if (fWorkspace == 0) return;
399  if (fDim == fWorkspace->NDim() && fType == fWorkspace->Type() )
400  return; // can use previously existing ws
401 
402  // otherwise clear workspace
403  fWorkspace->Clear();
404  // and create a new one
405  fWorkspace->Init(fDim);
406 }
407 
408 
409 
410 //----------- methods specific for VEGAS
411 
413 {
414  // returns the error sigma from the last iteration of the VEGAS algorithm
416  {
418  assert (ws != 0);
419  return ws->GetWS()->sigma;
420  }
421  else
422  {
423  std::cerr << "Parameter not mathcing integration type";
424  return 0;
425  }
426 
427 }
428 
429 
430 /**
431 */
433 {
434  // returns chi-squared per degree of freedom for the estimate of the integral
436  {
438  assert(ws != 0);
439  return ws->GetWS()->chisq;
440  }
441  else
442  {
443  std::cerr << "Parameter not mathcing integration type";
444  return 0;
445  }
446 }
447 
448 
449 
451 {
452  // internal method to check validity of GSL function pointer
453 
454  if (fFunction && fFunction->GetFunc() ) return true;
455  MATH_ERROR_MSG("GSLMCIntegrator::CheckFunction","Function has not been specified");
456  return false;
457 }
458 
459 const char * GSLMCIntegrator::GetTypeName() const {
460  if (fType == MCIntegration::kVEGAS) return "VEGAS";
461  if (fType == MCIntegration::kMISER) return "MISER";
462  if (fType == MCIntegration::kPLAIN) return "PLAIN";
463  return "UNDEFINED";
464 }
465 
467  IOptions * extraOpts = ExtraOptions();
471  opt.SetNCalls(fCalls);
472  opt.SetWKSize(0);
473  opt.SetIntegrator(GetTypeName() );
474  return opt;
475 }
476 
478  if (!fWorkspace) return 0;
479  return fWorkspace->Options();
480 }
481 
482 
483 } // namespace Math
484 } // namespace ROOT
485 
486 
487 
GSLMCIntegrator(MCIntegration::Type type=MCIntegration::kVEGAS, double absTol=-1, double relTol=-1, unsigned int calls=0)
constructor of GSL MCIntegrator using all the default options
double(* GSLMonteFuncPointer)(double *, size_t, void *)
unsigned int NCalls() const
maximum number of function calls
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
GSLMonteFunctionWrapper * fFunction
double Integral(const GSLMonteFuncPointer &f, unsigned int dim, double *a, double *b, void *p=0)
evaluate the Integral of a function f over the defined hypercube (a,b)
Type
enumeration specifying the integration types.
void SetWKSize(unsigned int size)
set workspace size
void SetFuncPointer(GSLMonteFuncPointer f)
int Status() const
return the Error Status of the last Integral calculation
void SetNCalls(unsigned int calls)
set maximum number of function calls
structures collecting parameters for MISER multidimensional integration
Definition: MCParameters.h:76
void SetFunction(const IMultiGenFunction &f)
method to set the a generic integration function
#define MATH_WARN_MSG(loc, str)
Definition: Error.h:47
double RelTolerance() const
absolute tolerance
structures collecting parameters for VEGAS multidimensional integration FOr implementation of default...
Definition: MCParameters.h:45
void SetRelTolerance(double relTolerance)
set the desired relative Error
double ChiSqr()
returns chi-squared per degree of freedom for the estimate of the integral in the Vegas algorithm ...
GSLRandomEngine Base class for all GSL random engines, normally user instantiate the derived classes ...
virtual ~GSLMCIntegrator()
destructor
GSLMCIntegrationWorkspace * fWorkspace
GSLRngWrapper class to wrap gsl_rng structure.
Definition: GSLRngWrapper.h:25
void SetParameters(const struct VegasParameters &p)
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
ROOT::Math::IntegratorMultiDimOptions Options() const
get the option used for the integration
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
std::string Integrator() const
name of multi-dim integrator
#define MATH_INFO_MSG(loc, str)
Definition: Error.h:44
Numerical multi dimensional integration options.
void SetFunction(const FuncType &f)
Fill gsl function structure from a C++ Function class.
Interface (abstract) class for multi numerical integration It must be implemented by the concrete Int...
double AbsTolerance() const
non-static methods for retrivieng options
void SetIntegrator(const char *name)
set multi-dim integrator name
MCIntegration::Type fType
void SetParameters(const VegasParameters &p)
set default parameters for VEGAS method
GSLRngWrapper * Engine()
internal method to return the engine Used by class like GSLMCIntegrator to set the engine ...
void SetOptions(const ROOT::Math::IntegratorMultiDimOptions &opt)
set the integration options
ROOT::Math::IOptions * ExtraOptions() const
get the specific options (for Vegas or Miser) in term of string- name
void SetType(MCIntegration::Type type)
set integration method
double Result() const
return the type of the integration used
void SetTypeName(const char *typeName)
set integration method using a name instead of an enumeration
GSLMCIntegrator & operator=(const GSLMCIntegrator &)
static ROOT::Math::IOptions * FindDefault(const char *name)
void SetGenerator(GSLRandomEngine &r)
set random number generator
const char * GetTypeName() const
return the name
double Sigma()
set parameters for PLAIN method
virtual ROOT::Math::IOptions * Options() const =0
retrieve option pointer corresponding to parameters create a new object to be managed by the user ...
Namespace for new Math classes and functions.
void SetMode(MCIntegration::Mode mode)
set integration mode for VEGAS method The possible MODE are : MCIntegration::kIMPORTANCE (default) : ...
Generic interface for defining configuration options of a numerical algorithm.
Definition: IOptions.h:30
wrapper to a multi-dim function withtout derivatives for Monte Carlo multi-dimensional integration al...
virtual MCIntegration::Type Type() const =0
void SetAbsTolerance(double absTolerance)
set the desired absolute Error
IOptions * ExtraOptions() const
return extra options
virtual void Clear()
free the workspace deleting the GSL pointer
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual bool Init(size_t dim)=0
initialize the workspace creating the GSL pointer if it is not there
void SetRelTolerance(double tol)
set the relative tolerance
double Error() const
return the estimate of the absolute Error of the last Integral calculation
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
void SetAbsTolerance(double tol)
non-static methods for setting options