numba_comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Playing with Numba and finite differences\n",
    "**David I. Ketcheson**\n",
    "\n",
    "**May 25, 2015**\n",
    "\n",
    "In this notebook, I test the speed of naive Python, vectorized\n",
    "numpy, numba, Cython, and Fortran for solving the heat equation\n",
    "on a regular Cartesian grid in 1- and 2-D with centered differences\n",
    "in space and the explicit Euler method in time.\n",
    "You can view this as an update to \n",
    "[Jake Vanderplas' blog post](https://jakevdp.github.io/blog/2013/06/15/numba-vs-cython-take-2/)\n",
    "from two years ago, but with a different kernel.\n",
    "\n",
    "Note that this operation is just a sparse matrix-vector multiply, but\n",
    "I'm not interested in implementing it that way because my real algorithms\n",
    "of interest are solvers for nonlinear hyperbolic PDEs.  There are \n",
    "better algorithms for this problem involving e.g. FFTs or implicit methods,\n",
    "but for the same reason I'm not primarlily interested in those approaches.\n",
    "\n",
    "I've run this on two machines and gotten some significantly different\n",
    "results, so YMMV.  The results below were computed on a Mac Pro with\n",
    "2 quad-core 2.66 GHz Intel Xeon processors and 16GB of RAM.\n",
    "\n",
    "**tl;dr: numba is competitive with Fortran, especially on grids of\n",
    "1000x1000 or larger; everything else is significantly slower.  My Cython\n",
    "skills may be lacking.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.18.2'"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "matplotlib.rcParams.update({'font.size': 15})\n",
    "import numpy as np\n",
    "from numba import jit\n",
    "import numba\n",
    "numba.__version__\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1D heat equation\n",
    "\n",
    "## Naive\n",
    "A naive implementation of the second-order centered-difference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 loops, best of 3: 117 µs per loop\n"
     ]
    }
   ],
   "source": [
    "def centered_difference_1D(u,dx=1):\n",
    "    D = np.empty_like(u)\n",
    "    for i in range(1,len(D)-1):\n",
    "        D[i] = (u[i+1] - 2*u[i] + u[i-1])/dx**2\n",
    "    return D\n",
    "\n",
    "N = 100\n",
    "u = np.random.rand(N)\n",
    "u[0] = 0; u[-1] = 0\n",
    "\n",
    "times = {}\n",
    "y = %timeit -o centered_difference_1D(u)\n",
    "times['Python'] = y.best\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numba jit\n",
    "Obviously, it's slow.  How much can we speed it up with numba?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def centered_difference_1D_numba(u,dx=1):\n",
    "    D = np.empty_like(u)\n",
    "    for i in range(1,len(D)-1):\n",
    "        D[i] = (u[i+1] - 2*u[i] + u[i-1])/dx**2\n",
    "    return D\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 107979.69 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "1000000 loops, best of 3: 1.33 µs per loop\n"
     ]
    }
   ],
   "source": [
    "y = %timeit -o centered_difference_1D_numba(u)\n",
    "times['numba'] = y.best\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not bad: a 100x speedup.\n",
    "\n",
    "## Time to solve the 1D Heat equation\n",
    "Now let's actually solve the heat equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python:\n",
      "1 loops, best of 3: 3.18 s per loop\n",
      "Numba:\n",
      "1 loops, best of 3: 116 ms per loop\n"
     ]
    }
   ],
   "source": [
    "def solve_heat_equation_1D(u,dx,dt,T):\n",
    "    \"\"\"Advance solution u by time T using the explicit Euler method.\"\"\"\n",
    "    n_steps = int(T/dt)\n",
    "    for j in range(n_steps):\n",
    "        u = u + dt*centered_difference_1D(u,dx)\n",
    "    return u\n",
    "\n",
    "@jit  # jit only makes a very small improvement here\n",
    "def solve_heat_equation_1D_numba(u,dx,dt,T):\n",
    "    \"\"\"Advance solution u by time T using the explicit Euler method.\"\"\"\n",
    "    n_steps = int(T/dt)\n",
    "    for j in range(n_steps):\n",
    "        u = u + dt*centered_difference_1D_numba(u,dx)\n",
    "    return u\n",
    "\n",
    "x = np.linspace(0,1,N)\n",
    "u = np.exp(-100*(x-0.5)**2)\n",
    "dx = x[1]-x[0]\n",
    "dt = 0.4 * dx**2\n",
    "\n",
    "times = {}\n",
    "print 'Python:'\n",
    "y = %timeit -o solve_heat_equation_1D(u,dx,dt,T=1)\n",
    "times['Python'] = y.best\n",
    "\n",
    "print 'Numba:'\n",
    "y = %timeit -o solve_heat_equation_1D_numba(u,dx,dt,T=1)\n",
    "times['numba'] = y.best\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2D centered difference kernel\n",
    "From here on, I won't time the basic Python implementation because\n",
    "it's so slow.  But I show the code for reference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def centered_difference_2D(u,dx=1.):\n",
    "    D = np.zeros_like(u)\n",
    "    for i in range(1,u.shape[0]-1):\n",
    "        for j in range(1,u.shape[1]-1):\n",
    "            D[i,j] = (u[i+1,j] + u[i,j+1] + u[i-1,j] + u[i,j-1] - 4*u[i,j])/dx**2\n",
    "    return D\n",
    "\n",
    "@jit\n",
    "def centered_difference_2D_numba(u,dx=1.):\n",
    "    D = np.zeros_like(u)\n",
    "    for i in range(1,u.shape[0]-1):\n",
    "        for j in range(1,u.shape[1]-1):\n",
    "            D[i,j] = (u[i+1,j] + u[i,j+1] + u[i-1,j] + u[i,j-1] - 4*u[i,j])/dx**2\n",
    "    return D\n",
    "\n",
    "def centered_difference_2D_vectorized(u,dx=1.):\n",
    "    D = np.zeros_like(u)\n",
    "    D[1:-1,1:-1] = (u[2:,1:-1]+u[1:-1,2:]+u[:-2,1:-1]+u[1:-1,:-2]-4*u[1:-1,1:-1])/dx**2\n",
    "    return D\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython\n",
    "I essentially borrowed this code from Jake Vanderplas' blog."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Cython extension is already loaded. To reload it, use:\n",
      "  %reload_ext Cython\n"
     ]
    }
   ],
   "source": [
    "%load_ext Cython\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "cimport cython\n",
    "import numpy as np\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.wraparound(False)\n",
    "def centered_difference_2D_cython(double[:, :] u, double dx):\n",
    "    cdef int M = u.shape[0]\n",
    "    cdef double[:, :] D = np.empty((M, M), dtype=np.float64)\n",
    "    for i in range(1,M-1):\n",
    "        for j in range(1,M-1):\n",
    "            D[i,j] = (u[i+1,j] + u[i,j+1] + u[i-1,j] + u[i,j-1] - 4*u[i,j])/dx**2\n",
    "    return np.asarray(D)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fortran\n",
    "I can call Fortran by using f2py.  Note that my compiler is gfortran,\n",
    "which is certainly not the best.  Also, I'm unsure of what optimization\n",
    "flags are being thrown by f2py."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting cdiff_fort.f90\n"
     ]
    }
   ],
   "source": [
    "%%file cdiff_fort.f90\n",
    "\n",
    "      subroutine centered_diff_fort(u,D,dx,m)\n",
    "          integer :: m\n",
    "          double precision, intent(in) :: u(m,m), dx\n",
    "          double precision, intent(out) :: D(m,m)\n",
    "          integer :: i,j\n",
    "          do j = 2,m-1 \n",
    "              do i = 2,m-1 \n",
    "                  D(i,j) = (u(i+1,j) + u(i,j+1) + u(i-1,j) + u(i,j-1) \n",
    "     &                   - 4*u(i,j))/dx**2\n",
    "              end do \n",
    "          end do \n",
    "      end subroutine centered_diff_fort\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compile the Fortran with f2py.\n",
    "# We'll direct the output into /dev/null so it doesn't fill the screen\n",
    "!f2py -c cdiff_fort.f90 -m centered_diff_fort > /dev/null\n",
    "# For some reason this fails in the notebook but works fine from a terminal\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time to compute the 2D centered difference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy vectorized\n",
      "10 loops, best of 3: 32 ms per loop\n",
      "numba\n",
      "The slowest run took 33.17 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "1 loops, best of 3: 6.14 ms per loop\n",
      "Cython\n",
      "10 loops, best of 3: 99.6 ms per loop\n",
      "Fortran\n",
      "100 loops, best of 3: 5.02 ms per loop\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGxBJREFUeJzt3XuUXWWd5vHvQ4JIGgmFCw2t3OSWwEyPo9BjREgWMSK2\njQMCGbR1IeCASwZvYIsOpLhoL8CQGQbp1iFcbGwM2IgiA5FEiiZBY7q1u6FJWi5BNG2kNQkQooDm\nN3+87yGbnVNnn7qlKvU+n7WyKue97PPuXVX7Ofvdl1JEYGZm5dlhtAdgZmajwwFgZlYoB4CZWaEc\nAGZmhXIAmJkVygFgZlaoxgCQdIikJZKek7RG0kWSOvaTtKOkKyTdL+k3kjZ3aPseSQ/mdv8i6eTB\nrIiZmQ1M0468B1gM/B44DrgY+BRwUcNy/wA4HdgILAPa3mwg6W3AN4AlwDuBO4GbJc3ufhXMzGww\n1OlGMEnnA+cC+0TExlx2HtALTImIZxvfQDobuCoitgobSYuACRHx9krZncCuEXHkANfFzMwGoGkK\n6FhgUWvnny0EdgZmDOWNJe0EzARuqVUtBKZLetVQlm9mZp01BcDBwKpqQUQ8CWzKdUOxP7BjffnA\nyjyug4a4fDMz66ApAHqADW3K1+e6oWj1ry9/fa3ezMxGgC8DNTMr1MSG+vXA5DblPWz5pD5Yrf71\n5ffU6gGQ5MeWmpkNQkSoXXnTEcAqYFq1QNJewCS2nrsfqMeAF+vLB6YCm4GfDHH5ZmbWQdMRwF3A\neZJ2qVwJNId0Evi+obxxRDwv6V7gJOArlao5wAP9XWLaX5KBjxLqOm2r7YWk3ojoHe1xjBfensNr\nrG/Ppn1i0xHAXwHPA7dJmiXpvwNzgSurl4ZKelTStbU3PlbSicAb8+v3SjpR0t6VZpcAMyXNlzRT\n0uWkS08v7nYFtxaj/G/uGBiDc9DMmnU8AoiIDZJmAVcDd5Dm5a8k3QhWNYGtw+QaYJ/WooBb89cP\nAV/Ny1+WQ+JS4CPA48ApEbF4kOtjZmZd6ngn8FjSOpRpngIa7fXpI93fNto0XqaAZkZE32iPY7zw\n9hxeY317Nu03HQDj1vgIADMbvKb9pu8DMDMrlAPAzKxQDgAzs0I5AMzMCuUAMDMrlAPAzKxQDgAz\ns0I5AMzMCuUAMDMrlAPAzKxQDgAzs0I5AMzMCuUAMDMrlAPAzKxQDgAzs0I5AMzMCuUAMDMrlAPA\nzKxQDgAzs0I5AMzMCuUAMDMrlAPAzKxQDgAzs0I5AMzMCuUAMDMrlAPAzKxQDgAzs0I5AMzMCuUA\nMDMrlAPAzKxQDgAzs0I1BoCkQyQtkfScpDWSLpLUTb/Jkq6XtE7SBkk3Sdq91manvLxHJW2S9Iik\nXkmvGMpKmZlZs4mdKiX1AIuBh4DjgAOAeaTguKBh2bfk9qcDAVwG3A4cVWlzRa7/HPBj4M3ApcBu\nwMcHtipmZjYQioj+K6XzgXOBfSJiYy47D+gFpkTEs/30mw4sA46KiKW57HBgOTA7IpbksqeAGyPi\nvErfecD7I2JKbZkBEBHqMN5IWWOgjtvKzMa/pv1m01TOscCi1s4/WwjsDMxo6Le2tfPPA1gBrM51\nVc/UXj/dMCYzMxsGTQFwMLCqWhARTwKbcl1/ptb7ZStzXcs1wJmS3ippF0lHAmcBVzcN3MzMhqbj\nOQCgB9jQpnx9rhtovw3Afq0XEdEr6TXA0kqbL0XEpQ3jMjOzIWoKgJHw0iS9pM8Dc4CzgX8G3ghc\nImldRMwdhbGZmRWjKQDWA5PblPfkuv6sA/bo1E/Sa4FPA2dFxIJcv1TSC8DVkv5PRPyqvgBJvZWX\nfRHR17AOZmbFkDQTmNlN26YAWAVMqy18L2AS7ef4q/2ObFM+Fbgt/38fYALwT7U2/5jHtQ+wVQBE\nRG/DmM3MipU/FPcBSOo4k9J0Evgu4BhJu1TK5pBOAt/X0G+KpCNaBZIOI83/35WLnshf31Tr++Za\nvZmZjYCm+wB2Ax4m3Qh2GbA/6Uaw+RFxYaXdo6TpmDMqZXcDB5LuI2jdCLY2ImZU2twCzAbmAg+S\nzgHMBe6OiP9WG4vvAxgQ3wdgVrqm/WbHAMgLmEa6LHM6af7+WqA3Kh0lrQbujYjTKmWTgfnA8aQj\njTuAcyJiXaXNLqQ7ik8AXgf8nDRFdElEPDeQFdnSxgGQOADMSjfkABgrHAAD5QAwK91Q7wQ2M7Nx\nygFgZlYoB4CZWaEcAGZmhXIAmJkVygFgZlYoB4CZWaEcAGZmhXIAmJkVygFgZlYoB4CZWaEcAGZm\nhXIAmJkVygFgZlYoB4CZWaEcAGZmhXIAmJkVygFgZlYoB4CZWaEcAGZmhXIAmJkVygFgZlYoB4CZ\nWaEcAGZmhXIAmJkVygFgZlYoB4CZWaEcAGZmhXIAmJkVygFgZlYoB4CZWaEaA0DSIZKWSHpO0hpJ\nF0nqpt9kSddLWidpg6SbJO3ept2rJX1Z0i8kbZK0UtIHBrtCZmbWnYmdKiX1AIuBh4DjgAOAeaTg\nuKBh2bfk9qcDAVwG3A4cVVn+rsDfAc8AZwO/Ag4Fdhz4qpiZ2UB0DADgLGAn4ISI2AgsyTvtXkmX\nR8Sz7TpJmg7MBo6KiKW5bA2wXNKsiFiSm36WtLOfGRHP57L7hrZKZmbWjaapnGOBRXnn37IQ2BmY\n0dBvbWvnDxARK4DVua7lQ8CCys7fzMy2kaYAOBhYVS2IiCeBTbmuP1Pr/bKVuQ5J+wF7AE9L+n+S\nnpf0lKR5kjwFZGY2wpoCoAfY0KZ8fa4baL8NlX5T8tfLgZ8BxwBfAD4CXNowLjMzG6KmcwAjIfJX\n5a8PRcSZ+f99kl4FfFbShZ4aMjMbOU1HAOuByW3Ke3Jdf9YBuzX0a329t9bmXtKJ5wMaxmZmZkPQ\ndASwCphWLZC0FzCJ9nP81X5HtimfCtyW//8Y8AJbjgReeov8NWhDUm/lZV9E9HUYh5lZUSTNBGZ2\n07YpAO4CzpO0S+VKoDmkk8CdLte8C7hA0hERsSwP6jBgv1xHRLwg6R7g6FrfWcBzwCPtFhwRvQ1j\nNjMrVv5Q3AcgaW6ntopo+0Gb3Hk34GHSjWCXAfuTbgSbHxEXVto9Svo0fkal7G7gQOBcttwItjYi\nZlTaHA4sBb4GfB34I+AS4OKI+IvaWCKvXP2Iodam//UpizpuKzMb/5r2mx0DIC9gGnA1MJ00b38t\n0BuVjpJWA/dGxGmVssnAfOB40rmGO4BzImJdbfnvAP6CdAfwL4GvAF+I2sAcAAPlADAr3ZADYKxw\nAAyUA8CsdE37TT8N1MysUA4AM7NCOQDMzArlADAzK5QDwMysUA4AM7NCOQDMzArlADAzK5QDwMys\nUA4AM7NCOQDMzArlADAzK5QDwMysUA4AM7NCOQDMzArlADAzK5QDwMysUA4AM7NCOQDMzArlADAz\nK5QDwMysUA4AM7NCOQDMzArlADAzK5QDwMysUA4AM7NCOQDMzArlADAzK5QDwMysUA4AM7NCOQDM\nzArlADAzK1RjAEg6RNISSc9JWiPpIknd9Jss6XpJ6yRtkHSTpN07tH+PpM2SVgx0JczMbOAmdqqU\n1AMsBh4CjgMOAOaRguOChmXfktufDgRwGXA7cFSb93klMB/4ZW5rZmYjrGMAAGcBOwEnRMRGYImk\nXYFeSZdHxLPtOkmaDswGjoqIpblsDbBc0qyIWFLrch7wc+Ax4D8MfnXMzKxbTVM5xwKL8s6/ZSGw\nMzCjod/a1s4fICJWAKtz3Usk7U0KgHMAdT90MzMbiqYAOBhYVS2IiCeBTbmuP1Pr/bKVua5qHrAw\nIv6xYSxmZjaMmqaAeoANbcrX57qB9tsA7Nd6Ielo0lTRgQ3jMDOzYTYal4EGgKSJwFXApRHx7/V6\nMzMbWU1HAOuByW3Ke3Jdf9YBezT0+zCwK3CjpN1y2SuAiZImA89FxO8axmdmZoPUFACrgGnVAkl7\nAZNoP8df7Xdkm/KpwG35/wcBrydd+lm3Hvgz4G/qFZJ6Ky/7IqKvwzjMzIoiaSYws6u2Ef3PuEj6\nDOkKnX1aVwJJOhfoBabUrg6q9nsL8ABwZEQsy2WHAT8E3h4R35O0P/C6ajfgM8C+wJnAqoh4qrLM\nAIiIfq8USm08g5So47Yys/Gvab/ZFAC7AQ+TbgS7DNifdNXO/Ii4sNLuUdKn8TMqZXeTTu6ey5Yb\nwdZGRL+Xj0q6ATg0Ig4f6IpsaeMASBwAZqVr2m92PAkcERuAWcAE4A5gLnBl/lo1oc2y5gD3AdcB\nNwIrgOMbxht4D25mtk10PAIYS3wEMFA+AjAr3ZCOAMzMbPxyAJiZFcoBYGZWKAeAmVmhHABmZoVy\nAJiZFcoBYGZWKAeAmVmhHABmZoVyAJiZFcoBYGZWKAeAmVmhmv4gjJkNg9ZDuSzxgwrHBgeA2Tbj\nDEi87x8rPAVkZlYoB4CZWaEcAGZmhXIAmJkVygFgZlYoB4CZWaEcAGZmhXIAmJkVygFgZlYoB4CZ\nWaEcAGZmhXIAmJkVygFgZlYoB4CZWaEcAGZmhXIAmJkVygFgZlYoB4CZWaG6CgBJh0haIuk5SWsk\nXSSpsa+kyZKul7RO0gZJN0navVK/g6TPSHogt/mVpEWSDhvKSpmZWbNuduI9wGLg98BxwMXAp4CL\nulj+LcBRwOnAqcDhwO2V+knAp4EHgPcBfwa8CCyV9KZuV8LMzAZOEZ3/ULWk84FzgX0iYmMuOw/o\nBaZExLP99JsOLAOOioiluexwYDkwOyKW5KOIV0XE05V+OwI/Ae6NiNMq5QEQEf3+RenUxn94O1HH\nbWXbln82q/yzua007Te7mQI6FljU2vlnC4GdgRkN/da2dv55ECuA1bmOiNhc3fnnsheBh4E9uxib\nmZkNUjcBcDCwqloQEU8Cm3Jdf6bW+2Urc11bknYC3kQ6CjAzsxHSTQD0ABvalK/PdQPtt6Gh3+eA\n3YCruxibmZkN0sRRet+2k6GS/gT4LPDJiHhk2w7JzKws3QTAemBym/KeXNefdcAe3fbLJ4gXAn8Z\nEVf1t1BJvZWXfRHR12EMZmZFkTQTmNlN224CYBUwrfYGe5Eu4Ww3x1/td2Sb8qnAbbXlHQTcCdwD\nnNNpMBHR2zhiM7NC5Q/FfQCS5nZq2805gLuAYyTtUimbQzoJfF9DvymSjmgV5Bu89st1rbI9gUXA\nI8Ap0XRdqpmZDYtu7gPYjXRZ5kPAZcD+wDxgfkRcWGn3KGlK5oxK2d3AgaT7CCL3XxsRM3L9zsD3\ngX2A95OmjVqej4gfV5bl+wAGxNdajyX+2azyz+a20rTfbJwCiogNkmaRrsq5gzR/fyXpRrCqCWx9\nRDEHmA9cl+vu4OVTPK8F/oj0m/GdWt8ngDc0jc/MzAan8QhgrPARwED5U9ZY4p/NKv9sbivDcSew\nmZmNQw4AM7NCOQDMzArlADAzK5QDwMysUA4AM7NCOQDMzArlADAzK5QDwMysUA4AM7NCOQDMzArl\nADAzK5QDwMysUA4AM7NCOQDMzArlADAzK5QDwMysUA4AM7NCOQDMzArlADAzK5QDwMysUA4AM7NC\nOQDMzArlADAzK5QDwMysUBNHewA2dkmK0R7DWBIRGu0xmA0nB4A1cAYk3vfb+OMpIDOzQvkIwMy2\nO56efLnBTk86AMxsO+UMSAY/PekpIDOzQjUGgKRDJC2R9JykNZIuktRNv8mSrpe0TtIGSTdJ2r1N\nu/dIelDSbyT9i6STB7syZmbWvY47ckk9wGLg98BxwMXAp4CLulj2LcBRwOnAqcDhwO215b8N+Aaw\nBHgncCdws6TZA1kJMzMbOEX0P48m6XzgXGCfiNiYy84DeoEpEfFsP/2mA8uAoyJiaS47HFgOzI6I\nJblsETAhIt5e6XsnsGtEHFlbZkDnkx2pjecFEw35unVvz6qhbU9vyyr/bA6v/rdn036zaSrnWGBR\na+efLQR2BmY09Fvb2vnnAawAVuc6JO0EzCQdKVQtBKZLelXD2MzMbAiaAuBgYFW1ICKeBDbluv5M\nrffLVuY6gP2BHdu0W5nHdVDD2MaovtEewDjTN9oDGGf6RnsA40zfaA9gSJoCoAfY0KZ8fa4baL8N\nlX49lbL6sqv125m+0R7AONM32gMYZ/pGewDjTN9oD2BIRuMyUE/cmZmNAU03gq0HJrcp72HLJ/V2\n1gF7NPRrfa0vv6dW/zLNdwCOhWe2dHOR1MgbnrslvT1bhr49x8K2hLGwPcfPzyZsz9uz6QhgFTCt\n9kZ7AZNoP8df7Te1TXn13MBjwIv15ec2m4GfNIzNzMyGoOkI4C7gPEm7VK4EmkM6CXxfQ78LJB0R\nEcsAJB0G7JfriIjnJd0LnAR8pdJ3DvBA/RJTP4rXzGx4Nd0HsBvwMPAQcBnpyp15wPyIuLDS7lGg\nLyLOqJTdDRxIuo8gcv+1ETGj0uYI0lmUq4FvAe8i3Wh2TEQsHp5VNDOzdjpOAUXEBmAWMAG4A5gL\nXJm/Vk1os6w5pKOE64AbgRXA8bXlLwNOBN4O3A28GzjFO38zs5HX8QhgPJP0XuCjwH8m3dj2U+A7\nwBcj4hdd9D8IeB/paOjpSvmppNDbJSI2jcDQrULSE8CtEXHeaI/Fxg5JvcCFbaoWR8Q7hrDcHYHP\nAd+MiH8a7HLGiiIfBy1pHvAx0o56HvAMcChwFuk8xQldLOYg0g/YdcDTDW1t5AS+tNjaexo4pk3Z\nUOxE+r1/HHAAbG8k/SnwCeC0iLihUnW/pK8AA30QnU9Om41Nv4uIHw7XwiS9svqyi/Y7R8Rvhuv9\nR0KJfw/gE8A/1Hb+AETE5ohYJOmHkq6v10u6QdKPJM0Avp2LV0vaLOnxWvM3SLpH0kZJKyUdX6tH\n0tmSHpH02/z147X6Xkn/LumNkn6QH8n9o/wU1VGXt8cKSbMl/XNe1/slHZLr983b5l3t+lVet9bz\njyX9vaRNeTn7StpT0rclPZsfFz6z/VB0gaS1ud1NknatVE6SdLWkVXkbPp5fj6nnTY2F7SnpCUlX\n9Lc9JU2Q9G+S6ucBkdQn6bYR2TjDTNLRkpYrPYZ+raQvSfqDSv3MvK3f0dpepItVnslNrs/1myXt\nXfnevE/SVyWtJ13YgqQPSloq6ddKj8f/nqQ318bT8Xs/UooKAKX5u+mkE86dXAucWPuB2AV4L7AA\n+BHp6iZIJ7bfQu0EN/A3pMdf/1fgEeDrkl5XWd6Hgatym3cDtwLzJP15bTmTSCfR/zK///PAbZJ2\n7mKVR1oAewOXA5cApwCvIT3Qr5u+VZNIlwPPy8vZG7iJ9LDAPtL2/TfgG7V1V25/NOnR458E/oT0\nPawueyJwAemx4xfk9rd2tZbbzljYnkGH7RkRvwduAD5YfTNJbwCOJP1+jBk5sCa2/uWyQ0n7gKdI\n071zSefzvtFmEQuAHwN/StoGR+fyS0i/928B1lbaf5E0zXQi8IVcti/w16RL3k8Bfkaacdiv0m8o\n3/vBi4hi/gFTSDeZfbih3a7ARuDUStlpwG+Bnvz63XlZe9f6nprLq313J930dmZ+vQOwBlhQ6/sl\n0rORXpFf9+Zlzay0+U+57B1jYHvekNdr/0rZe/L4DiL94G8G3tWm34rK69Z6Hlkp+0gu+5+Vsmm5\n7J2VsieAXwGTKmXvI/0Ni6n9jHsicERe1utHeztub9sTOKDNz+XFpEDZYbS3Y20b1P/NAr4O/Cv5\nIpjc/qRc/5b8emZ+Pa+23F1y+Qdr5a3vzd82jGuH/PO3Erig2+/9SG2noo4AKjqeNIyIZ0ifBk6t\nFJ8KfCsiOj0Co+q7leWtI33aaB0BvB7Yk60/gd5CCp//WCl7ISL6Kq9XVpYxFqyOiMcqrwc7vhci\n4v7K69Yyv9em7A9rfe+Jl19xdTvpyOCwVoGkD0j6cT6UfwFovddYe+rsmN+eEfEo8Hfk3w9JIh0R\n/HVEbB7gOEfS06QxV/8tB/6YdBVPdT9wG/A70geDqjsH+J5btZc0TdI3Ja3N7/EC6WnKB9aaDtf3\nvmulnQT+NWkKZe8u2i4A+iTtS7rP4W3kv2XQpfpTTl8AWieR9sxff1lr03pd/dOZ9TuiX0i/b1RP\nSI2mdusJAx9f/Y8LtZbz0vL7WfcghSuVdpskbSRvZ6XzLzcC1wCfIT2r6g+Bbw5inCNtzG/PbAFw\njaSPkqZB9iZdETeW/C4iflQvlDSF2u9eRPxe0q95+e8e9XZdeFn7fJ7pu8AvSOcff0raB13L1t/T\n4fred62oAIiIFyUtI80Dt7tGuNr2fkmPAB9iy5TNdzv1GYDWfQavqZW/Nn9dN0zvsy10uhrit/nr\nK2rlPQzfpZtiy3ZLBdIk0qF6azufBPwgIs6utOn0B41G0/awPSEdIV8FnEyaF/9BRPzrMI1hpP2C\nrddxAvBqtv7dG+h2rbefTjrynxURLz3fTOkpC3Xb/IrCEqeA/hdwmKQP1isk7SDpnZWi60iHuR8A\nvlo7ZBxKOv+cNF96cq38ZNJh64ODWOZo6fQL8hRpXvOlKxnyyfS3DvMYZldP2JNOcAbw9/n1K9ny\n/Wp5/zCPYbhsD9uTSJc33gycneu3umpuDFsOHC+puv87gfSBeGn7Li8Z6O996wT7Sz9/kt4K7NOm\n7Ta/n6WoIwCAiPiOpCuBBUrPIvo26YTvVNKNYI+z5SqhG4HPk4Ky/gPe+rRzlqSFwKaI6LTjfind\nI2Kz0p2KX86HnYtJf2LzLOD8iKjvrMayfj+15PX8FvAJST8lhdunSA8THM5PO78B7pR0BWlq5wrg\ntohoPXn2HuBLkj4L/JD0zKmj2y5p9G0P27NlAelndhPpxOr24lLSlT23S/or0hz7ZcDdEbG8U8c8\nbbYamCPpYdJRWacbwr5P2r/837w9X0+66mgNW3/PtvkRQHEBABAR50p6gPTp5WuklF5NCoMvVtr9\nUtJyYHM+8VVdxk8lnQucA/wP0qVdb2hVt3vbWv9rlW4s+Vj+9zPgkxHxv2t9xvJdrv2Nr1p2Nuly\nxGtIh9efJ51oO3SAy+k0hptJv2QLSFMV3yJd9dLyZdL35mOkT27fJV3Z8v0ulr8tbS/bMzWM+AdJ\na4B7o/b03jGg39+diHhY0rGkyzT/lnRt/9eAT7dZRjtnkfYT95Cm4/brr31EPCXppNz+dtJj7s8E\n/rzWfijfs0Er9llA3ZD0atKO+aMRsT0d4poNSv50e2tE1HeG7doeSpqunBUR94744GzYFXkE0CTP\nqx4KfJz06eDm0R2R2TbTzSMOdidNmV4CPOid//arxJPA3TiMND3wX0g3fPy2ob3ZeNHNlMBxpPso\nXsvL75Wx7YyngMzMCuUjADOzQjkAzMwK5QAwMyuUA8DMrFAOADOzQjkAzMwK9f8BJ96xiLwlpOIA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116131090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N=1000\n",
    "dx = 1.\n",
    "u = np.random.rand(N,N)\n",
    "\n",
    "times = {}\n",
    "#print 'Python'\n",
    "#y = %timeit -o centered_difference_2D(u,dx)\n",
    "#times['Python'] = y.best\n",
    "\n",
    "print 'numpy vectorized'\n",
    "y = %timeit -o centered_difference_2D_vectorized(u,dx)\n",
    "times['numpy'] = y.best\n",
    "\n",
    "print 'numba'\n",
    "y = %timeit -o centered_difference_2D_numba(u,dx)\n",
    "times['numba'] = y.best\n",
    "\n",
    "print 'Cython'\n",
    "y = %timeit -o centered_difference_2D_cython(u,dx)\n",
    "times['Cython'] = y.best\n",
    "\n",
    "import centered_diff_fort\n",
    "\n",
    "uF = np.zeros((N,N),order='F')\n",
    "uF[:,:] = np.random.rand(N,N)\n",
    "#print u.flags['F_CONTIGUOUS']\n",
    "dx = 1.\n",
    "print 'Fortran'\n",
    "y = %timeit -o centered_diff_fort.centered_diff_fort(uF,dx,N)\n",
    "times['Fortran'] = y.best\n",
    "\n",
    "plt.bar(range(len(times)), times.values(), align='center')\n",
    "plt.xticks(range(len(times)), times.keys());\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Timing versus N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 loops, best of 3: 68.9 µs per loop\n",
      "100000 loops, best of 3: 15 µs per loop\n",
      "1000 loops, best of 3: 232 µs per loop\n",
      "100000 loops, best of 3: 5.03 µs per loop\n",
      "1000 loops, best of 3: 521 µs per loop\n",
      "10000 loops, best of 3: 159 µs per loop\n",
      "100 loops, best of 3: 3.47 ms per loop\n",
      "10000 loops, best of 3: 73.1 µs per loop\n",
      "100 loops, best of 3: 3.25 ms per loop\n",
      "1000 loops, best of 3: 918 µs per loop\n",
      "10 loops, best of 3: 24 ms per loop\n",
      "1000 loops, best of 3: 481 µs per loop\n",
      "10 loops, best of 3: 34.4 ms per loop\n",
      "100 loops, best of 3: 6.23 ms per loop\n",
      "10 loops, best of 3: 99.2 ms per loop\n",
      "100 loops, best of 3: 5.15 ms per loop\n",
      "10 loops, best of 3: 139 ms per loop\n",
      "10 loops, best of 3: 25.5 ms per loop\n",
      "1 loops, best of 3: 403 ms per loop\n",
      "10 loops, best of 3: 21.4 ms per loop\n",
      "1 loops, best of 3: 1.28 s per loop\n",
      "1 loops, best of 3: 217 ms per loop\n",
      "1 loops, best of 3: 2.61 s per loop\n",
      "1 loops, best of 3: 182 ms per loop\n",
      "1 loops, best of 3: 4.73 s per loop\n",
      "1 loops, best of 3: 891 ms per loop\n",
      "1 loops, best of 3: 10.5 s per loop\n",
      "1 loops, best of 3: 761 ms per loop\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEkCAYAAACG1Y6pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXeYVOX1xz9fwIKCil0jCogiCPaCDa5YiYqKBRN7L1Fj\nohg16u7Yosbeu6iJFbvGHkZjYswvauxGE2s0GmOvidHz++PchWGY3Z3ZvVN293yeZ59l7n3nfd/Z\nO9xz3/Oe8z0yM4IgCIKgHvSq9wSCIAiCnksYoSAIgqBuhBEKgiAI6kYYoSAIgqBuhBEKgiAI6kYY\noSAIgqBudDsjJGmgpIckvSDpOUmn1HtOQRAEQWm6nRECvgEmm9kIYGVgTUkT6zynIAiCoAR96j2B\nrDGzd4F3039/I+kZYIn6zioIgiAoRbczQoVIWgDYCtio3nMJgiAIZqVh3HGShkq6WNIzkr6VNK2V\ndiPSPZ8vJL0tKSdpls8haQ5gKnCmmf212vMPgiAIKqeRVkIjgPHAY/i8ZhG1kzQAeBB4DpgADAVO\nx43pMQXtegO/Bp4wszOrPvMgCIKgQ6hRBEwlydLJSJoKzG9m44raHAkcBixlZp+nxyYDzcCiZvZZ\neuwyoJeZ7VHDjxAEQRBUSMO446w8azgeuK/FAKXcAPQFxgBIWgfYA1hV0lPpz4GZTzgIgiDoNI3k\njiuHYbg7bjpm9qakL4HlgLvN7Pc0kHENgiAIWqerGaEBwMcljn+UnisLSY3hgwyCIOiCmJmy6itW\nDEEQBEHd6GpG6CNg3hLHB6TnKsLMVO8fIFfvvip5Xzlt22rTkXOljmf5d2uEa9ddr1+jXruueP06\ne+3aOl/u8UrvseXQ1YzQS8DwwgOSBgJzpee6IvkG6KuS95XTtq02HTlXzpj1IN8g/VXyvnLattWm\n0nPljFcv8g3SX7nvK6dde21aO1/p8cyQNUiIdiFthGgfAUxm5hDtw5gRov15cV+t9N/yoXNA3szy\nGU09qBGSms2sud7zCConrl3XRFICTINs94QaxghJ6gtslr48FOiPGxfwqLevJM0HvIAnq54CLI0n\nq55pZsdWMJZBtn/IoLZISuLhoWsS167rUo17ZyMZoUHAq+nLlkkp/fdgM3szbTccOA9YC98Hugxo\ntgo+SBihIAiCyunWRqiWhBEKgiConGrcO7taYEKmSGpO/ZxBEARBG1TrXhkroSAIgqAsqnHv7GqK\nCTUllBV6NvGQEgTVJ4xQO8SNqGcSDyBBUBtiTyj2hIIgCNol9oQypFy/ZlriKFZCPZC49kEwKxEd\nFwRBEHQrwggFQRAEdSOMUBAEQVA3wggFQRAEdSOMUBAEQVA3up0RknShpH9I+q7ecwmCIAjaptsZ\nIeDXwCr1nkQQBEHQPt3OCJnZo2b2r3LaRrJqEARBeVTrXtntjFAlmFlzFNcCSbtJ+k7S+pIOk/R3\nSV9L+qukXQraDUrbNZXoozk9t2TBsSnpsfklXSHpfUmfSrpd0mJpm30lvSjpq/T3hKJ+p48p6QeS\nnknbvpEe613Q9uy07dAS81tM0v8kXZbV3y0IehLVulc2hBGSNFTSxekN5ltJ01ppN0LSQ5K+kPS2\npJykhvgM3YSTgB2BC4HDge+AKZLWLmpXqczGvcC8wDHApcD3gTskHYVX0b0SOAKYHZiaFjgsZgJw\nAXA7cBjwV6Ap7a+FS9Lfe5R4/6749z2MUBA0EI0iYDoCGA88hs9plpucpAHAg3hp7wnAULy0dy/8\n5hZ0ntmB1c3sfwCSpuLVbg8E/tCJfh83s4NaXkgC+AmwMLC8mX2eHv8t8DSwD3BUUR8rpHP7S/r6\nfEm3ALtJutjMHjez5yU9Buwq6WgzKwxO2QN4wcz+2InPEQRBxjTKKuJOM1vSzCYBL7TSZj9gDmCi\nmT1kZhcDOeCnkvoXNlR6l6sHElaLnypN/4IWAwRgZu8AL+MGvzOcVfT60fT3VS0GKB3vWeDTVsZ7\noMAAtXBq+nvrgmOXAIvhqy0AJI1J+7y88qkHQVBNGsIIWXkqquOB+wpvWsANQF9gbMuB1Of/JmCS\n3pJ0CTXEDNXip0rTf7XEsQ+BBTLu96P092sl2n7cyngvtnFscMGxG4BPgD0Lju0J/Ae4ut2ZBkFQ\nUxrFHVcOw3B33HTM7E1JX6bn7kqP7VWHuXUXvm3leIvRa+thodXvUhsPGe2NVzFm9rWkXwH7SloI\nNz7bAneY2Qcd7TcIgurQECuhMhmAPyUX81F6Lqg+H6a/5y9xbkiVxx7RxrHildYluFHcDfgBvloO\nV1wQdIZZA5QyoSuthDJHUnPBy3yEa7eNmX0m6V1gg8LjkoYAW1F6pZTV/tWGklY2s6fSMYVH8AHc\nVjTPZyX9CQ9G+BR4w8zuz2geQdCj2Ec6YDb4aX9YsBr9dyUj9BEe5lvMAGbsMVSEmTV3ZkI9iEL3\n2HnACZLuwcOlFwf2BZ4FVm/nvZ3hGeC3ks4H3gW2xI3h1Wb2eIn2lzAjHLs5ozkEQc9BWhNovsQ9\nDicCU05x93amdCUj9BIwvPCApIHAXOm5oHO0tmKxonOn4A8DOwMJ8Dy+4lgt/WnrveWOV4rb8Ui9\nI/E9wPeA44DjW2l/PXAm/v24spU2QRAUI62ORx6PxI3Plpj9Nz2X/XCNVt47zU2Z38zGFR0/ApgM\nLFWQV3IY/pS7aFHUXHtjtHzoHG244aLEc/1JE1dfBZrN7LgK3jcH8E88R2l8B8aNax/0LKTV8Pvp\ninji+hWY/WfGaSXANMi2vHdDrIQk9QU2S19+D+gvadv09d1m9hVwEXAwcIukU4Cl8Yz5MyoxQIWE\nO65bsyMwHzNUFIIgKIW0Kn4vXQX4BbBNofFpwczy1UjBbAgjBCwC3Jj+u2WVcmP678HAm2b2saQN\n8D2JO/F9oDMIf39QgKQtgKXw78XzFAUtBEGQIq2M/z9ZDTgZ2B6zr2s9jYYwQmb2OmWEi5vZixRF\nZnWGNDououK6F+fgwRJ/BvYqMxE6CHoO0oq48VkT3+PdAfc2tfO26qhoN9yeUC1o2RNqz68Z+wI9\nl7j2QbdDWgF3u62NG5+LyzE+M3dR3r2zErpSsmoQBEFQKdIoPODrPuD3wNKYnVWpAaoWYYSCIAi6\nI9LySDcCDwB/BIZidgZmX1beFUtJnJP5HAkjFARB0L2QRiBdD/wW3xtdGrPTMPui8q4YIXEV8CRQ\nlaCFHm2Eorx3EATdBmk40nV4Ls9TuPE5tYPGZ02J29K+0nIu+k22E07HisCEttvF5nTPJK590GWQ\nlsMLe26Eq4Sch9lnlXeD0j6OxFNjTgOuMOPLGW2yD0xoiBDtIAiCoEKkZYFjgU3wwpH7Y/Zp5d3Q\nG5gIHAHMiecMXW/GNxnOtlXCCAVBEHQlpGXwlc944GzggA4anzlwDcjDgQ9wGbO7zPguw9m2Sxih\nIAiCroA0FDga2BxPyh6K2SeVd0N/YB/gJ8BzwN7AI2aZlV2piB4dmBA0PpIGSfpOUlO95xIEdUEa\njHQ5Hmb9Om58jqvUAEksKHEcLgi8BrCFGZua8XC9DBCEEQq6Dj0vgibo2UgDkS7Cw6zfAZbBrBmz\nUhWm2+iGJSXOxqPcFgXWNmOSGU9lP+nK6XZGSNJISU9KelnS7ZL61XtOQRAEZSMtjnQe8DTwMTAM\ns2Mwq6h4p8RwiSvxcO3/AiPN2MeMV7KfdMfpdkYIL/lwlJktixe7O7yd9kEQBPVHWgTpTHyf5j/A\ncpgdgdm/K+uGNSRuAfLA34GhZkw2453M55wB3coISVoEGGRm96aHLge2aaN9JKsCknZL913Wl3SY\npL9L+lrSXyXtUtCu1f2Z9G/5naQlC45NSY/NL+kKSe9L+jRdoS6WttlX0ouSvkp/T2h9mvqBpGfS\ntm9IapLUu6jRcpIukPR8OtYXkv4sac+M/lxBkC3SgkinAi8CvYHlMTsUs3+V3wWS2FDiIeAmPMl0\niBknmFHRCqr1Mapzr+xu0XFLAP8oeP0WMLC1xlHUbhZOwvMELsSX7/sDUyT9zcz+UNCu0v2Ze/Fr\ncQywDF6c8A5JtwK7AZfhT34HA1MlLZuW9yhkAjAEryf1LrAlrgi8FF5evIWxwHrAHcBrwNzA9sCl\nkhYys5MrnHsQVAdpfuBQYD/gBmAFzP7R9puKu6A3sBWeYDoXro59bTVyfKpV1A4zq/sPMBS4GHgG\n+BaY1kq7EcBDwBfA23hce6+C86sBfyx43Rf4tEQ/5h+93Xm126Y7/OCG4DvgCaBPwfHFcb2oa9PX\ng9J2x5boozk9t2TBsSnpsXOL2p6eHn8D6FdwfFR6/KSCYy1jfgOsVNTPLem5NQuOzVVibsKfDD8u\n/Hxx7eOnLj8wn0Gzwb8NLjX33lTUB9jsYHuAvQT2ONhWYL0yn+ss45Z376zkp1FWQiPwxKvH8NXZ\nLE/akgYAD+L+0gm44Toddykekzb7B74aamFJZl4ZVR3lVJMoLmuqiqTMBWb2v+ljmL0jKdWN6hRn\nFb1+FM9RuMoKSrOb2bOSPm1lvAfM7C9Fx07FnwK3Bh5P+yiUGJkTXwkJVxIeCwzDK64GQW2R+uOr\n/UOAu4E1Mft7ZV3QD8/r+SnuvtsfyLtt65o0ihG608zuAJDXvZi/RJv9gDmAiemN6yFJ8wDNkk41\ns8/M7F1Jr0sab2b3AHsCN9fqQ0DVjEOteLXEsQ9pw6XZwX5bfNSvlWj7MbBAieMvtnFscMuBNBqy\nGXfBLVHiPQPammgQZI40N/Aj3PX2ILAOZi9X1gX9ccNzIB5wsJUZT2Q807rQEIEJlq7z2mE8cF/h\nkzPuR+2LP+G2sD9wYvoEvxz+tByUx7etHG8xrG1dp1YfaNq4vu2N1xGuxVdZdwE/xHW1NsSFHaFB\nvvNBD0Dqi/QTPEJtNWB9zHasxAClAQfb4w9cQ4F1zNiuuxggaJyVUDkMw58ipmNmb0r6Mj13V3rs\nWWCV2k+vR/Bh+rvUSnVIlcce0caxVwEkzYdLmlxlZgcUNpS0cXWnFwQp7gbeGxcE/ROwMWbPVN4N\nw/BAnEWAHcx4NNN5Nghd6alwAO6qKeYjwsVSE8zl4d8FNig8LmkIvjdTasWTla96Q0krF4wpZuSA\n3Zb+/jYdb6bvdRoOvleGcwmCWZFmR9oPeAXYGJiA2daVGiCJuSROwktx/wZYpbsaIOhaK6HMkdRc\n8DJvZvk6TaXRKXSPnQecIOke4HY8gm5f4Flg9Xbe2xmeAX4r6XxmhGhvAFxtZi1BCZ9Juh/YSdJX\nuNzJUrhY46u4SyQIskWaDdgFD5B6CdiW9DtZWTcID7o6Gw/SWsEaIME0zQ9KqtV/VzJCHwHzljg+\nADqWjGWRJ1RIa6sEKzp3Cn4ddsa/mM/jeTqrMetNvvi95Y5Xittx7asjcffre8BxwPFF7XbC66Fs\nAeyavuco4H/AFa30HQSVI/XB9x2PxdMNdsKsQysWiSG4MvbSwJ5mPJTZPDuJmeWV08PA92nOvv+G\nq6zaEh1nZuOKjj8MvG1mPyw4NhC/+FuY2d0VjBGVVYM2iWsftIqrdGyPJ0u/DxxDB70oEnPibuWD\n8UqmZ5jx34xm2mmUU0u11eOAuWlmJGRbWbUr7QndA2xSJEg6CfgSeLgjHYZsTxAEZSP1QtoWFxY9\nGA+XHtMJA7Qp7sZeCd/3ObnBDNA44HfA2fSa/Sze2uPOqozTCCshSX2BzdKXhwL9YfrC724z+yqN\nfHoBT1Y9BV+2ng6caWbHVjherISCNolrH0zHg2Am4Aot3+Dut3vp4M1TYiCewL0icJAZ92Q11SxQ\nTuvibu4lUJ/jWPc3/6XXbM3Av1l//XUh25VQo+wJLQLcmP675cLemP57MPCmmX0saQN8Y/xOfB/o\nDKiGlzIIgh6PG5/xuCuqD2587uyE8Zkdz2GbDJwL7GjG1xnNttMop9H4Z10GdBzr3vMhvedoMbyH\nAPdD9qW/G2IlVGtaVkL4k02rUXHxNNxziWvfg3HjsyF+Q+6P7/3cilmHb8AS6wPn4yohB5tRkVxP\nNVFOq+H3wlGgE1jnjrfo068ZFwI4BrjDksTSrYtpkO1KqEcboXDHBa0R176HIq2PG5+FcC/LjZ00\nPovh2wbrAD8Gbm8UnTfltCJufFYHTmKtW19k9vma8OqrxwI3WZLM9NnLvXdWNI8wQm23ixtRzySu\nfQ9DWhc3PkviN+ZrMWtNVqqM7uiDBy78HLgUONGML7KYamdRTsvjBnZd4BTWvP7PzLnI0XjqQw74\nlSXJ/0q+N4xQNoQRCtojrn0PQVoTNz7L4pvxV1OgJN+xLlkHuAAP3z7QjJc6Pc8MUE7DcNfiBsBp\nrHblI8w96Eg8v+8E4ApLkjaj88IIZUQYoaA94tp3c6RV8af+FfAb8BTMOhUeLbEQLpi8Ma54fWMj\nuN6U01B8b+f7wJmsfP69zDPicDzZ/GTgYkuSr8rqqwpGqCvlCWVO5AkFQQ9DWhnpdrzy7r3AMphd\n0hkDJNFbYj9cPeQjYLgZN9TbACmnQcrpMuCPwKuMOnkDxk5blnlG3Af8BRhqSXJWBQYoqcY8GyVE\nuy6EbE8Q9BCkFfB9kNF4nuEOmJV18227W1bHXW9fAxuY8Wxn++wsymkgLlW1PXAhy/18LItseCAe\n2XYebnw+qbRfq1J57x5thIIg6OZII/F9kHVxV9lOFFTf7Xi3zA+ciKvH/wy4pgFWPovh2oo7AZey\n9AHrssR2+wCPAJcBwyxJ/l3POZYijFAQBN0PaThufNYHfgnshlmno9MkeuHCuL/AqzaPMOuYgHJW\nKKeFcUO4OzCFJSatxdL77QY8CvwKWN6S5N06TrFNwggFdSP1Mf8W2N3MrqrzdILugDQMz3HZCM/P\n2YuZqzF3omtWxF1vfYDN6l3dVDktCByGlyr5NQtvNJrhR+2A1yG6GVjJkuStes6xHMIIBdORNBf+\nhd4Gr1raH6+m+gQuo/QrqzB3QtIgYDfgVjN7upVmdY8gCro40jJ4BNh4vJT7fngRxgy6Zl48ku4H\n6RiXmWUvX1P2fHIagEffHQDcxIA1RrPCKVvhK597gTUtSRpGkaE9wggFAEgaCtwNLAM8AJwE/BtY\nGH+qvBI3TD+rsOtB+JPpq7j6cBBkh7Q0bhg2x4vBHYRZxZvupbtGuOH5Ja7iv7wZddtTUU7z4qoL\nBwO302/Y2qx60SZ4FYHfAYklyQv1ml9H6XZGSNKFeEGzxc2sR4egl0uqYn4XbjAmmtltRU1+KalU\n0bqKhunEe4NgZqTBwNF4hd3zgKGYfZxd94xI+x0AbGvGY1n1XfFccuoD7I8b23uZY5F1GX39GPxh\n8S/AeEuSv9Rrfp2lO96kfw2sUu9JdDH2wjPGTy9hgAAwsz+b2UUAkp6W9IZKxGtK2k7Sd5J2krQb\nvucDcGV6/DtJ02Z9m3aX9LykryW9LmlyqXlI2krS7yV9LukzSY9KmlCi3euSpklaTtLdkj6V9LGk\nmyQtUv6fJmgopKWQLsFLt7+D5/k0Z2WAJPpJnIKvLm4FVq+zAdoANzRb0qffhoyd9gCjr78b2BbY\nzpJkQlc2QNANV0KWltctJ55dUjNtqGj3ILbF92UuKbP9JbgU/Ua4vHshewIfAzcBi+NuvaOAi3GX\nAXhp7kL2w8t5XJa+d2fgFEn/MLPrWhpJOgB/On0R99EL32+6TdK+ZnZpQZ8GfA/PjbgFLw++ErAv\nMA+wSZmfNWgENFPuy8XAsph9kF33CJiI7yc9DIwyo24RZcppEB5YsQpwKGMe+hr1uh74ANjDkqRD\nhTw7NacqJat2W9keSd+15o4L2Z6ZkfQB0MvMBpTZfl78KfQuM5tUcHwg8DpwoZkdmB5L8NXQbmZ2\ndVE/LefeAYZbupGcugffAP5mZmunxwYAb6VtV7E04klSf+ApfO9qoKX7AZJex8UotzezqQVjnodv\n6C5nZi+38Rl7xLVveKTv4bkvP8SFQE/D7P1sh2AZ/KFqIHCAWccqNWcyl5zmwvddDwTOYrUrbmPu\nwScDQ/FghN9YktTtpt2tZHskDZV0saRnJH1bwkXT0m6EpIckfSHpbUk5SY3rRpSsJj/ZMg9QdiRR\neqO/EdhS0vwFp3bHVyeXVzj+lVYQyWSeyf44HiTRwkbAXMA5VhBym77vHKAfXgOmkLcLDVBKy/ds\naIVzDGqJtBjS2Xj566+B5TD7WZYGSKKvxHHAY8CDwEr1MkDKScppe+AlYBgLjl2PsdP6M/fgafh3\ndpQlyd31NEDVop7uuBF4OOVj6Txm+eOmT78P4iW9J+A3jtNx43lM2mZP/KkB4AAzq5v/FoCu+fT8\nKR6OXQmX4El7OwNnp/tDuwNPmdlTFfb1aoljHwALFLwenP5+vkTbF4ratNcvRX0HjYLv1/0Md7NO\nAUZglqlbLE043Qo4Dd9bWsmMf2Q5RkXzyWkF/EFqAJptF8bcPwi/790LjGzkRNMsqKcRutPM7gCQ\nNBWYv0Sb/YA58Iitz4GHJM0DNEs61cw+M7PLKXryLrVhHrTJc8B6kgab2WvlvMHMHpP0HL4HdDYu\nD78ULo1SKR2u29KJfuM70khICwOHA3sA1wDLY/bPbIdA+MNsM/7Qu68ZD2Q5RkXzyWkBvIzEdkAT\n69z1F/rMfVY6t60sSf5Ur7nVkrq5tay8zajxwH02c8bzDXjZ2bGl3iDpMuBNwCS9JY+kCdqmxWW1\nV4XvuxQYKWl13Bh9hUcnFpKV+6Al+W5kiXMj0t+lVj5BIyMtiHQK7oaaExiF2Y+zNEASktgCX/U0\npz+r1ssAKac+yukAPMDGWOYnYxg7bQ36zH0zXgJ87Z5igKBMIyRpEUlDCl73krSvpLNKhcdmyDCY\nuSCUmb0JfJmemwUz28vMBppZ7/T3PlWcX3fhMuCvwGGtXU9Jq0rav+jwNbi/fjLu3rjZzD4tatPy\nANER91ehAXsA+AI4SFK/gnn1Bw7C97Tq9lQbVIi0ANJJ+PduHmBFzA7E7O3shkASmwF/wmsGnYAb\nn7qV2FZOY3EFku2YY+FNGTvtDRaf8CjwL2A5S5Kri0tqd3fKdcdNAV7BM3XBw2OPBP6G3xT2MrMr\ns58eA/CQ3WI+Ss8FGWBmX0naHFdMuE3S/bhP+gNgIVwEcmOKXG1m9nHqSt0JNxiXlej+edxAHCDp\nS+AT4D0zKxmIUsR0l5mZfSLpcPxJ8XFJU5gRoj0E2NcykmkJqohmkpy5GVgFszeyHQIBm+Irnrnw\n+9UtdZbaWRJXXhgNHMaYh75Ava7D76trW5K0GqnZ3SnXCK1MmkMiqTe+V/NzMztFUg6XkqiGEaoq\naZ5QCz06X8jM/i5pZTyPZhs8J6MfbvCfxG/2xa428O/FTng49SMl+v1a0g74U+hZ+B5fnhlRaq09\nkVrxOTO7UNI/8ZVXU3r4L8DWLfuLRe9vjW4XYdTw+F7uIfiD7B3AapS5/1j+EAh/WGrGV1c5YGqd\njU9f/Pv6Y+BcVr30BPoN/QUe+fkTS5Lf1Gtu5ZKmUiRV67+crRlJXwEbm9nvJK2BV+obaGZvpxO8\n28zm7vAk0sAEMxtXdPw94DwzO77o+OdAk5md3sHxIk8oIwq+D0ea2Sn1nk9WxLXPCHedHoivfu4D\njsPslWyHQHhgTA4PcMoBN5lVLeCl/Tnl1JL86hF4C6ybY+Txu+CBFycD51iSdKqceD2oRp5QuSuh\nt4Hl8Yz37wMv2Qzf7bz4vkA1eAkYXnggTYici6K9oqBuHAj8ly64Eg6qiCcc749HvD0MjMXsxWyH\nQLirOIcnK+eAG+ppfACU00g8YnRhNNuejLl/IK4sch89IOS6Uso1QpcDp0raEDdCRxWcWxOP8qgG\n9wCTJfUriJCbhAcmdDqpLGR7Ooa85MME/MFkR+BiM/tXfWcVNATSHMDe+D3icWBjzJ7JfhgS3Ogs\njoc5X2fG/7Iep6I5eYmFHLADcBzr3PVEGnINsLUlyeP1m13nqbtsj6RdgDVwiZQrWkKsJV0M/KHS\nomSpNMtm6ctD8WTJ5vT13elm+Xx4IuJzeF34pfFk1TPN7NhKxisaO9xxnUBeI+hVPODgHmAvy6hw\nWKMQ175CpNnxZOWf4yoHx2KWedE3iTH4jX4gcDzw6wYwPr3xFIXjgVtY+sBzWWKbybg+4ZHANd0l\n4q0a7ri6accV3Mhgxkax0n8PTkOxkZfpPQ9YC98kvwxoLjPPqLWxwwgFbRLXvkykPrhqxjF4tGwT\nVVAtkVgHNz5D8Jv9NfU2PgDKaV1c7eBz5ljop4y+MQGOAK4ATrAkKU5Z6NLUdE9I0pKVdNRiNCpo\n/zpl5CmZ+5E3qKTvcgl3XBB0EI+S3QGPUnwH2I0S0ZGdH4a1cOOzDB5hebUZ32Q9TqUopyXwlIX1\ngMmMeegT1OvXuCHuliHXNXfHSapk+Whm1jubKVWfWAkF7RHXvhVcPHgb3HX+Cb4C+i0Zu1Qk1sSN\nz3Dc+FxlRt2jyZTTnHik30+BC1nlopvoP+wkvB7XIV0h5Loz1Do6rjBzfh7c6r+AF3r6Fx6NMhFY\nDo+ACYKgu+J6jBPwIID/AocB91bB+KyOG5+ReC2qCQ1ifFo+/xnA08y/xvqMOmUnvBTJycDErhhy\n3QiUmyd0FfCVme1X4tzFwFxmtnMV5lcVYiUUtEdc+xQ3Ppvixmd24FjgjioYn1Vx47MSbnwuN+M/\nWY7RUZTTcDzReiCa7RDG3L8YPsf7gSN7Ush1PfOEtsaX4KW4mRkCmEEQdAfc+GyAG5958b2fWzDL\nNMpLYmXctbca8AtgW7Oq5R1WhHJq+dw7Ayew9u2PM9s8Z+EBVBO7esh1o1CuivbX+AZcKdalesmq\nVUVSc7U224KgyyKNwaWVLsAjU1fAbGqWBkhiRYlbcb3Ch4ClzTivEQyQcuqlnPbEE+L7M2S/hLHT\nVmS2eW4BLgTW6okGqFr3ynJXQhcAx0haALidGXtCWwH74EvTLoeZNdd7DkHQMEij8fDnpfEV0K8w\nyzQMWmIRXIR2HTz374dmfJXlGJ1BOa2Fh1x/w+wLbsVaN62HG+QrcZXrbhVyXQlmlq9GqbayjJCZ\nNUv6CK/iUGqtAAAgAElEQVR4WCjn/y5wmJmdVfqdQRA0PNKquNEZhUeiXYlZ5mHQEiOBu4BrgV3M\n+DLrMTqKclocDzAYB/yMMQ99iHpdjYdcr9MdQ64bhYqSVVMF7YHAorgBesvM6qrT1BEiMCFojx5x\n7aUV8GCANXBvxmWYVSUYQGITvP7UT8xKqrHXBeU0B67uPRm4lJXOuYF5R52Ah1z/xJLk7rpOsMGo\nRmBCRZVVzexbM3vdzP6Y/u5yBiiYFUmJpO/a+Fkjw7EGpXtxK2bVZ1Ah0nCkG/DorkeAoZidX0UD\ntD9ek2xigxmgzXBJsHUYsNoGjJ3Wi3lHPYjrUo4MA1Qbyt0TQtL3gM2B7+FleGfCzCJXqOtzLVAq\n2e7vJY51lEF4mO+rwNMZ9hu0hzQUj/baBNdg3AOzL6o3HL3TcTYB1jXL9HvUYdJ9n+PwkOuDGXP/\nQvj3/gFglCVJZqXFg/YpywhJ2hq4Hl85/QtmSh5r0XsLI9T1edLMrq1Gx5L6F1U+LWs5L98Jncuq\neLPs9vgD5LF4msXZwI+YtQx7xkPSD7gOL7uythkfVXO8clBOLYmwywMnsNatzzH7fGfi97UIua4T\n5brjTsJrYSxiZt8zs8EFP4PMbHAV51g2kgZKekjSC5Kek9Rtiqw1ApLGSHpA0seSvpT0hKQ9SrTL\nS3pN0mBJUyV9CHwiaVc8wxzgygJ337T0fS1uwV0l/UjSC8BXeHY+ktaQNEXSy5K+kPSppEclbVVi\nDlPSvuaRdKGk9yR9lbbPzL3Y0EgLIP0SeAb4GFgWs+NrYICWwGuPvQdsWm8DpJxWVk53ALcAd7L6\n1Zsydtp6zD7fzcBFwOgwQPWjXHfcQOAgM/uwmpPJgG+AyWb2pKTZgAckTTSzW+o9sS7C3JIWLDr2\ntZl9LmkLXLLpHbxa5GfAD4DLJA0xs6ML3mN4afCHgUdxOfuF8f2Hk/BaMxfjNyrwm1UhhwAL4KXD\n3wXeSo9vhW8YXw+8ASwI7ArcImlHM7uuxGe6D1+959L2PwXuljS4u5WfmI5XMz0k/bkJGIXZO7UZ\nmlXxNI6zgdPM6ldKXTmtgCfCjgZ+wSqXHEr/ZSbjYejnAAf25JDrhsHM2v3BfaU/KqdtI/3gX7SD\nSxy39KcZSNp4v9X7M9To75QA37Xycy2+Yn4D+BBYtOB9s+FG5n/A0ILj+fS9x7Ux1i5tnPs3sGCJ\n83OVONYXTyp8vuj4lLSv84qOb5se36edv0nXu/Ywh8HBBv80uNYKrkltvke2Fdj7YFvX8+9AM8vT\nzI0080+a+Qm3njWYadPOZdq0D5g27USmTZu/7teqC/6k/z8t6/8b5a6EfgJcK+kLPKLm4xLGrGFi\n/gHSxNqtgI1aa2NVSFZVPl+TJz9LkmqED1+MPzkX8i4uqTIQOMPMputkmdk3kk4FbgO2xDehp5/G\nV0wd4Woz+3fxwcLvWFrdtS++tzQN2LeoAm8LZxa9npb+HtrBuTUenjqxE77aew7YFLOaBX2kZbYP\nxe8T4834c63GnmkeOS2H731tAJzGiObJLDT2IODPzEg2fb8ec+sOWD2TVXGfMnihplIYUFEpB3mk\nzmS8WN3ywCNmtn6JdiOAc/El9cd4UbuctSEhIi8xPBWvwPrXSubVWapkHGrFK2b22+KDkial/3y+\nxHteSH8X7wu+bx3feyiZGChpYTyZcktgoaLTBswHFBuhV2dqZPZB+h9pgQ7OrXHwD7IV/jf5ENgJ\ns0drOwVmwxUQ1gRGm013ndZuDjkNxY3PpsAZLHPo4Sy++f7Ak/hKfpQlSU3ckUHllGuEZtl8zoAR\nwHjgsXQes6wgJA0AHsSf7ibgT6+n4+6hY9I2ewIHpm/ZH/g/4NfAE2ZW/BQc1I7OrIxneW8aJXc/\nXjrkLPzp9hPgW/z7+UNKBNpY6kcoQVd+WABpHC74OQf+MHcPrX/WKk2B+fCHva/xEOzP2nlLtuPn\nNAQ4GtgCOJtBexzJUjvvCfwFD0JY2ZKkomKbQe0pV7ZnShXGvtPM7gCQNBWYv0Sb/fD/ZBNTN8tD\nkuYBmiWdamafmdnlwOUtb5B0GfCpmR1WhTn3VFryO0aWODci/f1qiXOl6OiNcoX0J2dmucITkvbp\nYJ9dD2l1PLhjMP4gdgMZK1uXNw2G4OKj9wGHmlGzxHXltBTwc7ye2fksvuWKLHPITsBTwL14tNvf\najWfoHNUpJggaXFJ20jaO/29eEcHbuMJtZDxwH1Ffv4b8L2AsSXmtw7+VLyqpKfSnwOL2wUV8yTw\nJrC7pEVaDqYRiJPxjf7by+yr5VpW6g5rucnN9J2VNBIvNVLq+1S3yKzMcZWDm/EIxanAcMyuq5MB\nWhv4PXCeGYfUygAppyWU0wX49/F9Flh3FGOnfcgyhzwBrAqMtSTZJQxQ16LcZNXeuKT73sx8E/hO\n0iXAgW3t0XSCYbg7bjpm9qakL9NzdxWd+z0VGtagfczsu9SY3wr8X3rNPwcm4XsBJ5pZcTZ8a+6u\n5/Hw7gPS6/gJ8J6ZTWulfQsvpO89PA1KeBkP194H37NctcR7urbLDUBaEo/i3Bz4Jb7vUzfVaYkf\n4OHXu5pxT03GdHHRI4AdgcuYZ+RIVj53S+BPwBPAppYkob7RRSl3TygH7I7ne9yI53UsAmyPy198\nQLpHkzEDKBGJB3yUngtqhJndJWkD3Ac/Ga+y+QKwp5ldWdycVlYhZva1pB3wzfSzcHdrnhlRa629\n7ztJm+ERd7sCcwPPArvg1ThXKXcOXQIPwjgKL6h2IZ5oWur/Qo2mg/DN/92BDcx4tupj5rQIrty/\nG3AlfZcYxRrXbAL8AQ/Ln2hJ8n/VnkdQXcot7/0mcK6Z/bLEucPwXJwlOzyJdE/IzMYVHf8vXiri\nnKLjbwFX2cwJkpWM1/KhC/cW8maWL25n3V1JOShJ3a6973keigfb/Bo4EbPiZN4aT4k58ajUZYAt\nzahqOWvltBD+oLMn8CtmX/BU1rppLL4ifBs42pLk99WcQzCDtJhdkr5sgvqU916Y1sUmn8VXRdXg\nI7y0cDED0nOdohp5QkHQIaS+wAG4BuO9wGqYvVbfSYHEQrgb9p/A+lbFGkDKaQHcAO8LXEef/iuy\nzh1r4lGRnwD7WZLMkkIQVJf04TwPIKkp6/7LNUKv4BIt95c4NwmoVi7OS8DwwgOSBuKiiC9Vacwg\nqB1SH9zd1ISHnY/DrFQ+Vs2RWA7fd70RONqMqgRBKKcBeKLrAcBUes2+MuvdtwJwB+5SPQy415Kk\n67pXg1Yp1wgdD1wv3yS9Cd8TWhjfE1of2KE60+MeYHJRJvwkPI/k4c52LqmZEm64IKg6Ui9cQuh4\n3MW0HWZ/rO+kZiCxAZ7oeYQZxXt+2YyR07zAj4GDgdvQbKsz5v6huNGbG99nvj2MT2OQuuWy77fc\n/DZJG+NBCCvjmmHf4JEpTWb2QMUDu/ths/TloUB/3OcLcLeZfSVpPnzz+zm8Hv3SeLLqmWZ2bKVj\nFowdlVWDNqnatfek203wXJ/v8GCfB2udaNoWEnsBJwKTzNwNk2n/OfUHDsJXP78BjmfstMVxg7wI\nfh+40ZKk5uHnQduUe++sqM9Kv/tpuPaCwL+tE5VVJQ1iRoJjyyRaahMNNrM303bD8fDwtfB9oMuA\n5jLzjFobO4xQ0CZVufbSWrjKwSJ4lOEtDWZ8euHzmwhsZlZaPqnD/eck3GtyFp56cRxjp82HG5+h\nuPG51pLkf1mOG2RH3YxQqlLQz0rIwacJq5/ZzAXLGpqi6LhW3XFhhHoumV57aRS+slgJv9FejVlD\n3Wgl5gKuwTX5Jpoxi4Bsp/rPaX481HwUsDNjp32He1ZWxI3QFEuSb7IcM8iW1B03DepjhG4CPjaz\nvUucuxiY18yqtS+UObESCtojk2svDcFvtBsCJwMXYfZ1BtPLFInF8CCAF4B9zPhPpv3ntCnuwbiR\nFU47mwGrngasg6+6LrUkabi/SVCaaqyEyg1MWA8XBy3Fb/DqhEEQAEiL4e62SbgC/P40qKdAYkXg\nTryMx0lm2SX4Kqe58eTi7+Orn/fwUN/rgV0tSRqq/EtQH8o1QvMCX7Ry7j+EekEQgKu+T8bzXKYA\ny1GiLlKjILEZPs8Dzbgh075zGo279/4ArMDYaavjBuhQS5Jrshwr6NqUa4T+hmtXlcoTGs8MleUg\n6Hm4lt3BeJTn7cBKmNW8rk4lSByER+ZNMOOxzPrNaXZc3mcv4ABrsluUz++N7/tsa0nySFZjBd2D\nco3QOcBFqYzOlXj29OK4hteBtO6qa2jKyRMqCGIIgpmRZselZY7By5yvS42LKFaKRB88Om19YB0z\nMlNlUE7L46uft4GVGDvtX8rnT8UL761nSfJKVmMFtacR8oSOxgUV5yw4/BVwvJmdXIW5VY1qbK4F\nPQgvYbETvu/zN+AozJ6o76TaR2IevBRKL2B7Mz7JpN+cegGH4CurI4HLGTutL26QFsSFRj/IYqyg\nvtQ9TyhNHl0LrwXzAfCY1VHZt6OEEQo6hJeN3x0vK/A34HjMOq3cUW3S1c+OuJvsfuAgMzIJEU8L\nzE3BE9h3sSZ7Vfn8oni03V+BvSxJMo22C+pH3Y1QdyGMUFARru6xNy4u+gxufDLbR6kWqfFpWbG9\nBeSyUkBIE093waPfTgNOsyb7Vvn8SFxv7nLghJDc6V7UM0QbSSviJXVXA5YARpvZk5JOAn5nZjUp\ncBUENUPqh5eYPxT4I7AVZn+u76Tap8j4vAnsadZ5rcXp/XuphYtxlYMNrcmeBlA+vwnugvuxJcl1\nWY0XdG/KqkIqaTyu8LsIcBUzG6//4DpQQdA9kOZF+jkuK7U6sAlmWze6AZLoI7EbrjC/C258xmVs\ngLbAy7q8AqxeYID2w91yW4cBCiqh3JXQL4ApZra3XHq+sKbEX/CnxSDo2kjz46rOP8IV3Mdi9mJ9\nJ9U+EoWBEm+Q8coHpouOngmMAyZZk/0OQPl8b+BUXIx4XUuSSNcIKqJcI7QcXtOjFJ8C82cznc4h\n6WE8sbYX/hS7u5l1uvhd0M2RFgJ+CuwD3AaMxuxv9Z1U+xQZn9eB3c3IPA9HOa2He0B+C6xoTa7+\noHx+brz667zAWpYk8X8tqJhyjdD7eBmFB0ucG4H7nRuBzVuEVCWdjm8kH9la46gn1MNxeZ3D8Ii3\nG4BVMXu9rnMqgxoanzlw7budgf2sye6Yfi6fXxyX+3kG2N6S5L9Zjx80FtXKEyprTwi4DjhO0rrM\nKLuApGHAz/CnobpTYIB6Af1w49lW++YwQD0QaSDSecDz+IPYKMz2b3QDJDGbxO74ns9OuPHZoEoG\naEXg/4Bl8dVPoQFaEQ/UuBnYIwxQz6Ba98pyVbTnBKbiQoTvAoviWdGLAvcBE82sIb6Ikn6DR/C9\nAowzs1lyFCJEu4ciDcZXxtviqs5nYPZufSfVPunKZ2d85fMqHmr9u6qMlVNvfHXY8nO1Nc24SSif\n/z4egHCQJUmmenNB41P3PCFJG+Cy9AsCHwIPdrCq6lBc6HEtYHngETNbv0S7EbgK8WjgY/zGkTOz\nNisupiuhX+AlJmYJmggj1MOQlsXVPrbAa9qc1cjCoi3U0vgAKKchwNV41eTdrMnemOl8Pn8gnqYx\n0ZKk4fOkguypa55QOvBDwEMZjDsCFz59LJ3DLJZQrkj8IF7aewKek3A67kI8Jm2zJ65dB3CApQmE\nZvadpKtxyfigpyItj980N8IfZpamCyh8lDA+u1bZ+AgXHD0Jf3g7y5pmPOilEXBn4H/HdSxJXi3Z\nURB0gI6U954bF20cBrwHXGU28xNTGX2opTy3pKnA/GY2rqjNkbg7YCkz+zw9NhmvTLlocSXXVFJo\nDjN7L319LF4mfPcS48dKqDsjrYTfwNfDw4ovwOzT+k6qfVLjswtuOP+Or3wereqYOS0KXIonoO9k\nTfb8TOfz+X74nnBfXAW74Y14UD2qce9sNTBB0umSXi461h94ElfhnYSvSJ6WuzvKxsqzfOOB+1oM\nUMoN+H+GsSXaDwDulPS0pKfxDdVDK5lX0MWR1kC6Ay+0+AdgCGYnN7oBSgMO9sS11nYAdjFjoxoY\noIl4nt/TwJolDNASwO/wh83xYYCCatCWO259Zo16OwxYBtjLzK6Q51c8iAsj7pTx3IZRFBJuZm9K\n+jI9d1fRudeANTKeQ9AV8KjNY4DhwCnA9o1YRruYgpVPixr3LtU2PADKaV7cPbkWsLU1zaqDp3x+\nZVyE9Fzgl6EBF1SLtozQIFyqp5BtgBfN7AoAM3tf0ml4LkHWDMCDEYr5iKjkGkjCH5SOAZbC9zKu\nokGiNNsiNT674m63V4CdzPh9TcbOaRxeE+w3wErWZLNUTFY+vwUuQLq/JcnNtZhX0HNpywj1AaY/\nTUpaAH/SvKCo3Rt4qHaXI01WbSGSVrsCbnw2wY3Pgvhm+rWYfVPXeZWBxOzM2POptfFZAo8Q3BLY\ny5pmFRxWPi+8QuzhwOaWJH+qxdyCxiZNUk2q1X9bRugV/EmzJRpuM0B4XlAhC+Ph2lnzES4HUsyA\n9FynMbPmLPoJaoAbny1w49MXOBG4EbNv6zqvMqiz8RmKJ5Rvg+f3rGBNNkuBOeXzLRVXE2BtS5KK\ngo2C7kv6cJ4HkNTUZuMO0JYROhe4VNK8+MbkQcBreFGsQjbCw6iz5iV85TUdSQOBudJzQU/Ac762\nwfdNvgNOAG6lnVyxRiA1Pi1ut78CO5rxh5qMndMoPDF3Y9x7saw1lc6NUj4/D57O0BsPwc6k4moQ\nlEOrRsjMpsi1tQ7EVyRPAj8qVEaQtDBePz5XhbndA0yW1K8gQm4S8CVkoxAc2nENjKu1T8Jv4J+l\nv++mC1RhLGF8flhD47Mm7nZbAw9P38+aWo8OVD4/ELgbjyY8yJKk4d2aQX2olnZcXSqryitVbpa+\nPBToj+f/ANxtZl+leT8v4KusU3AB1dOBM83s2E6OH3lCjYrUIs55FC4RdTzwQBc1PrlaGJ802XR9\n/G+2DF5a4Qprsq/afF8+vxquGn4mcEZEwAXtUXfZnswGlQbhmeAwQy1B6b8Hm9mbabvhwHl4KOlH\nuGxPc5l5Rm2NH0ao0ZDmAHYDjsC/G8cDD3ch47MbbgRewo1P1WVtUuOzOW70BuARgr+2pvaDNJTP\nb4Unqe5tSXJbVScadBvqLtuTFeZqxe0qeJsXFNugWvMId1wD4KvivXEtweeAHTGrieuqs5QwPj+o\nkfHpDWyXjvstHiF4izW1H6SRRsD9FPgJnoDa0NVig8ahW7nj6k2shBoAqR9ekfdQ4HHgRMz+r76T\nKg+JIbjLcE/gRWq38pkd15Q7AvgXHiF4T6HKdZvvz+dnwwOO1sZDsBulDljQReg2K6GgB+PRlj/C\ny2g/DGyK2dP1nVT7SMyHrz52wSsNXw9MNOOJqo+d01y4wOhkfJ90T+B35RofAOXza+F7RZ/jZbgb\nWsoo6DmEEQpqg6ui/xiPtrwXWB+zF+o7qbZJlQ02xQ3PxsADwC+Be82oujJDKq/TYrB/j0vsVOQ+\nUz6/Cr6/Nir9faUlyf+ynmsQdJQwQkF1cX3BnwL7ALcDa2H2Sn0n1ToSwosi7oyLib6M19jZxyyb\nJOl255DTQsAhwL54qsK4YnHRdvvI50fiqROj8YCFiZYksxR4DIJ606ONUAQmVBHPMTsM2B24EVi1\nkctnSywF7Iivevrghme0GTWrnZNK6xyWzuEGYA1rsorGVz4/DE93GIe733a2JPky46kGPZAITMiQ\nCEyoIq5qcTh+Q/8VcCpm/6jvpEojMQ9e6nsXYCRwE258/mg2a6HFqs1jZmmdK4AzrMneqaiPfH4I\nrma/GZ73c64lyWdtvysIKiMCE4LGRRqMR21thyswj8Ds3fpOalYk+uD7OzvjNaumAecAd5tRU3dV\ngbTORri0zjKldN3a7MMVD36O/93PA5aJuj9BVyKMUNA5vKDhkXgJ9ouAZbHSGmX1It3nWRk3PD/A\nNRCvAQ40o6KbfibzqVBap2Qf+fyi+N99Z+ASYFlLkpp/liDoLGGEgo4hLY8/gW+EP4EPxawmG/fl\nIrEEM/Z55sINz3pm1DwwohVpnR3ak9aZpZ98fkHc3bkX7jocYUnScCvOICiXMEJBZUgr4YrW6+HS\n//s1Uvlsif7ARHyFsAowFU+K/b0ZNVfeLpDWOQqX1jmZMqV1Zuonn58PT+w9AA/0WMGSpCH32oKg\nEsIIBeUjjcf3e04HdsVmrcpZDyR6Axvihmdz4BHgYuBOM+pS5rtAWudIvATFicCt5UjrzNRPPt8f\nzxP6MXAnsJolyWsZTzcI6ka3jI6TdD6wv5mV1KeL6LgOIs0O9MYqcyFVC4kVccPzQ+AfuLvtejPe\nr9ucZpbWeQ83PvdWom4AoHx+LnzVMxl4EMhZkryc8XSDoCIiOq4MJK0HzA3th9hGnlCFFNSSqhcS\ni+NGZ2fcvXUNMM6svoUOC6R1DqOD0joAyufnwBN7j8Rr/IyzJKkoUTUIqkHkCZWBvBzAQ3ihvX/F\nSqh7IDE3sDVueNYAbsGNzyP12OeZaW4urXMAM6R1flGptA5MFxfdHd9vexo41pLkqSznGgSdJVZC\n7XMscJmZ/VsK+9KVSfd51scNzwR8VXAlsLUZdVcAUE6L4bpu+wG/waV1KtbCUz7fB4/gOxb4O7C9\nJckfs5xrEDQydTFCkobivu61gOWBR8xs/RLtRuDS86OBj/Gidjkzm+XpV9IKwBpm9nOFBeqySIzE\nDc+O+J7KNcDhZrxX14kxPdJtNHAQnuh6HR2Q1gFQPt8L2B6X2PkXsIclSSZl64OgK1GvldAI/D/x\nY+kcZvEJylWXH8QLnU0AhuJRWb2AY9I2e+KqzOCSKyMkvVbQx6vA6maVZaEHtUViEWbs8yyMy/1s\nYkZD7IUopzmBSbjxmQ84HzjAmqxiZYK0qNyWwHHAV2mfD0Zp7aCnUq/y3mop0S1pKjC/mY0ranMk\nvsm7lJl9nh6bjD85LmpmbepiSfou9oQaF4m++M14F7zI2u148mXejIrCmKtFKii6Px5w8BS+Kr/H\nmmZdibfblxufTfFyCr3xB6m7w/gEXYlusydk5Vm+8cB9LQYo5QbgFGAscFd7w3RwekGVkOgFjMEN\nz9bAn3DDs50ZjZFz5C63dfEVyob4qmyMNdlfO9xnPj8OOAGYF9/7udWSpK4BFUHQKDRyYMIw3B03\nHTN7U9KX6bk2jZCZ9a7i3IIKkBjOjH2ej/F9nqPNqEgpupoop764S/AgoC8uRbRXpZpuM/WZz6+D\nr3wG4iv46y1JGmKVFwSNQiMboQH4DauYj9JzQQMjsRBeFG4X4HvAtcAEMxqqlLdyWgoPsd4DeBwv\nqfBAR1xu0/vM51fDjc9wfO/n6qhmGgSlaWQjVHXSZNUWImm1k0jMCWyBG5718NXq0cBDZjTMTTh1\nuSXAwbh78CpgtDXZ3zvVbz6/Am50VsOVEra0JKl7gm8QdIY0STWpVv+NbIQ+wn3oxQxIz3UaM2vO\nop+eTLrPsw7ubtsWeBJ3t/3QjIYqqqac5gZ2wiMqe+GBBjtb00z7jpX3m88Px91tY/A9yx9YkjSE\ntFEQdJb04TwPIKkp6/4b2Qi9hLszpiOv2jlXeq7ThGxPx5FYBjc8OwNf4gEGK5jRcMrOymkInli6\nK/Aorm4wrVJJnVn6zeeXBprwqLfT8VyfhgiwCIKsqZZsTyMboXuAyZL6FUTITcJveJkk9cVKqDIk\n5sf3eXYGBuPJmtsAT9WyHHY5pC63jfBAg9G42sJq1mSvd7rvfH4Yru22OV6VdaglScOUswiCamBm\n+WroANRLMaEvsFn68ntAf0nbpq/vNldpvgj32d8i6RSg5anzjKKw7aB2rIDv9RwP3N9I+zwtKKf+\n+J7UQcB/cJfbJGuyTkv9pHs+P8flhM7FS2k3VCG/IOhq1CtZdRDQInXSMgGl/x5sZm+m7YbjobJr\n4ftAlwHNZeYZtTV+JKt2M5TTsrjLbSfgt7iRqFjFumTf+fwauPFZHTgDuMiSJB6Egh5HNe6d3UpF\nu1xa/pBAjtgT6rIop174fsxBwKr4Q8qF1mRvZdJ/Pj8Gj+5bDg84uCICDoKeSronNA3CCHWaWAl1\nbdLyCbvjK59P8VXP9dZkna6imsrrbIwbn0WBXwC/ilDrIOhGsj1B0BGU0wg8vHoH4D482u2xjFxu\nvXCh3KOBOfE8n5siyTQIqksYoaChUU698SCWg4BRwCXASGuyTCR/lM/3xksq/Bz4Gtd4uyO03YKg\nNoQRChoS5TQAL5F9APA+7nK7yZrsP5n0n8/PjgcxHIHX8zkMuC9UrYOgtvRoIxTJqo2HchqFr3q2\nw2V/drAm+1Nm/efzc+I6cT8DXgb2Bh4J4xMEbVOtZNUITAjqjnLqg+/HHAwsC1wIXGJNllk1VeXz\n/fBS3D8F/gycaEnyeFb9B0FPIAITgm6FcloQLxi3P/AP3OV2izVZZpFoyufnw4MZDsbDS8dbkjSU\nkncQ9GTCCAU1RzmtjLvctgZuA7a2Jnsy0zHy+YWAQ/DVz53AGEuSTDQHgyDIjjBCQU1QTrMBE3Hj\nsxRwAbCsNdn7mY6Tzy+OBxnshlfiXc2S5LUsxwiCIDvCCAVVRTktAuyDr0j+BpwJ3G5Nlmn+jfL5\nQXiwwSRgCjDKkuTtLMcIgiB7wggFVUE5rY6verYApgLjrcmeyXycGYrWW+Cit8MsSTJdXQVBUD26\nnRGS9DrwBdCyuf0DM4u9gBqgnGbHQ6sPwiVvzgcOsSb7MPOxXNH6KGAcHtAwNBStg6Dr0e2MEK7E\nPb5FibstIk8oG5TTYri7bR/geVxv7S5rsm8zH2tWReu9QtE6CKpPTyxq1xnKimGPonYdJy0atxa+\n6tkUuB7Y0Jrs+aqMN6ui9Q6haB0EtaNaRe26XbKqpNeAT9KXd+H1h/5X1CaSVTuIcpoT3/w/CJgP\nd4WyF2IAAAlnSURBVLldaU32ceZjhaJ1EDQU3aqekKShwGT8aXp54BEzW79EuxG4z3808DFeMyZn\nZiUFJiUtbmbvSJobuAb4k5mdXNQmjFAHUE5jgRuBp/Brco81lb4OnRpnVkXrk4AbQ9E6COpLd1NM\nGAGMBx5L5zGLNZQ0AHgQeA6/KQ0FTgd6AcekbfbEM+IB9jezPwKY2ReSLgf2re7H6FE8C4yxJvtr\nNTovULQ+Ci/NfSJweyhaB0H3pZ4rIbWU6ZY0FZjfzMYVtTkSTzxcysw+T49NBpqBRc3ss6L2cwF9\nzOxTSX2AS4F/mNkxRe1iJdRAlFC0PoFQtA6ChqNbrYSsPOs3HrivxQCl3IBvTI/F93wKWRS4WVIv\noDfwB/xpOmhAQtE6CIJGj44bhrvjpmNmb0r6Mj13V9G5V4GVaze9oCOkitb7Aofiitbbh6J1EPRM\nGt0IDcCDEYr5KD0XdCFC0ToIgmIa3QhVlTRZtYVIWq0SoWgdBF2XNEk1qVb/jW6EPgLmLXF8QHqu\nU0SyanUJResg6PqkD+d5AElNWfff6EboJWB44QFJA4G50nOdImR7qkOqaH04sAOhaB0E3YKeKttz\nDzBZUr+CCLlJwJfAw53tPFZC2ZIqWh+B53SFonUQdCOqJdtTNyMkqS+wWfrye0B/Sdumr+82s6/w\nG9nBwC2STgGWBpqAM4rCtoM6EorWQRB0lHomqw4CXk1ftkxC6b8Ht6hgSxoOnIfL+3yEy/Y0l5ln\n1NrYkayaASUUrS8KResg6L50K+24etLyhwRyxJ5QxZRQtL4iFK2DoHuT7glNgzBCnSZWQh0jdbud\nByxGKFoHQY+jW8n2BF2Sz/B9ulC0DoIgE2IlFARBEJRFNe6dvbLqKAiCIAgqpUcbIUnN1UrACoIg\n6E5U614Z7rggCIKgLMIdFwRBEHQrwggFQRAEdSOMUBAEQVA3wggFQRAEdSOMUBAEQVA3upURkjS3\npCmSXpL0nKT96z2nIAiCoHW6lRECTgdeMrPlzGwkMLWtxpEnFARBUB6RJ9QOkvoDLwNLmNm37bSN\nPKEgCIIKiTyhthkCvA+cI+kJSbdJWqrekwqCIAhapy5GSNJQSRdLekbSt5KmtdJuhKSHJH0h6W1J\nOUmtzbkPMBK4zcxWBW4HrqrSRwiCIAj+v717DbGiDuM4/v0lIZZSS1CWKEUXUyIqTUi7CFFhgUQY\nvqk3GuGLCCoiJFIJIoxSlMQMswiEJJHAbmIWRtGrvBVeKBS2XohkFlaS1D69mFk7bevu2d2Z8585\n5/eBA3vmcs4zPMw8Z3b+80wBUp0JTQXmAAeAQ/z7ZNUzJHUBnwB/A3OB54GnyB5E17vMQkm7Je0C\nzgd+jYjt+exNwLQyN8LMzEYmVRHaGhGTImI+sP8syywCRgMPRMSOiFhHVoCezK//EBFvRMSNEXFT\nRHwO7JM0PV//LmBfydthiXhASX05d9YoSRGK5kZDzAG2RcRvDdM2AWOAO86yziJglaQ9wBPAghEF\nalU2O3UANmyzUwdg1VHlgQmTgYONEyKiG/gjn/c/EXEgImZFxA0RMTsiDrUgzhEp8lfhcD9rKOs1\ns+xAywxnXlV/ORcdVzvmr6q5g/rlb6S5G2h+yn2vyo/37gJ+6Wf6iXzeiPUON0xNKm6k+HA/ayjr\nNbPsQMsMZ15/0yUtHTSQkhWZu5F8XpXzV9XcQf3yN9LcDTR/qNOLUuUzITMza3NVPhM6AVzQz/Su\nfN6w+SZVM7NqqPKZ0EFgSuMESROB8+hzrcjMzOqpykXoI+AeSWMbps0nG5iwM01IZmZWpFQdE8ZI\nmidpHjABuLj3vaQx+WKvAX8CWyTdKelRYCmwos+w7bJiXCvpR0k9ZX+XFUfSxLzLxv68k/ry1DHZ\n0EjaKWlP3lHlvfzGdasRSWuaPXYmaWAq6XLgcP62NwDlf1+RD8VG0hTgVeAWsutA64FlTd5nNNIY\nbyVriHo0Iqp8xmgNJI0HLouIXZLOBbYDqyNiS+LQrEmSxkXEyfzvV4DTEbE4cVjWJEm3AQuBhyNi\n1KDLt0sX7bJI6nERqi9Jq4HvI2J16lhsaPI+kWuBQxGxInU8NjhJo4EdwP3AsWaOnT64WtuSdBHZ\nzrAtdSw2NJI+BI6SNSVekzgca94SYH1E/NTsCm1VhErqzm0tUHTu8l9km4GVdeicUXdF5y8i7gXG\nA18Aq0oOv6MVlTtJ1wMzIuItDeEO1yrfJzQcvd25vyLbtoG6c39L1p37KrInsp4DPNeySK2vwnIn\naRSwEfg6IlaWHrlBCfteRPRIeht4p7ywjeJyNxOYKulIw3qHgZsj4vhZvz0i2uZFfo0r/3sz8Gk/\nyywGjgNjG6Y9DfwOjOv7eUBP6u3qhFeRuSMbwLIh9TZ10quo/AEXApc0zF8CvJl6+9r5VfRxs2F+\nU8fOtvoXVORbPoimunNLWg90AyHpB0mvFxqs/UcBubsdQNIssu7p05Q9a2q3pMcKD9j+o8B9rwvY\nKmmvpL3ANWTPEbOSFHnc7PvRzXx/u/07rhmTyU4rz4iIbkm93bnfz6c9kiA2G9hAubsW+CAivqTN\nrnW2kUH3vYg4AsxIEZwNqKnjZp/5gw7Phs7cWUvvzm2lce7qzfmrr9Jy14lFyMzMKqITi1Bp3bmt\ndM5dvTl/9VVa7jqxCLk7d305d/Xm/NVXabnrxCLk7tz15dzVm/NXX6Xlrq1Gx+UduO/L304AxuWd\nuiEbOXWKrDv342TduZcDV9LC7tzWP+eu3py/+kqdu7ZqYFqH7tzWP+eu3py/+kqdu7YqQmZmVi+d\neE3IzMwqwkXIzMyScREyM7NkXITMzCwZFyEzM0vGRcjMzJJxETIzs2RchMzMLBkXIbPEJC2T1CPp\n437mbZb0WYq4zFrBRcisOu6WNL2f6W5rYm3LRcisGn4GvgGeTR2IWSu5CJlVQwAvAHMlXZc6GLNW\ncREyq4YA3gW+w2dD1kFchMyqQXlL/BeBByVdnTogs1ZwETKrlo1AN7A4dSBmreAiZFYhEfEX8BLw\nkKRJqeMxK5uLkFn1bACOAc+QXStS2nDMyuMiZFYxEXEaeBlYAFyaOByzUrkImVXTOuAkMBPfrGpt\nzEXILL2gT6GJiFPAyjThmLWOslGhZmZmreczITMzS8ZFyMzMknERMjOzZFyEzMwsGRchMzNLxkXI\nzMyScREyM7NkXITMzCwZFyEzM0vmHy5QvBBt7mQQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11613eed0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "grids = (50, 200, 500, 1000, 2000, 5000, 10000)\n",
    "dx = 1.\n",
    "times = {}\n",
    "for method in ['Cython','numba','numpy','Fortran']:\n",
    "    times[method] = []\n",
    "\n",
    "for N in grids:\n",
    "    u = np.random.rand(N,N)\n",
    "    uF = np.zeros((N,N),order='F')\n",
    "    uF[:,:] = np.random.rand(N,N)\n",
    "    y = %timeit -o centered_difference_2D_vectorized(u,dx);\n",
    "    times['numpy'].append(y.best)\n",
    "    y = %timeit -o centered_difference_2D_numba(u,dx);\n",
    "    times['numba'].append(y.best)\n",
    "    y = %timeit -o centered_difference_2D_cython(u,dx);\n",
    "    times['Cython'].append(y.best)\n",
    "    y = %timeit -o centered_diff_fort.centered_diff_fort(uF,dx,N);\n",
    "    times['Fortran'].append(y.best)\n",
    "\n",
    "plt.loglog(grids,times['numpy'],\n",
    "         grids,times['numba'],\n",
    "         grids,times['Cython'],\n",
    "         grids,times['Fortran'])\n",
    "plt.legend(['numpy','numba','Cython','Fortran'],loc='best');\n",
    "plt.xlabel('N')\n",
    "plt.ylabel('Seconds');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fortran and Numba are the obvious winners here.  Obviously, I\n",
    "need to start using Numba more.\n",
    "\n",
    "\n",
    "## Time to solve the heat equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy vectorized\n",
      "1 loops, best of 3: 1.4 s per loop\n",
      "numba\n",
      "1 loops, best of 3: 825 ms per loop\n",
      "Cython\n",
      "1 loops, best of 3: 7.16 s per loop\n",
      "Fortran\n",
      "1 loops, best of 3: 587 ms per loop\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAAEICAYAAACd0wWxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFUhJREFUeJzt3X2UbFV95vHvw4vKHTGAMYoaApghQSaTrIxRUGOuECNI\n1KiIqKMLM8vIjBpf4kQzE/VGY9ag4mQSIGrAQIxCUEGyJCFipBOMCBokGnmREXwJ+DIIQRAQ4f7m\nj3Ma6tbte+tUd1X3vvd+P2v16j679jm1+1TVU7v22XVOqgpJUpt2WusGSJK2zJCWpIYZ0pLUMENa\nkhpmSEtSwwxpSWrYoJBO8sIkn09ya5J/TXJ6kr3n3ThJ2tFNDOkkzwbeD1wEPAN4PfAk4LwkmW/z\nJGnHlklfZklyFrBfVf3CSNnTgXOBA6vq6vk2UZJ2XEPHpL83tnxL/9uetCTN0ZCe9C8DHwNeStd7\nfhhwCvCDqnrK3FsoSTuwiSENkOQ/A6cCu/ZFnwaOrKpbtryWJGmlhhw4PBL4U+BdwHrgGGAv4Jwk\nTuGTpDkaMtzxReDyqnrRSNkBwFXAc6rqnL7M0+lJ0jJV1ZLH+Ib0hPcH/nlsY18G7uhvkyTNyS4D\n6nwV+PnRgiQHArv1t21iS+8GLUmyoao2rHU7thfuz9lxX87WtrA/J41CDAnpk4A/TnIDcD7wUOBN\nwHXAX6+4hZKkLZoY0lV1cpK7gf8GvIxujvRFwO9U1R1zbp8k7dCG9KSpqvcC751zW1bTwlo3YDuz\nsNYN2I4srHUDtjMLa92AlRo0T3rQhvpxlW1hTFqSWjEpO53nLEkNM6QlqWGGtCQ1zJCWpIYZ0pLU\nMENakhpmSEtSwwxpSWqYIS1JDTOkJalhhrQkNcyQlqSGGdKS1DBDWpIaNuRq4QtJNm7h53Gr0UhJ\n2lENuVr4gcDuo0XAW4CfA/auqo19Pc8nLUlTmpSdQy6fdeXYBu8H/AJwxmJAS5LmY9Dls8YcDuwB\nnDHtipOuiruj8VOHpEmWE9LHAN+oqk8t7y7N6Y75LGmyqWZ3JFkHPAM4az7NkSSNmnYK3tOBdSxj\nqEOSNL1pQ/oY4JqqumwejZEkbWrwmHSSHwGOAP7XhHobRhYXqmphWS2TpO1UkvXA+kF1J82THtno\nscD7gAOr6uolbp84T7qr44HDTpzdIWlidk4z3HEMcPlSAS1Jmo9BIZ3kR4FDgTPn2xxJ0qjBwx0T\nN+Rwx5Qc7pA02+EOSdIqM6QlqWGGtCQ1zJCWpIYZ0pLUMENakhpmSEtSwwxpSWqYIS1JDTOkJalh\nhrQkNcyQlqSGGdKS1DBDWpIaZkhLUsOGnvR/lyRvSHJNkjuTfCPJu+bdOEna0Q29EO1pwJOBDcBV\nwD7AgfNpkiRp0cSQTnI4cDTwH6vqqvk3SZK0aOLls5KcBexeVUdMqOfls6bi5bMkzebyWY8Frkly\nYpJbknw/yUeS7D3LhkqSNjekJ/0D4AfA5cAfAA8C3g58q6oOHqlnT3oq9qQlTc7OIQcOF1d8ZlXd\n3G/0m8DfJ3lyVV04k5ZKkjYzJKRvAr6yGNC9fwTuAh4NbBLSSTaMLC5U1cIK2yhJ25Uk64H1g+oO\nGO64EHhAVR0yUrYTcAfwmqo6uS9zuGMqDndIms2Bw48BP5PkwSNlTwJ2pRunliTNyZCe9O7AvwDX\nc9+Bw+OBK6rqqSP17ElPxZ60pBn0pKvqVuBQ4GbgTOBE4AK6L7hIkuZoYk968IbsSU/JnrSk2YxJ\nS5LWiCEtSQ0zpCWpYYa0JDXMkJakhhnSktQwQ1qSGmZIS1LDDGlJapghLUkNM6QlqWGGtCQ1zJCW\npIYZ0pLUMENakho2MaSTHJtk4xI/v7EaDZSkHdmQq4UvejLdxWcXXTfjtkiSxkwT0p+tqtvn1hJJ\n0mamGZP2Uk+StMqmCemvJPlhkqscj5ak1TFkuOMG4HeBS4GdgecD706yrqr+cJ6Nk6Qd3bKuFp7k\nTOCwqnrISJlXC5+KVwuXNDk7pzlwOOojwNFJfqKqvjZ2hxtGFheqamGZ9yFJ26Uk64H1g+ousyd9\nFHAWsN9iSNuTnpY9aUmTs3O53zg8CrhxvBctSZqticMdST4MXAx8qa//POBo4JXzbZokaciY9NXA\nS4Efp5sr/SXgRVX1gXk2TJK0zDHpJTfkmPSUHJOWNL8xaUnSKjCkJalhhrQkNcyQlqSGGdKS1DBD\nWpIaZkhLUsMMaUlqmCEtSQ0zpCWpYYa0JDXMkJakhhnSktQwQ1qSGmZIS1LDpgrpJI9IcluSjUnW\nzatRkqTOtD3pdwC34pn7JWlVDA7pJE8Cngq8k+4yWpKkORtyjUOS7Az8MfB7wPfm2iJJ0r2G9qSP\nA3YFTppjWyRJYyb2pJM8GHgL8MKquidxpEOSVsuQnvTbgIur6vx5N0aStKmt9qSTHAS8BHhSkj36\n4sWpd3skqaq6Y2ydDSOLC1W1MKO2StJ2Icl6YP2gulVbnk2X5NeAs7ey/ilV9Rt93QKoqi2Oh3R1\nnL3XyVb3laQdw6TsnDQmfRGbp/0RwOv739eusH2SpK3YakhX1XeBfxgtS7J//+dFVXX7vBomSVr+\nuTscs5CkVbDVMempNuSY9JQck5Y0OTs9C54kNcyQlqSGGdKS1DBDWpIaZkhLUsMMaUlqmCEtSQ0z\npCWpYYa0JDXMkJakhhnSktQwQ1qSGmZIS1LDDGlJapghLUkNmxjSSY5K8ukkNya5I8lVSf5nkl1X\no4GStCObdI1DgL2ATwDHA/8GPA7YADwMeOXcWiZJWt6VWZL8PvDyqtpzpMwrs0zFK7NImt+VWW4C\nHO6QpDkbMtwBQJKdgfsDP083zPHueTVKktQZHNLA94H79X9/EPjt2TdHkjRqmuGOg4EnAr8FHAn8\nyVxaJEm61+CedFVd3v/56SQ3AqcnOb6qrh2tl2TDyOJCVS2suJWStB1Jsh5YP6juMmd3/AfgC8Av\nV9Un+zJnd0zF2R2S5je74wn97+uWub4kaYCJwx1JzgcuAK4A7qEL6NcCZ1aVIS1JczRkTPpS4Fhg\nX+Bu4CvAG3AKniTN3bLGpJfckGPSU3JMWtL8xqQlSavAkJakhhnSktQwQ1qSGmZIS1LDDGlJapgh\nLUkNM6QlqWGGtCQ1zJCWpIYZ0pLUMENakhpmSEtSwwxpSWqYIS1JDZsY0kmOTnJekhuS3Jrkc0mO\nWY3GSdKObsiVWV4NXAv8JnAjcCTwwSQ/WlUnzrNxkrSjm3hlliR7VdVNY2UfAA6pqv1Hyrwyy1S8\nMoukGVyZZTyge5cDD19Z0yRJkyz3wOEhwNWzbIgkaXNDxqQ3keQw4JnAS2bfHEnSqKmuFp5kX+AS\n4FNV9Zyx2xyTnopj0pImZ+fgnnSSvYC/Aa4DXriVehtGFheqamHofUjSjiDJemD9oLpDetJJ1gGf\nAB5CN6vjxiXq2JOeij1pSTPoSSfZBfgQ8Cjg8UsFtCRpPoYMd5wMHAG8CnhIkoeM3HZZVd01l5ZJ\nkgZ9meU6YB9gvCtewH5V9fW+nsMdU3G4Q9Lk7JxqdsdK7ui+OoZ0x5CWNINvHEqS1o4hLUkNM6Ql\nqWGGtCQ1zJCWpIYZ0pLUMENakhpmSEtSwwxpSWqYIS1JDTOkJalhU18+S9peLZ5DQR3PLdMGQ1ra\nhDndMZ9b4XCHJDVsUEgn+ckk70nyhST3JLlw3g2TJA0f7ng03dVZLu7X8TOhJK2CoReiTfUVk3wY\n2KuqDh2r40n/p+JJ/1vj83OUz8/VMpOT/tesLt8iSZqKBw4lqWGGtCQ1zJCWpIYZ0pLUsJl/4zDJ\nhpHFhapamPV9SNK2LMl6YP2gutNO3HAK3qw4xak1Pj9H+fxcLZOyc1BPOsluwJH94iOA3ZMc1S+f\nV1V3rLShkqTNDf0yy77Atf3i4grp/96vqr5uT3pa9lRa4/NzlM/P1TKTnnRVfRUPMkrSqvNUpdsw\nz3+8KXt+2h4Z0ts8c7pjPmv75BCGJDXMkJakhhnSktQwQ1qSGuaBQ0lz4eyjTS139pEhLWmOzOnO\n8mcfOdwhSQ0zpCWpYYa0JDXMkJakhhnSktQwQ1qSGmZIS1LDBoV0kkcn+bsk309yfZLfS2LAS9Kc\nTfwyS5I9gU8A/wI8A/hJ4AS6gH/jXFsnSTu4Id84PA64P/DsqroN+LskDwI2JHl7Vd061xZK0g5s\nyJDFEcDf9gG96C+B3YBfmkurJEnAsJD+KeCq0YKq+jpwe3/bNmhhrRuwnVlY6wZsRxbWugHbmYW1\nbsCKDQnpPYF/W6L85v62bdDCWjdgO7Ow1g3YjiysdQO2Mwtr3YAVc4aGJDUsVVs/lWCSbwMnVtVb\nx8pvA95cVSf0y56TUJKWaUvnmx7Sk74KOHC0IMmPA+sYG6uWJM3WkCl4fwP89yQPHJnh8Ty6A4d/\nv1hpuVcdkCRt2ZDhjj2AK+i+zHI88Ci6L7P876p609xbKEk7sIkhDZDkQOBE4BC6WR2nABtqyMqS\npGUbNLujqq6sqsOqal1VPaKq3jzvgE7ynCSfTHJzkjuTXJ3khCR7D1z/gCQbkvzIWPmxSTYmWTef\nlmtckq8mecdat0Pt6F+bG5f4+fgKt7trv+2fnVVb11qTF6JNcgLwKuB9dEMr3wMOovuK+n7Aswds\n5gDgTf02bplPSzVQ4RVJtblbgKcuUbYS96d73V8L/PMKt9WE5kI6ydOB1wC/XlWnjdx0UZL3Ak+Z\ndpOzapukmbq7qi6d1caSPGB0cUD93arqjlnd/7y0+GWW1wD/NBbQAFTVxqr62ySXJvmz8duTnJbk\nsiS/BPxVX3xd/zHq2rHq+ye5IMltSa5M8qwltveKJNf0wy3XJHn12O0bkvy/JD+X5DP9qVwvS/LE\nZf/3M9bvk88meUqSL/T/70VJHt3fvm+/f5621Hojy4v/62OTfC7J7f129k2yd5K/SnJrki8lWb90\nU/LGJN/q6/1Ff6KuxRvXJTkxyVX9fry2X959XvtmWi3sy8Whoy3tyyQ7J7khyZuXaP9CkrPnsnNm\nLMmhSS5Jckf/f56U5N+N3L6+39e/sri/6I6bfa+v8mcjQyj7jDw2L0jy50luBs7tt/XiJJ9K8t0k\nN6UbZv1PY+3Z6mM/T02FdJJd6Q5Onj+h6inAUWMP2gOB5wCnApcBr+tvehZwcP971AeBjwK/BlwD\nnJnkESPbeynwR32dXwU+BJyQ5PVj21kHnA78SX//PwDOTrLbgH95NRSwD/B24K3A84EfoztJ1pB1\nR60D3ks3BPX8frt/AZxF9/3bZwE3AB8e+//T1z8U+C/Aa4Ej6R7H0W3vQnf628P734fS7fdWtLAv\ni63sy6q6BzgNePHonSXZH/hFutdHM/o3lV0Wf/qyg+gy4Dt0Q5tvBl4AfHiJTZwKfB54Ot0+OLQv\nfyvd6/5g4Fsj9d9JN6RyFPAHfdm+wPuB59Lt22/QfXLfb2S9lTz2K1NVzfwADwM2Ai+dUO9BwG3A\nsSNlvw7cCezZL/9qv619xtY9ti8fXXcv4IfAy/rlnYDrgVPH1j2J7jwm9+uXN/TbWj9S52f7sl9Z\n6/3Zt+e0/n971EjZM/s2HkD3BN0IPG2J9T47srz4v/7iSNl/7ct+d6TswL7s8JGyrwI3AutGyl4A\n3AP89BbavQvwhH5bj1zr/bgt7Uu6c76PPy/fQhf6O631fhzbB+M/hwFnAlfTzz7r6z+3v/3gfnl9\nv3zC2HYf2Je/eKx88bH5yIR27dQ/964E3jj0sZ/nvmqqJz1iqweZqup7dO+qx44UHwucW1U3D7yP\ne48iV9VNdO/aiz3pRwJ7s3kv7iy6N4ifGSm7q6oWRpavHNlGK66rqq+MLC+3jXdV1UUjy4vb/OQS\nZQ8fW/eCqrp9ZPmjdD3sxywWJHlRks/3H13vAhbv64Ap2zlPze/Lqvq/wD/Qvz6ShK5n/f6q2jhl\nO+fpFro2j/5cAjwWOKf6JOydDdxN98Y96rwp73Oz+kkOTHJOkm/193EX3Rk+//1Y1Vk99lNp7cDh\nd+mGC/YZUPdUYCHJvsDOwBPpzn091PiZ/e4CFg88LE7z+/ZYncXlvUbKNrnoQVXd1b0mGD2IsdaW\n+l9h+jaOX+BhcTv3bn8L/3/RvQkyUu/2dOd/2Rsg3TGB04GTgTcAN9GF0znLaOc8Nb8ve6cCJyd5\nOd1H/n3oZjq15O6qumy8MMnDGHvtVdU9Sb7Lpq89xusNsEn9/pjHx4Fv0h0P+xpdBp3C5o/prB77\nqTQV0lX1wyT/SDcmudVvM1bVRUmuAV7CfcMTK5pjOeKb/e8fGyt/aP/7phndz2rZ2pHuO/vf9xsr\n35PZTZsL9+27rqCbp/5A7tvXzwU+U1WvGKnT4kUltoV9Cd0nzT8CjqYbp/1MVV09ozbM2zfZ/H/c\nGXgwm7/2pt2v4/UPofsEfVhVfXnk/vZYYt01mSnW4nDHHwKPSfLi8RuS7JTk8JGi99F9pHsR8Odj\nH49W8i73r3Tjd0ePlR9N9xHti8vY5lra2hP5O3Rjbfcepe4Pwj5+xm14yuiBXroDYwV8rl9+APc9\nZoteOOM2zMK2sC+pbmrZGcAr+ts3mw3VsEuAZ2XTi10/m65T+akJ6077ul88KHvvcy/J44GfWKLu\nmsz1b6onDVBVH0vyLuDUJE+gm0p3G/DTdF9muZb7Zn+cDryN7s1m/Em42Gs4LslfArdX1dbC9d53\nyaramGQD8J7+I9Yn6C4VdhzwO1U1Hiat22IPoP9fzwVek+RrdG9Cv0V3Aq1Z9hzuAM5L983DhwPv\nAM6uqsUzKV4AnJTkfwCXAk/jviP1LdkW9uWiU+mes7fTHYzbVvw+3YyNjyZ5N92Y7/HA+VV1ydZW\n7IeIrgOel+QKuk83W/tSy8V0+fKn/f58JN1skuvZ/DFbk550cyENUFWvS/Jpul7AB+je7a6jC+x3\njtT7dpJLgI39wZLRbXwtyeuA3wReSTetZv/Fm5e627H1T0k3Of5V/c83gNdW1f8ZW6f1b9JtqY2j\nZa+gmw52Mt3HybfRHaA5aMrtbK0NZ9C9GE6l+2h+Lt2MhkXvoXt8XkXXC/o43ayFiwdsf7VsK/uy\nq1j1T0muBy6s9i4YvcXXTlVdkeQIuilyH6Gb+/wB4LeX2MZSjqPLiQvohp7221L9qvpOkuf29T8K\nfBl4GfD6sforecxWZNAJllqV5MF04fnyqtqWPs5Jy9L3Ej9UVeOBtVTdg+iG5g6rqgvn3jjNRZM9\n6Un6cb6DgFfTvcuesbYtklbNkK8770U3PPhW4IsG9LatxQOHQzyG7mPw4+gmrd85ob60vRjy0fcZ\ndHPMH8qm3yXQNmibHu6QpO3dttqTlqQdgiEtSQ0zpCWpYYa0JDXMkJakhhnSktSw/w+MNrcadlPe\nMQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10dc16b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def solve_heat_equation_2D(u,dx,dt,T):\n",
    "    \"\"\"Advance solution u to time T using the explicit Euler method.\"\"\"\n",
    "    n_steps = int(T/dt)\n",
    "    for j in range(n_steps):\n",
    "        u = u + dt*centered_difference_2D(u,dx)\n",
    "    return u\n",
    "\n",
    "def solve_heat_equation_2D_vectorized(u,dx,dt,T):\n",
    "    n_steps = int(T/dt)\n",
    "    for j in range(n_steps):\n",
    "        u = u + dt*centered_difference_2D_vectorized(u,dx)\n",
    "    return u\n",
    "\n",
    "@jit\n",
    "def solve_heat_equation_2D_numba(u,dx,dt,T):\n",
    "    n_steps = int(T/dt)\n",
    "    for j in range(n_steps):\n",
    "        u = u + dt*centered_difference_2D_numba(u,dx)\n",
    "    return u\n",
    "\n",
    "def solve_heat_equation_2D_cython(u,dx,dt,T):\n",
    "    n_steps = int(T/dt)\n",
    "    for j in range(n_steps):\n",
    "        u = u + dt*centered_difference_2D_cython(u,dx)\n",
    "    return u\n",
    "\n",
    "def solve_heat_equation_2D_fortran(u,dx,dt,T):\n",
    "    N = u.shape[0]\n",
    "    n_steps = int(T/dt)\n",
    "    for j in range(n_steps):\n",
    "        u = u + dt*centered_diff_fort.centered_diff_fort(uF,dx,N)\n",
    "    return u\n",
    "\n",
    "N = 200\n",
    "x = np.linspace(0,1,N)\n",
    "X, Y = np.meshgrid(x,x)\n",
    "u = np.exp(-100*((X-0.5)**2+(Y-0.5)**2))\n",
    "dx = x[1]-x[0]\n",
    "dt = 0.2 * dx**2\n",
    "T = 0.01\n",
    "\n",
    "times = {}\n",
    "#print 'Python'\n",
    "#y = %timeit -o solve_heat_equation_2D(u,dx,dt,T=T)\n",
    "#times['Python'] = y.best\n",
    "\n",
    "print 'numpy vectorized'\n",
    "y = %timeit -o solve_heat_equation_2D_vectorized(u,dx,dt,T=T)\n",
    "times['numpy'] = y.best\n",
    "\n",
    "print 'numba'\n",
    "y = %timeit -o solve_heat_equation_2D_numba(u,dx,dt,T=T)\n",
    "times['numba'] = y.best\n",
    "\n",
    "print 'Cython'\n",
    "y = %timeit -o solve_heat_equation_2D_cython(u,dx,dt,T=T)\n",
    "times['Cython'] = y.best\n",
    "\n",
    "print 'Fortran'\n",
    "uF = np.zeros((N,N),order='F')\n",
    "uF[:,:] = np.exp(-100*((X-0.5)**2+(Y-0.5)**2))\n",
    "y = %timeit -o solve_heat_equation_2D_fortran(u,dx,dt,T=T)\n",
    "times['Fortran'] = y.best\n",
    "\n",
    "\n",
    "plt.bar(range(len(times)), times.values(), align='center')\n",
    "plt.xticks(range(len(times)), times.keys());\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that numba and Fortran are both much closer to the numpy\n",
    "implementation when placed inside a time stepping loop; this\n",
    "doesn't make sense to me.  I tried inlining everything and applying\n",
    "numba's `jit`, but that gave something monstrously slow.  A possible\n",
    "reason is that numba is falling back to pyobject mode, but I'm\n",
    "not yet sufficiently expert with numba to spot why."
   ]
  }
 ],
 "metadata": {
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}