BUG: wrong output of inplace operators x += y, x *= y, ... for nonstandard strides

The following script works for numpy 1.4.1 but not for numpy > 1.6.

See also the discussion on
http://thread.gmane.org/gmane.comp.python.numeric.general/51379

import numpy
print('numpy version =', numpy.__version__)

x = numpy.array([1,2])
y = numpy.array([1,2,3,4]).reshape((2,2))

x_expected = x.copy()
x_expected[0] += y[0,0]
x_expected[1] += y[0,1]
x_expected[0] += y[1,0]
x_expected[1] += y[1,1]

u,v = numpy.broadcast_arrays(x, y)
# u[0,0] += v[0,0]    <===> x[0] += y[0,0]
# u[0,1] += v[0,1]    <===> x[1] += y[0,1]
# u[1,0] += v[1,0]    <===> x[0] += y[1,0]
# u[1,1] += v[1,1]    <===> x[1] += y[1,1]
u += v

print('x=\n', x)
print('x_expected\n', x_expected)

numpy.testing.assert_array_almost_equal(x_expected, x)

NumPy 1.4.1 output

$ python inplace_op_bug.py 
numpy version = 1.4.1
x=
[5 8]
x_expected
[5 8]

NumPy 1.8.0.dev-9597b1f output

$ python inplace_op_bug.py 
numpy version = 1.8.0.dev-9597b1f
x=
[4 6]
x_expected
[5 8]
Traceback (most recent call last):
  File "inplace_op_bug.py", line 24, in <module>
    numpy.testing.assert_array_almost_equal(x_expected, x)
  File "/Users/b45ch1/.virtualenvs/sci/lib/python2.7/site-packages/numpy-1.8.0.dev_9597b1f_20120920-py2.7-macosx-10.8-x86_64.egg/numpy/testing/utils.py", line 811, in assert_array_almost_equal
    header=('Arrays are not almost equal to %d decimals' % decimal))
  File "/Users/b45ch1/.virtualenvs/sci/lib/python2.7/site-packages/numpy-1.8.0.dev_9597b1f_20120920-py2.7-macosx-10.8-x86_64.egg/numpy/testing/utils.py", line 644, in assert_array_compare
    raise AssertionError(msg)
AssertionError: 
Arrays are not almost equal to 6 decimals

(mismatch 100.0%)
 x: array([5, 8])
 y: array([4, 6])

Edit: updated for python3 printing

@seberg

Legitimate issue in 1.9-devel. I'm thinking iterator problem, or iterator use problem. @seberg Thoughts?

This is an (ab)use of in-place array operations, and probably issues like this have been discussed a lot since 1.6... NumPy optimizes the iteration order and does buffering, we are not going to change that.

Basically this is the same as a.T += a giving corrupted results (ignoring the error). We could (and should!) add a big fat warning about it (and I am sure there are more issues open about that, in some cases we could even force C-order to do what is "expected"), but unless the in-place operation is known to be unbuffered and 1-d (and even there it is undocumented, so...), this is not well defined.

That's what I thought a first, but the rhs is independent of the left so I couldn't see how that would be a problem in this case. Fancy indexing, yes, here, I'm not so sure.

Or maybe broadcast isn't setting things up correctly, or maybe it can't. Hmm...

It's interesting that [4,6] is exactly y.sum(axis=0).

No, broadcasting is good, it is the iterator. It is more complicating, you need buffering to understand it (honestly I am not sure why this is buffered) and buffer reuse...

@rkern

@rkern, actually, just coincidence here. It really throws away all but the last row.

It is a buffering issue: you are updating the same location multiple times, but because access is buffered only the last update shows up, right? And unless you go the convoluted way of setting up fake shapes with broadcast_arrays it raises an error, so you really need to force it to be bitten by this one:

>>> x = np.array([1, 2])
>>> y = np.arange(1, 5).reshape(2, 2)
>>> x += y
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: non-broadcastable output operand with shape (2,) doesn't match the broadcast shape (2,2)
>>> x + y
array([[2, 4],
       [4, 6]])

@jaimefrio

