CoCalc Shared Filesscientific-python-lectures / Lecture-6A-Fortran-and-C.ipynb

# Using Fortran and C code with Python

J.R. Johansson (jrjohansson at gmail.com)

The latest version of this IPython notebook lecture is available at http://github.com/jrjohansson/scientific-python-lectures.

The other notebooks in this lecture series are indexed at http://jrjohansson.github.io.

In [1]:
%pylab inline
from IPython.display import Image

Populating the interactive namespace from numpy and matplotlib

The advantage of Python is that it is flexible and easy to program. The time it takes to setup a new calulation is therefore short. But for certain types of calculations Python (and any other interpreted language) can be very slow. It is particularly iterations over large arrays that is difficult to do efficiently.

Such calculations may be implemented in a compiled language such as C or Fortran. In Python it is relatively easy to call out to libraries with compiled C or Fortran code. In this lecture we will look at how to do that.

But before we go ahead and work on optimizing anything, it is always worthwhile to ask....

In [2]:
Image(filename='images/optimizing-what.png')


## Fortran

### F2PY

F2PY is a program that (almost) automatically wraps fortran code for use in Python: By using the f2py program we can compile fortran code into a module that we can import in a Python program.

F2PY is a part of NumPy, but you will also need to have a fortran compiler to run the examples below.

### Example 0: scalar input, no output

In [3]:
%%file hellofortran.f
C File  hellofortran.f
subroutine hellofortran (n)
integer n

do 100 i=0, n
print *, "Fortran says hello"
100     continue
end

Overwriting hellofortran.f

Generate a python module using f2py:

In [4]:
!f2py -c -m hellofortran hellofortran.f

running build running config_cc unifing config_cc, config, build_clib, build_ext, build commands --compiler options running config_fc unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options running build_src build_src building extension "hellofortran" sources f2py options: [] f2py:> /tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.c creating /tmp/tmpz2IPjB/src.linux-x86_64-2.7 Reading fortran codes... Reading file 'hellofortran.f' (format:fix,strict) Post-processing... Block: hellofortran Block: hellofortran Post-processing (stage 2)... Building modules... Building module "hellofortran"... Constructing wrapper function "hellofortran"... hellofortran(n) Wrote C/API module "hellofortran" to file "/tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.c" adding '/tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.c' to sources. adding '/tmp/tmpz2IPjB/src.linux-x86_64-2.7' to include_dirs. copying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpz2IPjB/src.linux-x86_64-2.7 copying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpz2IPjB/src.linux-x86_64-2.7 build_src: building npy-pkg config files running build_ext customize UnixCCompiler customize UnixCCompiler using build_ext customize Gnu95FCompiler Found executable /usr/bin/gfortran customize Gnu95FCompiler customize Gnu95FCompiler using build_ext building 'hellofortran' extension compiling C sources C compiler: x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC creating /tmp/tmpz2IPjB/tmp creating /tmp/tmpz2IPjB/tmp/tmpz2IPjB creating /tmp/tmpz2IPjB/tmp/tmpz2IPjB/src.linux-x86_64-2.7 compile options: '-I/tmp/tmpz2IPjB/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c' x86_64-linux-gnu-gcc: /tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.c In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4, from /tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.h:13, from /tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.c:17: /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp] #warning "Using deprecated NumPy API, disable it by " \ ^ x86_64-linux-gnu-gcc: /tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.c In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4, from /tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.h:13, from /tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.c:2: /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp] #warning "Using deprecated NumPy API, disable it by " \ ^ compiling Fortran sources Fortran f77 compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops Fortran f90 compiler: /usr/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops Fortran fix compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -Wall -fno-second-underscore -fPIC -O3 -funroll-loops compile options: '-I/tmp/tmpz2IPjB/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c' gfortran:f77: hellofortran.f /usr/bin/gfortran -Wall -Wall -shared /tmp/tmpz2IPjB/tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.o /tmp/tmpz2IPjB/tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.o /tmp/tmpz2IPjB/hellofortran.o -lgfortran -o ./hellofortran.so Removing build directory /tmp/tmpz2IPjB 

Example of a python script that use the module:

In [5]:
%%file hello.py
import hellofortran

hellofortran.hellofortran(5)

Overwriting hello.py
In [6]:
# run the script
!python hello.py

Fortran says hello Fortran says hello Fortran says hello Fortran says hello Fortran says hello Fortran says hello

### Example 1: vector input and scalar output

In [7]:
%%file dprod.f

subroutine dprod(x, y, n)

