No speed gains from Cython again?

cython, performance, python

Solution

The primary issue you are seeing is that Cython only optimizes single-item access for typed arrays. This means that each of the lines in your code where you are using vectorization from NumPy still involve creating and interacting with Python objects. The code you have there wasn't faster than the pure Python version because it wasn't really doing any of the computation differently. You will have to avoid this by writing out the looping operations explicitly. Here is a modified version of your code that runs significantly faster.

from numpy cimport ndarray as ar
from numpy cimport int32_t as int32
from numpy import empty
from numpy.random import randint
cimport cython
ctypedef int

# Notice the use of these decorators to tell Cython to turn off
# some of the checking it does when accessing arrays.
@cython.boundscheck(False)
@cython.wraparound(False)
def boots(int32 trial, ar[double, ndim=2] empirical, ar[double, ndim=2] expected):
    cdef:
        int32 length = empirical.shape[0], i, j, k
        int32 o
        ar[double, ndim=2] ret = empty((trial, 100))
        ar[int32] choices
        ar[double] m = empty(100), n = empty(100)
    for i in range(trial):
        # Still calling Python on this line
        choices = randint(0, length, length)
        # It was faster to compute m and n separately.
        # I suspect that has to do with cache management.
        # Instead of allocating new arrays, I just filled the old ones with the new values.
        o = choices[0]
        for k in range(100):
            m[k] = empirical[o,k]
        for j in range(1, length):
            o = choices[j]
            for k in range(100):
                m[k] += empirical[o,k]
        o = choices[0]
        for k in range(100):
            n[k] = expected[o,k]
        for j in range(1, length):
            o = choices[j]
            for k in range(100):
                n[k] += expected[o,k]
        # Here I simplified some of the math and got rid of temporary arrays
        for k in range(100):
            ret[i,k] = m[k] / n[k] - 1.
    ret.sort(axis=0)
    return ret[int(trial * 0.025)].reshape((10,10)), ret[int(trial * 0.975)].reshape((10,10))

If you want to have a look at which lines of your code involve Python calls, the Cython compiler can generate an html file showing which lines call Python. This option is called annotation. The way you use it depends on how you are compiling your cython code. If you are using the IPython notebook, just add the `--annotate` flag to the Cython cell magic.

You may also be able to benefit from turning on the C compiler optimization flags.

Problem

The following is my cython code, the purpose is to do a bootstrap. ``` def boots(int trial, np.ndarray[double, ndim=2] empirical, np.ndarray[double, ndim=2] expected): cdef int length = len(empirical) cdef np.ndarray[double, ndim=2] ret = np.empty((trial, 100)) cdef np.ndarray[long] choices cdef np.ndarray[double] m cdef np.ndarray[double] n cdef long o cdef int i cdef int j for i in range(trial): choices = np.random.randint(0, length, length) m = np.zeros(100) n = np.zeros(100) for j in range(length): o = choices[j] m.__iadd__(empirical[o]) n.__iadd__(expected[o]) empirical_boot = m / length expected_boot = n / length ret[i] = empirical_boot / expected_boot - 1 ret.sort(axis=0) return ret[int(trial * 0.025)].reshape((10,10)), ret[int(trial * 0.975)].reshape((10,10)) # test code empirical = np.ones((40000, 100)) expected = np.ones((40000, 100)) %prun -l 10 boots(100, empirical,expected) ``` It takes 11 seconds in pure python with fancy indexing, and no matter how hard I tuned in cython it stays the same. `np.random.randint(0, 40000, 40000)` takes 1 ms, so 100x takes 0.1s. `np.sort(np.ones((40000, 100))` takes 0.2s. Thus I feel there must be ways to improve `boots`.

Original source