Skip to content

Instantly share code, notes, and snippets.

@ketch

ketch/numba_comparison.ipynb

Last active August 29, 2015 14:21
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save ketch/ae87a94f4ef0793d5d52 to your computer and use it in GitHub Desktop.
Save ketch/ae87a94f4ef0793d5d52 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
numba_comparison.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- dom:TITLE: Playing with Numba and finite differences -->\n",
"# Playing with Numba and finite differences\n",
"<!-- dom:AUTHOR: David I. Ketcheson -->\n",
"<!-- Author: --> **David I. Ketcheson**\n",
"\n",
"Date: **May 26, 2015**\n",
"\n",
"**tl;dr: Fortran is fastest; numba and Cython are within a factor of two,\n",
"with numba somewhat faster than Cython; plain numpy is much slower.**\n",
"\n",
"\n",
"**Update**: *In the original version of this notebook (posted yesterday),\n",
"my Cython code omitted a type declaration for the loop counters, which\n",
"made it catastrophically slow. Many thanks to my tweeps for pointing that\n",
"out. With the added type declarations, Cython speed is similar to Numba and Fortran.\n",
"Also, some of you pointed out that I was needlessly re-allocating memory\n",
"every time the function was called, so I've now move the memory allocation\n",
"outside of the functions. This significantly sped up the numba and Fortran calls,\n",
"with the result that Cython is still about 2x slower than those (but much faster\n",
"than plain numpy).*\n",
"\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."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"'0.18.2'"
]
},
"execution_count": 1,
"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__"
]
},
{
"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": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"10000 loops, best of 3: 120 µs per loop\n"
]
}
],
"source": [
"def centered_difference_1D(u,D,dx=1):\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",
"D = np.empty_like(u)\n",
"u[0] = 0; u[-1] = 0\n",
"\n",
"times = {}\n",
"y = %timeit -o centered_difference_1D(u,D)\n",
"times['Python'] = y.best"
]
},
{
"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": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"@jit\n",
"def centered_difference_1D_numba(u,D,dx=1):\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"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The slowest run took 147324.47 times longer than the fastest. This could mean that an intermediate result is being cached \n",
"1000000 loops, best of 3: 665 ns per loop\n"
]
}
],
"source": [
"y = %timeit -o centered_difference_1D_numba(u,D)\n",
"times['numba'] = y.best"
]
},
{
"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": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Python:\n",
"1 loops, best of 3: 3.27 s per loop\n",
"Numba:\n",
"The slowest run took 4.46 times longer than the fastest. This could mean that an intermediate result is being cached \n",
"1 loops, best of 3: 86.7 ms per loop\n"
]
}
],
"source": [
"def solve_heat_equation_1D(u,D,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,D,dx)\n",
" return u\n",
"\n",
"@jit # jit only makes a very small improvement here\n",
"def solve_heat_equation_1D_numba(u,D,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,D,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,D,dx,dt,T=1)\n",
"times['Python'] = y.best\n",
"\n",
"print 'Numba:'\n",
"y = %timeit -o solve_heat_equation_1D_numba(u,D,dx,dt,T=1)\n",
"times['numba'] = y.best"
]
},
{
"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": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def centered_difference_2D(u,D,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,D,dx=1.):\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,D,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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Cython\n",
"I essentially borrowed this code from Jake Vanderplas' blog."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%load_ext Cython"
]
},
{
"cell_type": "code",
"execution_count": 8,
"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[:, :] D, double dx):\n",
" cdef int i, j\n",
" cdef int M = u.shape[0]\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)"
]
},
{
"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. When this routine is called below,\n",
"I'm careful to give it Fortran-ordered arrays."
]
},
{
"cell_type": "code",
"execution_count": 9,
"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(inout) :: 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"
]
},
{
"cell_type": "code",
"execution_count": 10,
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Time to compute the 2D centered difference"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"numpy vectorized\n",
"10 loops, best of 3: 32.1 ms per loop\n",
"numba\n",
"The slowest run took 40.73 times longer than the fastest. This could mean that an intermediate result is being cached \n",
"100 loops, best of 3: 3.24 ms per loop\n",
"Cython\n",
"100 loops, best of 3: 5.18 ms per loop\n",
"Fortran\n",
"100 loops, best of 3: 2.15 ms per loop\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEICAYAAACqMQjAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X20XFWd5vHvw0sDmUBy6aEbeuT9LcEeZ3RCjzGSROmI\nIM2LDWTQaYcGesksmdhKENGB3GjbvYLGjC5E2wGiNoKAIErTSSCRiwQQ08AaQBIFCeAEgwPJBUJE\nXvKbP/YuODmpU6fq3sp9SZ7PWnfde/Zb7Tq36vzqnL33KUUEZmZmzeww3B0wM7ORy0HCzMwqOUiY\nmVklBwkzM6vkIGFmZpUcJMzMrFJbQULSEZKWSXpJ0hpJcyXV1pU0TtJCSesk9Uu6StKepTJzJT0k\n6XlJL0haIem0UpkDJG1q8nN1Z0/XzMw6sVNdAUk9wFLgYeAE4BBgPinAXFRT/bpc/iwggHnATcDU\nQpndgSuBR4DXgVOB70l6PSJuKLV3HnBXYfvZuv6bmdnAqW4xnaQLgdnA/hGxIaedD/QCe0fEixX1\nJpMO6FMjYnlOOxK4F5gREctaPOZy4LmIODFvHwA8DhwfEf/SwfMzM7NBaOdy07HAkkaAyK4FdgOm\n1dRb2wgQABGxAlid81pZB+zcJF1t9NfMzLqknSBxOLCqmBARTwEbc16VCeV62cqctxlJO0kaL+nD\nwAzgG03qLpT0mqSnJc2XtGsb/TczswGqHZMAeoD+Junrc16n9fqBA4sJkt4J3J03NwGzIuJHhSIv\nA5cCtwIvAO8BLgAOBk6qfwpmZjYQ7QSJraE8EPIgMAkYDxwPfEXS8xHxXYCIWAvMKpT/iaRngMsk\nvS0iHhyKTpuZbW/aCRLrgXFN0ntyXpV1wF7t1IuIjcD9efPHksYB/wB8t0X7NwCXAW8nBRkk+Za2\nZmYDFBFbjPu2MyaxCphYTJC0LzCG5mMOxXpbjD1QPVZR9ADwFkk7tijjgGBmtpW1cyaxCDhf0tjC\nDKeZpIHrO2rqXSRpSkTcBSBpEmk8YlHNY04Bfh0Rr7coc0r+fV85o1k0HGkk9UZE73D3Y1vh/dk9\n3pfdNRr2Z6urMO0EiW+QxgNulDSPNFg8B/hycVqspMeAvog4GyAifirpVuA7kmbz5mK6OyPix7nO\n/qSFdNeQ1kGMBU4mBaFzCm3PIZ253ANsIC3Gmw3cEBEPt7MTzMysc7VBIiL6JR1Nml10M2k84cuk\nxXRFO7Ll5auZwAJSINgh1y8OQK8H1gCfAfYhzXz6OXBcRCwulFtFCgofJa3PeBK4BPhCXf/NzGzg\naldcjyaNU6ZRcrlpekT0DXc/thXen93jfdldo2F/tjp2OkiYmW3nWh07fatwMzOr5CBhZmaVHCTM\nzKySg4SZmVVykDAzs0oOEmZmVslBwszMKjlImJlZJQcJMzOr5CBhZmaVHCTMzKySg4SZmVUaru+4\nNrMm/BW8m/PNOoefg4TZiOM4kTg+jAS+3GRmZpUcJMzMrJKDhJmZVXKQMDOzSrVBQtIRkpZJeknS\nGklzJbVTb5ykhZLWSeqXdJWkPUtl5kp6SNLzkl6QtELSaQNpy8zMuq/l7CZJPcBS4GHgBOAQYD4p\nuFxU0/Z1ufxZpOka84CbgKmFMrsDVwKPAK8DpwLfk/R6RNzQYVtmZtZliqiebifpQmA2sH9EbMhp\n5wO9wN4R8WJFvcnAXcDUiFie044E7gVmRMSyFo+5HHguIk7stK1WX+ZtNhqk17CnwCbye3mItDp2\n1l02OhZY0ggQ2bXAbsC0mnprGwf1/OArgNU5r5V1wM5dasvMzAahLkgcDqwqJkTEU8DGnFdlQrle\ntjLnbUbSTpLGS/owMAP4xkDbMjOz7qlbcd0D9DdJX5/zOq3XDxxYTJD0TuDuvLkJmBURPxpIW2Zm\n1l3DcVuO8gXXB4FJwHjgeOArkp6PiO8OoC0zM+uiuiCxHhjXJL0n51VZB+zVTr2I2Ajcnzd/LGkc\n8A9AI0i03VaDpN7CZl9E9LXoq5nZdkfSdGB6Xbm6ILEKmFhqeF9gDM3HCYr1jmqSPgG4seYxHwD+\nWtKOEfH6QNqKiN6axzAz267lD899AJLmVJWrG7heBBwjaWwhbSZp4PqOmnp7S5rSSJA0iTSGsKjm\nMacAv84BYrBtmZnZINStkxhPWuj2MGkB28GkxXQLIuLiQrnHSJd1zi6kLQYOJa2zaCyAWxsR03L+\n/qSFdNcAjwNjgZOB/wacExHfbLetQjmvk7BRzeskirxOYqi0Ona2vNwUEf2SjgYuBW4mjQF8mbSY\nrmhHtjwrmQksIAWCHXL9WYX89cAa4DPAPqTZSj8HjouIxR22ZWZmW0HLM4nRxmcSNtr5TKLIZxJD\nZTArrs3MbDvmIGFmZpUcJMzMrJKDhJmZVXKQMDOzSg4SZmZWyUHCzMwqOUiYmVklBwkzM6vkIGFm\nZpUcJMzMrJKDhJmZVXKQMDOzSg4SZmZWyUHCzMwqOUiYmVklBwkzM6vkIGFmZpUcJMzMrFJbQULS\nEZKWSXpJ0hpJcyXV1pU0TtJCSesk9Uu6StKehfwdJH1a0t25zLOSlkiaVGrnAEmbmvxc3flTNjOz\ndu1UV0BSD7AUeBg4ATgEmE8KMBfVVL8ulz+L9O3u84CbgKk5fwzwKeBK4HM57VxguaR3RcT9pfbO\nA+4qbD9b138zMxu42iABnAPsAnwwIjYAyyTtAfRKuiQiXmxWSdJkYAYwNSKW57Q1wL2Sjo6IZcBG\n4MCIeL5QbxnwS1KwOLPU7C8i4medPUUzMxuodi43HQssyQGi4VpgN2BaTb21jQABEBErgNU5j4jY\nVAwQOe1V4BFgnyZtqo3+mplZl7QTJA4HVhUTIuIp0lnA4S3qTSjXy1bmvKYk7QK8g3Q2UbZQ0muS\nnpY0X9KudZ03M7OBa+dyUw/Q3yR9fc7rtF4/cGCLep8FxgOXFtJeztu3Ai8A7wEuAA4GTmrRlpmZ\nDUI7QWJriGaJkj4AfAb4ZEQ8+kbhiLXArELRn0h6BrhM0tsi4sGt2lszs+1UO0FiPTCuSXpPzquy\nDtir3XqSjiSNdXw9Ir7aRr9uAC4D3g5sFiQk9RY2+yKir432zMy2G5KmA9PryrUTJFYBE0uN70ua\nvtpszKFY76gm6ROAG0vtHQbcAtzG5mcMrTQ9GwGIiN422zAz2y7lD899AJLmVJVrZ+B6EXCMpLGF\ntJmkges7aurtLWlKIyEvkjsw5zXS9gGWAI8Cp0dE5cG/5JT8+742y5uZWYdUd0yWNJ40JfVh0mK4\ng0mL6RZExMWFco+RLu2cXUhbDBwKzObNxXRrI2Jazt8NuAfYH/gw6RJVw+8j4oFcbg7pzOUeYANp\nMd5s4JaIOLXweAEQEZ4qa6NSeg23+zlpWye/l4dIq2NnbZDIDUwkzS6aTBpPuBzoLX7ql7QauD0i\nziykjQMWACeTzlpuBmZFxLqcfwDwOOldUe7cExFxUC43kxQUDiWtz3gSuBr4Ql5XUftEzUYDB4ki\nB4mhMuggMVo4SNho5yBR5CAxVFodO30XWDMzq+QgYWZmlRwkzMyskoOEmZlVcpAwM7NKDhJmZlbJ\nQcLMzCo5SJiZWSUHCTMzq+QgYWZmlRwkzMyskoOEmZlVcpAwM7NKDhJmZlbJQcLMzCo5SJiZWSUH\nCTMzq+QgYWZmlRwkzMysUm2QkHSEpGWSXpK0RtJcSe3UGydpoaR1kvolXSVpz0L+DpI+LenuXOZZ\nSUskTeq0LTMz2zpaHuwl9QBLgdeBE4DPAecBc9to+zpgKnAWcAZwJHBTIX8M8CngbuBDwH8FXgWW\nS3pHh22ZmdlWoIiozpQuBGYD+0fEhpx2PtAL7B0RL1bUmwzcBUyNiOU57UjgXmBGRCzLZyO7R8Tz\nhXo7A78Ebo+IM9ttq1A/ACJCA9gXZsMuvYar35PbF/m9PERaHTvrLhsdCyxpBIjsWmA3YFpNvbWN\ng3p+8BXA6pxHRGwqBoic9irwCLBPJ22ZmdnWURckDgdWFRMi4ilgY86rMqFcL1uZ85qStAvwDtLZ\nxKDaMjOzwasLEj1Af5P09Tmv03r9NfU+C4wHLu1CW2ZmNkg7DcNjNr3gKukDwGeAT0bEo4Npy8zM\nuqMuSKwHxjVJ78l5VdYBe7VbLw9EXwt8PSK+Opi2cnu9hc2+iOhr0Vczs+2OpOnA9LpydUFiFTCx\n1PC+pOmrzcYJivWOapI+Abix1N5hwC3AbcCswbTVEBG9LfpmZrbdyx+e+wAkzakqVzcmsQg4RtLY\nQtpM0sD1HTX19pY0pZGQF8kdmPMaafsAS4BHgdOj+XzcttoyM7Puq1snMZ40JfVhYB5wMDAfWBAR\nFxfKPUa6rHN2IW0xcChpnUXk+msjYlrO3w24B9gf+DDpslLD7yPigXbbKpTzOgkb1bxOosjrJIZK\nq2NnyyCRK08kzTaaTBoDuBzoLX7ql7SawgK4nDYOWACcTDpjuRmYFRHrcv4BwOOkd0S5Y09ExEHt\nttXOEzUbDRwkihwkhsqggsRo4iBho52DRJGDxFAZzIprMzPbjjlImJlZJQcJMzOr5CBhZmaVHCTM\nzKySg4SZmVVykDAzs0oOEmZmVslBwszMKjlImJlZJQcJMzOr5CBhZmaVHCTMzKySg4SZmVVykDAz\ns0oOEmZmVslBwszMKjlImJlZJQcJMzOr1FaQkHSEpGWSXpK0RtJcSbV1JY2TtFDSOkn9kq6StGep\nzAxJ10h6QtImSXOatHNAziv/XN3+UzUzs07tVFdAUg+wFHgYOAE4BJhPCjAX1VS/Lpc/i/Tt7vOA\nm4CphTLHAH8K3AacTutvgT8PuKuw/Wxd/83MbOBqgwRwDrAL8MGI2AAsk7QH0Cvpkoh4sVklSZOB\nGcDUiFie09YA90o6OiKW5aLnR8TsnH9STV9+ERE/a6PPZmbWBe1cbjoWWJIDRMO1wG7AtJp6axsB\nAiAiVgCrc14jrdWZQ5k6KGtmZoPUTpA4HFhVTIiIp4CNOa/KhHK9bGXOG4iFkl6T9LSk+ZJ2HWA7\nZmbWhnYuN/UA/U3S1+e8Tuv1Awe28bhFLwOXArcCLwDvAS4ADgbqLlGZmdkAtRMktoZOLjEREWuB\nWYWkn0h6BrhM0tsi4sGu9s7MzID2gsR6YFyT9J6cV2UdsNcA6rXrBuAy4O3AZkFCUm9hsy8i+rrw\neGZm2wxJ04HpdeXaCRKrgImlxvcFxtB8zKFY76gm6ROAG9t43DqVZyMR0duF9s3Mtln5w3MfQLP1\naQ3tDFwvAo6RNLaQNpM0cH1HTb29JU1pJEiaRBqPWNTG49Y5Jf++rwttmZlZE6qbgSppPPAIaTHd\nPNJg8XxgQURcXCj3GOnSztmFtMXAocBs3lxMtzYiphXK7A8cmTevABYD1wMvRcSiXGYO6czlHmAD\naTHebOCWiDi10FYARISnytqolF7DHQ3ZbcPk9/IQaXXsrA0SuYGJpNlFk0njCZcDvcU1DpJWA7dH\nxJmFtHHAAuBk0lnLzcCsiFhXKHMGcGXeDN5cC/FERByUy8wkBYVDSeszngSuBr4QEa+280TNRgMH\niSIHiaEy6CAxWjhI2GjnIFHkIDFUWh07fRdYMzOr5CBhZmaVHCTMzKySg4SZmVVykDAzs0oOEmZm\nVslBwszMKjlImJlZJQcJMzOr5CBhZmaVHCTMzKySg4SZmVVykDAzs0oOEmZmVslBwszMKjlImJlZ\nJQcJMzOr5CBhZmaVHCTMzKxSbZCQdISkZZJekrRG0lxJ7dQbJ2mhpHWS+iVdJWnPUpkZkq6R9ISk\nTZLmDLQtMzPrvpYHe0k9wFLgdeAE4HPAecDcNtq+DpgKnAWcARwJ3FQqcwzwp8BtwEaqvwG+nbbM\nzKzLdqrJPwfYBfhgRGwAlknaA+iVdElEvNiskqTJwAxgakQsz2lrgHslHR0Ry3LR8yNids4/aZBt\nmZlZl9VdNjoWWJIDRMO1wG7AtJp6axsHdYCIWAGsznmNtKozh47bMjOz7qsLEocDq4oJEfEU6dLQ\n4S3qTSjXy1bmvE50sy0zM+tAXZDoAfqbpK/PeZ3W66+pt7XbMjOzDgzHFNh2LjENR1tmZlZSN3C9\nHhjXJL0n51VZB+w1gHpdaUtSb2GzLyL6OnxMM7NtmqTpwPS6cnVBYhUwsdTwvsAYmo8TFOsd1SR9\nAnBjXacG21ZE9Hb4GGZm25X84bkPoGqNGtRfbloEHCNpbCFtJmng+o6aentLmtJIkDQJODDndaKb\nbZmZWQfUahaqpPHAI8DDwDzgYGA+sCAiLi6Ue4x0WefsQtpi4FBgNmnsYB5pKuu0Qpn9SQvjAK4A\nFgPXAy9FxKJO2srlAiAi1NFeMBsh0mvYQ22J/F4eIq2OnS2DRK48EbgUmEwaA7gc6C2ucZC0Grg9\nIs4spI0DFgAnk85YbgZmRcS6QpkzgCvzZgCNDj4REQd10lbdEzUbDRwkihwkhsqggsRo4iBho52D\nRJGDxFBpdez0XWDNzKySg4SZmVVykDAzs0oOEmZmVslBwszMKjlImJlZJQcJMzOr5CBhZmaVHCTM\nzKySg4SZmVVykDAzs0oOEmZmVslBwszMKjlImJlZJQcJMzOr5CBhZmaVHCTMzKySg4SZmVVykDAz\ns0ptBQlJR0haJuklSWskzZVUW1fSOEkLJa2T1C/pKkl7Nil3oqSHJP1O0s8lnVbKP0DSpiY/V7f/\nVM3MrFM71RWQ1AMsBR4GTgAOAeaTAsxFNdWvy+XPIn27+zzgJmBqof13A98HvgacC3wAuEbS+oi4\nrdTeecBdhe1n6/pvZmYDp4hoXUC6EJgN7B8RG3La+UAvsHdEvFhRbzLpgD41IpbntCOBe4EZEbEs\npy0BdoyIPy/UvQXYIyKOytsHAI8Dx0fEv7ToawBEhOqeuNlIlF7Drd+T2w/5vTxEWh0727ncdCyw\npBEgsmuB3YBpNfXWNgJE7sAKYHXOQ9IuwHTSGUfRtcBkSbuXn0sb/TUzsy5pJ0gcDqwqJkTEU8DG\nnFdlQrletjLnARwM7Nyk3Mrct8NK6QslvSbpaUnzJe3aRv/NzGyAasckgB6gv0n6+pzXab1+4MBC\nGZqUW1/Kfxm4FLgVeAF4D3ABKcic1KIPZmY2CO0Eia2ho4uuEbEWmFVI+omkZ4DLJL0tIh5st63G\ntTdLfM3XzFppJ0isB8Y1Se/hzU/8zawD9qqp1/hdbr+nlN/MDcBlwNuBzYKEpN7CZl9E9G1e1XEi\ncXww215Jmk4aE26pnSCxCphYanxfYAzNxxyK9Y5qkj4BuDH//Svg1dz+naUym4Bftmi/8kgfEb0t\n6pmZbffyh+c+AElzqsq1M3C9CDhG0thC2kzSwPUdNfX2ljSlkSBpEmk8YlHu5O+B24FTS3VnAndX\nTa/NTsm/72vjOZiZ2QC0s05iPPAIaTHdPNJg8XxgQURcXCj3GOnSztmFtMXAoaR1Fo3FdGsjYlqh\nzBRSNLsU+CFwHGnR3DERsTSXmUM6c7kH2EBajDcbuCUiTi20VbtOwvPQizwPfaTx67PIr8+hMqh1\nEhHRDxwN7AjcDMwBvpx/F+3YpL2ZpLONK4FvAyuAk0vt30U6K/hzYDFwPHB6I0Bkq4D3At8CbgH+\nC3AJ8KG6/puZ2cDVnkmMJj6T6JQ/qY00fn0W+fU5VAa74trMzLZTDhJmZlbJQcLMzCo5SJiZWSUH\nCTMzqzRc926ybYTvhbU5z8axbY2DhHWB40Ti+GDbHl9uMjOzSg4SZmZWyUHCzMwqOUiYmVklBwkz\nM6vk2U1mtk3y9OzNDXR6toOEmW3DHCeSgU/P9uUmMzOr5CBhZmaVHCTMzKySg4SZmVWqDRKSjpC0\nTNJLktZImiupnXrjJC2UtE5Sv6SrJO3ZpNyJkh6S9DtJP5d02kDbMjOz7mp5sJfUAywFXgdOAD4H\nnAfMbaPt64CpwFnAGcCRwE2l9t8NfB9YBrwfuAW4RtKMTtsyM7OtICIqf4ALgeeAsYW084GXgN1b\n1JsMbALeXUg7MqcdXUhbAiwt1b0FuLPTtnJ6pKfU8jkFhH+I2n3Vzo/3p/fnSN2f3pft78tWx866\ny0bHAksiYkMh7VpgN2BaTb21EbG8EIxWAKtzHpJ2AaaTzhKKrgUmS9q93bbMzGzrqAsShwOrigkR\n8RSwMedVmVCul63MeQAHAzs3Kbcy9+uwDtoahfqGuwPbmL7h7sA2pG+4O7CN6RvuDgxKXZDoAfqb\npK/PeZ3W6y/U6ymkldumVK6urVGob7g7sI3pG+4ObEP6hrsD25i+4e7AoAzHFNgYoW2ZmVlJ3b2b\n1gPjmqT38OYn/mbWAXvV1Gv8LrffU8pfD/zbTvpQf2OvkfI1k+1MEtu6unMTNO/Phm1nfw7/voRu\n7M+RsC9hJOzPge7LujOJVcDE0gPtC4yh+ThBsV6z8YLi+MKvgFfL7ecym4Bf5u2qsYeqsQozM+uS\nujOJRcD5ksYWZjjNJA1c31FT7yJJUyLiLgBJk4ADcx4R8XtJtwOnAt8s1J0J3B0RL7bbVsNAb4Vr\nZmbNKc+RbZ4pjQceAR4G5pFmJM0HFkTExYVyjwF9EXF2IW0xcCgwmzR2MI80lXVaocwU0qjOpcAP\ngeNIi/WOiYilnbRlZmbd1zJIAEiaSDqITyaNAVwO9EahoqTVwO0RcWYhbRywADiZdFnrZmBWRKwr\ntX8i8HekIPB4bvu6Upm22jIzs+6qDRLbM0l/CXwMeDtpAeGTwD8DX4qI37RR/zDgQ6Qzr+cL6WcA\nV5JWsm/cCl23EklPANdHxPnD3RcbGST1Ahc3yVoaEe8bRLs7A58FfhAR/2eg7YwU/ma6CpLmAx8n\nHcznAy8AbwXOIY2HfLCNZg4jvQivBJ6vKWtbV+Ap07al54FjmqQNxi6k9/3jgIPEtkjSXwCfAM6M\niG8Vsu6U9E2gfAPC2ia71Tcz66rXIuJn3WpM0q7FzTbK7xYRv+vW428N/j6J5j4B3FcKEABExKaI\nWCLpZ5IWlvMlfUvS/ZKmAT/KyaslbZL0eKn4QZJuk7RB0kpJJzdp71xJj0p6Of/+21J+r6T/J+k/\nSvppvqX7/fkOuyNC3icrJM2Q9GB+vndKOiLnH5D3z3HN6hW2G8/1zyT9q6SNuZ0DJO0j6UeSXsy3\nnJ/evCu6SNLaXO4qSXsUMsdIulTSqrwfH8/buzdpa1iMhH0p6QlJX6zal5J2lPS0pDlN+t8n6cat\nsnO6TNJ7Jd2r9DUGayV9TdK/KeRPz/v6fY39RRq/fSEXWZjzN0nar/C/+ZCk70haT5qwg6SPSFou\n6Tmlr0T4saT/VOpPy//91uIgUaJ0PXEysLim6OXAKaUXzVjgL4ErgPtJs7EgDbi/M/8uupp0y/OT\ngEeB70n6d4X2/gb4ai5zPHA9MF/SBaV2xgDfBr6eH//3wI2SdmvjKQ+FAPYDLgE+D5wO/BHpZo7t\n1C0aQ5oyPT+3sx9wFelGkX2kffw08P3S81cu/17SLec/CXyA9H8str0TcBHp1vUX5fLXt/Ush8ZI\n2JdBi30ZEa8D3wI+UnwwSQcBR5HeHyNGDmo7NX5y2ltJx4Dfki4tzyGNL36/SRNXAA8Af0HaB+/N\n6Z8nve/fCawtlP8S6ZLWKcDf57QDgH8iLQk4Hfg16crFgYV6g/nfD9xgb228rf0Ae5MW8/1NTbk9\ngA3AGYW0M4GXgZ68fXxua79S3TNyerHunqTFhR/N2zsAa4ArSnW/Rrpv1R/k7d7c1vRCmf+Q0943\n3Psz9+db+bkdXEg7MffxMNIbZBNwXJN6Kwrbjed6VCHtv+e0/1lIm5jT3l9IewJ4FhhTSPsQ6btS\nJlT0eydgSm7rLcO9H0fTvgQOafK6/Bwp6Oww3PuxtA/KP0cD3wN+QZ7ck8ufmvPfmben5+35pXbH\n5vSPlNIb/5sbavq1Q37trQQuavd/v7X2k88kqrUc5IyIF0ifKs4oJJ8B/DAiWt2ypOjWQnvrSJ9a\nGmcSbwH2YctPsdeRAtS/L6S9EhF9he2VhTZGitUR8avC9kD7+EpE3FnYbrT54yZpf1Kqe1tsPpvs\nJtIZxqRGgqS/kvRAvnTwCtB4rMMYOUb8voyIx4CfkN8fkkQ6s/iniNjUYT+3pudJfS7+3Av8GWl2\nUvE4cCPwGumDQ9EtHT7mFuUlTZT0A0lr82O8QrrT9qGlot3637fNA9dbeo50uWa/NspeAfRJOgDY\nEXg3nX3HRfnutq8AjYGvffLvZ0plGtvFr299sVggIl5J70mKg2jDrdlzhc77+GJpu9HOG+1XPP8g\nBWEK5TZK2kDe10pjQt8GLgM+TboH2Z8APxhAP7emEb8vsyuAyyR9jHTJZT/STL+R5LWIuL+cKGlv\nSu+9iHhd0nNs/t6jXK4Nm5XPY163Ar8hjYc+SToGXc6W/9Nu/e/b5iBREhGvSrqLdE262RzqYtk7\nJT0K/DVvXh66tVWdDjTWYfxRKf2P8+/RtpCw1UyPl/PvPyil99C9aavizX2XEqQxpEsDjX19KvDT\niDi3UGYkruofDfsS0pn2V4HTSNfpfxoRv+hSH7a237Dlc9wR+EO2fO91ul/L5SeTriAcHRGNe9Y1\n7nhRNuQzJX25qbn/BUyS9JFyhqQdJL2/kHQl6ZT6r4DvlE5PBxPl/y/p+u1ppfTTSKfIDw2gzeHU\n6o30W9K11jdmaeRJAO/qch9mFCcakAZmA/jXvL0rb/7PGj7c5T50w2jYl0Sa2nkNcG7O32I24Ah2\nL3CypOIx8oOkD9bLm1d5Q6fv+8akgDdee5LeBezfpOyQr/XxmUQTEfHPkr4MXKF0f6kfkQapJ5AW\n0z3Om7Ofvg18gRRwy2+CxqemcyRdC2yMiFYH9zc+JUTEJqUVof+YT3GXkr4y9hzgwogoH8xGuspP\nQPm5/hD4hKQnSUHwPNKNJLv5yel3wC2Svki6jPRF4MaIaNxN+Dbga5I+A/yMdC+x9zZtaXiNhn3Z\ncAXpNbuRNBg8WvwdacbSTZK+QbrmPw9YHBH3tqqYL9GtBmZKeoR0dtdqUd09pOPL/8778y2k2VRr\n2PJ/NuSQUL7ZAAABZklEQVRnEg4SFSJitqS7SZ+CvkuK9qtJAeNLhXLPSLoX2JQH64ptPClpNjAL\n+B+kaW0HNbKbPWyp/uVKi3M+nn9+DXwyIr5SqjPSVxJX9bGYdi5pOuZlpNP5L5AGCN/aYTut+nAN\n6c14BenSyA9JM3oa/pH0//k46VPgraRZO/e00f5QGS37MhWMuE/SGtK93cpjIMOt8r0TEY9IOpY0\nRfUG0tqH7wKfatJGM+eQjhO3kS79HVhVPiJ+K+nUXP4m0tckfBS4oFR+MP+zAfO9mwZJ0h+SDt4f\ni4jRdDptNiD5U/L1EVE+YDYr+1bSpdGjI+L2rd456zqfSQxQvs77VuBvSZ8yrhneHpkNmXZuN7En\n6fLs54GHHCBGLw9cD9wk0mWI/0xaNPNyTXmzbUU7lx9OIK0x+WM2X0tko4wvN5mZWSWfSZiZWSUH\nCTMzq+QgYWZmlRwkzMyskoOEmZlVcpAwM7NK/x/wPdvx/ik7pAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10ad2acd0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"N=1000\n",
"dx = 1.\n",
"u = np.random.rand(N,N)\n",
"D = np.zeros_like(u)\n",
"\n",
"times = {}\n",
"#print 'Python'\n",
"#y = %timeit -o centered_difference_2D(u,D,dx)\n",
"#times['Python'] = y.best\n",
"\n",
"print 'numpy vectorized'\n",
"y = %timeit -o centered_difference_2D_vectorized(u,D,dx)\n",
"times['numpy'] = y.best\n",
"\n",
"print 'numba'\n",
"y = %timeit -o centered_difference_2D_numba(u,D,dx)\n",
"times['numba'] = y.best\n",
"\n",
"print 'Cython'\n",
"y = %timeit -o centered_difference_2D_cython(u,D,dx)\n",
"times['Cython'] = y.best\n",
"\n",
"import centered_diff_fort\n",
"\n",
"uF = np.zeros((N,N),order='F')\n",
"DF = 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,DF,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());"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Timing versus N"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The slowest run took 4.38 times longer than the fastest. This could mean that an intermediate result is being cached \n",
"10000 loops, best of 3: 66.6 µs per loop\n",
"100000 loops, best of 3: 6.98 µs per loop\n",
"100000 loops, best of 3: 18 µs per loop\n",
"100000 loops, best of 3: 4.14 µs per loop\n",
"1000 loops, best of 3: 508 µs per loop\n",
"10000 loops, best of 3: 119 µs per loop\n",
"1000 loops, best of 3: 205 µs per loop\n",
"The slowest run took 4.66 times longer than the fastest. This could mean that an intermediate result is being cached \n",
"10000 loops, best of 3: 55.8 µs per loop\n",
"100 loops, best of 3: 3.26 ms per loop\n",
"1000 loops, best of 3: 744 µs per loop\n",
"1000 loops, best of 3: 1.23 ms per loop\n",
"The slowest run took 4.50 times longer than the fastest. This could mean that an intermediate result is being cached \n",
"1000 loops, best of 3: 365 µs per loop\n",
"10 loops, best of 3: 34.3 ms per loop\n",
"100 loops, best of 3: 3.22 ms per loop\n",
"100 loops, best of 3: 5.16 ms per loop\n",
"100 loops, best of 3: 2.13 ms per loop\n",
"10 loops, best of 3: 131 ms per loop\n",
"100 loops, best of 3: 12.9 ms per loop\n",
"10 loops, best of 3: 20.7 ms per loop\n",
"100 loops, best of 3: 8.74 ms per loop\n",
"1 loops, best of 3: 1.27 s per loop\n",
"1 loops, best of 3: 102 ms per loop\n",
"10 loops, best of 3: 143 ms per loop\n",
"10 loops, best of 3: 54.4 ms per loop\n",
"1 loops, best of 3: 4.76 s per loop\n",
"1 loops, best of 3: 330 ms per loop\n",
"1 loops, best of 3: 524 ms per loop\n",
"1 loops, best of 3: 246 ms per loop\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEkCAYAAACG1Y6pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXeYXGX1xz/fJJTQEzrSApGQ0DuIwA0gEkB6BxEIvYmQ\noNTZAZUmVYqhg4iASBNEmrkigoWi9KII4QeCSBewwfn9ce4mk2k7szuzMzt7Ps+zz2bvvPe972Zm\n77nvKd8jMyMIgiAIWsGQVi8gCIIgGLyEEQqCIAhaRhihIAiCoGWEEQqCIAhaRhihIAiCoGWEEQqC\nIAhaRscZIUkXS/o/SZ+1ei1BEARBddRpdUKSvgi8ALxhZh1nZIMgCDqJjjNC3Uj6LIxQEARBexM3\n6SAIgqBltIURkjRa0hRJT0j6VNLUCuPGSbpf0keSXpOUl9QWv0MQBEFQP8NavYCMccAE4GF8TSU+\nQkkjgPuAp4CtgdHAWbghPbHfVhoEQRA0jLaICUmSZQuRdBMw0sw2LhpzLDAJWMrM/pkdmwx0AYuY\n2YeF8wGfRkwoCIKgvWmLm7TVZgknAHd3G6CMG4DhwEbdByRdBkwDTNKrki5p6GKDIAiChtEu7rha\nGIO746ZjZtMkfZy9dkd2bL8WrC0IgiDoBQPJCI0A3itz/N3stZqR1HofZBAEwQDFzNSoudrCHRcE\nQRAMTgaSEXoXmLfM8RHZa3VjZmr1F5Bv9Vz1nFfL2GpjevNaueON/H9rh/euU9+/dn3vBuL719f3\nrtrrtR7vzX22JwaSEXoOGFt4QNISwBzZa3UjqUtS0vel9Ym0Deaq57xaxlYb05vXarlmK0jbZL56\nzqtlbLUx9b5Wy/VaRdom89V6Xi3jehpT6fUejzfrXilrgxTtQqqkaH8LmMzMKdqTmJGi/c/iuapc\no/uXzgOpmaUNWHrQj0jqMrOuVq8jqJ947wYmmRGaCo2NCbWFEZI0HNgy+/FoYG7cuADcaWafSJoP\neAYvVj0dWBYvVj3HzE6q83oGjf2PDPoXSUk8PAxM4r0buDTj3tkuRmhp4KXsx+4FKfv3KDOblo0b\nC1wArIfHgS4DuqzOXyJ2QkEQBPXR0Tuh/iZ2QkEQBPXTjHvnQEpMaDhtkpgQBEHQ9gyaxIT+IHZC\nQRAE9dOMe+dAUkzod0JZYXATDylB0HwGtRGS1EUPiQlxIxqcxANIEMyMtGdTdDnDHdfDuDBCg5N4\n74PAkVgKOAXYDLQwRGJCEARB0GQkRkp8D3gMeBlYrhnXGdRGKLLjgiAIZkZiuMQxwPO4cMCKZpwE\nWr0p1wt3XPVx4ZIZnMR7Hww2JIYCe+Kut0eBY81m1uWM7LggCIKgoUgI2ByXQ/sQ2M2M3/TX9cMI\nBUEQDFIk1gTOABYDvgXcZka/usciJhQxoSAIBhkSy0j8GLgduAGP+9xazQCFYkIDiZhQ0BPx3ged\niMSCwAnAHsB5wNlmfFT7+aEdFwRBENSJxBwSxwPP4vf9cWacUo8BahZhhIIgCDoUiWES+wEvAqsA\n65lxuBl/b/HSphNGKEDS3pI+kzRe0iRJf5H0L0nPS9qrYNzS2bhcmTm6steWLDh2VXZspKQrJL0l\n6QNJt0laNBtzoKRnJX2Sfd+6aN7p15S0m6QnsrGvZMeGFow9Lxs7usz6FpX0P0mXNer/LQjaFQlJ\nbA08gaddb2fGzma82OKlldBxRkjSipIek/RCdrObq9VrGkB8F/cVXwwcA3wGXCXpC0Xj6g0k/gKY\nFzgRuBTYArhd0nF4J90r8cycWYGbsiaHxWwNXATcBkzCC+ly2XzdXJJ937fM+V/DP+9hhIKORmJd\n4AH873kyMN6M37d2VVUws476Ah4ENs/+fTpwcpkxln11AUmVuazVv08//Z/tjRucR4FhBccXA/4F\nXJf9vHQ27qQyc3Rlry1ZcOyq7Nj3i8aelR1/BZir4PhK2fHvFhzrvuZ/gVWL5rk5e22dgmO/AV4D\nhhSNfQF4qo7/k0Hx3sdX53yBLQd2E9irYPuCDW3s/CTd985GzttROyFJCwNLm9kvskOXAztUGm9m\nXdbg1t4S1h9fjVxzAReZ2f+6fzCz1/Gbd4l7q07OLfr5wez71Wb2z4LrPQl8UOF695rZH4uOnZF9\n367g2CXAovhuCwBJG2ZzXl7/0oOgvZFYWOIi4CHgEWCMGVeY8Wkjr9Poe2U3HWWEgMWB/yv4+VVg\nif5cgBnqj68mLf+lMsfeAeZv8LzvZt//WmbsexWu92yVY6MKjt0AvA9MLDg2Efg3cE2PKw2CAYLE\nXBJdwDO4x2KMGaeZ8XFrV1YfbWGEJI2WNCULOn8qaWqFceMk3S/pI0mvScpLKvwdoq6jb1R6cur+\nf622A6uovmHZXr4X16sbM/sXcC2wpaQFJc0D7AjcbmZv93beIGgXJGaROBjPePs8sKYZR5kxID/f\n7SLbMw6YADyMr6nkpiVpBHAf8BQepB6NxxaG4AFv8F3Q4gWnLcnMO6Ogb7yTfR9Z5rVlmnztcVWO\nFe+0LgEOxWNdHwDDCVdcMMDJNN62xxMOpgFbmvFYa1fVd9rFCP3MzG4HkHQT5W9yBwGzAdtncYT7\ns6fcLklnmNmHZvaGpJclTTCzu3A3zE/765fodMzsQ0lvAJsUHpe0DLAt5XdKjYpfbSppNTN7PLum\n8Aw+gFuL1vmkpN/jWXIfAK+Y2T0NWkcQ9DsSK+CZoHMAh5vRMZ/ntnDHVXHXFDIBuLswkI37/4cD\nGxUcOxj4jqQXgOWZEbwOek+he+wCYAVJd0k6SNLJ+A72Scq70RrlIn0C+KWkb0s6DLgH2Aa4xsx+\nV2b8JcAYYC08BTwIBhxZvc/hQIp/jlfvJAME7bMTqoUxuDtuOmY2TdLH2Wt3ZMeeBJrSfKnDqfQg\nYEWvnY7X/HwVT9l8Gt9xrJl9VTu31uuV4zY8U+9Y/P1+EzgZ731SjuuBc/AnxzBCwYBDYhH8szsS\nVzr4c4uX1BTaTsC02x1nZhsXHf8PMMnMzi86/iqe6ntCHdcIAdMBQla4+hLQZWYn13HebMDfgN+Z\n2YReXDfe+6BlZGoHU3AX3Clm/LfFSwKiqV3DkdRV8GParDz4oCXsAczHDBWFIGh7JObEE66+DOxo\n/dhcrhJZC4ekWfMPJCP0Lu4GKmYEM+pO6sLMuvqyoKD9kPQVYClcweFpipIWgqBdkVgDuA74HbCq\nGe+3eEnA9CLVFKCcbmRfGUhG6DlgbOEBSUvgPv/nyp7RA9lOKHZAncX5uNzQI8B+NSa9BEHLkBiK\na7wdBRxhxvUtXlJZmtXUbiAZobuAyZLmKsiQ2wX4GPhV65YVNBMze5k6sjjNbFTPo4KgPZBYkhlK\nHmuaMa2V62kFbZGYIGk4sGX249HA3Lg7BeBOM/tE0ny4PMVTeIbWsrjv9BwzO6nO60ViQlCVeO+D\nZiOxK75zPxs4s9Fab82gGYkJ7WKElmZG1Xv3gpT9e5SZTcvGjcXrVNbD40CX4VlTdf0S3f+RQJ4q\n7ri4EQ1e4r0PmoXEPPh9bG1gDzMebfGSaiJzx02FDjRC/U3shIKeiPc+aAYS6wM/BO4GJlkbtNeu\nh2bshNpCMaFVZN1Ak1avIwiCziYTHT0FlxE70oyDB6ABSpoyb+yEqo+Lp+HBSbz3QaOQGA38CBcA\n3seMN1q8pF4TO6EgCIIBQqb7NhHXVrwW2GLAGiBpKaSrmjH1oDZC4Y4LgqAZSMwP3AQcAYw34/tm\nTeuI3DykBZDOAR67tEn2YlAboWa09w6CYHAjsSnwJ+BlYB0znmrtinqBNBfSibgQwCzACvub7dWM\nSw1qIxS0P5KWlvRZM+RCgqCRSMwmcRZwFR77OdqMf7V4WfUhzYp0KN61dSywDmaHYdY0N+KgNkLh\njhtQDDxXRjBoyJrO/R5YGljFjHtbu6I6kYYg7Q48C3wF2AKz3TH7y4whIdvTcELANAiCvpC13D4M\nOAn4JnDlgIr9eIfiLwOnAv8B9sNsarmhZpb68MYyqI1QEARBbylqOvcFM15s8ZLqQ1oHOA1YFDgO\nuIUKNTvKS8DmzVjGoHbHBY6kvbO4y3hJkyT9RdK/JD0vaa+CcRXjM5lr8zNJSxYcuyo7NlLSFZLe\nkvSBpNskLZqNOVDSs5I+yb5vXXmZ2k3SE9nYVyTlJA0tGrS8pIskPZ1d6yNJj0ia2KD/riBA4ivA\n47ha+xcHlAGSxiLdjGfvXQusiNnN5QyQ8hqivLbDf88zmrGcQb0TilYOJXwXmB24GN+aHwxcJenP\nZvZQwbh63Q2/AF4FTgQ+j6et3i7pFmBvXAPw39nxmyQtl6lnF7I1sAyuufUGsA2Qw3sH7VswbiNg\nA+B24K/AnMDOwKWSFjSz0+pcexBMp6jp3E5mPNjiJdWOt77J4X9LZwJ7YPZJ2aF5DcX/bo4HPgFO\n5mQ+AH7Z6GUNaiMUMaESZgXWMrP/wfRW6y/hPu+Hqp3YA78zs8O7f8j8yt8AFgJW6G7NIemXeGrr\nAbh7oJCVs7X9Mfv5QvnT3N6SppjZ77LjPzSzKYUnSjoX/+P5lqTvdf9+QVAPWdO5HwF/oI2azvWI\nNBI4Fn9YuwRYDrP3yg7NaxZgz2z83/GuBvdYzozc9L/dhjKojVAzUF79EpS0XFMkZS4qvEGb2euS\nXgBG93Hec4t+fhA3QlcX9IbCzJ6U9EGF691bYIC6OQPYFtgO70aJmX3c/aKk2fGdkIB78V3SGLzj\nahDURNZ0bhJ+Q27bpnMlSHMAX8eb5f0UWAmz18sOzWt2YB88ueJFYH/gAcvNcNEpTRdtxjLDCDWY\nJhmH/uKlMsfeAZZo8Lzd7dj/Wmbse8D8ZY4/W+XY9EZ2kubCe1HtDCxe5pwR1RYaBIUUNJ0TA6Xp\nnDQLMBF3f/8GWB+zF8oOzWtO3PMwCXgM2NVy9tuZxqTpUOAgZvR4ayhhhIJCKjXV6jas1XZ5FT9L\nVfo99XS93nAd3iBxCvAA8HZ2nS3x3Vck4wQ1IbEL8H0GStM5aQiwI/BtYBqwDWaPlB2a1zzAocCR\n+N/JVpazx0vGpemq+N/Sf4AEGq/+0HFGSNLFeLHVYmZW9YYTiQl18072fWSZ15Zp8rXHVTn2EkDW\nfXcr3M13SOFASZs1d3lBp1DQdG4dXHS07I28rZA2xdOtAQ7B7L6yw/IaibvoDgXuAsZbzp4pGZem\nc+FNP/fE47NXMn78hs1YescZITxwmIOe1WojMaE+zOxDSW8AmxQel7QMHpspt+NpVIxsU0mrmfnT\nmjxCekz22q3Z90+z68308JGlg+/XwLUEHUpB07l7gNXbvuePtCZeaLo0nsl2E2aflQzLa2E8NrQf\ncDOwruXsz2WnTNNt8LbjvwJWsiT5OwBRrFobZvYgNCeLYxBT+J95AfBtSXcBtwGLAQcCTwJr9XBu\nX3gC+KWkC5mRor0JcE13ZlxmJO8B9pT0CV7bsBTu834JWLNBawk6DIlZ8BjKAcCBZtzW4iVVR1oO\nd7t9ETgZuByz/5YMy2txYDLwVfwBfTXLWdm4ltJ0Cdz9OBbYx5Kk4enY5eg4IxT0mkq7BCt67XRg\nXvxDneCZZvviN/jim3zxubVerxy3AS/gqaNjgDfxP75TisbtibslvgJ8LTvnOOB/wBUV5g4GMVnT\nuWvxpJjVzPhbi5dUGWkx3NOzA16vtA9mJbs15bUMnum2E3A5sILlrOzvpTQdBhyO76S+D+xqSdJv\nwqst66wqaTRuodcDVgAeMLPxZcaNw/9j1sU/JJcBeSuz5Sw677NKMaHorBr0RLz3nU+m+7YP/mB1\nCnCBGVXvKy1Dmgd/ADsAf5g6DbO3S4bltXw2bku86Pw8y9k/Kk6bpmvjiQfvAAdbkpTNopuxjMZ3\nVm3lTmgcMAHvOjiMMk/AkkYA9+EZGVvj9SNn4T7/E7MxE/FiSoBDzOzhpq88CIIBTdZ0bgqu4DG+\nbXv+eFxhV+B7eK3bqpi9WjIsr5Xxncx44Dzg65YrX5AKoDSdF/gOnk03GbjWkqQlO5JW7oTUnbqb\nVeaPNLONi8Yci+evL1VQVT8Zz1dfxMw+rDQ38GnshILeEu99ZyIxG3AIvlu4FjiubXv+uBfoAjwb\n9VDMflMyJK+1ceOzFv6APsVyMwrAS8anqXDDcy7wc+CbliTvVBpfuqQO2glVqR0pZAJwd2FVPXAD\nvn3eCLij+ARJl+G6TibpVeAuMzugAUsOgmCAkrnedsTjhc/iu5/2VM7wgusT8VjrycDFFElNKa8N\nceMzFr8f7mq58jpw089J01HAhcCSwC6WJG2he9fuiQljcHfcdMxsmqSPs9dKjJCZ7ddPawuCYAAg\n8QV8lzA7cIAZ97d4SeVxD84OeHFsisvsTC81ydopbIobqMXw1OwfWs7+U3XaNJ0FlxyahLv1zrYk\nqXpOf9LuRmgEnoxQzLuE/EoQBFXIst5Ow4tOjweubePEgzF4AtaiwJ6YPTDTy3mtg9fuzI3Hcm6w\nXM9CvErT9fHY16vAWpYk5aSyWkq7G6GmkikmdBPKCUHQAWRJByfi6fpnAV81o6qrqmVIc+IlBAfi\nxuWCwnqfTNvt28Bu+E7mOstVzwwGUJqOxN10W+ByVT/pbeJB1tY76c25tdDuRuhdvCalmBHMEMHs\nNaGYEASdg8TseKbsN4EbgXFm/L21q6qAu962wRMEHgJWLla4Vl6b4buYXwMrVku1nn6OJx7sgbvd\nbgLGWZL0qeVE9nCeZssuaWjZV9rdCD2HB96mI2/MNEf2Wp8I7bggGPhIDAF2wZsyPgFsYNb3+0PT\nkJbFXW9LA/tiNpMygfKaH9/BJcBBlrNf1DRtmi4HXISr0G9tSfL7Bq66e0fUcNrdCN0FTJY0V0GG\n3C7Ax7iuURAEgxiJDfGnfgF7m7XxfUEaDnwLFw89AzgXm5FUkCUe7Izvjm4EVrJc+TKUmaZN09my\neQ/HXXrftyQZMI0bW2aE5G/IltmPnwPmlrRj9vOd5m1nf4C3fL5Z0unAsrhkxdlFadu9ItxxQTAw\nkVgOj3mshsdUrm/bpAMAaSs8seARYLXigtNM4+0i/B63veVqK7pXmo7H75PPAKtZkpQUsjYKa5KA\naSuLVZdmRrOz7kUo+/coMxfZkzQWL9haD48DXQZ01VhnVOna3efmqeKOi4LFwUu89+2JxILASbiK\nwJnA+W1bbAogjcIVDMYAh2N2z0wv5zUE72L6bbyG51TL2b97nDZNF8RddhsBh1uS3N7opZdc091x\nU6GxxaotM0KtJBQTgp6I9769kBiO98GZhKtBn2JGj4H6luGt5Sfjaz4LOBub2bgor+WAS4HZgP0s\nZz1KBylNh+B6d6fiLSdyliR99grVSkcpJrQDkZjQWrInq18C+5jZ1S1eTtCGZEkHe+A7hUeA9cx4\nsbWr6gFpczzx4ElgDcxemenlvGbBjenRdAun5qzHrq1K0xVw19uswGaWJH9s9NKrXn+QJiY0lYgJ\nzYykOXCV3h1wgdm5cXXdR/FA6bVmPf+xFM25NLA3cIuZ/anCsMG3HQ96RGI8nnTwX2APM9pCZqYi\n0pLAOcAquOvtrpIhea2BhxT+DqxpOXu5x2nTdARe9/RVPCY+xZKk31uNNysmNKiNUDCDrLXGnbiq\n8L14uus/gIWALwFX4obpm3VOvTTuw38JqGSEgmA6EmPx7LEV8Kyvn5i18YOKNCu+qzkaj//sgdlM\ncSrlNQcuvPw13E33Q8tVj4VkfX4OwA3PrcCKliRvNnz9LWZQG6FwxzlZpuIduMHY3sxuLRpypryN\ncF86k0Z8JaiKxML4jXpHPOaxoxk9BulbirQpnjj1IrA2Zi+VDMlrPHAJ7lFY2XLWoyFRmm6G76re\nxF1vLX+Aa5Y7LhITehg3GILTkg7Hn+BOM7Pjahj/J2A+YOniLEVJO+FK53vhDznlupn+yszGF8aE\n8B5Rk/AU1TeAC83szDLX3hZ/klwFd+P9CTjDzG4vGvcy8FfgYDwwvAHwGb7LO8ys+o1gsLz37YDE\nHMBRwJHA1cB3zKi5vUBLkBbHhUbXAo7A7GclQ/KaD8/g2xw4xHKlY0rOSdPlcRfkGPzv4fZW9fkp\nRzMSE8r22wkGHTviN/RLahx/CbAE7qYrZiIuOvsTvKD4u9nxKbiW1554kLmQg3Cf94/wm9HfgNMl\n7VY4SNIhwM24AczjQd2RwK2S9i+a0/D6s6nAy2S6W8D2wDU1/p5BE5EYKrE38DywErCOGUe3tQGS\nZsF7mv0RV20ZV8EAbQc8jbeVX7EnA6Q0Hak0PReX6JmKu95uaycD1CwG9U6IqBMCQNLbwBAzq0mZ\nXNK8wOvAHWa2S8HxJfAb/sVmdlh2LMF3O3ub2TVF83S/9jowtrtJYeYefAX4s5l9ITs2AlcCfh1Y\nvaDJ4dzA43jsagkzez87/jLeN2VnM7up4JoX4E3Nljeziq2MB8t73yokNsWf+D8CJpnR/h2RpfG4\n620annjw55IheS2SjVkJ2N9yM6thl4z3NgsH4nHTnwInWZK81eilN4Jm1QkN6phQU7LjZhi45tLY\nG+Q8+O6jxkvb+5JuBHaTNNLMup9c98FjP5fXef0rC7vkmtknkn4HrFsw5ku4ZuD5hWoZZvahpPNx\nqZNN8T/kbl4rNEAZU3EjNBqoaISC5pDFfa4ElsOTXG5u66QDAGlR3KW7Pu4yvJViN7RL7uyDt464\nDNjTcla1iFZpujnu0nsN2MSS5MkmrL5hRHbcQGFgPj1/gKdj18MleKbPV4Hzspbq+wCPm9njdc5V\nEswF3saFGLsZlX0v1w3zmaIxPc1L0dxBPyCxGHA/num1rRlt01itLNIwXI/tOPzzvj9mH5UMy2uZ\n7PURwGaWs6r1O0rTsbhRG427n+8cDG63SoQRCgCeAjaQNMrMamp6ZWYPS3oKjwGdB2wCLIWn1tZL\ns2oeqs07EB8WBiwSS+IG6AozTm31enpE2gCX0XkT+CJmz5cMyWsYrohwLK5jd061RnNK0/nxdOvd\ncKHRbdupw2mrGNRGKFK0p3MTnj22H96BslYuxXdBa+HG6BM8uaCQRj3h/SX7viKZX7qAcdn3cjuf\noMVIjMJjf+ebcU6r11MVaWH8QWpjfJdyU7HrDUB5rYy7nT8E1rVcaXxo+liP+xyC/23dCIy1JGlf\nyaEKNCtFu6bsOEkLS1qm4Ochkg6UdK6krZuxsP7AzLrCAAHuw34emFTp/ZS0hqSDiw7/EPgXnjK9\nLfBTM/ugaEx3/KY37q/CP/578SD24ZLmKljX3LjL5MNsTNBGSHwez5I8s60NkDQM6TDcK/AmMBaz\nn5SJ/cyuvL4D3IdL6GxSyQApTaU03RKX79kCGG9JcthANEAwvbldw6l1J3QVXox1RPZzHt+C/hm/\nKexnZlc2fnlBf5AlAmyFKybcKuke/I/sbWBBYDywGUWuNjN7T9JNeNq14casmKdxA3GIpI+B94E3\nzax4N1OO6S6zLBniGNxF8jtJV2Wv7w0sAxxYmNwQtB6JccA9QM6s7mSV/kNaD2+j8B6wEWbPlB2W\n1xfxz/jTwCqWs4rJPJnO29m4i/oo4K7BHPepRq1GaDWyGhJJQ/G6juPN7HRJedwvGkZoAGNmf5G0\nGp4uugMejJ0Lb5/xGH6zL3a1gX8u9sTTqUvSUc3sX5J2xWuDzsUVg1NmuNQq/WFa8WtmdrGkv+E7\nr+42w38EtisuVq0yb0+vBQ1AYmXgbuAYM37Y6vWURVoQz2bbHK8ju76C620eXMFhO+Bwy9lPi8dM\nH5umC+AP6Tvhn/mLLUn+24TVdww11QlJ+gTYzMx+LWlt4Ld4TcZrmZ/wTjObs7lLbRyhmNA4Cj4P\nx5rZ6a1eT6OI9773SKyB76qPMOPGVq+nBM96OwCXCLoW6KLUjexD89oSuBjf0U22nL1bdlyazop3\nTD0O+DGQtyR5u9zYgUwrWzm8hosJ/hr3bT5nZq9lr80L7dFUKiuWvApYFJdoudPMKgpuRmJCQzgM\n+A+xEw4AiXWB24ADzSjWIGw90pdwTba/A5tgVrY2R3kthGd9rg3sbTn7ZdlxaSpgKzzl+s/AhpYk\nzzZj6a2m1a0cLgfOkIv1bYFb+27WAdrlP/2/wGQze0zSLMC9krY3s5vLDY5WDr0ja/mwNf5gsgcw\nxcz+3tpVBa1GYgO8WHhvM37e6vXMhLQcrtAwDle7vr2C6024e/l7uLzTRMvZx2WnTNOVcIO2GPB1\nS5KS1g2dREuLVc3sVEmv4U8FhzOzKOX8eKpuyzGzN3DxS8zsv5KeABZv7ao6koVwHbYPcY24Y1q7\nnKDVSGwCXA/sZsZ9rV7PdKT5cF3Cr+G1PDsVdzidPjSvpfCMt8WALS1nj5Qdl6YLASfjOoQn4/19\nIu7TSzpWO07S/Lim2JesqNAsYkJBT8R7XzsSE3D16x3NqKqV1m943Gc/PO5zO3AiFZTTlddQPJ5z\nEp7RdqblrMSoKE1nwx/Cv4WXJ5xsSVI2RtSp9GtMSN4lsGbMbFo947MmapOB9XC3zgNmNr7MuHF4\nq9x18RTKy4C8mX1WZe7Z8ALMc4oNUBAEjUNiG9wTsk3biJBKm+BusneAzbHKMjrKaxwebvgfsL7l\nyigjeNxnG9xF9yywviVJ3FcaRDV33Mt1zGPA0DqvPQ6YADycraPUP+vKyffhBWRb41pLZ+FFtidm\nYybiwXHw3jF/wFOJHzWz9i2OC4IBjsROuGL0FmaUdV31K/5g+z1gZTzl+pZycR8A5TUrXut4GL4D\nmmK50gdbpekquEFbCDjEkuSeJq1+0FLRHZcVL3YzD16o+AxwC55ZshDuE10eOMbMflzXhTN/R/bv\nm4CRZrZx0Zhj8Q/TUgXS/ZPxLfYi5YoTJV2GtyXYt8q1wx0XVCXe++pI7IHf8Dc3a3Hbdm8tcgIu\noPs94Nzi9trTh3riwRb4w+yLeLO5V0vGpenCeL+qbfD7zaWWJBV14QYLzXDH1VondDXwiZkdVOa1\nKcAcZvbVXi+ishF6APg/M9u94NiS+C5tazO7o2j8+nga+ZN4ijbA5WZ2QdG4MEJBVeK9r4zEvvgN\n+ktmlFWLjqkyAAAgAElEQVQX6KeFDAX2xZMDfg4cjycnlR/urrez8Tb2rmKQK5Ll8bjP1/Fkm6uB\nUyxJ3mvK+gcgrawT2g6voi/HT/H4SzMYAzNn2pjZtEz+ZQxwR9FrvyG6xQZB05A4CC/RGG/Wwn5M\nXrNyLt6GZEvMHqs4NK/58d1Mt3LHRcWJB1ncZzu8HfdTwBcsSaLfVD9QqxH6F66yXE4g8os0r1h1\nBJ6MUMy72WtBEPQTEkfiu4TErEWK5S6kfCawOr5bKatyDaC8ZsHjxCfgpQRjLWcl4qFK09XwuM/8\nwIGWJO2TYj4IqNUIXQScmKU938aMmNC2uPzFd5uzvOaSKSZ0E8oJQVABiW/iKc8bmVFXJmyDFjAP\nvgPbD3ep7VEp7gOgvCZk414FNracPVUyJk0XwXdGW+FahJdH3KeUTCkhadb8tRardkl6F2/HWyjn\n/wYwyczObcbi8B3PvGWOj8heawRhfIKgAhLCs8d2ww3Q6/28gKG4eO4puCDqSlgV9eq8SruWlo/7\nfANPeroSGGNJ8n5T1t8BZPfHNDNGueqj66fmpnZmdp6kC4AlgEVwA/SqmTWrKybAc8DYwgOZPtwc\n2Wt9ImR7gqAymQH6DvAV3ACVLfZs4gI2xOM+HwNbY+UVDACU10j8Brk77pnZ1nI2U9fSLO7zFXyH\n9DSwriVJxWZ0wcy0VLanYBGf4plpLzd8JeW5C5gsaa7uFG1gF/xD+au+Th4Cpk72hFNWoDFjXTP7\nfYOutTT+ZHuLmbU2tTeoSGaAzsJ7SY03o/8asUmj8JKQtfG4z41V4j7D8NYyJ+EJUuMsZ2+VjEvT\nsbhBWxI41JLk7iatvmNptYApkj6H+04/B8xe/LqZ1aUfJmk4sGX24+eAuSXtmP18p5l9gus4HQHc\nLOl0YFn8aefsAqPUa2InVMJ1UFZ48i9ljvWWpfEbxkvQ4vqSoCwSQ/Ai1DWBjc0a5vru6cJz4wWk\nB+IGYy/8PlB+eF5fxnc1f8M7nJYoYitN58PvGXviu7oLQ+etd7R0JyRpO1yccAielFC4zRWudlCv\niOXCML3XSPdTzo3Zv0cB07LOnZvgfxA/w+NAZ+Ppln0mdkIlPGZm1zVjYklzFxUX1/Rpln/q5zCz\nj5qxrmBmJIYCU/Ai9E3NKNtnp8EXHQLshRuJ+4GVmdEqpnR4XmPwXdoYXBH7Z2XiPt01RKfg2nEr\nWJKE0nsfaNZOCDPr8QvXS7odLyit6Zx2/iLr2lnLuFavtZ/+PxK8uPeoHsZtiKfpv4e7RB8F9i0z\nLgX+ij9M3IRreH2GKxl/VuZratE6voYLSj6Dp//nstfXxvtFvQB8hNeIPAhsW2YNV2VzzYM3JXsT\n+CQbv3a89+V+ZxsGdg3YVLC5+uW68EWDRwwesh7eF7oYQRfn0MVbdHE0XcxadtzUqeszdepjTJ36\nIFOnrtHq/9dO+qr13lnPV63uuCWAw83snRrHDwhiJ1TCnJIWKDr2LzP7p6Sv4JJNr+PSKB/iGVOX\nSVrGzE4oOMfw1uC/wm/6x+Ip/Q/gQePj8KftX2fjiwPeR+I1G5eQJcBkx7cFlsN35a8AC+AG62ZJ\ne1h56ai78d17Pht/FHCnpFHWAJdupyAxC95ldD5gSzPK9tBp4AWXwuM+6+Gq1D/uvsuVDM1P74Sa\nwz+DK1iutH+V0nTxbM4N8EzeH1uSdGabgBbQ6p3QvcChrbbC/W3NaxnTCV/M2IGU+7oOd8O+gu9o\nFik4bxbcyPwPGF1wPM3OPbnKtfaq8to/gAXKvD5HmWPD8UzJp4uOX5XNdUHR8R2z4wfEe9/9u9ps\nYLeC/Qxs9qZeD+Y0OMXgbYOclXlPZ1pbF5vSxZN08Uu6WLnsmKlTZ2fq1BOYOvVtpk79NlOnztnq\n/9NO/ar13lnPV607oW8A10n6CO+1XqJiYFa+++BgQ2naL09eliTN0DWbgleWF/IGHqBeAk8Ima7N\nZd448AzgVlzo8azCJeI7pt5wjVlpZXvhZyzr7jocjy1NBQ4syqLsplhJfWr2fXQv19ZRSAzHpbc+\nAXY2myne28gLCX8AOAt/cFkVKxUOnT48r89nY8fh9Ty3lYn7CN8dn433DlvLkqQ1Sg5Br6nVCD2R\nfb+iwuu9aeXQcprhjmuScegvXjSzklRtSbtk/3y6zDndApajio6/ZWa9DWqX1eyStBBe4b4NsGDR\ny4a7koqN0Ew3JTN7O8vwmb+Xa+sYJObEHyDeAvYyozlqAdIKeE+w+YE9MavY+E750k6olivthKo0\nXQE4D69Z3M+S5P5mLD2YQatTtCu2RRjIWKRoN5O+7IxLzs2y5O7Bs7bOBR4B3gc+xT+fu1NGvNbM\nKu1MB/LDQp+RmBsXAH4J2M+Mxhede4uFHPBVPCb3A8zKGros7lPYCXUFy5V2QlWajsjm2hXPfLs4\npHb6B2tliraZXdXwKwcDie46oRXLvDYu+16rG6S37sqVs6+8meULX5B0QC/nHJRIzIcXgj8BHGw2\nve1Joy4wBDc8p+J1ZytgpYkE04fnp3dCfRvY3HKlnVCzlOv9cQN0MzDOkqT/CmiDplGXYoKkxfBs\nlpF4kPphM+tfLakGEtlxNfMYMA3YR9IZZv6EKmkWvEX7Z7iwbS10u8vqdYd1P6nPtNuRtCIuwV/O\nuEVmVBESI/Ed5W+AIz1XoKEXWB2v6xsGbEsVpQ3lSzuhFsd9AJSmGwLn4zvfL1uSVGzXHTSPlrrj\n5CKCF+BPIoU3gc8kXQIcZlbaGrfdCXdcbZjZZ5IOw9Nj/5C95//EJZTWAb5jZsWqCpX27U/j6d2H\nZH2h3gfeNLOpFcZ380x27jFZUsILeLr2AfgT/RplzhnULrdiJBbE+3PdDXyzoQbIFfa/gycKHA9c\nSYV7gvIzdUI9E9jVcqWK2ErTJbPX18Mfdm6MlOvW0VJ3HL4F3gev97gRr+tYGNgZ72r4Nh5MDDoU\nM7sjU684Ab8hzIobholmdmXxcCrsQszsX5K6m4udC8yGp3RPLTi33HmfSdoSf3L+GjAn3kF3L2BV\nvL9MTWsYbGQ6cNvhWWTXALmGGSB/QD0Av0dcD4zFrKzMj/LTlRFOBe4EVrRcaSdUpelwXIHlCPzh\ndx9Lksi+7VBqbe89Dfi+mZ1Z5rVJwBFmtmQT1tcUor130BOd8t5LrIwb+wWBr5tVFaqtd/L1cSPx\nAXA4Zk9UHJrXqsCFeBbtoZazR0vGeMr1jvju5w/AJEuSVxq23qDPtLK990JUFpt8Et8VDTgiJhR0\nKhLz416KHfFdyiUNS8GWFsXTpzfGd8XXU+FpNnO9nYxnsx0PXGG5Ujed0nRlPOV6fnzn05N7Nuhn\nWp2i/SIu0XJPmdd2AZ5v2Ir6kYgJBZ1GJr9zEO4evx4Ya0Zj5LY8EeVwXHbpctz19mHZoXkJV64+\nHRcfHmc5e7tkXJrOj6da74inZ18SKdftSatjQqcA10taEq+ofxPfHe2M9xvZteErC4KgLiQ2xXcT\nr+M9gMoVF/d28k3xDLVXgfUxq/jgqbxWxF1vc+LN5Uoy5JSmw/CWDTk8zry8JUlHaVMGtVFrndCN\nkt7Dt9Xn4Zph/8VVlL9sZvc2b4lBEFRDYllc4mYlXKD19gYmHiyZzb0mLt91WxXX2zzMKE7NAZdY\nrrTzstJ0PG7Q3gI2sSQp6QMUDB7qae99D3BPlq69APAPa25r7yAIqpCpHhyHl058D9jVjJJU515O\nPjteu/MN3GBUbDCXud52ydZwN652UK676XDcoG2JG8ubI+U6qLVOaB5gLjN7PTM8bxa8thjwoVXw\nDfcnkn4FzIvXMr0E7GMV0kWz8V1EYkIwwMg6n+6JpzrfB6xsRuOKxqWtcI/Hn4A1MHu54tC8xuIZ\ncvMDO1vOHio7zrXersdrvVa2JHm/YesN+oVmJSbUmqL9E+A9M9u/zGtTgHnNrOVxocLunZLOAv5j\nZseWGRcp2kFV2vW9l1gH35kAHGHG7xo4+efxdO7RwBGY3V1xaF5z4ckP++Ju+ostV6oLl6VdH4DX\nhX0TuDJ2PwOXVqZobwAcXOG1nwM/aMxy+kaBARqCN1UbkFl7QVCMxGLAacAmeBO4HzVM802aE3fr\nHYg3hdsOs7ItHTLX2w544WsKrFSu4BRAaToSuBRYFtjAkuS5hqw36ChqNULz4u2Uy/FvYERjltN3\nJP0cD6K+iFdcB8GARWJ2PC5zNN5pdnkzGuP69nzbnfBYzq+BVTB7reLwvJbDWzJ8DtjTclVaMqTp\nBnin1luAPSxJGhOrCjqOEun7CvwZ2KrCaxOYobJcE5JGS5oi6QlJn0oqW5gmaZyk+yV9JOk1Sfls\nl1MRM9sC7zHyIO7XDoIBh4QktsOlkdYG1jbjuAYaoBWA+/EC0j0x26OSAVJecyiv7wAPAb8AVqtk\ngJSmQ5WmJ+GlHIdakhwZBiioRq07ofOBH0j6D3Al8DdgMVzD6zAqu+oqMQ43Xg9nayhVzpVG4EHX\np4CtcT/1WbjhPDEbMzG7PsAhZvYwTNcZuwYPhPaJbh9oEPQXEivisZlFgAPMuK+Bk8+LF4XuSc89\nfoQ3EDwXN0ArW66yar7SdAl89/MpsLolyYBV2A/6j5oSEwAknYD7jWcvOPwJcIqZnVbXRbOob/bv\nm4CRZrZx0Zhj8RTRpbpbNkuajP8BLVKcjSdpPmC2gjYDJwGjzGyfMtdveHAtCPpK1mbhZGYIA/+g\ngVI7Q/CHxu/i4qHH9dDjZ1n84XMZXOutquac0nQ7PDZ8LnCGJUmUb3QgrUxMwMy+LekCXFZ9flw5\n+2Eze6/ei1bpdlnIBODubgOUcQMuA7IR3hWykBHADZJmy35+kogJBQMAiWL1gLFmlEjc9OECG+Ne\nhE+AbXro8TMcz2I7jCxJwXLlkxRgptqfCcA2liS/bdi6g0FBXU3tMoNzV5PWUswYmNkNYWbTsh40\nYygyQmb2V9x3HgQDBomN8dilqwcYjVMPkMbiitRj8Yy6myqpHQAor63w3c+jeNzn1arTz1z7s2rU\n/gS9oWYjJGkVPIi5JrA4sK6ZPSbpu8CvzazRxmkEUG6X9S5tlI0XBL1BYhk8K21VPPPt1gZK7SyE\nu613wgtad8Ds3xWH5zUKd6ONBQ6ynJUTKp4xPmp/ggZSq2LCBOB2PDh5Ne426ObfuLJuf+2QGkam\nmNBNKCcETUdiLrw55IF4rc3uDZTaGQ4ciRu1a4HlsVLl6unD89OleY7M1rKz5SobK4jan8FIppSQ\nNGv+WndCpwJXmdn+koYxsxH6Iy4d32jexeuTihmRvdYowvgETSeT2tkD/1v6JbCKGRVrcuqcfAiw\nO95e+xFgPcxerHpKXl/G5XaeAta0XGVpnunnRO3PoCS7P6aZMcpVH10/tRqh5fEnpnJ8AIxszHJm\n4jncPTAdSUsAc2Sv9ZnoJxT0BxJr43GfocBOZjzcwMk3xBMDDK/3+XXV4XktCZyDuwGPsJzd2eMl\n0nQo7oo/BNjPkqQ4KSgYBDSrn1Ctxapv4dvvcowDpjVmOTNxF/BlSXMVHNsF+Bj4VSMuIKmrWaJ8\nQSCxqMRVwK3AFGDdhhkgaTmkW4BrcFfautUMkPKaXXkdCzwGPAGsWKMBWgLfuSV47U8YoEFKqzur\n/hg4WdLTMOOPSNIYPDB5RT0Xlfuut8x+/Bwwt6Qds5/vNJeM/wGeYn2zpNNxI5gDzi5K2+41sRMK\nmoHEbHicZTJwGS6180GDJl8AOAl3v50B7IZZRbdYgdbbGbgq9jqWs5oUTgpqf84DTo/an8FNs3ZC\ntapozw7cBGwBvIFXcr+Wfb8b2N4qCB5WmG9pvNUCzFBLUPbvUWY2LRvXLRO/Hh4HugzoqrHOqNr1\nu8/PEzGhoEFICFf3OAtPWz7ajD83aPLZ8QSgY/B6uTxW2rNnplPyWh13vc0HfKOngtPp581c+7O7\nJUnj3IfBgCXbCU2Fxhar1qyYkC1iE2BTvKndO8B9A7GraigmBI1GYgX8hr84cKQZVdOc65i4u2Hc\nqbgb7ZhqrbUBlNeiePr0lviu6fJyHU7Lnjuj9ucZ4EBLkrqL0YPOpaWKCdmF78dFDzuCaGoX9JVM\naqcL2A04BbjYjP82aPL18R3JLMA+9PA5zVKuuxW3LwfGWM5qKiDNan/2xzPsvgVcEbU/QSEtbWpX\ntJA5gYm4asGbwNVm9koT1tY0YicU9JUs7rMvboBuBk404x8NmnxZXJ5qbTwr7UeYVewdlMV9dsTj\nPo8Dk2uN+wAoTUfgtT+jgd0sSZ7tw+qDDqZfd0JZZ9KvmNlyBcfmxusQPo+74+YFjpK0tpm90KhF\n9RexEwrqRWIJvNB0PzzQv5kZf2rQ5COBE4C98Iy3vTD7uOopea2Bqx3MDexrOSvbFqXi+Wn6ReBH\neAbfnlH7E1SiFdlx4/EPZyGTcAO0n5ldIWlBXN/tJFwafkAR2XFBLWQJBxsDh+Kpyj8CxpvRmB2D\nNGs297H4rmoFMjX4iqfktRjuOtscb21yZa1xHyip/dnfkuRnvVx9MEhoVnZcNSO0NL7rKWQH4Fkz\nuyJb1FuSvofLzgdBRyExD74rOQT4DLgQ+FqDO5tuj7vengcSzJ6peoqrXB+Fx34uw+M+daV/K00X\nxw3pp8AaliSNUW4Igl5QzQgNgxmaVpLmxxUMLioa9wqeqh0EHUGW6XYosCu+0z8YeKBhAqN+kXXw\npIO5gIMwq9q4Lov77ITHfR4B1racvVTtnLLzpOm2eOFs1P4EbUE1I/Qi7pLrzobbEq/lubto3EJ4\nfGjAETGhoBuJWfAuoofiMlWXACs1TN9txoVG4enWX8TdaNdg1d1oymtNPO4zJ7C35er/vCpN58cb\n2m0GbBu1P0G9tCIm9H3gUnk74DfxQrm/Qkn9w5dwEcQBR8SEAslbaGdfL+Eut1vMqLn4usYLzYd3\nJp6I70ImYvZR1VM87vNd4Mt4wsJV9cR9YHrs5wC8MPsnwGpR+xP0hn6PCZnZVZIWxTsszotrTh1a\nqIwg71uyLf4BD4IBQZZosD6+69kc72a6hRlPNOFis+Aq8yfg7VBWxOxvVU/xuM/ReNznEnoR9wFQ\nmm6IN6l7D9jUkqTxv18Q9JG664Q6gagTGpxIzIlrrh0GDMd3PVeblW2e2NeLCXfvnYF7ECZjVtUI\nZHGfnbNzfg8cYzn7a92XdtHRM3BDOwn4SRSeBo2g5YoJnUbEhAYHEp/HM9z2Ah7EhUXvM6NiAWgf\nL7gm3jV1AeAIzH7R4yl5rYXHfYYDX7WcPVD3ZdO0u0ndN3DNxYmWJFXrjIKgVtpGMaETiJ1Q5yMx\nFE+mORRYDZexmWLGy0286JJ47c4muOL7lZj9r+opeX0Oj/t8Ca/buaYXcZ/uXdfZuGLCJEuSundQ\nQdATsRMKgh6QWAAP/h+MJ9RcCGzTsBba5S86D663dmB2vTGYVa0lUl5z4HGfI/GU6TGWq35O2XnS\ndCye6PA54ABLkqqp3kHQboQRCjoCibXwWM82uATNTmb8ockXHYaLfubwJoyrYPZ/VU/xuM+uwGnA\nb4E1ammtXTJPms6XXXdPXDH7IkuSxginBkE/EkYoGLBIzI63OTgUWBC4GO/h0xgh0coXHoIrHZwM\n/A2YgNnjPZ6W19p43Gc2YE/LVW/FXXaONB0C7IO7/W4HxlmSVO0rFATtTEcaIUkXAgebWdX25ZGY\nMDCRWBp3t+0DPIobg7vMaG71vxuf7fAdyL/xJIC76CGwqrwWx+M+mzAj7lN3UoTSdD085fq/wJaW\nJI/WO0cQ9JZWt/ceMEjaAK8s7zHjIopVBw4SQ/Dg/aHAF4BrgPXNeLEfLj4Er4fLAf/BhUZ/XoPx\nmQM3VF/Hd2ljLFd/a3ql6aK4vtwmeGfV6yLlOuhvWiFgOuCQNBsuibItno4bDHAk5gP2xlOsP8ID\n/7uZUVVtoEEXH4LHmHLA//BdzJ01GB/hTe5OAx4CVrdc/T23lKaz4QbsGLzfz/KWJI0RTw2CNqGj\njBDeUuIyM/tHMyx20H9IrILvenbCg/77AA81VES08sW7U55z+I76ROCOnowPgPJaB4/7DAN2t5w9\n2KslpOkW2TzPA+tZkjR/xxcELaAlRkjSaLxgcD1gBeABMxtfZtw4XMNuXVx65DIgb2W6TEpaGVjb\nzI5XWKABicSseMD/MLyVyA+AsWa80U8LELA13i2V7PvtNRqfJfBdeIJrxF3by7jP54FzgOWAIy1J\nfl7vHEEwkGjVTmgcMAF4OFtDyR+5pBG4jP5T+I1hNC59PwR/MkXSRPyGBS7OOE7SXwvmeAlYy8ze\nbtpvEvQZic/hNTb7A8/iN+HbzKha6NnABQj4Cm50hmTfb6vR+MyJP1Adjrc5OaiXcZ+5cXfffnj8\nZ3tLksaKqAZBG9ISxQRJsuzCkm4CRprZxkVjjsWDukuZ+R+1pMn4DWIR66kYUPqsUnZcKCa0nkxE\ndCPc5bYJcB1wkRlVm7o1eBECtsI/U8OYYXx63MEoryHMiPs8CHyrl3GfIcAe2Tz3AsdaklQVOA2C\nVtExiglWm+WbANzdbYAybsCfEjcC7ujpMr1cXtBEJOYGvoonGgzBEw0mmlG3SnQfFiFc0qcLmAVX\ngb+1FuMD0/v7nJ+du4vl7KFeLSNNC+fZwZLkt72ZJwgGMu2cmDAGd8dNx8ymSfo4e62qETKzoU1c\nW1AnEmPxXc/uwC9x91XaL4kGMxYhYAvc+MyGG59b6jA+C+H1PlvicZ+rexn3WQgvNt0Kd8FdZUnS\nHDHVIGhz2tkIjYCyEvvvZq8FbY7EMDyedyiegHIpsLIZVaVtmrAQ4TvrLlylOg/cXIfxmQX/HY7H\n65OWt5y9X/cy0rR0niSpe54g6CTa2Qg1nUwxoZtQTmgQEgvjAfYDgWm4y+2nDe9W2vNChDet68IL\nmPPAT2s1PgDKa1PcZfYqsKHl7NleLSVNZ54nSXo1TxD0N5lSQtKs+dvZCL2Ld3QtZkT2WqMI49MA\nskSD9fAn/S2Am4CtzfhjCxYjvCV2FzA3bnxuqtP4jMKzMVfB+/P8zHL1Z/EoTbvnWTWb5/ZQOwgG\nEtn9Mc2MUa7R87ezEXoOGFt4QNISwBzZa30mZHv6jsQceJbYYfgN/0LgMLOGPijUuhgBm+HGZ15m\nGJ+aNeWylOtv4dp0Z+MFp3W3gVCazgl8E0/A8HmSpHntJIKgyQxG2Z67gMmS5irIkNsF+Bj4VSMu\nEAKmvUdiNH6j/hpe73UscE/TupVWX4yATXGjMyL7/pM6jU93a+0zgV8Dq1queluGsvN4g7nueR4E\nVrUk6d8YWBA0gY7qrCppOJ5hBN7Ya25mVKnfaWafSJoPeAYvVj0dWBZ3a5xjZif18fpRJ9QLsm6l\nm+Mut7WAK4AfmNGaLp5ufDbBPzvz42raN9ZjfACU1yp4vGYe4IjetFgAUJp2zzMvcLglSa/mCYJ2\npRn3zlYZoaWBl7Ifuxeg7N+jzGxaNm4scAEea3gXl+3pqrHOqNr1u8/PEzuhmpHYGleruBC4wYxP\nWrQQARvjxmdB3Pjc0AvjMz9wCrAD3bqDdbbWBlCadhvAHXGf+aWWJM1tKxEE/Uy2E5oKHWCEWk3s\nhHqHhPq1rqd0AQLG48ZnYfzGf30vjM8w4IBsnhuAnOXsnbqXk6ZDC+a5EchZktQ9TxAMFDpGMaFd\niJhQfbTYAHUbn0Xx3cuPMatbW055bYS7zN4GNrGcPVn3HB73+RIe93kH+JIlyRP1zhMEA4mOigm1\nmtgJDSD8g98FfA43Ptf10vgsiRuNdfE45E97mXK9Fq7ztjiumnBzpFwHg4XYCTWY2Am1MdJGuPFZ\nAjc+P+ql8RnOjO6mFwD7WM4+rnueNF0O+DawPh5LvNKS5L/1zhMEA5XYCTWQ2Am1MdKGuPFZCjc+\n1/bS+AjvsHs28AgwqZcq14viyQY7ZHOdZ0lStxELgk4gdkJB5yJtgBufpfEdx7WY9WqnobxWAM4D\nFgEmWs5+WfccaTov3lb7IDwVfXlLkuhLFQQNZlAboXDHtQFSt3trGdz4/LAPxmc+3JDtjmfO/cBy\n9e2ilKaz43VQ38SV2le1JHm1N+sJgk4i3HENJNxxbYD0Bdz4jMaNzzV9MD5DgX1x992twImWs7fq\nmsPTrb+arelx4HhLkqd7s54g6FTCHRcMfKT18Bt9d6D/Gsx6ra6tvL4AfB/4BJhgOXu8rvM93for\neJ+gd4HdLEl61aQuCIL6CSMU9A9ufLrwhoTfAa7uo/FZDJdzGo/Hbn5cb8q10nT9bI55cdHSOyPd\nOgj6l0FthCIm1A9I6+LGZyxufK7qo/GZDTgSmAxcgjeY+2f1s4rmSNMV8Z3PKrgM0Y9CZicIqhMx\noQYSMaF+QFoHNz4r4Mbnyj4aH+Git+fgwrZHWc7+Utccabok7grcAi84vTjaKwRB7URMKGh/pLVx\n49O929gWs3/3acq8xuDGZxngcMvZL+o63wVGjwP2Bi4Clou22kHQHoQRChqDtBZufFbGjc92DTA+\n8+Dusr2BU4ELLFf7biprLPd1vKPpT4AVLUn+1pc1BUHQWMIIBX1DWhM3PqvihmL7BhifIXi69KnA\nL4CVLGdv1Hx+ms4CTMQN2K+B9SxJ/tyXNQVB0BzCCAW9Q1oDl7NZHY+v7IjV3wa7ZNq81sJTrgVs\nazn7fc3nerr1jngM6hVga0uSR/u6piAImkfHGSFJLwMfAd1um93M7LkKY7uI7Lj6kFbHjc+auPHZ\nuUHGZ2HcjTcBj99cYzmruVW40nSTbD1DgEMtSe7t65qCIJhBZMfViKS/Aht1d2etMCay43qDNAHv\nbns6cEmDjM8swOHAscBVwCmWsw9qPj9NV8ONz7LA8cBPLElqNl5BENROZMfVThiX5nA/MBqzhrT1\nVl6b4UKjLwMbWK78jrXsuWm6LK64sBEu13NZtFYIgoFHp+6EutNv7wC6rKgVQOyEWovyWgZvi7Ai\nXhVF16IAAAsDSURBVHh6Z61qB0rThfGEg12Bc4FzLUnqKlYNgqB3dNROSNJovOp9Pbyg8QEzG19m\n3Dg8UL0u8B7uDsqbVYwXrG9mr0uaE/gh3tDstCb8CkGdKK+5cLfbgcD3gF0sV1smndJ0Hvy9PBS4\nGm+t8I9mrTUIgv6hle64cXgQ+uFsHSVPwpJGAPcBTwFb44rLZ+HB5xOzMROBw7JTDjaz3wKY2UeS\nLsdveEELydQOdgXOAH4FrGI5e62mc9N0Nrynz3F4uvbqliR1N6cLgqA9aZk7TpIsu7ikm4CRZrZx\n0Zhj8affpcxcH0zSZLwuZREz+7Bo/BzAMDP7QNIw4FLg/8zsxKJx4Y7rJ5TXasD5wJy42sFvajrP\nWyt09wV6CjjOkuTJpi00CIIe6Sh3nNVm/SYAd3cboIwb8OysjfCYTyGLAD+VNAQYCjyE14wE/Yzy\nWgBPHNgW37VeYTnrUSQ0q/WZgLtQ/wnsZUny62auNQiC1tHu2XFjcHfcdMxsmqSPs9fuKHrtJWC1\n/lteUIzyGoa7z04CfgyMtZy9W9O5abou/oCxAO5+uz1aKwRBZ9PuRmgEnoxQzLvZa0EbobzG4663\nvwMbW86equm8NB2LF6quiRuva6K1QhAMDtrdCDWVTDGhm1BO6CXKayk8220t4CjgllpSrpWmi+Px\nva3xpIXdLUkaUoMUBEFjyJQSkmbN3+5G6F2862UxI7LXGkEYn16ivIbjXU2PwItO97Jcz4WsStOR\neCfTiXhjuuUsScrteIMgaDHZ/THNjFGu0fO3uxF6Du/IOR1JSwBzZK/1CTPr6uscg5Es5XoHfPfz\nO2A1y1WWSZp+XpoOxw3WJOBmYCVLktebudYgCBqDmaVS4xOK290I3QVMljRXQYbcLsDHeL1JnwgB\n0/pRXivicZ8FgL0t1/P/ndJ0GLAP/hT1W+CLliTPN3OdQRA0lo4TMJU0HG/XDHA0MDceHwC408w+\nkTQf3sr5KTxralm8WPUcMzupD9eOOqFekGm9XYu3yJ5iuZnlkErGe7r1dnjSwevAtyxJam7NEARB\ne9GMe2crjdDSwEvZj92LUPbvUd0q2JLGAhfg8j7v4rI9XTXWGVW6dve5eWInVDPKa3ZgLstZj3I5\nStMEr/WZHfgmcE+kWwfBwCXbCU2FDjFCrSR2Qs1DaboK3hF1eeAE4PporRAEnUFHKSa0AxETahxK\n01F4S4VNyZQSLEn+U/2sIAgGCh0XE2olsRNqHErTBfEdzx642vnZliQfVj8rCIKBSOyEGkzshHqP\n0nRuvDD1cOBHwDhLkr+3dlVBEDSL2Ak1kNgJ9R6l6azAAXgr7fuBkyxJXqp+VhAEnUDshIKWojRd\nGxclfQGYYEnyxxYvKQiCAc6gNkLhjqubV4GJliRpqxcSBEH/Eu64BhLuuCAIgvppxr1zSKMmCoIg\nCIJ6CSMUBEEQtIxBbYQkdTXLzxkEQdBJREyogURMKAiCoH4iJhQEQRB0FGGEgiAIgpYRRigIgiBo\nGWGEgiAIgpbRUUZI0pySrpL0nKSnJB3cw/jIjguCIKiByI6rAUk/AF42s9Oynxc0s7fKjIvsuCAI\ngjrpqPbejUbS3Liw5uJm9mkPY8MIBUEQ1EmkaFdnGeAt4HxJj0q6VdJSrV5UEARBUJmWGCFJoyVN\nkfSEpE8lTa0wbpyk+yV9JOk1SXlJldY8DFgRuNXM1gBuA65u0q8QBEEQNIBW7YTGAROAZ4HngRKf\noKQRwH3Ap8DWwMnA0UC+YMxESY9LegyYE3jfzO7NXr4BWKOZv0QQBEHQN1plhH5mZkua2S7AMxXG\nHATMBmxvZveb2RTcAB2VxX8ws8vNbDUzW93MHgCekLRmdv6XgCf+v727C5WiDuM4/v0hEYZSEpQl\nRtGLGRFBJaS9CEFRQXRhdFM3FuFFBF1ESKTeRBipKIkVYhEIRSJBrxIVQtFVpRSlFAqnLkQyAytJ\n4jxdzJxYt3XPnN2Z/c/M/j6wcHb+s7vP8DDz7Jz5z7MVb4cl4lmNzeXcWackRSiKzYa4G9gTEX90\nLHsLmA3cfobXrAI2S9oHPAmsHCpQq7PlqQOwgS1PHYDVR50nJiwCDnQuiIgJ4K987H8i4oeIWBYR\n10fE8og4OII4h1Lmt8JB32smryuybr91Bhmr6zfnsuNqY/7qmjtoXv6GzV2/8ZT7Xp1/3nse8HuP\n5cfzsaFNTTdMTSpvpvig7zWT1xVZt986g4z1Wi5p7bSBVKzM3A3zfnXOX11zB83L37C56zc+0+Vl\nqfOZkJmZtVydz4SOA+f2WD4vHxuYb1I1M6uHOp8JHQAWdy6QtBA4h65rRWZm1kx1LkIfAndJmtOx\n7EGyiQl704RkZmZlStUxYbakFZJWAAuAC6aeS5qdr/Yy8DewW9Idkh4D1gIbu6ZtVxXjNkm/SJqs\n+rOsPJIW5l02vs87qa9PHZPNjKS9kvblHVXeyW9ctwaRtLXosTNJA1NJlwKH8qdTASj/+7J8KjaS\nFgMvATeTXQfaDqwreJ/RsDHeQtYQ9UhE1PmM0TpImg9cHBFfSzoL+BjYEhG7E4dmBUmaGxEn8r83\nAKciYnXisKwgSbcCjwAPR8SsaddvSxftqkiadBFqLklbgJ8iYkvqWGxm8j6R24CDEbExdTw2PUln\nA58A9wNHixw7fXC11pJ0PtnOsCd1LDYzkj4AjpA1Jd6aOBwrbg2wPSJ+LfqCVhWhirpz2wiUnbv8\nG9kuYFMTOmc0Xdn5i4h7gPnA58DmisMfa2XlTtJ1wJKIeF0zuMO1zvcJDWKqO/eXZNvWrzv3d2Td\nua8ANpAV5GdHFql1Ky13kmYBO4GvImJT5ZEbVLDvRcSkpDeAN6sL2ygvd0uBayQd7njdIeCmiDh2\nxk+PiNY8yK9x5X/vAj7tsc5q4Bgwp2PZU8CfwNzu9wMmU2/XODzKzB3ZBJYdqbdpnB5l5Q84D7iw\nY3wN8Frq7Wvzo+zjZsd4oWNnq/4FFfmWT6NQd25J24EJICT9LOnVUoO105SQu9sAJC0j655+g7Lf\nmvpG0uOlB2ynKXHfmwe8K2m/pP3AVWS/I2YVKfO42f3WRT6/bf+OK2IR2WnlfyJiQtJUd+738mWP\nJojN+uuXu6uB9yPiC1p2rbNFpt33IuIwsCRFcNZXoeNm1/i007NhPHfWyrtzW2Wcu2Zz/pqrstyN\nYxEyM7OaGMciVFl3bqucc9dszl9zVZa7cSxC7s7dXM5dszl/zVVZ7saxCLk7d3M5d83m/DVXZblr\n1ey4vAP3vfnTBcDcvFM3ZDOnTpJ1536CrDv3euByRtid23pz7prN+Wuu1LlrVQPTJnTntt6cu2Zz\n/porde5aVYTMzKxZxvGakJmZ1YSLkJmZJeMiZGZmybgImZlZMi5CZmaWjIuQmZkl4yJkZmbJuAiZ\nmVkyLkJmiUlaJ2lS0kc9xnZJ+ixFXGaj4CJkVh93Srqxx3K3NbHWchEyq4ffgG+BZ1IHYjZKLkJm\n9RDAc8B9kq5NHYzZqLgImdVDAG8DP+KzIRsjLkJm9aC8Jf7zwAOSrkwdkNkouAiZ1ctOYAJYnToQ\ns1FwETKrkYj4B3gBeEjSJanjMauai5BZ/ewAjgJPk10rUtpwzKrjImRWMxFxCngRWAlclDgcs0q5\nCJnV0yvACWApvlnVWsxFyCy9oKvQRMRJYFOacMxGR9msUDMzs9HzmZCZmSXjImRmZsm4CJmZWTIu\nQmZmloyLkJmZJeMiZGZmybgImZlZMi5CZmaWjIuQmZkl8y8+DWPeWd3siQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10ada9310>"
]
},
"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",
" D = np.zeros((N,N))\n",
" uF = np.zeros((N,N),order='F')\n",
" DF = np.zeros((N,N),order='F')\n",
" uF[:,:] = np.random.rand(N,N)\n",
" y = %timeit -o centered_difference_2D_vectorized(u,D,dx);\n",
" times['numpy'].append(y.best)\n",
" y = %timeit -o centered_difference_2D_numba(u,D,dx);\n",
" times['numba'].append(y.best)\n",
" y = %timeit -o centered_difference_2D_cython(u,D,dx);\n",
" times['Cython'].append(y.best)\n",
" y = %timeit -o centered_diff_fort.centered_diff_fort(uF,DF,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');"
]
},
{
"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": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"numpy vectorized\n",
"1 loops, best of 3: 1.39 s per loop\n",
"numba\n",
"1 loops, best of 3: 728 ms per loop\n",
"Cython\n",
"1 loops, best of 3: 758 ms per loop\n",
"Fortran\n",
"1 loops, best of 3: 449 ms per loop\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGixJREFUeJzt3X28XFV97/HPl0QxaXg4WOVJHiLIQ2j13lfRQqmSJq8I\nEeFWykPLVV9gY8t9YVUU2nIVOIj2dQGBa324xgsKqJcSKGIpIg+BQ0EBUWvVQiqRJxulFkgwEBBN\nfvePtYbsTObMnjkzc+acle/79Tqvc2bttfas2XPmO3vWXnuPIgIzMyvLVsPugJmZ9Z/D3cysQA53\nM7MCOdzNzArkcDczK5DD3cysQLXhLmlvSUslfV/Sekm3d3MHkraS9G1JGyQdMfGumplZp2Z2UGce\nsBi4O9fvdmL8EmDX3M6T6s3MJkEnwzLXR8TuEXE8cH83K5c0AnwU+BCgCfTPzMwmoDbco7dTWM8F\n7gKW97AOMzPrUifDMhMi6bXAScBv4712M7NJNcjZMp8EPhkRDw3wPszMrIWB7LlL+mPgNYBnx5iZ\nDUHfw13SS4ALgPOBmZK2B7bNi+dI2iYi1lbqewaNmdkERUTLYW91c7xU0jXADhGxoE2d7YGn2qxm\nZUTsU6nvcDczm6Dxwn0QwzJrgT9g0zntOwNXAmcAt7VqNF4HpxJJoxExOux+lKKU7ekdlE1Nh9dy\nnenwv1n3f1cb7pJmsXHsfFdgG0nH5Ns3RMRzklYCYxGxJCLWA3c0rWPP/OcPIuK+LvpvNk0MO99H\n88+wTftcL0Yne+47Asvy343/4GX577nAY8AM6mfeDPu/38xsi1Eb7hHxCDXBHRFzO1jHjG46NkWN\nDbsDhRkbdgfKMX/YHSjN2LA70KuuDqgOpAN53KiEcTrbMqX/YX8wTeTX8iSpy05f8tfMrEAOdzOz\nAjnczcwK5HA3MyuQw93MrEAOdzOzAjnczcwK5HA3MyuQw93MrEAOdzOzAjnczcwK5HA3MyuQw93M\nrEAOdzOzAjnczcwK5HA3MyuQw93MrEAOdzOzAtWGu6S9JS2V9H1J6yXd3kGb10u6XNJDktZJWiHp\nLElb96fbZmbWTu0XZAPzgMXA3bl+J18WeRywB/Ax4EHgdcC5wGuBYybUUzMz61jtF2RLUuRKkq4B\ndoiIBTVtXh4RTzaVvRtYCuwRET+plPsLsm1a8xdkV/kLsidLz1+QHXXp37rNky2Kv5d/79Lt+szM\nrDuTeUD1YGAD8ONJvE8zsy3SpIS7pJ2ADwNXRMQTk3GfZmZbsoGHu6SXAsuAXwCnDvr+zMyss9ky\nEyZJwBXA/sAhEfF0m7qjlZtjETE2yL6ZmU03kuYD8zuq283x0k5ny1TqfwJYAiyKiG+OU8ezZWxa\n82yZKs+WmSx12TmwPXdJZwCnAMeOF+xmZjYYteEuaRZwRL65K7CNpMaJSDdExHOSVpKGUpbkNieQ\nTmC6DPippIMqq1zpg6pmZoPVyZ77jqQDorDxs+ey/Pdc4DFgBpsenF2Ul5+YfxoCOIk0Dm9mZgPS\n1Zj7QDrgMXeb5jzmXuUx98nS8xmqZmY2/TjczcwK5HA3MyuQw93MrEAOdzOzAjnczcwK5HA3MyuQ\nw93MrEAOdzOzAjnczcwK5HA3MyuQw93MrEAOdzOzAjnczcwK5HA3MyuQw93MrEAOdzOzAjnczcwK\n5HA3MytQbbhL2lvSUknfl7Re0u2drFjSdpK+IOkpSWskfUnSDr132czM6szsoM48YDFwd67f6TcB\nLwP2Bv40tzkPuA54U/fdNDOzbiiifVZLUuRKkq4BdoiIBTVtDga+AbwpIu7KZa8H7gUWRcTySt22\n3+BtNtWl/+FO93lKJ7+WJ0lddtYOy0Rd+re2GHi8Eex5PfcBD+dlZmY2QIM6oLofsKJF+QN5mZmZ\nDdCgwn0EWNOifE1eZmZmAzSMqZAenDQzG7BOZstMxFPAK1qUjwCrWzWQNFq5ORYRY/3vlpnZ9CVp\nPjC/k7qDCvcVwBtblO8HXNuqQUSMDqgvZmZFyDu9YwCSzm5Xd1DDMjcCO0k6pFEg6UBgbl5mZmYD\nVLvnLmkWcES+uSuwjaRj8u0bIuI5SStJQylLACLiHkk3A1dIOo2NJzHdGRG39f1RmJnZJjoZltmR\ndLYpbDwYuiz/PRd4DJjB5p8CjgcuBj6fl10PvLfH/pqZWQdqz1AdeAd8hqpNcz5DtcpnqE6Wns9Q\nNTOz6cfhbmZWIIe7mVmBHO5mZgVyuJuZFcjhbmZWIIe7mVmBHO5mZgVyuJuZFcjhbmZWIIe7mVmB\nHO5mZgVyuJuZFcjhbmZWIIe7mVmBHO5mZgVyuJuZFcjhbmZWIIe7mVmBasNd0jxJyyU9K2mVpHMk\nddLuDZJukfRk/rlF0hv6020zM2unbUhLGgFuBdYDRwEfAT4InFPTbo/cTsDbgXcAM4FbJO3ee7fN\nzKydmTXLTwa2Bo6OiGeA5ZK2BUYlnR8Ra8dpdwQwG3hbo46kbwJPAIuBpX3pvZmZtVQ3vLIYuCkH\ne8NVwCzg0DbtBPwaWFcpezaXaQL9NDOzLtSF+77AimpBRDxGCu1927S7GngauFDSKyS9ErgYeCov\nMzOzAaoblhkB1rQoX52XtRQRP5d0GPA14L25+GfAYRHx5EQ6amZmnRvIVEhJewI3APcCh5OGd74D\nfE3SboO4TzMz26huz301sF2L8pG8bDwfBH4JHBMR6wEk3QY8CJwGvK+5gaTRys2xiBir6ZuZ2RZF\n0nxgfid168J9BbB/08p3I82EWdGyRbIXcH8j2AEi4leS7gde3apBRIx20mEzsy1V3ukdA5B0dru6\ndcMyNwKHSZpTKTuedED1jjbtHgZ+S9KLbx6StgZ+C3ik5j7NzKxHdeH+WdLwyrWSFkr6M+Bs4KLq\n9EhJKyVdUmn3OWAX4DpJb5H0VuA6YMe8zMzMBqhtuEfEGmAhMAO4nhzs+XfVjOq6IuJfgEXAHOCL\nwOXAy4BFEfGDfnXezMxaU0QMtwNSAESET26yaSn9Dw/3dTR1yK/lSVKXnb4qpJlZgepmy0wJjXco\nS7xnZGZ1pkW4J873xLluZvU8LGNmViCHu5lZgRzuZmYFmkZj7tYvPkC9KR+gthI53LdYzvfEuW5l\n8rCMmVmBHO5mZgVyuJuZFcjhbmZWIIe7mVmBHO5mZgVyuJuZFcjhbmZWIIe7mVmBHO5mZgVyuJuZ\nFag23CXNk7Rc0rOSVkk6R1JHbwqSjpZ0n6R1kp6QdKOk2b1328zM2mkb0pJGgFuB9cBRwEeADwLn\n1K1Y0hLgy8ANwOHAEuBH+GJlZmYDp4jxrw4o6QzgNGCPiHgml50OjAI7RcTacdr9JvAw8P6IuLRt\nB2q+wXtjHV/FMOn92+W9Pau8Pfur9+1pnanLzrrhlcXATY1gz64CZgGHtml3HLABuLzzrpqZWb/U\nhfu+wIpqQUQ8BqzLy8bzu6QhmHdL+ndJL0i6R9LBPfXWzMw6UhfuI8CaFuWr87Lx7EQK/w8BpwNH\nAs8CX5f0ygn008zMujCoqZAC5gDviogrI+Im4A9JB2ZPGdB9mplZVjdzZTWwXYvykbysXbsNwFij\nICLWSvoOcECrBpJGKzfHImKsVT0zsy2VpPnA/E7q1oX7CmD/ppXvBsymaSy+yQOkTwXNR3G3IoX+\nZiJitKYvZmZbtLzTOwYg6ex2deuGZW4EDpM0p1J2POmA6h1t2l2ffy9oFEjaDvgd4Hs192lmZj2q\nm+e+PXA/8EPgPGAv4ELg4og4q1JvJWkoZUml7CukWTN/DTwJ/CWwH7BPRDxdqed57l3xvOz+8vbs\nL89znyw9zXOPiDXAQmAGaW/8bOCi/LtqRot1vR24Lte/GvglsKAa7GZmNhht99wnpQPec++S9zT7\ny9uzv7znPll6PUPVzMymIYe7mVmBHO5mZgVyuJuZFcjhbmZWIIe7mVmBHO5mZgVyuJuZFcjhbmZW\nIIe7mVmBHO5mZgVyuJuZFcjhbmZWIIe7mVmBHO5mZgVyuJuZFcjhbmZWIIe7mVmBZg67A2ZmVY2v\nj7P2Xz9ap3bPXdI8ScslPStplaRzJHW8xy9pK0nflrRB0hET7aiZbUnCPz1qu+cuaQS4FfghcBSw\nN3Ah6U3hzA7vYwmwK/3qsZmZ1arbAz8Z2Bo4OiKWR8RS4BzgA5K2qVt5fnP4KPAhwN+IbmY2SerC\nfTFwU0Q8Uym7CpgFHNrB+s8F7gKWT6x7ZmY2EXXhvi+woloQEY8B6/KycUl6LXAScBreazczm1R1\n4T4CrGlRvjova+eTwCcj4qGJdMzMzCZuIFMhJf0x8BrAs2PMzIagLtxXA9u1KB/JyzYj6SXABcD5\nwExJ2wPb5sVzJG0TEWtbtBut3ByLiLGavpmZbVEkzQfmd1Q3YvzZiZLuAFZFxAmVst2AR4EjI+KG\nFm22B55qc58rI2KfSv2A9pP1Ux3PokzU04kN4O25KW/P/vL27J/227IuO+v23G8ETpc0pzJj5njS\nAdU7xmmzFvgDNn12dgauBM4Abqu5TzMz61Hdnvv2wP2kk5jOA/YincR0cUScVam3kjSUsmSc9ewJ\nPAS8NSK+1rTMe+5d8Z5Rf3l79pe3Z/8McM89ItZIWgh8CrieNM5+ETDaVHUG9TNv/GyZmU2Stnvu\nk9IB77l3yXtG/eXt2V/env3T2567L/lrZlYgh7uZWYEc7mZmBXK4m5kVyOFuZlYgh7uZWYEc7mZm\nBXK4m5kVyOFuZlYgh7uZWYEc7mZmBXK4m5kVyOFuZlYgh7uZWYEc7mZmBXK4m5kVyOFuZlYgh7uZ\nWYEc7mZmBeoo3CXNk7Rc0rOSVkk6R1LbtpJeL+lySQ9JWidphaSzJG3dn66bmdl4ZtZVkDQC3Ar8\nEDgK2Bu4kPTGcGabpscBewAfAx4EXgecC7wWOKanXpuZWVu14Q6cDGwNHB0RzwDLJW0LjEo6PyLW\njtPuf0XEk5Xb/yTpeWCppN0i4ie9dd3MzMbTybDMYuCmHOwNVwGzgEPHa9QU7A3fy7936biHZmbW\ntU7CfV9gRbUgIh4D1uVl3TgY2AD8uMt2ZmbWhU7CfQRY06J8dV7WEUk7AR8GroiIJzptZ2Zm3ZuU\nqZCSXgosA34BnDoZ92lmtiXr5IDqamC7FuUjeVlbkgRcAewPHBIRT49Tb7Rycywixjrom5nZFkPS\nfGB+R3Ujom5ldwCrIuKEStluwKPAkRFxQ037TwBLgEUR8c0WywMgItRmHQHt+7nlUNtt1dEavD0r\nvD37y9uzf9pvy7rs7GRY5kbgMElzKmXHkw6o3tG2a9IZwCnA21sFu5mZDUYn4f5Z4JfAtZIWSvoz\n4Gzgour0SEkrJV1SuX0C6QSmK4CfSjqo8vOb/X0YZmZWVTvmHhFrJC0EPgVcTxpnvwgYbao6g03f\nLBaRPludmH9eXCVwEin0zcxsAGrH3AfeAY+5d8ljmv3l7dlf3p79M/gxdzMzm2Yc7mZmBXK4m5kV\nyOFuZlYgh7uZWYEc7mZmBXK4m5kVyOFuZlYgh7uZWYEc7mZmBXK4m5kVyOFuZlYgh7uZWYEc7mZm\nBXK4m5kVyOFuZlYgh7uZWYEc7mZmBaoNd0nzJC2X9KykVZLOkdRJu+0kfUHSU5LWSPqSpB36020z\nM2un7RdkSxoBbgV+CBwF7A1cSHpTOLNm3cty/T8lfSHiecB1wJt667KZmdVpG+7AycDWwNER8Qyw\nXNK2wKik8yNibatGkg4GFgFvioi7ctkq4F5JCyNief8egpmZNasbXlkM3JSDveEqYBZwaE27xxvB\nDhAR9wEP52VmZjZAdeG+L7CiWhARjwHr8rLx7NfcLnsgL5umxobdgcKMDbsDBRkbdgcKMzbsDvSs\nLtxHgDUtylfnZd22W1PTboobG3YHCjM27A4UZGzYHSjM2LA70LNhTIWMIdynmdkWpe6A6mpguxbl\nI3nZeJ4CXtFNO0k1oa/2iyfNOcPuQAfbqqO19L6KvvD27J/hb0vw9uynXrZl3Z77CmD/pjvbDZhN\n6zH1artWY+vjjcWbmVkf1e253wicLmlOZcbM8aQDqnfUtDtT0iER8Q0ASQcCc/OyF0XEVHiLNjMr\niiLG3+uXtD1wP+kkpvOAvUgnMV0cEWdV6q0ExiJiSaXs68BrgNPYeBLT4xHRbgqlmZn1QdthmYhY\nAywEZgDXA2cDF+XfVTNarOt40t7954HLgfuAt/XeZTMzq9N2z306kvRHwCnAfyWdbPUo8I/AxyPi\nZx203wc4gfTp5OlK+YmkN6o5EbFuAF23JpIeAa6OiNOH3RebGiSNAme1WHRrRLy5h/W+BPgQ8JWI\n+JeJrmcqqRtzn1YkXQi8jxTCFwK/AA4gXUZhLnB0B6vZh/TP83ng6Zq6NliBp87a5p4GDmtR1out\nSa/7hwCH+1Qi6UjgVOBdEXFZZdGdkj5HutZNV6vsV9/MrK9+HRHf6tfKJL2serOD+rMi4rl+3f+g\nlHQ991OB7zQFOwARsSEibpL0LUlfaF4u6TJJ35V0KPAPufhhSRskPdRU/dWSbpH0jKQHJG12HEHS\neyQ9KOn5/Pv9TctHJf2npP8i6Z58OeXvSvr9CT/6Psvb5D5JiyR9Pz/eOyXNy8v3zNvnLa3aVW43\nHusbJH1b0rq8nj0l7SzpHyStlfSvkua37orOlPR4rvelfPG6xsLZkj4laUXejg/l29sMatt0ayps\nS0mPSLpgvG0paYakn0pqPp6GpDFJ1w5k4/SZpAWS7pX0XH6cn5b0G5Xl8/O2fnNjewGfIn3KB/hC\nXr5B0u6V5+YESVdIWg18Na/rnZLukvSk0qXNb5P0O039afvcD1IR4a40XnYw8PWaqpcAxzQ92XOA\nPwIuBb5Lmt0D6eDvQWx+EPj/kS5d/IfAg8DfSdq1sr53A3+b67wVuBq4UNJfNa1nNulA8//J9/9L\n4FpJszp4yJMhgN2B84FzgT8BXkm6cFwnbatmA58jDZX9SV7vl0iXhR4jbeOfAtc0PX7l+gtIl47+\nAHAE6Xmsrnsm6RLUh+ffC0jbfaqYCtsyaLMtI2I9cBnwzuqdSXo18EbS62PKyG9GMxs/uewAUgb8\nnDQEezbp+Nk1LVZxKfDPwJGkbbAgl59Let0fBDxeqf9x0tDPMcDf5LI9gS8Cx5K27U9IIwVzK+16\nee57ExHT/gfYCdgAvLum3rbAM8CJlbJ3Ac8DI/n2W/O6dm9qe2Iur7bdAfgV8Of59lbAKuDSpraf\nJl1X56X59mhe1/xKndflsjcPe3vm/lyWH9telbL/lvu4D+kfewPwlhbt7qvcbjzWN1bK/kcu+3Cl\nbP9cdnil7BHgCWB2pewEYD2w3zj9ngkcktf1qmFvx+m0LUnfv9D8f/kR0pvFVsPejk3boPlnIfB3\nwL+RJ4rk+sfm5Qfl2/Pz7Qub1jsnl7+zqbzx3Px9Tb+2yv97DwBndvrcD3JbFbHnXtH24FtE/IL0\nLn5ipfhE4KsR0e5yClU3V9b3FGkvobHn/ipgZzbfa1xGemP57UrZCxExVrn9QGUdU8XDEfHjyu2J\n9vGFiLizcruxzttalO3S1PaW2HR20nWkPfoDGwWS3iHpn/NH7BeAxn3t02U/B2nKb8uIWAn8E/n1\nIUmkPfkvRsSGLvs5SE+T+lz9uRd4A2m2SzUHrgV+TXrDr7qhy/vcrL6k/SV9RdLj+T5eIF0t9zVN\nVfv13HellAOqT5KGNXbvoO6lwJikPUnz83+f7q4x33y1yxeAxgGZnfPv/2iq07hd/ZrBTb7oJCJe\nSK8lqgd3hq3VY4Xu+9j8pS6N9by4/nEef5DePKnUWyfpGfK2VjrmcTnwGeCvSdc12gX4ygT6OUhT\nfltmlwKfkXQKaWhid9LMsank1xHx3eZCSTvR9NqLiPWSnmTT1x7N9TqwSf18TOdm4Gek432PkjLo\nEjZ/Tvv13HeliHCPiF9J+gZpzLXVHNhq3TslPQicxMZhlJvbtelCYx79K5vKd8y/n+rT/UyWdjMH\nns+/X9pUPkL/pi+KjdsuFUizSR+hG9v6WOCeiHhPpc5UPAt6OmxLSJ9s/xY4jjQOfU9E/Fuf+jBo\nP2PzxzgDeDmbv/a63a7N9Q8mfWJfGBE/qtzf9i3aDmXmXUnDMv8bOFDSO5sXSNpK0uGVos+TPnq+\nA7ii6WNcL++q/04anzyuqfw40kfJH0xgncPU7gXwc9JY4otH/fPB6d/rcx8WVQ+Akw4YBvDtfPtl\nbHzOGv57n/vQD9NhWxJpit+VwHvy8s1ml01h9wJvk1TNtaNJO7F3tW7yom5f942D1S/+70n6PWCP\nFnWHcq5GEXvuABHxj5IuAi6VdAhpSuMzpCtRnkw6OaExm+Zy4GOkN7fmf97GXsrJkq4C1kVEu1B+\n8V05IjYonUG3NH8UvJX0dYQnA2dERHMITXXj7nHkx/pV4FRJj5LevD5IuqhcP/dUngNukHQBabjl\nAuDaiGhcXfQW4NOS/ifwLeAtbJz5MJVMh23ZcCnpf3Yd6SDldPFR0gyY6yR9ljSmfR7w9Yi4t13D\nPJT1MHC8pPtJn6bancx0Nylf/m/enq8izc5ZxebP2VD23IsJd4CIOE3SN0l7HV8mvbs+TAr6j1fq\n/Yeke4EN+SBSdR2PSjoNeC/wF6TpTa9uLG51t03tL1E6KeJ9+ecnwAci4hNNbab6mZfj9bFa9h7S\ntLzPkD72fox04OqALtfTrg9Xkl5El5KGEL5KmiHSsJT0/LyPtNd1M2kWyN0drH+yTJdtmSpGfEfp\nC+1vj4jmMf5hG/e1ExH3S1pMmqr496S5618G/rLFOlo5mZQTt5CGyOaOVz8ifi7p2Fz/OuBHwJ8D\nf9VUv5fnrCfFXVumE5JeTgrdUyJiOn3sNJuQvFd6dUQ0B12rugeQhhAXRsTtA++cDURRe+518jjm\nAcD7Se/qVw63R2aTppPT6ncgDWOeC/zAwT69lXRAtRMHkj6u/y7pZIXna+qblaKTj+hHkc4R2JFN\nzwWxaWiLHJYxMyvdlrbnbma2RXC4m5kVyOFuZlYgh7uZWYEc7mZmBXK4m5kV6P8DvRNHIsPPPGkA\nAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10adb7610>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def solve_heat_equation_2D(u,D,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,D,dx)\n",
" return u\n",
"\n",
"def solve_heat_equation_2D_vectorized(u,D,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,D,dx)\n",
" return u\n",
"\n",
"@jit\n",
"def solve_heat_equation_2D_numba(u,D,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,D,dx)\n",
" return u\n",
"\n",
"def solve_heat_equation_2D_cython(u,D,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,D,dx)\n",
" return u\n",
"\n",
"def solve_heat_equation_2D_fortran(u,D,dx,dt,T):\n",
" N = u.shape[0]\n",
" n_steps = int(T/dt)\n",
" for j in range(n_steps):\n",
" centered_diff_fort.centered_diff_fort(u,D,dx,N)\n",
" u = u + dt*D\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",
"D = np.zeros_like(u)\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,D,dx,dt,T=T)\n",
"#times['Python'] = y.best\n",
"\n",
"print 'numpy vectorized'\n",
"y = %timeit -o solve_heat_equation_2D_vectorized(u,D,dx,dt,T=T)\n",
"times['numpy'] = y.best\n",
"\n",
"print 'numba'\n",
"y = %timeit -o solve_heat_equation_2D_numba(u,D,dx,dt,T=T)\n",
"times['numba'] = y.best\n",
"\n",
"print 'Cython'\n",
"y = %timeit -o solve_heat_equation_2D_cython(u,D,dx,dt,T=T)\n",
"times['Cython'] = y.best\n",
"\n",
"print 'Fortran'\n",
"uF = np.zeros((N,N),order='F')\n",
"DF = 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(uF,DF,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());"
]
}
],
"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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment