taldcroft/perf_python.ipynb

## perf_python.ipynb
{
 "metadata": {
  "name": "perf_python"
 },
 "nbformat": 3,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "source": [
      "Wicked fast computation with Python",
      "=====================================",
      "",
      "- Python is a dynamic language with no built-in static typing.  ",
      "- Pure Python is generally slow for numeric computation.",
      "",
      "In this notebook I go through a series of techniques that are available",
      "for making numerical code fast within Python.  This is _heavily_ based on",
      "the examples found in [Performance Python](http://www.scipy.org/PerformancePython)",
      "and the [Speeding up Python][1] blog post by Travis Oliphant.",
      "",
      "[1]: http://technicaldiscovery.blogspot.com/2011/06/speeding-up-python-numpy-cython-and.html",
      "",
      "Before you make it faster",
      "--------------------------",
      "",
      "**Profiling!**",
      "",
      "Sometimes you'll be surprised where your code is spending time.",
      "",
      "Two places to start:",
      "",
      "- Standard library [profile](http://docs.python.org/2/library/profile.html) module: by function",
      "- [Line profiler](http://packages.python.org/line_profiler/): line-by-line",
      "",
      "Options not covered",
      "--------------------",
      "",
      "- Parallelization (e.g. [multiprocessing](http://docs.python.org/2/library/multiprocessing.html),",
      "  [IPython parallel processing](http://ipython.org/ipython-doc/dev/parallel/), MPI, etc etc.",
      "- [PyPy](pypy.org): alternate implementation of Python that uses ",
      "  [JIT](http://en.wikipedia.org/wiki/Just-in-time_compilation) compilation to achieve",
      "  high performance (geometric mean speedup of about 6 over CPython).",
      "  They are obsessed with [speed](http://speed.pypy.org/).",
      "",
      "The problem",
      "--------------",
      "",
      "Solve the 2D Laplace equation $\\nabla^2 u = 0$ with a specified boundary condition on $u$ using an iterative finite difference scheme (four point averaging, Gauss-Seidel or Gauss-Jordan). ",
      ""
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Pure Python solution"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def py_update(u, dx2, dy2):",
      "    \"\"\"Take a time step using straight Python loops and update `u` in place.\"\"\"",
      "    nx, ny = u.shape        ",
      "    dnr_inv = 0.5 / (dx2 + dy2)",
      "",
      "    for i in range(1, nx - 1):",
      "        for j in range(1, ny - 1):",
      "            u[i, j] = ((u[i - 1, j] + u[i + 1, j]) * dy2 +",
      "                       (u[i, j - 1] + u[i, j + 1]) * dx2) * dnr_inv"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Simple function for testing speed"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def calc(N=200, Niter=50, func=py_update, dx=0.1, dy=0.1):",
      "    dx2 = dx * dx",
      "    dy2 = dy * dy",
      "    u = np.zeros([N, N])",
      "    u[0] = 1",
      "    u[N - 1] = -1",
      "    for i in range(Niter):",
      "        func(u, dx2, dy2)",
      "    return u"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u = calc(func=py_update, N=20)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.imshow(u)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<matplotlib.image.AxesImage at 0x1045d26d0>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD5CAYAAAADZljUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvU2sdldZ//+91r7PeZ4Ho0ETWiutLfkBP4wBQtI4q2kC\nRRiIMFGIUULBGBJMDAxwohYHWgeO0IFRTOqkwQlCgjbEQQGZdAIJCSSQaGvTlP5ENPFPn/Ny73X9\nB+vtWtdea+21932fc5+nPdfJPuttv+/9Wd9rrb32vomZGdd2bdf2qjBz6B24tmu7tsuza+Cv7dpe\nRXYN/LVd26vIroG/tmt7Fdk18Nd2ba8iuwb+2q7tVWSbi1oxEV3Uqq/t2q5txmpP21cD/9RTT+H3\nf//3MY4jPvaxj+HTn/70ZJ4/Liz3NICH1270itvTeOUeG3B9fHeKfaZRtsqlH8cRn/jEJ/DUU0/h\nO9/5Dp588kl897vfXbl713Zt13ZZtgr4Z555Bm984xvxwAMP4OjoCB/84AfxxS9+cd/7dm3Xdm17\ntlXAv/DCC7jvvvti+t5778ULL7zQtewDazZ4h9gDh96BC7YHDr0DF2wPHHoHLsFWteF7O+Sefk1a\n/QNHBg8cDXgDAwQGmEFhAoMYIu4mhH4HGdbiMkQhXerD4HJ2yOx5yUDO83875u9d71W0Nx96B5St\n6RZuLfMWPU9hZirlt9KlODXiPmQiN4FEXOSTW5AJePZ8xLPnNm3n5W31GFcB//rXvx7PP/98TD//\n/PO49957J/M9/MBrpwtbC2ILYxmGLYy1PszTZF0lAAZgkQCX8VIamK8UAtCU8iKEXMlXq5nLa+XP\nlbXsTq0sarb2WU5ruVpZKb/IKvXHi2Etj+B86lJcpJkI1hCYDKwxsGRgDYl4ygMZPIDcO/nqd35Y\nOQMrgX/wwQfx/e9/H88++yx+7ud+Dp///Ofx5JNPTubju39ikmfsCDOOGKzFYEcM4+hCa30cGCxj\nGAlk2QHdM3H/FJ5YxBAiTYV8eUwd8VrerhXAKw32YD3QLwF8DuweAQaA4Mg2wyWT6ZusAcbBgT2a\nAeMwuNCYGIcZwIOBNcP0YPcN/GazwV/+5V/iV37lVzCOIz760Y/iF37hF6YzFoCncQszbjFsR2zG\nrZhGbEZgs2VsRsZmtKARwAgHdCku0wH6OQ/ApxkebK/oEvJY1oC/BHurAtjFU5gr6yk/tM1BvVSx\n13jVTbhDKGAmSmpei/codowPfqrFfZoHYLshbIfBTxtshw3Ohw0wDODNBjxsMA4DMCxDePVz+Pe+\n971473vf257p7tdM87Zb0PYcw3aLzfYcR1vjp3McbdlPFkdbwGzhYA6hjJfKJOBB+Y3KI4CtA5kJ\nUfEl6K6tlOJAX7jUE+iF/iKaBoeyXcHWed1q3RFK6CO3lE8gkR+AlmEtz8DRNohwKORtADsQzjeE\n843B+WaD880GtDkCNkfAZgPeHMFujkCbDbC5JOB7jO8qKPz2HOb8DOb8HJutwdE54ficcHzObtpa\nnwbMORzQpSmUDUjQ19z9AD35EB5mbrQGOE9jQVyGpbyl6r+vfoND2M7t6Up6jYL3xkugh8lI+Gtu\neaksQH3kw9Lky+wGODsyODsyGI4GmM0GODoCjo7BR8ewR0cYN8egoyPw5qhwJut2ocAXXfrzM5jz\nAcPZgOGMcHQOHJ8BN84sbpxb3DgzuHFGuHlOMOfswC5NAfRzuBO6RQ65rgCERaAZsLbg8XPuLACY\nVAg6D4V4T4hCet/Ngatg+wRcphcpt5p0HiC88ALgZKbpCLOGW6cHOKD1tJnmjUfA5pgwHBmY4wE4\nPgIfHYGPj2GPb2A8OoY5vgEcHbtpgV068DjbwJwZDKcGmzPg6BQ4PmPcOB1x82zEzTODW6cGN88A\ncwYgTOcq3PhQtpFkm16ngXhVI/DWdZCwFV5/aNOHOOoTZvJr4dI2/5IK4CqBvwRone6FvRf4Fuh6\nMiS88gC+ERWAB7/YFtdhiG8AHMNBfaziKm88Jgw3COZ4AB0P4BsbB/uNY4zHN3B+4wbo+Cboxg3g\n6AaW2MUCf9e0DU9nA+iEMJwCm1PG5sTi+HTEjdMtbp5ucevE4NYp4TWngDlFAl5CHk7i5ErBwR1c\nd0KCPViAWMJuktJbSqAv7Pxf7PLv0ubvUf2rZD0AL4kvddmXTAZJyWUoYTdayYfGZJCgvoEEt4yL\n9HgMmJsGdMOAbwywNzawN4+wvXGM4cYNDDdvwty4CbpxCzi+QsCXHsvR6QBzApjbwOaEcXQy4vhk\nxI2Tc9w8GXDrZMBrTgg/cQIMJwBO/STbPhL4uGI4tz5cAA163KkEs5XQe9AtuYk5h7801SoENNKl\nOFR8nxXBoW0NzDpey6u55S3gSzqh1T02xQXoAfKg9sZ4l16CvcEU9pB3BAd0x7S9AdAtAm4Y2Jsb\n2FsbjDePcH7zGOc3PfA3bwE3XwO+cRNL7PLb8CcG5jZjuM3Y3B5xdHuL49tb3Li9wc3bA26dGPzE\nbcJPHBOGYwAncCfrFFNl1xba8iUJABwNgwAeHnoSE/v8GaVvPf0D2hVAT7u+p6wVv0q2Bu6esjlX\nPaRrzmBR2SEAV2kJ+wT4Qm97FgYFvynCSnx7k4BbBL5lMN4aMN7aYHvrCGe3jrG5dQPm1k3QzdeA\nbr0GuHlresIbdrHA3yytfgB4AMGAMPjJwJCbBkMwRBgIGFRP52SgQ7BwtUKv/UbER6QOPq/8PAJ2\nROy0i9Mo4qIstvlRV3XZ2Qe0gYeKAwXAuV7WG78K1qvYkzi1l5lTd5kX2uSz0Aslj9Og0mKagL0p\npEP8GAlsPXnQ2af5FjDcJJhbhOEWOff+loG5NYBuGdDNAXRrcHwVGavbxQJfMtFIYn8y+MhBGJ+P\nMxDGEIerQcTuubmYSF65AHWYRkwf5Yk8sm57QepJTUaUSeAJeTrEDfI0ZBoiLcpKIcQ8tfLJ/Dp+\nVainPrhlugR9qZzEPFX4Negl8FWehJzUNOmFL4GtpwrwXFJ4nx8m3AD42E04AngDcHA9al7ujF06\n8AFa9ieM/WMJtuRhDy8MQAGeww9fHp+HauDn4FegQ4Luld4UgA99AKxDOY+Oy+UKeVB5Imjm96Tr\nmRdghRuwBnaWpmkZicgkvwD5JI/y8nj7UDkvzGsMQIMDPELvAZ9APwe5nDzwrNSdFfAhj29OYZ90\nBN4JwAOIJywqvCUHFoc3glwb2hABxCABv4uTV3en+jQHvJpo69Q9qDx5yDX4LB7pZYBL+FW+z84B\nVxVBnEeVyXK5niyNQrpSGUzsosCfufFIRZqVgAZcq7x29SXYlM83UXoS8+gKIIQecAl7VPgB+bP3\nGvClSbv0N3L4naJTUvljgI8JOHYesISeDdILXgvtMC49UYJ+E2B3xeyBZkPOC4hXIoEvw6xB1gl8\npvBigE4AP8A+Ad4TLCGP6QrwqOQVgS9UBs20MD1P0Q4BfAXyLE/P00qXoFZwo1A2gV2Wy3kU5GbI\n8yAUvwv0I7iK4Uau5pwpPKk8cs3cAH2AfeNF8k5y6YM7z+GESYAIEXTeOJVPFz1X+6T0zr2fKHxr\nWO42wS1H5gXQrVd7a4FB9CtIuNHKg1BqDbkuL5SVylVxO7+a2chfa7WbruXe03y+BDcGhfIa4JDl\n0nWnPMyg18BXQrSArw2dFQofwb5JmdKz6MXnIw+9hn1AeoRwRwAPJHVnErAnVecB4K3stGMx7InL\nCi+BP0fqqa/BL4EX7rwE3/jyCHN4004Dr/OBKdiMYiWA0ryqLGO05uI3ZX8mbxdrSnc5T7vwWVSB\nnrn4BbgBcR+EeJhHQ98TasA15Hr0XA1ynVfqtFOwx3hw4wPoG4rQ79J+Bw7m0iOdPAAg30nn2/W8\nIddrH10XApnQS5+gh5zk2PoNctA3mEIvQVev2VJ4YuDzmJGBX4tzGLOvwZdAVyqAaIWKQs/TdN+r\nDf3GMmutB/bGPEXwC+DqmSXgEuxwL5CIA4hvtcl7RsMOr+4QoLfiWfv9qCOuO+1uUBl2DzwiC8j6\nCwIjd47Ce9gzNz7k2QS7tQRjHNxuGeG+k8jXwEvYdVxC72FnBX2APYQS+AzwoOry1Vuh8hmoGm5d\nCWA6v85bBLnOuyjoW3CvUf4S5K2KQCu7nr8E9Fw6tNEV6NCwl95+O6qEJZc+PHsvqDwfe7CNCuVQ\nXnmsC+wwbfgAKCG255kpPpqz4RGdcOPrYyK94g/IwZaADyrc5mDrD2xMyoLac0pH+FW6CLSeUJkP\nlXkgQlTSIq9UOVCt0vD1Zq0OmJQpQLmQn7nhuqyUlnkFaMO2srxSWgIuy8S9k0FuCnkSbh03Kt2C\nXIfSpffqzqLdHmG/SQ54wUYa54s4sVlXbx+uDc/kHy1QVEAGsh5v6dLDMKrufLhQQdE19CGU7Xyp\n7BJwAX3xNVsFPpfKamDPwa/n12XA9Co32viL4z1WU/GOOOmy0rwV4FuQFyc5n1DF4kcrZKjHwhcq\ngKLC69dddV4GfB7nCLyvCI5INF8RveDp2BQstgMMvHF7yQBA7rm7DYNtQLAwsGx8yH5gDOcj38Jj\nNDnqJVw06fYMlbww5FZ22rU+paWhViP0Jsq/ZkIlXlL4Xnd9Lr5GInYAvmu+Gtg6vXSSrrCp5Ml7\nRA906QG+BL6P87GH+ph86OO+g84eU+yRtxtyH6kkgkX+9drAzBrYgQMpfILbA04O8hGDm2jAljcA\nAYYYhgAeGMYwzMCwG4bZwLW/jxjxAxfKfZtcYHlh5Xfxhko8XPwG8Lt8THMx7Luo9xVQ+K75WpDr\n9FrYLxP4EB4DOCKn3h5qHsi9rWkINnRcE2FLG4zkWLA0YPSMMIyrAPx0Ryi8s7TTDHcgY4TeROjd\n4zcLNg52DG4yg3W1YRgpxwBCY3QOeNkEqH0YU5ftCnyrIkBHfE6dX6nAL4G+2seDZbDrTrk54DvU\n3cUD7H4aCFZMbNynqa0HfusZCNA7YRRKv1LiD6TwBZWP0A9R6Qnwz+fZPRwfGBgsaGPA1nrXXjRi\na8CXboAW5Gug39Wtx4I4VPwqA79kmYty5TXUa4E3Km8OeN1p58vthtw0EOwQvjnv3PiRyAvexocm\n8mG9ODp2HENL7fLb8NEdITAHN8XAYoANtRoP2MK59AMR3OB2AowFbQBrLci/sM7h07MlyFtxDXwL\n+otU+DWw98C/pm3fY5cNfC/0cwpfA7wGfEndddnc4zhVAWQKvzER+NE46Edyrvs2wu6ZiOouVP5O\nUngAQt2Tykt3fsTgn7tTgn1w7NMA8EbADiy7IUIbXgNeg34N8L3wY0F8Tajjc2Vztm/IZXwO7lJZ\ny4vTwPfAPqfwMt37SK4AfVL38IMTfvLAb5E83ejSYzfYgQN32rFX+eSypIPcYuMOy4xu0E0Y6WQZ\ntCEYNmBwdG8mL9KUXDoNfEnlZWfd0nZ8DfRaHnDx4Ov4XNmcXTTwu7jxJfhbYM8B39tp1zH4hksK\nvzEYh/BrMg72IHYjQsedcOfJ5OysAP8gLj2Q2iCWKD6GCwcXoDd+NB0RMBr4t5cYzBzfmwew7iaQ\ncGvAd+2l74F+KfCtvFoZCuldYA92EdDvG/aSws+BTuiHfclIOzkuXig8Dyaq/EgO+i0NyaX3zVxL\npTb8HeTSB9hDDz2Th1102CWXHvFVRTswLDMMO20PY/CrKt66OXqUvUfhuZLXAn1pG35NuT7hc/El\ndpHAL4W+BXqvC3/RwG/geumHpO6uDe9+G240gwtpEAov4wH49Fx+rR3Wpce0lz536Tl+ODDAbgeG\nhfU/mQukbwlzf4ddDfhdXfpWG75Whpn4XPkc8Pt054E21Dq9Fvge8OcUvQb9voGvAl7IG+CewQ+i\nDe+n0Qxe3Qffht8I2IepS38nKTyACPuklx4Dthiw8S6945lBZnCDbmDdz+iSgz7CKy9a6+boAX6f\nLn0rDuwX9qXAl9JzdpnA1yqAXth7VH2fwLfgHwAeCBi80g+yl9634Y17HDdx6YXCszjAO6YNP1V4\ngnzWGMA3YBhijMY5/3HILVmEjgtrKX3yxx8/hX8ttdCdNC3wd1X4ud57zMTROf8rHfieyrxWwffC\n3wu7weyrsSzjQxhZl8Ksl167817dgzsfm7+RnXV2MIUHJPgGcRy96LgbYDHCwPgDD5+ytuQeYRhR\nFpXYP7oLb+OQAIFSJL/II9ru36GB74X8Mlx6YHfoLxr4FvRLgW9VApsANWUuu4uT+EINMBrCaMjD\n7Z+9U+iZN7FjLuuVRxpOmx7Jrb9swEHb8ICEXj+PD236AHiYRjIwRlQCMCAYB/cAwL9wAyYQM5gF\n6Ewu3roJxkLeRQO/K/SlOArxubJeK3mSPdC38q4C8Gtc+iF9jYY3JFx3eOAd+GPojTc57LIH3mZp\nCbluu9+hbXgJe6rN9Og7G9U9wD5ar/YmQQ8GyDDIEGjgTNkBpe41174G/b7a8C3g16p4LQ8L4mtt\nCfg94T6B1y5+Ce6lLr0q56wzjsSnqMh/sCKFAXY50MaaXMimQ2jzwTZJJNfbAV36cK/q3nrK3Jrw\nPH4Ew5CA3xgY66EnAwAgZjdBXndGuqv8KZMXu6ToAfJDAj/XPm+BfkjgddlVBH5O9btV3gEeXobh\nDWJbnSX0htwAG5oq/Eha2Q1kn1YNegCrVP7gnXbZ8/h4oKHTzjq4YWDYdWwY2Aj7aFw+uU/muO+J\ne7edmGFYbtWb/IJOC3qZdyjge133QwEfrHbfXTTsLchrwM8p/hKXXkIf4I6DatIbcTwgKnoAPEAf\n81CCnhTssh1/h7n0UAcRarXkzqdODAPvzgfFD9CzwcguBNyAHAzuVBj2Iu3PS/ohC9+G15DX2u6X\nCXwN+l7ADwU80Ib+qgDf494vcOmjWy8BV6++hvgoYScFe8hX9/60Zz614dfaAd+Hz9351Bup3XnV\nUx/gF7CPHngMbgSegXvI7XMBcP5Z61Zb/dDA90KPQt4hgQ+m70UNtc6TZZcNfK0CWKjyqb2ew24H\nk95zN7kLP3HnVRs+77jLh9PecW345JJMXXsNu2XfcUcJdmKn8kTCpScGw8K15gOlAIj9T1a5eHYz\naOhrFcAS4OcA7wG+tw0vw7m8VvyirAb7ZQG/FHbCctANYuecVbDHZ+1e4eNoUunCB/iRg586sGtu\n/R3n0of7eOrGa5eeyA3CGWFB2sVnAwKBrPE3sQOdghcRPwTI/ocu2H3IULfRW/FDAV9ru0Plo1LW\nG78oq4Fdi18W8C34F7r0GAgQ0LMRKm887IP/sIUAPQyqGUl85k2pe+jT2sfjuGBXpNNuOq4+V/k0\n8CY+n+dQGZB39xnw+k6wsHA/OMkjOdiJAMtlda+pfIjPAS5H4y2B3caT0g8+0Ib8KoBPhXgNcJ0+\npFu/UOFDBcCGYsgCdOt756NLD4pT+EBl1kFHUwZ0Zx0QvOMULrEr1WlXbscL0GEyhY9xOJfeEUQg\n/9yNyIKI3Pfuid0nsWrq3lL8uUnOVwO8pOwyD4Xy3jY8FpS34vu2OcB12UWrfAv+EF/aaWeEWx/a\n6wH6oPCyJ37BVBNGYB3swEEH3qSw3J7X7fpaL/6QTkP8RUACyE4+ikFxm27rlO2JyyACMIq2voKb\n5oC/CJe+pvD6RLbKW/GLtB7QdfqAwHMGN4kRc1DAp6GyukOuPJhmWAx6Dfpd7GC99IAEHeJgager\nYXedeQl8CTHyi60sfScHIBBCZx6NAPvPatEIgHgCuP7hiUkF0KPuu7TjoeJrQb9o4KkQ9sZbE9AH\n+Azw3FD44KIHBY/gh7hfXo+es2JQjfRQQzh9+638GC6/9/M2PLBe3YED9tLn+poPKJA1XHsaMAp4\no0SHKSivAt8Iv8L9ko373TqGi8fe/IKalyoAkrD3gl6DHpX8qwj93H1Xg7qU1wt+KC+100vpAtjc\nUvgIPdzzde+uR/U3FMPYNhdDZK3ucS+68kMF9Pr4+emAm3XQH9illwqPwgG2OvRcG558msDOHWeo\nm0OA76Fn37YPP1sVJhD8L93AdfD5L+NK0FmDzztAv6tLvw/Ql0LfC7lMl2AvlfVAH8pWAF9VdYP0\n5SQjetsl9IayOA9e4cl3zhHlg2gKKq/jVoGu4a+126XELbUr8FguufZAHXQDk7nwhFzlnQVykbt/\nHnb2lQHbsCWn7CTAT8/q/U9Ue5g16MwFd35Xde916Q8Bes1qFUAL8lrYA33I30HhW2rPcvIdcRwr\nAgF8aL8HuPUbcBH6ASVVr7vz7efud2gbPlgd9LwNYyfQB+CdS28mdV722S8DZO/HAwl0i7yvjwBj\nkJTd5Aoe41rZ59R9V6WHCGvxtRXBWivdeyWF12EP9Dqt83ogl+lFsCfQ8zRlZaOh3HUnKqs7mSL0\nSclTxzTDqGl/HXbAwX9qSt/XWt0D4gaUge4gDz30IwZ/Q1D5lBAA8cOTTHCfuyYbP4cXPo0HA1jL\nHnpk0AeFL6o7F8Il0GNyMq4u9DWoS3m1UOf1QK/TNcgLce3OcyUe5tWgWwW8/KWY0mg5p/Tpi7Pa\nlZ9O6bG0Hl1XEsX8pPTbauAfeOAB/NRP/RSGYcDR0RGeeeaZruWm7Y+2wtuo3wn0EAbYJ3WfSKTf\nqvCuOjtgDVk3Oo8YZB341gDGsuuQsXC/ZycgDq49CXUP8EvlXwT9UpVHJX7Z6g7Mg98LewjnINdl\nC4CvteV1nGU8/Mijgt7llTvpJnGE3vm6K19Sew29Vnd9S/TaauCJCE8//TR+5md+Zu0qAKHyKV56\nDp/j79JDnGPEIFcZVTwqe4Q+hAy2Boash59cfeASHnoBs4RauvIB4KUu/SvhsVwJ3lbZHOwh7I2v\nAV6quXLpWZSzIf+bhqESCD3y5D6g6sF3I+YcqNNXXAfVQ99qy5OCv9RLn9hYo+zBdnLpmXe5a9KT\ncP0MXrrz0/Y8i9whP/wAdtXF9ENu2H0ckxkwnHx2w+R+5MIQ4D+P1YI8c+17XPraPCWwrzL0Ndjn\n4j15Fwz8XCWQoIeH3g959ZWAjelcxWV7PcujsrrzJJ7GntSgDydhF+p2Uvh3vetdGIYBv/u7v4vf\n+Z3fWbS8VPWausfHbYWpSgSprUjQw3xswB5wsPHyTWJyH89gJvfLtULhY7wEew/cr6THci2Qe+K1\nvFpc5u0IvHTdp6AHVYd3333aq3vooBuDugtg5UAbDfmcay9hn7rziZHpyeu31cB/4xvfwD333IP/\n/M//xCOPPIK3vOUteOihh7J5vvbY12P8/od/Hvc/fH+2s1Poa269OwH5n1Z7uRbEtegtpM2HO8Kq\nI/P5RI5wDykZXxJW4esS8qtwPfycAxo8DoM01N+Eck6b3gF40nlr4mtsJeQs83cFPoyx8OlMqRXw\nHENS6TzMoU497zEMcFMOcQnmVlk+mo6KeWV33p0Eefmee/o5PPf0f5Svk7LVwN9zzz0AgNe97nX4\nwAc+gGeeeWYC/C8/9lBp0YifjNcfy01d+gS6+5Pt/OAVAOlR3TQEAJvzFH+G10YeQ158tG/8Y7wA\nMDz48fEeCSV3NIaPaTIDEJ2AxCTgD5OqMJCWr6kzQ+U3AJ+tHHpsDmIZX6H28vcFUkhF4FnE4UGe\nlHl1zuD283OMh/kpU/A0Ui53t2O7XQHdA3hd4eug54DrTm/g/ofvF2IKfO0z/4qarQL+5ZdfxjiO\n+Mmf/En8+Mc/xle+8hX88R//8aJ1JPS0g15T9qDoZYWfBzy/z8IScYgde7kVv7/N7H7dxt0vfl0G\n7kOZNsXBFN7XAZgdpEQpLibykHOhzN2BIe7KYzdJC+hOpZeVQ3H+HpuDvBLnTuBBYd0FyEM5ifIM\nfJ1HAmaknyczIu7zw3wadtYKLyqAOTe9H+zWJKGfwr7UVgH/0ksv4QMf+AAAYLvd4jd/8zfx7ne/\ne/F6SpCX1T4ouAbdpa1fUoaAvpf0qfKEBlGlBLkhuR/slJ3hnt2HpUMehxDepZegk68Q4PPItf19\nJRG/mV9w38l3IE5AL8FaUnokL2Ry+H5+Kp2WOSOxSAXiXriLlQa5yETpdUVQBF/Fw/4amkJPSeHT\nR1Iogs0VwGWPeumx25zKt57D1z9cufwy1WwV8G94wxvwrW99a0+7QNk0VXf5O/CmALpUd5PBXQI9\n3V/BpSf/WC79fFVQdrdWm/YsqDvnITy4FAFPoFME2EPmf+o6lLu43085r64Awm4HuDvc9eAdTMDm\nFEjeeiw7mwsgz+BulUuPIAOesnlYeAASWF3OSGUa7kkaKMMu2vP6mw0l2Hs76nKwyyPratNaO/BI\nu7o7H1x3CX4YcJMUvqTuAFTcbTHFQ8yIJR3gVuyb7CUQe0wKesgKwMPu4Y9KD0QvICo8XDliXpAw\nqfhpHRMFL1UColxCzoX5pQXwZQWg84qVQ0OpW15ADf6o2GEdXu1THAn+zBMQIWQ6AJ/gjt6cAj3k\nMUqwB9C9az8DdZ+6i3VOwE/eblYLZtesnD9nB3w9Vqbz9koJfKn2ugJIoVxjbjodYSdyz+KBqO6G\nyD2SI9cFGPcswE4V4MEefL+MB5nAEXTyNLCIlzrrJmrv83aBvOTGcyWcKwsb1GAvgX+Sr8DmFvwa\nbl9pTyoCSJhDXpp3Gs974SPsyrVn0CzsebwMe/Jcay/N5GIobY2bf5Bv2sm4aBXHA0zqHtz55NYn\n2OF1eLrW5NqXtuq2BbJg9tWFBzzG/VZiGx4BcgF7SAfIQVkaSKpOyF16Ny9Sx51Q+6Dq8VBEJVBS\n+qWQs1/PzkbqhtsB/hzyWpzSvBnsEuRSngI7C/N8kHgNm0xSe6RKoPQZtqUqr7/ipD9FXXbhg4+6\n3p0HDvxNu2Ba3aEOtqzoaWBO+hJ9cPn10/XpaWIgjbbzlDCc2rt2fXrQN3HnSYLuL1X4iasIJ4l5\nEviRSg16AXoCUnsfZaUvtdM1+BpyZkzOx1LTMMs2fA/ous0fAI1lUa3zuGybp0ogAKsBFh6jfPpC\n4n5T+RHAFuhUV/h1kCe3Xv+sVOIh3cW7QH8FvmmXDqLsxtegD0sZsSYDPYwmmdZ459KT6LALCp/3\nFlB2CYhXQc8eAAAgAElEQVQE4BR8gLwCyEBHrvJZnmjLOwgl9OJESZX3eaGYFOR6nklcnvxdTSq5\nT9ficb4i6L5chsKNTx14lDraEOrGpM4a4gz4CDRN8qS6Rnc7e/5ehr00zVUA+mWZ/FHd3PP4O1rh\np7DLeLoAQcU19KFmDEu1YM+3Kf0EDXnYlgt9lePBNmwT4JTvaRoQOYXbSMDk+wf+0VsGutvJOFFY\nJMRDOUSZNwl/Ne4jWQWw0IpQQ6l8Jd6qCCYqnlUA0mUX9wwJKCi/d0JZ7YMSky/LkBpgM4mX3fme\nSbfh8/UV9kUxIUMdX2IH+334WrwMvnw0p0GX6WUW4Jbw1/II7tv2tdslVBP1suAhOBfdwLquhPD4\nzjfQneeAGAfSPG6fkVFGlfikzC9HlfIlVnLfq+DHjU/LNPzSRY9lAnA9z1Spxb2jy4WAlOGfg70c\n753ysfLtX5RpQa/jS+2gvw8v4zXYy+15iRWDGifLRpd/uv0AJZDDr/Nq9S+jDLguyzr4fNrCRLcf\nviMQ5NMBenEGpIUKQcddBaLnC8eCCB7J+buvWDhvIl4Em6rzTAAXcTevhNnni572tJ02GD1TDfwe\n0JcCr9cxhb32ZRudfkW49DTJK09a5bXaly1UCRr8HtBzla6re6l8Mmno/YSg8iTUXd6+siKI+55b\ngp4TxLJCCPNp0Ne25ecUfTIPTZYpr4OU+y6UPVuHVvdC3g7g1wCfA753G/mXaVu98+JkKfDX2oG/\naReaqqXaeqru6X14V1ZCrKboBjYrk3DXQNfQLwK8B/rgxhfCHPTppc5UHrKBX1HzsAwARM9i2bWS\nFpQ2xkN+BX6XppjWdU0GfAyl6qtyAXu4X9YAX4NzKfS1SqAMe+17dXXw99XPenDgc6uBrqF3J6Ck\n7EHRS/DndWYZeBeneGkC2MtaXQluud4p9PDufIfSF85U0juRqdvposdOzhfjS++kEqxiJ+SeTtrv\naadmgE+KXoQ/QBDb7ersU+mq1qruNvQ9wJdg1+vTbfnahNlpNzvwc/ipbqYO6vU1dDCt6jIdYC6p\necC81Sm3M/Twag4FuYR9An5+9uRInAx6IG43zSvihXb9EotLa7dcxokKy6SNTpaT+QUXPuYBeS99\n9eyjUVYGvAf6HuBL6XUTClOUilV2wLH0Dm4SefkFKrXX50Gvue46XbotJPSldP9lQjMd82SFIN15\nALm6KyWX5QVVB7jcdo/bledjmZUBl253tqV8PpqWxXQJ+FiWx/MbX51Zqp3xZfdSC3KZN9fJV0r3\nw15X9rXQH2ws/RT0UAHIi9iGXqt3yJPQz7Xd5yBnrGiry20UVD6LB3UP8VBecOvD2QpnsaXq2byl\nTr+dvkeIiVteBNznsbo3a5WB3MusVz6rUHIYctArZTtAPxfvcfP1/HZWRvQxlDu417hnBwDeS5KP\nS9Al5KQONr8YBCD/yKUGXLvvEFsMYdpCarOXoF8FuoI+xkvwl8plHoCKLuaXXINdrACmKr/GpmdU\npCtKr5fL0hp4BblcVwYE5Wd/DexhexroFuwt4Gtl+dS+w/Rd1D6X/XaAl2fSLqd4wCwcKAT0UCfZ\n/So8K9hLCt+CX55eNxpf502h18vNQS/nj3FS5ayWoekyMgTalz9Lk0zXmzdLLa9WCh11k70sp7Ot\nC4+hdpNLd7+kgD2wA2iW94BfAnxJBSCh75OM8vlYYwfutHPIy8Gs7PNrFwgzsNdAkWEJ5lreEsDn\nL9mMus9UFvLMzWhsY95Sut/aVUaeblUMk3mpXHWVjjKHW98fu4E+B34pv7dXv7W+0v6X7pzpNVhm\nB315JiGfQy+VPeSULlxJ2QEU1V2eOpk3B/3S2ydsI4Q9FVBXHpUBXwJ/T3rO2hDPqX/fvLUj5cl8\nLdBRyFteCfQMqJnr3Oud5P6k4yxDv9auwI9J5u35ALsEHZieiJobL4GeyytNNejDHrYgnwN8Lfwt\nbPYBuLwKS8qme7O0AghbqJe14C/DXm/HA3XAdVkvsHMqvqYCkOdieh52s4M9lnNuu1R3l5+c+qnu\nT0/GtM0OYDXsEvq5ecLe1vLnwr558zoeaANchnk/7faazQPutloqb/VCtN36GtiAviphm0vgXzr1\n9Ogvg7yk8JjE19gVGWlHCL8NoxUeqJ2IFLYAL5XvMoW9beWFdH9YvsQ6XNYqrsO9L+j13tTy10Pe\n245vXYV50Et5u0C/FPZ2JRDOQ37V1l7BAz6HlyoPpIvLCArvcqfgt1Rdl8sKIMTzfdm9Y64G+jr4\n+4HXt4DG7+oA77ZaKl8C+fRM9F2dMP8+wNfz9IBeg1+eCx3KY8zPR34ul9qVaMOzSMmLKmEHyi3Z\nKegMNGCvddatBz3sTRvw5VBPw3SO9NlY5rpfNPClsjbkeq/yeGmZOvAyXurcW5buLVvSm99boeTn\nKxz/blfv4L30U/cuXVp50Bp+DXoCGqBqWdmlD+tfkqfLS2W1sA55G3x9tnrhL9mu0LdgL82zVMnn\nl5m7KjK+H/BLeUvgXwI6Z+ma0i+3A78P7yy/daYHpmFPy+d5tUqglGbkqq7TYbtXH/a6NtZsP7dO\nX5VRq5Z2V/n6VUnlbZDXlun0Lu312jY1FYzdrxhwcJfeWa2lGS5jTU2kYs+BHvJkWkLutrgO7qXQ\nl8JWK66uddOz04fh8nml9ah7ad59uPI14Nt5ffG1y+xzyte9P9CDXQHgazdPWdVrNge6zlsLb2+F\nsCRcBnm/a38VbKkrvwT4EJYgT2Xr4z3z9kLc8gL0+pJRZVpvV+C79LLlrsuWWw/ocgoq77bahnaN\nqq+DPZ2Rfhx2b5dflOVVGKD3tO6jEGrVWyvctQJo5dWW3a+yTxV+X3xciTY84Bz4acn8LRyA1ZAD\nU9BlmQZdh6XbJ5Xt1m6fQh+2lR9/DZM7Rd2DtauxNvByGT1/L5ytedauS87bB28f9Lnlqp6Dv+5q\nH/QDGCTSudIn9Z/G2pYgD+/HlcoYAfkQSrVPlUG5UuitDHrC6dFNb++U7m/LXxWbd+N1ejqvXo9L\nLwN73Txy32p3wDzkedmSF6/zc5Sf03V2Rdrw8lm8hF9WCdM2vVTpqcm35csTx/Uub3fvCnmKl44q\nxeVZutOAn6+6gpVv7vbyy9z8EC6Zd35dYd/3p/L5VdTg7w79FfjiTUgnLQ8HHQ5VpvdxW6ftcdzu\nLsCvCXW8BH6aL6Xn3fnLasfPX4fyEdX2dE7R2kd+eaHc15IiLwO/dB7L4E/Pyxo7+Eg7iXA4YWtB\nDyD31KNS3bV7H7a5S9hbFva8dNndfKV0G4GSUbVkbVn9eui9q+XXVX26/BroW2XrAZf72g+6Trfz\nIfL1OdjNroBLDySkg3tfB70GP8eKYnpagkfhpvxt+7TOsBZMykphq6w1Ty1vejSlozi8a19T4by8\nL3++jTp3Fupgl/LWVQjT7cmrF1S6Bn9Ybkkak7jcB2nLr/QBv2mXQx5MuvsBYhmWLMEu86aA5y/Z\nkjjdNfDXK3dvnj4SoHYZ2y58Gf79W+uGWwL7dH215Wply2Av5bXL0nbn1L6s8nk6rGdpHGpbtfPR\na1dC4YMrLtFee1DSrUcEvQS+LHFLTsFfB3Bvudzr9tEua78fAvg67D1HNl1feT5dwfS5+Hk8d5l7\nlmnfAW3ww3pKlU25AipVMmVb4+If7COWeah1nJD7APMqD7GucCok6BxPpyuXp7Wm6jLeA++S+Pzl\nlEfVf6tP592v9ba7lzQsWkpenr+n6lsXbwEv4/MqX46XtjXXnNB34BrQgx280266+wnrGuR18HN8\nNejJpU9eBGVzpzgm8f1DP29zKJUrj8MAv1u1s6ZtepHgp31qX8Vwd+0CeLnMrQNZqOPL7Yq59Jj8\nTy3uOXXXyEqkpRsvB9Votz4sDbFNGe+BPxxRe75Sus/6b8+LsdrRrlP3NabX3wNvuyw/a3Vlz+ed\nwl6uAGrbXuZ57McOPrRWf/lGgwssBb/UVMgVXq4LMR62mC5oC9A5eJcsu8Rat+au6+41WSWmrZfh\nv5z9aMM9l26XAbWznkO+DPbe7aa43o915/igH7EEJPTIVNjNpVvXy9Yt8+QnM9NcrJZo1fu7p2t5\nS6x269VukX2bPkN6y9N5Ls6WQt6btwS8Oux5uGbb9e2vt4N/0y6lAY14yYXWefWbK+/qkxqf+xT1\nW7UH2H3N0wfJ/K3p1n1x1vJbSnt10bYGoiX5deD64G4peiu/3C+it73crkCnHZC0VoOPmOrvuNMV\niES8DPvU10hr0lsowxuW75m391JN55zmlPXgcmyqevmeXC3o+8uAMmh90Mv1z6l7zz7p7e1qV+Cb\ndkCtJqu59GXY5Rqnrf5pb4F2/lml832Yy1uT32u6iSNLXH6evnib9ysuE/jW9pbuxzK1zcGu9WXs\nYx+m+7PODt5pF6x0iCXQS/lh+WntywLtNuxTdc+3WbI1gNfK1gFS6iK7DOh7GxIXA33rCHetaOrA\nAz3Q6/g+9inf3m52JR7LAeVbN9hc513ruTxnl6Z0wXSPQnnbNdu1vGT9N0ifI7hfm1fAy7K1bnKw\n0rWZP6P1TrVyeh/nZX/X84q04dvWOmlzj+nqDm/qNajBHtbfY2vn280VvWzHubQHwQ63J/s6Cz2P\nfOXc07LdmzXl+0iuZzf4m8A/+uij+PKXv4y77roL3/72twEAP/rRj/Abv/EbeO655/DAAw/gH/7h\nH/Da1752xabrPZB9S7dVX1rN8aZm+XRbS2zpMks7nsIcbluHs6sAvLQ5cdjPNtIay2XLzoXcr3Lf\nlNxa711bNv0VqMw+8pGP4KmnnsryHn/8cTzyyCP43ve+h3e+8514/PHHV20YcLucH07uCslpuiw1\nJwsTJ65MtnMaMXTPK5dpTXqe0nI967ka0+YK7EN+HueuydLrWb6fCLZ4r7XvTX1f1/I1Fzkv66yp\n8A899BCeffbZLO9LX/oSvvrVrwIAPvzhD+Phhx9eAf30STXQrh2XuEfr1HKZ7XNE25oBJHoNbp8u\nz66asvfYRQyEalmt0zml60+F0rKy41k/QF6+74vb8C+99BLuvvtuAMDdd9+Nl156afFG64MI5sGf\nX/eu8BzG7tT9vhNszYjIyzH9cJqLoMtfVr7QNvycERGI6jfi1x77eozf//DP4/6H79drKMC+voez\ntexVBKg9bk7v3+H392LtYqHbN/Rrn760HuXm4y3y/+mhsow7e+7p5/Dc0//RtQ+Lgb/77rvxgx/8\nAD/7sz+LF198EXfddVd13l9+7KHOteawL+3x7AWnv4NvfUWz1JZUUtfWtqWPSJdC3+5cW74vuYbn\ncenaS+hLdv/D92di+rXP/Gt1u4uBf9/73ocnnngCn/70p/HEE0/g/e9//9JVRONCvDVCaUl6d+jT\nfpTLd3sGXJvvGvq6rXW9cz9yHfStR3at5/lL9jl/+4PiWqSa542A5dYE/kMf+hC++tWv4oc//CHu\nu+8+/Mmf/An+4A/+AL/+67+Oz33uc/Gx3DqjGLIKAan2fQo4F18C/b6fty4pvwY+2VJYlpS3KoG5\n+XWevEa1cRZzFUkpLJUm5NdVfk3gn3zyyWL+v/zLv6zamLPSDUxZqCHtBXgJ9KV0yi/v7z5GVfU8\nQXi1Q98D+poRjkvUvWfeEuhrlD2tI+CcQ1+Hezn4B315Rp6WFEo4qZK/n0pB7kd5/1JZC+6ljwFL\nHXLcKH81Ab/mPYQlZfJqlparQa/jGuxSer5SygGXDnuu7QAV51puB/mIZTkvh7sXcDlva75antx+\neT+nwNf6GZaAv0v/xCvVWm7vrvlr1HwO/hr4c/FkZcDl/LJnPmwdcZ7l4F+Bt+XkLtfhbYG+NNTx\n0r7MXep8vdO8fN3LoS95Aa9sm4e0lb+rCy/TNbhDmHes9fayz1t5PW4tU/DX2QE/cZWntUtbA1vH\nd4M+35epikMtU7r8chldeeU2Vxm8GlW9ZNMrcPGQ6/SSsKzMfaC3+gJSPMG+izsPHPxtOa2k8+Bf\nJPQ1Jy/vT1gC/f7TrwbbFdq5eXrKWqFW+Tno5+DX/QAhL2wxV/k0xxq1P9g37crus8svtctLky4L\n6VpYdsF1POyD3N8a8PuH/tUO/BK496Xmpfgc6KW8taqu9yG59S4V/u/jTrgC78Onwyq5wK2pNE/I\n02GPypcufz7/HPR94C6rAMJ5ejVYu1G0FvZS9d0qq1/dnip/vhLoV3yZp++RdVJw0Mdy8pIgnurl\nk1tfGb7pPKE8bLcFfB7O3Q5AP9hrK4BXsu3qhq+dbxrvB78Gehn2kNvnik8Rr5X125Vow+upf/55\n6NtwhnQpLKv7LtD3VwZufbV0y65a5dB7c89Vc/sBWa+3FE9XuBdyDXpYi9x+Kl9r+7myBwd+anNq\nXvcEgKWVQH9Yhr4O+1Loa/0KpXTJrrbr33OzzrvapbK5+JoKoHzH1MtK0JdgX2dhWzTJW2NXAHhS\nU0n125VAC/rppZHpEJ8L5y850IZ99wpAnq+pXTVl11a/4fOGXci7KMDbSt9TrdeBLyt6vwufW9in\nsGfyWAhroT94p13ZpddKvgz68uWQ7fl58Gvz5bdB2asA1sIfjn2KQQmLVv7Vs5oPMufPzMO/PK9+\nJcqwt9Otu3CtlRoEIc5YD/0BgS+dEnnK02mVl2Ep+DnkpUphaRdNuUJBti3EdCnM88Jxt2/9NvBX\nHfZgJT+kDXitemwB3VdWuzpaFqZ3ENCG3cD2nIyqEfSXb/Icl7vuqh/45RmZTqDXwF8H/bzSz8f7\nQG/dQvI4a3npvEwhLzu+0zy9vkPatEovV2kpb72695TtB/ayfEjILcxOFUAZ9hTegW341g2ag68h\nrU02fi20r1IoX9KerplSP8E87D2VANADvC6/mrAHK0G/Dvgy/GtDDXkoK1XhPbC3ILcw3dBLLS+P\n19sN+gOOtAvx/JLyJF6HdAr7/GQLl3W+O6aU16fwdcDlsbr15+cmR6WldeV5ro7pI3F5Ka7nmYN/\nV9jzsHbVSndGGXiIvAC2hjxUCC3wSexFQlvDvpsd+G252uXM694lyq3B7lf7ucvaVvfdVT4ce7Bl\nwF/ttnzpSFPcpadHtF91z6vZKezlBlt+57XVvRTXgLeQlaPl2eOfQx/2MMWXKv0VeCzXcq5qp7kO\ntvWnewn85fWjud1WBTA9pn7gZXweeD8/p7heF+JNpMOLLys2Nuhige+BfnpVen3JurobMGyMJ9iX\nKL0+4gC6O6Nh62kuedZ77aC99EDNoepRdZOF2q23ar4exe8BvAb5cuD97cJzsKs4p7icX8evjqnb\nmFMcANxXzmtHUYCedlH2aVxeybDMGlkALEhAXq6MbPypp7b6c9yCBD0HfF07/kr00tcgqYFu4+lL\np5GL8eXQr60EuoFnnQdktzr3wd5y6dd15+zfygqt9pRnoKdyhRFBopq7XlJ+DXztjuu/6rXtlVXe\nwKLezk9HmgAPU3Lhp0q/xA7wiavabaAB6TndJgI+BV1DbhTwy34DbCnsk/gEdF0JAGXwBdzcVnUW\n81wNY6/IKd2uBAqgs5i3pO4swS+pd7kdnkMOldZ32XLga6F282U87ElSdT2lSmDtVT54Gz5cghBv\nwV6Dt6TwU+in8fY21lUK+XH4OJeOrwA7y1sVU7VvKbxox18phQ/Axju0BLrKJ/ZKnpajcKIm0Psz\nylwv6wQcQNE/ZDGPviPStiDWMVX3WiiPXO9ZAH+fdqDHcro7R4KfTmUOvz7V/dDPAa9d/2VTvh/h\n+GKc22XxfAjwI7wyD1AKL+cV8/j41TBxpaMSu3yg0maPTdMQ9+ugBH2u5lrtOW6jDHsrX8JNxbtE\nLwNMAS/tWxl6ggFi6PYsvUCbpsTGrnbANrw8CHkJULkUNdhqrrwOp/DXKgBdgaydwPlxhOMOlUA8\nDwFsblcEWRwVVc+8gsNacM1jHEBfo84tENz1qPgSenB0+0uqPtXLfujDHcNZOvSZ53dGOtY24NM8\n18WX+vZlFcXZXrqjkPCvt4N/tbYEeV4JaBBLedMeew14L/D7rAAATNNa8TXkHmIdyrJw7iaKf4Vg\nD+ZAFdD7dBf04Kj4RAV1j+l8uSn4yyYGw4BQVvW0zpKiz+Vp2F1cAz+tknJ1Xw/9wdvwztqQ53Wv\nKcJnM9hrKt9276eP9lpt/XYFAJ7eapO8zHWvhErtXRyot+/D+bwKFlxtRNrLoAfFFm3xuTCuW6b9\nOpmFW798Ml7dHfRa1dN6436LUMZL6u7iULA7+CXw02orj6+1g78eO3Xl+y5TAFsDXoa8BnztcV69\nEugFH/o4uJzfhLyi/inu0pN2/FWReRL+CIv2e4y7NMXZCxC1YA93S/Qa5BkO4RrYXfs9gU5ZHFH5\ne+DO40FoUIBdAl9qnEgr5fXYwV+PlfViXgEsu1TlTrv5ymAf4JfKMHMMAfAm7IXKAOGcSXee0/ms\nw06Nwgsqi4CmHWsqvAyDCy+gj+58o0deq6+7o+ojMPTDWwd7UvgE+nTdMuyNA0Hhw5RgD4qfq/y0\nCtvFDva2nAZdt+GXgK6hbgHfgrxVETiF1s/yp7dN3Hfpvmft9pAOcHuFDi465+kAvlTyKvBXRdmV\nZT3zFPJ4Uk4EpA46Tvmy3R/An7TnOeudz+8mql41hnXlVFN4UmGqFOK2w3GIeEnZy/dvErr0vh1U\nuB/YgYO/LSfr5nll1268hrbl5leVngugcwn4wmg9TrBruAOAeR5y+EuKXlN4OaCG7xzoiYA4GjzW\nXbFR77swGETu+hMxmHzYSkOk/TojFrIiIMDA5FeOgppzqgxiXmq3u30KsHvFJ84AdnFklQAK6bIx\n8pdqy6reii+1g/fSl0DXz9d74K5Npfk17Bn4PAWfY34BdKYshIxnKo9iWRl6N29elvJSGsg6764Y\n9BFGRL8eWQM+ywvpEHche5Aj9NIbUGHcGuVhqFAk2C5uy8CTV3FmMPlccpATWzC5HvwEvWvX+70v\nwK8bPr4ChPFlcmjPFP59wQ4c/Dl8Dr1uu9QUfin4xXlZqTmXwbdsIEGNsCvQUz7E/MhdcwU+gOS+\nZ5AjzRuVHhHouJzIy9NXw1iA7qDNlT3MJSsA9q48R5e+kRfXKSH34Ic2Pvk2vAaevLrLdAH0iUtP\nBGYXYgJ1Dn+wUJelcXUOdF+FIDQupuKX1qDBX2sHAF7XVnWFb3Wg7TplFYiAuxRnJgGsgl1WABOo\nkbfNodvtyONi3gzymB9OnE77ILj9VwX6yHFQ5jbsYnQNEF15EYcHnDTYCnzVT0BkvJp7dz9CPgWe\nQGBKgDvFJzCHZgdnZSWVl8ZwnXRTjU4NAt3DkPq1NPj5yV1TARzo9+EpS7dUXcJe62lfPSnAY1qq\nPQuYWbjxcTIeQBLzJXe92imXQS/yW2VABnnaRjyRV4b1YA5GyjIIod2d5yc3XoCOHPoJ7CIuVT1T\neHKwhvZ9mkjAn/LctVSAk4VhDzqRgB7NyYjDyyGfhtq1R2Nae50P/GOSpdqsDnsf4APGpSqvYec8\njzPoKcuDht8f4AR83WaX4EplLqXFOsQJzNv1V9ClB5C/vBd75mWbHvEQojpH6FMFEOAOlUCoIGKd\nEdIxLtNCzYPCF/O8y07kmwxC5UG+wy7lGeT3cH5QuSXgka2h1n6fq0jW2sF/W648TdvqffD3wy5B\nbwNvwJYU8B52K0CX6i6AjnEFcmyfi5Mh3fkMajVfjMfykHc1XXoE1YX//CIV5hEKL6F3qowM/Ah6\ncqJiXly3cP0l8DAsIKdM8WEIsp8+9aDnoIfe+zqIU/gT8KGqcKi7bkMNfpqzpzJZYldipF1N6SX0\nc+69hWnCPmKYwi7Vng2sBN4m6MEUoYdVym6NiFMGnWybxzRSOs8T+aiUxbjIh4Z+rxdoP+ZluAo7\n4F12EpBCgC/cdA0/pmWx8vDbICPc9/AcnxhkKKm9YVeBEzuVN0nFgxsfR9yx9wIqIJYuQQ588BlK\nsJfa8dorXm9XYqRdzZ2fA3xumkBeyJNt9gnsVgNvoqpDpR3sJgdZAT/JhygD8nmg5wvxggtfAv+q\nmPDey7An9QcBMNqFzycJf1yPn1d6AHIZsgwY13EHFoAzudC340EENl7jrYiHnnkE0G1sm/sjqMKf\no8qVv7kOu5C3++W9EiPttAtTg16DvwT+DHYJeMG1z6ag4NZEBQ9wsyVRAYRyzEwllz2li/FmBdCI\nH9pkC4PyvBSnFJ/APc2LlYDsDQsdeiFPlhn4kXnu0Rw4wG6dwsMrPlyaLYONV3b2Kh86ZUPHHWlA\nWyaRDc/9y396EE5YMm8q6JO4zK7ESDsNfU3pl4AulbwWL8Iu1d3H2TolhyUBvUlpofYTwC2meTWQ\n10J/FWGXVlL2Wl54hlWZSpVAXM76uBX5FlHNg0vvoPewM7t2PVza+AE2bCj1yjP883cPPeAqeVmh\nVQ893eES6eDIp7/8Md90ad1qWwf9wV36kqK3gB4LkwR55AD0gDG200OZj0/AJlhrvIqHjjqTVF0q\ne4QcHmhfbtGGu1aGRnqXsn1Z+27usyXAd0BfhL1VSRg4wA3HedkQyPgyAtgYkLFgYwBjwWTAxq2b\nDcEaC0OANQSY0BQQO629FXFsAWDpv+Z4W7jPYoTKwKqc3PsF1l/mA3+1VrsurfeaXJmDN59GP1kW\noAu4S2VsE+wBbuuVOgNdKXmAPwNfwr4r9NhD2e4XZp3VKoAS3DJ+kcAbUW5IpeGBhoOfHOxkjOu0\nMwYwgDUMYxzsrj8gAU9EYv99XKh/OC4CYxQqn16UJYG26yeowz53ouftwCPtwiTBNpAj7fLn5nXF\njyrPBpYD2Br0IXPdOSq7Dv2juALsXFJ0DX0L7lI+KvHLBH6fnoG+Hy8S+Nq8RsWjsju4g3ITIYMe\nHnYYE0GHYQe7YRhjYZniK7oR8hgWjgtS2QcBfA5/wNyI6kC263OBXGcHGWk3fb6o1byvR36UKu9h\ndtq1Zd8AABwRSURBVKEC304rAh4pA52j0geoTZwnws4J9EzxS8DX4L4TgS+VzYnMHOSl+AI3vTmf\nzjNyIifAHmQNPTz0zv1n8IAYl+BTbB6E5RX4ohuf/fFJ0ENowJkDn8PuRM75BvoyyJql367IwJt6\nO74G/Qip8IPogRdqHkC3EvxBQU7OlR/zPB5rql5Q+axNjzbgVxH42jI966rdc4cCvpaO7jwXoYdv\ny8MAGBL81g6ggSPkZNil3bvPbiEDJwayqRDOHYWA/f2aAy/D0HZ3+Bvolr5279dc6oMPvEkDC+qw\nl6EPo+pEDzwn6DPIJfh2yJQ8qPgkLtISag2/y0PbpW9BD+wO+EUp/Frgay79HPgXAbp22Yli513o\nuEvz+DLroOag4r4nnwXsFL0BcufJ+MOIn7Px3WsMhOf8yZ0PgKcw1/n80VzuEe/W+jro67Fux6fQ\nl2Evj6grK3xS8gD8KOJB4TES2CJz3VnkQSv7iPTobcyhz1R+DvI7AfhWBVADupZXg7uUtw/gJ3GK\nyg7DUZXzisCnB7j5LANDAt7BDvdSjY+7891wqwcgfRGAhbo7rEMXXVJ1G+/zUu98gn+9NYF/9NFH\n8eUvfxl33XUXvv3tbwMAHnvsMfzt3/4tXve61wEA/uzP/gzvec97Fm20/Jyx/XhOd9gldR9Ex11o\nvyu33ib4x9ErvAc2wj5iUgkkyDEFO7r3EGn0qfxVA34O9tI656BfC3wv9E3AVTrrsBPwh/Hzoo0f\nRthhIK/sAXIRDgxiXznIcx+Ph92y2fWQ6i6Bd5hrV3764yjpRO0CfRP4j3zkI/i93/s9/PZv/3Y6\nFiJ88pOfxCc/+cnVG/VrKsJe7qHX6UGouwCffSde1oZ3brwdPfQ+DGBjlLAjVgRZXHfMZRWBKq/B\n3VJ/VOJry5ZaC3Adr4Gs00viFwW8evyGMGw3uPKTR3bsKoWBwMzuug7wCu85Ng5kGpBc9nAAYbvW\nb4fDDMGll/cwC8BTe73nV5CA9dA3gX/ooYfw7LPPTvJ5D99Qqrnz5QMt9dLL0XP5sNkxU/Y0jeMA\nO/qRc6OHNrjpPs4qPYUbZdDDPD2gX2Xge+Jrwa7F9wV8Ky97Bs9e1SmDngwBA7vhzh52ZvjReC6O\nIZ0PE9KE6QhAdS2mwNvCPW28Yz/9fmISSH0Sl9mqNvxnP/tZ/P3f/z0efPBB/MVf/AVe+9rXrtp4\nundlO6VvfLxsu4/Zc3bp1hfA99BHSPUkVD/L01O1EsA86HcC8DXYg61x2efiNaCXAl95Bp+Unsrq\nb5Aq8dhG96Ab5Oc+gJ8dM/tmQljIhczOJQg6nmDX0LcHnoWD2kVuFwP/8Y9/HH/0R38EAPjDP/xD\nfOpTn8LnPve54rxfe+zrMX7/wz+P+x++v3owtZ+Q0q+xupFzor0ensUHmLcOaDsa96ht6zviwrSF\ngjtMBchLwLdAH9EH+VV/NLfGpd8F9BAuUfWlKj9x4St5A8rXtZAvv5PgvmrMsAJ0ef6IgdFruKXB\ng8+wJJQ/5luEfquaSy9P4nNPP4fnnv4P9Nhi4O+6664Y/9jHPoZf/dVfrc77y489NLM2qewNRY9v\ntQnIs+GzAXgD3nrQI+zOfcfWAx3CDGqqKD1ywOdgX6Pu+4S8pM5zpudd4tLv6sLrvN6pp83eVHjU\nK4DStS1VAvE6Upzybxyyv9zuxJFhWBowkm+xkwCehNJTiYHpY7rkGQP3P3w/7n/4/nhKv/aZf0XN\nFgP/4osv4p577gEAfOELX8Bb3/rWpasAALXzudJnoGeuuim47h76MQGfgb+Vyk6AVPga4KWyJW79\nEpUHlsPcUwaRN38xlqWBi4X+AjvtZuEPal6aKk23ADrBKz3DK72vQ5hAQ4CbMcYPZlr/gc0BhpK7\nP9KAwfsD+Y+gStd+nTWB/9CHPoSvfvWr+OEPf4j77rsPn/nMZ/D000/jW9/6FogIb3jDG/DXf/3X\nizcqexqnbn2hcy57fXXaC2/ZAc+jgd0G6Am8dcrOXuGjO79Fv7q/0lz6WlkpvwT+vt14mbfEnV8K\nfmuS89TUvFixe9CzD5q6IucsuB0gP/beEsMaoeomufZjgB8D5Ii78q8dJYaWWhP4J598cpL36KOP\nrtqQtryjTtZkBXfeT2N4ni4H1XBqvzslD9BTrvIe/qTwJde+MC1x52uQt1Qf2C/oSxS+B/JS3hr1\nnosvneaevc+pfG3qhV1cU/ceTd6x5lrrbrYIvHFt/Ax2YlhjXfudbdGlZ6T2POJ21tmB3paTsOfQ\np0kCn78Q456lD3Go7JgB70EXISTsWyqregv+Je78mnY8KvFd1H0tzHP5wS4D+qWufS29ZAou/Qbd\nwMvfA0yfOWNYXxGEO9yyn4yN7rs14f62otMuufOlzuzpyeu3A70tJ6cc+vQb8KGDzj+CU0Nj43P1\nkOddeanmXJgc9OhT9ldCp92u4OuyfbnvpfKlyj6Xv1TdtUsf2u41+BnZV4nTMFvXY59+2cYCvBEv\n3lg3Hh9pmBnR4PMNDA2Zsucdd4jhGugP/MWbdABa2eN4edlxFz9akTrpZFyDjvMQRwSdJfA9nXZL\n3PldXfp9q720XV19aRcJfY+qL4W/BXevS1+A3w3KAdLv+7mH9Szu7PB9ezDBDu6BmxOyoOjCnYdP\nI39fpPw8fv4ylezAH7GcuvXh+XsagDMo6AX4o3Dvo/uOCHuMbwl8DuHWY1lPfS/sV9Wl3wX+ml20\nS78P8OfUfURd3avuu8wL4LF33OXBOPBB6W5PQPt2PAnYWUI/Iv8mRGLjQjvtLsoS7NKlL7wGKx/D\nRXV3oI8B9jBtjYAdAnIXx7lXd18RrHLpe9S+B/KSW49Gek0ZFsRLVivfp6LreXeFe4nCl9RdPpIr\nufE6n91pIiZx6sOPbfi7PP6CLlz7PfzR4IH3bXoWyj9R+AQ7YrjODvwBjHyUXXXQjXrzLb4ME2FP\nwEfXPcAeAI9pTBW+B/zLcOll/CpAX7OLgn6Jqi+pBEqAl9Rdu/Q18EMeu51Pr5YEZYfHHtnbs5my\njy4+GguyBmQGpO/cDIKBacfdLnaATruCG6+Hz8Jgi414rz286WYmnXNxkmBruEt5F+3SXyTwPfPJ\nsDfeY4dW+DVuf03RS8CXBtvovBE59Nn1Jaf8yoMLP0+WRogSLBkQ+SG11noxszAUXHvXpI0dd6TB\nXw7/4RSe5aM457qHF2EC7Fu7wda772MYI781sFsCnxvwueh511DPxeeUvQf4Hpd+DnxgGchL4uiM\nr7ElMMt4q3zt1KPwLVdelwdlD3APlfyee4NTnGNlkMackG/XjzT4n7ty38YfKX3UZQyws3DtiVZd\nwsO04f0veVhx4Jb9hyxYQB+es/uXYcYwqObcwJ57VQ8Qt9S8V+HXDrxZ0mHXAn6Jq75G3S8C+H2C\nf1HQ9yi7Br4FupxnrmmnpjAMN/0GjcFIAxA/jMkgM0QOxmwgjnfvyTcjVnj3lw98+NmeCL18QcYf\nKDvgt3bj1D100m0H2C1F2PmcwGe+c26Nwu/q0rcUvgZ+CX5U4nNAt/J0vs7T8SW2T1c+xPcBfQ38\npc/hNeSlCqAxtr6Ux1Y1Z8m9KOM+jGlAxmAcBrBlkN1gSwM24SkVxHP5+As469ryBxl440LybXeK\nvfHBrd+yc+mTwqdXXmMb3sOOcwBnWK70vQq/BPgW1LUyYB7wHiUvlfXG19hluvRz5bsofAt43Wtv\nkFcCNdC5kgd/3xO5z12TgfXAYxiAEeCBYWzoxfffeSA/AMfDHn6NV/7cVa8dzqX3bg2rN+KiS88b\nbO0A61U+9MRzbMM72PmMgDNy0Peqew34lvI33LSmS9+TRiXeC3ornMtba5fp0rfKeqBfC3wJ/ppL\n33E/sD8OJg89ESBgx+A8AWLGNnyyjVKPvVN4Ez3kO8Olh3Bt2LXj48sxLFx63mC0G/H4LXTYmdRh\nd+5BP8XVU/gS2LU8VOK7AD+Xp+NL7CJc+hDuAvwSl76WXwO+Vgm0PDw5+Z1xsDtlh2G3ni27T2t5\n15+sxUgb33GXt+HjsNp1Hv1hHsuFkOWbcHFyvfMhDD/9FF9/jR12xrXfz+Cgryn8XKddD/wdPbCz\nLn0tb41bXwp7y1rxHtuX+16atwZ3KW8J+HOQ9zyeq70fXwO8cC+4d2ucO28MgY2BHeBgHwEeCTYA\nb/TvLlBy67H+efxBXfr4lpx06a1U+EF8qko+gxe99MGdD+34XTrt5lz6Hte+BHQtf6nKl8LeslZ8\nje3TlQ/xGuA6vRT8OdhlWa+660E6cxMA94s35EIPOw0AbwDjP5dO1r1KO1rvzhsx6o4SN6B10B9w\nLH1w54VL72EfbXLpOfwM1CgG3Mge+jM4lT/F/hS+NBKvBTqjDHwP6Lu68nN5PeVrbB/ueynvIoBf\n4tYvVfjOCp8RtkHAQLADgbYE2hjQ1kFPowc+KrwfVWp8+z08h7+THsvFwf/ysZx26Tk9g3ewE+A/\nbuEG2hjfaecV/hS7K3wrr+a6917wmuoDU6h3Ab4W37fK79OV18tdpMq3QA95NTXvbcPXJoKHncED\nAYNxkG8ZtAFg/U9RWwKsdbAb7877nnrLrv1/Rz2WczcxAUz577jF32gXL8icG+DcAGcUe+P5jIBT\nAp960OU0p+oybLbbuU/htbLPAs99wLfi2FN8F4VfC/RcvKX0Mj470bzK19z9yfN2yoHXg3N6YA+i\nQYi/Viu/i8/GuN+yG9i5+MapepjCsFw2vsMu/J7dCpU/CPDx55mCcqswxk8NcELg2wScEHDbTycA\nbvspxLXCb1GGPXssx4UBN4U85jLYllUaZbhbrjuHk4IpjKW82Th3zleI9xiJZZbG50w+apI3c4gz\npjD3VAZhuXCNMsgprxQGIH6bfgTc78tRDntIlzxCPZ0DOEJ2b7hHb/54Qx67vPgzaEM+sUwHW6Hy\nh2nDM3KXOX6YguIHJ7H1Sh4hL4EvwqDwJbh1egsHdrioVkJeyA8Aa9Aj5DyFHTLdyIcMeYEL3wH2\nkrJeK4Eo19UCvBd+CWgJ2ibsNM0Py1hKnW2GhKvtl5NlI3LQjUz7SmADB3O4l2Vcwh48yuzrOFDw\nJ+hhCbzxcG8IvIHblrxn7iiX3iL+cKN7fTUPY8fcSYKaTyiqOgtlj3Gp8D1T9DKCeotwVGkNeQRc\ngi9VWqn7RMEbih4SVeBFpQBd1lNeSC+1Gtxz6VnTNzGl7RGSByChB3LVR6Eszh8gJ/HYjqbQG0L8\npdmRkLf3RToIVgAywL4RYYB+JPdjlOERHWvQ3b4xw4neEYGPSPxisZ/gXPsoMAvtcG34Ee5gwjvr\nZ5R3xJ2bAuwhTmWXvgf0sE0NswRcxwPktVDCjwXxcD5kKIEtlWXpwjL6XLfStbyWSSZZ5dXSIa+1\nLRbzSEjDC+Xk/2nI5fYCwDHu84Oyx4pBgE8qzH5SWsAtflc+VgwB8gD8thIPQpSpOwlFzysA+N8+\nhAc//XAd0i/dDjPns2IHU3hYIH5/7hwJ+tPwuE3CTcAJknvv4c/a8gF43ZaqpcN+sAA8xjnPD9Xw\nBPhCHkphWEchTwTTZTCdN4uqK176kc/aTbFW5TXUel0lsHvzMkWnFOq8AH4GvK4AVOUwAbwVFx5A\nFqp8CbWeZFlo7wvYwz3DljKld6HozwqqDsRn+Nmgn4V2mM9Uxw47JDfe9767ybhQtNM5g5y8O09J\n5aXCq/6BSbw4EkpALkM5QcwXIS6kfTQpbyGU84ggbiPGJxGxnkLZnIruyzTktW3MKfvEAuA0VfgY\np7ROqejNtIA5S/uZqTAZAsjkfQSyKRBeqglQa+hlngee/TfqY3teqzvDeSSjcU+vkAbauE47OHbs\nnfI+vIQsc+nhx8V7yE8IuG3AUcW1slPu1p9i+vy8FsqOuNZkrah2/c5H+AuVwARoVnEU5hHxLKor\ngk7ge+HvKa8ZqXQJ8jXbYAEgtMJr2EOFgJSPRjzCLkMU8sIUQLciLjwLgh8Hj6mqlyoCqfAICu9D\nuPY9/M9VhV+0SfsDB3tYv81vrSV2WJc+9Mh7l55PnRvPJwS+bZKC34avAJAru6wAzjCFujWiLrrY\nmELMdponw1LeBO4S4AX4Y6CvnoS6BnQlv2b7UvgW4D3baam+eP0zQh8IkyovXfYQIREH8vWUwEYl\nj0wCHXDyS/7ASKw3wLwVcf0ufQi3iICHH60gZoRPYsnbKsIe9im48b7XPvNQF9pBO+3S44s0sCa6\n8ll7HcKth+iso7wN3z02XsCUtbc97BFgWZWWgLeFMqg8kQ5xDbyMz0GOwrJLwd+39W4y8FedXyi7\nhLyUH9cnoI6rkTVQAfAa8DBiWQ89jPe3SYRplgzy0oczwhQ77YK6e7BZufMMuJ45P4Q2PIPfUHrs\nl3kLy+yACk9+uKx4CSZTeAK/TJj0xmeVgJhmR8+JKY5y00orz3oAf2aAPOt8qPlblYAMe+KtPFV2\nAO6jacUP1twnAbIEk1U624CsEPTGC+vScOsycFJ1wM8jiZSVClLHXe+rtPGWEKpuk+LHWwT+0Zsf\nd48Nxe83sm+/x3t4oV0o8ONZOW88BcZThj1l8ImfblvY2wzctm562QIvU3lUnUyHSQNvUQY9KHzW\ntmAVluIV4Kt5JfhraXTGFwB/SFu1+QqgkzzkcdbLylDOW3tRXqfDULzSWzWqLLraYpbaGHwNaE0b\nGPAvxYPJgk2Y3NBbOzDGDcMOgN0A49A8qRO7UOB//NK0qr99m3H7ZYuT2xYnL29x+vI5zl8+w/bl\nDcbbG4wvG9iXjVfyIQF9KkL5Savam22aVYhwAngrbAFdgh4qH5X4EtB1urcyuJOsA+6ayhfXIfNq\nU60S6H3Fjuq3gt4V+Wj4TKwimNAXHrew9mXY8TbG7Qm252c4OzvD5vQcw8mI4cSCbjHwY8DeqrlT\nZbtg4Kd5JycO+pOTEae3R5zdPsfZ7TOc395gPBlgb3t3/jYDp8P0BRkNvQS+BnsR+hLQSyFvTVBh\nLY5CvFW2pqK4E2wJ4CXI9Twyb+nUertGxr1p2PUuMJL7f64W187mCGA7gseXMW5vY3t+gvOzU2xO\nz3B2ssVwewtzewRuWuAWYG9ikV0s8P9vWvucngInJxYnJxanp1ucnWxxfnKG7cmA8YRgT8k9ijth\n9zw+vOseH9uhrPAl0GvQA5XCVh4KZTXIa+72GtiXzHsngi6tpt69yl5Kt5oHvdDXQmUl8IOFTuSw\nuF5GPqbejuDxNuz2BOP5KbZnpzg/Pcdwcg5zewt62YJvMvgmY7xxxRX+7Ixxeso4PR1xejbi7PQc\n56cG21OD8RSwZwCfWuDUAmcmwV0LSwqvRXtiuyh3r9u+JERnWuf1zH8nWQv4UqiXK6XnmgqlPFOJ\nq7R+GaZ0r4UygrtX5eaksosXbvh8hD0/gT07xXh6iu3JKc5vnMHcOAfdGME3RtgbFvYGsD0ubLNh\nl96GPzsHzs8szs5GnJ1tcXZucH5msD0jbM8Y4xnDno/A2daNp699nLL0umvNO5+oOwqFc1D3uOkl\nGOcUfQ7qWl7vcneStdzzXkWvpXuaDEuUX2xD3wpWxI0o64AdRwDOLPjsFPb0DOPxGc6PT0HHZ8Dx\nFny8hT22sMeM8ZhxfnSVFP7/TfPOz9lPFufnW2y3hPNzwvacMZ5b2PMRfL4Fb88d8LUx8ZN32zHv\nnUceWso8V4YF4VweZvL2Me+dZHNtcx0vpVvz9DQR5oBX8+jLqzvsAvRyE9KFDy/fnEMMy7XgozOM\nR+fYHp2Djs5BR2fA0Tn4aIQ9GjEeWZwfAUdHhcNv2KUr/HYLbLcW2+2I7ZawHQnbLTBuLcbtCLvd\ngrdbYLtJ4+27BtNgXrAzmwN5CcQ9LvoS9/siyu4E63HNW/P3zDvnOcw1AXTcmz71VsxW8gDCs/nC\nc3weGHZzDrvZYtxsQZstsNmCY57F+Yax2QCbzRVS+P+v0Ia3I2McrZ9GjBYxbset/xz1uXvAaGna\nITeqMMR7xDnmrYW4F+i1bnrL5ua902EPtkbB16yn1R+woC+hdEuU6grpeYbPaOnH/DFuwcMIa0aM\nwwgMI+zg4tthi2EYMQwWw8AwV/05PFvAWoa1ow8tLI/pW3Z+YmsA/1GA7kfm6AjzvanE58qW5PWU\nLZnnIpa9StaCeZma9VcWvZ1+Ki/ALe8xEqFcLIzZCd0AWzQeADDYuN+OZ7KwxsL4tIkTwxiAzBVS\n+JcLbXhmBrP1LwmEOOVTqQe0t5mdbWwmPVtwUZC+UuA8tC2tAJYs39lcaEGvF9Ht+0a3gCUGEYOJ\nY5zI+lBOsweZ2aU/h3fG8K+sXdu13fm2pHtmwSr7VrOM+MLogWu7tmt7pdoBgH/28jd5afbsoXfg\ngu3ZQ+/ABduzh96BC7dr4Pdqzx56By7Ynj30DlywPXvoHbhwu3bpr+3aXkV2Dfy1XduryIi59G3j\nPax45S9jXNu1XdvuVsP6wh7LXVA9cm3Xdm072LVLf23X9iqya+Cv7dpeRXZpwD/11FN4y1vegje9\n6U348z//88va7KXZAw88gLe97W14xzvegV/6pV869O7sZI8++ijuvvtuvPWtb415P/rRj/DII4/g\nzW9+M9797nfjf/7nfw64h7tZ6fgee+wx3HvvvXjHO96Bd7zjHXjqqacOuIcXZ5cC/DiO+MQnPoGn\nnnoK3/nOd/Dkk0/iu9/97mVs+tKMiPD000/jm9/8Jp555plD785O9pGPfGRywz/++ON45JFH8L3v\nfQ/vfOc78fjjjx9o73a30vERET75yU/im9/8Jr75zW/iPe95z4H27mLtUoB/5pln8MY3vhEPPPAA\njo6O8MEPfhBf/OIXL2PTl2qvlI7Khx56CD/90z+d5X3pS1/Chz/8YQDAhz/8YfzjP/7jIXZtL1Y6\nPuCVc/1adinAv/DCC7jvvvti+t5778ULL7xwGZu+NCMivOtd78KDDz6Iv/mbvzn07uzdXnrpJdx9\n990AgLvvvhsvvVT42MEdbp/97Gfx9re/HR/96Efv6CZLyy4F+FfDM/lvfOMb+OY3v4l//ud/xl/9\n1V/h61//+qF36cKMiF5x1/TjH/84/v3f/x3f+ta3cM899+BTn/rUoXfpQuxSgH/961+P559/Pqaf\nf/553HvvvZex6Uuze+65BwDwute9Dh/4wAfu+Ha8trvvvhs/+MEPAAAvvvgi7rrrrgPv0X7trrvu\nihXZxz72sVfc9Qt2KcA/+OCD+P73v49nn30WZ2dn+PznP4/3ve99l7HpS7GXX34Z//u//wsA+PGP\nf4yvfOUrWQ/wK8He97734YknngAAPPHEE3j/+99/4D3ar7344osx/oUvfOEVd/2i8SXZP/3TP/Gb\n3/xm/j//5//wn/7pn17WZi/F/u3f/o3f/va389vf/nb+xV/8xTv++D74wQ/yPffcw0dHR3zvvffy\n3/3d3/F//dd/8Tvf+U5+05vexI888gj/93//96F3c7Xp4/vc5z7Hv/Vbv8Vvfetb+W1vexv/2q/9\nGv/gBz849G5eiF3YWPpru7Zru3p2PdLu2q7tVWTXwF/btb2K7Br4a7u2V5FdA39t1/Yqsmvgr+3a\nXkV2Dfy1XduryP5/z0jBdzySmPIAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Simple timing tests for pure Python",
      "",
      "This is slow ...  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%time x = calc(func=py_update, N=200)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 15.91 s, sys: 0.13 s, total: 16.04 s",
        "Wall time: 16.16 s"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Vectorize with NumPy",
      "",
      "[NumPy](http://numpy.scipy.org/) is the cornerstone of most numerical Python applications.  It provides:",
      "",
      "* A powerful N-dimensional array object",
      "* Sophisticated (broadcasting) functions for vectorized operations",
      "* Tools for integrating C/C++ and Fortran code",
      "* Useful linear algebra, Fourier transform, and random number capabilities",
      "",
      "A NumPy array is a Python object that stores a statically-typed data array and supports functional operations on the data in memory at the C-level.",
      "",
      "Operations that can be cleanly vectorized will be fast.  In this case there is a simple way to express the update step as a single vector expression."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def np_update(u, dx2, dy2):",
      "    \"\"\"Takes a time step using a numpy expressions.\"\"\"",
      "    dnr_inv = 0.5/(dx2 + dy2)",
      "",
      "    u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1]) * dy2 + ",
      "                     (u[1:-1,0:-2] + u[1:-1, 2:]) * dx2) * dnr_inv"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit calc(func=np_update, N=200)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 27.6 ms per loop"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Magic with weave",
      "",
      "[SciPy](http://scipy.org/) provides the [weave](http://www.scipy.org/Weave) module which let's you insert inline snippets of C code.  Importantly it provides two more things:",
      "",
      "- Wonderful and magical glue that transports NumPy arrays and Python values into the C-layer.  ",
      "- Macros that signicantly ease writing code that does array access.",
      "",
      "Inline code gets written out as a temporary C-file and compiled once at run time.  This get cached for future use."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "from scipy import weave",
      "from scipy.weave import converters",
      "def inline_update(u, dx2, dy2):        ",
      "    \"\"\"Takes a time step using inlined C code -- this version uses",
      "    blitz arrays.\"\"\"        ",
      "    nx, ny = u.shape",
      "    dnr_inv = 0.5 / (dx2 + dy2)",
      "",
      "    code = \"\"\"",
      "            for (int i=1; i<nx-1; ++i) {",
      "               for (int j=1; j<ny-1; ++j) {",
      "                   u(i,j) = ((u(i-1,j) + u(i+1,j))*dy2 +",
      "                             (u(i,j-1) + u(i,j+1))*dx2)*dnr_inv;",
      "               }",
      "           }",
      "           return_val = 0;",
      "           \"\"\"",
      "    err = weave.inline(code, ['u', 'dx2', 'dy2', 'dnr_inv', 'nx','ny'],",
      "                       type_converters = converters.blitz)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Inline compilation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%time u = calc(func=inline_update, N=20)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s",
        "Wall time: 0.00 s"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%time u = calc(func=inline_update, N=20)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s",
        "Wall time: 0.00 s"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Run for record"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit calc(func=inline_update, N=200)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 18 ms per loop"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Slicing and weaving",
      "",
      "Weave supports array slice syntax:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def blitz_update(u, dx2, dy2):        ",
      "    \"\"\"Takes a time step using a numpy expression that has been",
      "    blitzed using weave.\"\"\"        ",
      "    # The actual iteration",
      "    dnr_inv = 0.5/(dx2 + dy2)",
      "    expr = \"u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 + \"\\",
      "           \"(u[1:-1,0:-2] + u[1:-1, 2:])*dx2)*dnr_inv\"",
      "    weave.blitz(expr, check_size=0)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit calc(func=blitz_update, N=200)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 17.9 ms per loop"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Numexpr",
      "",
      "The [numexpr](http://code.google.com/p/numexpr/wiki/UsersGuide) package is also worth checking out.  It is similar to `weave` but with a focus on highly-optimized computation of a single expression.  Attention is paid to details of CPU caching and use of multiple cores where possible.  No separate C compilation step is required."
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Cython",
      "",
      "[Cython](http://www.cython.org) gives you the combined power of Python and C to let you:",
      "",
      "- Write Python code that calls back and forth from and to C or C++ code natively at any point.",
      "- **Easily tune readable Python code into plain C performance by adding static type declarations.**",
      "- Use combined source code level debugging to find bugs in your Python, Cython and C code.",
      "- Interact efficiently with large data sets, e.g. using multi-dimensional NumPy arrays.",
      "- Quickly build your applications within the large, mature and widely used CPython ecosystem.",
      "- Integrate natively with existing code and data from legacy, low-level or high-performance libraries and applications."
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "#### cy_laplace_slow.pyx"
     ]
    },
    {
     "cell_type": "plaintext",
     "source": [
      "cimport cython",
      "cimport numpy as np",
      "",
      "@cython.boundscheck(False)",
      "@cython.wraparound(False)",
      "def timestep(double[:, ::1] u, double dx2, double dy2):",
      "    cdef double mult = 0.5 / (dx2 + dy2)",
      "",
      "    for i in range(1, u.shape[0] - 1):",
      "        for j in range(1, u.shape[1] - 1):",
      "            u[i, j] = ((u[i+1, j] + u[i-1, j]) * dy2 + (u[i, j+1] + u[i, j-1]) * dx2) * mult",
      "",
      "    return None"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "#### Build setup.py:",
      "",
      "    from distutils.core import setup",
      "    from distutils.extension import Extension",
      "    from Cython.Distutils import build_ext",
      "",
      "    import numpy",
      "",
      "    ext = Extension(\"cy_laplace_slow\", [\"cy_laplace_slow.pyx\"],",
      "                     include_dirs = [numpy.get_include()])",
      "                    ",
      "    setup(ext_modules=[ext],",
      "          cmdclass = {'build_ext': build_ext})",
      "",
      "#### Build command",
      "",
      "    % python setup.py build_ext --inplace",
      "",
      "This creates the file `cy_laplace_slow.so` which is an *importable Python module*.  Once the build is complete then the module is installed and always available for import, unlike with `scipy.weave`.",
      "",
      "An important detail is that the first part of the build process is creating a C file `cy_laplace_slow.c`.  For packages that get distributed one normally includes the C file so that users only need a C compiler, not Cython, to build."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import cy_laplace_slow",
      "def cython_update(u, dx2, dy2):",
      "    \"\"\"Takes a time step using a cython module that",
      "    implements the loop in Cython.\"\"\"",
      "    cy_laplace_slow.timestep(u, dx2, dy2)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit calc(func=cython_update, N=200)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 400 ms per loop"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "markdown",
     "source": [
      "#### Rats, not so fast, what happened?",
      "",
      "`cython --annotate` to the rescue!  This generates an HTML page that shows where the performance bottlenecks are occurring.  You really need to get all the yellow (pure Python) *out* of inner loops."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import cy_laplace",
      "def cython_update(u, dx2, dy2):",
      "    \"\"\"Takes a time step using a cython module that",
      "    implements the loop in Cython.\"\"\"",
      "    cy_laplace.timestep(u, dx2, dy2)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit calc(func=cython_update, N=200)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 17.7 ms per loop"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "markdown",
     "source": [
      "### ctypes",
      "",
      "Yet another option for connecting Python to C code is the standard library [ctypes][1] module.  NumPy has a special [ndarray.ctypes][2] object that simplifies using `ctypes` with NumPy arrays.",
      "",
      "I use this interface for the [core integration routine][3] in the [Xija][4] thermal modeling code.  It was quite easy to set up and works reliably without any additional library dependencies.",
      "",
      "[1]: http://docs.python.org/2/library/ctypes.html",
      "[2]: http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.ctypes.html",
      "[3]: https://github.com/sot/xija/blob/master/xija/core.c",
      "[4]: http://cxc.cfa.harvard.edu/mta/ASPECT/tool_doc/xija/"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Fortran with f2py",
      "",
      "This is included within `scipy` and takes care of gluing Python / NumPy to Fortran codes.",
      "",
      "**Fortran is still the gold standard for speed!**",
      "",
      "https://github.com/certik/laplace_test",
      "",
      "https://github.com/scipy/speed/tree/master/laplace",
      "",
      "    Method        Time Rel. speed",
      "    ------------  ---- ----------",
      "    NumPy         2.03  1.00",
      "    Cython        1.25  0.61",
      "    Fortran loop  0.47  0.23",
      "    Fortran array 0.19  0.09",
      "",
      "#### Example code",
      "",
      "    subroutine for_update1(u, dx2, dy2, nx, ny)",
      "    real(8), intent(inout) :: u(nx,ny)",
      "    real(8), intent(in) :: dx2, dy2",
      "    integer, intent(in) :: nx, ny",
      "    integer :: i, j",
      "    do j = 2, ny-1",
      "        do i = 2, nx-1",
      "            u(i, j) = ((u(i+1, j) + u(i-1, j)) * dy2 + &",
      "                       (u(i, j+1) + u(i, j-1)) * dx2) / (2*(dx2+dy2))",
      "        end do",
      "    end do",
      "    end subroutine",
      "",
      "    subroutine for_update2(u, dx2, dy2, nx, ny)",
      "    real(8), intent(inout) :: u(nx,ny)",
      "    real(8), intent(in) :: dx2, dy2",
      "    integer, intent(in) :: nx, ny",
      "    u(2:nx-1,2:ny-1) = ((u(3:,2:ny-1)+u(:ny-2,2:ny-1))*dy2 + &",
      "            (u(2:nx-1,3:) + u(2:nx-1,:ny-2))*dx2) / (2*(dx2+dy2))",
      "    end subroutine",
      "",
      "`f2py` generally just works (on linux and my desktop iMac it was fine), *BUT* on my mac laptop I had problems.",
      "",
      "Takeaway: if you have legacy Fortran or absolutely need that last factor of 2-5 in speed then spend the time to make it work."
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Numba",
      "",
      "[Numba](https://github.com/numba/numba) is aiming to be the next-generation tool for fast numerical Python.  It uses the",
      "remarkable [LLVM](http://llvm.org) compiler infrastructure to compile Python syntax to machine code."
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "<iframe src=\"http://www.slideshare.net/slideshow/embed_code/13809448\" width=\"427\" height=\"356\" frameborder=\"0\" marginwidth=\"0\" marginheight=\"0\" scrolling=\"no\" style=\"border:1px solid #CCC;border-width:1px 1px 0;margin-bottom:5px\" allowfullscreen webkitallowfullscreen mozallowfullscreen> </iframe> <div style=\"margin-bottom:5px\"> <strong> <a href=\"http://www.slideshare.net/teoliphant/numba\" title=\"Numba\" target=\"_blank\">Numba</a> </strong> from <strong><a href=\"http://www.slideshare.net/teoliphant\" target=\"_blank\">Travis Oliphant</a></strong> </div>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Add 3 lines of code (2 of which are boilerplate imports)",
      "# to convert pure Python to a compiled C-speed routine",
      "#",
      "from numba.decorators import jit",
      "from numba import float64",
      "",
      "@jit(argtypes=[float64[:,:], float64, float64])",
      "def numba_update(u, dx2, dy2):",
      "    \"\"\"Takes a time step using straight forward Python loops.\"\"\"",
      "    nx, ny = u.shape        ",
      "    dnr_inv = 0.5 / (dx2 + dy2)",
      "",
      "    for i in range(1, nx - 1):",
      "        for j in range(1, ny - 1):",
      "            u[i, j] = ((u[i - 1, j] + u[i + 1, j]) * dy2 +",
      "                       (u[i, j - 1] + u[i, j + 1]) * dx2) * dnr_inv"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit calc(func=numba_update, N=200)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 18.5 ms per loop"
       ]
      }
     ],
     "prompt_number": 36
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Numba II",
      "",
      "Consider an example that really cannot be vectorized: Welford's method for calculating variance.  This is an iterative (or streaming) method that is [numerically more stable][1] than the traditional $s^2 = \\frac{1}{N-1}\\sum (x_i - \\overline x)^2$.",
      "",
      "I googled on \"welford numpy\" and found [this implementation][2] which is shown in the next cell.",
      "",
      "[1]: http://www.johndcook.com/blog/2008/09/26/comparing-three-methods-of-computing-standard-deviation/",
      "[2]: http://projects.scipy.org/numpy/attachment/ticket/1098/test_var.py"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def welford(vals):",
      "    # welford's method (Welford 1962, see also Knuth, AOCP Vol. 2 page 232)",
      "    m = vals[0]",
      "    s = 0.0",
      "    for k in xrange(1, len(vals)):",
      "        oldm = m",
      "        x = vals[k]",
      "        newm = oldm + (x - oldm) / (k + 1)",
      "        s = s + (x - newm) * (x - oldm)",
      "        m = newm",
      "    return s / (k+1)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 40
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vals = np.random.uniform(size=1e6)",
      "%time var = welford(vals)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 5.44 s, sys: 0.07 s, total: 5.52 s",
        "Wall time: 5.49 s"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "markdown",
     "source": [
      "**OK** now let's just go for it and numba-fy the Python `welford` function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "numba_welford = jit(argtypes=[float64[:]], restype=float64)(welford)",
      "%timeit numba_welford(vals)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 12.4 ms per loop"
       ]
      }
     ],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit np.var(vals)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 18.7 ms per loop"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "markdown",
     "source": [
      "Summary",
      "--------",
      "",
      "Python provides a strong suite of options for doing computationally intensive numerical analysis.  In combination with batteries-included standard library and the full scientific stack (NumPy, SciPy, Matplotlib), this makes Python a great choice for many analysis tasks."
     ]
    }
   ]
  }
 ]
}