Adding the erf special function to numpy?

[normalhuman]

EDIT: the original issue here was resolved a long time ago, Intel removed their erf addition which wasn't in numpy. Issue retitled to reflect the discussion about adding np.erf

When np.erf is given a complex number, it returns the number itself instead of the value of the error function. For example, scipy.special.erf(1+2j) returns (-0.5366435657785664-5.0491437034470374j) but np.erf returns (1+2j) itself.

Reproducing code example:

import numpy as np
print(np.erf(1))     # prints 0.8427007929497149  which is correct 
print(np.erf(1+2j))   # prints (1+2j) which is incorrect

Error message:

None, but the return value is wrong for complex arguments. It is simply the argument itself. If complex inputs are not supported, an error should be thrown instead.

Numpy/Python version information:

1.15.4 3.7.1 (default, Oct 23 2018, 19:19:42) 
[GCC 7.3.0]

[eric-wieser]

np.erf is not present in any of our releases.

If it exists on your system, then someone has modified your numpy installation.

[normalhuman]

Indeed, the only other mention of np.erf I see in this issue tracker says

can you say if your numpy version has np.random_intel or np.erf defined? (it is just that there are patched numpy versions from Intel, and if yours is one we may have to push it to Intel)

What is the right tracker for such an issue, if there is one?

[mattip]

@oleksandr-pavlyk ping

[mhvk]

@normalhuman - aside: if your code is meant to be portable, I'd strongly recommend using scipy.special since other people will not have np.erf... (this is the main reason we do not think it is a good idea to inject new functions in the numpy namespace).

[oleksandr-pavlyk]

@normalhuman Yes, please use scipy.special.erf for portability. np.erf in patched NumPy from Intel was meant for real arguments only, like in computation of Black-Scholes formula.

I will make sure to disable np.erf for complex inputs.

[njsmith]

@oleksandr-pavlyk please consider removing it from the numpy namespace entirely, and moving it into a separate package namespace. Changing the API of a project like this is extremely poor etiquette, because it creates confusion among users (like we saw here), and sets us up for more breakage later. (What are you going to do if numpy does add an erf function, that's incompatible with yours?) Putting this functionality in a seperate "Intel add-ons" package avoids both issues.

[normalhuman]

@oleksandr-pavlyk For context, np.erf is automatically used by SymPy's lambdify utility which converts a SymPy object (a mathematical formula) into a Python callable: lambdify(x, erf(x), "numpy") uses np.erf if it's available, because the function name matches. Naming the function np.erf_intel, similar to np.random_intel, would avoid the issue.

[oleksandr-pavlyk]

@njsmith Would it be OK to keep numpy.core.umath.erf? Would numpy community consider adding np.erf supported for reals? This is part of cmath header in C++11.

[mhvk]

Personally, I think erf is common enough (also part of FORTRAN standard) that it makes sense to support it generally rather than point people to scipy.special. With that, it would be less of a problem to have a compiler/processor specific implementation. I'm a bit less sure about only supporting real, when there is an implementation in scipy.special that also supports complex. That said, the complex one does not seem to be part of the C++11 or FORTRAN standards.

Let me ping @rgommers to see what he thinks about moving/duplicating scipy.special.erf in numpy.

[rkern]

The line that we currently have drawn for special functions in the numpy namespace is "is it in C99?" I'm not opposed to updating that, but I'd like to draw a similar line referencing some other specific external standard and be prepared to add everything new in that standard at once, not piecemeal.

Disclosure: I work for Enthought, which is also in the business of distributing optimized builds of numpy, so please consider this next recommendation in that light. However, I promise that it is offered in good faith, only informed by our past mistakes when we overzealously added additional functionality to numpy and had to revert them because they caused problems.

@oleksandr-pavlyk Please don't add stuff to numpy.core.umath, either. The same arguments apply. When we add one, it might work incompatibly to yours. Also, code that harvests the available ufuncs, like SymPy's lambdify, would be perfectly reasonable to look directly at numpy.core.umath and find yours there.

[rgommers]

Let me ping @rgommers to see what he thinks about moving/duplicating scipy.special.erf in numpy.

I agree with @rkern's thoughts about scope.

In addition, in this particular case it doesn't seem like effort well spent, since there'll not be many use cases where erf is needed but people really need to avoid a dependency on SciPy.

like SymPy's lambdify, would be perfectly reasonable to look directly at numpy.core.umath and find yours there.

this I'm not sure about - I know SymPy and other libs do dig in there, but core.umath isn't public I think (at least in our Python API, in our C API it is - see e.g. the docs and NEP 15) and we're free to change it in principle.

[njsmith]

The line that we currently have drawn for special functions in the numpy namespace is "is it in C99?"

I don't think that's correct... do you mean C89? Traditionally we've supported non-C99 compilers (like old MSVC), and C99 does have erf (plus lots of other things we don't – erfc, gamma and log-gamma, etc.).

[mhvk]

@rkern - agreed with your point that we should support what a particular standard provides. Our documentation indeed seems to suggest we support C99, which, as @njsmith notes includes erf, erfc, lgamma and tgamma [1]. So, perhaps that argues that we should include those?

p.s. Nevertheless, agreed with all here that one should not insert any such functions in any of the numpy namespaces without there being a generic implementation; astropy does inspect np.core.umath and was bitten by this already (which led to a very hard to debug problem since the error could not be reproduced locally; see astropy/astropy#7058 (comment))

[1] https://en.wikipedia.org/wiki/C_mathematical_functions

[rkern]

Sorry, our policy has more or less been "if it's in C99, we'll let it in if someone wants to do the work; otherwise, keep it in scipy.special". So strike my request to upgrade everything at once. Piecemeal is how we do this, apparently. So, "erf is in C99" is a perfectly good reason to ask for it, but someone still has to do the work.

Our support for C99 compilers is neither here nor there. We don't rely on the implementation being provided by the C99 compiler. This was just an arbitrary, objective policy for deciding which special functions should get into numpy instead of scipy.special so we didn't have to have subjective arguments about which special functions are important or popularly used.

[rkern]

Basically, the policy takes a negative form: "if a special function is not in C99, it doesn't go into numpy".

[seberg]

Adding erf sounds fine, although I am not 100% sure that we should add it if we do not (or cannot) add the complex versions as well (assuming there is not really a precedence for that already).

Of course it would slow down development and be more involved. But I think it should be reasonable if Intel puts in special code that will be used instead of the default one when compiling with the hard dependency of MKL. Similar to SIMD inner loops that I think we have to replace others when available. That is assuming that it does not add a huge amount of new complexity probably.

[pv]

real erf is very simple, complex is much more complex to get right (the impl in scipy is in c++ which makes copypasting it to numpy harder)

…

On December 28, 2018 1:47:16 PM UTC, Sebastian Berg ***@***.***> wrote: Adding `erf` sounds fine, although I am not 100% sure that we should add it if we do not (or cannot) add the complex versions as well (assuming there is not really a precedence for that already). Of course it would slow down development and be more involved. But I think it should be reasonable if Intel puts in special code that will be used instead of the default one when compiling with the hard dependency of MKL. Similar to SIMD inner loops that I think we have to replace others when available. That is assuming that it does not add a huge amount of new complexity probably.

[charris]

real erf is very simple, complex is much more complex

True. What we have done in the past for such functions is copy from FreeBSD, much of whose implementation seems to have come from Silicon Graphics back in the day. I don't know if that would work here...

[mhvk]

I think just having real erf would already be useful and solve the problem at hand. It is also most in the spirit of only supporting what the standard C libraries support. It obviously would need to include a pointer to scipy.special for the complex case.