%matplotlib inline

import numpy as np
import numexpr as ne
import numba
import math
import random
import matplotlib.pyplot as plt
import scipy as sp
import sys

%load_ext Cython

"""
Explanation: <h1 align="center">Scientific Programming in Python</h1>
<h2 align="center">Topic 5: Accelerating Python with Cython: Writting C in Python </h2>

Notebook created by Martín Villanueva - martin.villanueva@usm.cl - DI UTFSM - May2017.
End of explanation
"""


def primes_python(n):
    primes = [False, False] + [True] * (n - 2)
    i= 2
    while i < n:
        # We do not deal with composite numbers.
        if not primes[i]:
            i += 1
            continue 
        k= i+i
        # We mark multiples of i as composite numbers.
        while k < n:
            primes[k] = False
            k += i 
        i += 1
    # We return all numbers marked with True.
    return [i for i in range(2, n) if primes[i]]

primes_python(20)

"""
Explanation: Table of Contents

1.- Cython Basic Usage
2.- Advanced usage
3.- Pure C in Python

<div id='cython' />
1.- Cython Basic Usage
Cython is both a Superset of Python and a Python Library that lets you combine C and Python in various ways. There are two main use-cases:
1. Optimizing your Python code by statically compiling it to C.
2. Wrapping a C/C++ library in Python.
In order to get it properly working, you need Cython and a C compiler:
1. Cython: conda install cython
2. C compiler: Install GNU C compiler with your package manager (Unix/Linux) or install Xcode (OSX).

We will introduce the basic Cython usage by impementing the  Eratosthenes Sieve Algorithm, which is an algorithm to find all prime numbers smaller than a given number.
End of explanation
"""


tp = %timeit -o primes_python(10000)

"""
Explanation: Let's evaluate the performance for the first version:
End of explanation
"""


%%cython
def primes_cython1(n):
    primes = [False, False] + [True] * (n - 2)
    i= 2
    while i < n:
        # We do not deal with composite numbers.
        if not primes[i]:
            i += 1
            continue 
        k= i+i
        # We mark multiples of i as composite numbers.
        while k < n:
            primes[k] = False
            k += i 
        i += 1
    # We return all numbers marked with True.
    return [i for i in range(2, n) if primes[i]]

tc1 = %timeit -o primes_cython1(10000)

"""
Explanation: And now we write our first Cython version, by just adding %%cython magic in the first line of the cell:
End of explanation
"""


%%cython
def primes_cython2(int n):
    # Note the type declarations below
    cdef list primes = [False, False] + [True] * (n - 2)
    cdef int i = 2
    cdef int k = 0
    # The rest of the functions is unchanged
    while i < n:
        # We do not deal with composite numbers.
        if not primes[i]:
            i += 1
            continue 
        k= i+i
        # We mark multiples of i as composite numbers.
        while k < n:
            primes[k] = False
            k += i 
        i += 1
    # We return all numbers marked with True.
    return [i for i in range(2, n) if primes[i]]

tc2 = %timeit -o primes_cython2(10000)

print("Cython version 1 speedup: {0}".format(tp.best/tc1.best))
print("Cython version 2 speedup: {0}".format(tp.best/tc2.best))

"""
Explanation: We achieve x2 speed improvement doing (practically) nothing!.
When we add %%cython at the beginning of the cell, the code gets compiled by Cython into a C extension. Then, this extension is loaded, and the compiled function is readily available in the interactive namespace. 
Lets help the compiler by explicitly defining the type of the variables with the cdef macro/keyword:
End of explanation
"""


@numba.jit(nopython=True)
def primes_numba(n):
    primes = [False, False] + [True] * (n - 2)
    i= 2
    while i < n:
        # We do not deal with composite numbers.
        if not primes[i]:
            i += 1
            continue 
        k= i+i
        # We mark multiples of i as composite numbers.
        while k < n:
            primes[k] = False
            k += i 
        i += 1
    # We return all numbers marked with True.
    res = []
    for i in range(2,n):
        if primes[i]: res.append(i)
    return res

tn = %timeit -o primes_numba(10000)

"""
Explanation: Then: In general, Cython will be the most efficient when it can compile data structures and operations directly to C by __making as few CPython API calls as possible__. Specifying the types of the variables often leads to greater speed improvements.
Just for curiosity let's see the performance Numba's JIT achieves:
End of explanation
"""


%%cython -a
def primes_cython1(n):
    primes = [False, False] + [True] * (n - 2)
    i= 2
    while i < n:
        # We do not deal with composite numbers.
        if not primes[i]:
            i += 1
            continue 
        k= i+i
        # We mark multiples of i as composite numbers.
        while k < n:
            primes[k] = False
            k += i 
        i += 1
    # We return all numbers marked with True.
    return [i for i in range(2, n) if primes[i]]

