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.

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

Feature/issue 1010 - Replace boost::math:: functions with cmath (Part 1) #1067

Merged
merged 10 commits into from
Dec 10, 2018
Merged

Feature/issue 1010 - Replace boost::math:: functions with cmath (Part 1) #1067

merged 10 commits into from
Dec 10, 2018

Conversation

andrjohns
Copy link
Collaborator

@andrjohns andrjohns commented Nov 18, 2018

Summary

Addresses part 1 of issue #1010, replacing boost::math:: functions that don't need input checking with their c++11 equivalents.

The following functions have been replaced:

  • asinh
  • cbrt
  • erf
  • erfc
  • expm1
  • isinf
  • isnan
  • round
  • trunc

Tests

No major changes to tests, only replacing boost::math::foo() with stan::math::foo() where relevant

Side Effects

N/A

Checklist

@andrjohns andrjohns changed the title Feature/issue 1010 std functs Feature/issue 1010 - Replace boost::math:: functions with cmath Nov 19, 2018
Copy link
Contributor

@bob-carpenter bob-carpenter left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks so much for taking on this massive PR. There are several easy, but systematic changes that need to be made across the PR.

@@ -14,7 +15,7 @@ namespace math {
template <typename T>
inline fvar<T> gamma_p(const fvar<T> &x1, const fvar<T> &x2) {
using boost::math::digamma;
using boost::math::lgamma;
using stan::math::lgamma;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be removed---no need to import stan::math::lgamma when you're in the stan::math namespace.

Copy link
Collaborator Author

@andrjohns andrjohns Nov 25, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a good point. There are a few functions in the rest of the library with unnecessary using statements, looks like there's already an issue for it: #426. I'll clean-up the statements in this pull and have a look at the rest in a separate pull

@@ -11,7 +11,7 @@ namespace math {
template <typename T>
inline fvar<T> gamma_q(const fvar<T>& x1, const fvar<T>& x2) {
using boost::math::digamma;
using boost::math::tgamma;
using stan::math::tgamma;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not going to flag all these, but all the rest of these using stan::math::X should be removed.

@@ -60,7 +61,7 @@ inline fvar<T> gamma_p(double x1, const fvar<T> &x2) {
if (is_inf(x1))
return fvar<T>(u, std::numeric_limits<double>::quiet_NaN());

T der2 = exp(-x2.val_ + (x1 - 1.0) * log(x2.val_) - boost::math::lgamma(x1));
T der2 = exp(-x2.val_ + (x1 - 1.0) * log(x2.val_) - stan::math::lgamma(x1));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should not need to be qualified with stan::math::, because without a using std::lgamma, the standard version isn't a candidate.

template <typename T>
inline fvar<T> hypot(const fvar<T>& x1, double x2) {
using std::sqrt;
template <typename T, typename T1>
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[optional]
I would prefer these all to be numbered, so this would be T1 and T2 rather than T and T1.

template <typename T>
inline fvar<T> hypot(double x1, const fvar<T>& x2) {
using std::sqrt;
template <typename T, typename T1>
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[optional]
Ordering is backwards here. I think it should be

template <typename T1, typename T2>
inline fvar<typename boost::promote_args<T1, T2>::type>
hypot(const T1& x1, const fvar<T2>& x2) { ... }

This way, the type names line up wit hthe that way, the type names line up with the variable names and everything stays in order. This kind of thing may seem minor, but it's a big deal for program readability.

@@ -104,7 +104,8 @@ inline var fdim(const var& a, const var& b) {
* @return The positive difference between the first and second
* arguments.
*/
inline var fdim(double a, const var& b) {
template <typename Ta>
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[optional]
Just use <typename T> --- there's nothing to confuse it with.

@@ -45,7 +48,7 @@ inline double lgamma(double x) {
* @return natural logarithm of the gamma function applied to
* argument
*/
inline double lgamma(int x) { return boost::math::lgamma(x, boost_policy_t()); }
inline double lgamma(int x) { return lgamma(static_cast<double>(x)); }
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think the cast will be necessary with C++11. See (4) in: https://en.cppreference.com/w/cpp/numeric/math/lgamma

@@ -52,7 +51,7 @@ namespace math {
*/
template <typename T>
inline typename boost::math::tools::promote_args<T>::type lmgamma(int k, T x) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When can T be anything but double (or int promoted to double)?

@@ -30,11 +31,11 @@ inline double Phi(double x) {
if (x < -37.5)
return 0;
else if (x < -5.0)
return 0.5 * boost::math::erfc(-INV_SQRT_2 * x);
return 0.5 * stan::math::erfc(-INV_SQRT_2 * x);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll stop writing these, too, but all these qualifications are not necessary.

*
* @param x Value to test.
* @return <code>1</code> if the value is infinite.
*/
inline int is_inf(double x) { return boost::math::isinf(x); }
inline bool is_inf(double x) { return std::isinf(x); }
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[comment]
Thanks --- bool is the correct return type here.

@bob-carpenter
Copy link
Contributor

I get that we have to sort out the thread-safe thing, but why persist in Boost's convention of returning NaN when the input is zero when std::lgamma correctly returns infinity?

My PR for just lgamma, which I have closed without merging, just went with the std::lgamma values. I agree with Ben for two reasons---returning infinity makes the most sense mathemaitcally and we should be following the standard library conventions wherever possible.

@andrjohns
Copy link
Collaborator Author

The difficulty with the standard library and errors is that it never actually throws. There's a clearer writeup here, but essentially the function returns -nan and the errno variable is set to EDOM or ERANGE depending on whether a range or domain error occurs. The user then has to check this errno variable after calling the function to determine whether an error occurred.

As far as my google-fu can tell, there's no way to 'enable' throwing on domain/range errors int he standard library functions.

@bob-carpenter
Copy link
Contributor

As far as my google-fu can tell, there's no way to 'enable' throwing on domain/range errors int he standard library functions.

I also want to continue to follow the standard library's throw protocol. This is why we explicitly specified Boost error policies before---to make sure they did not throw.

If you want to throw, the way to do it is:

...
double y = std::lgamma(x);
if (...error code is on...)
  throw ...;
return y;

The tricky part there is thread safety for error codes.

@bob-carpenter
Copy link
Contributor

@andrjohns : Please ping me when the changes are ready to review or if you need help getting them implemented.

@andrjohns
Copy link
Collaborator Author

@bob-carpenter just confirming the implementation for throwing here before I update all functions.

Using acosh as an example:

inline double acosh(double x) {
  errno = 0;
  double y = std::acosh(x);
  if (errno != 0)
    throw std::domain_error(strerror(errno));
  return y;
}

The errno variable has to be reset to 0 before each function call, otherwise an error code from a previous call could remain. In terms of thread-safety, since c++11 each thread has its own errno variable that gets written/called (https://en.cppreference.com/w/cpp/error/errno).

The alternative to this is manually checking for valid inputs (consistent with the standard library) and throwing prior to calling the function:

inline double acosh(double x) {
  check_greater_or_equal("acosh", "x", x, 1.0);
  return std::acosh(x);
}

Is there a preference?

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Nov 30, 2018 via email

@andrjohns
Copy link
Collaborator Author

Ah I understand what you mean now. The main reason that I was confused is that the current tests expect that the function throws, and so I was trying to maintain the current behaviour.

Won't implementing the standard library math functions in this way lead to inconsistencies with the boost math functions (i.e. the boost functions will throw while the standard functions will silently fail)?

@bgoodri
Copy link
Contributor

bgoodri commented Nov 30, 2018 via email

@andrjohns
Copy link
Collaborator Author

Right, but validating the input beforehand (i.e. check_greater_or_equal("acosh", "x", x, 1.0);) results in throwing, and Bob's saying that it should not throw at all?

@bgoodri
Copy link
Contributor

bgoodri commented Nov 30, 2018 via email

@sakrejda
Copy link
Contributor

@andrjohns I'm glad you find this as confusing as I do (I've had a few very similar conversations). It's confusing for everyone. I don't think there's a clear answer for what to do when switching from boost::lgamma to std::lgamma because 1) we don't want to change the behavior of current programs; and 2) we are (currently) following what the std library does when using std functions. We can have both so somebody has to pick one.

@syclik is (is going to be?) the lead for the math library so it makes sense for him to decide what to do with this.

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Nov 30, 2018

Ah I understand what you mean now. The main reason that I was confused is that the current tests expect that the function throws, and so I was trying to maintain the current behaviour.

I don't think so. Here's the current version (stan/math/prim/scal/fun/acosh.hpp):

inline double acosh(double x) {
  if (unlikely(is_nan(x)))
    return x;
  else
    return boost::math::acosh(x, boost_policy_t());
}

Won't implementing the standard library math functions in this way lead to inconsistencies with the boost math functions (i.e. the boost functions will throw while the standard functions will silently fail)?

No, we use policies on the boost functions so that they don't throw. Here's the policy (stan/math/prim/scal/fun/boost_policy.hpp):

typedef boost::math::policies::policy<
    boost::math::policies::overflow_error<
        boost::math::policies::errno_on_error>,
    boost::math::policies::pole_error<boost::math::policies::errno_on_error> >
    boost_policy_t;

Throwing from Stan functions is fine (at least for now) because we can control the exception type and message.

The math library policy right now is to follow the throw behavior of the C++ standard library.

@syclik is (is going to be?) the lead for the math library so it makes sense for him to decide what to do with this.

Is. He hasn't expressed an opinion yet.

It's confusing for everyone. I don't think there's a clear answer for what to do when switching from boost::lgamma to std::lgamma because 1) we don't want to change the behavior of current programs; and 2) we are (currently) following what the std library does when using std functions. We can have both so somebody has to pick one.

Is it clear now?

@bob-carpenter
Copy link
Contributor

And just to be super explicit, here's the reverse-mode implementation:

inline var acosh(const var& a) {
  return var(new acosh_vari(stan::math::acosh(a.val()), a.vi_));
}

It calls the basic math version so it also won't throw.

@andrjohns
Copy link
Collaborator Author

I don't think so. Here's the current version (stan/math/prim/scal/fun/acosh.hpp):

It's not defined in the function header, but the throwing is expected nonetheless. This is from acosh_test.cpp in the vanilla math library:

TEST(MathFunctions, acosh_exception) {
  using stan::math::acosh;
  EXPECT_THROW(acosh(0.5), std::domain_error);
}

No, we use policies on the boost functions so that they don't throw. Here's the policy (stan/math/prim/scal/fun/boost_policy.hpp):

That policy sets the behaviour for pole and overflow errors, not domain errors. If you update the policy to add the same behaviour for domain errors:

typedef boost::math::policies::policy<
    boost::math::policies::overflow_error<
        boost::math::policies::errno_on_error>,
    boost::math::policies::pole_error<boost::math::policies::errno_on_error>,
    boost::math::policies::domain_error<boost::math::policies::errno_on_error> >
    boost_policy_t;

Then re-run the test for acosh:

test/unit/math/prim/scal/fun/acosh_test --gtest_output="xml:test/unit/math/prim/scal/fun/acosh_test.xml"
Running main() from lib/gtest_1.8.1/src/gtest_main.cc
[==========] Running 4 tests from 1 test case.
[----------] Global test environment set-up.
[----------] 4 tests from MathFunctions
[ RUN      ] MathFunctions.acosh
[       OK ] MathFunctions.acosh (0 ms)
[ RUN      ] MathFunctions.acosh_exception
test/unit/math/prim/scal/fun/acosh_test.cpp:17: Failure
Expected: acosh(0.5) throws an exception of type std::domain_error.
  Actual: it throws nothing.
[  FAILED  ] MathFunctions.acosh_exception (0 ms)
[ RUN      ] MathFunctions.acosh_inf_return
[       OK ] MathFunctions.acosh_inf_return (1 ms)
[ RUN      ] MathFunctions.acosh_nan
[       OK ] MathFunctions.acosh_nan (0 ms)
[----------] 4 tests from MathFunctions (1 ms total)

[----------] Global test environment tear-down
[==========] 4 tests from 1 test case ran. (1 ms total)
[  PASSED  ] 3 tests.
[  FAILED  ] 1 test, listed below:
[  FAILED  ] MathFunctions.acosh_exception

 1 FAILED TEST
test/unit/math/prim/scal/fun/acosh_test --gtest_output="xml:test/unit/math/prim/scal/fun/acosh_test.xml" failed

@bob-carpenter
Copy link
Contributor

Hmm. That's bad. That's a bug in the traits and the tests, in my opinion, given the goal of matching the standard library throw behavior. At the same time, I've been thinking that in Stan 3 we want to change behavior and have everything throw. As is, we're in this intermediate situation where things like sqrt and log that have been implemented in std:: all along don't throw, e.g.,

inline double sqrt(int x) { return std::sqrt(static_cast<double>(x)); }

but the functions like lgamma and acosh that were in C99 (and hence some C++ compilers but not all of them) were implemented with Boost in such a way that they didn't match the C99 (and now C++11) behavior.

The question is now what we want to do about it. Let's see if @syclik has an opinion. I'll try to ping him via email.

My inclination is to not try to go and change all of our behavior back to matching the standard library only to change it up again.

@andrjohns
Copy link
Collaborator Author

That sqrt code also brings up another point with the use of c++11 - the static_cast for integers is no longer needed. According to the spec, all integral types can be used directly with the std:: functions. Using the sqrt ref as an example:

Additional overloads are provided in this header (cmath) for the integral types: These overloads effectively cast x to a double before calculations (defined for T being any integral type).

@syclik
Copy link
Member

syclik commented Dec 2, 2018

@andrjohns, sorry about the confusion. You're digging right into code that touches a lot of edge cases. I'll try to summarize what people are asking me to decide on, although I don't think there's a decision here.

I believe what people are asking is for Stan's math functions that are implemented in the standard library, what behavior should the functions have with the error conditions?

The answer is: first follow the existing Math library's implementation. Then fall back to the standard library's implementation. Then Boost. If the implementation throws, that's fine, but we would like to throw with a reasonable message. We currently check prior to calling a function to form an error message, but we don't have to do it that way (given the code is clean, the patterns are reasonable, and there isn't some sort of performance hit by doing it that way).

Hopefully that explains it a bit. To get into motivation, it's as Bob mentioned: changing behavior at these edge cases does change the behavior of the Stan programs. We won't do that at the Math library level unless we have a good reason and we announce that we will do this. At some point we should, but I don't think this PR should be changing a lot of behavior. I'd rather roll all of that up together in one release that is behavior-changing to happen all at once. This is why testing for known edge cases matter... turns out that Boost has changed its strategy for some functions over they years. We want to be consistent independent of what the standard library implementation on different compilers do and what Boost does.

Before this PR gets merged, it really needs to be tested on Windows. I know the standard library is supposed to implement a specification, but sometimes implementations don't behave the same way in either unspecified areas or edge cases.

@bob-carpenter, does that answer what you were asking? I'm saying we should match the current Math implementation first, no matter if that's inconsistent with either the standard library or Boost (or both). And then use the standard library's rules.

@syclik
Copy link
Member

syclik commented Dec 2, 2018

That sqrt code also brings up another point with the use of c++11 - the static_cast for integers is no longer needed. According to the spec, all integral types can be used directly with the std:: functions. Using the sqrt ref as an example:

@andrjohns, feel free to remove the static_cast. As long as the behavior is the same, that's great. (We'll have to empirically verify that it's the same because doc isn't always correct.)

@bob-carpenter
Copy link
Contributor

@bob-carpenter, does that answer what you were asking? I'm saying we should match the current Math implementation first, no matter if that's inconsistent with either the standard library or Boost (or both). And then use the standard library's rules.

Yes. Thanks.

@bob-carpenter
Copy link
Contributor

That sqrt code also brings up another point with the use of c++11 - the static_cast for integers is no longer needed. According to the spec, all integral types can be used directly with the std:: functions.

Simplifying the code is great. The only downside is that it's more code changes to implement and to review.

@andrjohns
Copy link
Collaborator Author

Before this PR gets merged, it really needs to be tested on Windows

@syclik I'll run the tests locally on windows before I push anything. Is the Rtools 3.5 compiler the right target (gcc 4.9.3, I believe)? Any other compilers I should test against?

@bob-carpenter
Copy link
Contributor

Yes, the Rtools compiler should be sufficient.

@syclik --- isn't this being tested as part of continuous integration? I know there was some traffic saying Windows tests were broken, but I thought this issue was the main cause. If we don't have Windows tests turned on, that seems like a higher priority than anything else for the math lib.

@syclik
Copy link
Member

syclik commented Dec 3, 2018 via email

@syclik
Copy link
Member

syclik commented Dec 3, 2018 via email

@andrjohns
Copy link
Collaborator Author

Testing under windows was a good idea. For no reason that I can determine acosh (and only acosh) is returning a nan for an infinity input when it should be returning infinity:

int main() {
  double inf = std::numeric_limits<double>::infinity();
  std::cout << "acosh(inf) = " << std::acosh(inf) << std::endl;
  return 0;
}

On Linux:

andrew@VMLinux:~$ g++ --version
g++ (Ubuntu 8.2.0-7ubuntu1) 8.2.0

andrew@VMLinux:~$ g++ -std=c++11 test.cpp

andrew@VMLinux:~$ ./a.out 
acosh(inf) = inf

On Windows:

C:\Users\Andrew Johnson\Documents>g++ --version
g++ (x86_64-posix-seh, Built by MinGW-W64 project) 4.9.3

C:\Users\Andrew Johnson\Documents>g++ -std=c++11 test.cpp

C:\Users\Andrew Johnson\Documents>a.exe
acosh(inf) = nan

On Linux, I also tested with g++-4.8, g++-6, and g++-7 and all returned the correct output.

@sakrejda
Copy link
Contributor

sakrejda commented Dec 4, 2018

I would, on Windows only, pre-check for Inf and return Inf. I don't know that we have a policy for this sort of thing.

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Dec 4, 2018 via email

@sakrejda
Copy link
Contributor

sakrejda commented Dec 4, 2018

Our goal's to have Stan behave the same way on all platforms, if that's what you mean by policy.

I was trying to remember if platform-specific ifdef blocks were permited or if we should just add the pre-call check to all platforms.

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Dec 4, 2018 via email

@andrjohns
Copy link
Collaborator Author

So for this case should a platform-specific ifdef or a blanket (all platform) input check be used?

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Dec 5, 2018 via email

@andrjohns
Copy link
Collaborator Author

Sounds good to me, I'll make the changes and ping you when they're ready

@syclik
Copy link
Member

syclik commented Dec 6, 2018 via email

@andrjohns
Copy link
Collaborator Author

andrjohns commented Dec 6, 2018

That should be fine. I'll do one PR for the cases that are a straight swap from boost::math::foo() to stan::math::foo(), and another for the functions that need error-checking/edge-case handling

@syclik
Copy link
Member

syclik commented Dec 6, 2018 via email

@andrjohns andrjohns changed the title Feature/issue 1010 - Replace boost::math:: functions with cmath Feature/issue 1010 - Replace boost::math:: functions with cmath (Part 1) Dec 9, 2018
@andrjohns
Copy link
Collaborator Author

andrjohns commented Dec 9, 2018

@bob-carpenter This is ready for review again. To simplify the review process, the changes have been broken into three parts (covered in more detail in this comment). This pull replaces the boost functions that don't need input checking and throwing.

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Dec 9, 2018 via email

Copy link
Contributor

@bob-carpenter bob-carpenter left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. This looks great.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants