Cython

Christian Elsasser (che@physik.uzh.ch)

The content of this lecture might be distributed under CC by-sa.

An introductionnary example - Fibonacci numbers

In [1]:
def fib0(n):
    a,b=1,1
    for i in range(n):
        a,b=a+b,a
    return a
In [2]:
fib0(5)
Out[2]:
13
In [3]:
fib0(6)
Out[3]:
21
In [4]:
%timeit fib0(10000)
100 loops, best of 3: 3.39 ms per loop
In [5]:
import cython
In [6]:
%load_ext Cython
In [7]:
%%cython
def helloWorld():
    print("Hello, World!")
In [8]:
helloWorld
Out[8]:
<function _cython_magic_a9ecb934c282c3ed435f053056cba526.helloWorld>
In [9]:
helloWorld()
Hello, World!
In [10]:
%%cython
def fib1(n):
    a,b=1,1
    for i in range(n):
        a,b=a+b,a
    return a
In [11]:
fib1(5)
Out[11]:
13
In [12]:
%timeit fib1(10000)
100 loops, best of 3: 2.72 ms per loop
In [13]:
%%cython
def fib2(int n):
    cdef int a,b,i
    a,b=1,1
    for i in range(n):
        a,b=a+b,a
    return a
In [14]:
fib2(5)
Out[14]:
13
In [15]:
%timeit fib2(10000)
100000 loops, best of 3: 3.43 µs per loop

Second Example - Numerical integration

In [16]:
%%cython

from math import sin,exp

def f0(double x):
    return sin(x)*exp(-x)

def integrate0(double a,double b,int N):
    cdef double dx,s
    cdef int i
    dx = (b-a)/N
    s = 0.0
    for i in range(N):
        s+=f0(a+(i+0.5)*dx)
    return s*dx
In [17]:
from math import pi
In [18]:
integrate0(0.0,pi,1000000)
Out[18]:
0.5216069591323289
In [19]:
%timeit integrate0(0.0,pi,1000000)
1 loop, best of 3: 263 ms per loop
In [20]:
%%cython

from math import sin,exp

cdef double f1(double x):
    return sin(x)*exp(-x)

def integrate1(double a,double b,int N):
    cdef double dx,s
    cdef int i
    dx = (b-a)/N
    s = 0.0
    for i in range(N):
        s+=f1(a+(i+0.5)*dx)
    return s*dx
In [21]:
%timeit integrate1(0.0,pi,1000000)
10 loops, best of 3: 186 ms per loop
In [22]:
f1(5)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-22-c6005772367a> in <module>()
----> 1 f1(5)

NameError: name 'f1' is not defined
In [23]:
%%cython

from libc.math cimport sin,exp

cpdef double f2(double x):
    return sin(x)*exp(-x)

def integrate2(double a,double b,int N):
    cdef double dx,s
    cdef int i
    dx = (b-a)/N
    s = 0.0
    for i in range(N):
        s+=f2(a+(i+0.5)*dx)
    return s*dx
In [24]:
%timeit integrate2(0.0,pi,1000000)
10 loops, best of 3: 25.4 ms per loop
In [25]:
import numpy as np

def f3(x):
    return np.sin(x)*np.exp(-x)

def integrate3(a,b,N):
    dx = (b-a)/N
    x = (np.arange(N)+0.5)*dx+a
    return np.sum(f3(x))*dx
In [26]:
%timeit integrate3(0.0,pi,1000000)
10 loops, best of 3: 46.7 ms per loop

Cython and STL containers

In [27]:
%%cython
#distutils: language = c++
from libcpp.vector cimport vector
def manipulateVector(vector[double] v):
    v.push_back(42.0) # Equivalent to append for list in python
    return v
In [28]:
manipulateVector([1,2,3,4])
Out[28]:
[1.0, 2.0, 3.0, 4.0, 42.0]

Cython and exceptions

In [29]:
%%cython
cpdef int raiseError():
    raise RuntimeError("A problem!")
    return 1
In [30]:
raiseError()
Exception ignored in: '_cython_magic_fa7a72687a7c423eabc376188b792d18.raiseError'
RuntimeError: A problem!
Out[30]:
0
In [31]:
%%cython
cpdef int raiseErrorBetter() except *:
    raise RuntimeError("A problem!")
    return 1
In [32]:
try:
    raiseErrorBetter()
except RuntimeError:
    print("Error caught")
Error caught
In [33]:
%%cython
cpdef int raiseException(int n) except -1:
    # Assuming C function returning an integer; if -1 --> Error appeared
    return n
In [34]:
raiseException(1)
Out[34]:
1
In [35]:
raiseException(-1)
---------------------------------------------------------------------------
SystemError                               Traceback (most recent call last)
<ipython-input-35-b5db4e46ee95> in <module>()
----> 1 raiseException(-1)

SystemError: <built-in function raiseException> returned NULL without setting an error

Back to integration - Cython and Classes

In [36]:
%%cython

from libc.math cimport sin,exp

cdef class Integrand(object):
    cpdef double evaluate(self,double x):
        raise NotImplementedError()
        
cdef class SinExpFunction(Integrand):
    cpdef double evaluate(self,double x):
        return sin(x)*exp(-x)
    
def integrate(Integrand f, double a, double b, int N):
    cdef int i
    cdef double dx,s
    dx = (b-a)/N
    s = 0.0
    for i in range(N):
        s += f.evaluate(a+(i+0.5)*dx)
    return s*dx
In [37]:
f = SinExpFunction()
In [38]:
import math
In [39]:
%timeit integrate(f,0,math.pi,10000)
1000 loops, best of 3: 265 µs per loop
In [40]:
class Poly(Integrand):
    def evaluate(self,x):
        return 3*x-x**3
In [41]:
g = Poly()
integrate(g,0,math.pi,10000)
Out[41]:
-9.5478660351052
In [42]:
%timeit integrate(g,0,math.pi,10000)
100 loops, best of 3: 3.18 ms per loop