Skip to content

Commit

Permalink
Merge pull request #1080 from boostorg/issue184_2024
Browse files Browse the repository at this point in the history
Reinstate root finding protection against huge jumps.
  • Loading branch information
jzmaddock authored Feb 6, 2024
2 parents 75dcb3e + 04c2c24 commit 03ea9c8
Show file tree
Hide file tree
Showing 18 changed files with 160 additions and 121 deletions.
3 changes: 2 additions & 1 deletion include/boost/math/distributions/binomial.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -258,7 +258,7 @@ namespace boost

typedef typename Policy::discrete_quantile_type discrete_quantile_type;
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
return detail::inverse_discrete_quantile(
result = detail::inverse_discrete_quantile(
dist,
comp ? q : p,
comp,
Expand All @@ -267,6 +267,7 @@ namespace boost
RealType(1),
discrete_quantile_type(),
max_iter);
return result;
} // quantile

}
Expand Down
5 changes: 2 additions & 3 deletions include/boost/math/distributions/chi_squared.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -327,9 +327,8 @@ RealType chi_squared_distribution<RealType, Policy>::find_degrees_of_freedom(
RealType result = r.first + (r.second - r.first) / 2;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
" either there is no answer to how many degrees of freedom are required"
" or the answer is infinite. Current best guess is %1%", result, Policy());
policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
" either there is no answer to how many degrees of freedom are required or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
}
return result;
}
Expand Down
19 changes: 6 additions & 13 deletions include/boost/math/distributions/detail/generic_mode.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,7 @@ typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::val
// Oops we don't know how to handle this, or even in which
// direction we should move in, treat as an evaluation error:
//
return policies::raise_evaluation_error(
function,
"Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type());
return policies::raise_evaluation_error(function, "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type()); // LCOV_EXCL_LINE
}
do
{
Expand Down Expand Up @@ -82,11 +80,9 @@ typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::val
max_iter).first;
if(max_iter >= policies::get_max_root_iterations<policy_type>())
{
return policies::raise_evaluation_error<value_type>(
function,
"Unable to locate solution in a reasonable time:"
" either there is no answer to the mode of the distribution"
" or the answer is infinite. Current best guess is %1%", result, policy_type());
return policies::raise_evaluation_error<value_type>(function, // LCOV_EXCL_LINE
"Unable to locate solution in a reasonable time: either there is no answer to the mode of the distribution" // LCOV_EXCL_LINE
" or the answer is infinite. Current best guess is %1%", result, policy_type()); // LCOV_EXCL_LINE
}
return result;
}
Expand Down Expand Up @@ -135,11 +131,8 @@ typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::
max_iter).first;
if(max_iter >= policies::get_max_root_iterations<policy_type>())
{
return policies::raise_evaluation_error<value_type>(
function,
"Unable to locate solution in a reasonable time:"
" either there is no answer to the mode of the distribution"
" or the answer is infinite. Current best guess is %1%", result, policy_type());
return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
" either there is no answer to the mode of the distribution or the answer is infinite. Current best guess is %1%", result, policy_type()); // LCOV_EXCL_LINE
}
return result;
}
Expand Down
5 changes: 2 additions & 3 deletions include/boost/math/distributions/detail/generic_quantile.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -84,9 +84,8 @@ typename Dist::value_type generic_quantile(const Dist& dist, const typename Dist
value_type result = ir.first + (ir.second - ir.first) / 2;
if(max_iter >= policies::get_max_root_iterations<forwarding_policy>())
{
return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:"
" either there is no answer to quantile"
" or the answer is infinite. Current best guess is %1%", result, forwarding_policy());
return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
" either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, forwarding_policy()); // LCOV_EXCL_LINE
}
return result;
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,7 @@ typename Dist::value_type
while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b))
{
if(count == 0)
return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type());
return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type()); // LCOV_EXCL_LINE
a = b;
fa = fb;
b *= multiplier;
Expand All @@ -242,7 +242,7 @@ typename Dist::value_type
return 0;
}
if(count == 0)
return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type());
return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type()); // LCOV_EXCL_LINE
b = a;
fb = fa;
a /= multiplier;
Expand Down Expand Up @@ -276,6 +276,11 @@ typename Dist::value_type
//
std::pair<value_type, value_type> r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type());
max_iter += count;
if (max_iter >= policies::get_max_root_iterations<policy_type>())
{
return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
" either there is no answer to quantile or the answer is infinite. Current best guess is %1%", r.first, policy_type()); // LCOV_EXCL_LINE
}
BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
return (r.first + r.second) / 2;
}
Expand Down
13 changes: 4 additions & 9 deletions include/boost/math/distributions/find_scale.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -80,9 +80,7 @@ namespace boost