@jaimefrio, yes the buffer gets reused since I do not think the iterator has a notion for "there may be a reason that I can't reuse a buffer". It may be possible to inhibit buffer reuse or something like that if we also add warnings for overlapping memory.

You could also add a "just use the out array as initialized" to the reduction array, or add a specific path in the ufunc machinery for this: If "first input is output" -> make it readwrite and call reduction iteration logic (which only expects two operands).

Personally I think the general warnings are more important before trying to make this work again, and if you start trying to consider general memory overlaps, things may get ugly anyway...

We now have enough logic to detect the memory overlap, and buffering works as it should, the issue here seems to be something else: what do we expect when "broadcasting" in this fashion? IMO in the face of ambiguity we should raise an error.

In more detail: Note u.strides is (0, 8) even though its shape is (2, 2), so the data in u is actually multiple copies of x.

u, v = np.broadcast_arrays(x, y)
print(u.strides, u.shape)

The add ufunc uses writeback mechanisms to create an output array for u += v, but the writeback resolution cannot properly write 4 data elements to 2 slots. We should detect the improper output shape/strides and raise an error much like .x += y will do, it outputs

ValueError: non-broadcastable output operand with shape (2,) doesn't match the broadcast shape (2,2)

@mattip

@mattip do you suggest we always check if an output array has self overlap and warn if it does? Luckily I guess in most cases this is trivial, but it seems like it could be expensive to solve in odd cases.

Frankly, as cool as this trickery can be (I used to use it in the old days too), I do not think it is normally necessary, since that is why we have things like reduce and accumulate for.

So, personally, I do not really care anymore. Nowadays (ab)using self overlap by writing to it will never really do what you expect, functions to create self overlapping arrays warn either strongly in docs about it, or set the array to not being writeable. I think we got rid of this kind (ab)use for good. If someone really wants to do it, they probably have to go to C, luckily most of the abuses are probably not necessary (there used to be some 1-D array abuses that might even still be explicitly supported, at least in indexing assignment they exist).

I guess it is too bad that even now the result can depend on memory order probably, but meh :). I am against promoting self overlap tricks, and this stuff has been unreliable for so many years, that I hope it is a phenomenon that mostly grew out.

do you suggest we always check if an output array has self overlap and warn if it does? Luckily I guess in most cases this is trivial, but it seems like it could be expensive to solve in odd cases.

Not warn, I think we should error as the results are indeterminate. It is hard to imagine a case where this is desired behavior, if desired the user should be doing u += v[0] or so.

Sorry, right error makes sense. Now overlap check might be expensive in weird cases, but I guess not in anything normal, so adding an error sounds like the right idea.

I still believe that the correct return value is x_expected = [5 8].

Strides are really a very nice and elegant way to work with multi-dimensional arrays and it is a pity that they don't work as expected in numpy >= 1.6

From my understanding, the issue is that

import numpy
x = numpy.array([1,2])
y = numpy.array([1,2,3,4]).reshape((2,2))
u, v = numpy.broadcast_arrays(x, y)

print(x.shape, x.strides, x.data)
print(y.shape, y.strides, y.data)
print(v.shape, v.strides, v.data)

returns

((2,), (8,), <read-write buffer for 0x7fab783d2990, size 16, offset 0 at 0x7fab7835d8f0>)
((2, 2), (16, 8), <read-write buffer for 0x7fab782eb580, size 32, offset 0 at 0x7fab7835d8b0>)
((2, 2), (16, 8), <read-write buffer for 0x7fab782eb7b0, size 32, offset 0 at 0x7fab7835d9b0>)

but

In [21]: u.shape
Out[21]: (2, 2)

In [22]: u.strides
Out[22]: (0, 8)

In [23]: u.data
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-23-706b8128eb61> in <module>()
----> 1 u.data

AttributeError: cannot get single-segment buffer for discontiguous array

I would have expected that u.data is a contiguous memory block

<read-write buffer for 0x7fab783d2990, size 16, offset 0 at 0x7fab7835d8f0>

Is it possible in numpy to create a numpy array from an existing buffer?
I tried

 u = numpy.ndarray(shape=(2,2), strides=(0,8), buffer=x.data, dtype=numpy.int64)

but get the same error