double precision x(n), y
y = 1.0

do 100 i=1, n
y = y * x(i)
100    continue
end

Overwriting dprod.f
In [8]:
!rm -f dprod.pyf
!f2py -m dprod -h dprod.pyf dprod.f

Reading fortran codes... Reading file 'dprod.f' (format:fix,strict) Post-processing... Block: dprod {} In: :dprod:dprod.f:dprod vars2fortran: No typespec for argument "n". Block: dprod Post-processing (stage 2)... Saving signatures to file "./dprod.pyf"

The f2py program generated a module declaration file called dsum.pyf. Let's look what's in it:

In [9]:
!cat dprod.pyf


The module does not know what Fortran subroutine arguments is input and output, so we need to manually edit the module declaration files and mark output variables with intent(out) and input variable with intent(in):

In [10]:
%%file dprod.pyf
python module dprod ! in
interface  ! in :dprod
subroutine dprod(x,y,n) ! in :dprod:dprod.f
double precision dimension(n), intent(in) :: x
double precision, intent(out) :: y
integer, optional,check(len(x)>=n),depend(x),intent(in) :: n=len(x)
end subroutine dprod
end interface
end python module dprod

Overwriting dprod.pyf

Compile the fortran code into a module that can be included in python:

In [11]:
!f2py -c dprod.pyf dprod.f

running build running config_cc unifing config_cc, config, build_clib, build_ext, build commands --compiler options running config_fc unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options running build_src build_src building extension "dprod" sources creating /tmp/tmpWyCvx1/src.linux-x86_64-2.7 f2py options: [] f2py: dprod.pyf Reading fortran codes... Reading file 'dprod.pyf' (format:free) Post-processing... Block: dprod Block: dprod Post-processing (stage 2)... Building modules... Building module "dprod"... Constructing wrapper function "dprod"... y = dprod(x,[n]) Wrote C/API module "dprod" to file "/tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.c" adding '/tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.c' to sources. adding '/tmp/tmpWyCvx1/src.linux-x86_64-2.7' to include_dirs. copying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpWyCvx1/src.linux-x86_64-2.7 copying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpWyCvx1/src.linux-x86_64-2.7 build_src: building npy-pkg config files running build_ext customize UnixCCompiler customize UnixCCompiler using build_ext customize Gnu95FCompiler Found executable /usr/bin/gfortran customize Gnu95FCompiler customize Gnu95FCompiler using build_ext building 'dprod' extension compiling C sources C compiler: x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC creating /tmp/tmpWyCvx1/tmp creating /tmp/tmpWyCvx1/tmp/tmpWyCvx1 creating /tmp/tmpWyCvx1/tmp/tmpWyCvx1/src.linux-x86_64-2.7 compile options: '-I/tmp/tmpWyCvx1/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c' x86_64-linux-gnu-gcc: /tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.c In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4, from /tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.h:13, from /tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.c:18: /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp] #warning "Using deprecated NumPy API, disable it by " \ ^ /tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.c:111:12: warning: ‘f2py_size’ defined but not used [-Wunused-function] static int f2py_size(PyArrayObject* var, ...) ^ x86_64-linux-gnu-gcc: /tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.c In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4, from /tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.h:13, from /tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.c:2: /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp] #warning "Using deprecated NumPy API, disable it by " \ ^ compiling Fortran sources Fortran f77 compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops Fortran f90 compiler: /usr/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops Fortran fix compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -Wall -fno-second-underscore -fPIC -O3 -funroll-loops compile options: '-I/tmp/tmpWyCvx1/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c' gfortran:f77: dprod.f /usr/bin/gfortran -Wall -Wall -shared /tmp/tmpWyCvx1/tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.o /tmp/tmpWyCvx1/tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.o /tmp/tmpWyCvx1/dprod.o -lgfortran -o ./dprod.so Removing build directory /tmp/tmpWyCvx1 

#### Using the module from Python

In [12]:
import dprod

In [13]:
help(dprod)

Help on module dprod: NAME dprod FILE /home/rob/Desktop/scientific-python-lectures/dprod.so DESCRIPTION This module 'dprod' is auto-generated with f2py (version:2). Functions: y = dprod(x,n=len(x)) . DATA __version__ = '$Revision:$' dprod = <fortran object> VERSION
In [14]:
dprod.dprod(arange(1,50))

6.082818640342675e+62
In [15]:
# compare to numpy
prod(arange(1.0,50.0))