%%cython -a
def primes_cython2(int n):
    # Note the type declarations below
    cdef list primes = [False, False] + [True] * (n - 2)
    cdef int i = 2
    cdef int k = 0
    # The rest of the functions is unchanged
    while i < n:
        # We do not deal with composite numbers.
        if not primes[i]:
            i += 1
            continue 
        k= i+i
        # We mark multiples of i as composite numbers.
        while k < n:
            primes[k] = False
            k += i 
        i += 1
    # We return all numbers marked with True.
    return [i for i in range(2, n) if primes[i]]

"""
Explanation: Numba wins this time! but: This is not the final form of Cython... 
Inspecting Cython bottlenecks with annotations
We can inspect the C code generated by Cython with the -a argument. Let's inspect the code used above.
The non-optimized lines will be shown in a gradient of yellow (white lines are faster, yellow lines are slower), telling you which lines are the least efficiently compiled to C. By clicking on a line, you can see the generated C code corresponding to that line.
End of explanation
"""


# Matrices to use
A = np.random.random((1000,3))
B = np.random.random((500,3))

def dist(a, b):
    return np.sqrt(np.sum((a-b)**2))

def distance_matrix_python(A, B):
    m = A.shape[0]
    n = B.shape[0]
    D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i],B[j])
    return D

%timeit distance_matrix_python(A,B)

%%cython -a
import numpy as np

def dist(a, b):
    return np.sqrt(np.sum((a-b)**2))

def distance_matrix_cython0(A, B):
    m = A.shape[0]
    n = B.shape[0]
    D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i],B[j])
    return D

%timeit distance_matrix_cython0(A,B)

"""
Explanation: Alternative usage of Cython: Outside the notebook
If you want to use Cython outside the notebook (the way it was thought...), you have to do the work of the magic:
1. Write the function into a .pyx file.
2. Cythonize it with cython filename.pyx generating the filename.c file.
3. Compile it with GCC: 
gcc -shared -fPIC -fwrapv -O3 -fno-strict-aliasing -I/home/mavillan/anaconda3/include/python3.6m -o primes.so primes.c
<div id='cython++' />
2.- Advanced usage
In this section we will consider the example of computing a distance matrix: Given the matrices $A_{m,3}$ and $B_{n,3}$ (each row is a 3D-position), the distance matrix has entries $D_{i,j} = d(A[i],B[j])$.
NumPy Arrays
You can use NumPy from Cython exactly the same as in regular Python, but by doing so you are losing potentially high speedups because Cython has support for fast access to NumPy arrays.
End of explanation
"""


%%cython -a
import numpy as np
cimport numpy as cnp

ctypedef cnp.float64_t float64_t

def dist(cnp.ndarray[float64_t, ndim=1] a, cnp.ndarray[float64_t, ndim=1] b):
    return np.sqrt(np.sum((a-b)**2))

def distance_matrix_cython1(cnp.ndarray[float64_t, ndim=2] A, cnp.ndarray[float64_t, ndim=2] B):
    cdef:
        int m = A.shape[0]
        int n = B.shape[0]
        int i,j
        cnp.ndarray[float64_t, ndim=2] D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i], B[j])
    return D

%timeit -n 10 distance_matrix_cython1(A,B)

%%cython -a
import numpy as np
cimport numpy as cnp

ctypedef cnp.float64_t float64_t
from libc.math cimport sqrt

def dist(cnp.ndarray[float64_t, ndim=1] a, cnp.ndarray[float64_t, ndim=1] b):
    cdef:
        int i = 0
        int n = a.shape[0]
        float ret = 0
    for i in range(n):
        ret += (a[i]-b[i])**2
    return sqrt(ret)


def distance_matrix_cython2(cnp.ndarray[float64_t, ndim=2] A, cnp.ndarray[float64_t, ndim=2] B):
    cdef:
        int m = A.shape[0]
        int n = B.shape[0]
        int i,j
        cnp.ndarray[float64_t, ndim=2] D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i], B[j])
    return D

%timeit -n 10 distance_matrix_cython2(A,B)

"""
Explanation: Now let's improve this naive Cython implementation by statically defining the types of the variables:
End of explanation
"""


%%cython -a
import numpy as np
cimport numpy as cnp

ctypedef cnp.float64_t float64_t
from libc.math cimport sqrt

def dist(float64_t[::1] a, float64_t[::1] b):
    cdef:
        int i = 0
        int n = a.shape[0]
        float ret = 0
    for i in range(n):
        ret += (a[i]-b[i])**2
    return sqrt(ret)


