Skip to content

Commit

Permalink
Merge pull request #714 from boostorg/tanh_sinh_performance
Browse files Browse the repository at this point in the history 
Improve tanh_sinh.
  • Loading branch information
@jzmaddock
jzmaddock committed Nov 3, 2021
2 parents afaa03a + 075ff0d commit 535fcc2
Showing 1 changed file with 130 additions and 21 deletions.
151 changes: 130 additions & 21 deletions include/boost/math/quadrature/detail/tanh_sinh_detail.hpp
View file
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

0 comments on commit 535fcc2

Please sign in to comment.