Fast 2 dimensional array of floats for python (access/write)
Question: For my project use, I need to store certain amount (~100x100) of floats in two
dimensional array. And during the function calculation I need to read and
write to the array and since the function is really the bottleneck (consuming
98% of time) I really would need it to be fast.

I did some experiments with numpy and cython:

    
    
    import numpy
    import time
    cimport numpy
    cimport cython
    
    cdef int col, row
    
    DTYPE = numpy.int
    ctypedef numpy.int_t DTYPE_t
    cdef numpy.ndarray[DTYPE_t, ndim=2] matrix_c = numpy.zeros([100 + 1, 100 + 1], dtype=DTYPE)
    
    time_ = time.time()
    for l in xrange(5000):
        for col in xrange(100):
            for row in xrange(100):
                matrix_c[<unsigned int>row + 1][<unsigned int>col + 1] = matrix_c[<unsigned int>row][<unsigned int>col]
    print "Numpy + cython time: {0}".format(time.time() - time_)
    

but I found out that in spite of all my attempts, the version using python
lists, is still significantly faster.

Code using lists:

    
    
    matrix = []
    for i in xrange(100 + 1):
        matrix.append([])
        for j in xrange(100 + 1):
            matrix[i].append(0)
    
    time_ = time.time()
    for l in xrange(5000):
        for col in xrange(100):
            for row in xrange(100):
                matrix[row + 1][col + 1] = matrix[row][col]
    print "list time: {0}".format(time.time() - time_)
    

And results:

    
    
    list time: 0.0141758918762
    Numpy + cython time: 0.484772920609
    

Have I done something wrong? If not, is there anything that would help me to
improve the results?


Answer: Here's my version of the code you have. There are three functions, dealing
with integer arrays, 32 bit floating point arrays and double precision
floating point arrays, respectively.

    
    
    from numpy cimport ndarray as ar
    cimport numpy as np
    import numpy as np
    cimport cython
    import time
    
    @cython.boundscheck(False)
    @cython.wraparound(False)
    def access_write_int(ar[int,ndim=2] c, int n):
        cdef int l, col, row, h=c.shape[0], w=c.shape[1]
        time_ = time.time()
        for l in range(n):
            for row in range(h-1):
                for col in range(w-1):
                    c[row+1,col+1] = c[row,col]
        print "Numpy + cython time: {0}".format(time.time() - time_)
    
    @cython.boundscheck(False)
    @cython.wraparound(False)
    def access_write_float(ar[np.float32_t,ndim=2] c, int n):
        cdef int l, col, row, h=c.shape[0], w=c.shape[1]
        time_ = time.time()
        for l in range(n):
            for row in range(h-1):
                for col in range(w-1):
                    c[row+1,col+1] = c[row,col]
        print "Numpy + cython time: {0}".format(time.time() - time_)
    
    @cython.boundscheck(False)
    @cython.wraparound(False)
    def access_write_double(ar[double,ndim=2] c, int n):
        cdef int l, col, row, h=c.shape[0], w=c.shape[1]
        time_ = time.time()
        for l in range(n):
            for row in range(h-1):
                for col in range(w-1):
                    c[row+1,col+1] = c[row,col]
        print "Numpy + cython time: {0}".format(time.time() - time_)
    

To call these functions from Python I run this

    
    
    import numpy as np
    from numpy.random import rand, randint
    
    print "integers"
    c = randint(0, high=20, size=(101,101))
    access_write_int(c, 5000)
    print "32 bit float"
    c = rand(101, 101).astype(np.float32)
    access_write_float(c, 5000)
    print "double precision"
    c = rand(101, 101)
    access_write_double(c, 5000)
    

The following changes are important:

  1. Avoid slicing the array by accessing it using indices of the form `[i,j]` instead of `[i][j]`

  2. Define the variables `l`, `col`, and `row`, as integers so that the for loops run in C.

  3. Use the function decorators `@cython.boundscheck(False)` and '@cython.wraparound(False)` to turn off boundschecking and wraparound indexing for the key portion of the program. This allows for out of bounds memory accesses, so you should do this _only_ when you are certain that your indices are what they should be.

  4. Swap the two innermost `for` loops so that you access your array according to how it is arranged in memory. This makes a much bigger difference for larger arrays. The arrays given by `np.zeros` `np.random.rand`, etc. are usually C contiguous, so rows are stored in contiguous blocks and it is faster to vary the index along the rows in the outer `for` loop and not the inner one. If you want to keep the for loops as they are, consider taking the transpose of your array before you run the function on it so that the columns are in contiguous blocks instead.