def distance_matrix_cython3(float64_t[:,::1] A, float64_t[:,::1] B):
    cdef:
        int m = A.shape[0]
        int n = B.shape[0]
        int i,j
        float64_t[:,::1] D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i], B[j])
    return D

%timeit -n 10 distance_matrix_cython3(A,B)

"""
Explanation: Typed Memory Views
Typed memoryviews allow efficient access to memory buffers, such as those underlying NumPy arrays, without incurring any Python overhead. Memoryviews are similar to the current NumPy array buffer support (np.ndarray[np.float64_t, ndim=2]), but they have more features and cleaner syntax.
They can handle a wider variety of sources of array data. For example, they can handle C arrays and the Cython array type (Cython arrays).
Syntaxis: dtype[:,::1] where ::1 indicates the axis where elements are contiguous.
End of explanation
"""


%%cython -a -c=-fPIC -c=-fwrapv -c=-O3 -c=-fno-strict-aliasing
import numpy as np
cimport numpy as cnp

ctypedef cnp.float64_t float64_t
from libc.math cimport sqrt

def dist(float64_t[::1] a, float64_t[::1] b):
    cdef:
        int i = 0
        int n = a.shape[0]
        float ret = 0
    for i in range(n):
        ret += (a[i]-b[i])**2
    return sqrt(ret)


def distance_matrix_cython4(float64_t[:,::1] A, float64_t[:,::1] B):
    cdef:
        int m = A.shape[0]
        int n = B.shape[0]
        int i,j
        float64_t[:,::1] D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i], B[j])
    return D

%timeit -n 10 distance_matrix_cython4(A,B)

"""
Explanation: Compiler optimization
With the -c option we can pass the compiler (gcc) optimization options. Below we use the most common of them:
End of explanation
"""


%%cython -a -c=-fPIC -c=-fwrapv -c=-O3 -c=-fno-strict-aliasing
#!python
#cython: cdivision=True, boundscheck=False, nonecheck=False, wraparound=False, initializedcheck=False

import numpy as np
cimport numpy as cnp

ctypedef cnp.float64_t float64_t
from libc.math cimport sqrt

def dist(float64_t[::1] a, float64_t[::1] b):
    cdef:
        int i = 0
        int n = a.shape[0]
        float ret = 0
    for i in range(n):
        ret += (a[i]-b[i])**2
    return sqrt(ret)


def distance_matrix_cython5(float64_t[:,::1] A, float64_t[:,::1] B):
    cdef:
        int m = A.shape[0]
        int n = B.shape[0]
        int i,j
        float64_t[:,::1] D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i,:], B[j,:])
    return D

%timeit -n 10 distance_matrix_cython5(A,B)

"""
Explanation: Compiler directives
 Compiler directives are instructions which affect the behavior of Cython code. 
cdivision (True / False)
* If set to False, Cython will adjust the remainder and quotient operators C types to match those of Python ints (which differ when the operands have opposite signs) and raise a ZeroDivisionError when the right operand is 0. This has up to a 35% speed penalty. If set to True.
boundscheck (True / False)
* If set to False, Cython is free to assume that indexing operations ([]-operator) in the code will not cause any IndexErrors to be raised. Lists, tuples, and strings are affected only if the index can be determined to be non-negative (or if wraparound is False). Conditions which would normally trigger an IndexError may instead cause segfaults or data corruption if this is set to False. Default is True.
nonecheck (True / False)
* If set to False, Cython is free to assume that native field accesses on variables typed as an extension type, or buffer accesses on a buffer variable, never occurs when the variable is set to None. Otherwise a check is inserted and the appropriate exception is raised. This is off by default for performance reasons. Default is False.
wraparound (True / False)
* In Python arrays can be indexed relative to the end. For example A[-1] indexes the last value of a list. In C negative indexing is not supported. If set to False, Cython will neither check for nor correctly handle negative indices, possibly causing segfaults or data corruption. Default is True.
initializedcheck (True / False)
* If set to True, Cython checks that a memoryview is initialized whenever its elements are accessed or assigned to. Setting this to False disables these checks. Default is True.
For all the compilation directives see here.
End of explanation
"""


%%cython -a -c=-fPIC -c=-fwrapv -c=-O3 -c=-fno-strict-aliasing
#!python
#cython: cdivision=True, boundscheck=False, nonecheck=False, wraparound=False, initializedcheck=False

import numpy as np
cimport numpy as cnp

ctypedef cnp.float64_t float64_t
from libc.math cimport sqrt

cdef float64_t dist(float64_t[::1] a, float64_t[::1] b):
    cdef:
        int i = 0
        int n = a.shape[0]
        float ret = 0
    for i in range(n):
        ret += (a[i]-b[i])**2
    return sqrt(ret)


