1 #ifndef STAN_MATH_PRIM_ARR_FUNCTOR_integrate_1d_HPP 2 #define STAN_MATH_PRIM_ARR_FUNCTOR_integrate_1d_HPP 7 #include <boost/math/quadrature/exp_sinh.hpp> 8 #include <boost/math/quadrature/sinh_sinh.hpp> 9 #include <boost/math/quadrature/tanh_sinh.hpp> 51 inline double integrate(
const F& f,
double a,
double b,
52 double relative_tolerance) {
57 bool used_two_integrals =
false;
61 auto f_wrap = [&](
double x) {
62 return f(x, std::numeric_limits<double>::quiet_NaN());
64 boost::math::quadrature::sinh_sinh<double> integrator;
65 Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1, &levels);
67 boost::math::quadrature::exp_sinh<double> integrator;
74 auto f_wrap = [&](
double x) {
75 return f(-(x + b), std::numeric_limits<double>::quiet_NaN());
77 Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
80 boost::math::quadrature::tanh_sinh<double> integrator_right;
81 auto f_wrap = [&](
double x) {
82 return f(-x, std::numeric_limits<double>::quiet_NaN());
84 Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
86 + integrator_right.integrate(f_wrap, -b, 0, relative_tolerance,
87 &error2, &L2, &levels);
88 used_two_integrals =
true;
91 boost::math::quadrature::exp_sinh<double> integrator;
93 auto f_wrap = [&](
double x) {
94 return f(x + a, std::numeric_limits<double>::quiet_NaN());
96 Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
99 boost::math::quadrature::tanh_sinh<double> integrator_right;
100 auto f_wrap = [&](
double x) {
101 return f(x, std::numeric_limits<double>::quiet_NaN());
103 Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
105 + integrator_right.integrate(f_wrap, a, 0, relative_tolerance,
106 &error2, &L2, &levels);
107 used_two_integrals =
true;
110 auto f_wrap = [&](
double x,
double xc) {
return f(x, xc); };
111 boost::math::quadrature::tanh_sinh<double> integrator;
112 if (a < 0.0 && b > 0.0) {
113 Q = integrator.integrate(f_wrap, a, 0.0, relative_tolerance, &error1, &L1,
115 + integrator.integrate(f_wrap, 0.0, b, relative_tolerance, &error2,
117 used_two_integrals =
true;
119 Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1,
124 static const char*
function =
"integrate";
125 if (used_two_integrals) {
126 if (error1 > relative_tolerance * L1) {
127 domain_error(
function,
"error estimate of integral below zero", error1,
129 " exceeds the given relative tolerance times norm of " 130 "integral below zero");
132 if (error2 > relative_tolerance * L2) {
133 domain_error(
function,
"error estimate of integral above zero", error2,
135 " exceeds the given relative tolerance times norm of " 136 "integral above zero");
139 if (error1 > relative_tolerance * L1) {
141 function,
"error estimate of integral", error1,
"",
142 " exceeds the given relative tolerance times norm of integral");
194 template <
typename F>
196 const F& f,
const double a,
const double b,
197 const std::vector<double>& theta,
const std::vector<double>& x_r,
198 const std::vector<int>& x_i, std::ostream& msgs,
199 const double relative_tolerance
200 =
std::sqrt(std::numeric_limits<double>::epsilon())) {
201 static const char*
function =
"integrate_1d";
206 domain_error(
function,
"Integration endpoints are both", a,
"",
"");
210 std::bind<double>(f, std::placeholders::_1, std::placeholders::_2,
211 theta, x_r, x_i, &msgs),
212 a, b, relative_tolerance);
int isinf(const stan::math::var &a)
Return 1 if the specified argument is positive infinity or negative infinity and 0 otherwise...
void check_less_or_equal(const char *function, const char *name, const T_y &y, const T_high &high)
Check if y is less or equal to high.
fvar< T > sqrt(const fvar< T > &x)
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...
void 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.
double integrate_1d(const F &f, const double a, const double b, const std::vector< double > &theta, const std::vector< double > &x_r, const std::vector< int > &x_i, std::ostream &msgs, const double relative_tolerance=std::sqrt(std::numeric_limits< double >::epsilon()))
Compute the integral of the single variable function f from a to b to within a specified relative tol...