{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# A Python Tour of Data Science: High Performance Computing\n",
"\n",
"[Michaël Defferrard](http://deff.ch), *PhD student*, [EPFL](http://epfl.ch) [LTS2](http://lts2.epfl.ch)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exercise: is Python slow ?\n",
"\n",
"That is one of the most heard complain about Python. Because [CPython](https://en.wikipedia.org/wiki/CPython), the Python reference implementation, [interprets the language](https://en.wikipedia.org/wiki/Interpreted_language) (i.e. it compiles Python code to intermediate bytecode which is then interpreted by a virtual machine), it is inherentably slower than [compiled languages](https://en.wikipedia.org/wiki/Compiled_language), especially for computation heavy tasks such as number crunching.\n",
"\n",
"There are three ways around it that we'll explore in this exercise:\n",
"1. Specialized libraries.\n",
"1. Compile Python to machine code.\n",
"1. Implement in a compiled language and call from Python.\n",
"\n",
"In this exercise we'll compare many possible implementations of a function. Our goal is to compare the execution time of our implementations and get a sense of the many ways to write efficient Python code. We test seven implementations:\n",
"\n",
"Along the exercise we'll use the function\n",
"$$accuracy(\\hat{y}, y) = \\frac1n \\sum_{i=0}^{i=n-1} 1(\\hat{y}_i = y_i),$$\n",
"where $1(x)$ is the [indicator function](https://en.wikipedia.org/wiki/Indicator_function) and $n$ is the number of samples. This function computes the accuracy, i.e. the percentage of correct predictions, of a classifier. A pure Python implementation is given below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def accuracy_python(y_pred, y_true):\n",
" \"\"\"Plain Python implementation.\"\"\"\n",
" num_correct = 0\n",
" for y_pred_i, y_true_i in zip(y_pred, y_true):\n",
" if y_pred_i == y_true_i:\n",
" num_correct += 1\n",
" return num_correct / len(y_true)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below we test and measure the execution time of the above implementation. The [%timeit](http://ipython.readthedocs.io/en/stable/interactive/magics.html?highlight=timeit#magic-timeit) function provided by IPython is a useful helper to measure the execution time of a line of Python code. As we'll see, the above implementation is very inefficient compared to what we can achieve."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"c = 10 # Number of classes.\n",
"n = int(1e6) # Number of samples.\n",
"\n",
"y_true = np.random.randint(0, c, size=n)\n",
"y_pred = np.random.randint(0, c, size=n)\n",
"print('Expected accuracy: {}'.format(1/c))\n",
"print('Empirical accuracy: {}'.format(accuracy_python(y_pred, y_true)))\n",
"\n",
"%timeit accuracy_python(y_pred, y_true)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 1 Specialized libraries\n",
"\n",
"Specialized libraries, which provide efficient compiled implementations of the heavy computations, is an easy way to solve the performance problem. That is for example NumPy, which uses efficient BLAS and LAPACK implementations as a backend. SciPy and scikit-learn fall in the same category.\n",
"\n",
"Implement below the accuracy function using:\n",
"1. Only functions provided by [NumPy](http://www.numpy.org/). The idea here is to *vectorize* the computation.\n",
"2. The implementation of [scikit-learn](http://scikit-learn.org/), our machine learning library.\n",
"\n",
"Then test that it provides the correct result and measure it's execution time. How much faster are they compared to the pure Python implementation ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def accuracy_numpy(y_pred, y_true):\n",
" \"\"\"Numpy implementation.\"\"\"\n",
" return np.sum(y_pred == y_true) / y_true.size\n",
"\n",
"def accuracy_sklearn(y_pred, y_true):\n",
" \"\"\"Scikit-learn implementation.\"\"\"\n",
" from sklearn import metrics\n",
" return metrics.accuracy_score(y_pred, y_true)\n",
"\n",
"assert np.allclose(accuracy_numpy(y_pred, y_true), accuracy_python(y_pred, y_true))\n",
"assert np.allclose(accuracy_sklearn(y_pred, y_true), accuracy_python(y_pred, y_true))\n",
"\n",
"%timeit accuracy_numpy(y_pred, y_true)\n",
"%timeit accuracy_sklearn(y_pred, y_true)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 2 Compiled Python\n",
"\n",
"The second option of choice, when the algorirthm does not exist in our favorite libraries and we have to implement it, it to implement in Python and compile it to machine code.\n",
"\n",
"Below you'll compile Python with two frameworks.\n",
"1. [Numba](http://numba.pydata.org) is a [just-in-time (JIT)](https://en.wikipedia.org/wiki/Just-in-time_compilation) compiler for Python, using the [LLVM](http://llvm.org) compiler infrastructure.\n",
"1. [Cython](http://cython.org), which requires type information, [transpiles](https://en.wikipedia.org/wiki/Source-to-source_compiler) Python to C then compiles the generated C code.\n",
"\n",
"While these two approaches offer maximal compatibility with the CPython and NumPy ecosystems, another approach is to use another Python implementation such as [PyPy](http://pypy.org), which features a just-in-time compiler and supports multiple back-ends (C, CLI, JVM). Alternatives are [Jython](http://www.jython.org), which runs Python on the Java platform, and [IronPython](http://ironpython.net) / [PythonNet](http://pythonnet.github.io) for the .NET platform."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from numba import jit\n",
"# Decorator, same as accuracy_numba = jit(accuracy_python)\n",
"\n",
"@jit\n",
"def accuracy_numba(y_pred, y_true):\n",
" \"\"\"Plain Python implementation, compiled by LLVM through Numba.\"\"\"\n",
" num_correct = 0\n",
" for y_pred_i, y_true_i in zip(y_pred, y_true):\n",
" if y_pred_i == y_true_i:\n",
" num_correct += 1\n",
" return num_correct / len(y_true)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%load_ext Cython"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%cython\n",
"cimport numpy as np\n",
"cimport cython\n",
"\n",
"@cython.boundscheck(False) # Turn off bounds-checking for entire function.\n",
"@cython.wraparound(False) # Turn off negative index wrapping for entire function.\n",
"def accuracy_cython(np.ndarray[long, ndim=1] y_pred, np.ndarray[long, ndim=1] y_true):\n",
" \"\"\"Python implementation with type information, transpiled to C by Cython.\"\"\"\n",
" cdef int num_total = len(y_true)\n",
" cdef int num_correct = 0\n",
" cdef int n = y_pred.size\n",
" for i in range(n):\n",
" if y_pred[i] == y_true[i]:\n",
" num_correct += 1\n",
" return num_correct / num_total"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Evaluate below the performance of those two implementations, while testing their correctness. How do they compare with plain Python and specialized libraries ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"assert np.allclose(accuracy_numba(y_pred, y_true), accuracy_python(y_pred, y_true))\n",
"assert np.allclose(accuracy_cython(y_pred, y_true), accuracy_python(y_pred, y_true))\n",
"\n",
"%timeit accuracy_numba(y_pred, y_true)\n",
"%timeit accuracy_cython(y_pred, y_true)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3 Using C from Python\n",
"\n",
"Here we'll explore our third option to make Python faster: implement in another language ! Below you'll\n",
"1. implement the accuracy function in C,\n",
"1. compile it, e.g. with the GNU compiler collection ([GCC](https://gcc.gnu.org/)),\n",
"1. call it from Python."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%file function.c\n",
"\n",
"double accuracy(long* y_pred, long* y_true, int n) {\n",
" int num_correct = 0;\n",
"\n",
" for(int i = 0; i < n; i++) {\n",
" if(y_pred[i] == y_true[i])\n",
" num_correct++;\n",
" }\n",
"\n",
" return (double) num_correct / n;\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The below cell describe a shell script, which will be executed by IPython as if you typed the commands in your terminal. Those commands are compiling the above C program into a dynamic library with GCC. You can use any other compiler, text editor or IDE to produce the C library. Windows users, you may want to use Microsoft toolchain."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%script sh\n",
"FILE=function\n",
"gcc -c -O3 -Wall -std=c11 -pedantic -fPIC -o $FILE.o $FILE.c\n",
"gcc -o lib$FILE.so -shared $FILE.o\n",
"file lib$FILE.so"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The below cell finally create a wrapper around our C library so that we can easily use it from Python."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import ctypes\n",
"\n",
"libfunction = np.ctypeslib.load_library('libfunction', './')\n",
"libfunction.accuracy.restype = ctypes.c_double\n",
"libfunction.accuracy.argtypes = [\n",
" np.ctypeslib.ndpointer(dtype=np.int),\n",
" np.ctypeslib.ndpointer(dtype=np.int),\n",
" ctypes.c_int\n",
"]\n",
"\n",
"def accuracy_c(y_pred, y_true):\n",
" n = y_pred.size\n",
" return libfunction.accuracy(y_pred, y_true, n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Evaluate below the performance of your C implementation, and test its correctness. How does it compare with the others ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"assert np.allclose(accuracy_c(y_pred, y_true), accuracy_python(y_pred, y_true))\n",
"%timeit accuracy_c(y_pred, y_true)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4 Using Fortran from Python\n",
"\n",
"Same idea as before, with Fortran ! [Fortran](https://en.wikipedia.org/wiki/Fortran) is an imperative programming language developed in the 50s, especially suited to numeric computation and scientific computing. While you probably won't write new code in Fortran, you may have to interface with legacy code (especially in large and old corporations). Here we'll resort to the [f2py](https://docs.scipy.org/doc/numpy-dev/f2py/) utility provided by the Numpy project for the (almost automatic) generation of a wrapper."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%file function.f\n",
"\n",
"! The content of this cell is written to the \"function.f\" file in the current directory.\n",
"\n",
" SUBROUTINE DACCURACY(YPRED, YTRUE, ACC, N)\n",
"\n",
"CF2PY INTENT(OUT) :: ACC\n",
"CF2PY INTENT(HIDE) :: N\n",
" INTEGER*4 YPRED(N)\n",
" INTEGER*4 YTRUE(N)\n",
" DOUBLE PRECISION ACC\n",
" INTEGER N, NCORRECT\n",
"\n",
" NCORRECT = 0\n",
"\n",
" DO 10 J = 1, N\n",
" IF (YPRED(J) == YTRUE(J)) THEN\n",
" NCORRECT = NCORRECT + 1\n",
" END IF\n",
" 10 CONTINUE\n",
"\n",
" ACC = REAL(NCORRECT) / N\n",
" END"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The below command compile the Fortran code and generate a Python wrapper."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"!f2py -c -m function function.f >> /dev/null"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import function\n",
"def accuracy_fortran(y_pred, y_true):\n",
" return function.daccuracy(y_pred, y_true)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Evaluate below the performance of your Fortran implementation, and test its correctness. How does it compare with the others ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"assert np.allclose(accuracy_fortran(y_pred, y_true), accuracy_python(y_pred, y_true))\n",
"%timeit accuracy_fortran(y_pred, y_true)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 5 Analysis\n",
"\n",
"Plot a graph with `n` as the x-axis and the execution time of the various methods on the y-axis. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%capture\n",
"# Suppress outputs.\n",
"\n",
"N = 8 # Number of samples: from 1 to 10^N.\n",
"\n",
"methods = ['accuracy_python', 'accuracy_sklearn', 'accuracy_fortran', 'accuracy_numpy',\n",
" 'accuracy_numba', 'accuracy_cython', 'accuracy_c']\n",
"\n",
"times = np.empty((len(methods), N+1))\n",
"for n in range(0, N+1):\n",
" y_true = np.random.randint(2, size=10**n)\n",
" y_pred = np.random.randint(2, size=len(y_true))\n",
" for i, method in enumerate(methods):\n",
" res = %timeit -o eval(method)(y_pred, y_true)\n",
" times[i, n] = res.best"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"plt.figure(figsize=(15, 5))\n",
"for i, method in enumerate(methods):\n",
" plt.semilogy(range(N+1), times[i,:], '.-', label=method)\n",
"plt.legend(loc='upper left')\n",
"plt.ylabel('execution time [s]')\n",
"plt.xlabel('log_10 number of samples')\n",
"plt.xticks(range(N+1))\n",
"plt.show();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To reflect on what we've done, answer the following questions:\n",
"1. Which implementation is the fastest ?\n",
"1. Which solution was the easiest to work with ?\n",
"1. Which solution would you use in which situation ?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It turns out the compiled versions by Numba and Cython are almost as fast as C ! Although much easier to write. In this case, they are even faster than Fortran, NumPy and scikit-learn. This gives us the best of both worlds: an interpreted language for rapid development, which can then be compiled for efficient execution in production."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}