# 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 [ ]: