Automatic Differentiation
 
Loading...
Searching...
No Matches
integrate_1d.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_HPP
2#define STAN_MATH_PRIM_FUNCTOR_INTEGRATE_1D_HPP
3
8#include <boost/math/quadrature/exp_sinh.hpp>
9#include <boost/math/quadrature/sinh_sinh.hpp>
10#include <boost/math/quadrature/tanh_sinh.hpp>
11#include <cmath>
12#include <functional>
13#include <ostream>
14#include <vector>
15
16namespace stan {
17namespace math {
51template <typename F>
52inline double integrate(const F& f, double a, double b,
53 double relative_tolerance) {
54 static constexpr const char* function = "integrate";
55 double error1 = 0.0;
56 double error2 = 0.0;
57 double L1 = 0.0;
58 double L2 = 0.0;
59 size_t levels;
60
61 auto one_integral_convergence_check = [](auto& error1, auto& rel_tol,
62 auto& L1) {
63 if (error1 > rel_tol * L1) {
64 [error1]() STAN_COLD_PATH {
66 function, "error estimate of integral", error1, "",
67 " exceeds the given relative tolerance times norm of integral");
68 }();
69 }
70 };
71
72 auto two_integral_convergence_check
73 = [](auto& error1, auto& error2, auto& rel_tol, auto& L1, auto& L2) {
74 if (error1 > rel_tol * L1) {
75 [error1]() STAN_COLD_PATH {
77 function, "error estimate of integral below zero", error1, "",
78 " exceeds the given relative tolerance times norm of "
79 "integral below zero");
80 }();
81 }
82 if (error2 > rel_tol * L2) {
83 [error2]() STAN_COLD_PATH {
85 function, "error estimate of integral above zero", error2, "",
86 " exceeds the given relative tolerance times norm of "
87 "integral above zero");
88 }();
89 }
90 };
91
92 // if a or b is infinite, set xc argument to NaN (see docs above for user
93 // function for xc info)
94 auto f_wrap = [&f](double x) { return f(x, NOT_A_NUMBER); };
95 if (std::isinf(a) && std::isinf(b)) {
96 boost::math::quadrature::sinh_sinh<double> integrator;
97 double Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
98 &levels);
99 one_integral_convergence_check(error1, relative_tolerance, L1);
100 return Q;
101 } else if (std::isinf(a)) {
102 boost::math::quadrature::exp_sinh<double> integrator;
108 if (b <= 0.0) {
109 double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1,
110 &L1, &levels);
111 one_integral_convergence_check(error1, relative_tolerance, L1);
112 return Q;
113 } else {
114 boost::math::quadrature::tanh_sinh<double> integrator_right;
115 double Q = integrator.integrate(f_wrap, a, 0.0, relative_tolerance,
116 &error1, &L1, &levels)
117 + integrator_right.integrate(
118 f_wrap, 0.0, b, relative_tolerance, &error2, &L2, &levels);
119 two_integral_convergence_check(error1, error2, relative_tolerance, L1,
120 L2);
121 return Q;
122 }
123 } else if (std::isinf(b)) {
124 boost::math::quadrature::exp_sinh<double> integrator;
125 if (a >= 0.0) {
126 double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1,
127 &L1, &levels);
128 one_integral_convergence_check(error1, relative_tolerance, L1);
129 return Q;
130 } else {
131 boost::math::quadrature::tanh_sinh<double> integrator_left;
132 double Q = integrator_left.integrate(f_wrap, a, 0, relative_tolerance,
133 &error1, &L1, &levels)
134 + integrator.integrate(f_wrap, relative_tolerance, &error2,
135 &L2, &levels);
136 two_integral_convergence_check(error1, error2, relative_tolerance, L1,
137 L2);
138 return Q;
139 }
140 } else {
141 auto f_wrap = [&f](double x, double xc) { return f(x, xc); };
142 boost::math::quadrature::tanh_sinh<double> integrator;
143 if (a < 0.0 && b > 0.0) {
144 double Q = integrator.integrate(f_wrap, a, 0.0, relative_tolerance,
145 &error1, &L1, &levels)
146 + integrator.integrate(f_wrap, 0.0, b, relative_tolerance,
147 &error2, &L2, &levels);
148 two_integral_convergence_check(error1, error2, relative_tolerance, L1,
149 L2);
150 return Q;
151 } else {
152 double Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1,
153 &L1, &levels);
154 one_integral_convergence_check(error1, relative_tolerance, L1);
155 return Q;
156 }
157 }
158}
159
173template <typename F, typename... Args,
174 require_all_st_arithmetic<Args...>* = nullptr>
175inline double integrate_1d_impl(const F& f, double a, double b,
176 double relative_tolerance, std::ostream* msgs,
177 const Args&... args) {
178 static constexpr const char* function = "integrate_1d";
179 check_less_or_equal(function, "lower limit", a, b);
180 if (unlikely(a == b)) {
181 if (std::isinf(a)) {
182 throw_domain_error(function, "Integration endpoints are both", a, "", "");
183 }
184 return 0.0;
185 } else {
186 return integrate(
187 [&](const auto& x, const auto& xc) { return f(x, xc, msgs, args...); },
188 a, b, relative_tolerance);
189 }
190}
191
238template <typename F>
239inline double integrate_1d(const F& f, double a, double b,
240 const std::vector<double>& theta,
241 const std::vector<double>& x_r,
242 const std::vector<int>& x_i, std::ostream* msgs,
243 const double relative_tolerance
244 = std::sqrt(EPSILON)) {
245 return integrate_1d_impl(integrate_1d_adapter<F>(f), a, b, relative_tolerance,
246 msgs, theta, x_r, x_i);
247}
248
249} // namespace math
250} // namespace stan
251
252#endif
#define STAN_COLD_PATH
#define unlikely(x)
require_all_t< std::is_arithmetic< scalar_type_t< std::decay_t< Types > > >... > require_all_st_arithmetic
Require all of the scalar types satisfy std::is_arithmetic.
void check_less_or_equal(const char *function, const char *name, const T_y &y, const T_high &high, Idxs... idxs)
Throw an exception if y is not less than high.
static constexpr double NOT_A_NUMBER
(Quiet) not-a-number value.
Definition constants.hpp:56
static constexpr double EPSILON
Smallest positive value.
Definition constants.hpp:41
void throw_domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw a domain error with a consistently formatted message.
return_type_t< T_a, T_b, Args... > integrate_1d_impl(const F &f, const T_a &a, const T_b &b, double relative_tolerance, std::ostream *msgs, const Args &... args)
Return the integral of f from a to b to the given relative tolerance.
return_type_t< T_a, T_b, T_theta > integrate_1d(const F &f, const T_a &a, const T_b &b, const std::vector< T_theta > &theta, const std::vector< double > &x_r, const std::vector< int > &x_i, std::ostream *msgs, const double relative_tolerance)
Compute the integral of the single variable function f from a to b to within a specified relative tol...
double integrate(const F &f, double a, double b, double relative_tolerance)
Integrate a single variable function f from a to b to within a specified relative tolerance.
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
bool isinf(const stan::math::var &a)
Return 1 if the specified argument is positive infinity or negative infinity and 0 otherwise.
Definition std_isinf.hpp:16
Adapt the non-variadic integrate_1d arguments to the variadic integrate_1d_impl interface.