Introduction to Cython.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intro to Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why Cython\n",
    "\n",
    "![DevTime](whycython.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Outline:\n",
    "\n",
    "* Speed up Python code\n",
    "* Interact with NumPy arrays\n",
    "* Release GIL and get parallel performance\n",
    "* Wrap C/C++ code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: speed up your Python code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to integrate the function $f(x) = x^4 - 3x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    y = x**4 - 3*x\n",
    "    return y\n",
    "    \n",
    "def integrate_f(a, b, n):\n",
    "    dx = (b - a) / n\n",
    "    dx2 = dx / 2\n",
    "    s = f(a) * dx2\n",
    "    for i in range(1, n):\n",
    "        s += f(a + i * dx) * dx\n",
    "    s += f(b) * dx2\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's time this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit integrate_f(-100, 100, int(1e5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not too bad, but this can add up. So how can we speed this up? Use Cython!\n",
    "\n",
    "**Cython is an optimising static compiler for both the Python programming language and the extended Cython programming language** (based on Pyrex). It makes writing C extensions for Python as easy as Python itself.\n",
    "\n",
    "Cython gives you the combined power of Python and C to let you\n",
    "\n",
    "- write Python code that calls back and forth from and to C or C++ code natively at any point.\n",
    "- easily tune readable Python code into plain C performance by adding static type declarations, also in Python syntax.\n",
    "- use combined source code level debugging to find bugs in your Python, Cython and C code.\n",
    "- interact efficiently with large data sets, e.g. using multi-dimensional NumPy arrays.\n",
    "- quickly build your applications within the large, mature and widely used CPython ecosystem.\n",
    "- integrate natively with existing code and data from legacy, low-level or high-performance libraries and applications.\n",
    "\n",
    "In brief, Cython code is sort of Python code with the best parts of C. Cython code needs to be compiled into C code, which is then compiled as a shared library when your packages is installed. \n",
    "\n",
    "Jupyter \"makes it easy\" to write and compile Cython code through the magic below. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first use our plain-old Python code, and compile it to C without doing anything to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "def f2(x):\n",
    "    y = x**4 - 3*x\n",
    "    return y\n",
    "    \n",
    "def integrate_f2(a, b, n):\n",
    "    dx = (b - a) / n\n",
    "    dx2 = dx / 2\n",
    "    s = f2(a) * dx2\n",
    "    for i in range(1, n):\n",
    "        s += f2(a + i * dx) * dx\n",
    "    s += f2(b) * dx2\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit integrate_f2(-100, 100, int(1e5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a little bit faster, which is nice since all we did was to call Cython on the exact same code. But can we do better?\n",
    "\n",
    "An easy speed up is to annotate the type of objects you are using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "def f3(double x):\n",
    "    y = x**4 - 3*x\n",
    "    return y\n",
    "    \n",
    "def integrate_f3(double a, double b, int n):\n",
    "    dx = (b - a) / n\n",
    "    dx2 = dx / 2\n",
    "    s = f3(a) * dx2\n",
    "    for i in range(1, n):\n",
    "        s += f3(a + i * dx) * dx\n",
    "    s += f3(b) * dx2\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit integrate_f3(-100, 100, int(1e5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final bit of \"easy\" Cython optimization is \"declaring\" the variables inside the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "def f4(double x):\n",
    "    y = x**4 - 3*x\n",
    "    return y\n",
    "    \n",
    "def integrate_f4(double a, double b, int n):\n",
    "    cdef:\n",
    "        double dx = (b - a) / n\n",
    "        double dx2 = dx / 2\n",
    "        double s = f4(a) * dx2\n",
    "        int i = 0\n",
    "    for i in range(1, n):\n",
    "        s += f4(a + i * dx) * dx\n",
    "    s += f4(b) * dx2\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit integrate_f4(-100, 100, int(1e5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's another way to declare types of elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "def f4(double x):\n",
    "    y = x**4 - 3*x\n",
    "    return y\n",
    "    \n",
    "def integrate_f4(double a, double b, int n):\n",
    "    cdef double dx = (b - a) / n\n",
    "    cdef double dx2 = dx / 2\n",
    "    cdef double s = f4(a) * dx2\n",
    "    cdef int i = 0\n",
    "    for i in range(1, n):\n",
    "        s += f4(a + i * dx) * dx\n",
    "    s += f4(b) * dx2\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4X speedup with so little effort is pretty nice. What else can we do?\n",
    "\n",
    "Cython has a nice \"-a\" flag (for annotation) that can provide clues about why your code is slow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "\n",
    "def f4(double x):\n",
    "    y = x**4 - 3*x\n",
    "    return y\n",
    "    \n",
    "def integrate_f4(double a, double b, int n):\n",
    "    cdef:\n",
    "        double dx = (b - a) / n\n",
    "        double dx2 = dx / 2\n",
    "        double s = f4(a) * dx2\n",
    "        int i = 0\n",
    "    for i in range(1, n):\n",
    "        s += f4(a + i * dx) * dx\n",
    "    s += f4(b) * dx2\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a lot of yellow still! How do we reduce this?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: change the `f4` declaration to C\n",
    "\n",
    "Rename the functions to `f5` and `integrate_f5` and finish optimizing!\n",
    "\n",
    "Tip: Python function calls can be expensive – in Cython doubly so because one might need to convert to and from Python objects to do the call. In our example above, the argument is assumed to be a C double both inside f() and in the call to it, yet a Python float object must be constructed around the argument in order to pass it.\n",
    "\n",
    "Cython provides a syntax for declaring a C-style function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit integrate_f5(-100, 100, 10**5) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate_f5(-100, 100, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: work with NumPy arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above is a very small subset of Python. Most scientific application deal not with single values, but with arrays of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def mean3filter(arr):\n",
    "    arr_out = np.empty_like(arr)\n",
    "    for i in range(1, arr.shape[0] - 1):\n",
    "        arr_out[i] = np.sum(arr[i-1:i+2]) / 3\n",
    "    arr_out[0] = (arr[0] + arr[1]) / 2\n",
    "    arr_out[-1] = (arr[-1] + arr[-2]) / 2\n",
    "    return arr_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.random.rand(10**5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit mean3filter(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import cython\n",
    "import numpy as np\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "def mean3filter2(double[::1] arr):\n",
    "    cdef double[::1] arr_out = np.empty_like(arr)\n",
    "    cdef int i\n",
    "    for i in range(1, arr.shape[0]-1):\n",
    "        arr_out[i] = np.sum(arr[i-1 : i+2]) / 3\n",
    "    arr_out[0] = (arr[0] + arr[1]) / 2\n",
    "    arr_out[-1] = (arr[-1] + arr[-2]) / 2\n",
    "    return np.asarray(arr_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit mean3filter2(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rubbish! How do we fix this?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: use `%%cython -a` to speed up the code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Cython in production code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `setup.py` to build your Cython files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "```python\n",
    "from distutils.core import setup\n",
    "from distutils.extension import Extension\n",
    "from Cython.Distutils import build_ext\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "setup(\n",
    "  cmdclass = {'build_ext': build_ext},\n",
    "  ext_modules = [\n",
    "    Extension(\"prange_demo\", [\"prange_demo.pyx\"],\n",
    "              include_dirs=[np.get_include()],\n",
    "              extra_compile_args=['-fopenmp'],\n",
    "              extra_link_args=['-fopenmp', '-lgomp']),\n",
    "  ]\n",
    ")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Write a Cython module with a setup.py to run the mean-3 filter, then import from the notebook.\n",
    "- Write a mean3.pyx\n",
    "- Write a corresponding setup.py\n",
    "- Run python setup.py build_ext -i\n",
    "\n",
    "NOTE: Don't add extra compiler or link args!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'mean3'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-0f875b93d1ce>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mmean3\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmean3filter3\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mrin\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mrout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mmean3filter3\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrin\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'mean3'"
     ]
    }
   ],
   "source": [
    "from mean3 import mean3filter3\n",
    "rin = np.random.rand(10)\n",
    "rout = np.zeros(10)\n",
    "mean3filter3(rin, rout)\n",
    "print(rout) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean3filter3?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr_out = np.zeros_like(r)\n",
    "%timeit mean3filter3(r, arr_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3: write parallel code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Warning:**: Dragons afoot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "import cython\n",
    "from cython.parallel import prange\n",
    "import numpy as np\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "def mean3filter4(double[::1] arr, double[::1] out):\n",
    "    cdef int i, j, k = arr.shape[0]-1\n",
    "    with nogil:\n",
    "        for i in prange(1, k-1, schedule='static',\n",
    "                        chunksize=(k-2) // 2, num_threads=2):\n",
    "            for j in range(i-1, i+1):\n",
    "                out[i] += arr[j]\n",
    "            out[i] /= 3\n",
    "    out[0] = (arr[0] + arr[1]) / 2\n",
    "    out[-1] = (arr[-1] + arr[-2]) / 2\n",
    "    return np.asarray(out) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rin = np.random.rand(10**8)\n",
    "rout = np.empty_like(rin)\n",
    "\n",
    "%timeit mean3filter4(rin, rout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# comparing to single-threaded...\n",
    "rin = np.random.rand(10**8)\n",
    "rout = np.empty_like(rin)\n",
    "\n",
    "%timeit mean3filter3(rin, rout)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 4: interact with C/C++ code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "# distutils: language=c++\n",
    "import cython\n",
    "from libcpp.vector cimport vector\n",
    "\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "def build_list_with_vector(double[::1] in_arr):\n",
    "    cdef vector[double] out\n",
    "    cdef int i\n",
    "    for i in range(in_arr.shape[0]):\n",
    "        out.push_back(in_arr[i])\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "build_list_with_vector(np.random.rand(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: C++ int graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "#distutils: language=c++\n",
    "from cython.operator cimport dereference as deref, preincrement as inc\n",
    "\n",
    "from libcpp.vector cimport vector\n",
    "from libcpp.map cimport map as cppmap\n",
    "\n",
    "cdef class Graph:\n",
    "    cdef cppmap[int, vector[int]] _adj\n",
    "    \n",
    "    cpdef int has_node(self, int node):\n",
    "        return self._adj.find(node) != self._adj.end()\n",
    "    \n",
    "    cdef void add_node(self, int new_node):\n",
    "        cdef vector[int] out\n",
    "        if not self.has_node(new_node):\n",
    "            self._adj[new_node] = out\n",
    "            \n",
    "    def add_edge(self, int u, int v):\n",
    "        self.add_node(u)\n",
    "        self.add_node(v)\n",
    "        self._adj[u].push_back(v)\n",
    "        self._adj[v].push_back(u)\n",
    "        \n",
    "    def __getitem__(self, int u):\n",
    "        return self._adj[u]\n",
    "    \n",
    "    cdef vector[int] _degrees(self):\n",
    "        cdef vector[int] deg\n",
    "        cdef int first = 0\n",
    "        cdef vector[int] edges\n",
    "        cdef cppmap[int, vector[int]].iterator it = self._adj.begin()\n",
    "        while it != self._adj.end():\n",
    "            deg.push_back(deref(it).second.size())\n",
    "            it = inc(it)\n",
    "        return deg\n",
    "            \n",
    "    def degrees(self):\n",
    "        return self._degrees()\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g0 = Graph() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g0.add_edge(1, 5)\n",
    "g0.add_edge(1, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g0[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g0.has_node(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g0.degrees()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "g = nx.barabasi_albert_graph(100000, 6)\n",
    "with open('graph.txt', 'w') as fout:\n",
    "    for u, v in g.edges_iter():\n",
    "        fout.write('%i,%i\\n' % (u, v))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit list(g.degree())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myg = Graph()\n",
    "def line2edges(line):\n",
    "    u, v = map(int, line.rstrip().split(','))\n",
    "    return u, v\n",
    "\n",
    "edges = map(line2edges, open('graph.txt'))\n",
    "\n",
    "for u, v in edges:\n",
    "    myg.add_edge(u, v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit mydeg = myg.degrees()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numba and interacting with SciPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def addarr(x, y):\n",
    "    result = np.zeros_like(x)\n",
    "    for i in range(x.size):\n",
    "        result[i] = x[i] + y[i]\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How long does this take in pure Python?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = int(1e6)\n",
    "a = np.random.rand(n)\n",
    "b = np.random.rand(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit -r 1 -n 1 addarr(a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "About half a second on my machine. Let's try with Numba using its JIT decorator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numba\n",
    "\n",
    "addarr_nb = numba.jit(addarr) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit -r 1 -n 1 addarr_nb(a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first time it runs, it's only a tiny bit faster. That's because of the nature of JITs: they only compile code *as it is being run*, in order to use object type information of the objects passed into the function. (Note that, in Python, the arguments `a` and `b` to `addarr` could be anything: an array, as expected, but also a list, a tuple, even a `Banana`, if you've defined such a class, and the meaning of the function body is different for each of those types.)\n",
    "\n",
    "Let's see what happens the next time we run it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit -r 1 -n 1 addarr_nb(a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whoa! Now the code takes 5ms, about 100 times faster than the pure Python version. And the NumPy equivalent?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit -r 1 -n 1 a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only marginally faster than Numba, even though NumPy addition is implemented in highly optimised C code. \n",
    "\n",
    "Be careful! Both Numba and numpy are actively developped, and some operations may be faster in specific version of Numba versus specific version of numpy, but not others."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Write a mean filter using Numba, and time it, using `-r 1 -n 1` as above. How does it compare to Cython?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# decorate here...\n",
    "def mean3filter_nb(arr, arr_out):\n",
    "    pass  # ... and fill this in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rin = np.random.rand(10**8)\n",
    "rout = np.zeros_like(rin)\n",
    "\n",
    "%timeit -r 1 -n 1 mean3filter_nb(rin, rout)  # jit warmup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rin = np.random.rand(10**8)\n",
    "rout = np.zeros_like(rin)\n",
    "\n",
    "%timeit -r 1 -n 1 mean3filter_nb(rin, rout)  # faster!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SciPy LowLevelCallable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See:\n",
    "- https://ilovesymposia.com/2017/03/12/scipys-new-lowlevelcallable-is-a-game-changer/\n",
    "- https://ilovesymposia.com/2017/03/15/prettier-lowlevelcallables-with-numba-jit-and-decorators/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numba\n",
    "from numba import cfunc, carray\n",
    "from numba.types import intc, CPointer, float64, intp, voidptr\n",
    "from scipy import LowLevelCallable\n",
    "\n",
    "def jit_filter_function(filter_function):\n",
    "    jitted_function = numba.jit(filter_function, nopython=True)\n",
    "    \n",
    "    @cfunc(intc(CPointer(float64), intp, CPointer(float64), voidptr))\n",
    "    def wrapped(values_ptr, len_values, result, data):\n",
    "        values = carray(values_ptr, (len_values,), dtype=float64)\n",
    "        result[0] = jitted_function(values)\n",
    "        return 1\n",
    "    \n",
    "    return LowLevelCallable(wrapped.ctypes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit_filter_function\n",
    "def fmin(values):\n",
    "    result = np.inf\n",
    "    for v in values:\n",
    "        if v < result:\n",
    "            result = v\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = np.random.random((2048, 2048))\n",
    "footprint = np.ones((3, 3))\n",
    "%timeit ndi.generic_filter(image, fmin, footprint=footprint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cython can also be used to create LowLevelCallables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython --name=test9\n",
    "\n",
    "from libc.stdint cimport intptr_t\n",
    "from numpy.math cimport INFINITY\n",
    "\n",
    "cdef api int erosion_kernel(double* input_arr_1d, intptr_t filter_size,\n",
    "                            double* return_value, void* user_data):\n",
    "    cdef:\n",
    "        double[:] input_arr\n",
    "        ssize_t i\n",
    "        double result = INFINITY\n",
    "    for i in range(filter_size):\n",
    "        if input_arr_1d[i] < return_value[0]:\n",
    "            result = input_arr_1d[i]\n",
    "    return_value[0] = result\n",
    "    return 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import test9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import LowLevelCallable, ndimage\n",
    "import sys\n",
    "\n",
    "def erosion_fast(image, footprint):\n",
    "    out = ndimage.generic_filter(\n",
    "            image,\n",
    "            LowLevelCallable.from_cython(test9, name='erosion_kernel'),\n",
    "            footprint=footprint\n",
    "    )\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit erosion_fast(image, footprint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Write a mean filter using `generic_filter`, and time it. How does it compare to the above?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write a filter function here...\n",
    "\n",
    "def mean3filter_ndi(arr_in, arr_out):\n",
    "    pass  # ... then fill this in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rin = np.random.rand(10**8)\n",
    "rout = np.zeros_like(rin)\n",
    "\n",
    "%timeit -r 1 -n 1 mean3filter_ndi(rin, rout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rin = np.random.rand(10**8)\n",
    "rout = np.zeros_like(rin)\n",
    "\n",
    "%timeit -r 1 -n 1 mean3filter3(rin, rout)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is about 2-3x slower than Cython and 6x slower than Numba, probably accounting for some overhead of allowing generic filters. (But to be honest I don't know why this is so much slower!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Concluding remarks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rout.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some pros and cons about Cython and Numba:\n",
    "\n",
    "- Cython pros:\n",
    "  * very wide support\n",
    "  * easy to distribute compiled code to most users\n",
    "  * quite developed optimizing workflow (i.e. cython -a)\n",
    "- Cython cons:\n",
    "  * need to use a new language\n",
    "\n",
    "\n",
    "- Numba pros:\n",
    "  * quite easy to use, especially if you're starting from Cython code\n",
    "  * often eye-popping, face-melting performance\n",
    "- Numba cons:\n",
    "  * hard to install outside of conda\n",
    "  * hard to optimise. If it's slow, you have to guess (though they are helpful on mailing list)\n",
    "  * many parts of Python still unsupported, e.g. dicts.\n",
    "  * project still young and some people are paranoid that it could disappear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}