Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <gsl/gsl_integration.h>
7 : #include <limits>
8 : #include <memory>
9 : #include <vector>
10 :
11 : /// \cond
12 : namespace PUP {
13 : class er;
14 : } // namespace PUP
15 : /// \endcond
16 :
17 : /// \ingroup NumericalAlgorithmsGroup
18 : /// Numerical integration algorithms
19 : namespace integration {
20 :
21 : namespace detail {
22 :
23 : // The GSL functions require the integrand to have this particular function
24 : // signature. In particular, any extra parameters to the functions must be
25 : // passed as a void* and re-interpreted appropriately.
26 : template <typename IntegrandType>
27 : double integrand(const double x, void* const params) {
28 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
29 : auto* const function = reinterpret_cast<IntegrandType*>(params);
30 : return (*function)(x);
31 : }
32 :
33 : class GslQuadAdaptiveImpl {
34 : public:
35 : explicit GslQuadAdaptiveImpl(size_t max_intervals) noexcept;
36 :
37 : GslQuadAdaptiveImpl() = default;
38 : GslQuadAdaptiveImpl(const GslQuadAdaptiveImpl&) = delete;
39 : GslQuadAdaptiveImpl& operator=(const GslQuadAdaptiveImpl&) = delete;
40 : GslQuadAdaptiveImpl(GslQuadAdaptiveImpl&&) noexcept = default;
41 : GslQuadAdaptiveImpl& operator=(GslQuadAdaptiveImpl&& rhs) noexcept = default;
42 : ~GslQuadAdaptiveImpl() noexcept = default;
43 :
44 : void pup(PUP::er& p) noexcept; // NOLINT(google-runtime-references)
45 :
46 : protected:
47 : template <typename IntegrandType>
48 : gsl_function* gsl_integrand(IntegrandType&& integrand) const noexcept {
49 : gsl_integrand_.function = &detail::integrand<IntegrandType>;
50 : gsl_integrand_.params = &integrand;
51 : return &gsl_integrand_;
52 : }
53 :
54 : struct GslIntegrationWorkspaceDeleter {
55 : void operator()(gsl_integration_workspace* workspace) const noexcept;
56 : };
57 :
58 : size_t max_intervals_ = 0;
59 : std::unique_ptr<gsl_integration_workspace, GslIntegrationWorkspaceDeleter>
60 : workspace_{};
61 :
62 : private:
63 : void initialize() noexcept;
64 :
65 : mutable gsl_function gsl_integrand_{};
66 : };
67 :
68 : void check_status_code(int status_code);
69 :
70 : void disable_gsl_error_handling() noexcept;
71 :
72 : } // namespace detail
73 :
74 : /// Each type specifies which algorithm from GSL should be used. It should be
75 : /// chosen according to the problem.
76 1 : enum class GslIntegralType {
77 : /// gsl_integration_qag()
78 : StandardGaussKronrod,
79 : /// gsl_integration_qags()
80 : IntegrableSingularitiesPresent,
81 : /// gsl_integration_qagp()
82 : IntegrableSingularitiesKnown,
83 : /// gsl_integration_qagi()
84 : InfiniteInterval,
85 : /// gsl_integration_qagiu()
86 : UpperBoundaryInfinite,
87 : /// gsl_integration_qagil()
88 : LowerBoundaryInfinite
89 : };
90 :
91 : /*!
92 : * \brief A wrapper around the GSL adaptive integration procedures
93 : *
94 : * All templates take upper bounds to the absolute and relative error;
95 : * `tolerance_abs` and `tolerance_rel(default = 0.0)`, respectively. To compute
96 : * to a specified absolute error, set `tolerance_rel` to zero. To compute to a
97 : * specified relative error, set `tolerance_abs` to zero. For details on the
98 : * algorithm see the GSL documentation on `gsl_integration`.
99 : *
100 : * Here is an example how to use this class. For the function:
101 : *
102 : * \snippet Test_GslQuadAdaptive.cpp integrated_function
103 : *
104 : * the integration should look like:
105 : *
106 : * \snippet Test_GslQuadAdaptive.cpp integration_example
107 : */
108 : template <GslIntegralType TheIntegralType>
109 1 : class GslQuadAdaptive;
110 :
111 : /*!
112 : * \brief Use Gauss-Kronrod rule to integrate a 1D-function
113 : *
114 : * The algorithm for "StandardGaussKronrod" uses the QAG algorithm to employ an
115 : * adaptive Gauss-Kronrod n-points integration rule. Its function takes a
116 : * `lower_boundary` and `upper_boundary` to the integration region, an upper
117 : * bound for the absolute error `tolerance_abs` and a `key`. The latter is an
118 : * index to the array [15, 21, 31, 41, 51, 61], where each element denotes how
119 : * many function evaluations take place in each subinterval. A higher-order rule
120 : * serves better for smooth functions, whereas a lower-order rule saves time for
121 : * functions with local difficulties, such as discontinuities.
122 : */
123 : template <>
124 : class GslQuadAdaptive<GslIntegralType::StandardGaussKronrod>
125 1 : : public detail::GslQuadAdaptiveImpl {
126 : public:
127 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
128 : template <typename IntegrandType>
129 0 : double operator()(IntegrandType&& integrand, const double lower_boundary,
130 : const double upper_boundary, const double tolerance_abs,
131 : const int key, const double tolerance_rel = 0.) const {
132 : double result = std::numeric_limits<double>::signaling_NaN();
133 : double error = std::numeric_limits<double>::signaling_NaN();
134 : detail::disable_gsl_error_handling();
135 : const int status_code = gsl_integration_qag(
136 : gsl_integrand(std::forward<IntegrandType>(integrand)), lower_boundary,
137 : upper_boundary, tolerance_abs, tolerance_rel, this->max_intervals_, key,
138 : this->workspace_.get(), &result, &error);
139 : detail::check_status_code(status_code);
140 : return result;
141 : }
142 : };
143 :
144 : /*!
145 : * \brief Integrates a 1D-function with singularities
146 : *
147 : * The algorithm for "IntegrableSingularitiesPresent" concentrates new,
148 : * increasingly smaller subintervals around an unknown singularity and makes
149 : * successive approximations to the integral which should converge towards a
150 : * limit. The integration region is defined by `lower_boundary` and
151 : * `upper_boundary`.
152 : */
153 : template <>
154 : class GslQuadAdaptive<GslIntegralType::IntegrableSingularitiesPresent>
155 1 : : public detail::GslQuadAdaptiveImpl {
156 : public:
157 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
158 : template <typename IntegrandType>
159 0 : double operator()(IntegrandType&& integrand, const double lower_boundary,
160 : const double upper_boundary, const double tolerance_abs,
161 : const double tolerance_rel = 0.) const {
162 : double result = std::numeric_limits<double>::signaling_NaN();
163 : double error = std::numeric_limits<double>::signaling_NaN();
164 : detail::disable_gsl_error_handling();
165 : const int status_code = gsl_integration_qags(
166 : gsl_integrand(std::forward<IntegrandType>(integrand)), lower_boundary,
167 : upper_boundary, tolerance_abs, tolerance_rel, this->max_intervals_,
168 : this->workspace_.get(), &result, &error);
169 : detail::check_status_code(status_code);
170 : return result;
171 : }
172 : };
173 :
174 : /*!
175 : * \brief Integrates a 1D-function where singularities are known
176 : *
177 : * The algorithm for "IntegrableSingularitiesKnown" uses user-defined
178 : * subintervals given by a vector of doubles `points`, where each element
179 : * denotes an interval boundary.
180 : */
181 : template <>
182 : class GslQuadAdaptive<GslIntegralType::IntegrableSingularitiesKnown>
183 1 : : public detail::GslQuadAdaptiveImpl {
184 : public:
185 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
186 : template <typename IntegrandType>
187 0 : double operator()(IntegrandType&& integrand,
188 : const std::vector<double>& points,
189 : const double tolerance_abs,
190 : const double tolerance_rel = 0.) const {
191 : double result = std::numeric_limits<double>::signaling_NaN();
192 : double error = std::numeric_limits<double>::signaling_NaN();
193 : detail::disable_gsl_error_handling();
194 : // The const_cast is necessary because GSL has a weird interface where
195 : // the number of singularities does not change, but it doesn't take
196 : // the argument as a `const double*`. However, the first thing
197 : // `gsl_integration_qagp` does internally is forward to another function
198 : // that does take `points` by `const double*`. If `gsl_integration_qagp`
199 : // were to change the size of `points` this code would be severely broken
200 : // because `std::vector` allocates with `new`, while GSL would likely use
201 : // `malloc` (or some other C allocator). Mixing (de)allocators in such a way
202 : // is undefined behavior.
203 : const int status_code = gsl_integration_qagp(
204 : gsl_integrand(std::forward<IntegrandType>(integrand)),
205 : const_cast<double*>(points.data()), // NOLINT
206 : points.size(), tolerance_abs, tolerance_rel, this->max_intervals_,
207 : this->workspace_.get(), &result, &error);
208 : detail::check_status_code(status_code);
209 : return result;
210 : }
211 : };
212 :
213 : /*!
214 : * \brief Integrates a 1D-function over the interval \f$ (-\infty, +\infty) \f$
215 : *
216 : * The algorithm for "InfiniteInterval" uses the semi-open interval
217 : * \f$ (0, 1] \f$ to map to an infinite interval \f$ (-\infty, +\infty) \f$. Its
218 : * function takes no parameters other than limits `tolerance_abs` and optional
219 : * `tolerance_rel` for the absolute error and the relative error.
220 : */
221 : template <>
222 : class GslQuadAdaptive<GslIntegralType::InfiniteInterval>
223 1 : : public detail::GslQuadAdaptiveImpl {
224 : public:
225 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
226 : template <typename IntegrandType>
227 0 : double operator()(IntegrandType&& integrand, const double tolerance_abs,
228 : const double tolerance_rel = 0.) const {
229 : double result = std::numeric_limits<double>::signaling_NaN();
230 : double error = std::numeric_limits<double>::signaling_NaN();
231 : detail::disable_gsl_error_handling();
232 : const int status_code = gsl_integration_qagi(
233 : gsl_integrand(std::forward<IntegrandType>(integrand)), tolerance_abs,
234 : tolerance_rel, this->max_intervals_, this->workspace_.get(), &result,
235 : &error);
236 : detail::check_status_code(status_code);
237 : return result;
238 : }
239 : };
240 :
241 : /*!
242 : * \brief Integrates a 1D-function over the interval \f$ [a, +\infty) \f$
243 : *
244 : * The algorithm for "UpperBoundaryInfinite" maps the semi-open interval
245 : * \f$ (0, 1] \f$ to a semi-infinite interval \f$ [a, +\infty) \f$, where
246 : * \f$ a \f$ is given by `lower_boundary`.
247 : */
248 : template <>
249 : class GslQuadAdaptive<GslIntegralType::UpperBoundaryInfinite>
250 1 : : public detail::GslQuadAdaptiveImpl {
251 : public:
252 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
253 : template <typename IntegrandType>
254 0 : double operator()(IntegrandType&& integrand, const double lower_boundary,
255 : const double tolerance_abs,
256 : const double tolerance_rel = 0.) const {
257 : double result = std::numeric_limits<double>::signaling_NaN();
258 : double error = std::numeric_limits<double>::signaling_NaN();
259 : detail::disable_gsl_error_handling();
260 : const int status_code = gsl_integration_qagiu(
261 : gsl_integrand(std::forward<IntegrandType>(integrand)), lower_boundary,
262 : tolerance_abs, tolerance_rel, this->max_intervals_,
263 : this->workspace_.get(), &result, &error);
264 : detail::check_status_code(status_code);
265 : return result;
266 : }
267 : };
268 :
269 : /*!
270 : * \brief Integrates a 1D-function over the interval \f$ (-\infty, b] \f$
271 : *
272 : * The algorithm for "LowerBoundaryInfinite" maps the semi-open interval
273 : * \f$ (0, 1] \f$ to a semi-infinite interval \f$ (-\infty, b] \f$, where
274 : * \f$ b \f$ is given by `upper_boundary`.
275 : */
276 : template <>
277 : class GslQuadAdaptive<GslIntegralType::LowerBoundaryInfinite>
278 1 : : public detail::GslQuadAdaptiveImpl {
279 : public:
280 : using detail::GslQuadAdaptiveImpl::GslQuadAdaptiveImpl;
281 : template <typename IntegrandType>
282 0 : double operator()(IntegrandType&& integrand, const double upper_boundary,
283 : const double tolerance_abs,
284 : const double tolerance_rel = 0.) const {
285 : double result = std::numeric_limits<double>::signaling_NaN();
286 : double error = std::numeric_limits<double>::signaling_NaN();
287 : detail::disable_gsl_error_handling();
288 : const int status_code = gsl_integration_qagil(
289 : gsl_integrand(std::forward<IntegrandType>(integrand)), upper_boundary,
290 : tolerance_abs, tolerance_rel, this->max_intervals_,
291 : this->workspace_.get(), &result, &error);
292 : detail::check_status_code(status_code);
293 : return result;
294 : }
295 : };
296 : } // namespace integration
|
