8000 BUG: Incorrect result from `add.at()` · Issue #23457 · numpy/numpy · GitHub
[go: up one dir, main page]

Skip to content

BUG: Incorrect result from add.at() #23457

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
WarrenWeckesser opened this issue Mar 26, 2023 · 1 comment · Fixed by #23467
Closed

BUG: Incorrect result from add.at() #23457

WarrenWeckesser opened this issue Mar 26, 2023 · 1 comment · Fixed by #23467

Comments

@WarrenWeckesser
Copy link
Member
WarrenWeckesser commented Mar 26, 2023

Describe the issue:

I'm using a 64 bit Linux computer with the main development branch of NumPy. When the dtype of the indices parameter of np.add.at is not a 64 bit integer, and the length of indices is greater than 8192, and the third parameter b is an array with the same length as indices, I get an incorrect result from np.add.at(a, indices, b). In the example below, instead of adding the value from b at index 8192, it appears to wrap around and add the value from index 0 again. So instead of getting the expected value 115, a[0] is 25.

Git bisect points the finger at #23136.

Reproduce the code example:

import numpy as np


# The result is incorrect for n > 8192
n = 8193

# The result is incorrect for dtypes int16 and int32; int64 is OK.
indices = np.zeros(n, dtype=np.int16)

b = np.zeros(n)
b[0] = 10
b[1] = 5
b[8192:] = 100

a = np.zeros(1)
np.add.at(a, indices, b)

print(f'{a[0]     = }')
print(f'expected = {b.sum()}')

Error message:

a[0]     = 25.0
expected = 115.0

Runtime information:

In [33]: import sys, numpy; print(numpy.__version__); print(sys.version)
1.25.0.dev0+1008.ga37978a106
3.10.8 (main, Oct 24 2022, 08:32:26) [GCC 11.2.0]

In [34]: numpy.show_runtime()
[{'numpy_version': '1.25.0.dev0+1008.ga37978a106',
  'python': '3.10.8 (main, Oct 24 2022, 08:32:26) [GCC 11.2.0]',
  'uname': uname_result(system='Linux', node='pop-os', release='6.2.0-76060200-generic', version='#202302191831~1678319661~22.04~4d98339 SMP PREEMPT_DYNAMIC Thu M', machine='x86_64')},
 {'simd_extensions': {'baseline': ['SSE', 'SSE2', 'SSE3'],
                      'found': ['SSSE3',
                                'SSE41',
                                'POPCNT',
                                'SSE42',
                                'AVX',
                                'F16C',
                                'FMA3',
                                'AVX2'],
                      'not_found': ['AVX512F',
                                    'AVX512CD',
                                    'AVX512_KNL',
                                    'AVX512_KNM',
                                    'AVX512_SKX',
                                    'AVX512_CLX',
                                    'AVX512_CNL',
                                    'AVX512_ICL']}},
 {'architecture': 'Zen',
  'filepath': '/usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so',
  'internal_api': 'openblas',
  'num_threads': 24,
  'prefix': 'libopenblas',
  'threading_layer': 'pthreads',
  'user_api': 'blas',
  'version': '0.3.20'}]

Context for the issue:

I ran into this while experimenting with a method to use numpy.add.at to multiply a sparse COO matrix by a regular 1d NumPy array. It worked for small arrays, but then broke for larger ones.

@seberg
Copy link
Member
seberg commented Mar 26, 2023

Ah, fix should be:

diff --git a/numpy/core/src/umath/ufunc_object.c b/numpy/core/src/umath/ufunc_object.c
index a5e8f4cbe4..97a5e4125d 100644
--- a/numpy/core/src/umath/ufunc_object.c
+++ b/numpy/core/src/umath/ufunc_object.c
@@ -5935,6 +5935,9 @@ trivial_at_loop(PyArrayMethodObject *ufuncimpl, NPY_ARRAYMETHOD_FLAGS flags,
 
         res = ufuncimpl->contiguous_indexed_loop(
                 context, args, inner_size, steps, NULL);
+        if (args[2] != NULL) {
+            args[2] += *inner_size;
+        }
     } while (res == 0 && iter->outer_next(iter->outer));
 
     if (res == 0 && !(flags & NPY_METH_NO_FLOATINGPOINT_ERRORS)) {

Don't want to writing test tonight though (or checking). If the integer is intp compatible, the loop will only have one iteration (normally).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants
0