14_Parallel.ipynb

{
 "metadata": {
  "celltoolbar": "Slideshow",
  "name": "",
  "signature": "sha256:37d9700aed0ef121cd6b2da5f2e33a7e8cd30fdeb970ba067691bb6edcce5cec"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Parallel programming\n",
      "===================="
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Because my boss tells me to\n",
      "---------------------------\n",
      "\n",
      "Much like optimisation, parallel processing is one of those things that can easily make everything you do very complicated and incomprehensible to anyone else giving you little in return (or worse).\n",
      "\n",
      "Parallelising will not reduce the amount of work your CPU has to do. There is no magic that makes your code more efficient. For a problem that is entirely CPU bound, the best case scenraio is that you will get your results quicker by parallelising it, but if you code things badly it will take longer and use more resources. For codes where the CPU is idle a lot, during IO for example, parallel codes can take advantage of idle CPU time, but many scientific codes will have their major component that being CPU intensive.\n",
      "\n",
      "Nevertheless, if we can take advantage of the availabilty of large computing resources then we can possibly achieve results that would take years in days or hours."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Speedup\n",
      "-------\n",
      "\n",
      "The normal measures for how well something is parallelised is either **speedup** or **efficiency**: \n",
      "\n",
      "$S_p = \\frac{T_1}{T_p}$; \n",
      "\n",
      "$E_p = \\frac{S_p}{p} = \\frac{T_1}{pT_p}$\n",
      "\n",
      "Where $T_1$ is the time to execute on a single processor, and $T_p$ is the time to run on $p$ processors. These refer to **wall clock time** or **execution time**, the **CPU time** will be the total time doing computation over all CPUs added together."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Not everything parallelises well. **Amdahl's Law** is a limit to what you can do with parallelisation:\n",
      "\n",
      "$S_{max} = (F_{serial} + \\frac{F_{parallel}}{p})^{-1}$\n",
      "\n",
      "* *If you only parallelise half of your CPU bound code you can theoretically only get a maximum speedup of 2*.\n",
      "    + If my code does two things I need to parallelise both.\n",
      "* Profile before doing any optimisation to see if it is worth it.\n",
      "    + Make single processor performace the best it can be first."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "def amdahls_law(serial, p):\n",
      "    return (serial + (1-serial)/p)**-1\n",
      "\n",
      "x = np.arange(1, 64)\n",
      "plt.plot(x, [amdahls_law(0.5, item) for item in x], label=\"50% parallel\")\n",
      "plt.plot(x, [amdahls_law(0.2, item) for item in x], label=\"80% parallel\")\n",
      "plt.plot(x, [amdahls_law(0.05, item) for item in x], label=\"95% parallel\")\n",
      "plt.legend(); plt.xlabel(\"number of processes\"); plt.ylabel(\"speedup\")"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 1,
       "text": [
        "<matplotlib.text.Text at 0x40d4d50>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEPCAYAAABLIROyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVOX+wPHPsAmBiAsgiguhoigIikuaigu44i43W1Sw\n/JVXzbq3utbtpplbardcstUly63yqrmllbiUSyJuuIYQ7iKC7NvM+f1xYmRTQJkZBr7v1+t5nWXO\nzPmC8jxznvOc56tRFEVBCCFEtWNh6gCEEEKYhjQAQghRTUkDIIQQ1ZQ0AEIIUU1JAyCEENWUNABC\nCFFNGawBCA8Px9XVFR8fn0L7Fy9eTKtWrWjTpg1vvPGGoU4vhBCiFFaG+uCwsDAmT57MmDFj9Pv2\n7NnDli1bOHnyJNbW1iQkJBjq9EIIIUphsCuAbt26Ubt27UL7li1bxrRp07C2tgbA2dnZUKcXQghR\nCqPeA7h48SL79u2jc+fOBAYGcvToUWOeXgghRAEG6wIqSV5eHklJSRw6dIjff/+d0NBQLl26ZMwQ\nhBBC/MWoDYC7uzvDhw8HoEOHDlhYWJCYmEjdunULHdesWTNiYmKMGZoQQpg9T09P/vjjjzIfb9Qu\noKFDh/LLL78AcOHCBXJycopV/gAxMTEoimK25Z133jF5DBK/6eOojvGbc+xVIf7yfnE22BXA6NGj\n2bt3L4mJiTRq1Ih3332X8PBwwsPD8fHxwcbGhq+++spQpxdCCFEKgzUAa9euLXH/6tWrDXVKIYQQ\n5SBPAhtAYGCgqUN4JBK/aZlz/OYcO5h//OWlURSl0iWE0Wg0VMKwhBCiUitv3WnUUUBCCOOqU6cO\nSUlJpg5DVLDatWtz586dR/4cuQIQogqTv6Wq6X7/ruX995Z7AEIIUU1JAyCEENWUNABCCFFNSQMg\nhBAGFhgYyJdffgnAypUr6datW5neN27cON5++22DxSUNgBDCZAIDA7Gzs6NmzZrUrFmTVq1aFXr9\n559/pmXLltjb29OrVy/i4+P1r61Zs4YGDRrg4eFBRESEfn9MTAxdu3atVDe/NRoNGo3GaO8rK2kA\nhBAmo9FoWLp0KampqaSmpnL27Fn9a7dv32bEiBHMmjWLpKQkAgIC+Nvf/gaoMwtPmzaNqKgolixZ\nwuTJk/XvmzJlCh9++KFBK86C8vLyDPr5hmzIpAEQQpjU/Sq4jRs30qZNG0aMGIGNjQ3Tp0/nxIkT\nXLhwgcTERBo2bIirqyu9e/fWTyv/3Xff0ahRIzp06PDAc65cuZKuXbsyefJknJycaNWqlX6iSoAV\nK1bg7e2No6Mjnp6efPbZZ/rXIiIicHd35/3338fNzY3x48eTnJzMoEGDcHFxoU6dOoSEhHD16tUy\n/fznzp0jKCiIunXr0rJlS7799tsyva8iSAMghDCpadOm4ezszJNPPsnevXv1+6Ojo2nbtq1++7HH\nHqNZs2ZER0fj4uJCYmIiV69eZffu3bRp04a0tDRmzZrFnDlzynTeI0eO0KxZMxITE5kxYwbDhw/X\nPzTn6urKtm3bSElJYcWKFbzyyitERUXp33vz5k2SkpKIj4/n008/RafTMX78eOLj44mPj8fOzo5J\nkyaVGkN6ejpBQUE8++yzJCQksG7dOiZOnFjoSsiQpAEQoprTaB69PKx58+YRGxvLtWvXmDBhAiEh\nIcTGxgJq5ejo6FjoeEdHR1JTU9FoNCxbtoyRI0fywQcf8Pnnn/Of//yHKVOmcPz4cXr16kW/fv2I\njo6+77ldXFx4+eWXsbS0JDQ0FC8vL7Zt2wbAgAED8PDwAKB79+4EBwezf/9+/XstLCyYMWMG1tbW\n2NraUqdOHYYNG4atrS0ODg68+eabhRqz+9m6dSseHh6MHTsWCwsL/Pz8GD58uNGuAmQqCCGqOVPe\nK+3YsaN+fcyYMaxdu5Zt27YxadIkHBwcSElJKXT83bt3qVmzJgC9evXi4MGDAJw4cYJjx46xYMEC\nmjZtyq+//kp8fDzPP/+8/piiGjZsWGi7SZMmXL9+HYAdO3YwY8YMLl68iE6nIyMjA19fX/2xzs7O\n2NjY6LczMjJ45ZVX+PHHH/VXEWlpaSiK8sB7EX/++SeHDx8ulD89Ly+PMWPG3P+XVoHkCkAIUSm1\nbt2aEydO6LfT09OJiYmhdevWhY5TFIXJkyezaNEiEhIS0Gq1NGrUiICAAE6ePHnfzy/aR//nn3/S\noEEDsrOzGTFiBK+//jq3bt0iKSmJAQMGFLpXUbRSX7hwIRcuXODIkSPcvXuXvXv36pO0PEjjxo3p\n0aMHSUlJ+pKamsrSpUtL/f1UBGkAhBAmcffuXX788UeysrLIy8vjm2++Yf/+/fTr1w+AYcOGcfr0\naTZu3EhWVhYzZszAz8+PFi1aFPqcL774gvbt2+Pr60vdunXJzMzk7Nmz7NmzB09Pz/ue/9atWyxa\ntIjc3Fy+/fZbzp07x4ABA8jJySEnJ4d69ephYWHBjh072LVr1wN/lrS0NOzs7KhVqxZ37txhxowZ\nZfodDBw4kAsXLvD111+Tm5tLbm4uv//+O+fOnQMMOwIIpAtICGEiubm5vP3225w7dw5LS0tatWrF\n5s2badasGQD16tXj+++/Z9KkSTz77LN07tyZdevWFfqM27dvs2jRIn03j5WVFUuWLKFXr17Y2dmx\nYsWK+56/U6dOXLx4EWdnZ+rXr8/333+v74pZtGgRoaGhZGdnExISwpAhQwq9t+gVwNSpU3n66aep\nV68eDRs25NVXX2XLli0lnrfg2P6aNWuya9cuXn31VV599VV0Oh1+fn588MEHxY41BJkNVIgqTP6W\nSrZy5Uq+/PLLQjd2zYnMBiqEEOKRGKwBCA8Px9XVFR8fn2KvLVy4EAsLiwpJaCCEEOVl6K4Vc2Gw\nBiAsLIydO3cW23/58mV2795NkyZNDHVqIYR4oLFjx7Jv3z5Th2FyBmsAunXrVmhsa75XX32V999/\n31CnFUIIUUZGvQewefNm3N3dCz1QIYQQwjSMNgw0IyOD2bNns3v3bv0+GZ0ghAHk5MCvv8KOHaaO\nRFRyRmsAYmJiiIuL00/udOXKFdq3b8+RI0dwcXEpdvz06dP164GBgQQGBhopUiHMUHy8WuHv2AF7\n9oCXF/z1QJWouiIiIgrlQigvgz4HEBcXR0hICKdOnSr2moeHB5GRkdSpU6d4UDJ2WYgHy8uDgwdh\n2za13LgBfftC//4QHAzOzoD8LVVVlf45gNGjR9OlSxcuXLhAo0aNij2RJ0OwhCin5GRYtw6efhpc\nXeHll8HaGj7/XG0Avv4annlGX/mLyqPapYRcu3Yt165dIzs7m8uXLxMWFlbo9UuXLpX47V8IUUBs\nLHz4IfTqBY0bwzffQGAgnDwJx47BzJnQuTNYWpo60ody5coVQkJCqFu3Lm5ubkyePBmtVqt/XVJC\nSkpIIaoPRYHjx2H6dPDzg06d4NQp9dv+9evwww8wYQIUmcrYXE2ZMoV69epx/fp1jh8/zt69e/n4\n448BSQmZT1JCClGV6XTqqJ1//AMefxxGjIDUVFiyRK30v/wShgwBe3tTR1rhoqOj+dvf/oaNjQ2u\nrq7069ePM2fOAJIS0hikARDCFLRadbTO3/8O7u7w4ovg4ACbN8Mff8DChfDkk2bbtVNWffv2Zc2a\nNWRmZnL16lV27Nihnw5aUkIankwHLYSxaLWwfz9s2AAbN6rdOCNHQkQEFJnj3pg0Mx69q0R55+G6\nKaZPn06fPn1wdHREq9Uybtw4/dTL6enpOBe5oV1SSkhbW9tiKSFnzpyJjY0NCxcuLJZAJl9+SkiA\n0NBQFi5cyLZt23j22WcZMGCA/riCKSH9/f2Bwikh89NCDhs2TP+eN998k169epX68xdMCQkUSgn5\nn//8pxy/yYcjDYAQhqTTqcM1162D774DNzcIDVW7fB6QrMSYHrbyfuTzKgp9+/Zl1KhRHD58mNTU\nVMLDw3njjTeYN2+epIQ0AukCEqKi5d/IfeMN8PBQb9q6usK+ferInX/9q9JU/qZ0+/ZtIiMjmTRp\nEtbW1tSpU4dx48axfft2QFJCGoM0AEJUlLg4mD0bWreGoUPBwkIdtXP6NPz739C8uakjrFTq1auH\nm5sby5YtQ6vVkpyczKpVq/T9/kOHDpWUkAYeyioNgBCPIjlZfRCre3cICIDLl9Xt2FiYMwd8fUEe\neiyRRqNh48aN/PDDD9SrV4/mzZtTo0YN/vvf/wJqN8v333/PW2+9RZ06dTh69Oh9U0LOnDkTKJwS\ncuLEiSxevPi+5y+YEvLtt9/Wp4SsWbOmPiVknTp1WLt2bZlSQmZmZlKvXj26dOlC//7979v1U1JK\nyHXr1tGwYUPc3NyYNm0aOTk5xY41BEkJKUR5abWwezesWgXbt0OfPjBmjDr3To0apo6uEPlbKpmk\nhFTJTWAhyurCBVi+HFavVkfwjB2rjtWvW9fUkQnxUKQBEOJB0tPV0TtffAEXL8Jzz6nf/r29TR2Z\neASSElIlXUBClCQqCj79VB2z36ULjB8Pgwapk6+ZEflbqpqkC0iIipaeDmvXqhX/rVvwwgvqPDxV\nZN4dIYqSKwAhzpyBjz+GNWvU0TwTJqhz61eBaRjkb6lqkisAIR5Fbq46787SpXD+vPpt/+RJdV4e\nIaoJaQBE9ZKQAJ99BsuWqTNv/v3vMGwYFHisX4jqQhoAUT2cPAkffaROwjZihDp+v8DcLkJUR/Ik\nsKi6dDo1X26vXmqu3McfV8fyf/GFVP7CqKpdSkghTCYrS02i0qYNvP22OoQzNhbeekvy5VYyZ8+e\npVevXjg5OdG8eXM2bdqkfy0uLg4LCwtq1qypL7NmzdK/LikhH51Bu4DCw8PZtm0bLi4unDp1CoDX\nXnuNrVu3YmNjg6enJytWrKBWrVqGDENUF0lJ6mieJUugXTt12bOnzMVTSeXl5TFkyBAmTpzIzz//\nTEREBCEhIURFRdG8wMR5KSkpxSrBgikhjx49yuTJk/V1jClSQlpZGa4qNduUkGFhYezcubPQvuDg\nYKKjozlx4gQtWrQoc/YeIe7r2jV47TVo1kx9Wvenn+51/UjlX2mdO3eO69evM3XqVDQaDT179qRr\n166sXr260HE6na7YeyUlZMUwaAPQrVu3QokOAIKCgrCwUE/bqVMnrly5YsgQRFX2xx/qmP02bSAn\nR316d+VKdTpmYZZ0Oh2nT58utK9JkyY0atSI8PBwEhMTAXWmUEkJ+ehMeg9g+fLlhVKvCVEm58+r\ns2927gz166s3dj/6CBo3NnVk5kmjefTyELy8vHBxcWH+/Pnk5uaya9cu9u3bR2ZmJqBW8kePHiU+\nPp7IyEhSU1N55plnADUlY35KyA8++KBYSshevXrRr18/oqOj73v+/JSQlpaWhIaG4uXlxbZt2wAY\nMGAAHh4eQOGUkPkKpoS0tbWlTp06DBs2DFtbWxwcHHjzzTfZu3dvqb+DgikhLSwsCqWENAaTDQOd\nNWsWNjY2PP300yW+Pn36dP16YGAggYGBxglMVF7R0fDee/Dzz/Dyy7B4Mcj9o0dnopul1tbWbNq0\nicmTJzNv3jw6dOhAaGgotra2ANjb29OuXTtArayXLFmCm5sb6enp2NvbS0pI1O6ogjfAy8skDcDK\nlSvZvn07P//8832PKdgAiGru7FmYPl1Nnv7qq+qDXH/lhRXmzcfHp1AF1qVLF8LCwh74nqL3BPJT\nQi5ZsqRQSkgXF5dyp4QcMmSIPiXk119/zZAhQ7C0tGTYsGFlTgnp4uLC8ePHadeuXakNQH5KyNIy\njt1P0S/HZc1Els/oXUA7d+5k/vz5bN68Wd/SC1GiP/5Qu3p69FBH9cTEqHl2pfKvMk6dOkVWVhYZ\nGRksWLCAmzdvMm7cOEDtoz9//jw6nY7ExESmTJlCz5499Unh80lKyIdn0AZg9OjRdOnShfPnz9Oo\nUSOWL1/O5MmTSUtLIygoCH9/fyZOnGjIEIQ5unwZnn9e7eNv1kxtCN54AxwcTB2ZqGCrV6+mQYMG\nuLq6smfPHnbv3o31X1NuX7p0if79++Po6IiPjw92dnasXbu20PslJeSjkdlAReWRmKjm0V2xAl58\nEf75TygyikyUj/wtlUxSQqrkSWBheunpMGsWeHlBRgacPq1uS+UvhEFJAyBMR6tVp2xo0UKdrO3g\nQfVJXjc3U0cmqjhJCamSLiBhGrt3q108jo6wcCF07GjqiKok+VuqmiQhjDBPZ8+qFf/58/D+++pc\n/PJNTAiTkC4gYRzJyfDKK2rKxT591DSMw4dL5S+ECUkDIAxLq1Xn32/ZUr3Be+aM2hBIBi4hTE66\ngIThHDoEkyZBjRpqBq6/HusXxlO7dm252VkFFZ1k82HJTWBR8e7cgX/9C7ZuVfv5n3lGunqEMAJ5\nDkCYjqLAqlXg7a128Zw5A88+K5W/EJWUdAGJinH2rPr0bnq6+s0/IMDUEQkhSiFXAOLR5OTAzJnQ\nrRuMHAmHD0vlL4SZkCsA8fCOHFEnbWvUCI4dk4QsQpgZuQIQ5ZeRoc7LP3jwvZu9UvkLYXakARDl\n89tv4OcHN26ok7Y9/bTc5BXCTEkXkCibrCz4z39g9WpYulR9ilcIYdakARClO3oUxo6FVq3UWTud\nnU0dkRCiAkgXkLg/rVadl3/gQPj3v+Hbb6XyF6IKkSsAUbI//1Qf4rK2hshIcHc3dURCiAomVwCi\nuLVroUMHdZTPTz9J5S9EFWWwBiA8PBxXV1d8fHz0++7cuUNQUBAtWrQgODiY5ORkQ51ePIy0NBgz\nBmbMgJ074bXXwEK+IwhRVRnsrzssLIydO3cW2jd37lyCgoK4cOECvXv3Zu7cuYY6vSivkyfVJ3it\nrNQuH5m5U4gqz2ANQLdu3YpNWbplyxbGjh0LwNixY9m0aZOhTi/KSlHg88+hd2946y1Yvhzs7U0d\nlRDCCIx6E/jmzZu4uroC4Orqys2bN415elFUair83//BqVOwf7+atEUIUW2YbBSQRqN5YKKK6dOn\n69cDAwMJDAw0fFDVyblzaj7erl3VCdwee8zUEQkhyikiIoKIiIiHfr9BE8LExcUREhLCqVOnAGjZ\nsiURERHUr1+f69ev07NnT86dO1c8KEkIY1j/+x9MmABz5qiTuQkhqoRKnRBm8ODBrFq1CoBVq1Yx\ndOhQY55eaLXw5pswdaqaolEqfyGqNYNdAYwePZq9e/dy+/ZtXF1deffddxkyZAihoaHEx8fTtGlT\nNmzYgJOTU/Gg5Aqg4t25A6NHQ14erFsnT/QKUQWVt+6UnMDVwblzEBICQ4bA3LnqUE8hRJVTqbuA\nhAns2gU9eqhDPBcskMpfCKEntUFVpSiweLF6o/f77+HJJ00dkRCikpEGoCrKzYXJk+HXX+HgQWja\n1NQRCSEqIWkAqprUVAgNVbN0/fYb1Kxp6oiEEJVUme4B5OTkcOLECU6dOkVOTo6hYxIP6/p1tb+/\nUSPYskUqfyHEA5XaAGzbtg1PT0+mTJnCpEmT8PT0ZPv27caITZTH2bPQpQuMGAGffio3e4UQpSp1\nGKiXlxfbtm2jWbNmAMTExDBgwADOnz9vuKBkGGj57NsHo0apo3yee87U0QghTKS8dWepXxMdHR31\nlT/A448/jqOj48NFJyreli3qE71r1kCfPqaORghhRkq9AnjxxReJj48nNDQUgG+//ZbGjRsTFBQE\nwPDhwys+KLkCKJvVq9WkLVu3qnP5CyGqtQp/EnjcuHH6DwZQFKXQLJ4rVqx4iDBLCUoagNItXgzv\nvw8//gje3qaORghRCchUEFWdosDMmeq3/927ZYy/EEKvwu8BhIWFFTsBwPLly8sZmnhkigL//Kea\nqH3/fqhf39QRCSHMWKkNwMCBA/WVfmZmJv/73/9o0KCBwQMTRSiKOo3zb79BRAQUSbcphBDlVe4u\nIJ1OR9euXTl48KChYpIuoKIURZ3a4ehR2LkTSphCWwghKrwLqKgLFy6QkJBQ3reJh6XTwd//DseP\nqzd8a9UydURCiCqi1AbAwcFB3wWk0WhwdXVl3rx5Bg9MoFb+L74IZ86olb88fyGEqEAyCqiyUhSY\nOBGio2HbNpnXRwhRqgrrAoqMjCw03r+odu3alS8yUXaKoj7gdeyYOuJHKn8hhAHc9wogMDAQjUZD\nZmYmkZGR+Pr6AnDy5EkCAgLkJrAhzZgBGzfCnj1Qp46poxFCmIkKSwkZERHBnj17aNCgAceOHSMy\nMpLIyEiioqIeeRjonDlzaN26NT4+Pjz99NNkZ2c/0udVKQsWwNq1aipHqfyFEAZU6nTQ586dw8fH\nR7/dpk0bzp49+9AnjIuL4/PPP+fYsWOcOnUKrVbLunXrHvrzqpRly2DpUrXbx9XV1NEIIaq4UkcB\n+fr68vzzz/Pss8+iKApr1qyhbdu2D31CR0dHrK2tycjIwNLSkoyMDBo2bPjQn1dlrFsHs2bB3r3g\n7m7qaIQQ1UCpo4AyMzNZtmwZ+/fvB6B79+689NJL2NraPvRJP/vsM/7xj39gZ2dH3759Wb16deGg\nqts9gF9+gaeegp9/hgJXW0IIUR4GmQwuIyOD+Ph4WrZs+UjBgZpQJiQkhP3791OrVi1GjRrFyJEj\neeaZZ+4FpdHwzjvv6LcDAwMJDAx85HNXSsePQ3AwbNgAVfVnFEIYREREBBEREfrtGTNmVGwDsGXL\nFl577TWys7OJi4sjKiqKd955hy1btjxUwOvXr2f37t188cUXAKxevZpDhw6xdOnSe0FVlyuAuDh4\n8kn44AM1kbsQQjyCChsFlG/69OkcPnyY2n9NPubv78+lS5ceOsCWLVty6NAhMjMzURSFn376Ce/q\nOJ99YiL06wevvy6VvxDCJEptAKytrXEqMvmYhUWpb7uvtm3bMmbMGAICAvTPFkyYMOGhP88sZWZC\nSAgMGQJTppg6GiFENVVqF1B4eDi9e/dm7ty5bNy4kUWLFpGbm8snn3xiuKCqcheQosAzz6jz/KxZ\nA4/QmAohREEV3gW0ePFioqOjqVGjBqNHj8bR0ZEPP/zwkYKs1mbOhJgYWLFCKn8hhEmVeTK49PR0\n7O3tDR0PUIWvADZsUDN6HT4Mbm6mjkYIUcVU+BXAb7/9hre3t34I6IkTJ5g4ceLDR1hd/f67Oq//\n5s1S+QshKoVSG4CpU6eyc+dO6tWrB6g3cffu3WvwwKqUK1dg2DD4/HPw9zd1NEIIAZShAQBo3Lhx\noW0rq3InEqu+MjNh6FCYNEldCiFEJVFqTd64cWN+/fVXAHJycli0aBGtWrUyeGBVgqKo3T6envDG\nG6aORgghCim1AVi2bBkvv/wyV69epWHDhgQHBxd6alc8wGefqTd8Dx+GByTXEUIIU5CUkIZy+LD6\nsNeBA9CihamjEUJUAxU+Cih/8rZ69erh7OzMkCFDHmkqiGrh1i0YNUq96SuVvxCikiq1AXj66acJ\nDQ3l+vXrXLt2jVGjRjF69GhjxGae8vLUqZ2fe06d6kEIISqpUruAfH19OXnyZKF9bdu25cSJE4YL\nypy7gKZNg8hI2LEDLC1NHY0Qohopb91Z6k3g/v37M2fOHP23/vXr19O/f3/u3LkDQB3JW3vP7t2w\nejUcOyaVvxCi0iv1CqBp06Zo7jOCRaPRGOR+gFleAdy8qT7k9fXX0KuXqaMRQlRDFX4T+P333+fE\niRPExsYSFhaGn58f3333HbGxsXIzOJ9OB2PGQHi4VP5CCLNRagMwc+ZMHB0dOXDgAL/88gvjx4+X\nuYCKWrgQ0tJg+nRTRyKEEGVWagNg+Vdf9tatW3nhhRcYNGgQOTk5Bg/MbBw+DPPnq3P7yxQZQggz\nUmoD0LBhQyZMmMD69esZOHAgWVlZ6HQ6Y8RW+d29C6NHwyefQJMmpo5GCCHKpdSbwOnp6ezcuRNf\nX1+aN2/O9evXOXXqFMHBwYYLylxuAoeFQY0aagMghBAmVt66U6aCeFhbtsArr8CJE+DgYOpohBCi\n4kcBGUJycjIjR46kVatWeHt7c+jQIVOE8fBu34YXX4SVK6XyF0KYLZNcAYwdO5YePXoQHh5OXl4e\n6enp1KpV615QlfkKQFEgNBSaNlVv/gohRCVR6buA7t69i7+//wOfIajUDcDatfDee+p0D7a2po5G\nCCH0Kn0XUGxsLM7OzoSFhdGuXTteeOEFMjIyjB3Gw7l2DaZOha++kspfCGH2jD5wPS8vj2PHjrFk\nyRI6dOjA1KlTmTt3Lu+++26h46YXeKgqMDCQwMBA4wZalKLA88/DxInQvr1pYxFCCCAiIoKIiIiH\nfr/Ru4Bu3LjBE088QWxsLAAHDhxg7ty5bN269V5QlbEL6Ouv4YMP1Ae/rK1NHY0QQhRT6buA6tev\nT6NGjbhw4QIAP/30E61btzZ2GOVz+zb8859qghep/IUQVYRJRgGdOHGC559/npycHDw9PVmxYkXl\nHgU0dizUq6fO+SOEEJVUpR8FVBaVqgH46Se17//0aRnzL4So1Cp9F5BZychQH/j6+GOp/IUQVY5c\nATzItGkQF6eO/RdCiEpOuoAqysmT0KcPnDoFrq6mjUUIIcpAuoAqgk4HEybA7NlS+QshqixpAEqy\nahVYWKgpHoUQooqSLqCi7t6Fli3hhx8gIMA0MQghxEOQewCP6h//gJQU9aEvIYQwI9IAPIqzZ6F7\nd4iOBhcX459fCCEegdwEfliKAlOmwL//LZW/EKJakAYg3+bN6nTPEyeaOhIhhDAK6QICyMyE1q3h\ns8/Usf9CCGGGpAvoYSxcCP7+UvkLIaoVuQK4dQtatYLff4fHHzfOOYUQJqVTdOTp8vRFq9PeW1e0\nxfbn73vQdr9m/bDQmPY7dXnrTqNnBKt0Zs6EZ5+Vyl+IAhRFQatoydHmkKPNIVebqy51ucW2c7W5\nhfYX3FdwmafLK7QvT5dXbH/+eqHlX/sL7nvUAmBlYYW1hTWWFpZYWVhhZWGFpabAuoVlse2CxxTc\ntrSwJNgz2OQNQHlV7yuAmBjo1Ekd/unsbPjzCVGEoijk6fLIysvSl2xttrrMyy62na3NLnGZo80p\ntC9/O0f8gamFAAAZo0lEQVSbU+j1/Aq9LMVSY4m1pTU2ljZYW6hLG0sbrC2tsbawLvRa0fVCy7/W\n8yvc/P1WFlbF1gsek79e8LWC+wpWzAX3Fd0uuC+/0ja3irqs5DmA8njqKWjTRh36KcRf8nR5pOek\nk56bTkZuRoklMzdTXeZl6rcz8zLvLQusZ+VlkZmbqa/g8/flFwA7KztsrWypYVUDOys7aljVULct\na9x//a/tGpY1sLG0KbZuY2lTbLvgvoIVesF9+ZWypYWlif8lRHlJA1BWR4/C4MFw8SLY2xv2XMIg\nFEUhW5tNanYqqTmppGSnkJqdSlpOGqk5fy3/2taX3Hvr+ZV8wfX0nHS0ihZ7a3ses34Mext1mV/s\nrOywt7HHzsqu0D47a7sSl7ZWtoX22VrZFitWFtITKyqGNABloSjqiJ/QUPi//zPcecQD5WpzSc5K\nJikrieSsZH1JylS372bf5W7WXXWZfZeU7BTuZqnL/ALgWMMRxxqO1KxRk5o2NfVLBxsH/bJgsbex\nx97aXr/uYOOAvbW9fr+NpQ0ajcbEvx0hyk9uApfFrl1w5QqMH2/qSKoERVG4m32X2xm3i5XEjETu\nZN4hMVNd5pekrCQyczNxsnWitl1tatvWxsnWSV9q1ahFLdtauDm4Ucu2ln67Vo1ahSp8WytbU//4\nQpgtkzUAWq2WgIAA3N3d+eGHH4x3Yp0O3ngD5swBq+rZ/pWFTtGRkJ7A9bTr3Ei7wY20G9xMu8nN\ndLXcSr/FrfRbJKQncDvjNrZWtjjbO1PvsXr3il096tjVoYlTE+ra1aWOXR19qW1Xm5o2NeWbthAm\nZLIa8KOPPsLb25vU1FTjnnjtWrCzg2HDjHveSiQzN5MrKVeIvxvPlZQrXE29ytWUq+oy9SrXUq+R\nkJ6g/wZe36E+9R3q42rvSoOaDfCv74+LvQsu9i76Sl++iQthfkzSAFy5coXt27fz1ltv8cEHHxjv\nxFotvPuumuS9Cn/zzMjNIDYpltjkWOKS44hNiiXubhxxyXFcvnuZlOwUGjo2pHGtxrg7utOwZkNa\n1mtJ78d707BmQxrUbICrgys2ljam/lGEEAZkkgbglVdeYf78+aSkpBj3xOvWqTN99upl3PMaQHZe\nNn/c+YPziec5d/scf9z5g5ikGP648weJGYk0dWqKR20PPJw8aOrUlM7unWni1ITGtRrjYu9SZcdB\nCyHKzugNwNatW3FxccHf35+IiIj7Hjd9+nT9emBgIIGBgY92Yq1Wfep3yRKz+vaflZfFudvnOH3r\nNNG3ojmdcJqzCWe5knKFJk5N8KrrhVddLzq7d+ZZ32dpVqcZDWs2lDHcQlQDERERD6xHS2P0YaBv\nvvkmq1evxsrKiqysLFJSUhgxYgRfffXVvaAMMQx0zRq162f//krbACRmJBJ1I4rjN44TdSOKqOtR\nxCbH4lnbkzYubWjj0obWzq3xdvbm8dqPY21pbeqQhRCViFk9B7B3714WLFhQbBRQhTcAWq063fPi\nxRAUVHGf+wiy8rKIuh7F4auHOXL1CIevHiYhPQG/+n741/fH380fv/p+eDt7S1+8EKJMzO45AKMM\nA9ywAerWNel0z6nZqfx2+Tf2/bmPffH7OHb9GF51vejUsBNBjwfx7+7/pmW9ltI3L4Qwmqr/JLBW\nq873s2iRUb/95+nyOHzlMD/G/MiumF2cvnWa9g3a06NJD7o36U5n98442DgYLR4hRNVndlcABvft\nt1C7tlG+/SdmJLLl/Ba2XtzKL7G/0NSpKX09+zKn9xyeaPSEjJUXQlQqVfsKQKsFHx/48EMIDn70\nzyvB1ZSrbDq3iY3nNnL02lGCHg9isNdggj2Dqe9Q3yDnFEKIksgVQEGbNoGjY4V3/aTlpPHdme9Y\ndWIVJ26cYFCLQUzuOJlgz2Aes36sQs8lhBCGUnUbAEWB+fPh9dcrZNinoijs/XMvK4+vZNO5TXRv\n0p1JHSYxqMUgaljVqICAhRDCuKpuA3DgACQmwpAhj/QxWXlZrDm1hg8PfYhW0TLefzzz+szD1cG1\nggIVQgjTqLoNwPz58OqrYPlwT8TeSr/Fx79/zCdHP6GdWzsWBi+kz+N9ZPZKIUSVUTVvAp87Bz16\nQFycOvNnOSRlJvH+r+/z2bHPGNlqJFM7T6WVc6uHj0UIIYxEbgIDLFwIEyeWq/LPyM1g0eFFLDy4\nkGEth3HixRO4O7obMEghhDCtqtcA3LgB338PFy6U6XBFUVhxfAVv73mbro26ciDsAF71vAwcpBBC\nmF7VawAWL4bRo6FevVIPjbkTw4StE0jJTmHT3zbRoWEHIwQohBCVQ9W6B5CWBk2bwuHD4Ol538O0\nOi0fHvqQOQfm8K8n/8XUzlOxsqh6baEQonqp3vcAvvwSAgMfWPmfu32OMf8bg72NPYeeP0SzOs2M\nF58QQlQiVecKQKtVK/7166FTpxIP+eH8D4RvCWdG4AxeDHhRZt4UQlQp1fcKYOtWcHMrsfJXFIXZ\n+2fz8dGP2fLUFp5o9IQJAhRCiMql6jQAH3+sDv0sIi0njbDNYcTfjefI80do6NjQBMEJIUTlUzX6\nQC5cgOPHYdSoQruvp16n6/KuONg4sHfcXqn8hRCigKrRACxbBuHhYHtvvv2E9AR6f9WbUd6jWD54\nuczFL4QQRZj/TeD0dGjcGCIj1SGgqNM59FzVk5AWIczsNdNwgQohRCVS3pvA5n8FsHYtdO2qr/xT\nslPo900/env05t2e75o2NiGEqMRM0gBcvnyZnj170rp1a9q0acOiRYse7oMUBZYuhb//HVDn8xm0\nZhDt6rdjQfACmblTCCEewCRdQDdu3ODGjRv4+fmRlpZG+/bt2bRpE61aqbNulvky5rffYOxYOH+e\nPHQMXDMQNwc3lg9ZLmP8hRDVjll0AdWvXx8/Pz8AHBwcaNWqFdeuXSv/B338Mbz0ElhYMGvfLHSK\nji8HfymVvxBClIHJbwLHxcXRo0cPoqOjcXBwUIMqSyt26xZ4ecGlSxxIjWbkhpFE/V8UbjXdjBC1\nEOYt/89LUe6VgtslHVPW1wsuH/T+ko4vy2eVtq8sn1naZ5W0XdrxrVpVSPbZR2JWTwKnpaUxcuRI\nPvroI33ln2/69On69cDAQAIDAwu/+YsvYMQIku00PLvqWT4P+bzaV/6Kos6IkZurlrw8tRRcL1i0\n2uLrJS1LKzpd8fWCy/xSdLssrylKyeslHZP/h150X1mXRddL2i7ra+Up+f92ZT3ufseXtL/gvoLv\nL0ijuVcKbuevl7SvtNcLLh/0/pKOL8tnlbavLJ9Z2meVtP2g40+cAGtrjCoiIoKIiIiHfr/JrgBy\nc3MZNGgQ/fv3Z+rUqYWDKq0V0+nA0xNlwwae+nMBLo+5sHjAYgNHXDpFgcxMdVLStDR1hGp6urqe\nmQkZGWrJX8/KulcyM9VldnbxkpOjloLrubnFl7m5YGGh/ifML1ZWarG2VrNj5u+ztCy8zF8vaTt/\nn4VF4X35JX9/wdeLrhcsRfdpNMWP1WiKH5N/XEnHFKzEir5e0vJB733Qdmn7H7ZA+Y673/El7S+4\nr+i6qFrM4gpAURTGjx+Pt7d3scq/TPbuhZo1WWl5ijMJZ1j5/MoKjU+rVXuYbt2ChIR75fZtSE6G\npKR7y7t3ISVFLampagXr4KAWe/t7y8ceU4ud3b2lnZ367JqTk7qsUePesmixsSlerK3vLfOLhdz+\nEEKUkUmuAA4cOED37t3x9fXVD9WcM2cO/fr1U4MqrRUbM4ZbLdxpXeNz9ozdQxuXNuU6f3a2mi44\nJgb++ENdxsfDtWtw9apa8Ts5gasrODurpV49tdSpA7Vrq6/nLx0doWZNtRj7ElAIIfKV9wrA5DeB\nS/LAHyIlBaVxY0JmeNHvieeY1HHSAz/rzh2IilLLsWNqiY2FRo3U2aObNVOXTZpAgwbQsCHUr69+\nsxZCCHNiFl1Aj2TDBm528CbOOoWXAl4q9rJWC7//Dtu2wfbtcPEi+PmBvz/07QvTpkHLlvJNXQgh\nzK4BUFasYJ7fTWYEvo+lhaV+/8mT8N//qmkB6teHAQPU7S5d1BuYQgghCjOvW4bnz5N9/gy/tnZk\nWKthgPptf8gQ9dt9y5Zw9CicOgXz5kH37lL5CyHE/ZhV9ahbsZx1fpa80+c9Io9a8PbbEB0Nb7wB\n69apo2qEEEKUjflcAWi1ZK34nB+7NUB3fgADB8Lw4eoonkmTpPIXQojyMpsrAO3OHfxhm4GXxwc8\n/7yGrVuhY0dTRyWEEObLbIaBxgZ35Cv7JD4+eIEd2zW0a2ei4IQQopKqksNAcxJuUGd/JKtabOGn\n3Rp8fEwdkRBCmD+zuAL4OvQpbCL34rP1On+lDBBCCFFElbwCaDD4KaxHPC2VvxBCVCCzuAIQQghR\nOrPICCaEEML0pAEQQohqShoAIYSopqQBEEKIakoaACGEqKakARBCiGpKGgAhhKimTNIA7Ny5k5Yt\nW9K8eXPmzZtnihCEEKLaM3oDoNVqmTRpEjt37uTMmTOsXbuWs2fPGjsMg4qIiDB1CI9E4jctc47f\nnGMH84+/vIzeABw5coRmzZrRtGlTrK2teeqpp9i8ebOxwzAoc/9PJPGbljnHb86xg/nHX15GbwCu\nXr1Ko0aN9Nvu7u5cvXrV2GEIIUS1Z/QGQKPRGPuUQgghSqIY2cGDB5W+ffvqt2fPnq3MnTu30DGe\nnp4KIEWKFClSylE8PT3LVR8bfTbQvLw8vLy8+Pnnn2nQoAEdO3Zk7dq1tJK5noUQwqiMng/AysqK\nJUuW0LdvX7RaLePHj5fKXwghTKBS5gMQQghheJXuSWBze0gsPDwcV1dXfAokKr5z5w5BQUG0aNGC\n4OBgkpOTTRjh/V2+fJmePXvSunVr2rRpw6JFiwDziT8rK4tOnTrh5+eHt7c306ZNA8wn/nxarRZ/\nf39CQkIA84q/adOm+Pr64u/vT8eOHQHzij85OZmRI0fSqlUrvL29OXz4sFnEf/78efz9/fWlVq1a\nLFq0qNyxV6oGwBwfEgsLC2Pnzp2F9s2dO5egoCAuXLhA7969mTt3romiezBra2v++9//Eh0dzaFD\nh1i6dClnz541m/htbW3Zs2cPx48f5+TJk+zZs4cDBw6YTfz5PvroI7y9vfUj5Mwpfo1GQ0REBFFR\nURw5cgQwr/hffvllBgwYwNmzZzl58iQtW7Y0i/i9vLyIiooiKiqKyMhIHnvsMYYNG1b+2B95WE8F\n+u233wqNEJozZ44yZ84cE0ZUNrGxsUqbNm30215eXsqNGzcURVGU69evK15eXqYKrVyGDBmi7N69\n2yzjT09PVwICApTTp0+bVfyXL19Wevfurfzyyy/KoEGDFEUxr/8/TZs2VW7fvl1on7nEn5ycrHh4\neBTbby7x5/vxxx+VJ598UlGU8sdeqa4AqspDYjdv3sTV1RUAV1dXbt68aeKIShcXF0dUVBSdOnUy\nq/h1Oh1+fn64urrqu7PMKf5XXnmF+fPnY2Fx70/RnOLXaDT06dOHgIAAPv/8c8B84o+NjcXZ2Zmw\nsDDatWvHCy+8QHp6utnEn2/dunWMHj0aKP/vvlI1AFXxITGNRlPpf660tDRGjBjBRx99RM2aNQu9\nVtnjt7Cw4Pjx41y5coV9+/axZ8+eQq9X5vi3bt2Ki4sL/v7+903kXZnjB/j111+Jiopix44dLF26\nlP379xd6vTLHn5eXx7Fjx5g4cSLHjh3D3t6+WJdJZY4fICcnhx9++IFRo0YVe60ssVeqBqBhw4Zc\nvnxZv3358mXc3d1NGNHDcXV15caNGwBcv34dFxcXE0d0f7m5uYwYMYLnnnuOoUOHAuYVf75atWox\ncOBAIiMjzSb+3377jS1btuDh4cHo0aP55ZdfeO6558wmfgA3NzcAnJ2dGTZsGEeOHDGb+N3d3XF3\nd6dDhw4AjBw5kmPHjlG/fn2ziB9gx44dtG/fHmdnZ6D8f7uVqgEICAjg4sWLxMXFkZOTw/r16xk8\neLCpwyq3wYMHs2rVKgBWrVqlr1grG0VRGD9+PN7e3kydOlW/31ziv337tn6UQ2ZmJrt378bf399s\n4p89ezaXL18mNjaWdevW0atXL1avXm028WdkZJCamgpAeno6u3btwsfHx2zir1+/Po0aNeLChQsA\n/PTTT7Ru3ZqQkBCziB9g7dq1+u4feIi/XQPfnyi37du3Ky1atFA8PT2V2bNnmzqcUj311FOKm5ub\nYm1trbi7uyvLly9XEhMTld69eyvNmzdXgoKClKSkJFOHWaL9+/crGo1Gadu2reLn56f4+fkpO3bs\nMJv4T548qfj7+ytt27ZVfHx8lPfff19RFMVs4i8oIiJCCQkJURTFfOK/dOmS0rZtW6Vt27ZK69at\n9X+v5hK/oijK8ePHlYCAAMXX11cZNmyYkpycbDbxp6WlKXXr1lVSUlL0+8obuzwIJoQQ1VSl6gIS\nQghhPNIACCFENSUNgBBCVFPSAAghRDUlDYAQQlRT0gAIIUQ1JQ2AMGuBgYFERkYa/DyLFi3C29ub\n5557zuDnEsJYjJ4RTIiK9CjztOTl5WFlVbY/gWXLlunTmJZX/qM2lXlOGVE9yRWAMLi4uDhatWrF\nhAkTaNOmDX379iUrKwso/A3+9u3beHh4ALBy5UqGDh1KcHAwHh4eLFmyhAULFtCuXTueeOIJkpKS\n9J+/evVq/P398fHx4ffffwfUqQnCw8Pp1KkT7dq1Y8uWLfrPHTx4ML179yYoKKhYrB988AE+Pj74\n+Pjw0UcfAfDiiy9y6dIl+vXrx4cffljo+JUrVzJkyBB69uxJixYtePfdd/U/s5eXF2PHjsXHx4fL\nly/z2muv4ePjg6+vLxs2bNB/xrx58/D19cXPz0+f1CYmJob+/fsTEBBA9+7dOX/+PADffvstPj4+\n+Pn50aNHDwCio6Pp1KkT/v7+tG3blpiYGAC+/vpr/f4XX3wRnU6HVqtl3Lhx+jiK/jyimjHos8pC\nKGq+BCsrK+XEiROKoihKaGio8vXXXyuKoiiBgYFKZGSkoiiKkpCQoDRt2lRRFEVZsWKF0qxZMyUt\nLU1JSEhQHB0dlU8//VRRFEV55ZVXlA8//FBRFEXp0aOHMmHCBEVRFGXfvn36vAzTpk3TnyMpKUlp\n0aKFkp6erqxYsUJxd3cv8RH5o0ePKj4+PkpGRoaSlpamtG7dWjl+/LiiKOq894mJicXes2LFCsXN\nzU25c+eOkpmZqbRp00Y5evSoEhsbq1hYWCiHDx9WFEVRvvvuOyUoKEjR6XTKzZs3lcaNGyvXr19X\ntm/frnTp0kXJzMzUx6ooitKrVy/l4sWLiqIoyqFDh5RevXopiqIoPj4+yrVr1xRFUZS7d+8qiqIo\nkydPVr755htFURQlNzdXyczMVM6cOaOEhIQoeXl5iqIoysSJE5WvvvpKiYyMVIKCgvTxJycnl+Wf\nUFRR0gUkjMLDwwNfX18A2rdvT1xcXKnv6dmzJ/b29tjb2+Pk5KRPmejj48PJkycBtVslfzKsbt26\nkZKSwt27d9m1axc//PADCxYsACA7O5v4+Hg0Gg1BQUE4OTkVO9+BAwcYPnw4dnZ2AAwfPpx9+/bR\ntm3bB8YZHBxM7dq19e85cOAAQ4cOpUmTJvo0ib/++itPP/00Go0GFxcXevTowe+//87evXsJDw/H\n1tYWACcnJ9LS0jh48GChKX5zcnIA6Nq1K2PHjiU0NJThw4cD8MQTTzBr1iyuXLnC8OHDadasGT//\n/DORkZEEBAQA6mR5rq6uhISEcOnSJaZMmcLAgQMJDg4u9d9BVF3SAAijqFGjhn7d0tJS3wVkZWWF\nTqcD0O8r6T0WFhb6bQsLC/Ly8u57rvy+9o0bN9K8efNCrx0+fBh7e/v7vk8pMDWWoiil9tsXfV1R\nFH1yl6LnUe4z7VbR/TqdDicnJ6Kiooodu2zZMo4cOcK2bdto3749kZGRjB49ms6dO7N161YGDBjA\np59+CsDYsWOZPXt2sc84efIkO3fu5JNPPmHDhg18+eWXD/wZRdUl9wCESeRXek2bNuXo0aMAfPfd\nd+V6b/76+vXrAfUbvJOTE46OjvTt21ef5B7QV6b3q4RBvYLYtGkTmZmZpKens2nTJrp161ZqLLt3\n7yYpKYnMzEw2b95M165di52nW7durF+/Hp1OR0JCAvv27aNTp04EBQWxYsUKMjMzAUhKSsLR0REP\nDw/970NRFP0VT0xMDB07dmTGjBk4Oztz5coVYmNjadq0KZMnT2bIkCGcOnWK3r17891335GQkACo\nidrj4+NJTEwkLy+P4cOHM3PmTI4dO1am37momuQKQBhF0W/K+dv//Oc/CQ0N5bPPPmPgwIH6/UWz\nGRVdL3icra0t7dq1Iy8vj+XLlwPw9ttvM3XqVHx9fdHpdDz++ONs2bLlgVmS/P39GTdunL7b5oUX\nXtB3/9zvPRqNho4dOzJixAiuXLnCc889R7t27YiLiyv0nmHDhnHw4EHatm2LRqNh/vz5uLi40Ldv\nX44fP05AQAA2NjYMHDiQ9957j2+++YaXXnqJ9957j9zcXEaPHo2vry+vv/46Fy9eRFEU+vTpg6+v\nL/PmzWP16tVYW1vj5ubGW2+9hZOTE++99x7BwcHodDqsra35+OOPsbW1JSwsTH/VVRkTngvjkemg\nhXgEK1euJDIyksWLF5s6FCHKTbqAhHgElT1nrBAPIlcAQghRTckVgBBCVFPSAAghRDUlDYAQQlRT\n0gAIIUQ1JQ2AEEJUU9IACCFENfX/1gw3qwQPh7sAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x40ca910>"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Benefits\n",
      "--------\n",
      "\n",
      "There are still good reasons to parallelise your code:\n",
      "\n",
      "* Chips are getting more cores, not faster\n",
      "    + Clusters are getting bigger with faster interconnects\n",
      "    + GPGPUs give you a many cores on a single computer\n",
      "* Your problem easily scales perfectly (**embarrassingly parallel**) and can utilise all the resources effectively\n",
      "* Your problem would otherwise take years\n",
      "* An algorithm has been specifically designed for parallel computing\n",
      "* Take advantage of distributed hardware resources such as memory\n",
      "* To take advantage of **concurrency** in asynchronous problems with slow responses\n",
      "* You know that you will be able to maintain a parallel implementation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Zero code parallelisation\n",
      "-------------------------\n",
      "\n",
      "If you can manually chop up your data into segments that you can process entirely individually then you can run each of those instances of your code concurrently without having to recode anything. This can help to think about some concepts of parallel computing:\n",
      "\n",
      "* How do you divide the data up evenly so that each process takes a similar amount of time?\n",
      "    + Long running processes can leave idle CPUs until everything is finished\n",
      "    + Short processes where there are many more than the cores you have available will fill in gaps\n",
      "* How do the processes communicate with each other?\n",
      "    + How are the data and parameters supplied to each instance?\n",
      "    + How are the data recombined at the end?\n",
      "\n",
      "Maybe the easiest solution to your problem is to have a bunch of serial jobs run with a shell script (or `subprocess`) or submitted through a queuing system."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "![parallel tasks](files/Images/Par001.png)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Parallel programming\n",
      "--------------------\n",
      "\n",
      "Two very general classes of parallelisation are **shared memory** and **distributed memory** computing."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "![shared and distributed memory](files/Images/Par002.png)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Shared memory\n",
      "-------------\n",
      "\n",
      "Shared memory is usually achieved through **lightweight threads** and is very well standardised in compiled languages through **OpenMP**. It is designed to run on a single machine (or a shared memory cluster).\n",
      "\n",
      "Shared memory can be _difficult_ to do effectively manually as you need to avoid threads working on the same data or depending on the results of other threads. \n",
      "\n",
      "Compiling your code with `-openmp -ftree-parallelize-loops=n` in GCC 4.4 or later (where n is the number of cores you will be using), or `-openmp -parallel` with intel compilers will do some automatic code analysis and produce parallel binaries. In practice you can also to tell the compiler what to parallelise.\n",
      "\n",
      "* `!$omp parallel` is used in fortran\n",
      "* `#pragma omp parallel` is used in C\n",
      "* Consider private (loop variable) and public (everything) scopes\n",
      "* Consider workload and static dynamic schedulers\n",
      "* Setting `OMP_NUM_THREADS` at runtime to tweak performance\n",
      "\n",
      "There is also a lot of information available on optmiising your **OpenMP directives**, ask me if you'd like the notes from the HPCVL parallel workshop."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Distributed memory\n",
      "------------------\n",
      "\n",
      "The standard for distributed memory computing is **MPI** and it is intended to run on any number of cores across any number of nodes on most high performance computing (HPC) facilities. It is generally more more complicated than OpenMP as you have to deal with **communication** via **message passing**. Running MPI works by having multiple instances of the same program running in different places and communicating on the same computer or over a network in order to decide what work each should do.\n",
      "\n",
      "A basic MPI program can be implemented with just the following commands (but there are over a hundred):\n",
      "\n",
      "* `MPI_INIT` sets up everything needed to start; must be executed before enything else\n",
      "* `MPI_COMM_SIZE` finds out how many processes are in the MPI group\n",
      "* `MPI_COMM_RANK` determines the identifier of this particular process\n",
      "* `MPI_SEND` sends a *message*\n",
      "* `MPI_RECV` recieves a *message*\n",
      "* `MPI_FINALIZE` closes all the MPI stuff.\n",
      "\n",
      "Each running copy of the code is a completely independent process so the distribution of data amongst the processes has to be considered:\n",
      "\n",
      "* **Decomposition** has only parts of the data with each process but there needs to be processing to slice up the data and ensure that there is overlap between regions if necessary.\n",
      "* **Duplication** has all the data on every node which can make code faster with less communication but uses more memory.\n",
      "\n",
      "Often the `RANK 0` process will do reading, writing and data distribution and management as it is always guaranteed to exist (there are parallel IO codes but they are for heavy data processing like CERN)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Parallel programming in Python\n",
      "------------------------------\n",
      "\n",
      "Each Python process is strictly bound to a single thread by the global interface lock (GIL). This means that a Python process cannot use multi-core machines, but everything is **thread safe**, so if you use Python's `threading` library you won't see speedup with CPU bound tasks but sharing between threads is easy.\n",
      "\n",
      "There are a number of ways to achieve true **concurrency** either through libraries or multiple processes. Since version 2.6, Python's `multiprocessing` module, for example, has given an API for easy multiple process management that sidesteps the GIL."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file f_group_hist.py\n",
      "\n",
      "from random import random, choice\n",
      "\n",
      "FUNC = [\"H\", \"C\", \"N\", \"O\", \"F\", \"I\", \"P\", \"S\"]\n",
      "\n",
      "class FunctionalGroup(object):\n",
      "    \"\"\"Acts like a particle with a random history.\"\"\"\n",
      "    id_count = 0\n",
      "    def __init__(self, age=2000):\n",
      "        \"\"\"Initialise a random name and history.\"\"\"\n",
      "        self.position_history = []\n",
      "        self.generate_trajectory(length=age)\n",
      "        FunctionalGroup.id_count += 1\n",
      "        self.id = \"<Func|%-5s|%i>\" % (\n",
      "            \"\".join(choice(FUNC) for _i in range(5)),\n",
      "            FunctionalGroup.id_count)\n",
      "    def generate_trajectory(self, length=2000):\n",
      "        \"\"\"Fill it with random positions beacuse this is not real.\"\"\"\n",
      "        for _timestep in range(length):\n",
      "            self.position_history.append([random(), random(), random()])\n",
      "    def __repr__(self):\n",
      "        \"\"\"Help distinguish between instances.\"\"\"\n",
      "        return self.id\n",
      "\n",
      "def euclidean_dist(point_one, point_two):\n",
      "    \"\"\"Calculate the distance between two 3D points.\"\"\"\n",
      "    return ((point_one[0] - point_two[0])**2 +\n",
      "            (point_one[1] - point_two[1])**2 +\n",
      "            (point_one[2] - point_two[2])**2)**0.5\n",
      "\n",
      "def consecutive_pairs(iterable):\n",
      "    \"\"\"Generator that yields tuples containing the previous two items.\"\"\"\n",
      "    iter_object = iter(iterable)\n",
      "    # Take the first one so we always have two to give\n",
      "    previous = iter_object.next()\n",
      "    for item in iter_object:\n",
      "        yield (previous, item)\n",
      "        previous = item\n",
      "\n",
      "def distance_travelled(particle):\n",
      "    \"\"\"Total distance travelled for a particle over over a .position_history.\"\"\"\n",
      "    return sum(euclidean_dist(from_position, to_position) \n",
      "               for from_position, to_position \n",
      "               in consecutive_pairs(particle.position_history))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Overwriting f_group_hist.py\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from f_group_hist import FunctionalGroup, distance_travelled"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "my_groups = [FunctionalGroup() for _null in range(512)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Serial\n",
      "------\n",
      "\n",
      "For our example we will apply a distance calcaultor to a set of functional groups that have been moving around. In the first instance, we can do this serially."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def serial_distancer(groups):\n",
      "    \"\"\"Construct a dictionary with all the distances travelled.\"\"\"\n",
      "    return {group.id: distance_travelled(group) for group in groups}"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit serial_distancer(my_groups)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 1.57 s per loop\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Multiprocessing\n",
      "---------------\n",
      "\n",
      "The `multiprocessing` module spawns new processes that act similar to threads, but they are a bit more heavyweight and each can run on a separate core (outside the GIL). You need to do a bit of management for processes, but there are classes like `Pool` and `Queue` that help take care of distribution of work and **locking**."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import multiprocessing\n",
      "import os\n",
      "\n",
      "# A simple example of a pool\n",
      "\n",
      "def worker(args):\n",
      "    \"\"\"Give some info on the system.\"\"\"\n",
      "    print(\"PID = {}; args = {}\".format(os.getpid(), args))\n",
      "    return os.getpid(), args"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# A Pool manages all your processes\n",
      "pool = multiprocessing.Pool(processes=4)\n",
      "result = pool.map(worker, [1, 2, \"string\", [\"li\", \"st\"], 2.2, worker])\n",
      "result"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "PID = 21009; args = string\n",
        "PID = 21007; args = 1\n",
        "PID = 21008; args = 2\n",
        "PID = 21010; args = ['li', 'st']\n",
        "PID = 21008; args = 2.2\n",
        "PID = 21007; args = <function worker at 0x420d410>\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "[(21007, 1),\n",
        " (21008, 2),\n",
        " (21009, 'string'),\n",
        " (21010, ['li', 'st']),\n",
        " (21008, 2.2),\n",
        " (21007, <function __main__.worker>)]"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# run on our Pool of 4 workers\n",
      "%timeit pool.map(distance_travelled, my_groups)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 3.49 s per loop\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Queue\n",
      "-----\n",
      "Using  `Pool.map` with our function didn't show any benefit. The distribution of work is probably not optimised, instead all the data goes to every process, which is slow in this case. We can also manually create `Process`es and a `Queue` to demonstrate how we can distribute work and use shared objects for more complex output."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import multiprocessing\n",
      "from multiprocessing import Queue\n",
      "\n",
      "def mp_distancer(groups, nprocs):\n",
      "    def worker(groups, out_q):\n",
      "        \"\"\"\n",
      "        Worker to be invoked in a process. Pass the groups to be\n",
      "        measured and the results dictionary is pushed to a queue.\n",
      "        \"\"\"\n",
      "        outdict = {}\n",
      "        for group in groups:\n",
      "            outdict[group.id] = distance_travelled(group)\n",
      "        out_q.put(outdict)\n",
      "\n",
      "    # Each process gets a chunk and the queue to put results into\n",
      "    out_q = multiprocessing.Queue()\n",
      "    # Divide up the work into equally sized chunks\n",
      "    chunksize = int(math.ceil(len(groups) / float(nprocs)))\n",
      "    # A list stores references to our processes\n",
      "    procs = []\n",
      "\n",
      "    for proc_idx in range(nprocs):\n",
      "        work = groups[chunksize*proc_idx:chunksize*(proc_idx+1)]\n",
      "        # Create each process using the function to execute and \n",
      "        # the tuple of arguments with work and a shared queue\n",
      "        p = multiprocessing.Process(target=worker, args=(work, out_q))\n",
      "        procs.append(p)\n",
      "        p.start()\n",
      "\n",
      "    # Collect all results into a single dict. We know how many to expect.\n",
      "    result_dict = {}\n",
      "    for _i in range(nprocs):\n",
      "        # The get will block until there is an item to get\n",
      "        result_dict.update(out_q.get())\n",
      "\n",
      "    # Wait for all worker processes to finish\n",
      "    # Queue must be empty before joining to prevent deadlocks\n",
      "    for p in procs:\n",
      "        p.join()\n",
      "\n",
      "    return result_dict"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mp_distancer(my_groups, 2) == serial_distancer(my_groups)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit mp_distancer(my_groups, 1)\n",
      "%timeit mp_distancer(my_groups, 2)\n",
      "%timeit mp_distancer(my_groups, 8)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 1.76 s per loop\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 907 ms per loop\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 317 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Fianlly we see some speedup, athough it is far from *perfect scaling*. In every case there is overhead from process setup and communication which makes the single core version slower than the serial version. You must decide whether it is worth your time to make things parallel. \n",
      "\n",
      "The processes can share memory, although in this case you need to be careful to lock variables as you manipulate them."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "MPI\n",
      "----\n",
      "\n",
      "Although `multiprocessing` can be made to work across nodes on clusters, Python also has wrappers around standard MPI libraries. With MPI, many instances of your code run completely separately and communicate **messages** to each other. In this example we use `scatter` and `collect` but a similar thing could be achieved with `send` and `recv`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file mpi_distancer.py\n",
      "\n",
      "#!/usr/bin/env python\n",
      "\n",
      "\"\"\"\n",
      "mpi_distancer.py\n",
      "\n",
      "Simple implementation of a history distance calculation\n",
      "with data scatter and gather via MPI.\n",
      "\n",
      "\"\"\"\n",
      "\n",
      "import math\n",
      "import time\n",
      "\n",
      "# The import takes care of MPI_Init and setup\n",
      "from mpi4py import MPI\n",
      "\n",
      "# Local definitions that we'll use\n",
      "from f_group_hist import FunctionalGroup, distance_travelled\n",
      "\n",
      "# comm is the communicator that handles our collection of MPI tasks\n",
      "comm = MPI.COMM_WORLD\n",
      "\n",
      "# These identify how many processes are running and which one this is\n",
      "rank, size = comm.rank, comm.size\n",
      "print \"Running as rank %d out of %d\" % (rank, size)\n",
      "\n",
      "if rank == 0:\n",
      "    # treat 0 as the root node as it is always there\n",
      "    start_time = time.time()\n",
      "    # Instantiate all the data; \n",
      "    # if data input is time consuming you can distribute this step too\n",
      "    my_groups = [FunctionalGroup() for _i in range(511)]\n",
      "    print(\"Generated data in {} seconds\".format(time.time() - start_time))\n",
      "    chunksize = int(math.ceil(len(my_groups) / float(size)))\n",
      "    # We need to partition the groups ourselves (not always necessary)\n",
      "    my_groups = [my_groups[chunksize*i:chunksize*(i+1)] for i in range(size)]\n",
      "else:\n",
      "    # All other processes don't care\n",
      "    my_groups = None\n",
      "\n",
      "# scatter sends out pickled python objects to all the other processes\n",
      "my_groups = comm.scatter(my_groups, root=0)\n",
      "\n",
      "if rank == 0:\n",
      "    print(\"Sent data after {} seconds\".format(time.time() - start_time))\n",
      "    # using NumPy arrays with Scatter will be faster\n",
      "\n",
      "print(\"Rank {} has {} items; first is {}\".format(rank, len(my_groups), my_groups[0]))\n",
      "\n",
      "# Here's the processing bit that we do\n",
      "results = {group.id: distance_travelled(group) for group in my_groups}\n",
      "\n",
      "print(\"Rank {} done\".format(rank))\n",
      "\n",
      "# Wait for everybody to synchronize _here_\n",
      "comm.Barrier()\n",
      "\n",
      "# collects the results as a list with results\n",
      "# for each node in a separate item\n",
      "results = comm.gather(results, root=0)\n",
      "\n",
      "if rank == 0:\n",
      "    # Combine all our data together as before, but only on the root node\n",
      "    final_results = {}\n",
      "    for result in results:\n",
      "        final_results.update(result)\n",
      "\n",
      "    # All done\n",
      "    print(\"Total time: {} seconds\".format(time.time() - start_time))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Overwriting mpi_distancer.py\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# mpiexec starts as many pythons as we need\n",
      "!mpiexec -n 1 python mpi_distancer.py"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Running as rank 0 out of 1\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Generated data in 1.97901010513 seconds\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Sent data after 8.92737817764 seconds\r\n",
        "Rank 0 has 511 items; first is <Func|FFFSP|1>\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Rank 0 done\r\n",
        "Total time: 10.5981180668 seconds\r\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!mpiexec -n 4 python mpi_distancer.py"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Running as rank 1 out of 4\r\n",
        "Running as rank 2 out of 4\r\n",
        "Running as rank 0 out of 4\r\n",
        "Running as rank 3 out of 4\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Generated data in 1.98770403862 seconds\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Rank 1 has 128 items; first is <Func|NFNFS|129>\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Rank 2 has 128 items; first is <Func|PSFNO|257>\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Rank 3 has 127 items; first is <Func|FSFPH|385>\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Rank 1 done\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Rank 2 done\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Rank 3 done\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Sent data after 6.26479887962 seconds\r\n",
        "Rank 0 has 128 items; first is <Func|OPSOH|1>\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Rank 0 done\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Total time: 6.71548199654 seconds\r\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We see some speedup as the calculation part scales very well, but calcaultion setup and distribution of big chunks of data can take a significant amount of time and reduce the efficiency of the parallelisation.\n",
      "\n",
      "MPI parallelisation is best when the communication is kept to a minimum."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Unlocking the GIL with Cython and OpenMP\n",
      "----------------------------------------\n",
      "\n",
      "As we have seen, we can incorporate code in other languages in Python. Depending on how you compile that external code, you can take advantage of many-core machines. Since Cython translates to C, there are features that let you leverage OpenMP parallelisation from your Cython modules."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext cythonmagic"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython -f -c-fopenmp --link-args=-fopenmp -c-g\n",
      "# enable extra compiler flags for openmp compilation\n",
      "\n",
      "import numpy\n",
      "cimport cython\n",
      "cimport numpy\n",
      "# use prange for automatic work distribution\n",
      "from cython.parallel cimport prange\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "def cy_euclid_par(numpy.ndarray[numpy.float64_t, ndim=2] coordinates,\n",
      "                  int num_threads=1):\n",
      "    \"\"\"Calculate euclidean distance of all coordinates.\"\"\"\n",
      "    cdef int idx, size\n",
      "    cdef numpy.ndarray[numpy.float64_t, ndim=1] output \n",
      "\n",
      "    output = numpy.empty(coordinates.shape[0])\n",
      "    size = coordinates.shape[0]\n",
      "    \n",
      "    # code block releases the gil\n",
      "    with nogil:\n",
      "        # this block runs separately for each OpenMP thread\n",
      "        # and prange gives idx different values per thread\n",
      "        for idx in prange(size, schedule='static',\n",
      "                          num_threads=num_threads):\n",
      "            # communicate with shared buffers (not Python objects)\n",
      "            output[idx] = (coordinates[idx, 0]**2 +\n",
      "                           coordinates[idx, 1]**2 +\n",
      "                           coordinates[idx, 2]**2)**0.5\n",
      "    return output"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n",
      "positions = numpy.random.random((3000, 3))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit cy_euclid_par(positions)\n",
      "%timeit cy_euclid_par(positions, num_threads=2)\n",
      "%timeit cy_euclid_par(positions, num_threads=4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 586 \u00b5s per loop\n",
        "1000 loops, best of 3: 302 \u00b5s per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "10000 loops, best of 3: 157 \u00b5s per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "A few notes:\n",
      "\n",
      "+ Speedup is good in this simple case as all the heavy processing in inside the loop\n",
      "+ Can't interact with Python objects within `nogil` regions\n",
      "    + Everything needs to be typed and can get complicated\n",
      "    + Cython gives helpful errors on compilation to prevent this\n",
      "    + use `with gil` blocks for python interaction\n",
      "+ Need to use NumPy `array`s for a shared memory buffer (Python types don't work)\n",
      "+ `prange` splits the work up on a schedule (these are OpenMP concepts - see figure):\n",
      "    + `static` cuts work into equal chunks; no overhead; good where each job takes the same time.\n",
      "    + `dynamic` work is given to threads when they request it in chunks of one; most overhead; good where jobs take different times or threads are different speeds\n",
      "    + `guided` threads are initially given a large chunk with decreasing size chunks as requested by threads; some overhead; may or may not be better than `dynamic`\n",
      "+ Work can also be distributed manually\n",
      "\n",
      "![scheduling work](files/Images/distribute_work.png)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "GPU\n",
      "---\n",
      "\n",
      "Python has modules for CUDA and PyOpenCL for GPU massively parallel computing. These work a lot like `weave` that we saw previously in that you compile the functions that run on the GPU but you also need to transfer your data to and from the GPU efficiently. They look to be fairly well documented, but that is an exercise left to the reader."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pycuda.autoinit\n",
      "import pycuda.driver as drv\n",
      "import numpy\n",
      "\n",
      "from pycuda.compiler import SourceModule\n",
      "mod = SourceModule(\"\"\"\n",
      "__global__ void multiply_them(float *dest, float *a, float *b)\n",
      "{\n",
      "  const int i = threadIdx.x;\n",
      "  dest[i] = a[i] * b[i];\n",
      "}\n",
      "\"\"\")\n",
      "\n",
      "multiply_them = mod.get_function(\"multiply_them\")\n",
      "\n",
      "a = numpy.random.randn(400).astype(numpy.float32)\n",
      "b = numpy.random.randn(400).astype(numpy.float32)\n",
      "\n",
      "dest = numpy.zeros_like(a)\n",
      "multiply_them(\n",
      "        drv.Out(dest), drv.In(a), drv.In(b),\n",
      "        block=(400,1,1), grid=(1,1))\n",
      "\n",
      "print a[0], b[0], a[0]*b[0], dest[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-0.243131 0.078306 -0.0190386 -0.0190386\n"
       ]
      }
     ],
     "prompt_number": 20
    }
   ],
   "metadata": {}
  }
 ]
}