BUG: (1.25.0) ufunc.at returns wrong results when indices is an array containing negative ints

[rsokl]

Describe the issue:

According to docs, numpy.add.at(x, indices, v) should be equivalent to x[indices] += v (other than when indices contains repeat elements). As of numpy 1.25.0 this is no longer the case when indices is an array (or list) containing negative values.

For example, here is the issue demonstrated via numpy.add.at

>>> x1 = np.array([1.0])
>>> x2 = np.array([1.0])
>>> index = np.array([-1])
>>> v = 1.0

>>> np.add.at(x1, index, v)
>>> x2[index] += v
>>> x1, x2  # x1 and x2 should match, but do not for numpy 1.25.0
(array([1.]), array([2.]))

And via numpy.multiply.at

>>> x1 = np.array([1.0])
>>> x2 = np.array([1.0])
>>> index = np.array([-1])
>>> v = 0.0

>>> np.multiply.at(x1, index, v)
>>> x2[index] *= v
>>> x1, x2  # x1 and x2 should match, but do not for numpy 1.25.0
(array([1.]), array([0.]))

Reproduce the code example:

The following code demonstrates that the discrepancy occurs specifically when indices is an array (not a scalar) of negative ints.

import numpy as np
print(f"{np.__version__=}\n")

print("out1 = np.array([0.0])")
print("out2 = np.array([0.0])")
print()

for n in [0, -1]:
    for index in [n, np.array([n])]:        
        out1 = np.array([0.0])
        out2 = np.array([0.0])
        np.add.at(out1, index, 1.0)
        out2[index] += 1
        print(f"{index=!r}\n----------")
        print(f"np.add.at(out1, index, 1.0) -> {out1!r}")
        print(f"out2[index] += 1 -> {out2}")
        print()
    print()

Error message:

No response

Runtime information:

1.25.0
3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 17:59:51) [MSC v.1935 64 bit (AMD64)]
[{'numpy_version': '1.25.0',
  'python': '3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 17:59:51) '
            '[MSC v.1935 64 bit (AMD64)]',
  'uname': uname_result(system='Windows', node='Bohm', release='10', version='10.0.22621', machine='AMD64')},
 {'simd_extensions': {'baseline': ['SSE', 'SSE2', 'SSE3'],
                      'found': ['SSSE3',
                                'SSE41',
                                'POPCNT',
                                'SSE42',
                                'AVX',
                                'F16C',
                                'FMA3',
                                'AVX2',
                                'AVX512F',
                                'AVX512CD',
                                'AVX512_SKX',
                                'AVX512_CLX',
                                'AVX512_CNL',
                                'AVX512_ICL'],
                      'not_found': []}},
 {'filepath': 'C:\\Users\\rsokl\\miniconda3\\envs\\npt\\Library\\bin\\libblas.dll',
  'internal_api': 'mkl',
  'num_threads': 4,
  'prefix': 'libblas',
  'threading_layer': 'intel',
  'user_api': 'blas',
  'version': '2022.1-Product'}]
None

Context for the issue:

This is a major bug, as ufuncs are silently returning incorrect results for common usage patterns. E.g., this is causing mygrad's CI to fail (and ultimately segfault - but I am not 100% sure that this issue is the precise cause for the segfault).

P.S. I caught this bug using the property based testing library Hypothesis 😄

[mattip]

Hmm. It seems this is a result of the "fastest path" work I did in #23136. The "fast path" and "slow path" that use the iter API are not affected. Apparently we don't yet have a test with negative indices.

[mattip]

It is coming from this line

numpy/numpy/core/src/umath/ufunc_object.c

Line 5822 in c143f7b

args[1] = (char *) iter->outer_ptrs[0];

which directly accesses the iter->outer_ptr. Is there a missing call to something that will unwrap the indexes if they are passed in as a ndarray and not a list of ints?

[mattip]

Thanks @rsokl