BUG: ufuncs that call into Python do not ignore FPEs e.g. for inf - inf

[Steven-R-Hall]

Describe the issue:

I'm writing a program that uses the casadi module to automatically differentiate Python expressions. The program constructs a numpy ndarray (with dtype='O') of objects with the casadi SX data type. Say the array is called 'x'. When performing arithmetic operations on the array, sometimes a Numpy warning is generated. For example,

>>> 2**31  /  x

generates the warning

RuntimeWarning: invalid value encountered in true_divide

and

>>> 2**31 + x

generates the warning

RuntimeWarning: invalid value encountered in add

Despite the warning, the correct value is produced. However, the following code runs without warning:

>>> (2**31-1)  /  x
>>> (2**31-1)  +  x

My expectation is that

>>> numpy.true_divide(y, x)

should run without a warning whenever y is a float or int, x is a numpy.ndarray with dtype='O', and y/x[i] can be calculated without warning for all values of i in range. But that is not what happens here. Curiously, it fails only when the integer or float is at least 2 **31, even though no such behavior exists in casadi (or at least not obviously).

The sample code below is the simplest code that demonstrates the issue. It seems that the line calculating y2 should never generate a warning if the line calculating y1 does not.

Reproduce the code example:

import numpy as np
import casadi as cd

x = np.zeros([1], dtype='O')
x[0] = cd.SX.sym('r')  # only element of array x is a casadi symbol

y1 = 2 ** 31 / x[0]  # no warning or error generated
y2 = 2 ** 31 / x  # generates warning, but result should be np.array([y1], dtype='O')
y3 = (2 ** 31 - 1) / x  # ... but this works fine with no warning or error

Error message:

true_divide_problem.py:8: RuntimeWarning: invalid value encountered in true_divide
  y2 = 2 ** 31 / x  # generates warning, but result should be np.array([y1], dtype='O')

NumPy/Python version information:

1.21.4 3.9.9 (main, Nov 30 2021, 10:15:14)

[seberg]

This is a bit confusing, because the warning seems to be caused by casadi and not NumPy, but casadi would ignore it (not the Python warning, but the CPU one).
The way floating point warnings work, is that the CPU has flags for these warnings that can be read out (or "raised"). Most Python code will just ignore them completely, but NumPy does not.

For example, code like this:

def myfunc(x):
    return x - x

myfunc(np.inf)  # does not warn

import numpy as np
vec_myfunc = np.frompyfunc(myfunc, nin=1, nout=1)

vec_myfunc(np.inf)  # does warn

In that case np.inf - np.inf warns (on the CPU/instruction level as it should as per IEEE). But Python ignores the warning. When the np.inf - np.inf is "wrapped" into a ufunc though, NumPy does not ignore it anymore.

The same seems to happen here, casadi runs into such an "invalid" operation (i.e. a NaN is created internally somewhere). But casadi normally just ignores it (which is probably correct). NumPy picks it up anyway though.

NumPy could work around this by storing and restoring those FPE (floating point exception) flags on the CPU. Considering that Python calls are not super fast, that would probably not be a problem. casadi could probably avoid creating a NaN internally here also.

[Steven-R-Hall]

It appears you are right. Here is my reworked minimal code that generates the issue:

import numpy as np
import casadi as cd

def myfunc(x):
    return x / cd.SX.sym('r')

vec_myfunc = np.frompyfunc(myfunc, nin=1, nout=1)

y1 = 2 ** 31 / cd.SX.sym('r')  # does not warn
y2 = vec_myfunc(2 ** 31)  # does warn

which generates

true_divide_problem_2.py:10: RuntimeWarning: invalid value encountered in myfunc (vectorized)
y2 = vec_myfunc(2 ** 31)  # does warn

Just to understand the problem further,

1. How does numpy generate an error in your example? myfunc doesn't generate a warning in the python interpreter, so how is one generated in numpy? I suppose that's a question about how np.frompyfunc vectorizes functions.

2. What is the right way to set warnings so that vec_myfunc generates a warning if and only if my_func would generate a warning for the corresponding scalar input? I assume that that would be

   np.seterr(invalid='ignore')

