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)
1.74 ms ± 53.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
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_be1788cc2bbd7f7fd1a4044281992e09.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
Out[11]:
<function _cython_magic_4a4ac9ad4d385768791792f92e54ab70.fib1>
In [12]:
fib1(5)
Out[12]:
13
In [13]:
%timeit fib1(10000)
1.59 ms ± 31.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [14]:
%%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 [15]:
fib2
Out[15]:
<function _cython_magic_1a42b2c15adc8f8d5baea688dbf17820.fib2>
In [16]:
fib2(5)
Out[16]:
13
In [17]:
%timeit fib2(10000)
2.24 µs ± 17 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Second Example - Numerical integration

In [18]:
%%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 [19]:
from math import pi
In [20]:
integrate0(0.0,pi,1000000)
Out[20]:
0.5216069591323289
In [21]:
%timeit integrate0(0.0,pi,1000000)
101 ms ± 3.51 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [22]:
%%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 [23]:
%timeit integrate1(0.0,pi,1000000)
69.7 ms ± 1.32 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [24]:
f1(5)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-24-c5ed7dd0edc1> in <module>
----> 1 f1(5)

NameError: name 'f1' is not defined
In [25]:
%%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 [26]:
%timeit integrate2(0.0,pi,1000000)
12.2 ms ± 150 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [27]:
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 [28]:
%timeit integrate3(0.0,pi,1000000)
19.5 ms ± 872 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Cython and STL containers

In [29]:
%%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 [30]:
manipulateVector([1,2,3,4])
Out[30]:
[1.0, 2.0, 3.0, 4.0, 42.0]

Cython and exceptions

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

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

Back to integration - Cython and Classes

In [38]:
%%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 [39]:
f = SinExpFunction()
In [40]:
import math
In [41]:
%timeit integrate(f,0,math.pi,1000000)
14.9 ms ± 65.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [42]:
class Poly(Integrand):
    def evaluate(self,x):
        return 3*x-x**3
In [43]:
g = Poly()
integrate(g,0,math.pi,10000)
Out[43]:
-9.5478660351052
In [44]:
%timeit integrate(g,0,math.pi,1000000)
157 ms ± 1.44 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [ ]: