18#ifndef STAN_MATH_PRIM_FUNCTOR_HCUBATURE_HPP
19#define STAN_MATH_PRIM_FUNCTOR_HCUBATURE_HPP
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};
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};
55static constexpr std::array<double, 4>
gwd7{
56 1.2948496616886969327061143267787e-01, 2.797053914892766679014677714229e-01,
57 3.8183005050511894495036977548818e-01,
58 4.1795918367346938775510204081658e-01};
73inline void combination(Eigen::Matrix<int, Eigen::Dynamic, 1>& c,
const int dim,
74 const int p,
const int x) {
78 for (; k < x; r =
choose(dim - c[0], p - 1), k = k + r) {
82 for (
int i = 1; i < p - 1; i++) {
84 for (; k < x; r =
choose(dim - c[i], p - (i + 1)), k = k + r) {
89 c[p - 1] = (p > 1) ? c[p - 2] + x - k : x;
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++) {
110 for (
size_t j = 0; j < k; j++) {
111 p.coeffRef(c.coeff(j) - 1, i) = lambda;
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));
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++) {
140 index.push_back(
false);
141 p(c.array() - 1.0, current_col).setConstant(lambda);
143 for (
int j = 1; j != inner_iter_len; j++, current_col++) {
144 p.col(current_col) = p.col(current_col - 1);
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) {
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;
155 for (
int h = 0; h != first_zero + 1; h++) {
158 p(c.segment(0, first_zero + 1).array() - 1, current_col).array() *= -1;
177 std::array<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>, 4>,
178 Eigen::Matrix<double, 5, 1>, Eigen::Matrix<double, 4, 1>>
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);
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);
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));
217template <
typename F,
typename T_a,
typename T_b,
typename ParsPairT>
219 const ParsPairT& pars_pair) {
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);
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];
233 for (
auto i = 0; i != 7; i++) {
235 [](
auto&& integrand,
auto&& c,
auto&& delta,
auto&&... args) {
236 return integrand(c + delta, args...) + integrand(c - delta, args...);
238 pars_pair, integrand, c, deltax[i]);
241 Idash += fx *
gwd7[i / 2];
244 delta_t V =
fabs(delta);
247 return std::make_pair(I,
fabs(I - Idash));
268template <
typename F,
typename GenzMalik,
typename T_a,
typename T_b,
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);
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++) {
284 return std::make_tuple(Scalar(0.0), Scalar(0.0), 0);
286 c[i] = (a[i] + b[i]) / 2;
287 deltac[i] = (b[i] - a[i]) / 2.0;
290 for (std::size_t i = 0; i != dim; i++) {
293 Eigen::Matrix<Scalar, 5, 1> f = Eigen::Matrix<Scalar, 5, 1>::Zero();
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));
302 [](
auto&& integrand,
auto&& c,
auto&& p2,
auto&&... args) {
303 return integrand(c + p2, args...) + integrand(c - p2, args...);
305 pars_tuple, integrand, c, p2);
306 auto p3 = deltac.cwiseProduct(points[1].
col(i));
308 [](
auto&& integrand,
auto&& c,
auto&& p3,
auto&&... args) {
309 return integrand(c + p3, args...) + integrand(c - p3, args...);
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);
316 for (
auto i = 0; i != points[2].cols(); i++) {
318 [](
auto&& integrand,
auto&& cc,
auto&&... args) {
319 return integrand(cc, args...);
321 pars_tuple, integrand, c + deltac.cwiseProduct(points[2].
col(i)));
323 for (
auto i = 0; i != points[3].cols(); i++) {
325 [](
auto&& integrand,
auto&& cc,
auto&&... args) {
326 return integrand(cc, args...);
328 pars_tuple, integrand, c + deltac.cwiseProduct(points[3].
col(i)));
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);
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) {
342 maxdivdiff = divdiff[i];
343 }
else if ((
fabs(delta) <= deltaf) && (deltac[i] > deltac[kdivide])) {
347 return std::make_tuple(I, E, kdivide);
361template <
typename T_a,
typename T_b>
363 template <
typename Vec1,
typename Vec2>
365 :
a_(
std::forward<Vec1>(a)),
366 b_(
std::forward<Vec2>(b)),
369 Eigen::Matrix<T_a, Eigen::Dynamic, 1>
a_;
370 Eigen::Matrix<T_b, Eigen::Dynamic, 1>
b_;
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;
415 const int maxEval = max_eval <= 0 ? 1000000 : max_eval;
420 std::array<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>, 4>,
421 Eigen::Matrix<double, 5, 1>, Eigen::Matrix<double, 4, 1>>
425 std::tie(result, err)
429 std::tie(result, err, kdivide)
433 = (dim == 1) ? 15 : 1 + 4 * dim + 2 * dim * (dim - 1) + std::pow(2, dim);
435 auto evals_per_box = numevals;
439 if ((error <=
fmax(reqRelError *
fabs(val), reqAbsError))
440 || (numevals >= maxEval)) {
443 numevals += 2 * evals_per_box;
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);
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))
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());
463 eig_vec_b mb = Eigen::Map<const eig_vec_b>(box.b_.data(), box.b_.size());
472 std::tie(result_1, err_1)
474 std::tie(result_2, err_2)
478 integrand, genz_malik, dim, ma, box.b_, pars);
480 integrand, genz_malik, dim, box.a_, mb, pars);
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];
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;
497 for (
auto&& box : ms) {
500 for (
auto err_i : err_vec) {
auto col(T_x &&x, size_t j)
Return the specified column of the specified kernel generator expression using start-at-1 indexing.
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
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],...
static constexpr std::array< double, 8 > wd7
static constexpr std::array< double, 4 > gwd7
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.
fvar< T > fmax(const fvar< T > &x1, const fvar< T > &x2)
Return the greater of the two specified arguments.
constexpr decltype(auto) apply(F &&f, Tuple &&t, PreArgs &&... pre_args)
static constexpr double INFTY
Positive infinity.
fvar< T > fabs(const fvar< T > &x)
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
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...