So originally I labeled this issue as "Unexpected behavior of ndarray with dtype='O'". I guess after looking into it a bit it is consistent with numpy trying to be more precise about handling floating point numbers consistently with respect to the IEEE standard. The behavior still seems a little odd, though, because in this case the result isn't expected to be an int or a float, but just some object, which could be any valid python object. I have no idea whats going on inside casadi (the casadi module is a wrapper around C++ code), but in the test case above there's nothing that should generate a warning. Indeed, as far as I can tell, no floating point operation at all should be generated by the example, since the true_divide in myfunc is not an arithmetic divide at all, but a symbolic operation.

[seberg]

The NumPy ufunc machinery does this, since it has to do it especially for Numeric arrays. In principle NumPy could skip doing it here (because it should never happen), but that wouldn't even the absolute correct fix.
The absolute correct thing NumPy could do is explicitly ignore the warnings here if they originate in Python code, but not if they originate elsewhere (which may be a bit of a weird back and forth).

Using np.errstate() is the correct way to silence these warnings locally, yes.

As for casadi, there are two possible causes:

• Most likely, there is just some FPE triggering operation somewhere. That might even be something like the 2**32 being converted to float32 and then back again to an integer...
• The compiler is doing optimizations that lead to sloppy handling (because a lot of programs just don't care).

Since casadi is a library, they might want to avoid whatever triggers it.

[seberg]

I have tried a workaround for it:

diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
index 3a8a54913..7115ccd67 100644
--- a/numpy/core/src/umath/loops.c.src
+++ b/numpy/core/src/umath/loops.c.src
@@ -206,27 +206,37 @@ PyUFunc_FF_F_As_DD_D(char **args, npy_intp const *dimensions, npy_intp const *st
  **                         GENERIC OBJECT lOOPS                             **
  *****************************************************************************/
 
+#include <fenv.h>
+
 /*UFUNC_API*/
 NPY_NO_EXPORT void
 PyUFunc_O_O(char **args, npy_intp const *dimensions, npy_intp const *steps, void *func)
 {
+    fenv_t fenv;
+    fegetenv(&fenv);
+
     unaryfunc f = (unaryfunc)func;
     UNARY_LOOP {
         PyObject *in1 = *(PyObject **)ip1;
         PyObject **out = (PyObject **)op1;
         PyObject *ret = f(in1 ? in1 : Py_None);
         if (ret == NULL) {
+            fesetenv(&fenv);
             return;
         }
         Py_XDECREF(*out);
         *out = ret;
     }
+    fesetenv(&fenv);
 }

In most cases, the speed difference is neglegible due to buffering. The most extreme case I found was:

arr = np.ones(10000000, dtype=object)
where=np.ones_like(arr, dtype=bool)
where[::2] = False

%timeit np.absolute(arr, where=where)

Changing from

277 ms ± 2.42 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

on the current version to:

1.37 s ± 1.17 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

with the fgetenv dance on each element (where a calculation is done). Floating point environment mutation is not super quick unfortunately...

EDIT: Of course the above is probably the most extreme case possible by a very far amount, so it may well be OK in practice!

EDIT2: There may also be speedups available with some dance to avoid setting the flags when they are not used anyway?

[seberg]

The following diff is only ~10% unless warnings are actually being silenced. So this is still slow:

arr = np.full(10000000, fill_value=np.inf, dtype=object)
where=np.ones_like(arr, dtype=bool)
where[::2] = False

%timeit np.subtract(arr, arr, where=where)

but checking errors is much faster than setting them, so doing it only if necessary really saves a lot of time.

New improved diff

diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
index 3a8a54913..19e0f6314 100644
--- a/numpy/core/src/umath/loops.c.src
+++ b/numpy/core/src/umath/loops.c.src
@@ -206,21 +206,31 @@ PyUFunc_FF_F_As_DD_D(char **args, npy_intp const *dimensions, npy_intp const *st
  **                         GENERIC OBJECT lOOPS                             **
  *****************************************************************************/
 
+#include <fenv.h>
+
 /*UFUNC_API*/
 NPY_NO_EXPORT void
 PyUFunc_O_O(char **args, npy_intp const *dimensions, npy_intp const *steps, void *func)
 {
+    fexcept_t excepts, after;
+    fegetexceptflag(&excepts, FE_ALL_EXCEPT);
+
     unaryfunc f = (unaryfunc)func;
     UNARY_LOOP {
         PyObject *in1 = *(PyObject **)ip1;
         PyObject **out = (PyObject **)op1;
         PyObject *ret = f(in1 ? in1 : Py_None);
         if (ret == NULL) {
-            return;
+            goto finish;
         }
         Py_XDECREF(*out);
         *out = ret;
     }
+  finish:
+    fegetexceptflag(&after, FE_ALL_EXCEPT);
+    if (after != excepts) {
+        fesetexceptflag(&excepts, FE_ALL_EXCEPT);
+    }
 }
 
 /*UFUNC_API*/
@@ -263,6 +273,9 @@ PyUFunc_O_O_method(char **args, npy_intp const *dimensions, npy_intp const *step
 NPY_NO_EXPORT void
 PyUFunc_OO_O(char **args, npy_intp const *dimensions, npy_intp const *steps, void *func)
 {
+    fexcept_t excepts;
+    fegetexceptflag(&excepts, FE_ALL_EXCEPT);
+
     binaryfunc f = (binaryfunc)func;
     BINARY_LOOP {
         PyObject *in1 = *(PyObject **)ip1;
@@ -270,11 +283,14 @@ PyUFunc_OO_O(char **args, npy_intp const *dimensions, npy_intp const *steps, voi
         PyObject **out = (PyObject **)op1;
         PyObject *ret = f(in1 ? in1 : Py_None, in2 ? in2 : Py_None);
         if (ret == NULL) {
-            return;
+            goto finish;
         }
         Py_XDECREF(*out);
         *out = ret;
     }
+
+  finish:
+    fesetexceptflag(&excepts, FE_ALL_EXCEPT);
 }
 
 NPY_NO_EXPORT void

[Steven-R-Hall]

@seberg, thank you for this. Now that I understand the situation, I'm not sure whether this is a bug or not, and whether all Numpy users should see a slowdown. As you noted, there must be something going on in the casadi code that generates flags which casadi ignores.

In my use case, the arrays are of modest size (maybe 10000 elements or fewer), and are operated on only a few times (order 10s). So the slowdown would be relatively unimportant for me. But for others, a 10% hit might be significant.

[seberg]

I it a bug, but a very low priority one. Arguably, it almost never gets triggered (and if it does, like here, the triggering side is also doing something non-ideal).

Doing the fix I said is a simple judgment call to me. Getting it 100% right is probably not going to happen (at least not soon), it is just too complex. Whether getting it more right here is worth it, I don't have a strong opinion on.

In practice, the slowdown should not be noticable in the vast majority of cases and serious slowdowns will only happen in some extreme cases (errors require silencing on very small chunked/unbuffered inputs, which probably means the use of where= which is super rare in itself).
OTOH, it also makes the code more complex for very little gain ;).

[Steven-R-Hall]

The reason I said it's maybe not a bug is this: if myfunc is pure python, I think the Numpy philosophy dictates that if CPU flags are set there should be a warning, even though Python doesn't warn, which is what it does. Even in my case, myfunc is a mixture of Python and an extension, so I don't think Numpy even has the ability to tell whether the flags are the result of a pure Python invalid result, or a byproduct of an extension where the invalid flag means something entirely different (perhaps).

I'll leave it to the Numpy team to figure this out, but thanks again so much for the quick and informative responses.

[seberg]

We discussed this today a bit in the triage meeting, and we agreed that it seems just not worth the hassle. The number of people running into this seems low to non-existing.

Let me close the issue, if someone runs into the issue we can reconsider. I do think in principle it would be cleaner to try to clean it up, but the win just seems super small.

[wshanks]

I am not asking for any action on this, but I wanted to note a context where this effect came up for me. In the uncertainties package, there is a function that calls numpy.vectorize on a function that does a < 0 comparison to its argument. When the argument is np.nan, this triggers a floating point exception and produces the numpy invalid value warning. This case was covered in more detail in this stackoverflow question. I note it here for the purpose of tracking if this issue comes up enough to do something about because I think comparison to numpy.nan inside of numpy.vectorize might be a common case where this warning shows up.