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

Reinstate root finding protection against huge jumps. #1080

Merged
merged 2 commits into from
Feb 6, 2024
Merged

Reinstate root finding protection against huge jumps. #1080

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a57fb0a
Reinstate root finding protection against huge jumps.
jzmaddock Feb 5, 2024
04c2c24
BOOST_ASSERT->BOOST_MATH_ASSERT.
jzmaddock Feb 6, 2024
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
3 changes: 2 additions & 1 deletion include/boost/math/distributions/binomial.hpp
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
9 changes: 7 additions & 2 deletions include/boost/math/distributions/detail/inv_discrete_quantile.hpp
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
Loading