how is y += x computed when y.strides = (0, 8) and x.strides=(16, 8) ?
Hi, I'm using numpy 1.6.1 on Ubuntu 12.04.1 LTS. A code that used to work with an older version of numpy now fails with an error. Were there any changes in the way inplace operations like +=, *=, etc. work on arrays with nonstandard strides? For the script:  start of code  import numpy x = numpy.arange(6).reshape((3,2)) y = numpy.arange(2) print 'x=\n', x print 'y=\n', y u,v = numpy.broadcast_arrays(x, y) print 'u=\n', u print 'v=\n', v print 'v.strides=\n', v.strides v += u print 'v=\n', v # expectation: v = [[6,12], [6,12], [6,12]] print 'u=\n', u print 'y=\n', y # expectation: y = [6,12]  end of code  I get the output  start of output  x= [[0 1] [2 3] [4 5]] y= [0 1] u= [[0 1] [2 3] [4 5]] v= [[0 1] [0 1] [0 1]] v.strides= (0, 8) v= [[4 6] [4 6] [4 6]] u= [[0 1] [2 3] [4 5]] y= [4 6]  end of output  I would have expected that v += u performs an elementbyelement += v[0,0] += u[0,0] # increments y[0] v[0,1] += u[0,1] # increments y[1] v[1,0] += u[1,0] # increments y[0] v[1,1] += u[1,1] # increments y[1] v[2,0] += u[2,0] # increments y[0] v[2,1] += u[2,1] # increments y[1] yielding the result y = [6,12] but instead one obtains y = [4, 6] which could be the result of v[2,0] += u[2,0] # increments y[0] v[2,1] += u[2,1] # increments y[1] Is this the intended behavior? regards, Sebastian
Poor Sebastian, you make the mistake of asking difficult questions. I noticed that it should be [6, 10] not [6, 12], and in fact is with numpy1.4.1; while I observe the [4, 6] result with numpy1.6.1. Logs follow: numpy1.4.1 in Python2.6.5 on Mac (intel 64bit) with Python + numpy built from sources dualarch. The result in terms of output does not depend on the architecture chosen for run. The other is numpy1.6.1 with Python2.7.2. numpy1.4.1 (64bit Python 2.6.5): Python 2.6.5 (r265:79063, Jul 18 2010, 12:14:53) [GCC 4.2.1 (Apple Inc. build 5659)] on darwin Type "help", "copyright", "credits" or "license" for more information.
import numpy print numpy.__version__ 1.4.1 import numpy
x = numpy.arange(6).reshape((3,2)) y = numpy.arange(2)
print 'x=\n', x x= [[0 1] [2 3] [4 5]] print 'y=\n', y y= [0 1]
u,v = numpy.broadcast_arrays(x, y)
print 'u=\n', u u= [[0 1] [2 3] [4 5]] print 'v=\n', v v= [[0 1] [0 1] [0 1]] print 'v.strides=\n', v.strides v.strides= (0, 8)
v += u
print 'v=\n', v # expectation: v = [[6,12], [6,12], [6,12]] v= [[ 6 10] [ 6 10] [ 6 10]] print 'u=\n', u u= [[0 1] [2 3] [4 5]] print 'y=\n', y # expectation: y = [6,12] y= [ 6 10]
And numpy1.6.1 (64bit Python2.7.2): Python 2.7.2 (default, Mar 15 2012, 15:42:23) [GCC 4.2.1 (Apple Inc. build 5664)] on darwin Type "help", "copyright", "credits" or "license" for more information. [fiddling with sys.maxint edited out]
import numpy
x = numpy.arange(6).reshape((3,2)) y = numpy.arange(2)
print 'x=\n', x x= [[0 1] [2 3] [4 5]] print 'y=\n', y y= [0 1]
u,v = numpy.broadcast_arrays(x, y)
print 'u=\n', u u= [[0 1] [2 3] [4 5]] print 'v=\n', v v= [[0 1] [0 1] [0 1]] print 'v.strides=\n', v.strides v.strides= (0, 8)
v += u
print 'v=\n', v # expectation: v = [[6,12], [6,12], [6,12]] v= [[4 6] [4 6] [4 6]] print 'u=\n', u u= [[0 1] [2 3] [4 5]] print 'y=\n', y # expectation: y = [6,12] y= [4 6]
Maybe this helps bisecting it. Friedrich. P.S.: I took the time scrolling through the tickets, with an empty set resulting (first three pages by date or so). This does not mean such a ticket does not exist. Also the docs are rather quiet about this (e.g. http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.__iadd__.h..., or http://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html#arithmeticand...). My guess is that v is duplicated during the inplace addition (including strides), and the results of the calculation (using the original v) are written to that copy. After that, the copy's data would be [4, 6] with the strides as for v. So copying over to v's "original" data or letting it point to the copy and freeing the remaining copy would explain the behaviour after all. Maybe there's a sensible explanation for why it should be done that way. It is and remains a guess after all. I find the problem difficult to term after all. Maybe "Inplace operation not inplace for zerostrided first operand?" I don't think it is easy to get a good notion of an "inplace operation" for a zerostrided first operand. It might be that your definition does just not match up. But I think inplace operations should be expected to be independent on the order of execution of the elementwise operations, even if these elements share data (as in v in this case, the first operand). This criterion is fulfilled by your expectation and numpy1.4.1 but not by numpy1.6.1. I noticed that it's not necessary to duplicate v's strides in the hypothesis noted above. The neglection of the other elementoperations would then happen when copying the elementwise results over to v's storage. A remedy might be to let the backcopy respect the operation it is used for. So if the backcopying happens because an addition is to be inplace, it could use inplace addition instead of assigment. As a result, the elementwise operations on the "repeat()'ed" copy would be "chained together" by the operation under use. The noncommutative operations like subtraction or division would use addition or mulitiplication, resp., because the subtraction or division character can be thought of as an inversion of the second operand only before the whole inplace operations takes place, leading to a commutative inplace operation. Hence the "chained" operands in the backcopy operation form an associative and commutative expression after all. Nevertheless, a duplication of the first operands would occur in our case, since the same data from the first operand appears in several element's expressions in the copy buffer. The result would then be, in our example here, neither [6, 10] nor [4, 6], but [6, 12] instead. What makes actually more sense if you observe the symmetry regarding the quotient of the second column and the first column (the second column seen as a series, as it is to be collapes onto the same datum). It seems it comes down to giving the assigment a sensible meaning in case the first operand (i.e., the target) is partially zerostrided or contains less data elements that its view does. The second case might appear for stride tricks. It seems, it comes down to deciding in that case how to collapse the value of the second operand (the source of the assigment) with the values of the first operand (the target). And it seems from this consideration here, that the substitution from "a += b" to "a = a + b" is not generally valid if the substitute needs reshaping of the first operand "a" to match the shape of "b". This repetition is what makes the substitute generate a different result than the breakdown to elementwise operation in the substituted form. It seems, that inplace operations and parallelisation by array notions do not commute. To illustrate this noncummativity of inplace operations and parallelism, note that the implementation of "a += b" which yielded the 'expected result' (as in numpy1.4.1) goes first by elements, and then into the inplace operation (the iteration happens on a higher code level than the operation). In the other case, which is done by numpy1.6.1 in an equivalent way at least, the operations is executed first on array operands, and then an iterations is employed. Unfortunately, for the second approach, due to the nature of computers, an "execution of an operation on arrays" is abstract, and needs concrete implementation on the element level, s.t. the "a = a + b" substitute leads to an elementwise assigment, which is not accounted for by the algorithmic idea, and leads to the divergence between the two approaches' results. I like the equivalence criterion of "a += b" <=> "a = a + b", not caring for dtype downcasts. These downcasts are actually a Python idiosyncrasy and leads to the known paradox only when neglecting the character of Python assignment. The important notion of the substitutional equivalence of "a += b" <=> "a = a + b" seems to be to me, that the r.h.s. of the equivalence refers to additions, and explains the inplace operation by addition. Unfortunately, in the case of a zerostrided first operand, the assignment in between becomes illdefined. But if then we then stand by "explaining inplace addition by pure addition" it would seem natural to me, to chain the operands resulting from the substitute "a + b" together by just that "+" operation. I don't know what symbol to put for that, maybe a "=+", how strange this might ever look. So the usual Python teaching of replacing "+=" by "=" and "+" seems to need a more precise definition of that "=" used in case of arrays as operands. For scalars, of course an ordinary assigment does very well. In numpy, I would hence state as a formula containing the preceding sentences: "a += b <=> a =+ a + b". It would be possible to implement this by filling the first operand by the unity of the respective operation (0 for addition and 1 for multiplication), carrying out the elementwise operation on broadcasted arrays, and feeding it back to the original first (left) operand by chaining using the resp. operation, which could in turn be done by scalar inplace operation using the unity data as a start. Funny is, e.g. the following hypothetical Python session (I turn the arrows upwards for this): ^^^ a = numpy.ones(()) ^^^ b = numpy.zeros((10, 1)) ^^^ a array(1.0) ^^^ b array([[ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.]]) ^^^ u, v = numpy.broadcast_arrays(a, b) ^^^ u array([[ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.]]) ^^^ v array([[ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.]]) ^^^ u += v ^^^ u array([[10.], [10.], [10.], [10.], [10.], [10.], [10.], [10.], [10.], [10.]]) ^^^ v array([[ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.]]) ^^^ a array(10.0) ^^^ b array([[ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.]]) ^^^ Funnily too, the operations "a += b" (in my case), or "y += x" would be perfectly defined even if "a + b" or "y + x" (latter in Sebastian's case) require broadcasting of the first operand. Currently they yield the following result: (numpy1.5.1, as above:)
a += b Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: invalid return array shape
(numpy1.6.1, as above:)
y += x Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: nonbroadcastable output operand with shape (2) doesn't match the broadcast shape (3,2)
The result for Sebstian would then be [6, 12] for the y variable as noted above already. Actually it's interesting that any arbitrarily added additional parallelisation by forcing the first operand to be broadcasted yields changes in it even if the second operand is the zero element w.r.t. the operation considered. But I think that's rather by design then. Noone was asking for invariance under arbitrary repetition, no? I don't know how this would work with multithreading. Probably not very well as it would require locking on the thenshared target operand on the l.h.s. of, e.g., "a += b". It's not easy to tell apart the cases where all elements of "a" are independent on each other. But actually it's possible at all only in restriction of the dtype to an atomic data type. What's always the case if object arrays are considered as pointer arrays. But I really know too little on this multithreading optimisation as it's really not my point of interest and rather an additional idea which probably can be made consistent with the idea proposed here if some work is done on it. This was the postscriptum, apparently. F. Am 31.08.2012 um 11:31 schrieb Sebastian Walter:
Hi,
I'm using numpy 1.6.1 on Ubuntu 12.04.1 LTS. A code that used to work with an older version of numpy now fails with an error.
Were there any changes in the way inplace operations like +=, *=, etc. work on arrays with nonstandard strides?
For the script:
 start of code 
import numpy
x = numpy.arange(6).reshape((3,2)) y = numpy.arange(2)
print 'x=\n', x print 'y=\n', y
u,v = numpy.broadcast_arrays(x, y)
print 'u=\n', u print 'v=\n', v print 'v.strides=\n', v.strides
v += u
print 'v=\n', v # expectation: v = [[6,12], [6,12], [6,12]] print 'u=\n', u print 'y=\n', y # expectation: y = [6,12]
 end of code 
I get the output
 start of output  x= [[0 1] [2 3] [4 5]] y= [0 1] u= [[0 1] [2 3] [4 5]] v= [[0 1] [0 1] [0 1]] v.strides= (0, 8) v= [[4 6] [4 6] [4 6]] u= [[0 1] [2 3] [4 5]] y= [4 6]
 end of output 
I would have expected that
v += u
performs an elementbyelement +=
v[0,0] += u[0,0] # increments y[0] v[0,1] += u[0,1] # increments y[1] v[1,0] += u[1,0] # increments y[0] v[1,1] += u[1,1] # increments y[1] v[2,0] += u[2,0] # increments y[0] v[2,1] += u[2,1] # increments y[1]
yielding the result
y = [6,12]
but instead one obtains
y = [4, 6]
which could be the result of
v[2,0] += u[2,0] # increments y[0] v[2,1] += u[2,1] # increments y[1]
Is this the intended behavior?
regards, Sebastian _______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
Hey, No idea if this is simply not support or just a bug, though I am guessing that such usage simply is not planned. However, this also has to do with buffering, so unless the behaviour is substantially changed, I would not expect even predictable results. I have used things like a[1:] += a[:1] (cumsum is better here, and clear as to what it does), but this is different, as here the same data is read after being written. And here an example showing the buffer is involved: In [14]: np.setbufsize(1024) Out[14]: 1024 In [19]: x = numpy.arange(6*1000).reshape(1,2) In [20]: y = numpy.arange(2) In [21]: u,v = numpy.broadcast_arrays(x, y) In [22]: v += u In [23]: v Out[23]: array([[21348, 21355], [21348, 21355], [21348, 21355], ..., [21348, 21355], [21348, 21355], [21348, 21355]]) In [24]: np.setbufsize(1000000) Out[24]: 1024 # note it gives old bufsize... In [25]: x = numpy.arange(6*1000).reshape(1,2) In [26]: y = numpy.arange(2) In [27]: u,v = numpy.broadcast_arrays(x, y) In [28]: v += u In [29]: v Out[29]: array([[5998, 6000], [5998, 6000], [5998, 6000], ..., [5998, 6000], [5998, 6000], [5998, 6000]]) On Thu, 20120906 at 01:55 +0200, Friedrich Romstedt wrote:
Poor Sebastian, you make the mistake of asking difficult questions.
I noticed that it should be [6, 10] not [6, 12], and in fact is with numpy1.4.1; while I observe the [4, 6] result with numpy1.6.1. Logs follow:
numpy1.4.1 in Python2.6.5 on Mac (intel 64bit) with Python + numpy built from sources dualarch. The result in terms of output does not depend on the architecture chosen for run. The other is numpy1.6.1 with Python2.7.2.
numpy1.4.1 (64bit Python 2.6.5):
Python 2.6.5 (r265:79063, Jul 18 2010, 12:14:53) [GCC 4.2.1 (Apple Inc. build 5659)] on darwin Type "help", "copyright", "credits" or "license" for more information.
import numpy print numpy.__version__ 1.4.1 import numpy
x = numpy.arange(6).reshape((3,2)) y = numpy.arange(2)
print 'x=\n', x x= [[0 1] [2 3] [4 5]] print 'y=\n', y y= [0 1]
u,v = numpy.broadcast_arrays(x, y)
print 'u=\n', u u= [[0 1] [2 3] [4 5]] print 'v=\n', v v= [[0 1] [0 1] [0 1]] print 'v.strides=\n', v.strides v.strides= (0, 8)
v += u
print 'v=\n', v # expectation: v = [[6,12], [6,12], [6,12]] v= [[ 6 10] [ 6 10] [ 6 10]] print 'u=\n', u u= [[0 1] [2 3] [4 5]] print 'y=\n', y # expectation: y = [6,12] y= [ 6 10]
And numpy1.6.1 (64bit Python2.7.2):
Python 2.7.2 (default, Mar 15 2012, 15:42:23) [GCC 4.2.1 (Apple Inc. build 5664)] on darwin Type "help", "copyright", "credits" or "license" for more information. [fiddling with sys.maxint edited out]
import numpy
x = numpy.arange(6).reshape((3,2)) y = numpy.arange(2)
print 'x=\n', x x= [[0 1] [2 3] [4 5]] print 'y=\n', y y= [0 1]
u,v = numpy.broadcast_arrays(x, y)
print 'u=\n', u u= [[0 1] [2 3] [4 5]] print 'v=\n', v v= [[0 1] [0 1] [0 1]] print 'v.strides=\n', v.strides v.strides= (0, 8)
v += u
print 'v=\n', v # expectation: v = [[6,12], [6,12], [6,12]] v= [[4 6] [4 6] [4 6]] print 'u=\n', u u= [[0 1] [2 3] [4 5]] print 'y=\n', y # expectation: y = [6,12] y= [4 6]
Maybe this helps bisecting it.
Friedrich.
P.S.: I took the time scrolling through the tickets, with an empty set resulting (first three pages by date or so). This does not mean such a ticket does not exist. Also the docs are rather quiet about this (e.g. http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.__iadd__.h..., or http://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html#arithmeticand...).
My guess is that v is duplicated during the inplace addition (including strides), and the results of the calculation (using the original v) are written to that copy. After that, the copy's data would be [4, 6] with the strides as for v. So copying over to v's "original" data or letting it point to the copy and freeing the remaining copy would explain the behaviour after all. Maybe there's a sensible explanation for why it should be done that way.
It is and remains a guess after all.
I find the problem difficult to term after all. Maybe "Inplace operation not inplace for zerostrided first operand?" I don't think it is easy to get a good notion of an "inplace operation" for a zerostrided first operand. It might be that your definition does just not match up. But I think inplace operations should be expected to be independent on the order of execution of the elementwise operations, even if these elements share data (as in v in this case, the first operand). This criterion is fulfilled by your expectation and numpy1.4.1 but not by numpy1.6.1.
I noticed that it's not necessary to duplicate v's strides in the hypothesis noted above. The neglection of the other elementoperations would then happen when copying the elementwise results over to v's storage.
A remedy might be to let the backcopy respect the operation it is used for. So if the backcopying happens because an addition is to be inplace, it could use inplace addition instead of assigment. As a result, the elementwise operations on the "repeat()'ed" copy would be "chained together" by the operation under use. The noncommutative operations like subtraction or division would use addition or mulitiplication, resp., because the subtraction or division character can be thought of as an inversion of the second operand only before the whole inplace operations takes place, leading to a commutative inplace operation. Hence the "chained" operands in the backcopy operation form an associative and commutative expression after all. Nevertheless, a duplication of the first operands would occur in our case, since the same data from the first operand appears in several element's expressions in the copy buffer.
The result would then be, in our example here, neither [6, 10] nor [4, 6], but [6, 12] instead. What makes actually more sense if you observe the symmetry regarding the quotient of the second column and the first column (the second column seen as a series, as it is to be collapes onto the same datum). It seems it comes down to giving the assigment a sensible meaning in case the first operand (i.e., the target) is partially zerostrided or contains less data elements that its view does. The second case might appear for stride tricks.
It seems, it comes down to deciding in that case how to collapse the value of the second operand (the source of the assigment) with the values of the first operand (the target). And it seems from this consideration here, that the substitution from "a += b" to "a = a + b" is not generally valid if the substitute needs reshaping of the first operand "a" to match the shape of "b". This repetition is what makes the substitute generate a different result than the breakdown to elementwise operation in the substituted form. It seems, that inplace operations and parallelisation by array notions do not commute.
To illustrate this noncummativity of inplace operations and parallelism, note that the implementation of "a += b" which yielded the 'expected result' (as in numpy1.4.1) goes first by elements, and then into the inplace operation (the iteration happens on a higher code level than the operation). In the other case, which is done by numpy1.6.1 in an equivalent way at least, the operations is executed first on array operands, and then an iterations is employed. Unfortunately, for the second approach, due to the nature of computers, an "execution of an operation on arrays" is abstract, and needs concrete implementation on the element level, s.t. the "a = a + b" substitute leads to an elementwise assigment, which is not accounted for by the algorithmic idea, and leads to the divergence between the two approaches' results.
I like the equivalence criterion of "a += b" <=> "a = a + b", not caring for dtype downcasts. These downcasts are actually a Python idiosyncrasy and leads to the known paradox only when neglecting the character of Python assignment. The important notion of the substitutional equivalence of "a += b" <=> "a = a + b" seems to be to me, that the r.h.s. of the equivalence refers to additions, and explains the inplace operation by addition. Unfortunately, in the case of a zerostrided first operand, the assignment in between becomes illdefined. But if then we then stand by "explaining inplace addition by pure addition" it would seem natural to me, to chain the operands resulting from the substitute "a + b" together by just that "+" operation. I don't know what symbol to put for that, maybe a "=+", how strange this might ever look. So the usual Python teaching of replacing "+=" by "=" and "+" seems to need a more precise definition of that "=" used in case of arrays as oper ands. For scalars, of course an ordinary assigment does very well. In numpy, I would hence state as a formula containing the preceding sentences: "a += b <=> a =+ a + b".
It would be possible to implement this by filling the first operand by the unity of the respective operation (0 for addition and 1 for multiplication), carrying out the elementwise operation on broadcasted arrays, and feeding it back to the original first (left) operand by chaining using the resp. operation, which could in turn be done by scalar inplace operation using the unity data as a start.
Funny is, e.g. the following hypothetical Python session (I turn the arrows upwards for this):
^^^ a = numpy.ones(()) ^^^ b = numpy.zeros((10, 1)) ^^^ a array(1.0) ^^^ b array([[ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.]]) ^^^ u, v = numpy.broadcast_arrays(a, b) ^^^ u array([[ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.], [ 1.]]) ^^^ v array([[ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.]]) ^^^ u += v ^^^ u array([[10.], [10.], [10.], [10.], [10.], [10.], [10.], [10.], [10.], [10.]]) ^^^ v array([[ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.]]) ^^^ a array(10.0) ^^^ b array([[ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.], [ 0.]]) ^^^
Funnily too, the operations "a += b" (in my case), or "y += x" would be perfectly defined even if "a + b" or "y + x" (latter in Sebastian's case) require broadcasting of the first operand. Currently they yield the following result:
(numpy1.5.1, as above:)
a += b Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: invalid return array shape
(numpy1.6.1, as above:)
y += x Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: nonbroadcastable output operand with shape (2) doesn't match the broadcast shape (3,2)
The result for Sebstian would then be [6, 12] for the y variable as noted above already.
Actually it's interesting that any arbitrarily added additional parallelisation by forcing the first operand to be broadcasted yields changes in it even if the second operand is the zero element w.r.t. the operation considered. But I think that's rather by design then. Noone was asking for invariance under arbitrary repetition, no?
I don't know how this would work with multithreading. Probably not very well as it would require locking on the thenshared target operand on the l.h.s. of, e.g., "a += b". It's not easy to tell apart the cases where all elements of "a" are independent on each other. But actually it's possible at all only in restriction of the dtype to an atomic data type. What's always the case if object arrays are considered as pointer arrays. But I really know too little on this multithreading optimisation as it's really not my point of interest and rather an additional idea which probably can be made consistent with the idea proposed here if some work is done on it.
This was the postscriptum, apparently. F.
Am 31.08.2012 um 11:31 schrieb Sebastian Walter:
Hi,
I'm using numpy 1.6.1 on Ubuntu 12.04.1 LTS. A code that used to work with an older version of numpy now fails with an error.
Were there any changes in the way inplace operations like +=, *=, etc. work on arrays with nonstandard strides?
For the script:
 start of code 
import numpy
x = numpy.arange(6).reshape((3,2)) y = numpy.arange(2)
print 'x=\n', x print 'y=\n', y
u,v = numpy.broadcast_arrays(x, y)
print 'u=\n', u print 'v=\n', v print 'v.strides=\n', v.strides
v += u
print 'v=\n', v # expectation: v = [[6,12], [6,12], [6,12]] print 'u=\n', u print 'y=\n', y # expectation: y = [6,12]
 end of code 
I get the output
 start of output  x= [[0 1] [2 3] [4 5]] y= [0 1] u= [[0 1] [2 3] [4 5]] v= [[0 1] [0 1] [0 1]] v.strides= (0, 8) v= [[4 6] [4 6] [4 6]] u= [[0 1] [2 3] [4 5]] y= [4 6]
 end of output 
I would have expected that
v += u
performs an elementbyelement +=
v[0,0] += u[0,0] # increments y[0] v[0,1] += u[0,1] # increments y[1] v[1,0] += u[1,0] # increments y[0] v[1,1] += u[1,1] # increments y[1] v[2,0] += u[2,0] # increments y[0] v[2,1] += u[2,1] # increments y[1]
yielding the result
y = [6,12]
but instead one obtains
y = [4, 6]
which could be the result of
v[2,0] += u[2,0] # increments y[0] v[2,1] += u[2,1] # increments y[1]
Is this the intended behavior?
regards, Sebastian _______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
_______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
On Thu, Sep 6, 2012 at 1:41 AM, Sebastian Berg <sebastian@sipsolutions.net> wrote:
Hey,
No idea if this is simply not support or just a bug, though I am guessing that such usage simply is not planned.
I think that's right... currently numpy simply makes no guarantees about what order ufunc loops will be performed in, or even if they will be performed in any strictly sequential order. In ordinary cases this lets it make various optimizations, but it means that you can't count on any specific behaviour for the unusual case where different locations in the output array are stored in overlapping memory. Fixing this would require two things: (a) Some code to detect when an array may have internal overlaps (sort of like np.may_share_memory for axes). Not entirely trivial. (b) A "fallback mode" for ufuncs where if the code in (a) detects that we are (probably) dealing with one of these arrays, it processes the operations in some predictable order without buffering. I suppose if someone wanted to come up with these two pieces, and it didn't look like it would cause slowdowns in common cases, the code in (b) avoided creating duplicate code paths that increased maintenance burden, etc., then probably noone would object to making these arrays act in a better defined way? I don't think most people are that worried about this though. Your original code would be much clearer if it just used np.sum... n
I'd like to have a look at the implementation of iadd in numpy, but I'm having a real hard time to find the corresponding code. I'm basically stuck at https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/number.... Could someone give me a pointer where to find it? Respectively, could someone point me to some documentation where the (folder/file) structure of the numpy sources is explained? Sebastian On Thu, Sep 6, 2012 at 2:58 PM, Nathaniel Smith <njs@pobox.com> wrote:
On Thu, Sep 6, 2012 at 1:41 AM, Sebastian Berg <sebastian@sipsolutions.net> wrote:
Hey,
No idea if this is simply not support or just a bug, though I am guessing that such usage simply is not planned.
I think that's right... currently numpy simply makes no guarantees about what order ufunc loops will be performed in, or even if they will be performed in any strictly sequential order. In ordinary cases this lets it make various optimizations, but it means that you can't count on any specific behaviour for the unusual case where different locations in the output array are stored in overlapping memory.
Fixing this would require two things: (a) Some code to detect when an array may have internal overlaps (sort of like np.may_share_memory for axes). Not entirely trivial. (b) A "fallback mode" for ufuncs where if the code in (a) detects that we are (probably) dealing with one of these arrays, it processes the operations in some predictable order without buffering.
I suppose if someone wanted to come up with these two pieces, and it didn't look like it would cause slowdowns in common cases, the code in (b) avoided creating duplicate code paths that increased maintenance burden, etc., then probably noone would object to making these arrays act in a better defined way? I don't think most people are that worried about this though. Your original code would be much clearer if it just used np.sum...
n _______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
On Wed, Oct 17, 2012 at 11:38 AM, Sebastian Walter <sebastian.walter@gmail.com> wrote:
I'd like to have a look at the implementation of iadd in numpy, but I'm having a real hard time to find the corresponding code.
I'm basically stuck at https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/number....
n_ops is essentially a map of function pointers set in the umath module (see PyArray_SetNumericOps, used in umath module). IOW, n_ops.add is a pointer to umath.add, itself set up through a generated file __umath_generated.c: f = PyUFunc_FromFuncAndData(add_functions, ...); PyDict_SetItemString(dictionary, "add", f); At that point, you need to look into the ufunc machinery: for double all along, the add type resolver should be pretty simple and in the end call DOUBLE_add.
Could someone give me a pointer where to find it? Respectively, could someone point me to some documentation where the (folder/file) structure of the numpy sources is explained?
There is sadly not much explanation on the code structure for the C part. src/multiarray contains the code for the multiarray extension (array, dtype, broadcasting, iteration) and src/umath contains the code for the umath extension (ufunc machinery + core loops implementation). David
participants (5)

David Cournapeau

Friedrich Romstedt

Nathaniel Smith

Sebastian Berg

Sebastian Walter