Compound assignment operators in Python's NumPy library The "vectorizing" of fancy indexing by Python's NumPy library sometimes gives unexpected results.

For example:

import numpy
a = numpy.zeros((1000,4), dtype='uint32')
b = numpy.zeros((1000,4), dtype='uint32')
i = numpy.random.random_integers(0,999,1000)
j = numpy.random.random_integers(0,3,1000)

a[i,j] += 1
for k in xrange(1000):
    b[i[k],j[k]] += 1


It gives different results in the arrays 'a' and 'b' (i.e. the appearance of tuple (i,j) appears as 1 in 'a' regardless of repeats, whereas repeats are counted in 'b').

This is easily verified as follows:

numpy.sum(a)
883
numpy.sum(b)
1000


It is also notable that the fancy indexing version is almost two orders of magnitude faster than the for loop.

Is there an efficient way for NumPy to compute the repeat counts as implemented using the for loop in the provided example?
Ben,

The following modification to your code is a more accurate reflection of my larger problem:



    from timeit import timeit
    import numpy

    # Dimensions
    ROWS = 10000000
    SAMPLES = 10000
    COLS = 4

    # Number of times to run timeit
    TEST_COUNT = 100 

    # Matrices
    a = numpy.zeros((ROWS, COLS), dtype='uint32')
    b = numpy.zeros((ROWS, COLS), dtype='uint32')
    c = numpy.zeros((ROWS, COLS), dtype='uint32')

    # Indices
    i = numpy.random.random_integers(0, ROWS - 1, SAMPLES)
    j = numpy.random.random_integers(0, COLS - 1, SAMPLES)

    # Fancy
    def fancy():
        a[i,j] += 1
    print "Fancy:", timeit("fancy()", setup="from __main__ import fancy", number=TEST_COUNT)

    # Loop
    def loop():
        for k in xrange(SAMPLES):
            b[i[k],j[k]] += 1
    print "Loop:", timeit("loop()", setup="from __main__ import loop", number=TEST_COUNT)

    # Flat
    def flat():
        flatIndices = i*COLS + j
        sums = numpy.bincount(flatIndices)
        c.flat[:len(sums)] += sums

    print "Flat:", timeit("flat()", setup="from __main__ import flat", number=TEST_COUNT)

    # Check that for_loop and custom are equivalent
    print "Equivalent:", (c == b).all()



With the exception that there are probably few samples that are actually repeated in this version. Let's ignore that difference for now. On my machine, your code gave the following timings:



    Fancy: 0.109203100204
    Loop: 7.37937998772
    Flat: 126.571173906
    Equivalent: True



This slowdown is largely attributable to the large sparse array output by bincount().