## Complete example

In this example, I will show how the addition function and the filter function implemented previously using the other approaches can be implemented using ctypes. First, the C-code which implements the algorithms contains the functions zadd, dadd, sadd, cadd, and dfilter2d. The zadd function is

/* Add arrays of contiguous data */

typedef struct {double real; double imag;} cdouble;

typedef struct {float real; float imag;} cfloat;

void zadd(cdouble *a, cdouble *b, cdouble *c, long n) {

c->real = a->real + b->real; c->imag = a->imag + b->imag; a++; b++; c++;

with similar code for cadd, dadd, and sadd that handles complex float, double, and float data-types, respectively:

void cadd(cfloat *a, cfloat *b, cfloat *c, long n) {

c->real = a->real + b->real; c->imag = a->imag + b->imag;

void dadd(double *a, double *b, double *c, long n) {

The code.c file also contains the function dfilter2d:

/* Assumes b is contiguous and a has strides that are multiples of sizeof(double)

void dfilter2d(double *a, double *b, int *astrides, int *dims) {

int i, j, M, N, S0, S1; int r, c, rm1, rp1, cp1, cm1;

M = dims[0]; N = dims[1]; S0 = astrides[0]/sizeof(double); S1=astrides[1]/sizeof(double); for (i = 1; i<M-1; i + +) {

r = i*S0; rp1 = r+S0; rm1 = r-S0; for (j=1; j <N-1; j++) {

c = j*S1; cp1 = j+S1; cm1 = j-S1; b[i*N+j] = a[r+c] +

(a [rp1 + c] + a[rm1 + c] + a[r+cp1] + a[r+cm1])*0.5 + (a[rp1+cp1] + a[rp1+cm1] +

A possible advantage this code has over the Fortran-equivalent code is that it takes arbitrarily strided (i.e. non-contiguous arrays) and may also run faster depending on the optimization capability of your compiler. But, it is a obviously more complicated than the simple code in filter.f. This code must be compiled into a shared library. On my Linux system this is accomplished using gcc -o code.so -shared code.c

Which creates a shared-library named code.so in the current directory. On Windows don't forget to either add __declspec(dllexport) in front of void on the line preceeding each function definition, or write a code.def file that lists the names of the functions to be exported.

A suitable Python interface to this shared library should be constructed. To do this create a file named interface.py with the following lines at the top:

import numpy as N import os

.path = os .path.dirname ('_file__' ) lib = N.ctypeslib. load-library('code', _path) _typedict = {'zadd' : complex, 'sadd' : N.single, 'cadd' : N.csingle, 'dadd' : float} for name in _typedict.keys () : val = getattr(lib, name) val.restype = None -type = -typedict[name]

val.argtypes = [N.ctypeslib.ndpointer(-type, flags='aligned, contiguous'), N.ctypeslib.ndpointer(_type, flags='aligned, contiguous'), N.ctypeslib.ndpointer(-type, flags='aligned, contiguous,'\

This code loads the shared library named code.<ext> located in the same path as this file. It then adds a return type of void to the functions contained in the library. It also adds argument checking to the functions in the library so that ndarrays can be passed as the first three arguments along with an integer (large enough to hold a pointer on the platform) as the fourth argument.

Setting up the filtering function is similar and allows the filtering function to be called with ndarray arguments as the first two arguments and with pointers to integers (large enough to handle the strides and shape of an ndarray) as the last two arguments.

lib.dfilter2d.restype=None lib.dfilter2d.argtypes = [N.ctypeslib.ndpointer(float, ndim=2, flags='aligned'), N.ctypeslib.ndpointer(float, ndim=2, flags='aligned, contiguous,'\ 'writeable'), ctypes.POINTER(N.ctypeslib.c_intp) , ctypes.POINTER(N.ctypeslib.c_intp) ]

Next, define a simple selection function that chooses which addition function to call in the shared library based on the data-type:

def select(dtype):

return lib.sadd, single elif dtype.char in ['F']:

return lib.cadd, csingle elif dtype.char in ['DG']:

return lib.dadd, float return func, ntype

Finally, the two functions to be exported by the interface can be written simply as def add(a, b):

requires = ['CONTIGUOUS', 'ALIGNED'] a = N.asanyarray(a) func, dtype = select(a.dtype) a = N.require(a, dtype, requires) b = N.require(b, dtype, requires) c = N.empty.like(a) func(a,b,c,a.size) return c and def filter2d(a):

a = N.require(a, float, ['ALIGNED']) b = N.zeros-like(a)

lib.dfilter2d(a, b, a.ctypes.strides, a.ctypes.shape) return b

0 0