MAINT: Translate binom to C++ #19471

steppi · 2023-11-03T17:58:25Z

Reference issue

This PR closes no issues, but continues the work discussed in the RFC issue #19404.

What does this implement/fix?

This PR translates the real valued binomial coefficient function binom from Cython into C++. Again, the translation is as literal as possible.

Additional information

Since binom depends on cephes I needed to work out how to include cephes functions in C++ code. This was not immediately possible because cephes has a header cephes_names.h which defines aliases of the form beta -> cephes_beta etc. So far so good, but because of how they are implemented, the aliases bleed into all code in the same translation unit, for instance, expanding expm1 when it appears in boost to cephes_expm1. Rather than trying to sort that out directly, I've created a new header cephes.hh for inclusion in C++ files, which puts cephes functions in the scipy::special::cephes namespace without a prefix in the name, and undefs all of the aliases

For now, I've kept the binom implementation within orthogonal_eval.pxd and it is still used internally there.

As is, this binom implementation is not header only because it depends on cephes. We'll need to discuss the possibility of making cephes or parts of cephes header only in the future.

This PR should unblock the work that is stalled on gh-19287.

cc @izaid, @tirthasheshpatel, @rlucas7, @ilayn

steppi · 2023-11-03T19:28:36Z

The unrelated failure seems to be due to the recent NumPy integer dtype changes. #19466, #19462, numpy/numpy#24224.

rlucas7 · 2023-11-04T00:51:22Z

The unrelated failure seems to be due to the recent NumPy integer dtype changes. #19466, #19462, numpy/numpy#24224.

it seems so yeah, there's a discussion on #19466 that seems to (maybe?) address the issue for sobol.

I can take a closer look this weekend once I update #19287 to get merged lambert code incorporated. Thanks for getting this up so quickly @steppi

ev-br · 2023-11-04T08:14:50Z

Let me devil advocate a bit. Since we pull stuff from boost already, can https://beta.boost.org/doc/libs/1_37_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/factorials/sf_binomial.html

replace our implementation?

steppi · 2023-11-04T11:29:13Z

Let me devil advocate a bit. Since we pull stuff from boost already, can https://beta.boost.org/doc/libs/1_37_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/factorials/sf_binomial.html

replace our implementation?

I’d prefer to be able to use boost if we could. We’re avoiding boost for now until they make the changes necessary to get their stuff to run on Cuda.

izaid · 2023-11-04T12:42:00Z

scipy/special/_binom.h

+#include "cephes.hh"
+
+
+using std::numeric_limits;


I don't think we should import things like this in header files. Best to keep it fully qualified.

If you really want to do using std::numeric_limits;, you can actually do that within either the functions or namespaces it is used in. But I'd stick with not doing it.

Ah, good catch. Agreed on not wanting to have unexpected side effects when including a header. There are some imports in _lambertw.h that can be removed in another PR.

scipy/scipy/special/_lambertw.h

Lines 32 to 34 in 4d790c1

using namespace std::complex_literals;

using std::numeric_limits;

scipy/special/_binom.h

izaid · 2023-11-04T12:43:16Z

scipy/special/_binom.h