if (result <= 0)
{ // If policy isn't to throw, return the scale <= 0.
policies::raise_evaluation_error<typename Dist::value_type>(function,
"Computed scale (%1%) is <= 0!" " Was the complement intended?",
result, Policy());
policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE
}
return result;
} // template <class Dist, class Policy> find_scale
Expand Down Expand Up @@ -140,9 +138,7 @@ namespace boost
// ( z - location) / (quantile(complement(Dist(), q))
if (result <= 0)
{ // If policy isn't to throw, return the scale <= 0.
policies::raise_evaluation_error<typename Dist::value_type>(function,
"Computed scale (%1%) is <= 0!" " Was the complement intended?",
result, Policy());
policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE
}
return result;
} // template <class Dist, class Policy, class Real1, class Real2, class Real3> typename Dist::value_type find_scale
Expand Down Expand Up @@ -192,9 +188,8 @@ namespace boost
// ( z - location) / (quantile(complement(Dist(), q))
if (result <= 0)
{ // If policy isn't to throw, return the scale <= 0.
policies::raise_evaluation_error<typename Dist::value_type>(function,
"Computed scale (%1%) is <= 0!" " Was the complement intended?",
result, policies::policy<>()); // This is only the default policy - also Want a version with Policy here.
policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", // LCOV_EXCL_LINE
result, policies::policy<>()); // This is only the default policy - also Want a version with Policy here. LCOV_EXCL_LINE
}
return result;
} // template <class Dist, class Real1, class Real2, class Real3> typename Dist::value_type find_scale
Expand Down
23 changes: 16 additions & 7 deletions include/boost/math/distributions/inverse_gaussian.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -388,11 +388,16 @@ inline RealType quantile(const inverse_gaussian_distribution<RealType, Policy>&
// digits used to control how accurate to try to make the result.
// To allow user to control accuracy versus speed,
int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.
using boost::math::tools::newton_raphson_iterate;
result =
newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType, Policy>(dist, p), guess, min, max, get_digits, m);
return result;
newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType, Policy>(dist, p), guess, min, max, get_digits, max_iter);
if (max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
" either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
}
return result;
} // quantile

template <class RealType, class Policy>
Expand Down Expand Up @@ -459,11 +464,15 @@ inline RealType quantile(const complemented2_type<inverse_gaussian_distribution<
// int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
// digits used to control how accurate to try to make the result.
int get_digits = policies::digits<RealType, Policy>();
std::uintmax_t m = policies::get_max_root_iterations<Policy>();
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
using boost::math::tools::newton_raphson_iterate;
result =
newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType, Policy>(c.dist, q), guess, min, max, get_digits, m);
return result;
result = newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType, Policy>(c.dist, q), guess, min, max, get_digits, max_iter);
if (max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
" either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
}
return result;
} // quantile

template <class RealType, class Policy>
Expand Down
24 changes: 18 additions & 6 deletions include/boost/math/distributions/kolmogorov_smirnov.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -379,10 +379,16 @@ inline RealType quantile(const kolmogorov_smirnov_distribution<RealType, Policy>

RealType k = detail::kolmogorov_smirnov_quantile_guess(p) / sqrt(n);
const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.

return tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor<RealType, Policy>(dist, p),
k, RealType(0), RealType(1), get_digits, m);
RealType result = tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor<RealType, Policy>(dist, p),
k, RealType(0), RealType(1), get_digits, max_iter);
if (max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
" either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
}
return result;
} // quantile

template <class RealType, class Policy>
Expand All @@ -403,11 +409,17 @@ inline RealType quantile(const complemented2_type<kolmogorov_smirnov_distributio
RealType k = detail::kolmogorov_smirnov_quantile_guess(RealType(1-p)) / sqrt(n);

const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.

return tools::newton_raphson_iterate(
RealType result = tools::newton_raphson_iterate(
detail::kolmogorov_smirnov_complementary_quantile_functor<RealType, Policy>(dist, p),
k, RealType(0), RealType(1), get_digits, m);
k, RealType(0), RealType(1), get_digits, max_iter);
if (max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
" either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
}
return result;
} // quantile (complemented)

template <class RealType, class Policy>
Expand Down
34 changes: 12 additions & 22 deletions include/boost/math/distributions/non_central_beta.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -109,9 +109,7 @@ namespace boost
}
if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
{
return policies::raise_evaluation_error(
"cdf(non_central_beta_distribution<%1%>, %1%)",
"Series did not converge, closest value was %1%", sum, pol);
return policies::raise_evaluation_error("cdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
}
}
return sum;
Expand Down Expand Up @@ -190,9 +188,7 @@ namespace boost
}
if(static_cast<std::uintmax_t>(i - k) > max_iter)
{
return policies::raise_evaluation_error(
"cdf(non_central_beta_distribution<%1%>, %1%)",
"Series did not converge, closest value was %1%", sum, pol);
return policies::raise_evaluation_error("cdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
}
last_term = term;
}
Expand All @@ -206,9 +202,7 @@ namespace boost
}
if(static_cast<std::uintmax_t>(count + k - i) > max_iter)
{
return policies::raise_evaluation_error(
"cdf(non_central_beta_distribution<%1%>, %1%)",
"Series did not converge, closest value was %1%", sum, pol);
return policies::raise_evaluation_error("cdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
}
pois *= i / l2;
beta -= xterm;
Expand Down Expand Up @@ -324,7 +318,7 @@ namespace boost
{
if(count == 0)
{
b = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol);
b = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol); // LCOV_EXCL_LINE
return std::make_pair(a, b);
}
//
Expand Down Expand Up @@ -362,7 +356,7 @@ namespace boost
}
if(count == 0)
{
a = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol);
a = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol); // LCOV_EXCL_LINE
return std::make_pair(a, b);
}
//
Expand Down Expand Up @@ -508,12 +502,14 @@ namespace boost

if(max_iter >= policies::get_max_root_iterations<Policy>())
{
// LCOV_EXCL_START
return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
" either there is no answer to quantile of the non central beta distribution"
" or the answer is infinite. Current best guess is %1%",
policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
function), Policy());
// LCOV_EXCL_STOP
}
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
Expand Down Expand Up @@ -583,9 +579,7 @@ namespace boost
}
if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
{
return policies::raise_evaluation_error(
"pdf(non_central_beta_distribution<%1%>, %1%)",
"Series did not converge, closest value was %1%", sum, pol);
return policies::raise_evaluation_error("pdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
}
}
return sum;
Expand Down Expand Up @@ -817,10 +811,8 @@ namespace boost
typedef typename Policy::assert_undefined_type assert_type;
static_assert(assert_type::value == 0, "Assert type is undefined.");

return policies::raise_evaluation_error<RealType>(
function,
"This function is not yet implemented, the only sensible result is %1%.",
std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity?
return policies::raise_evaluation_error<RealType>(function, "This function is not yet implemented, the only sensible result is %1%.", // LCOV_EXCL_LINE
std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity? LCOV_EXCL_LINE
}

template <class RealType, class Policy>
Expand All @@ -830,10 +822,8 @@ namespace boost
typedef typename Policy::assert_undefined_type assert_type;
static_assert(assert_type::value == 0, "Assert type is undefined.");

return policies::raise_evaluation_error<RealType>(
function,
"This function is not yet implemented, the only sensible result is %1%.",
std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity?
return policies::raise_evaluation_error<RealType>(function, "This function is not yet implemented, the only sensible result is %1%.", // LCOV_EXCL_LINE
std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity? LCOV_EXCL_LINE
} // kurtosis_excess

template <class RealType, class Policy>
Expand Down
Loading

0 comments on commit 03ea9c8

Please sign in to comment.