1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Base/Math/IntegratorMCMiser.h
//! @brief Defines and implements template class IntegratorMCMiser.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
#ifdef SWIG
#error no need to expose this header to Swig
#endif // SWIG
#ifndef BORNAGAIN_BASE_MATH_INTEGRATORMCMISER_H
#define BORNAGAIN_BASE_MATH_INTEGRATORMCMISER_H
#include <gsl/gsl_monte_miser.h>
#include <memory>
//! Alias template for member function with signature double f(double)
template <typename T>
using miser_integrand = double (T::*)(const double*, size_t, const void*) const;
//! Template class to use Monte Carlo MISER integration of class member functions.
//!
//! Wraps an integrator from GNU Scientific Library.
//! Since this class holds a persistent workspace, we need at least one instance per thread.
//! Standard usage for integration inside a class T:
//! - Create a handle to an integrator:
//! 'auto integrator = make_integrator_miser(this, mem_function, dimension)'
//! - Call: 'integrator.integrate(lmin, lmax, data, n_points)'
template <typename T> class IntegratorMCMiser {
public:
//! structure holding the object and possible extra parameters
struct CallBackHolder {
const T* m_object_pointer;
const miser_integrand<T> m_member_function;
const void* m_data;
};
//! to integrate p_member_function, which must belong to p_object
IntegratorMCMiser(const T* p_object, miser_integrand<T> p_member_function, size_t dim);
IntegratorMCMiser(IntegratorMCMiser&&) = default;
~IntegratorMCMiser();
//! perform the actual integration over the ranges [min_array, max_array]
double integrate(double* min_array, double* max_array, const void* params,
size_t n_points) const;
private:
//! static function that can be passed to gsl integrator
static double StaticCallBack(double* d_array, size_t dim, void* v)
{
auto* p_cb = static_cast<CallBackHolder*>(v);
auto mf = static_cast<miser_integrand<T>>(p_cb->m_member_function);
return (p_cb->m_object_pointer->*mf)(d_array, dim, p_cb->m_data);
}
const T* m_object;
const miser_integrand<T> m_member_function;
const size_t m_dim; //!< dimension of the integration
gsl_monte_miser_state* m_gsl_workspace;
gsl_rng* m_random_gen;
};
//! Alias template for handle to a miser integrator
template <typename T> using P_integrator_miser = std::unique_ptr<IntegratorMCMiser<T>>;
//! Template function to create an integrator object
template <typename T>
P_integrator_miser<T> make_integrator_miser(const T* object, miser_integrand<T> mem_function,
size_t dim)
{
P_integrator_miser<T> P_integrator(new IntegratorMCMiser<T>(object, mem_function, dim));
return P_integrator;
}
// ************************************************************************************************
// Implementation
// ************************************************************************************************
template <typename T>
IntegratorMCMiser<T>::IntegratorMCMiser(const T* p_object, miser_integrand<T> p_member_function,
size_t dim)
: m_object(p_object)
, m_member_function(p_member_function)
, m_dim(dim)
, m_gsl_workspace{nullptr}
{
m_gsl_workspace = gsl_monte_miser_alloc(m_dim);
const gsl_rng_type* random_type;
gsl_rng_env_setup();
random_type = gsl_rng_default;
m_random_gen = gsl_rng_alloc(random_type);
}
template <typename T> IntegratorMCMiser<T>::~IntegratorMCMiser()
{
gsl_monte_miser_free(m_gsl_workspace);
gsl_rng_free(m_random_gen);
}
template <typename T>
double IntegratorMCMiser<T>::integrate(double* min_array, double* max_array, const void* params,
size_t n_points) const
{
CallBackHolder cb = {m_object, m_member_function, params};
gsl_monte_function f;
f.f = StaticCallBack;
f.dim = m_dim;
f.params = &cb;
double result, error;
gsl_monte_miser_integrate(&f, min_array, max_array, m_dim, n_points, m_random_gen,
m_gsl_workspace, &result, &error);
return result;
}
#endif // BORNAGAIN_BASE_MATH_INTEGRATORMCMISER_H
|
