Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
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
51namespace ROOT {
52namespace 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
73GSLMCIntegrator::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
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
102GSLMCIntegrator::GSLMCIntegrator(const char * type, double absTol, double relTol, unsigned int calls):
103 fType(MCIntegration::kDEFAULT),
104 fDim(0),
105 fCalls(calls),
106 fAbsTol(absTol),
107 fRelTol(relTol),
108 fResult(0),fError(0),fStatus(-1),
109 fExtGen(false),
110 fWorkspace(0),
111 fFunction(0)
112{
113 // constructor of GSL MCIntegrator. Vegas MC is set as default integration type if type == 0
115
116 //set random number generator
117 fRng = new GSLRngWrapper();
118 fRng->Allocate();
119 // use the default options for the needed extra parameters
122 if (opts != 0) SetParameters( VegasParameters(*opts) );
123 }
124 else if (fType == MCIntegration::kMISER) {
126 if (opts != 0) SetParameters( MiserParameters(*opts) );
127 }
128
129}
130
131
132
134{
135 // delete workspace
136 if (fWorkspace) delete fWorkspace;
137 if (fRng != 0 && !fExtGen) delete fRng;
138 if (fFunction != 0) delete fFunction;
139 fRng = 0;
140
141}
142
143
144// disable copy ctrs
145
146
149{}
150
152
153
154
155
156
158{
159 // method to set the a generic integration function
162 fDim = f.NDim();
163}
164
165void GSLMCIntegrator::SetFunction( GSLMonteFuncPointer f, unsigned int dim, void * p )
166{
167 // method to set the a generic integration function
170 fFunction->SetParams ( p );
171 fDim = dim;
172}
173
174
175
176double GSLMCIntegrator::Integral(const double* a, const double* b)
177{
178 // evaluate the Integral of a over the defined interval (a[],b[])
179 assert(fRng != 0);
180 gsl_rng* fr = fRng->Rng();
181 assert(fr != 0);
182 if (!CheckFunction()) return 0;
183
184 // initialize by creating the right WS
185 // (if dimension and type are different than previous calculation)
186 DoInitialize();
187
189 {
191 assert(ws != 0);
192 fStatus = gsl_monte_vegas_integrate( fFunction->GetFunc(), (double *) a, (double*) b , fDim, fCalls, fr, ws->GetWS(), &fResult, &fError);
193 }
194 else if (fType == MCIntegration::kMISER)
195 {
197 assert(ws != 0);
198 fStatus = gsl_monte_miser_integrate( fFunction->GetFunc(), (double *) a, (double *) b , fDim, fCalls, fr, ws->GetWS(), &fResult, &fError);
199 }
200 else if (fType == MCIntegration::kPLAIN)
201 {
203 assert(ws != 0);
204 fStatus = gsl_monte_plain_integrate( fFunction->GetFunc(), (double *) a, (double *) b , fDim, fCalls, fr, ws->GetWS(), &fResult, &fError);
205 }
206 /**/
207 else
208 {
209
210 fResult = 0;
211 fError = 0;
212 fStatus = -1;
213 std::cerr << "GSLIntegrator - Error: Unknown integration type" << std::endl;
214 throw std::exception();
215 }
216
217 return fResult;
218
219}
220
221
222double GSLMCIntegrator::Integral(const GSLMonteFuncPointer & f, unsigned int dim, double* a, double* b, void * p )
223{
224 // evaluate the Integral for a function f over the defined interval (a[],b[])
225 SetFunction(f,dim,p);
226 return Integral(a,b);
227}
228
229
230/* to be added later
231 double GSLMCIntegrator::Integral(GSLMonteFuncPointer f, void * p, double* a, double* b)
232 {
233
234 }
235
236*/
237//MCIntegration::Type GSLMCIntegrator::MCType() const {return fType;}
238
239/**
240 return the Result of the last Integral calculation
241*/
242double GSLMCIntegrator::Result() const { return fResult; }
243
244/**
245 return the estimate of the absolute Error of the last Integral calculation
246*/
247double GSLMCIntegrator::Error() const { return fError; }
248
249/**
250 return the Error Status of the last Integral calculation
251*/
252int GSLMCIntegrator::Status() const { return fStatus; }
253
254
255// setter for control Parameters (getters are not needed so far )
256
257/**
258 set the desired relative Error
259*/
260void GSLMCIntegrator::SetRelTolerance(double relTol){ this->fRelTol = relTol; }
261
262/**
263 set the desired absolute Error
264*/
265void GSLMCIntegrator::SetAbsTolerance(double absTol){ this->fAbsTol = absTol; }
266
268 // delete previous exist generator
269 if (fRng && !fExtGen) delete fRng;
270 fRng = r.Engine();
271 fExtGen = true;
272}
273
275{
276 // create workspace according to the type
277 fType=type;
278 if (fWorkspace != 0) {
279 if (type == fWorkspace->Type() ) return;
280 delete fWorkspace; // delete because is a different type
281 fWorkspace = 0;
282 }
283 //create Workspace according to type
286 }
287 else if (type == MCIntegration::kMISER) {
289 }
290 else {
291 // default: use VEGAS
293 MATH_WARN_MSG("GSLMCIntegration","Invalid integration type : use Vegas as default");
295 }
297 }
298}
299
301{
302 // set the integration type using a string
303 std::string typeName = (type!=0) ? type : "VEGAS";
304 if (type == 0) MATH_INFO_MSG("GSLMCIntegration::SetTypeName","use default Vegas integrator method");
305 std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );
306
307 MCIntegration::Type integType = MCIntegration::kVEGAS; // default
308
309 if (typeName == "PLAIN") {
310 integType = MCIntegration::kPLAIN;
311 }
312 else if (typeName == "MISER") {
313 integType = MCIntegration::kMISER;
314 }
315 else if (typeName != "VEGAS") {
316 MATH_WARN_MSG("GSLMCIntegration::SetTypeName","Invalid integration type : use Vegas as default");
317 }
318
319 // create the fWorkspace object
320 // if it exists already with the same type it will not be re-created
321 SetType(integType);
322}
323
324
326{
327 // set integration mode for VEGAS method
329 {
331 assert(ws != 0);
332 if(mode == MCIntegration::kIMPORTANCE) ws->GetWS()->mode = GSL_VEGAS_MODE_IMPORTANCE;
333 else if(mode == MCIntegration::kSTRATIFIED) ws->GetWS()->mode = GSL_VEGAS_MODE_STRATIFIED;
334 else if(mode == MCIntegration::kIMPORTANCE_ONLY) ws->GetWS()->mode = GSL_VEGAS_MODE_IMPORTANCE_ONLY;
335 }
336
337 else std::cerr << "Mode not matching integration type";
338}
339
341{
342 // set integration options
343 SetTypeName(opt.Integrator().c_str() );
346 fCalls = opt.NCalls();
347
348 //std::cout << fType << " " << MCIntegration::kVEGAS << std::endl;
349
350 // specific options
351 ROOT::Math::IOptions * extraOpt = opt.ExtraOptions();
352 if (extraOpt) {
354 VegasParameters p(*extraOpt);
355 SetParameters(p);
356 }
357 else if (fType == MCIntegration::kMISER ) {
358 MiserParameters p(fDim); // if possible pass dimension
359 p = (*extraOpt);
360 SetParameters(p);
361 }
362 else {
363 MATH_WARN_MSG("GSLMCIntegrator::SetOptions","Invalid options set for the chosen integration type - ignore them");
364 }
365 }
366}
367
368
370{
371 // set method parameters
373 {
375 assert(ws != 0);
376 ws->SetParameters(p);
377 }
378 else
379 MATH_ERROR_MSG("GSLIntegrator::SetParameters"," Parameters not matching integration type");
380}
381
383{
384 // set method parameters
386 {
388 assert(ws != 0);
389 ws->SetParameters(p);
390 }
391 else
392 MATH_ERROR_MSG("GSLIntegrator::SetParameters"," Parameters not matching integration type");
393}
394
395
397{
398 // initialize by setting integration type
399
400 if (fWorkspace == 0) return;
401 if (fDim == fWorkspace->NDim() && fType == fWorkspace->Type() )
402 return; // can use previously existing ws
403
404 // otherwise clear workspace
405 fWorkspace->Clear();
406 // and create a new one
408}
409
410
411
412//----------- methods specific for VEGAS
413
415{
416 // returns the error sigma from the last iteration of the VEGAS algorithm
418 {
420 assert (ws != 0);
421 return ws->GetWS()->sigma;
422 }
423 else
424 {
425 std::cerr << "Parameter not mathcing integration type";
426 return 0;
427 }
428
429}
430
431
432/**
433*/
435{
436 // returns chi-squared per degree of freedom for the estimate of the integral
438 {
440 assert(ws != 0);
441 return ws->GetWS()->chisq;
442 }
443 else
444 {
445 std::cerr << "Parameter not mathcing integration type";
446 return 0;
447 }
448}
449
450
451
453{
454 // internal method to check validity of GSL function pointer
455
456 if (fFunction && fFunction->GetFunc() ) return true;
457 MATH_ERROR_MSG("GSLMCIntegrator::CheckFunction","Function has not been specified");
458 return false;
459}
460
461const char * GSLMCIntegrator::GetTypeName() const {
462 if (fType == MCIntegration::kVEGAS) return "VEGAS";
463 if (fType == MCIntegration::kMISER) return "MISER";
464 if (fType == MCIntegration::kPLAIN) return "PLAIN";
465 return "UNDEFINED";
466}
467
469 IOptions * extraOpts = ExtraOptions();
473 opt.SetNCalls(fCalls);
474 opt.SetWKSize(0);
476 return opt;
477}
478
480 if (!fWorkspace) return 0;
481 return fWorkspace->Options();
482}
483
484
485} // namespace Math
486} // namespace ROOT
487
488
489
#define MATH_INFO_MSG(loc, str)
Pre-processor macro to report messages which can be configured to use ROOT error or simply an std::io...
Definition Error.h:77
#define MATH_ERROR_MSG(loc, str)
Definition Error.h:83
#define MATH_WARN_MSG(loc, str)
Definition Error.h:80
ROOT::R::TRInterface & r
Definition Object.C:4
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define a(i)
Definition RSha256.hxx:99
int type
Definition TGX11.cxx:121
void SetAbsTolerance(double tol)
non-static methods for setting options
double RelTolerance() const
absolute tolerance
void SetRelTolerance(double tol)
set the relative tolerance
double AbsTolerance() const
non-static methods for retrieving options
void SetWKSize(unsigned int size)
set workspace size
IOptions * ExtraOptions() const
return extra options
virtual void Clear()
free the workspace deleting the GSL pointer
virtual ROOT::Math::IOptions * Options() const =0
retrieve option pointer corresponding to parameters create a new object to be managed by the user
virtual MCIntegration::Type Type() const =0
virtual bool Init(size_t dim)=0
initialize the workspace creating the GSL pointer if it is not there
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)
double Error() const
return the estimate of the absolute Error of the last Integral calculation
GSLMonteFunctionWrapper * fFunction
GSLMCIntegrator(MCIntegration::Type type=MCIntegration::kVEGAS, double absTol=-1, double relTol=-1, unsigned int calls=0)
constructor of GSL MCIntegrator.
const char * GetTypeName() const
return the name
GSLMCIntegrator & operator=(const GSLMCIntegrator &)
double Result() const
return the type of the integration used
void SetOptions(const ROOT::Math::IntegratorMultiDimOptions &opt)
set the integration options
void SetParameters(const VegasParameters &p)
set default parameters for VEGAS method
void SetGenerator(GSLRandomEngine &r)
set random number generator
ROOT::Math::IntegratorMultiDimOptions Options() const
get the option used for the integration
ROOT::Math::IOptions * ExtraOptions() const
get the specific options (for Vegas or Miser) in term of string- name
GSLMCIntegrationWorkspace * fWorkspace
void SetMode(MCIntegration::Mode mode)
set integration mode for VEGAS method The possible MODE are : MCIntegration::kIMPORTANCE (default) : ...
void SetTypeName(const char *typeName)
set integration method using a name instead of an enumeration
void SetRelTolerance(double relTolerance)
set the desired relative Error
virtual ~GSLMCIntegrator()
destructor
double Sigma()
set parameters for PLAIN method
double ChiSqr()
returns chi-squared per degree of freedom for the estimate of the integral in the Vegas algorithm
int Status() const
return the Error Status of the last Integral calculation
void SetAbsTolerance(double absTolerance)
set the desired absolute Error
void SetFunction(const IMultiGenFunction &f)
method to set the a generic integration function
void SetType(MCIntegration::Type type)
set integration method
wrapper to a multi-dim function withtout derivatives for Monte Carlo multi-dimensional integration al...
void SetFunction(const FuncType &f)
Fill gsl function structure from a C++ Function class.
void SetFuncPointer(GSLMonteFuncPointer f)
GSLRandomEngine Base class for all GSL random engines, normally user instantiate the derived classes ...
GSLRngWrapper class to wrap gsl_rng structure.
Documentation for the abstract class IBaseFunctionMultiDim.
Definition IFunction.h:62
Generic interface for defining configuration options of a numerical algorithm.
Definition IOptions.h:31
Numerical multi dimensional integration options.
static ROOT::Math::IOptions * FindDefault(const char *name)
find specific options - return 0 if not existing
void SetIntegrator(const char *name)
set multi-dim integrator name
void SetNCalls(unsigned int calls)
set maximum number of function calls
unsigned int NCalls() const
maximum number of function calls
std::string Integrator() const
name of multi-dim integrator
Interface (abstract) class for multi numerical integration It must be implemented by the concrete Int...
Type
enumeration specifying the integration types.
double(* GSLMonteFuncPointer)(double *, size_t, void *)
Class for adapting any multi-dimension C++ functor class to C function pointers used by GSL MonteCarl...
Namespace for new Math classes and functions.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
structures collecting parameters for MISER multidimensional integration
structures collecting parameters for VEGAS multidimensional integration FOr implementation of default...
void ws()
Definition ws.C:66