6.0828186403426752e+62
In [16]:
dprod.dprod(arange(1,10), 5) # only the 5 first elements

120.0

Compare performance:

In [17]:
xvec = rand(500)

In [18]:
timeit dprod.dprod(xvec)

1000000 loops, best of 3: 882 ns per loop
In [19]:
timeit xvec.prod()

100000 loops, best of 3: 4.45 µs per loop

### Example 2: cummulative sum, vector input and vector output

The cummulative sum function for an array of data is a good example of a loop intense algorithm: Loop through a vector and store the cummulative sum in another vector.

In [20]:
# simple python algorithm: example of a SLOW implementation
# Why? Because the loop is implemented in python.
def py_dcumsum(a):
b = empty_like(a)
b[0] = a[0]
for n in range(1,len(a)):
b[n] = b[n-1]+a[n]
return b


Fortran subroutine for the same thing: here we have added the intent(in) and intent(out) as comment lines in the original fortran code, so we do not need to manually edit the fortran module declaration file generated by f2py.

In [21]:
%%file dcumsum.f
c File dcumsum.f
subroutine dcumsum(a, b, n)
double precision a(n)
double precision b(n)
integer n
cf2py  intent(in) :: a
cf2py  intent(out) :: b
cf2py  intent(hide) :: n

b(1) = a(1)
do 100 i=2, n
b(i) = b(i-1) + a(i)
100    continue
end

Overwriting dcumsum.f

We can directly compile the fortran code to a python module:

In [22]:
!f2py -c dcumsum.f -m dcumsum

running build running config_cc unifing config_cc, config, build_clib, build_ext, build commands --compiler options running config_fc unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options running build_src build_src building extension "dcumsum" sources f2py options: [] f2py:> /tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c creating /tmp/tmpfvrMl6/src.linux-x86_64-2.7 Reading fortran codes... Reading file 'dcumsum.f' (format:fix,strict) Post-processing... Block: dcumsum Block: dcumsum Post-processing (stage 2)... Building modules... Building module "dcumsum"... Constructing wrapper function "dcumsum"... b = dcumsum(a) Wrote C/API module "dcumsum" to file "/tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c" adding '/tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.c' to sources. adding '/tmp/tmpfvrMl6/src.linux-x86_64-2.7' to include_dirs. copying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpfvrMl6/src.linux-x86_64-2.7 copying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpfvrMl6/src.linux-x86_64-2.7 build_src: building npy-pkg config files running build_ext customize UnixCCompiler customize UnixCCompiler using build_ext customize Gnu95FCompiler Found executable /usr/bin/gfortran customize Gnu95FCompiler customize Gnu95FCompiler using build_ext building 'dcumsum' extension compiling C sources C compiler: x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC creating /tmp/tmpfvrMl6/tmp creating /tmp/tmpfvrMl6/tmp/tmpfvrMl6 creating /tmp/tmpfvrMl6/tmp/tmpfvrMl6/src.linux-x86_64-2.7 compile options: '-I/tmp/tmpfvrMl6/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c' x86_64-linux-gnu-gcc: /tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4, from /tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.h:13, from /tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c:18: /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp] #warning "Using deprecated NumPy API, disable it by " \ ^ /tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c:111:12: warning: ‘f2py_size’ defined but not used [-Wunused-function] static int f2py_size(PyArrayObject* var, ...) ^ x86_64-linux-gnu-gcc: /tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.c In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17, from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4, from /tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.h:13, from /tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.c:2: /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp] #warning "Using deprecated NumPy API, disable it by " \ ^ compiling Fortran sources Fortran f77 compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops Fortran f90 compiler: /usr/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops Fortran fix compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -Wall -fno-second-underscore -fPIC -O3 -funroll-loops compile options: '-I/tmp/tmpfvrMl6/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c' gfortran:f77: dcumsum.f /usr/bin/gfortran -Wall -Wall -shared /tmp/tmpfvrMl6/tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.o /tmp/tmpfvrMl6/tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.o /tmp/tmpfvrMl6/dcumsum.o -lgfortran -o ./dcumsum.so Removing build directory /tmp/tmpfvrMl6 
In [23]:
import dcumsum

In [24]:
a = array([1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0])

In [25]:
py_dcumsum(a)

array([ 1., 3., 6., 10., 15., 21., 28., 36.])
In [26]:
dcumsum.dcumsum(a)

array([ 1., 3., 6., 10., 15., 21., 28., 36.])
In [27]:
cumsum(a)

array([ 1., 3., 6., 10., 15., 21., 28., 36.])

