BUG: arr.fill() does not work with USE_SETITEM scalars

[proyan]

TLDR:

When creating a user-defined data-type using boost-python, np.fill() uses a memory manipulation scheme which doesn't correspond to the PyObject ABI created by boost-python. On the other hand, direct assignment using mat[:] accesses the data using PyArray_Data, and returns the correct pointer.
As a result, there is garbage value assignment in fill(), and sane value assignment in direct access

Detailed:

Hello NumPy developers,
First of all, I want to thank you all for your hardwork.
We are a team of robotics engineers/researchers, and one of our packages includes extension of the Eigen3 library to python using numpy.ndarray as a container. (see eigenpy library)

One of the recent features proposed in this library is the definition of user-types as numpy scalars, (https://github.com/jcarpent/eigenpy/blob/devel/include/eigenpy/user-type.hpp#L235). The unit-test for this feature, (https://github.com/jcarpent/eigenpy/blob/devel/unittest/python/test_user_type.py), does not correctly perform some allocation in the fill() instruction, and thus the test fails at

  mat.fill(mat.dtype.type(20.))

On the other hand,

  mat[:] = mat.dtype.type(20.)

succeeds without any issue.
I have tried to compare the behaviour of the mat.fill() command with mat[:] and the numpy user-defined dtype numpy.core._rational_tests.rational in gdb, to understand better how the wrong address is being propagated in the following sample script:

import user_type
import numpy as np
from numpy.core._rational_tests import rational
rows = 1
cols = 1

mat = user_type.create_double(rows,cols)
a = user_type.CustomDouble(10.)
x = rational(1)
y = np.array([x])
print("memory address of a,",hex(id(a)))
print("memory address of a.data tuple,",hex(id(a.__array_interface__["data"])))
print("memory address of a.__array_struct__,",a.__array_struct__)
print("memory address of a.data, ,",hex(a.__array_interface__["data"][0]) )
print("memory address of CustomDouble,",hex(id(np.dtype(a))))
print("memory address of mat,",hex(id(mat)))
print("memory address of mat.data tuple,",hex(id(mat.__array_interface__["data"])))
print("memory address of mat.data ,",hex(mat.__array_interface__["data"][0]))
mat[:] = a //First test
mat.fill(a) //Second test
y.fill(x)      // numpy implementation

The output of this script is:

memory address of a, 0x7ffff71b8c20
memory address of a.data tuple, 0x7ffff7142580
memory address of a.__array_struct__, <capsule object NULL at 0x7ffff73132a0>
memory address of a.data, , 0xd996c0
memory address of CustomDouble, 0xb4bb20
memory address of mat, 0x7ffff11676f0
memory address of mat.data tuple, 0x7ffff7142580
memory address of mat.data , 0xd36050

In all three tests, source variables (a and x) are being copied to destination variables (mat and y)

In the sample run that I did,

First Test

mat[:] = a

The variable a is copied to the destination pointer through the copyswapn function defined here, with the arguments

Breakpoint 3, eigenpy::internal::SpecialMethods<CustomType<double>, 256>::copyswapn (dst=0xd36050, dstride=8, src=0xd996c0, sstride=0, n=1, swap=0, array=0x7ffff0257b70)
    at /home/rbudhira/devel/src/sot/eigenpy/include/eigenpy/user-type.hpp:136

The backtrace shows that the correct data pointer inside a is accessed through PyArray_Data here:

numpy/numpy/core/src/multiarray/array_assign_array.c

Line 272 in 3409bd3

return PyArray_AssignRawScalar(

This test thus passes.

Second Test

mat.fill(a)

mat.fill however, is not able to correctly acquire the pointer for the location of the source data. The arguments given to copyswapn function are:

Breakpoint 3, eigenpy::internal::SpecialMethods<CustomType<double>, 256>::copyswapn (dst=0xd36050, dstride=8, src=0x7ffff71b8c30, sstride=0, n=1, swap=0, array=0x7ffff0257a50)
    at /home/rbudhira/devel/src/sot/eigenpy/include/eigenpy/user-type.hpp:136

which don't correspond to the correct location of the source data.
The backtrace shows that this data pointer is created here:

numpy/numpy/core/src/multiarray/scalarapi.c

Line 171 in 3409bd3

memloc = (uintptr_t)scalar;

which assumes that the data begins immediately after PyObject_HEAD. This is obviously not the case.

This command fails.

Third test:

y.fill(x)

The numpy example data-type rational works, since it is following the ABI that is assumed by np.fill(). In this case, the data pointer created by both PyArray_Data and memory manipulation matches.

Surprisingly, fillwithscalar functional, which is defined in both our custom user type, and numpy's rational type, is never used.

Could you please explain how boost-python classes could be used as numpy user-defined data types, and how to fix the support of ndarray.fill?

Thanks,Rohan

[proyan]

This issue corresponds to the user-type fix stack-of-tasks/eigenpy#240

[seberg]

The backtrace shows that this data pointer is created here: memloc = (uintptr_t)scalar;

I don't understand, that seems correct and continues calculating the right offset, or is that logic wrong? Can you first check that descr.alignment is actually correct. You are setting here:
https://github.com/jcarpent/eigenpy/blob/devel/include/eigenpy/user-type.hpp#L238 and that offset looks a lot like it may include the PyObject * head, which would be incorrect

Surprisingly, fillwithscalar functional, which is defined in both our custom user type, and numpy's rational type, is never used.

Lucky us :)! one less weird function to worry about...

---

That said... the fill function does like a train-wreck that does not properly support USE_GETITEM and USE_SETITEM. I personally, think the recommendation to subclass np.generic is not all that helpful. (I get that in theory, it allows us to directly access the memory, but in practice it feels like we need a better solution...). Anyway, that probably does not help you much :/

One thing, that I must make sure though: Are you aware of my work on NEP 40-43? I know it is a bit slower than we hoped, but it is getting close and the last big steps are about to be finalized (there are quite a few small ones after that, but still). If you are interested, I would be happy to talk about it, or show you how to experiment with the new DType system.

[proyan]

The memory manipulation that is being carried out is:

memloc = (uintptr_t(scalar) + sizeof(PyObject) + align -1)/align)*align

