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

1d newton #1000

Open
wants to merge 20 commits into
base: develop
Choose a base branch
from
Open

1d newton #1000

Show file tree
Hide file tree
Changes from 13 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
d30f4c5
basic structure of solution
ryanelandt Jun 27, 2023
2cb0f12
created helper class for bracket
ryanelandt Jun 28, 2023
7bb9a80
added unit test new behavior
ryanelandt Jun 28, 2023
4da9504
passes unit tests locally, but needs refactor
ryanelandt Jul 6, 2023
4e6c6f6
update newton bounds
ryanelandt Jul 12, 2023
dd5693f
added special cases
ryanelandt Jul 12, 2023
acd0241
added more tests
ryanelandt Jul 12, 2023
89ac04c
Refactored into state and solvers
ryanelandt Jul 12, 2023
c690172
Merge pull request #1 from boostorg/develop
ryanelandt Jul 12, 2023
849e9d5
clarity improvements in comments and code
ryanelandt Jul 12, 2023
7a4c561
fix inspect issues
ryanelandt Jul 12, 2023
90582fc
fix CI errors: division by zero and array include
ryanelandt Jul 12, 2023
42b0bbc
direct init array<bool, 2> to fix clang 6 error
ryanelandt Jul 12, 2023
2cea888
fix another windows 1/0 compile error
ryanelandt Jul 22, 2023
e94fa65
Merge remote-tracking branch 'origin/develop'
ryanelandt Jul 31, 2023
6e29567
Merge branch 'develop' into 1d_newton
ryanelandt Jul 31, 2023
632a259
add BitSpace for smarter bisection
ryanelandt Aug 2, 2023
33319e7
Remove void_t because it requires C++17
ryanelandt Aug 2, 2023
2756324
refactor approach
ryanelandt Aug 4, 2023
07271f3
added () around min and max
ryanelandt Aug 4, 2023
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
4 changes: 2 additions & 2 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 @@ -382,7 +382,7 @@ inline RealType quantile(const kolmogorov_smirnov_distribution<RealType, Policy>
std::uintmax_t m = 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), boost::math::tools::max_value<RealType>(), get_digits, m);
k, RealType(0), RealType(1), get_digits, m);
ryanelandt marked this conversation as resolved.
Show resolved Hide resolved
} // quantile

template <class RealType, class Policy>
Expand All @@ -407,7 +407,7 @@ inline RealType quantile(const complemented2_type<kolmogorov_smirnov_distributio

return tools::newton_raphson_iterate(
detail::kolmogorov_smirnov_complementary_quantile_functor<RealType, Policy>(dist, p),
k, RealType(0), boost::math::tools::max_value<RealType>(), get_digits, m);
k, RealType(0), RealType(1), get_digits, m);
} // quantile (complemented)

template <class RealType, class Policy>
Expand Down
17 changes: 6 additions & 11 deletions include/boost/math/special_functions/beta.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1021,18 +1021,15 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de
BOOST_MATH_ASSERT((p_derivative == 0) || normalised);

if(!(boost::math::isfinite)(a))
return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
ryanelandt marked this conversation as resolved.
Show resolved Hide resolved
if(!(boost::math::isfinite)(b))
return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
if(!(boost::math::isfinite)(x))
return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
if (!(0 <= x && x <= 1))
return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);

if(p_derivative)
*p_derivative = -1; // value not set.

if((x < 0) || (x > 1))
return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);

if(normalised)
{
if(a < 0)
Expand Down Expand Up @@ -1422,18 +1419,16 @@ T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
// start with the usual error checks:
//
if (!(boost::math::isfinite)(a))
return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
if (!(boost::math::isfinite)(b))
return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
if (!(boost::math::isfinite)(x))
return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
if (!(0 <= x && x <= 1))
return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);

if(a <= 0)
return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
if(b <= 0)
return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol);
if((x < 0) || (x > 1))
return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
//
// Now the corner cases:
//
Expand Down
23 changes: 12 additions & 11 deletions include/boost/math/special_functions/detail/ibeta_inverse.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,23 +25,16 @@ namespace boost{ namespace math{ namespace detail{
template <class T>
struct temme_root_finder
{
temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {}
temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {
const T x_extrema = 1 / (1 + a);
BOOST_MATH_ASSERT(0 < x_extrema && x_extrema < 1);
}

boost::math::tuple<T, T> operator()(T x)
{
BOOST_MATH_STD_USING // ADL of std names

T y = 1 - x;
if(y == 0)
{
T big = tools::max_value<T>() / 4;
return boost::math::make_tuple(static_cast<T>(-big), static_cast<T>(-big));
}
if(x == 0)
{
T big = tools::max_value<T>() / 4;
return boost::math::make_tuple(static_cast<T>(-big), big);
}
T f = log(x) + a * log(y) + t;
T f1 = (1 / x) - (a / (y));
return boost::math::make_tuple(f, f1);
Expand Down Expand Up @@ -410,6 +403,14 @@ T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol)
T lower = eta < mu ? cross : 0;
T upper = eta < mu ? 1 : cross;
T x = (lower + upper) / 2;

// Special cases due to lack of numerical precision
if (cross == 1) {
return 1;
} else if (cross == 0) {
return 0;
}

x = tools::newton_raphson_iterate(
temme_root_finder<T>(u, mu), x, lower, upper, policies::digits<T, Policy>() / 2);
#ifdef BOOST_INSTRUMENT
Expand Down
Loading