Benchmark the different implementations:

In [28]:
a = rand(10000)

In [29]:
timeit py_dcumsum(a)

100 loops, best of 3: 4.83 ms per loop
In [30]:
timeit dcumsum.dcumsum(a)

100000 loops, best of 3: 12.2 µs per loop
In [31]:
timeit a.cumsum()

10000 loops, best of 3: 27.4 µs per loop

## ctypes

ctypes is a Python library for calling out to C code. It is not as automatic as f2py, and we manually need to load the library and set properties such as the functions return and argument types. On the otherhand we do not need to touch the C code at all.

In [32]:
%%file functions.c

#include <stdio.h>

void hello(int n);

double dprod(double *x, int n);

void dcumsum(double *a, double *b, int n);

void
hello(int n)
{
int i;

for (i = 0; i < n; i++)
{
printf("C says hello\n");
}
}

double
dprod(double *x, int n)
{
int i;
double y = 1.0;

for (i = 0; i < n; i++)
{
y *= x[i];
}

return y;
}

void
dcumsum(double *a, double *b, int n)
{
int i;

b[0] = a[0];
for (i = 1; i < n; i++)
{
b[i] = a[i] + b[i-1];
}
}

Overwriting functions.c

Compile the C file into a shared library:

In [33]:
!gcc -c -Wall -O2 -Wall -ansi -pedantic -fPIC -o functions.o functions.c
!gcc -o libfunctions.so -shared functions.o


The result is a compiled shared library libfunctions.so:

In [34]:
!file libfunctions.so

libfunctions.so: ELF 64-bit LSB shared object, x86-64, version 1 (SYSV), dynamically linked, BuildID[sha1]=d68173ae6a804f703472af96f413b81a189db4b8, not stripped

Now we need to write wrapper functions to access the C library: To load the library we use the ctypes package, which included in the Python standard library (with extensions from numpy for passing arrays to C). Then we manually set the types of the argument and return values (no automatic code inspection here!).

In [35]:
%%file functions.py

import numpy
import ctypes

_libfunctions = numpy.ctypeslib.load_library('libfunctions', '.')

_libfunctions.hello.argtypes = [ctypes.c_int]
_libfunctions.hello.restype  =  ctypes.c_void_p

_libfunctions.dprod.argtypes = [numpy.ctypeslib.ndpointer(dtype=numpy.float), ctypes.c_int]
_libfunctions.dprod.restype  = ctypes.c_double

_libfunctions.dcumsum.argtypes = [numpy.ctypeslib.ndpointer(dtype=numpy.float), numpy.ctypeslib.ndpointer(dtype=numpy.float), ctypes.c_int]
_libfunctions.dcumsum.restype  = ctypes.c_void_p

def hello(n):
return _libfunctions.hello(int(n))

def dprod(x, n=None):
if n is None:
n = len(x)
x = numpy.asarray(x, dtype=numpy.float)
return _libfunctions.dprod(x, int(n))

def dcumsum(a, n):
a = numpy.asarray(a, dtype=numpy.float)
b = numpy.empty(len(a), dtype=numpy.float)
_libfunctions.dcumsum(a, b, int(n))
return b

Overwriting functions.py
In [36]:
%%file run_hello_c.py

import functions

functions.hello(3)

Overwriting run_hello_c.py
In [37]:
!python run_hello_c.py

C says hello C says hello C says hello
In [38]:
import functions


### Product function:

In [39]:
functions.dprod([1,2,3,4,5])

120.0

### Cummulative sum:

In [40]:
a = rand(100000)

In [41]:
res_c = functions.dcumsum(a, len(a))

In [42]:
res_fortran = dcumsum.dcumsum(a)

In [43]:
res_c - res_fortran

array([ 0., 0., 0., ..., 0., 0., 0.])

### Simple benchmark

In [44]:
timeit functions.dcumsum(a, len(a))

1000 loops, best of 3: 286 µs per loop
In [45]:
timeit dcumsum.dcumsum(a)

10000 loops, best of 3: 119 µs per loop
In [46]:
timeit a.cumsum()

1000 loops, best of 3: 261 µs per loop

## Cython

A hybrid between python and C that can be compiled: Basically Python code with type declarations.

In [47]:
%%file cy_dcumsum.pyx

cimport numpy

def dcumsum(numpy.ndarray[numpy.float64_t, ndim=1] a, numpy.ndarray[numpy.float64_t, ndim=1] b):
cdef int i, n = len(a)
b[0] = a[0]
for i from 1 <= i < n:
b[i] = b[i-1] + a[i]
return b

