Optimize python for Connected Component Labeling Area of Subsets
Question: I have a binary map on which I do Connected Component Labeling and get
something like this for a 64x64 grid - <http://pastebin.com/bauas0NJ>

Now I want to group them by label, so that I can find their area and their
center of mass. This is what I do:

    
    
    #ccl_np is the computed array from the previous step (see pastebin)
    #I discard the label '1' as its the background
    unique, count = np.unique(ccl_np, return_counts = True)
    xcm_array = []
    ycm_array = []
    
    for i in range(1,len(unique)):
        subarray = np.where(ccl_np == unique[i])
        xcm_array.append("{0:.5f}".format((sum(subarray[0]))/(count[i]*1.)))
        ycm_array.append("{0:.5f}".format((sum(subarray[1]))/(count[i]*1.)))
    
    final_array = zip(xcm_array,ycm_array,count[1:])
    

I want a fast code (as I will be doing this for grids of size 4096x4096) and
was told to check out numba. Here's my naive attempt :

    
    
    unique, inverse, count = np.unique(ccl_np, return_counts = True, return_inverse = True)
    xcm_array = np.zeros(len(count),dtype=np.float32)
    ycm_array = np.zeros(len(count),dtype=np.float32)
    inverse = inverse.reshape(64,64)
    
    @numba.autojit
    def mysolver(xcm_array, ycm_array, inverse, count):
        for i in range(64):
            for j in range(64):
                pos = inverse[i][j]
                local_count = count[pos]
                xcm_array[pos] += i/(local_count*1.)
                ycm_array[pos] += j/(local_count*1.)
    
    
    mysolver(xcm_array, ycm_array, inverse, count)
    
    final_array = zip(xcm_array,ycm_array,count)
    

To my surprise, using numba was slower or at best equal to the speed of the
previous way. What am I doing wrong ? Also, can this be done in Cython and
will that be faster ?

I am using the included packages in the latest Anaconda python 2.7
distribution.


Answer: I believe the issue might be that you are timing jit'd code incorrectly. The
first time you run the code, your timing includes the time it takes numba to
compile the code. This is called warming up the jit. If you call it again,
that cost is gone.

    
    
    import numpy as np
    import numba as nb
    
    unique, inverse, count = np.unique(ccl_np, return_counts = True, return_inverse = True)
    xcm_array = np.zeros(len(count),dtype=np.float32)
    ycm_array = np.zeros(len(count),dtype=np.float32)
    inverse = inverse.reshape(64,64)
    
    def mysolver(xcm_array, ycm_array, inverse, count):
        for i in range(64):
            for j in range(64):
                pos = inverse[i][j]
                local_count = count[pos]
                xcm_array[pos] += i/(local_count*1.)
                ycm_array[pos] += j/(local_count*1.)
    
    @nb.jit(nopython=True)
    def mysolver_nb(xcm_array, ycm_array, inverse, count):
        for i in range(64):
            for j in range(64):
                pos = inverse[i,j]
                local_count = count[pos]
                xcm_array[pos] += i/(local_count*1.)
                ycm_array[pos] += j/(local_count*1.)
    

Then the timings with `timeit` which runs the code multiple times. First the
plain python version:

    
    
    In [4]:%timeit mysolver(xcm_array, ycm_array, inverse, count)
    10 loops, best of 3: 25.8 ms per loop
    

and then with numba:

    
    
    In [5]: %timeit mysolver_nb(xcm_array, ycm_array, inverse, count)
    The slowest run took 3630.44 times longer than the fastest. This could mean         that an intermediate result is being cached 
    10000 loops, best of 3: 33.1 µs per loop
    

The numba code is ~1000 times faster.

