Automatic Differentiation
 
Loading...
Searching...
No Matches
hcubature.hpp
Go to the documentation of this file.
1// For the code in hcubature.hpp, the package JuliaMath/HCubature.jl
2// written in Julia by Steven G. Johnson served as a template.
3// It comes with the following MIT "Expat" license:
4//
5// Copyright (c) 2017: Steven G. Johnson.
6//
7// Permission is hereby granted, free of charge, to any person obtaining
8// a copy of this software and associated documentation files
9// (the "Software"), to deal in the Software without restriction,
10// including without limitation the rights to use, copy, modify, merge,
11// publish, distribute, sublicense, and/or sell copies of the Software,
12// and to permit persons to whom the Software is furnished to do so,
13// subject to the following conditions:
14//
15// The above copyright notice and this permission notice shall be
16// included in all copies or substantial portions of the Software.
17
18#ifndef STAN_MATH_PRIM_FUNCTOR_HCUBATURE_HPP
19#define STAN_MATH_PRIM_FUNCTOR_HCUBATURE_HPP
20
28#include <queue>
29#include <tuple>
30
31namespace stan {
32namespace math {
33
34namespace internal {
35
36static constexpr std::array<double, 8> xd7{
37 -9.9145537112081263920685469752598e-01,
38 -9.4910791234275852452618968404809e-01,
39 -8.6486442335976907278971278864098e-01,
40 -7.415311855993944398638647732811e-01,
41 -5.8608723546769113029414483825842e-01,
42 -4.0584515137739716690660641207707e-01,
43 -2.0778495500789846760068940377309e-01};
44
45static constexpr std::array<double, 8> wd7{
46 2.2935322010529224963732008059913e-02,
47 6.3092092629978553290700663189093e-02,
48 1.0479001032225018383987632254189e-01,
49 1.4065325971552591874518959051021e-01,
50 1.6900472663926790282658342659795e-01,
51 1.9035057806478540991325640242055e-01,
52 2.0443294007529889241416199923466e-01,
53 2.0948214108472782801299917489173e-01};
54
55static constexpr std::array<double, 4> gwd7{
56 1.2948496616886969327061143267787e-01, 2.797053914892766679014677714229e-01,
57 3.8183005050511894495036977548818e-01,
58 4.1795918367346938775510204081658e-01};
59
73inline void combination(Eigen::Matrix<int, Eigen::Dynamic, 1>& c, const int dim,
74 const int p, const int x) {
75 int r = 0;
76 int k = 0;
77 c[0] = 0;
78 for (; k < x; r = choose(dim - c[0], p - 1), k = k + r) {
79 c[0]++;
80 }
81 k = k - r;
82 for (int i = 1; i < p - 1; i++) {
83 c[i] = c[i - 1];
84 for (; k < x; r = choose(dim - c[i], p - (i + 1)), k = k + r) {
85 c[i]++;
86 }
87 k = k - r;
88 }
89 c[p - 1] = (p > 1) ? c[p - 2] + x - k : x;
90}
91
101template <typename Scalar>
102inline Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> combos(
103 const int k, const Scalar lambda, const int dim) {
104 Eigen::Matrix<int, Eigen::Dynamic, 1> c(k);
105 const auto choose_dimk = choose(dim, k);
106 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> p
107 = Eigen::MatrixXd::Zero(dim, choose_dimk);
108 for (size_t i = 0; i < choose_dimk; i++) {
109 combination(c, dim, k, i + 1);
110 for (size_t j = 0; j < k; j++) {
111 p.coeffRef(c.coeff(j) - 1, i) = lambda;
112 }
113 }
114 return p;
115}
116
127template <typename Scalar>
128inline Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> signcombos(
129 const int k, const Scalar lambda, const int dim) {
130 Eigen::Matrix<int, Eigen::Dynamic, 1> c(k);
131 const auto choose_dimk = choose(dim, k);
132 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> p
133 = Eigen::MatrixXd::Zero(dim, choose_dimk * std::pow(2, k));
134 int current_col = 0;
135 const auto inner_iter_len = std::pow(2, k);
136 std::vector<bool> index;
137 index.reserve(inner_iter_len * (choose_dimk + 1));
138 for (int i = 1; i != choose_dimk + 1; i++) {
139 combination(c, dim, k, i);
140 index.push_back(false);
141 p(c.array() - 1.0, current_col).setConstant(lambda);
142 current_col += 1;
143 for (int j = 1; j != inner_iter_len; j++, current_col++) {
144 p.col(current_col) = p.col(current_col - 1);
145 int first_zero
146 = std::distance(std::begin(index),
147 std::find(std::begin(index), std::end(index), false));
148 const std::size_t index_size = index.size();
149 if (first_zero == index_size) {
150 index.flip();
151 p(c.segment(0, index.size()).array() - 1, current_col).array() *= -1;
152 index.push_back(true);
153 p(c[index.size() - 1] - 1, current_col) = -lambda;
154 } else {
155 for (int h = 0; h != first_zero + 1; h++) {
156 index[h].flip();
157 }
158 p(c.segment(0, first_zero + 1).array() - 1, current_col).array() *= -1;
159 }
160 }
161 index.clear();
162 }
163 return p;
164}
165
176inline std::tuple<
177 std::array<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>, 4>,
178 Eigen::Matrix<double, 5, 1>, Eigen::Matrix<double, 4, 1>>
179make_GenzMalik(const int dim) {
180 std::array<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>, 4> points;
181 Eigen::Matrix<double, 5, 1> weights;
182 Eigen::Matrix<double, 4, 1> weights_low_deg;
183 double twopn = std::pow(2, dim);
184 weights[0]
185 = twopn * ((12824.0 - 9120.0 * dim + 400.0 * dim * dim) * 1.0 / 19683.0);
186 weights[1] = twopn * (980.0 / 6561.0);
187 weights[2] = twopn * ((1820.0 - 400.0 * dim) * 1.0 / 19683.0);
188 weights[3] = twopn * (200.0 / 19683.0);
189 weights[4] = 6859.0 / 19683.0;
190 weights_low_deg[0] = twopn * ((729 - 950 * dim + 50 * dim * dim) * 1.0 / 729);
191 weights_low_deg[1] = twopn * (245.0 / 486);
192 weights_low_deg[2] = twopn * ((265.0 - 100.0 * dim) * 1.0 / 1458.0);
193 weights_low_deg[3] = twopn * (25.0 / 729.0);
194 points[0] = combos(1, std::sqrt(9.0 * 1.0 / 70.0), dim);
195 double l3 = std::sqrt(9.0 * 1.0 / 10.0);
196 points[1] = combos(1, l3, dim);
197 points[2] = signcombos(2, l3, dim);
198 points[3] = signcombos(dim, std::sqrt(9.0 * 1.0 / 19.0), dim);
199 return std::make_tuple(std::move(points), std::move(weights),
200 std::move(weights_low_deg));
201}
202
217template <typename F, typename T_a, typename T_b, typename ParsPairT>
218inline auto gauss_kronrod(const F& integrand, const T_a a, const T_b b,
219 const ParsPairT& pars_pair) {
220 using delta_t = return_type_t<T_a, T_b>;
221 const delta_t c = 0.5 * (a + b);
222 const delta_t delta = 0.5 * (b - a);
223 auto f0 = math::apply([](auto&& integrand, auto&& c,
224 auto&&... args) { return integrand(c, args...); },
225 pars_pair, integrand, c);
226
227 auto I = f0 * wd7[7];
228 auto Idash = f0 * gwd7[3];
229 std::array<delta_t, 7> deltax;
230 for (int i = 0; i < 7; ++i) {
231 deltax[i] = delta * xd7[i];
232 }
233 for (auto i = 0; i != 7; i++) {
234 auto fx = math::apply(
235 [](auto&& integrand, auto&& c, auto&& delta, auto&&... args) {
236 return integrand(c + delta, args...) + integrand(c - delta, args...);
237 },
238 pars_pair, integrand, c, deltax[i]);
239 I += fx * wd7[i];
240 if (i % 2 == 1) {
241 Idash += fx * gwd7[i / 2];
242 }
243 }
244 delta_t V = fabs(delta);
245 I *= V;
246 Idash *= V;
247 return std::make_pair(I, fabs(I - Idash));
248}
249
268template <typename F, typename GenzMalik, typename T_a, typename T_b,
269 typename ParsTupleT>
270inline auto integrate_GenzMalik(const F& integrand, const GenzMalik& genz_malik,
271 const int dim,
272 const Eigen::Matrix<T_a, Eigen::Dynamic, 1>& a,
273 const Eigen::Matrix<T_b, Eigen::Dynamic, 1>& b,
274 const ParsTupleT& pars_tuple) {
275 auto&& points = std::get<0>(genz_malik);
276 auto&& weights = std::get<1>(genz_malik);
277 auto&& weights_low_deg = std::get<2>(genz_malik);
278 using delta_t = return_type_t<T_a, T_b>;
279 Eigen::Matrix<delta_t, Eigen::Dynamic, 1> c(dim);
280 Eigen::Matrix<delta_t, Eigen::Dynamic, 1> deltac(dim);
282 for (auto i = 0; i != dim; i++) {
283 if (a[i] == b[i]) {
284 return std::make_tuple(Scalar(0.0), Scalar(0.0), 0);
285 }
286 c[i] = (a[i] + b[i]) / 2;
287 deltac[i] = (b[i] - a[i]) / 2.0;
288 }
289 delta_t v = 1.0;
290 for (std::size_t i = 0; i != dim; i++) {
291 v *= deltac[i];
292 }
293 Eigen::Matrix<Scalar, 5, 1> f = Eigen::Matrix<Scalar, 5, 1>::Zero();
294 f.coeffRef(0)
295 = math::apply([](auto&& integrand, auto&& c,
296 auto&&... args) { return integrand(c, args...); },
297 pars_tuple, integrand, c);
298 Eigen::Matrix<Scalar, Eigen::Dynamic, 1> divdiff(dim);
299 for (auto i = 0; i != dim; i++) {
300 auto p2 = deltac.cwiseProduct(points[0].col(i));
301 auto f2i = math::apply(
302 [](auto&& integrand, auto&& c, auto&& p2, auto&&... args) {
303 return integrand(c + p2, args...) + integrand(c - p2, args...);
304 },
305 pars_tuple, integrand, c, p2);
306 auto p3 = deltac.cwiseProduct(points[1].col(i));
307 auto f3i = math::apply(
308 [](auto&& integrand, auto&& c, auto&& p3, auto&&... args) {
309 return integrand(c + p3, args...) + integrand(c - p3, args...);
310 },
311 pars_tuple, integrand, c, p3);
312 f.coeffRef(1) += f2i;
313 f.coeffRef(2) += f3i;
314 divdiff[i] = fabs(f3i + 12.0 * f.coeff(0) - 7.0 * f2i);
315 }
316 for (auto i = 0; i != points[2].cols(); i++) {
317 f.coeffRef(3) += math::apply(
318 [](auto&& integrand, auto&& cc, auto&&... args) {
319 return integrand(cc, args...);
320 },
321 pars_tuple, integrand, c + deltac.cwiseProduct(points[2].col(i)));
322 }
323 for (auto i = 0; i != points[3].cols(); i++) {
324 f.coeffRef(4) += math::apply(
325 [](auto&& integrand, auto&& cc, auto&&... args) {
326 return integrand(cc, args...);
327 },
328 pars_tuple, integrand, c + deltac.cwiseProduct(points[3].col(i)));
329 }
330
331 Scalar I = v * weights.dot(f);
332 Scalar Idash = v * weights_low_deg.dot(f.template head<4>());
333 Scalar E = fabs(I - Idash);
334
335 int kdivide = 0;
336 Scalar deltaf = E / (std::pow(10, dim) * v);
337 Scalar maxdivdiff = 0.0;
338 for (auto i = 0; i != dim; i++) {
339 Scalar delta = divdiff[i] - maxdivdiff;
340 if (delta > deltaf) {
341 kdivide = i;
342 maxdivdiff = divdiff[i];
343 } else if ((fabs(delta) <= deltaf) && (deltac[i] > deltac[kdivide])) {
344 kdivide = i;
345 }
346 }
347 return std::make_tuple(I, E, kdivide);
348}
349
361template <typename T_a, typename T_b>
362struct Box {
363 template <typename Vec1, typename Vec2>
364 Box(Vec1&& a, Vec2&& b, return_type_t<T_a, T_b> I, int kdivide)
365 : a_(std::forward<Vec1>(a)),
366 b_(std::forward<Vec2>(b)),
367 I_(I),
368 kdiv_(kdivide) {}
369 Eigen::Matrix<T_a, Eigen::Dynamic, 1> a_;
370 Eigen::Matrix<T_b, Eigen::Dynamic, 1> b_;
372 int kdiv_;
373};
374
375} // namespace internal
376
403template <typename F, typename T_a, typename T_b, typename ParsTuple,
404 typename TAbsErr, typename TRelErr>
405inline auto hcubature(const F& integrand, const ParsTuple& pars, const int dim,
406 const Eigen::Matrix<T_a, Eigen::Dynamic, 1>& a,
407 const Eigen::Matrix<T_b, Eigen::Dynamic, 1>& b,
408 const int max_eval, const TAbsErr reqAbsError,
409 const TRelErr reqRelError) {
411 using eig_vec_a = Eigen::Matrix<T_a, Eigen::Dynamic, 1>;
412 using eig_vec_b = Eigen::Matrix<T_b, Eigen::Dynamic, 1>;
413 using namespace boost::math::quadrature;
414
415 const int maxEval = max_eval <= 0 ? 1000000 : max_eval;
416 Scalar result;
417 Scalar err;
418 auto kdivide = 0;
419 std::tuple<
420 std::array<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>, 4>,
421 Eigen::Matrix<double, 5, 1>, Eigen::Matrix<double, 4, 1>>
422 genz_malik;
423
424 if (dim == 1) {
425 std::tie(result, err)
426 = internal::gauss_kronrod(integrand, a[0], b[0], pars);
427 } else {
428 genz_malik = internal::make_GenzMalik(dim);
429 std::tie(result, err, kdivide)
430 = internal::integrate_GenzMalik(integrand, genz_malik, dim, a, b, pars);
431 }
432 auto numevals
433 = (dim == 1) ? 15 : 1 + 4 * dim + 2 * dim * (dim - 1) + std::pow(2, dim);
434
435 auto evals_per_box = numevals;
436 Scalar error = err;
437 Scalar val = result;
438
439 if ((error <= fmax(reqRelError * fabs(val), reqAbsError))
440 || (numevals >= maxEval)) {
441 return result;
442 }
443 numevals += 2 * evals_per_box;
444 using box_t = internal::Box<T_a, T_b>;
445 std::vector<box_t> ms;
446 ms.reserve(numevals);
447 ms.emplace_back(a, b, result, kdivide);
448 auto get_largest_box_idx = [](auto&& box_vec) {
449 auto max_it = std::max_element(box_vec.begin(), box_vec.end());
450 return std::distance(box_vec.begin(), max_it);
451 };
452 std::vector<Scalar> err_vec;
453 err_vec.reserve(numevals);
454 err_vec.push_back(err);
455 while ((numevals < maxEval)
456 && (error > fmax(reqRelError * fabs(val), reqAbsError))
457 && fabs(val) < INFTY) {
458 auto err_idx = get_largest_box_idx(err_vec);
459 auto&& box = ms[err_idx];
460 auto w = (box.b_[box.kdiv_] - box.a_[box.kdiv_]) / 2;
461 eig_vec_a ma = Eigen::Map<const eig_vec_a>(box.a_.data(), box.a_.size());
462 ma[box.kdiv_] += w;
463 eig_vec_b mb = Eigen::Map<const eig_vec_b>(box.b_.data(), box.b_.size());
464 mb[box.kdiv_] -= w;
465 int kdivide_1{0};
466 int kdivide_2{0};
467 Scalar result_1;
468 Scalar result_2;
469 Scalar err_1;
470 Scalar err_2;
471 if (dim == 1) {
472 std::tie(result_1, err_1)
473 = internal::gauss_kronrod(integrand, ma[0], box.b_[0], pars);
474 std::tie(result_2, err_2)
475 = internal::gauss_kronrod(integrand, box.a_[0], mb[0], pars);
476 } else {
477 std::tie(result_1, err_1, kdivide_1) = internal::integrate_GenzMalik(
478 integrand, genz_malik, dim, ma, box.b_, pars);
479 std::tie(result_2, err_2, kdivide_2) = internal::integrate_GenzMalik(
480 integrand, genz_malik, dim, box.a_, mb, pars);
481 }
482 box_t box1(std::move(ma), box.b_, result_1, kdivide_1);
483 box_t box2(box.a_, std::move(mb), result_2, kdivide_2);
484 result += result_1 + result_2 - box.I_;
485 err += err_1 + err_2 - err_vec[err_idx];
486 ms[err_idx].I_ = 0;
487 err_vec[err_idx] = 0;
488 ms.push_back(std::move(box1));
489 ms.push_back(std::move(box2));
490 err_vec.push_back(err_1);
491 err_vec.push_back(err_2);
492 numevals += 2 * evals_per_box;
493 }
494 result = 0.0;
495 err = 0.0;
496
497 for (auto&& box : ms) {
498 result += box.I_;
499 }
500 for (auto err_i : err_vec) {
501 err += err_i;
502 }
503 return result;
504} // hcubature
505
506} // namespace math
507} // namespace stan
508#endif
auto col(T_x &&x, size_t j)
Return the specified column of the specified kernel generator expression using start-at-1 indexing.
Definition col.hpp:23
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > combos(const int k, const Scalar lambda, const int dim)
Compute a matrix [p] of all [dim]-component vectors with [k] components equal to [lambda] and other c...
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > signcombos(const int k, const Scalar lambda, const int dim)
Compute a matrix [p] of all [dim]-component vectors with [k] components equal to [±lambda] and other ...
auto integrate_GenzMalik(const F &integrand, const GenzMalik &genz_malik, const int dim, const Eigen::Matrix< T_a, Eigen::Dynamic, 1 > &a, const Eigen::Matrix< T_b, Eigen::Dynamic, 1 > &b, const ParsTupleT &pars_tuple)
Compute the integral of the function to be integrated (integrand) from a to b for more than one dimen...
auto gauss_kronrod(const F &integrand, const T_a a, const T_b b, const ParsPairT &pars_pair)
Compute the integral of the function to be integrated (integrand) from a to b for one dimension.
static constexpr std::array< double, 8 > xd7
Definition hcubature.hpp:36
std::tuple< std::array< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >, 4 >, Eigen::Matrix< double, 5, 1 >, Eigen::Matrix< double, 4, 1 > > make_GenzMalik(const int dim)
Compute the points and weights corresponding to a [dim]-dimensional Genz-Malik cubature rule.
void combination(Eigen::Matrix< int, Eigen::Dynamic, 1 > &c, const int dim, const int p, const int x)
Get the [x]-th lexicographically ordered set of [p] elements in [dim] output is in [c],...
Definition hcubature.hpp:73
static constexpr std::array< double, 8 > wd7
Definition hcubature.hpp:45
static constexpr std::array< double, 4 > gwd7
Definition hcubature.hpp:55
auto hcubature(const F &integrand, const ParsTuple &pars, const int dim, const Eigen::Matrix< T_a, Eigen::Dynamic, 1 > &a, const Eigen::Matrix< T_b, Eigen::Dynamic, 1 > &b, const int max_eval, const TAbsErr reqAbsError, const TRelErr reqRelError)
Compute the [dim]-dimensional integral of the function from to within specified relative and absol...
int choose(int n, int k)
Return the binomial coefficient for the specified integer arguments.
Definition choose.hpp:29
fvar< T > fmax(const fvar< T > &x1, const fvar< T > &x2)
Return the greater of the two specified arguments.
Definition fmax.hpp:23
constexpr decltype(auto) apply(F &&f, Tuple &&t, PreArgs &&... pre_args)
Definition apply.hpp:52
static constexpr double INFTY
Positive infinity.
Definition constants.hpp:46
fvar< T > fabs(const fvar< T > &x)
Definition fabs.hpp:15
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
STL namespace.
Box(Vec1 &&a, Vec2 &&b, return_type_t< T_a, T_b > I, int kdivide)
Eigen::Matrix< T_b, Eigen::Dynamic, 1 > b_
return_type_t< T_a, T_b > I_
Eigen::Matrix< T_a, Eigen::Dynamic, 1 > a_
Compute the integral of the function to be integrated (integrand) from a to b for more than one dimen...