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 <cctype> // 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(nullptr),
82 fFunction(nullptr)
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 != nullptr) SetParameters( VegasParameters(*opts) );
94 }
95 else if (fType == MCIntegration::kMISER) {
97 if (opts != nullptr) 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(nullptr),
111 fFunction(nullptr)
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 != nullptr) SetParameters( VegasParameters(*opts) );
123 }
124 else if (fType == MCIntegration::kMISER) {
126 if (opts != nullptr) SetParameters( MiserParameters(*opts) );
127 }
128
129}
130
131
132
134{
135 // delete workspace
136 if (fWorkspace) delete fWorkspace;
137 if (fRng != nullptr && !fExtGen) delete fRng;
138 if (fFunction != nullptr) delete fFunction;
139 fRng = nullptr;
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
160 if(fFunction == nullptr) fFunction = new GSLMonteFunctionWrapper();
162 fDim = f.NDim();
163
164 // now we can initialize the workspace
165 DoInitialize();
166}
167
168void GSLMCIntegrator::SetFunction( GSLMonteFuncPointer f, unsigned int dim, void * p )
169{
170 // method to set the a generic integration function
171 if(fFunction == nullptr) fFunction = new GSLMonteFunctionWrapper();
173 fFunction->SetParams ( p );
174 fDim = dim;
175
176 // now we can initialize the workspace
177 DoInitialize();
178}
179
180
181
182double GSLMCIntegrator::Integral(const double* a, const double* b)
183{
184 // evaluate the Integral of a over the defined interval (a[],b[])
185 assert(fRng != nullptr);
186 gsl_rng* fr = fRng->Rng();
187 assert(fr != nullptr);
188 if (!CheckFunction()) return 0;
189
190 // initialize by creating the right WS
191 // (if dimension and type are different than previous calculation)
192 DoInitialize(); // this is still needed if type is changed after setting the function
193
195 {
197 assert(ws != nullptr);
198 fStatus = gsl_monte_vegas_integrate( fFunction->GetFunc(), (double *) a, (double*) b , fDim, fCalls, fr, ws->GetWS(), &fResult, &fError);
199 }
200 else if (fType == MCIntegration::kMISER)
201 {
203 assert(ws != nullptr);
204 fStatus = gsl_monte_miser_integrate( fFunction->GetFunc(), (double *) a, (double *) b , fDim, fCalls, fr, ws->GetWS(), &fResult, &fError);
205 }
206 else if (fType == MCIntegration::kPLAIN)
207 {
209 assert(ws != nullptr);
210 fStatus = gsl_monte_plain_integrate( fFunction->GetFunc(), (double *) a, (double *) b , fDim, fCalls, fr, ws->GetWS(), &fResult, &fError);
211 }
212 /**/
213 else
214 {
215
216 fResult = 0;
217 fError = 0;
218 fStatus = -1;
219 std::cerr << "GSLIntegrator - Error: Unknown integration type" << std::endl;
220 throw std::exception();
221 }
222
223 return fResult;
224
225}
226
227
228double GSLMCIntegrator::Integral(const GSLMonteFuncPointer & f, unsigned int dim, double* a, double* b, void * p )
229{
230 // evaluate the Integral for a function f over the defined interval (a[],b[])
231 SetFunction(f,dim,p);
232 return Integral(a,b);
233}
234
235
236/* to be added later
237 double GSLMCIntegrator::Integral(GSLMonteFuncPointer f, void * p, double* a, double* b)
238 {
239
240 }
241
242*/
243//MCIntegration::Type GSLMCIntegrator::MCType() const {return fType;}
244
245/**
246 return the Result of the last Integral calculation
247*/
248double GSLMCIntegrator::Result() const { return fResult; }
249
250/**
251 return the estimate of the absolute Error of the last Integral calculation
252*/
253double GSLMCIntegrator::Error() const { return fError; }
254
255/**
256 return the Error Status of the last Integral calculation
257*/
258int GSLMCIntegrator::Status() const { return fStatus; }
259
260
261// setter for control Parameters (getters are not needed so far )
262
263/**
264 set the desired relative Error
265*/
266void GSLMCIntegrator::SetRelTolerance(double relTol){ this->fRelTol = relTol; }
267
268/**
269 set the desired absolute Error
270*/
271void GSLMCIntegrator::SetAbsTolerance(double absTol){ this->fAbsTol = absTol; }
272
274 // delete previous exist generator
275 if (fRng && !fExtGen) delete fRng;
276 fRng = r.Engine();
277 fExtGen = true;
278}
279
281{
282 // create workspace according to the type
283 fType=type;
284 if (fWorkspace != nullptr) {
285 if (type == fWorkspace->Type() ) return;
286 delete fWorkspace; // delete because is a different type
287 fWorkspace = nullptr;
288 }
289 //create Workspace according to type
292 }
293 else if (type == MCIntegration::kMISER) {
295 }
296 else {
297 // default: use VEGAS
299 MATH_WARN_MSG("GSLMCIntegration","Invalid integration type : use Vegas as default");
301 }
303 }
304}
305
307{
308 // set the integration type using a string
309 std::string typeName = (type!=nullptr) ? type : "VEGAS";
310 if (type == nullptr) MATH_INFO_MSG("GSLMCIntegration::SetTypeName","use default Vegas integrator method");
311 std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );
312
313 MCIntegration::Type integType = MCIntegration::kVEGAS; // default
314
315 if (typeName == "PLAIN") {
316 integType = MCIntegration::kPLAIN;
317 }
318 else if (typeName == "MISER") {
319 integType = MCIntegration::kMISER;
320 }
321 else if (typeName != "VEGAS") {
322 MATH_WARN_MSG("GSLMCIntegration::SetTypeName","Invalid integration type : use Vegas as default");
323 }
324
325 // create the fWorkspace object
326 // if it exists already with the same type it will not be re-created
327 SetType(integType);
328}
329
330
332{
333 // set integration mode for VEGAS method
335 {
337 assert(ws != nullptr);
338 assert(ws->GetWS() != nullptr);
339 if(mode == MCIntegration::kIMPORTANCE) ws->GetWS()->mode = GSL_VEGAS_MODE_IMPORTANCE;
340 else if(mode == MCIntegration::kSTRATIFIED) ws->GetWS()->mode = GSL_VEGAS_MODE_STRATIFIED;
341 else if(mode == MCIntegration::kIMPORTANCE_ONLY) ws->GetWS()->mode = GSL_VEGAS_MODE_IMPORTANCE_ONLY;
342 }
343
344 else std::cerr << "Mode not matching integration type";
345}
346
348{
349 // set integration options
350 SetTypeName(opt.Integrator().c_str() );
353 fCalls = opt.NCalls();
354
355 // specific options
356 ROOT::Math::IOptions * extraOpt = opt.ExtraOptions();
357 if (extraOpt) {
359 VegasParameters p(*extraOpt);
361 }
362 else if (fType == MCIntegration::kMISER ) {
363 MiserParameters p(fDim); // if possible pass dimension
364 p = (*extraOpt);
366 }
367 else {
368 MATH_WARN_MSG("GSLMCIntegrator::SetOptions","Invalid options set for the chosen integration type - ignore them");
369 }
370 }
371}
372
373
375{
376 // set method parameters
378 {
380 assert(ws != nullptr);
381 ws->SetParameters(p);
382 }
383 else
384 MATH_ERROR_MSG("GSLIntegrator::SetParameters"," Parameters not matching integration type");
385}
386
388{
389 // set method parameters
391 {
393 assert(ws != nullptr);
394 ws->SetParameters(p);
395 }
396 else
397 MATH_ERROR_MSG("GSLIntegrator::SetParameters"," Parameters not matching integration type");
398}
399
400
402{
403 // initialize by setting integration type
404
405 if (fWorkspace == nullptr) return;
406 if (fDim > 0 && fDim == fWorkspace->NDim() && fType == fWorkspace->Type() )
407 return; // can use previously existing ws
408
409 // otherwise clear workspace (if existing)
410 if (fWorkspace->NDim() > 0 ) fWorkspace->Clear();
411 // and create a new one
413}
414
415
416
417//----------- methods specific for VEGAS
418
420{
421 // returns the error sigma from the last iteration of the VEGAS algorithm
423 {
425 assert (ws != nullptr);
426 return ws->GetWS()->sigma;
427 }
428 else
429 {
430 std::cerr << "Parameter not matching integration type";
431 return 0;
432 }
433
434}
435
436
437/**
438*/
440{
441 // returns chi-squared per degree of freedom for the estimate of the integral
443 {
445 assert(ws != nullptr);
446 return ws->GetWS()->chisq;
447 }
448 else
449 {
450 std::cerr << "Parameter not matching integration type";
451 return 0;
452 }
453}
454
455
456
458{
459 // internal method to check validity of GSL function pointer
460
461 if (fFunction && fFunction->GetFunc() ) return true;
462 MATH_ERROR_MSG("GSLMCIntegrator::CheckFunction","Function has not been specified");
463 return false;
464}
465
466const char * GSLMCIntegrator::GetTypeName() const {
467 if (fType == MCIntegration::kVEGAS) return "VEGAS";
468 if (fType == MCIntegration::kMISER) return "MISER";
469 if (fType == MCIntegration::kPLAIN) return "PLAIN";
470 return "UNDEFINED";
471}
472
477 opt.SetNCalls(fCalls);
478 opt.SetWKSize(0);
480 return opt;
481}
482
483/// return a new option object which is managed by user
484std::unique_ptr<ROOT::Math::IOptions> GSLMCIntegrator::ExtraOptions() const {
485 if (!fWorkspace) return nullptr;
486 return fWorkspace->Options();
487}
488
490 // set the additional options
492}
493
494
495} // namespace Math
496} // namespace ROOT
497
498
499
#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
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define a(i)
Definition RSha256.hxx:99
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char mode
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
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 void SetOptions(const ROOT::Math::IOptions &)=0
set options
virtual MCIntegration::Type Type() const =0
virtual std::unique_ptr< ROOT::Math::IOptions > Options() const =0
retrieve option pointer corresponding to parameters create a new object to be managed by the user
virtual bool Init(size_t dim)=0
initialize the workspace creating the GSL pointer if it is not there
void SetExtraOptions(const ROOT::Math::IOptions &opt)
Set the extra options for Vegas and Miser.
double Result() const override
return the type of the integration used
GSLMonteFunctionWrapper * fFunction
GSLMCIntegrator(MCIntegration::Type type=MCIntegration::kVEGAS, double absTol=-1, double relTol=-1, unsigned int calls=0)
constructor of GSL MCIntegrator.
ROOT::Math::IntegratorMultiDimOptions Options() const override
get the option used for the integration
double Integral(const GSLMonteFuncPointer &f, unsigned int dim, double *a, double *b, void *p=nullptr)
evaluate the Integral of a function f over the defined hypercube (a,b)
const char * GetTypeName() const
return the name
GSLMCIntegrator & operator=(const GSLMCIntegrator &)
void SetFunction(const IMultiGenFunction &f) override
method to set the a generic integration function
~GSLMCIntegrator() override
destructor
void SetRelTolerance(double relTolerance) override
set the desired relative Error
void SetOptions(const ROOT::Math::IntegratorMultiDimOptions &opt) override
set the integration options
void SetParameters(const VegasParameters &p)
set default parameters for VEGAS method
void SetGenerator(GSLRandomEngine &r)
set random number generator
void SetAbsTolerance(double absTolerance) override
set the desired absolute Error
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
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
double Error() const override
return the estimate of the absolute Error of the last Integral calculation
int Status() const override
return the Error Status of the last Integral calculation
void SetType(MCIntegration::Type type)
set integration method
std::unique_ptr< ROOT::Math::IOptions > ExtraOptions() const
get the specific options (for Vegas or Miser) in term of string- name.
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.
void SetParameters(const struct VegasParameters &p)
Documentation for the abstract class IBaseFunctionMultiDim.
Definition IFunction.h:61
Generic interface for defining configuration options of a numerical algorithm.
Definition IOptions.h:28
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
std::string Integrator() const override
name of multi-dim integrator
unsigned int NCalls() const
maximum number of function calls
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.
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
Structure collecting parameters for MISER multidimensional integration.
Structures collecting parameters for VEGAS multidimensional integration For implementation of default...