+	    }
+	    if (k > 1E8*std::fabs(n)) {
+		// avoid loss of precision
+		num = cephes::Gamma(1 + n) / std::fabs(k)


Curious why is it Gamma with a capital 'G'?

No clue, but cephes gamma function has been capital G in SciPy from the start: 52c64a9.

izaid · 2023-11-04T12:43:48Z

scipy/special/_extra_special.h

 #include "_lambertw.h"

+inline double binom(double n, double k) {


Do we want this to be npy_double or just double?

I don't think it matters for double. It seems like the convention in other parts of the code base is to just have double. Check out specfun_wrappers for instance.

scipy/special/cephes.hh

izaid · 2023-11-04T12:46:05Z

scipy/special/cephes.hh

+#undef kolmogp
+#undef kolmogc
+#undef kolmogci
+#undef owens_t


Yeah... I've actually messed about with cephes in SciPy before. It has some things going on that are definitely not best practices. Our long term aim should be to not use it like this, and just adapt it to current best practices for C / C++.

Yeah, I plan to make another mailing list post and RFC issue about what to do with cephes.

izaid · 2023-11-04T12:49:00Z

@steppi Just read through and left some comments! Generally looks quite good, well done. A few things I think we should change.

Yeah, we'll need to do something to integrate cephes a bit better. I don't think it's so hard, more annoying, but what we have now is certainly not ideal.

rlucas7

This passes all tests for me in stats and special locally, and has a +1 from @izaid so merging to unblock my work on stirling2 Thanks @steppi + @izaid

steppi · 2023-11-12T00:48:58Z

Thanks @rlucas7!

steppi added 3 commits November 3, 2023 12:52

Translate binomial coefficient to C++

f44db2b

Fix typos in use of floor

940e434

Avoid polluting global namespace when include cephes.hh

2f449a5

steppi requested a review from person142 as a code owner November 3, 2023 17:58

steppi added scipy.special C/C++ Items related to the internal C/C++ code base maintenance Items related to regular maintenance tasks labels Nov 3, 2023

izaid reviewed Nov 4, 2023

View reviewed changes

scipy/special/_binom.h Outdated Show resolved Hide resolved

izaid reviewed Nov 4, 2023

View reviewed changes

scipy/special/cephes.hh Outdated Show resolved Hide resolved

izaid reviewed Nov 4, 2023

View reviewed changes

dschmitz89 mentioned this pull request Nov 4, 2023

ENH: special: add ufunc for logarithm of the binomial coefficient #19473

Open

steppi added 3 commits November 4, 2023 10:15

Use more standard namespace name

7eabaf1

Remove using statement which shouldn't be in header

26d6730

Remove unnecessary else

dbb7265

steppi mentioned this pull request Nov 6, 2023

RFC: Thoughts on making our vendored Cephes a cleaner header only library. #19481

Closed

izaid approved these changes Nov 6, 2023

View reviewed changes

rlucas7 approved these changes Nov 11, 2023

View reviewed changes

rlucas7 merged commit bf819ba into scipy:main Nov 11, 2023
21 of 23 checks passed

j-bowhay added this to the 1.12.0 milestone Nov 12, 2023

rlucas7 mentioned this pull request Nov 14, 2023

MAINT: Translate loggamma from Cython into C++ #19496

Merged

steppi mentioned this pull request Dec 12, 2023

MAINT: Add binom to new C++ special lib along with its cephes dependencies #19681

Merged

steppi mentioned this pull request Jan 8, 2024

MAINT: Translate Digamma from Cython into C++ #19834

Merged

steppi deleted the binom_cpp branch March 7, 2024 19:33

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

MAINT: Translate binom to C++ #19471

MAINT: Translate binom to C++ #19471

steppi commented Nov 3, 2023

steppi commented Nov 3, 2023 •

edited

rlucas7 commented Nov 4, 2023

ev-br commented Nov 4, 2023

steppi commented Nov 4, 2023

izaid Nov 4, 2023

steppi Nov 4, 2023

izaid Nov 4, 2023

steppi Nov 4, 2023

izaid Nov 4, 2023

steppi Nov 4, 2023

izaid Nov 4, 2023

steppi Nov 4, 2023

izaid commented Nov 4, 2023

rlucas7 left a comment

steppi commented Nov 12, 2023

	using namespace std::complex_literals;

	using std::numeric_limits;

		#include "_lambertw.h"

		inline double binom(double n, double k) {

MAINT: Translate binom to C++ #19471

MAINT: Translate binom to C++ #19471

Conversation

steppi commented Nov 3, 2023

Reference issue

What does this implement/fix?

Additional information

steppi commented Nov 3, 2023 • edited

rlucas7 commented Nov 4, 2023

ev-br commented Nov 4, 2023

steppi commented Nov 4, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

izaid commented Nov 4, 2023

rlucas7 left a comment

Choose a reason for hiding this comment

steppi commented Nov 12, 2023

steppi commented Nov 3, 2023 •

edited