Adding types¶

To add types we use custom Cython syntax, so we are now breaking Python sourcecompatibility. Consider this code (read the comments!) :

# tag: numpy# You can ignore the previous line.# It's for internal testing of the cython documentation.importnumpyasnp# "cimport" is used to import special compile-time information# about the numpy module (this is stored in a file numpy.pxd which is# distributed with Numpy).# Here we've used the name "cnp" to make it easier to understand what# comes from the cimported module and what comes from the imported module,# however you can use the same name for both if you wish.cimportnumpyascnp# It's necessary to call "import_array" if you use any part of the# numpy PyArray_* API. From Cython 3, accessing attributes like# ".shape" on a typed Numpy array use this API. Therefore we recommend# always calling "import_array" whenever you "cimport numpy"cnp.import_array()# We now need to fix a datatype for our arrays. I've used the variable# DTYPE for this, which is assigned to the usual NumPy runtime# type info object.DTYPE=np.int64# "ctypedef" assigns a corresponding compile-time type to DTYPE_t. For# every type in the numpy module there's a corresponding compile-time# type with a _t-suffix.ctypedefcnp.int64_tDTYPE_t# "def" can type its arguments but not have a return type. The type of the# arguments for a "def" function is checked at run-time when entering the# function.## The arrays f, g and h is typed as "np.ndarray" instances. The only effect# this has is to a) insert checks that the function arguments really are# NumPy arrays, and b) make some attribute access like f.shape[0] much# more efficient. (In this example this doesn't matter though.)defnaive_convolve(cnp.ndarrayf,cnp.ndarrayg):ifg.shape[0]%2!=1org.shape[1]%2!=1:raiseValueError("Only odd dimensions on filter supported")assertf.dtype==DTYPEandg.dtype==DTYPE# The "cdef" keyword is also used within functions to type variables. It# can only be used at the top indentation level (there are non-trivial# problems with allowing them in other places, though we'd love to see# good and thought out proposals for it).## For the indices, the "int" type is used. This corresponds to a C int,# other C types (like "unsigned int") could have been used instead.# Purists could use "Py_ssize_t" which is the proper Python type for# array indices.cdefintvmax=f.shape[0]cdefintwmax=f.shape[1]cdefintsmax=g.shape[0]cdefinttmax=g.shape[1]cdefintsmid=smax//2cdefinttmid=tmax//2cdefintxmax=vmax+2*smidcdefintymax=wmax+2*tmidcdefcnp.ndarrayh=np.zeros([xmax,ymax],dtype=DTYPE)cdefintx,y,s,t,v,w# It is very important to type ALL your variables. You do not get any# warnings if not, only much slower code (they are implicitly typed as# Python objects).cdefints_from,s_to,t_from,t_to# For the value variable, we want to use the same data type as is# stored in the array, so we use "DTYPE_t" as defined above.# NB! An important side-effect of this is that if "value" overflows its# datatype size, it will simply wrap around like in C, rather than raise# an error like in Python.cdefDTYPE_tvalueforxinrange(xmax):foryinrange(ymax):s_from=max(smid-x,-smid)s_to=min((xmax-x)-smid,smid+1)t_from=max(tmid-y,-tmid)t_to=min((ymax-y)-tmid,tmid+1)value=0forsinrange(s_from,s_to):fortinrange(t_from,t_to):v=x-smid+sw=y-tmid+tvalue+=g[smid-s,tmid-t]*f[v,w]h[x,y]=valuereturnh

After building this and continuing my (very informal) benchmarks, I get:

In [21]:importconvolve2In [22]:%timeit -n2 -r3 convolve2.naive_convolve(f, g)2 loops, best of 3: 828 ms per loop

Efficient indexing¶

There’s still a bottleneck killing performance, and that is the array lookupsand assignments. The[]-operator still uses full Python operations –what we would like to do instead is to access the data buffer directly at Cspeed.

What we need to do then is to type the contents of thendarray objects.We do this with a special “buffer” syntax which must be told the datatype(first argument) and number of dimensions (“ndim” keyword-only argument, ifnot provided then one-dimensional is assumed).

These are the needed changes:

...defnaive_convolve(cnp.ndarray[DTYPE_t,ndim=2]f,cnp.ndarray[DTYPE_t,ndim=2]g):...cdefcnp.ndarray[DTYPE_t,ndim=2]h=...

