Python: converting iterated assignment with atomic assignment using numpy does not work when matrix height> 256

I am working using numpy 1.6.2 and python 2.7 .
Given a matrix N x M x D

A

and a matrix I

that contains a list of indices.
I have to fill the matrix of zeros with the ACopy

sum of the element A

according to the definitions found in I

(see code).

Here is my code:

ACopy = zeros(A.shape)
for j in xrange(0, size(A, 0)):
  i = I[j]
  ACopy[j, i, :] = A[j, i, :] + A[j, i + 1, :]

      

        
        
        
      

    

Index Matrix:

I = array([2, 0, 3, 2, 1])

      

        
        
        
      

    

A

:

A = array([[[ 0,  1,  2],
        [ 3,  4,  5],
        [ 6,  7,  8],
        [ 9, 10, 11],
        [12, 13, 14]],

       [[15, 16, 17],
        [18, 19, 20],
        [21, 22, 23],
        [24, 25, 26],
        [27, 28, 29]],

       [[30, 31, 32],
        [33, 34, 35],
        [36, 37, 38],
        [39, 40, 41],
        [42, 43, 44]],

       [[45, 46, 47],
        [48, 49, 50],
        [51, 52, 53],
        [54, 55, 56],
        [57, 58, 59]],

       [[60, 61, 62],
        [63, 64, 65],
        [66, 67, 68],
        [69, 70, 71],
        [72, 73, 74]]])

      

        
        
        
      

    

I'm trying to improve my code by avoiding the for loop in this way:

r = r_[0:len(I)]
ACopy[r, I, :] = A[r, I, :] + A[r, I + 1, :]

      

        
        
        
      

    

I noticed that the output matrices are ACopy

different and I cannot figure out why. Any idea?

Thanks everyone!

EDIT . I am calculating a lot of matrices and I am trying to use np.array_equals(ACopy1,ACopy2)

where ACopy1

is the output of the first method and ACopy2

is the output of the second method. Sometimes the matrices are the same, but not every time. The output of the two methods should be the same or is there any burgundy case?

EDIT2 . I noticed that this weird behavior only happens when the height is over 256 . Here's my test suite:

from numpy import *

w = 5
h = 257

for i in xrange(1000):
  Z = random.rand(w, h, 5)

  I = (random.rand(w) * h - 1).astype(uint8)
  r = r_[0:w]
  ZCopy = zeros(Z.shape)
  ZCopy2 = zeros(Z.shape)

  for j in xrange(0, size(Z, 0)):
    i = I[j]
    ZCopy[j, i, :] = Z[j, i, :] + Z[j, i + 1, :]

  ZCopy2[r, I, :] = Z[r, I, :] + Z[r, I + 1, :]

  if (ZCopy - ZCopy2).any() > 0:
    print(ZCopy, ZCopy2, I)
    raise ValueError