Our CustomDouble

alignment = 8
sizeof(PyObject) = 16

Input: scalar: 0x7ffff71b8c20
Output: (void*)memloc: 0x7ffff71b8c30 = (scalar + 16)

Numpy rational:

alignment = 4
sizeof(PyObject) = 16
Input: scalar: 0x7f9696e45710
Output: (void*)memloc: 0x7f9696e45720 = (scalar+16)

In both cases, the memory manipulation is working in the same manner. However, the desired location for the data pointer in our CustomDouble is 0xd996c0, which is accessed through PyArray_Data when doing mat[:]=a , and not through memory manipulation. I am guessing this might be because of a different ABI used by boost python structures, rather than the PyRational structure that was defined in _rational_tests.c

[proyan]

think the recommendation to subclass np.generic is not all that helpful. (I get that in theory, it allows us to directly access the memory, but in practice it feels like we need a better solution...). Anyway, that probably does not help you much :/

Without inheriting from PyGenericArrType_Type, the command np.array([user_type.CustomDouble(10.)]) creates an array with dtype=np.object_ rather than an array with dtype=CustomDouble.

However, since the Generic implementation is not used in such a case, the fill and [:] assignments both seem to work!

If you could propose a way to register the new user-type as a scalar, without using np.generic, we might be able to avoid the check here

numpy/numpy/core/src/multiarray/convert.c

Line 381 in c3d1107

else if (PyArray_IsScalar(obj, Generic)) {

and proceed to fill the scalar properly

One thing, that I must make sure though: Are you aware of my work on NEP 40-43? I know it is a bit slower than we hoped, but it is getting close and the last big steps are about to be finalized (there are quite a few small ones after that, but still). If you are interested, I would be happy to talk about it, or show you how to experiment with the new DType system.

I would be very happy to talk about extending our work using the new features, would it work with the currently available numpy releases as well?

[seberg]

Without inheriting from PyGenericArrType_Type, the command np.array([user_type.CustomDouble(10.)]) creates an array with dtype=np.object_ rather than an array with dtype=CustomDouble.

Good point, the new API does this correctly, but the old API doesn't have the right hooks. (This is even inside NumPy already, but not public API...)

However, the desired location for the data pointer in our CustomDouble is 0xd996c0

I don't understand how this is plausible? That memory location seems completely different, does your object store a pointer to the memory? PyArrayGeneric_Type assumes that the memory address is that of an aligned C-struct starting with the PyObject. If that is not the case you simply can't inherit from np.generic.

would it work with the currently available numpy releases as well?

No, probably next release, although realistically in some "experimental" form. But getting buy-in might help make sure it happens more smoothly...

[seberg]

@proyan or am I just misunderstanding here: Your dtype actually does set USE_GETITEM and USE_SETITEM (can you see what dtype.flags is, if you are not sure)? And the problem is that arr.fill() does not honor that contract...

[proyan]

I don't understand how this is plausible? That memory location seems completely different, does your object store a pointer to the memory? PyArrayGeneric_Type assumes that the memory address is that of an aligned C-struct starting with the PyObject. If that is not the case you simply can't inherit from np.generic.

I see. I assumed as much from the addresses that I obtained. The python class is defined using boost::python here:
https://github.com/jcarpent/eigenpy/blob/2edf3ead8931af8401c7637a12ed2d699dedc909/unittest/user_type.cpp#L172
and
and it is being registered as a new dtype here:
https://github.com/jcarpent/eigenpy/blob/devel/include/eigenpy/user-type.hpp#L235

---

I may be mistaken, but I am guessing that boost python uses the structure of PyArrayObject instead of an aligned C-struct. But perhaps you would know better.

---

Your dtype actually does set USE_GETITEM and USE_SETITEM (can you see what dtype.flags is, if you are not sure)? And the problem is that arr.fill() does not honor that contract...

Yes, we define getitem, setitem, copyswap and copyswapn.

However, when using mat.fill(), setitem is not used. After PyArray_FillWithScalar here:

numpy/numpy/core/src/multiarray/convert.c

Line 381 in c3d1107

else if (PyArray_IsScalar(obj, Generic)) {

and a few other calls, the control passes to:

https://github.com/jcarpent/eigenpy/blob/2edf3ead8931af8401c7637a12ed2d699dedc909/include/eigenpy/user-type.hpp#L133

with the wrong source pointer address.

---

can you see what dtype.flags is, if you are not sure)?

In [1]: mat.dtype.flags                                                                                                                                                                                            
Out[1]: 112

[proyan]

By the way, thanks a lot for the quick responses! It is quite refreshing to see real-time activity on such a giant project.

[seberg]

Yes, we define getitem, setitem, copyswap and copyswapn.

The important part here is the flags, and they indicate NPY_USE_GETITEM and NPY_USE_SETITEM. The typical code paths never call scalar_value if the NPY_USE_SETITEM and NPY_USE_GETITEM are set. arr.fill() and its C-version fails to check this, however.

So it is indeed a bug in PyArray_FillWithScalar, which should probably just use the new PyArray_Pack, except that might modify some crazy corner cases., but maybe we can get away with it...

[proyan]

Currently, there are two ways forward that I have seen

Not inheriting from PyGenericArrType_Type

Array creation returns array of dtype object

In [14]: a = np.array([user_type.CustomDouble(10.)])                                                                                                                                                               
In [15]: a                                                                                                                                                                                                         
Out[15]: 
array([value: 10], dtype=object)

But fill works

Inheriting from PyGenericArrType_Type

Array creation returns array of correct dtype

In [14]: a = np.array([user_type.CustomDouble(10.)])                                                                                                                                                               
In [15]: a                                                                                                                                                                                                         
Out[15]: 
array([value: 10], dtype=CustomDouble)

But fill doesn't work

---

Do you see a way out in either of the above two lanes? Or perhaps a third workaround that I might use while ensuring compatibility with the existing numpy releases?

[seberg]

Note that array creation will work as long as you use dtype=customdtype.

Do you see a way out in either of the above two lanes? Or perhaps a third workaround that I might use while ensuring compatibility with the existing numpy releases?

No, I don't see a way right now. Aside from making sure the data is stored directly on the scalar in a C-struct manner (because, why not?). arr.fill is simply buggy, you would have to fix NumPy or live with arr.fill segfaulting until the bug is fixed. The current custom dtype API somewhat works, but it has a lot of weird limitation, there is a reason for replacing it...

There is internal API to allow customizing this in the future, but it is not public and intentionally a bit limited.

[proyan]

Thanks for the clarification and the quick response @seberg.

[seberg]

I really would like to apply the diff below. But it would change two things: First arr.fill(array) could in some weird cases behave differently (mainly unpacking object arrays, but that is an easy special case to add). Some np.float64(np.nan) fills and integers will error start erroring out, and that tiny test change with the void dtype shows that there is another small change there...

(Basically, this would align arr.fill(value) with arr1d[0] = value, except modiyfing all elements.

It is also interesting, that the current code actually uses safe casting for the case of array.fill(other_array), but otherwise always uses force casting.

diff --git a/numpy/core/src/multiarray/convert.c b/numpy/core/src/multiarray/convert.c
index 29a2bb0e8..30af810d7 100644
--- a/numpy/core/src/multiarray/convert.c
+++ b/numpy/core/src/multiarray/convert.c
@@ -19,6 +19,7 @@
 #include "array_assign.h"
 
 #include "convert.h"
+#include "array_coercion.h"
 
 int
 fallocate(int fd, int mode, off_t offset, off_t len);
@@ -357,151 +358,28 @@ PyArray_ToString(PyArrayObject *self, NPY_ORDER order)
 NPY_NO_EXPORT int
 PyArray_FillWithScalar(PyArrayObject *arr, PyObject *obj)
 {
-    PyArray_Descr *dtype = NULL;
-    npy_longlong value_buffer[4];
-    char *value = NULL;
-    int retcode = 0;
+    npy_longlong value_buffer_stack[4];
+    char *value_buffer_heap = NULL;
+    char *buf = (char *)value_buffer_stack;
 
-    /*
-     * If 'arr' is an object array, copy the object as is unless
-     * 'obj' is a zero-dimensional array, in which case we copy
-     * the element in that array instead.
-     */
-    if (PyArray_DESCR(arr)->type_num == NPY_OBJECT &&
-                        !(PyArray_Check(obj) &&
-                          PyArray_NDIM((PyArrayObject *)obj) == 0)) {
-        value = (char *)&obj;
-
-        dtype = PyArray_DescrFromType(NPY_OBJECT);
-        if (dtype == NULL) {
-            return -1;
-        }
+    if ((size_t)PyArray_ITEMSIZE(arr) > sizeof(value_buffer_stack)) {
+        value_buffer_heap = PyMem_MALLOC(PyArray_ITEMSIZE(arr));
+        buf = value_buffer_heap;
     }
-    /* NumPy scalar */
-    else if (PyArray_IsScalar(obj, Generic)) {
-        dtype = PyArray_DescrFromScalar(obj);
-        if (dtype == NULL) {
-            return -1;
-        }
-        value = scalar_value(obj, dtype);
-        if (value == NULL) {
-            Py_DECREF(dtype);
-            return -1;
-        }
-    }
-    /* Python boolean */
-    else if (PyBool_Check(obj)) {
-        value = (char *)value_buffer;
-        *value = (obj == Py_True);
-
-        dtype = PyArray_DescrFromType(NPY_BOOL);
-        if (dtype == NULL) {
-            return -1;
-        }
+    if (PyArray_DESCR(arr)->flags & NPY_NEEDS_INIT) {
+        memset(buf, '\0', PyArray_DESCR(arr)->elsize);
     }
-    /* Python integer */
-    else if (PyLong_Check(obj)) {
-        /* Try long long before unsigned long long */
-        npy_longlong ll_v = PyLong_AsLongLong(obj);
-        if (error_converting(ll_v)) {
-            /* Long long failed, try unsigned long long */
-            npy_ulonglong ull_v;
-            PyErr_Clear();
-            ull_v = PyLong_AsUnsignedLongLong(obj);
-            if (ull_v == (unsigned long long)-1 && PyErr_Occurred()) {
-                return -1;
-            }
-            value = (char *)value_buffer;
-            *(npy_ulonglong *)value = ull_v;
-
-            dtype = PyArray_DescrFromType(NPY_ULONGLONG);
-            if (dtype == NULL) {
-                return -1;
-            }
-        }
-        else {
-            /* Long long succeeded */
-            value = (char *)value_buffer;
-            *(npy_longlong *)value = ll_v;
 
-            dtype = PyArray_DescrFromType(NPY_LONGLONG);
-            if (dtype == NULL) {
-                return -1;
-            }
-        }
-    }
-    /* Python float */
-    else if (PyFloat_Check(obj)) {
-        npy_double v = PyFloat_AsDouble(obj);
-        if (error_converting(v)) {
-            return -1;
-        }
-        value = (char *)value_buffer;
-        *(npy_double *)value = v;
-
-        dtype = PyArray_DescrFromType(NPY_DOUBLE);
-        if (dtype == NULL) {
-            return -1;
-        }
-    }
-    /* Python complex */
-    else if (PyComplex_Check(obj)) {
-        npy_double re, im;
-
-        re = PyComplex_RealAsDouble(obj);
-        if (error_converting(re)) {
-            return -1;
-        }
-        im = PyComplex_ImagAsDouble(obj);
-        if (error_converting(im)) {
-            return -1;
-        }
-        value = (char *)value_buffer;
-        ((npy_double *)value)[0] = re;
-        ((npy_double *)value)[1] = im;
-
-        dtype = PyArray_DescrFromType(NPY_CDOUBLE);
-        if (dtype == NULL) {
-            return -1;
-        }
-    }
-
-    /* Use the value pointer we got if possible */
-    if (value != NULL) {
-        /* TODO: switch to SAME_KIND casting */
-        retcode = PyArray_AssignRawScalar(arr, dtype, value,
-                                NULL, NPY_UNSAFE_CASTING);
-        Py_DECREF(dtype);
-        return retcode;
+    /* NOTE: PyArray_Pack effectively uses unsafe casting. */
+    if (PyArray_Pack(PyArray_DESCR(arr), buf, obj) < 0) {
+        PyMem_FREE(value_buffer_heap);
+        return -1;
     }
-    /* Otherwise convert to an array to do the assignment */
-    else {
-        PyArrayObject *src_arr;
-
-        /**
-         * The dtype of the destination is used when converting
-         * from the pyobject, so that for example a tuple gets
-         * recognized as a struct scalar of the required type.
-         */
-        Py_INCREF(PyArray_DTYPE(arr));
-        src_arr = (PyArrayObject *)PyArray_FromAny(obj,
-                        PyArray_DTYPE(arr), 0, 0, 0, NULL);
-        if (src_arr == NULL) {
-            return -1;
-        }
-
-        if (PyArray_NDIM(src_arr) != 0) {
-            PyErr_SetString(PyExc_ValueError,
-                    "Input object to FillWithScalar is not a scalar");
-            Py_DECREF(src_arr);
-            return -1;
-        }
-
-        retcode = PyArray_CopyInto(arr, src_arr);
 
-        Py_DECREF(src_arr);
-        return retcode;
-    }
+    int res = PyArray_AssignRawScalar(
+            arr, PyArray_DESCR(arr), buf, NULL, NPY_UNSAFE_CASTING);
+    PyMem_FREE(value_buffer_heap);
+    return res;
 }
 
 /*
diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py
index 5f0a725d2..fc2358f33 100644
--- a/numpy/core/tests/test_multiarray.py
+++ b/numpy/core/tests/test_multiarray.py
@@ -67,7 +67,7 @@ def _aligned_zeros(shape, dtype=float, order="C", align=None):
     # data pointer --- so we use and allocate size+1
     buf = buf[offset:offset+size+1][:-1]
     data = np.ndarray(shape, dtype, buf, order=order)
-    data.fill(0)
+    data.fill(np.uint8(0))
     return data

[proyan]

Thanks for the fix @seberg. I agree PyArray_Pack should solve the issue. I could validate it on our package once there is a PR.

[seberg]

@proyan the main branch now uses the normal scalar item setting (see gh-20924). So I expect this should be fixed, but it is not part of the NumPy 1.23 release.

Closing this, but please comment/reopen/open a new issue if there is still something to be done!

[proyan]

ref @jcarpent