Overwriting cy_dcumsum.pyx

A build file for generating C code and compiling it into a Python module.

In [48]:
%%file setup.py

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

setup(
cmdclass = {'build_ext': build_ext},
ext_modules = [Extension("cy_dcumsum", ["cy_dcumsum.pyx"])]
)

Overwriting setup.py
In [49]:
!python setup.py build_ext --inplace

running build_ext cythoning cy_dcumsum.pyx to cy_dcumsum.c warning: /usr/local/lib/python2.7/dist-packages/Cython/Includes/numpy.pxd:869:17: Non-trivial type declarators in shared declaration (e.g. mix of pointers and values). Each pointer declaration should be on its own line. warning: /usr/local/lib/python2.7/dist-packages/Cython/Includes/numpy.pxd:869:24: Non-trivial type declarators in shared declaration (e.g. mix of pointers and values). Each pointer declaration should be on its own line. building 'cy_dcumsum' extension x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC -I/usr/include/python2.7 -c cy_dcumsum.c -o build/temp.linux-x86_64-2.7/cy_dcumsum.o In file included from /usr/include/python2.7/numpy/ndarraytypes.h:1761:0, from /usr/include/python2.7/numpy/ndarrayobject.h:17, from /usr/include/python2.7/numpy/arrayobject.h:4, from cy_dcumsum.c:352: /usr/include/python2.7/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp] #warning "Using deprecated NumPy API, disable it by " \ ^ In file included from /usr/include/python2.7/numpy/ndarrayobject.h:26:0, from /usr/include/python2.7/numpy/arrayobject.h:4, from cy_dcumsum.c:352: /usr/include/python2.7/numpy/__multiarray_api.h:1629:1: warning: ‘_import_array’ defined but not used [-Wunused-function] _import_array(void) ^ In file included from /usr/include/python2.7/numpy/ufuncobject.h:327:0, from cy_dcumsum.c:353: /usr/include/python2.7/numpy/__ufunc_api.h:241:1: warning: ‘_import_umath’ defined but not used [-Wunused-function] _import_umath(void) ^ x86_64-linux-gnu-gcc -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -D_FORTIFY_SOURCE=2 -g -fstack-protector --param=ssp-buffer-size=4 -Wformat -Werror=format-security build/temp.linux-x86_64-2.7/cy_dcumsum.o -o /home/rob/Desktop/scientific-python-lectures/cy_dcumsum.so
In [50]:
import cy_dcumsum

In [51]:
a = array([1,2,3,4], dtype=float)
b = empty_like(a)
cy_dcumsum.dcumsum(a,b)
b

array([ 1., 3., 6., 10.])
In [52]:
a = array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])

In [53]:
b = empty_like(a)
cy_dcumsum.dcumsum(a, b)
b

array([ 1., 3., 6., 10., 15., 21., 28., 36.])
In [54]:
py_dcumsum(a)

array([ 1., 3., 6., 10., 15., 21., 28., 36.])
In [55]:
a = rand(100000)
b = empty_like(a)

In [56]:
timeit py_dcumsum(a)

10 loops, best of 3: 50.1 ms per loop
In [57]:
timeit cy_dcumsum.dcumsum(a,b)

1000 loops, best of 3: 263 µs per loop

### Cython in the IPython notebook

When working with the IPython (especially in the notebook), there is a more convenient way of compiling and loading Cython code. Using the %%cython IPython magic (command to IPython), we can simply type the Cython code in a code cell and let IPython take care of the conversion to C code, compilation and loading of the function. To be able to use the %%cython magic, we first need to load the extension cythonmagic:

In [58]:
%load_ext cythonmagic

In [62]:
%%cython

cimport numpy

def cy_dcumsum2(numpy.ndarray[numpy.float64_t, ndim=1] a, numpy.ndarray[numpy.float64_t, ndim=1] b):
cdef int i, n = len(a)
b[0] = a[0]
for i from 1 <= i < n:
b[i] = b[i-1] + a[i]
return b

In [63]:
timeit cy_dcumsum2(a,b)

1000 loops, best of 3: 265 µs per loop

## Versions

In [64]:
%reload_ext version_information

%version_information ctypes, Cython

SoftwareVersion
Python2.7.6 (default, Mar 22 2014, 22:59:56) [GCC 4.8.2]
IPython1.1.0
OSposix [linux2]
ctypes1.1.0
Cython0.20.2
Tue Aug 26 23:37:29 2014 JST