Usage:

In [18]:importconvolve3In [19]:%timeit -n3 -r100 convolve3.naive_convolve(f, g)3 loops, best of 100: 11.6 ms per loop

Note the importance of this change.

Gotcha: This efficient indexing only affects certain index operations,namely those with exactlyndim number of typed integer indices. So ifv for instance isn’t typed, then the lookupf[v,w] isn’toptimized. On the other hand this means that you can continue using Pythonobjects for sophisticated dynamic slicing etc. just as when the array is nottyped.

Tuning indexing further¶

The array lookups are still slowed down by two factors:

1. Bounds checking is performed.

2. Negative indices are checked for and handled correctly. The code above isexplicitly coded so that it doesn’t use negative indices, and it(hopefully) always access within bounds. We can add a decorator to disablebounds checking:

   ...cimportcython@cython.boundscheck(False)# turn off bounds-checking for entire function@cython.wraparound(False)# turn off negative index wrapping for entire functiondefnaive_convolve(cnp.ndarray[DTYPE_t,ndim=2]f,cnp.ndarray[DTYPE_t,ndim=2]g):...

Now bounds checking is not performed (and, as a side-effect, if you ‘’do’’happen to access out of bounds you will in the best case crash your programand in the worst case corrupt data). It is possible to switch bounds-checkingmode in many ways, seeCompiler directives for moreinformation.

Also, we’ve disabled the check to wrap negative indices (e.g. g[-1] givingthe last value). As with disabling bounds checking, bad things will happenif we try to actually use negative indices with this disabled.

The function call overhead now starts to play a role, so we compare the lattertwo examples with larger N:

In [11]:%timeit -n3 -r100 convolve4.naive_convolve(f, g)3 loops, best of 100: 5.97 ms per loopIn [12]:N=1000In [13]:f=np.arange(N*N,dtype=np.int64).reshape((N,N))In [14]:g=np.arange(81,dtype=np.int64).reshape((9,9))In [17]:%timeit -n1 -r10 convolve3.naive_convolve(f, g)1 loops, best of 10: 1.16 s per loopIn [18]:%timeit -n1 -r10 convolve4.naive_convolve(f, g)1 loops, best of 10: 597 ms per loop

(Also this is a mixed benchmark as the result array is allocated within thefunction call.)

What typing does not do¶

The main purpose of typing things asndarray is to allow efficientindexing of single elements, and to speed up access to a small number ofattributes such as.shape. Typing does not allow Cython to speedup mathematical operations on the whole array (for example, adding two arraystogether). Typing does not allow Cython to speed up calls to Numpy globalfunctions or to methods of the array.

More generic code¶

It would be possible to do:

defnaive_convolve(object[DTYPE_t,ndim=2]f,...):

i.e. useobject rather thancnp.ndarray. Under Python 3.0 thiscan allow your algorithm to work with any libraries supporting the bufferinterface; and support for e.g. the Python Imaging Library may easily be addedif someone is interested also under Python 2.x.

There is some speed penalty to this though (as one makes more assumptionscompile-time if the type is set tocnp.ndarray, specifically it isassumed that the data is stored in pure strided mode and not in indirectmode).

Buffer options¶

The following options are accepted when creating buffer types:

• ndim - an integer number of dimensions.

• mode - a string from:

  • "c" - C contiguous array,

  • "fortran" - Fortran contiguous array,

  • "strided" - non-contiguous lookup into a single block of memory,

  • "full" - any valid buffer, including indirect arrays.

• negative_indices - boolean value specifying whether negative indexing is allowed, essentiallya per-variable version of the compiler directivecython.wraparound.

• cast - boolean value specifying whether to allow the user to view the array as a differenttype. The sizes of the source and destination type must be the same. In C++ this would beequivalent toreinterpret_cast.

In all cases these parameters must be compile-time constants.

As an example of how to specify the parameters:

cdefcnp.ndarray[double,ndim=2,mode="c",cast=True]some_array

cast can be used to get a low-level view of an array with non-native endianness:

cdefcnp.ndarray[cnp.uint32,cast=True]values=np.arange(10,dtype='>i4')

although correctly interpreting the cast data is the user’s responsibility.