In [53]: u
Out[53]: 
array([[1, 2],
       [1, 2]])

In [54]: u.data
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-54-706b8128eb61> in <module>()
----> 1 u.data

AttributeError: cannot get single-segment buffer for discontiguous array

What is the reason that a discontiguous array gets created?

That depends on your definition of contiguous. The numpy definition is that you can iterate through all elements by just doing one itemsize step. That is clearly violated, because you have to do a zero step.

The point is that what you feel is obvious, is dark magic to the majority of users. We want to protect users who get puzzled by u += u.T and not users who figure out how the (hidden!) C-backend inside numpy must be doing its iteration to achieve what they want. Also the specific reason why this went away: NumPy tries to optimize for speed in normal cases and not for ease of understanding the hidden C-level to use it in a way not quite intended.

All of your points are valid, as I said, I have used these things myself. But, we cannot optimize numpy for such marginal needs of those who are able to wrap their head around this and prefer that to the other points. The second point is that I believe this is a "big enough hammer" thing. If you have a thing that you really must do using this "obvious trick", we should have a more obvious way to do it, because stride tricks like this are not and should not obvious to the majority of users (feel free to debunk these claims).

Now, if someone suggests a solution that does limit the majority of use cases, I am not opposed. But I truly believe that the usage of such tricks can either already be avoided, or we should rather talk about other solutions for the problem. Just because the hammer looks so nice to a technically very adept part of numpy users, is IMO not enough reason to go to any lengths to support it.

About the last question, sorry. Yes of course you can get it. For example one way is to use arr.base, at some point you must reach a contiguous buffer. The code above seems to work for me in numpy 14.1, but it is likely a difference between python 2.7 and python 3 or so, because it might be that in some versions numpy cannot return a memoryview or so (which supports the non-contiguous buffer). I am surprised you (probably?) cannot put the array directly into np.ndarray as the buffer object, but it is probably the exact same reason.

I would have expected that u.data is a contiguous memory block

u holds basically two copies of the two data elements in x, this is what stride[0]==0 means

u	0	1
0	`x[0]`	`x[1]`
1	`x[0]`	`x[1]`

This overlap is non-contiguous.

Now if I add to this a true 2x2 ndarray with 4 values, what should the values of u[0,0] and u[0,1] be? They point to the same memory location and cannot be different, so a choice must be made. The change you noticed between the two NumPy versions is a change in the preference, but I think the correct answer (PR #11667) is not to allow ambiguous behaviour, as much as possible.

I am going to close this issue. There is no TODO here that I can see happening. NumPy now has overlap detection, and if anyone wants in-place magic aside from accumulate, the solution would be to petition for a new feature.

But overall, I think the actual solution is to write a cython, or other function working with the buffers directly (and doing whatever sanity checks make sense). It is simply too corner case and requires too much low-level understanding to expose nicely to Python IMO.

When this issue was new, writing something new in cython (or similar) was less easy, and these tricks were far more common. I remember being happy about hacking them myself. But the times have changed and it seems for the better :).

Navigation Menu

Repository navigation

BUG: wrong output of inplace operators x += y, x *= y, ... for nonstandard strides #2705

NumPy 1.4.1 output

NumPy 1.8.0.dev-9597b1f output

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

BUG: wrong output of inplace operators x += y, x *= y, ... for nonstandard strides #2705

Description

NumPy 1.4.1 output

NumPy 1.8.0.dev-9597b1f output

Activity

charris commented on Feb 21, 2014

seberg commented on Feb 21, 2014

charris commented on Feb 21, 2014

charris commented on Feb 21, 2014

rkern commented on Feb 21, 2014

seberg commented on Feb 21, 2014

seberg commented on Feb 21, 2014

jaimefrio commented on Feb 21, 2014

seberg commented on Feb 21, 2014

mattip commented on Aug 2, 2018

seberg commented on Aug 2, 2018

mattip commented on Aug 2, 2018

seberg commented on Aug 2, 2018

b45ch1 commented on Aug 3, 2018

seberg commented on Aug 3, 2018

seberg commented on Aug 3, 2018

mattip commented on Aug 3, 2018

seberg commented on May 20, 2021

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions