From 11134978df32ccc0d77ab3599cc1a3f9be5d2c73 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Mon, 25 Oct 2021 18:29:14 +0100 Subject: [PATCH 1/3] tanh_sinh: improve thrash detection. Significantly reduces the number of iterations in many cases. --- .../quadrature/detail/tanh_sinh_detail.hpp | 76 +++++++++++++++++-- 1 file changed, 68 insertions(+), 8 deletions(-) diff --git a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp index 7cef821aae..e38d735fda 100644 --- a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp +++ b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp @@ -359,15 +359,75 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin // // If the error is increasing, and we're past level 4, something bad is very likely happening: // - if ((err > last_err) && (k > 4) && (++thrash_count > 1)) + if ((err * 1.5 > last_err) && (k > 4)) { - // We could raise an evaluation_error, but since we likely have some sort of result, just return the last one - // (ie before the error started going up) - I1 = I0; - L1_I1 = L1_I0; - --k; - err = last_err; - break; + bool terminate = false; + if ((++thrash_count > 1) && (last_err < 1e-3)) + // Probably just thrashing, abort: + terminate = true; + else if(thrash_count > 2) + // OK, terrible error, but giving up anyway! + terminate = true; + else if (last_err < boost::math::tools::root_epsilon()) + // Trying to squeeze precision that probably isn't there, abort: + terminate = true; + else + { + // Take a look at the end points, if there's significant new area being + // discovered, then we're not able to get close enough to the endpoints + // to ever find the integral: + std::size_t j1 = (std::min)(max_left_index, abscissa_row.size() - 1); + std::size_t j2 = (std::min)(max_right_index, abscissa_row.size() - 1); + Real x = abscissa_row[j1]; + Real xc = x; + Real w = weight_row[j1]; + if (j1 >= first_complement_index) + { + // We have stored x - 1: + BOOST_MATH_ASSERT(x < 0); + x = 1 + xc; + } + else + { + BOOST_MATH_ASSERT(x >= 0); + xc = x - 1; + } + + result_type ym = f(-x, xc); + + x = abscissa_row[j2]; + xc = x; + w = weight_row[j2]; + if (j2 >= first_complement_index) + { + // We have stored x - 1: + BOOST_MATH_ASSERT(x < 0); + x = 1 + xc; + } + else + { + BOOST_MATH_ASSERT(x >= 0); + xc = x - 1; + } + result_type yp = f(x, -xc); + + result_type term = (yp + ym) * w; + + if (abs(term / sum) > err) + terminate = true; + } + + if (terminate) + { + // We could raise an evaluation_error, but since we likely have some sort of result, just return the last one + // (ie before the error started going up) + I1 = I0; + L1_I1 = L1_I0; + --k; + err = last_err; + break; + } + // Fall through and keep going, assume we've discovered a new feature of f(x).... } // // Termination condition: From 035188f9748aec6fb0779612d830b1449c18c95f Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Tue, 26 Oct 2021 16:29:15 +0100 Subject: [PATCH 2/3] tanh_sinh: detect and allow non-finite values at end points. --- .../quadrature/detail/tanh_sinh_detail.hpp | 106 +++++++++++++++--- 1 file changed, 91 insertions(+), 15 deletions(-) diff --git a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp index e38d735fda..7d2ea43ce6 100644 --- a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp +++ b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp @@ -195,6 +195,14 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin using std::sqrt; using boost::math::constants::half; using boost::math::constants::half_pi; + + // + // The type of the result: + typedef decltype(std::declval()(std::declval(), std::declval())) result_type; + + Real h = m_t_max / m_inital_row_length; + result_type I0 = half_pi() * f(0, 1); + Real L1_I0 = abs(I0); // // We maintain 4 integer values: // max_left_position is the logical index of the abscissa value closest to the @@ -222,22 +230,42 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin while (max_right_position && fabs(m_abscissas[0][max_right_position]) < right_min_complement) --max_right_position; // + // Check for non-finite values at the end points: + // + result_type yp, ym; + do + { + yp = f(-1 - m_abscissas[0][max_left_position], m_abscissas[0][max_left_position]); + if ((boost::math::isfinite)(yp)) + break; + --max_left_position; + } while (m_abscissas[0][max_left_position] < 0); + do + { + ym = f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]); + if ((boost::math::isfinite)(ym)) + break; + --max_right_position; + } while (m_abscissas[0][max_right_position] < 0); + + if ((max_left_position == 0) && (max_right_position == 0)) + { + return policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", ym, Policy()); + } + + I0 += yp * m_weights[0][max_left_position] + ym * m_weights[0][max_right_position]; + L1_I0 += abs(yp * m_weights[0][max_left_position]) + abs(ym * m_weights[0][max_right_position]); + // // Assumption: left_min_complement/right_min_complement are sufficiently small that we only // ever decrement through the stored values that are complements (the negative ones), and // never ever hit the true abscissa values (positive stored values). // BOOST_MATH_ASSERT(m_abscissas[0][max_left_position] < 0); BOOST_MATH_ASSERT(m_abscissas[0][max_right_position] < 0); - // - // The type of the result: - typedef decltype(std::declval()(std::declval(), std::declval())) result_type; - Real h = m_t_max / m_inital_row_length; - result_type I0 = half_pi()*f(0, 1); - Real L1_I0 = abs(I0); for(size_t i = 1; i < m_abscissas[0].size(); ++i) { - if ((i > max_right_position) && (i > max_left_position)) + if ((i >= max_right_position) && (i >= max_left_position)) break; Real x = m_abscissas[0][i]; Real xc = x; @@ -249,9 +277,8 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin } else xc = x - 1; - result_type yp, ym; - yp = i <= max_right_position ? f(x, -xc) : 0; - ym = i <= max_left_position ? f(-x, xc) : 0; + yp = i < max_right_position ? f(x, -xc) : 0; + ym = i < max_left_position ? f(-x, xc) : 0; I0 += (yp + ym)*w; L1_I0 += (abs(yp) + abs(ym))*w; } @@ -311,12 +338,34 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin ++max_right_position; ++max_right_index; } + // + // We also check that our endpoints don't hit singularities: + // + do + { + yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]); + if ((boost::math::isfinite)(yp)) + break; + max_left_position -= 2; + --max_left_index; + } while (abscissa_row[max_left_index] < 0); + do + { + ym = f(1 + abscissa_row[max_right_index], -abscissa_row[max_right_index]); + if ((boost::math::isfinite)(ym)) + break; + --max_right_index; + max_right_position -= 2; + } while (abscissa_row[max_right_index] < 0); + + sum += yp * weight_row[max_left_index] + ym * weight_row[max_right_index]; + absum += abs(yp * weight_row[max_left_index]) + abs(ym * weight_row[max_right_index]); for(size_t j = 0; j < weight_row.size(); ++j) { // If both left and right abscissa values are out of bounds at this step // we can just stop this loop right now: - if ((j > max_left_index) && (j > max_right_index)) + if ((j >= max_left_index) && (j >= max_right_index)) break; Real x = abscissa_row[j]; Real xc = x; @@ -333,8 +382,8 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin xc = x - 1; } - result_type yp = j > max_right_index ? 0 : f(x, -xc); - result_type ym = j > max_left_index ? 0 : f(-x, xc); + yp = j >= max_right_index ? 0 : f(x, -xc); + ym = j >= max_left_index ? 0 : f(-x, xc); result_type term = (yp + ym)*w; sum += term; @@ -393,7 +442,7 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin xc = x - 1; } - result_type ym = f(-x, xc); + ym = f(-x, xc); x = abscissa_row[j2]; xc = x; @@ -409,7 +458,7 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin BOOST_MATH_ASSERT(x >= 0); xc = x - 1; } - result_type yp = f(x, -xc); + yp = f(x, -xc); result_type term = (yp + ym) * w; @@ -440,6 +489,33 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin // if (err <= abs(tolerance*L1_I1)) { + // + // A quick sanity check: have we at some point narrowed our boundaries as a result + // of non-finite values? If so let's check that the area isn't on an increasing + // trajectory at our new end point, and increase our error estimate by the last + // good value as an estimate for what we may have discarded. + // + if ((max_left_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_left_index + 1]) > left_min_complement)) + { + yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]) * weight_row[max_left_index]; + ym = f(-1 - abscissa_row[max_left_index - 1], abscissa_row[max_left_index - 1]) * weight_row[max_left_index - 1]; + if (abs(yp) > abs(ym)) + { + return policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Inetgration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", I1, Policy()); + } + err += abs(yp); + } + if ((max_right_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_right_index + 1]) > right_min_complement)) + { + yp = f(-1 - abscissa_row[max_right_index], abscissa_row[max_right_index]) * weight_row[max_right_index]; + ym = f(-1 - abscissa_row[max_right_index - 1], abscissa_row[max_right_index - 1]) * weight_row[max_right_index - 1]; + if (abs(yp) > abs(ym)) + { + return policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Inetgration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", I1, Policy()); + } + err += abs(yp); + } + break; } From 075ff0dc6cc87624c6ed367c2916ddd8ac7ab292 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 28 Oct 2021 12:05:12 +0100 Subject: [PATCH 3/3] Better estimate error due to truncation at endpoints. --- .../quadrature/detail/tanh_sinh_detail.hpp | 57 +++++-------------- 1 file changed, 15 insertions(+), 42 deletions(-) diff --git a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp index 7d2ea43ce6..b36ccf54ad 100644 --- a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp +++ b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp @@ -309,6 +309,7 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin h *= half(); result_type sum = 0; Real absum = 0; + Real endpoint_error = 0; auto const& abscissa_row = this->get_abscissa_row(k); auto const& weight_row = this->get_weight_row(k); std::size_t first_complement_index = this->get_first_complement_index(k); @@ -360,6 +361,17 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin sum += yp * weight_row[max_left_index] + ym * weight_row[max_right_index]; absum += abs(yp * weight_row[max_left_index]) + abs(ym * weight_row[max_right_index]); + // + // We estimate the error due to truncation as the value contributed by the two most extreme points. + // In most cases this is tiny and can be ignored, if it is significant then either the area of the + // integral is so far our in the tails that our exponent range can't reach it (example x^-p at double + // precision and p ~ 1), or our function is truncated near epsilon, and we have had to narrow our endpoints. + // In this latter case we may over-estimate the error, but this is the best we can do. + // In any event, we do not add endpoint_error to the error estimate until we terminate the main loop, + // otherwise it can make things appear to be non-converged, when in reality, they are as converged as they + // will ever be. + // + endpoint_error = absum; for(size_t j = 0; j < weight_row.size(); ++j) { @@ -425,44 +437,7 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin // Take a look at the end points, if there's significant new area being // discovered, then we're not able to get close enough to the endpoints // to ever find the integral: - std::size_t j1 = (std::min)(max_left_index, abscissa_row.size() - 1); - std::size_t j2 = (std::min)(max_right_index, abscissa_row.size() - 1); - Real x = abscissa_row[j1]; - Real xc = x; - Real w = weight_row[j1]; - if (j1 >= first_complement_index) - { - // We have stored x - 1: - BOOST_MATH_ASSERT(x < 0); - x = 1 + xc; - } - else - { - BOOST_MATH_ASSERT(x >= 0); - xc = x - 1; - } - - ym = f(-x, xc); - - x = abscissa_row[j2]; - xc = x; - w = weight_row[j2]; - if (j2 >= first_complement_index) - { - // We have stored x - 1: - BOOST_MATH_ASSERT(x < 0); - x = 1 + xc; - } - else - { - BOOST_MATH_ASSERT(x >= 0); - xc = x - 1; - } - yp = f(x, -xc); - - result_type term = (yp + ym) * w; - - if (abs(term / sum) > err) + if (abs(endpoint_error / sum) > err) terminate = true; } @@ -473,7 +448,7 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin I1 = I0; L1_I1 = L1_I0; --k; - err = last_err; + err = last_err + endpoint_error; break; } // Fall through and keep going, assume we've discovered a new feature of f(x).... @@ -503,7 +478,6 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin { return policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Inetgration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", I1, Policy()); } - err += abs(yp); } if ((max_right_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_right_index + 1]) > right_min_complement)) { @@ -513,9 +487,8 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin { return policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Inetgration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", I1, Policy()); } - err += abs(yp); } - + err += endpoint_error; break; }