accelerating_python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Accelerating your Python code interactively"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lesson you will learn\n",
    "* how to speed up your Python code with **Numba**\n",
    "* how to speed up your Python code with **Cython**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are various ways to speed up the execution of code in a Jupyter notebook.  \n",
    "We will start scratching the surface by looking at Numba and Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other ways to increase performance include\n",
    "\n",
    "- **Numexpr** - a package that speeds up complex array operations (sometimes slow in numpy due to creation of many temporary arrays)\n",
    "- wrapping C libraries with **ctypes**\n",
    "- **CFFI** (C Foreign Function Interface for Python)\n",
    "- **SWIG** (Simplified Wrapper and Interface Generator)\n",
    "- **ipyparallel** (distributing parallel tasks on multicore computer)\n",
    "- creating a new magic command which compiles C/C++/Fortran code, runs it and returns result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that Cython and CFFI will be covered in more detail in the [MMA session](https://github.com/coderefinery/mma)! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy offers fast vector computations of arrays, but some algorithms can't be vectorized and require explicit loops.  \n",
    "However, Python loops are slow. `Numba` can speed up Python through a *just-in-time* (JIT) compiler, compiling Python code directly to machine code. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at an example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random walk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first import everything"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"red\"> *Exercise 3.1*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Load the random_walk.py example into a cell below\n",
    "* split up the functions over cells (either via Edit menu or keyboard shortcut). \n",
    "* The initialization of `n` and the call to `walk()` don't need to be in a function, and you can remove the `__name__` stuff\n",
    "* plot the random walk trajectory\n",
    "* time the execution of `walk()`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### <font color=\"green\"> *Solution*\n",
    "\n",
    "> The solution can be found in the solutions.ipynb notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Just-in-time* compilation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now import the Numba package, and then add the decorator `@jit` immediately before the function definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numba import jit, vectorize, float64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def step_numba():\n",
    "    return 1. if random.random() > .5 else -1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def walk_numba(n):\n",
    "    x = np.zeros(n)\n",
    "    dx = 1. / n\n",
    "    for i in range(n - 1):\n",
    "        x_new = x[i] + dx * step_numba()\n",
    "        if x_new > 5e-3:\n",
    "            x[i + 1] = 0.\n",
    "        else:\n",
    "            x[i + 1] = x_new\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `nopython=True` argument activates nopython mode, which means that the code is translated directly to machine code, bypassing the CPython interpreter. The nopython mode is faster than python mode, but is more limited and for example lists and dictionaries are not supported. For optimal performance, try to stick with nopython!  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`n` is still 10000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t2 = %timeit -o walk_numba(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How much faster was the code with just-in-time compilation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print (t1.best/t2.best)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quite a decent speedup!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cython is a superset of Python which supports calling C functions and declaring C types on variables and class attributes.\n",
    "Cython allows you to\n",
    "- wrap C/C++ libraries into Python\n",
    "- optimize Python code by statically compiling with C\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will optimize Python by \"cythonizing\" it. For wrapping C/C++ libraries, see http://docs.cython.org/en/latest/src/tutorial/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can also give compiler flags: `%%cython --compile-args=-fopenmp --link-args=-fopenmp`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Mandelbrot fractal "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "size = 200\n",
    "iterations = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pure python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mandelbrot_python(m, size, iterations):\n",
    "    for i in range(size):\n",
    "        for j in range(size):\n",
    "            c = -2 + 3./size*j + 1j*(1.5-3./size*i)\n",
    "            z = 0\n",
    "            for n in range(iterations):\n",
    "                if np.abs(z) <= 10:\n",
    "                    z = z*z + c\n",
    "                    m[i, j] = n\n",
    "                else:\n",
    "                    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1 -o \n",
    "m = np.zeros((size, size))\n",
    "mandelbrot_python(m, size, iterations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_py = _"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### First cython attempt\n",
    "\n",
    "First just add the cython magic. The `-a` option activates annotation, which shows \"how much Python\" there's left in the cythonized code (less yellow means less Python and faster execution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "import numpy as np\n",
    "\n",
    "def mandelbrot_cython(m, size, iterations):\n",
    "    for i in range(size):\n",
    "        for j in range(size):\n",
    "            c = -2 + 3./size*j + 1j*(1.5-3./size*i)\n",
    "            z = 0\n",
    "            for n in range(iterations):\n",
    "                if np.abs(z) <= 10:\n",
    "                    z = z*z + c\n",
    "                    m[i, j] = n\n",
    "                else:\n",
    "                    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1 -o\n",
    "m = np.zeros((size, size), dtype=np.int32)\n",
    "mandelbrot_cython(m, size, iterations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_cy0 = _"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Small speedup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Second attempt\n",
    "\n",
    "Now add type information, use *memory views* for NumPy arrays, and replace `np.abs()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "import numpy as np\n",
    "\n",
    "def mandelbrot_cython(int[:,::1] m, \n",
    "                      int size, \n",
    "                      int iterations):\n",
    "    cdef int i, j, n\n",
    "    cdef complex z, c\n",
    "    for i in range(size):\n",
    "        for j in range(size):\n",
    "            c = -2 + 3./size*j + 1j*(1.5-3./size*i)\n",
    "            z = 0\n",
    "            for n in range(iterations):\n",
    "                if z.real**2 + z.imag**2 <= 100:\n",
    "                    z = z*z + c\n",
    "                    m[i, j] = n\n",
    "                else:\n",
    "                    break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have less yellow!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1 -o\n",
    "m = np.zeros((size, size), dtype=np.int32)\n",
    "mandelbrot_cython(m, size, iterations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_cy = _"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Huge speedup!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"red\"> *Exercise 3.2*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Take the following Python code which returns prime numbers\n",
    "2. Run it for `kmax=100` and time it\n",
    "3. Add the most simple cythonization, using the annotation option\n",
    "4. Run again and time, compare to pure Python\n",
    "5. Cythonize it properly! (Hint: replace `p=[]` also). Use annotation, and compare to simple cython\n",
    "6. Run again and time, compare to pure Python and simple Cython\n",
    "6. Try just-in-time compilation using `Numba` and compare to Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def primes(kmax):  \n",
    "    p = []\n",
    "    result = []  \n",
    "    if kmax > 1000:\n",
    "        kmax = 1000\n",
    "    k = 0\n",
    "    n = 2\n",
    "    while k < kmax:\n",
    "        i = 0\n",
    "        while i < k and n % p[i] != 0:\n",
    "            i = i + 1\n",
    "        if i == k:\n",
    "            p.append(n)\n",
    "            k = k + 1\n",
    "            result.append(n)\n",
    "        n = n + 1\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### <font color=\"green\"> *Solution*\n",
    "\n",
    "> The solution can be found in the solutions.ipynb notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}