Use C division to prevent requiring the GIL during loop #5

The following code excerpt from a cython extension is an example where the division remains as python division rather than C division. The issue here is that to do so it must reacquire the GIL, this happens inside of a heavily nested loop and is therefore not performant for threaded use.

src

...
 # release the GIL
    with nogil:

        # Now run the proper convolution
        for i in range(wkx, nx - wkx):
            for j in range(wky, ny - wky):
                top = 0.
                bot = 0.
                for ii in range(i - wkx, i + wkx + 1):
                    for jj in range(j - wky, j + wky + 1):
                        val = f[ii, jj]
                        ker = g[<unsigned int>(nkx - 1 - (wkx + ii - i)),
                                <unsigned int>(nky - 1 - (wky + jj - j))]
                        if not npy_isnan(val):
                            top += val * ker
                            bot += ker
                if normalize_by_kernel:
                    if bot == 0:
                        conv[i, j] = f[i, j]
                    else:
                        conv[i, j] = top / bot    #<--------------------------------------------------
                else:
                    conv[i, j] = top
    # GIL acquired again here
    return conv

This can be seen using cython's inspection feature, i.e. cython -a boundary_non.pyx

After specifying the use of C division with the decorator @cython.cdivision(True) the GIL is no longer involved.

astropy#3949 is the original issue related to releasing the GIL for this code. It relates to its performant usage with DASK, c.f. previous authors benchmark.

Test

• astropy 6117643

Utilizing the Jupyter notebook test code from the following benchmark.

import numpy as np
import dask.array as da
from astropy.convolution import convolve, Gaussian2DKernel

Nz,Ny,Nx = (64,1024,2048)
shape = (Nz,Ny,Nx)
big_numpy_array = np.random.rand(*shape)

cs = (1,Ny,Nx)
big_dask_array = da.from_array(big_numpy_array, chunks=cs)

def filter_func_ap(b):
    ker = Gaussian2DKernel(1)#, x_size=111, y_size=111)
    return convolve(b.squeeze(), ker, boundary=None)[np.newaxis,:,:]

%time big_dask_array.map_blocks(filter_func_ap).sum().compute()

%timeit filter_func_ap(big_numpy_array[0])

• With python division (@cython.cdivision(False)) the 1st test ran with an average wallclock time of 2.6s (over 10 runs). The second test, an average of 216ms/loop.
• With C division (@cython.cdivision(True)) the 1st test ran with an average wallclock time of 2.47s (over 10 runs). The second test, an average of 206ms/loop.

The use of C division and thus not having to re-acquire the GIL for this op shaves ~10ms off per loop. This is a 4.9% performance increase (per loop). The performance increase reported in the benchmark to release the GIL, showed a performance increase per loop of 6.6%. However, I believe this to be bogus. When testing the with gil: version myself, the time per loop is ~206ms. Effectively there isn't a performance increase per loop between tests of NOT releasing the GIl + python division AND when releasing the GIL + C division. This makes perfect sense as the former test where the GIL is not released does not, therefore, have to reacquire it within the nested loop for the division op. It does, however, have to be reacquired for the division op when releasing the GIL and using python division.

To Do

Investigate any numerical differences present in the results between Python and C division.
(None found).