Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve tanh_sinh. #714

Merged
merged 3 commits into from
Nov 3, 2021
Merged

Improve tanh_sinh. #714

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1113497
tanh_sinh: improve thrash detection.
jzmaddock Oct 25, 2021
035188f
tanh_sinh: detect and allow non-finite values at end points.
jzmaddock Oct 26, 2021
075ff0d
Better estimate error due to truncation at endpoints.
jzmaddock Oct 28, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
151 changes: 130 additions & 21 deletions include/boost/math/quadrature/detail/tanh_sinh_detail.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,6 +195,14 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) 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<F>()(std::declval<Real>(), std::declval<Real>())) result_type;

Real h = m_t_max / m_inital_row_length;
result_type I0 = half_pi<Real>() * 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
Expand Down Expand Up @@ -222,22 +230,42 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) 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<F>()(std::declval<Real>(), std::declval<Real>())) result_type;

Real h = m_t_max / m_inital_row_length;
result_type I0 = half_pi<Real>()*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;
Expand All @@ -249,9 +277,8 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) 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;
}
Expand Down Expand Up @@ -282,6 +309,7 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
h *= half<Real>();
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);
Expand Down Expand Up @@ -311,12 +339,45 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) 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]);
//
// 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)
{
// 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;
Expand All @@ -333,8 +394,8 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) 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;

Expand All @@ -359,15 +420,38 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) 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<Real>())
// 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:
if (abs(endpoint_error / 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 + endpoint_error;
break;
}
// Fall through and keep going, assume we've discovered a new feature of f(x)....
}
//
// Termination condition:
Expand All @@ -380,6 +464,31 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) 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());
}
}
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 += endpoint_error;
break;
}

Expand Down