Stan Math Library  2.20.0
reverse mode automatic differentiation
integrate_1d.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_ARR_FUNCTOR_integrate_1d_HPP
2 #define STAN_MATH_PRIM_ARR_FUNCTOR_integrate_1d_HPP
3 
7 #include <boost/math/quadrature/exp_sinh.hpp>
8 #include <boost/math/quadrature/sinh_sinh.hpp>
9 #include <boost/math/quadrature/tanh_sinh.hpp>
10 #include <functional>
11 #include <limits>
12 #include <ostream>
13 #include <vector>
14 
15 namespace stan {
16 namespace math {
50 template <typename F>
51 inline double integrate(const F& f, double a, double b,
52  double relative_tolerance) {
53  double error1 = 0.0;
54  double error2 = 0.0;
55  double L1 = 0.0;
56  double L2 = 0.0;
57  bool used_two_integrals = false;
58  size_t levels;
59  double Q = 0.0;
60  if (std::isinf(a) && std::isinf(b)) {
61  auto f_wrap = [&](double x) {
62  return f(x, std::numeric_limits<double>::quiet_NaN());
63  };
64  boost::math::quadrature::sinh_sinh<double> integrator;
65  Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1, &levels);
66  } else if (std::isinf(a)) {
67  boost::math::quadrature::exp_sinh<double> integrator;
73  if (b <= 0.0) {
74  auto f_wrap = [&](double x) {
75  return f(-(x + b), std::numeric_limits<double>::quiet_NaN());
76  };
77  Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
78  &levels);
79  } else {
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());
83  };
84  Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
85  &levels)
86  + integrator_right.integrate(f_wrap, -b, 0, relative_tolerance,
87  &error2, &L2, &levels);
88  used_two_integrals = true;
89  }
90  } else if (std::isinf(b)) {
91  boost::math::quadrature::exp_sinh<double> integrator;
92  if (a >= 0.0) {
93  auto f_wrap = [&](double x) {
94  return f(x + a, std::numeric_limits<double>::quiet_NaN());
95  };
96  Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
97  &levels);
98  } else {
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());
102  };
103  Q = integrator.integrate(f_wrap, relative_tolerance, &error1, &L1,
104  &levels)
105  + integrator_right.integrate(f_wrap, a, 0, relative_tolerance,
106  &error2, &L2, &levels);
107  used_two_integrals = true;
108  }
109  } else {
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,
114  &levels)
115  + integrator.integrate(f_wrap, 0.0, b, relative_tolerance, &error2,
116  &L2, &levels);
117  used_two_integrals = true;
118  } else {
119  Q = integrator.integrate(f_wrap, a, b, relative_tolerance, &error1, &L1,
120  &levels);
121  }
122  }
123 
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,
128  "",
129  " exceeds the given relative tolerance times norm of "
130  "integral below zero");
131  }
132  if (error2 > relative_tolerance * L2) {
133  domain_error(function, "error estimate of integral above zero", error2,
134  "",
135  " exceeds the given relative tolerance times norm of "
136  "integral above zero");
137  }
138  } else {
139  if (error1 > relative_tolerance * L1) {
140  domain_error(
141  function, "error estimate of integral", error1, "",
142  " exceeds the given relative tolerance times norm of integral");
143  }
144  }
145  return Q;
146 }
147 
194 template <typename F>
195 inline double integrate_1d(
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";
202  check_less_or_equal(function, "lower limit", a, b);
203 
204  if (a == b) {
205  if (std::isinf(a))
206  domain_error(function, "Integration endpoints are both", a, "", "");
207  return 0.0;
208  } else {
209  return integrate(
210  std::bind<double>(f, std::placeholders::_1, std::placeholders::_2,
211  theta, x_r, x_i, &msgs),
212  a, b, relative_tolerance);
213  }
214 }
215 
216 } // namespace math
217 } // namespace stan
218 
219 #endif
int 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
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)
Definition: sqrt.hpp:13
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...

     [ Stan Home Page ] © 2011–2018, Stan Development Team.