def distance_matrix_cython6(float64_t[:,::1] A, float64_t[:,::1] B):
    cdef:
        int m = A.shape[0]
        int n = B.shape[0]
        int i,j
        float64_t[:,::1] D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i,:], B[j,:])
    return D

%timeit -n 10 distance_matrix_cython6(A,B)

"""
Explanation: Pure C functions
With the cdef keywork we can realy define C function, as we shown below. In such functions all variable types should be defined and should have a return type, and can't be called directly in Python, i.e, only can be called by functions defined in the same module.
There is a midpoint between def and cdef which automatically creates a Python function with the same name, so the function can be called directly.
End of explanation
"""


%%cython  -c=-fPIC -c=-fwrapv -c=-O3 -c=-fno-strict-aliasing
#!python
#cython: cdivision=True, boundscheck=False, nonecheck=False, wraparound=False, initializedcheck=False

cimport numpy as cnp
ctypedef cnp.float64_t float64_t

cdef float64_t test1(float64_t a, float64_t b):
    return a+b

test1(1.,1.)

%%cython  -c=-fPIC -c=-fwrapv -c=-O3 -c=-fno-strict-aliasing
#!python
#cython: cdivision=True, boundscheck=False, nonecheck=False, wraparound=False, initializedcheck=False

cimport numpy as cnp
ctypedef cnp.float64_t float64_t

cpdef float64_t test2(float64_t a, float64_t b):
    return a+b

test2(1,1)

"""
Explanation: Example of cdef and cpdef
End of explanation
"""


%%cython -a -c=-fPIC -c=-fwrapv -c=-O3 -c=-fno-strict-aliasing
#!python
#cython: cdivision=True, boundscheck=False, nonecheck=False, wraparound=False, initializedcheck=False

import numpy as np
cimport numpy as cnp

ctypedef cnp.float64_t float64_t
from libc.math cimport sqrt

cdef inline float64_t dist(float64_t[::1] a, float64_t[::1] b):
    cdef:
        int i = 0
        int n = a.shape[0]
        float ret = 0
    for i in range(n):
        ret += (a[i]-b[i])**2
    return sqrt(ret)

def distance_matrix_cython7(float64_t[:,::1] A, float64_t[:,::1] B):
    cdef:
        int m = A.shape[0]
        int n = B.shape[0]
        int i,j
        float64_t[:,::1] D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i,:], B[j,:])
    return D

%timeit -n 10 distance_matrix_cython7(A,B)

"""
Explanation: Function inlining
In computing, inline expansion, or inlining, is a manual or compiler optimization that replaces a function call site with the body of the called function. 
As a rule of thumb: Some inlining will improve speed at very minor cost of space, but excess inlining will hurt speed, due to inlined code consuming too much of the instruction cache, and also cost significant space.
End of explanation
"""


@numba.jit(nopython=True)
def dist(a, b):
    n = a.shape[0]
    ret = 0
    for i in range(n):
        ret += (a[i]-b[i])**2
    return math.sqrt(ret)

@numba.jit(nopython=True)
def distance_matrix_numba(A, B):
    m = A.shape[0]
    n = B.shape[0]
    D = np.empty((m,n))
    for i in range(m):
        for j in range(n):
            D[i,j] = dist(A[i,:], B[j,:])
    return D

%timeit -n 10 distance_matrix_numba(A,B)

"""
Explanation: What about Numba?
End of explanation
"""


%%cython -a -c=-fPIC -c=-fwrapv -c=-O3 -c=-fno-strict-aliasing
#!python
#cython: cdivision=True, boundscheck=False, nonecheck=False, wraparound=False, initializedcheck=False

cdef class A(object):
    def   d(self): return 0
    cdef  int c(self): return 0
    cpdef int p(self): return 0

    def test_def(self, long num):
        while num > 0:
            self.d()
            num -= 1

    def test_cdef(self, long num):
        while num > 0:
            self.c()
            num -= 1

    def test_cpdef(self, long num):
        while num > 0:
            self.p()
            num -= 1

%%timeit n = 1000000
a1 = A()
a1.test_def(n)

%%timeit n = 1000000
a1 = A()
a1.test_cdef(n)

%%timeit n = 1000000
a1 = A()
a1.test_cpdef(n)

"""
Explanation: 4.- Other advanced things you can do with Cython
We have seen that with Cython we can implement our algorithms achieving C performance. Moreover it is very versatile and we can do some other advanced thing with it:
Object-oriented programming: Classes and methods
To support object-oriented programming, Cython supports writing normal Python classes exactly as in Python.
End of explanation
"""


