{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tools for high-performance computing applications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "J.R. Johansson (jrjohansson at gmail.com)\n",
    "\n",
    "The latest version of this [IPython notebook](http://ipython.org/notebook.html) lecture is available at [http://github.com/jrjohansson/scientific-python-lectures](http://github.com/jrjohansson/scientific-python-lectures).\n",
    "\n",
    "The other notebooks in this lecture series are indexed at [http://jrjohansson.github.io](http://jrjohansson.github.io)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## multiprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python has a built-in process-based library for concurrent computing, called `multiprocessing`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiprocessing\n",
    "import os\n",
    "import time\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def task(args):\n",
    "    print((\"PID =\", os.getpid(), \", args =\", args))\n",
    "    \n",
    "    return os.getpid(), args"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PID = 28995 , args = test\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(28995, 'test')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "task(\"test\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pool = multiprocessing.Pool(processes=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PID = 29006 , args = 1\n",
      "PID = 29009 , args = 4\n",
      "PID = 29007 , args = 2\n",
      "PID = 29008 , args = 3\n",
      "PID = 29006 , args = 6\n",
      "PID = 29009 , args = 5\n",
      "PID = 29007 , args = 8\n",
      "PID = 29008 , args = 7\n"
     ]
    }
   ],
   "source": [
    "result = pool.map(task, [1,2,3,4,5,6,7,8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(29006, 1),\n",
       " (29007, 2),\n",
       " (29008, 3),\n",
       " (29009, 4),\n",
       " (29009, 5),\n",
       " (29006, 6),\n",
       " (29008, 7),\n",
       " (29007, 8)]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The multiprocessing package is very useful for highly parallel tasks that do not need to communicate with each other, other than when sending the initial data to the pool of processes and when and collecting the results. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IPython parallel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IPython includes a very interesting and versatile parallel computing environment, which is very easy to use. It builds on the concept of ipython engines and controllers, that one can connect to and submit tasks to. To get started using this framework for parallel computing, one first have to start up an IPython cluster of engines. The easiest way to do this is to use the `ipcluster` command,\n",
    "\n",
    "    $ ipcluster start -n 4\n",
    "\n",
    "Or, alternatively, from the \"Clusters\" tab on the IPython notebook dashboard page. This will start 4 IPython engines on the current host, which is useful for multicore systems. It is also possible to setup IPython clusters that spans over many nodes in a computing cluster. For more information about possible use cases, see the official documentation [Using IPython for parallel computing](http://ipython.org/ipython-doc/dev/parallel/).\n",
    "\n",
    "To use the IPython cluster in our Python programs or notebooks, we start by creating an instance of `IPython.parallel.Client`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.parallel import Client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "cli = Client()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the 'ids' attribute we can retreive a list of ids for the IPython engines in the cluster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 1, 2, 3]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cli.ids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of these engines are ready to execute tasks. We can selectively run code on individual engines:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getpid():\n",
    "    \"\"\" return the unique ID of the current process \"\"\"\n",
    "    import os\n",
    "    return os.getpid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "28995"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# first try it on the notebook process\n",
    "getpid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30181"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# run it on one of the engines\n",
    "cli[0].apply_sync(getpid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[30181, 30182, 30183, 30185]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# run it on ALL of the engines at the same time\n",
    "cli[:].apply_sync(getpid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this cluster of IPython engines to execute tasks in parallel. The easiest way to dispatch a function to different engines is to define the function with the decorator:\n",
    "\n",
    "    @view.parallel(block=True)\n",
    "\n",
    "Here, `view` is supposed to be the engine pool which we want to dispatch the function (task). Once our function is defined this way we can dispatch it to the engine using the `map` method in the resulting class (in Python, a decorator is a language construct which automatically wraps the function into another function or a class).\n",
    "\n",
    "To see how all this works, lets look at an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "dview = cli[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "@dview.parallel(block=True)\n",
    "def dummy_task(delay):\n",
    "    \"\"\" a dummy task that takes 'delay' seconds to finish \"\"\"\n",
    "    import os, time\n",
    "\n",
    "    t0 = time.time()\n",
    "    pid = os.getpid()\n",
    "    time.sleep(delay)\n",
    "    t1 = time.time()\n",
    "    \n",
    "    return [pid, t0, t1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate random delay times for dummy tasks\n",
    "delay_times = numpy.random.rand(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to map the function `dummy_task` to the random delay time data, we use the `map` method in `dummy_task`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[30181, 1395044753.2096598, 1395044753.9150908],\n",
       " [30182, 1395044753.2084103, 1395044753.4959202],\n",
       " [30183, 1395044753.2113762, 1395044753.6453338],\n",
       " [30185, 1395044753.2130392, 1395044754.1905618]]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dummy_task.map(delay_times)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do the same thing again with many more tasks and visualize how these tasks are executed on different IPython engines:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def visualize_tasks(results):\n",
    "    res = numpy.array(results)\n",
    "    fig, ax = plt.subplots(figsize=(10, res.shape[1]))\n",
    "    \n",
    "    yticks = []\n",
    "    yticklabels = []\n",
    "    tmin = min(res[:,1])\n",
    "    for n, pid in enumerate(numpy.unique(res[:,0])):\n",
    "        yticks.append(n)\n",
    "        yticklabels.append(\"%d\" % pid)\n",
    "        for m in numpy.where(res[:,0] == pid)[0]:\n",
    "            ax.add_patch(plt.Rectangle((res[m,1] - tmin, n-0.25),\n",
    "                         res[m,2] - res[m,1], 0.5, color=\"green\", alpha=0.5))\n",
    "        \n",
    "    ax.set_ylim(-.5, n+.5)\n",
    "    ax.set_xlim(0, max(res[:,2]) - tmin + 0.)\n",
    "    ax.set_yticks(yticks)\n",
    "    ax.set_yticklabels(yticklabels)\n",
    "    ax.set_ylabel(\"PID\")\n",
    "    ax.set_xlabel(\"seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "delay_times = numpy.random.rand(64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAADSCAYAAAAYCPc3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHfdJREFUeJzt3X9wFPX9x/HXHQFiSpofChcN1jtQhIRL7vhhxoI1MQRE\nKqLlh4DkKvSHdSaMjkiLdNra4RvCSCmpaB07FNLaAZk6jVEq7aVNwEgtmh8HHYNYyg2gScaAIQkE\nI7n7/sF4YwgJv3K32eT5mMmMt7f72ff7czfh5W521xIMBoMCAABAn2Y1ugAAAABcGqENAADABAht\nAAAAJkBoAwAAMAFCGwAAgAlEGV1AuFksFqNLAAAAuGzd3dhjQBxpCwaDA/bn5z//ueE1GP0z0OeA\n/unf6Bron/7p//J/ejIgQhsAAIDZEdoAAABMgNDWz2VmZhpdguEG+hzQf6bRJRiK/jONLsFQ9J9p\ndAm9yhK81AlUk7NYLJc8RwwAANAX9JRbONIGAABgAoQ2AAAAEyC0AQAAmAChDQAAwAQIbQAAACbQ\n7x9jJUlP7HpCTWebIrY/X71P6Unpphk3kvpDDz0J92dv1PyZ+XMzc+2Xcrm9RWIOrmYfvnqfJBn+\n+fAdMVZ8dLw23rvR6DJMYUCEtqazTbLH2yO2v4qjFWHZX7jGjaT+0ENPwv3ZGzV/Zv7czFz7pVxu\nb5GYg6vZR8XRCkky/PPhO2Isf5Pf6BJMg9OjAAAAJkBoAwAAMAFCGwAAgAkQ2gAAAEwgbKHt7Nmz\nysjIkMvlUkpKilatWiVJOnnypHJycjRmzBhNnz5dTU1NoeVZWVmKjY1VXl5ep7G2bNkip9Op9PR0\nzZw5UydOnJAkbd26VcOHD5fb7Zbb7dbvf//7cLUDAABgqLCFtujoaJWVlammpkb79+9XWVmZKioq\nVFBQoJycHB06dEjZ2dkqKCgIrb9mzRqtX7++0zjt7e1asWKFdu/eLZ/Pp7S0NG3atEnS+YeqLly4\nUNXV1aqurtbSpUvD1Q4AAIChwnp6NCYmRtL54NXR0aGEhASVlJTI4/FIkjwej4qLi0PrTpkyRUOH\nDu00RlRUlBISEtTa2qpgMKhTp04pOTlZkhQMBhUMBsPZAgAAQJ8Q1tAWCATkcrlks9mUlZWl1NRU\nNTQ0yGazSZJsNpsaGho6bWOxWDoXaLWqsLBQ48ePV3Jysg4ePBg6omaxWPTaa68pLS1N8+bN0/Hj\nx8PZDgAAgGHCenNdq9WqmpoanTp1SjNmzFBZWVmn9y0WS5eQdqHm5mYtX75cPp9PDodDeXl5Wrt2\nrVavXq37779fixYt0uDBg/Xyyy/L4/HoH//4R5cxarbXyB/tlyTZXXbZXfbeahEAAOCqlZeXq7y8\n/LLWjcgTEeLi4jRr1ixVVlbKZrOpvr5eSUlJqqur04gRI3rctra2Vg6HQw6HQ5I0b948rVu3TpKU\nmJgYWm/ZsmVauXLlRcdwPezq83eEBgAAA09mZqYyMzNDr5999tlu1w3b6dHGxsbQlaFtbW3yer1y\nu92aPXu2ioqKJElFRUWaM2dOp+0u/Bu1UaNG6eDBg2psbJQkeb1epaSkSJLq6upC65WUlISWAwAA\n9DdhO9JWV1cnj8ejQCCgQCCgJUuWKDs7W263W/Pnz9fmzZtlt9u1Y8eO0DZ2u10tLS1qb29XcXGx\nvF6vxo4dq/z8fGVlZclqtcput2vr1q2SpOeff14lJSWKiorS9ddfH1oOAADQ34QttDmdTlVVVXVZ\nnpiYqNLS0otu4/f7L7o8NzdXubm5XZbn5+crPz//muoEAAAwA56IAAAAYAKENgAAABMgtAEAAJgA\noQ0AAMAECG0AAAAmQGgDAAAwgYg8EcFo8dHx8jf5I7a/2CGxYdlfuMaNpP7QQ0/C/dkbNX9m/tzM\nXPulXG5vkZiDq9lH7JBYSTL88+E7Yqz46HijSzANS/DCRxD0MxaLpctTFgAAAPqinnILp0cBAABM\ngNAGAABgAoQ2AAAAEyC0AQAAmAChDQAAwAQIbQAAACZAaAMAADABQhsAAIAJENoAAABMgNAGAABg\nAoQ2AAAAEyC0AQAAmAChDQAAwAQIbQAAACZAaAMAADABQhsAAIAJENoAAABMgNAGAABgAoQ2AAAA\nE4gyuoBIeGLXE2o629Tt+756n9KT0iNYUfhcay9Xun1vz91Xx+tp7L74mYWjpqsd0+zfg6tRcbRC\nU78x1dAaelN3cxofHa+N9240oCIARhsQoa3pbJPs8fZu3684WtHj+2Zyrb1c6fa9PXdfHa+nsfvi\nZxaOmq52TLN/D65GyYclhtfQm7qbU3+TP+K1AOgbOD0KAABgAoQ2AAAAEyC0AQAAmAChDQAAwATC\nFtrOnj2rjIwMuVwupaSkaNWqVZKkkydPKicnR2PGjNH06dPV1NQUWp6VlaXY2Fjl5eV1GmvLli1y\nOp1KT0/XzJkzdeLECUnSSy+9pLS0NLndbt15553y+XzhagcAAMBQYQtt0dHRKisrU01Njfbv36+y\nsjJVVFSooKBAOTk5OnTokLKzs1VQUBBaf82aNVq/fn2ncdrb27VixQrt3r1bPp9PaWlp2rRpkyRp\n8eLF2r9/v6qrq/XMM8/oqaeeClc7AAAAhgrr6dGYmBhJ54NXR0eHEhISVFJSIo/HI0nyeDwqLi4O\nrTtlyhQNHTq00xhRUVFKSEhQa2urgsGgmpublZycLEmKjY0Nrdfa2qobbrghnO0AAAAYJqz3aQsE\nApowYYIOHz6sH/3oR0pNTVVDQ4NsNpskyWazqaGhodM2Foul02ur1arCwkKNHz9ew4YN05gxY/TC\nCy+E3n/xxRe1YcMGnT59Wnv37g1nOwAAAIYJa2izWq2qqanRqVOnNGPGDJWVlXV632KxdAlpF2pu\nbtby5cvl8/nkcDiUl5entWvXavXq1ZKkxx9/XI8//ri2bdumpUuXdtmHJNVsr5E/2i9Jsrvssrvs\nvdIfAADAtSgvL1d5efllrRuRJyLExcVp1qxZqqyslM1mU319vZKSklRXV6cRI0b0uG1tba0cDocc\nDockad68eVq3bl2X9RYsWKDHHnvsomO4Hnb1qzulAwCA/iEzM1OZmZmh188++2y364btb9oaGxtD\nV4a2tbXJ6/XK7XZr9uzZKioqkiQVFRVpzpw5nbYLBoOdXo8aNUoHDx5UY2OjJMnr9SolJUWS9NFH\nH4XW27lzp9LS0sLVDgAAgKHCdqStrq5OHo9HgUBAgUBAS5YsUXZ2ttxut+bPn6/NmzfLbrdrx44d\noW3sdrtaWlrU3t6u4uJieb1ejR07Vvn5+crKypLVapXdbtfWrVslSS+88IJKS0s1ePBgDR8+XFu2\nbAlXOwAAAIYKW2hzOp2qqqrqsjwxMVGlpaUX3cbv9190eW5urnJzc7ss37hx4zXVCAAAYBY8EQEA\nAMAECG0AAAAmQGgDAAAwAUIbAACACRDaAAAATIDQBgAAYAIReSKC0eKj4+Vv8nf7fuyQ2B7fN5Nr\n7eVKt+/tufvqeD2N3Rc/s3DUdLVjmv17cDUSr0s0vIbe1N2cxkfHR74YAH2CJXjhIwj6GYvF0uUp\nCwAAAH1RT7mF06MAAAAmQGgDAAAwgR5D29atWzVhwgTFxMQoJiZGkyZNCj3sHQAAAJHT7YUIRUVF\nKiws1IYNG+R2uxUMBlVdXa2nn35aFovlos8CBQAAQHh0eyFCRkaGtm/fLofD0Wm53+/XggUL9O9/\n/zsiBV4rLkQAAABmcVUXIrS0tHQJbJJkt9vV0tLSe9UBAADgkroNbdHR0d1u1NN7AAAA6H3dnh69\n7rrrdOutt150o8OHD+vMmTNhLay3cHoUAACYRU+5pdsLEWpra8NWEAAAAK4MT0QAAADoI67qSNuw\nYcNksVi6HbC5ubl3qgMAAMAlcaQNAACgj7iqI21tbW166aWXdPjwYTmdTi1btkxRUd2uDgAAgDDq\n9kjb/PnzNWTIEE2dOlVvvfWW7Ha7CgsLI13fNeNIGwAAMIuecku3oc3pdOrAgQOSpHPnzmny5Mmq\nrq4OX5VhQmgDAABmcVVPRPjqqVBOiwIAABir2yNtgwYNUkxMTOh1W1ubrrvuuvMbmejqUY60AQAA\ns7iqCxE6OjrCVhAAAACuTLenRwEAANB3DIg/Vvtu8XeNLuGifPU+pSelX/M6ZtZTf5fq3YxzY8aa\nu3OtvUR6LiKxvyvZR1+cv778/TSitr48H1+Kj47Xxns3Gl0GImRAhDZ7vN3oEi6q4mjFJWu7nHXM\nrKf+LtW7GefGjDV351p7ifRcRGJ/V7KPvjh/ffn7aURtfXk+vuRv8htdAiKI06MAAAAmQGgDAAAw\nAUIbAACACRDaAAAATCBsoe3s2bPKyMiQy+VSSkqKVq1aJUk6efKkcnJyNGbMGE2fPl1NTU2h5VlZ\nWYqNjVVeXl6nsbZs2SKn06n09HTNnDlTJ06ckCRt2LBBqampSk9P17Rp03T06NFwtQMAAGCosIW2\n6OholZWVqaamRvv371dZWZkqKipUUFCgnJwcHTp0SNnZ2SooKAitv2bNGq1fv77TOO3t7VqxYoV2\n794tn8+ntLQ0bdq0SZI0YcIEVVZWyufzae7cuVq5cmW42gEAADBUWE+PfvkYrPb2dnV0dCghIUEl\nJSXyeDySJI/Ho+Li4tC6U6ZM0dChQzuNERUVpYSEBLW2tioYDKq5uVnJycmSpMzMTEVHR0uSMjIy\ndPz48XC2AwAAYJiwhrZAICCXyyWbzaasrCylpqaqoaFBNptNkmSz2dTQ0NBpG4vF0rlAq1WFhYUa\nP368kpOTVVtbq6VLl3bZ1+bNm3XfffeFrxkAAAADhfXmularVTU1NTp16pRmzJihsrKyTu9bLJYu\nIe1Czc3NWr58uXw+nxwOh/Ly8rR27VqtXr06tM4rr7yiqqoq/frXv77oGOVby0P/bXfZZXfZr7on\nAACA3lJeXq7y8vLLWjciT0SIi4vTrFmzVFlZKZvNpvr6eiUlJamurk4jRozocdva2lo5HA45HA5J\n0rx587Ru3brQ+6WlpcrPz9eePXs0ePDgi46R+d3MXusFAACgt2RmZiozMzP0+tlnn+123bCdHm1s\nbAxdGdrW1iav1yu3263Zs2erqKhIklRUVKQ5c+Z02i4YDHZ6PWrUKB08eFCNjY2SJK/Xq5SUFElS\ndXW1HnvsMb3xxhu64YYbwtUKAACA4cJ2pK2urk4ej0eBQECBQEBLlixRdna23G635s+fr82bN8tu\nt2vHjh2hbex2u1paWtTe3q7i4mJ5vV6NHTtW+fn5ysrKktVqld1u19atWyVJK1eu1OnTpzV37lxJ\n0i233BK6sAEAAKA/CVtoczqdqqqq6rI8MTFRpaWlF93G7/dfdHlubq5yc3O7LPd6vddUIwAAgFnw\nRAQAAAATILQBAACYAKENAADABAhtAAAAJkBoAwAAMIGI3FzXaP4mv9ElXFTskNhL1nY565hZT/1d\nqnczzo0Za+7OtfYS6bmIxP6uZB99cf768vfTiNr68nx8KT463ugSEEGW4IV3s+1nLBZLlxv2AgAA\n9EU95RZOjwIAAJgAoQ0AAMAECG0AAAAmQGgDAAAwAUIbAACACRDaAAAATIDQBgAAYAKENgAAABMg\ntAEAAJgAoQ0AAMAECG0AAAAmQGgDAAAwAUIbAACACRDaAAAATIDQBgAAYAKENgAAABMgtAEAAJgA\noQ0AAMAECG0AAAAmEGV0AZHwxK4n1HS2yegyIs5X71N6UrrRZURMuPrtzXEvHMtMn9GV1mqm3i6l\nP/XypXD31Ffn7Frq6qs9XS6z1m/m35u9bUCEtqazTbLH240uI+IqjlYMqL7D1W9vjnvhWGb6jK60\nVjP1din9qZcvhbunvjpn11JXX+3pcpm1fjP/3uxtnB4FAAAwAUIbAACACRDaAAAATIDQBgAAYAKE\nNgAAABMIW2g7e/asMjIy5HK5lJKSolWrVkmSTp48qZycHI0ZM0bTp09XU1NTaHlWVpZiY2OVl5fX\naawtW7bI6XQqPT1dM2fO1IkTJyRJe/bs0YQJEzR48GC99tpr4WoFAADAcGELbdHR0SorK1NNTY32\n79+vsrIyVVRUqKCgQDk5OTp06JCys7NVUFAQWn/NmjVav359p3Ha29u1YsUK7d69Wz6fT2lpadq0\naZMk6ZZbblFRUZEWLVoUrjYAAAD6hLCeHo2JiZF0Pnh1dHQoISFBJSUl8ng8kiSPx6Pi4uLQulOm\nTNHQoUM7jREVFaWEhAS1trYqGAyqublZycnJks6HNqfTKauVs7wAAKB/C2vaCQQCcrlcstlsysrK\nUmpqqhoaGmSz2SRJNptNDQ0NnbaxWCydC7RaVVhYqPHjxys5OVm1tbVaunRpOMsGAADoc8L6RASr\n1aqamhqdOnVKM2bMUFlZWaf3LRZLl5B2oebmZi1fvlw+n08Oh0N5eXlau3atVq9efdl11GyvkT/a\nL0myu+yyu+xX2goAAECv89f45a/xX9a6EXmMVVxcnGbNmqXKykrZbDbV19crKSlJdXV1GjFiRI/b\n1tbWyuFwyOFwSJLmzZundevWdVmvp/Dnetg1YB95AQAA+q4LDybtLtrd7bphOz3a2NgYujK0ra1N\nXq9Xbrdbs2fPVlFRkSSpqKhIc+bM6bRdMBjs9HrUqFE6ePCgGhsbJUler1cpKSldtrlwOwAAgP4k\nbEfa6urq5PF4FAgEFAgEtGTJEmVnZ8vtdmv+/PnavHmz7Ha7duzYEdrGbrerpaVF7e3tKi4ultfr\n1dixY5Wfn6+srCxZrVbZ7XZt3bpVkvTee+/poYce0meffaY333xTv/jFL3TgwIFwtQQAAGCYsIU2\np9OpqqqqLssTExNVWlp60W38fv9Fl+fm5io3N7fL8smTJ+vYsWPXVCcAAIAZcK8MAAAAEyC0AQAA\nmAChDQAAwAQIbQAAACZAaAMAADCBiNxc12jx0fHyN/mNLiPiYofEDqi+w9Vvb4574Vhm+oyutFYz\n9XYp/amXL4W7p746Z9dSV1/t6XKZtX4z/97sbZZgP78rrcVi4ca7AADAFHrKLZweBQAAMAFCGwAA\ngAkQ2gAAAEyA0AYAAGAChDYAAAATILT1c+Xl5UaXYLiBPgf0X250CYai/3KjSzAU/ZcbXUKvIrT1\nc/3tC3s1Bvoc0H+50SUYiv7LjS7BUPRfbnQJvYrQBgAAYAKENgAAABMYEE9EAAAAMIvuolm/f/Zo\nP8+kAABggOD0KAAAgAkQ2gAAAEyA0AYAAGAC/Tq07dq1S2PHjtVtt92mdevWGV1ORC1dulQ2m01O\np9PoUgxx7NgxZWVlKTU1VePHj9dvfvMbo0uKqLNnzyojI0Mul0spKSlatWqV0SUZoqOjQ263W/ff\nf7/RpUSc3W5XWlqa3G637rjjDqPLibimpibNnTtX48aNU0pKit59912jS4qYDz/8UG63O/QTFxc3\n4H4Hrl27VqmpqXI6nVq0aJE+//xzo0vqFf326tGOjg7dfvvtKi0tVXJysiZPnqxt27Zp3LhxRpcW\nEW+//baGDRum3NxcHThwwOhyIq6+vl719fVyuVxqbW3VxIkTVVxcPGA+f0k6c+aMYmJidO7cOU2d\nOlXr16/X1KlTjS4rojZs2KDKykq1tLSopKTE6HIiyuFwqLKyUomJiUaXYgiPx6O7775bS5cu1blz\n53T69GnFxcUZXVbEBQIBJScna9++fbr55puNLici/H6/7rnnHtXW1mro0KFasGCB7rvvPnk8HqNL\nu2b99kjbvn37dOutt8put2vw4MF6+OGH9frrrxtdVsTcddddSkhIMLoMwyQlJcnlckmShg0bpnHj\nxumTTz4xuKrIiomJkSS1t7ero6NjwP3jffz4cf31r3/V9773vQF7FflA7fvUqVN6++23tXTpUklS\nVFTUgAxsklRaWqrRo0cPmMAmSV//+tc1ePBgnTlzRufOndOZM2eUnJxsdFm9ot+Gto8//rjTl3Tk\nyJH6+OOPDawIRvH7/aqurlZGRobRpURUIBCQy+WSzWZTVlaWUlJSjC4pop588kk999xzslr77a+5\nHlksFk2bNk2TJk3S7373O6PLiagjR45o+PDhevTRRzVhwgR9//vf15kzZ4wuyxDbt2/XokWLjC4j\nohITE/XUU0/pG9/4hm666SbFx8dr2rRpRpfVK/rtbzNuqgtJam1t1dy5c1VYWKhhw4YZXU5EWa1W\n1dTU6Pjx49qzZ0+/ewZfT958802NGDFCbrd7wB5teuedd1RdXa233npLL7zwgt5++22jS4qYc+fO\nqaqqSo8//riqqqr0ta99TQUFBUaXFXHt7e164403NG/ePKNLiajDhw9r48aN8vv9+uSTT9Ta2qo/\n/elPRpfVK/ptaEtOTtaxY8dCr48dO6aRI0caWBEi7YsvvtB3vvMdPfLII5ozZ47R5RgmLi5Os2bN\n0vvvv290KRGzd+9elZSUyOFwaOHChfrnP/+p3Nxco8uKqBtvvFGSNHz4cD344IPat2+fwRVFzsiR\nIzVy5EhNnjxZkjR37lxVVVUZXFXkvfXWW5o4caKGDx9udCkR9f777+ub3/ymrr/+ekVFRemhhx7S\n3r17jS6rV/Tb0DZp0iR99NFH8vv9am9v16uvvqrZs2cbXRYiJBgMatmyZUpJSdETTzxhdDkR19jY\nqKamJklSW1ubvF6v3G63wVVFTn5+vo4dO6YjR45o+/btuueee/SHP/zB6LIi5syZM2ppaZEknT59\nWn//+98H1JXkSUlJuvnmm3Xo0CFJ5/+uKzU11eCqIm/btm1auHCh0WVE3NixY/Xuu++qra1NwWBQ\npaWl/ebPQ/rtY6yioqK0adMmzZgxQx0dHVq2bNmAunJw4cKF2r17t06cOKGbb75Zv/zlL/Xoo48a\nXVbEvPPOO3rllVdCtzyQzl8Cfu+99xpcWWTU1dXJ4/EoEAgoEAhoyZIlys7ONroswwy0P5doaGjQ\ngw8+KOn8qcLFixdr+vTpBlcVWc8//7wWL16s9vZ2jR49Wlu2bDG6pIg6ffq0SktLB9zfM0pSenq6\ncnNzNWnSJFmtVk2YMEE/+MEPjC6rV/TbW34AAAD0J/329CgAAEB/QmgDAAAwAUIbAACACRDaAAAA\nTIDQBgAAYAKENgAAABMgtAFAhJSXl+v+++83ugwAJkVoAwAAMAFCG4AB4/Tp05o1a5ZcLpecTqd2\n7NihyspKZWZmatKkSbr33ntVX18vSfrvf/+radOmyeVyaeLEiTpy5Igk6emnn5bT6VRaWpp27Ngh\n6fwRtMzMTM2bN0/jxo3TI488Etrnrl27NG7cOE2cOFF/+ctfQst3794tt9stt9utCRMmqLW1NYIz\nAcCM+u1jrADgQrt27VJycrJ27twpSWpubtbMmTNVUlKi66+/Xq+++qpWr16tzZs3a/HixXrmmWf0\nwAMPqL29XR0dHXrttdfk8/m0f/9+ffrpp5o8ebK+9a1vSZJqamr0wQcf6MYbb9SUKVO0d+/e0ONz\nysrKNHr0aC1YsCD0SK1f/epXevHFF3XnnXfqzJkzGjp0qGHzAsAcONIGYMBIS0uT1+vVT37yE1VU\nVOjo0aP6z3/+o2nTpsntduv//u//9PHHH6u1tVWffPKJHnjgAUnSkCFDdN111+mdd97RokWLZLFY\nNGLECN1999167733ZLFYdMcdd+imm26SxWKRy+XSkSNHdPDgQTkcDo0ePVqS9Mgjj+jLJwdOmTJF\nTz75pJ5//nl99tlnGjRokGHzAsAcONIGYMC47bbbVF1drZ07d+qnP/2psrKylJqaqr1793Zar6Wl\npdsxLnxc85dHzr56pGzQoEE6d+5clwfVf3XbH//4x/r2t7+tnTt3asqUKfrb3/6m22+//ap7A9D/\ncaQNwIBRV1en6OhoLV68WCtWrNC+ffvU2Niod999V5L0xRdf6IMPPlBsbKxGjhyp119/XZL0+eef\nq62tTXfddZdeffVVBQIBffrpp9qzZ4/uuOOOLkFOOh/mxo4dK7/fr//973+SpG3btoXeP3z4sFJT\nU7Vy5UpNnjxZH374YQRmAICZcaQNwIBx4MABPf3007JarRoyZIh++9vfatCgQVq+fLlOnTqlc+fO\n6cknn1RKSor++Mc/6oc//KF+9rOfafDgwfrzn/+sBx98UP/617+Unp4ui8Wi5557TiNGjFBtbW2X\no2rS+aNvL7/8smbNmqWYmBjdddddOn36tCSpsLBQZWVlslqtGj9+vGbOnBnp6QBgMpbgxf4XEQAA\nAH0Kp0cBAABMgNAGAABgAoQ2AAAAEyC0AQAAmAChDQAAwAQIbQAAACZAaAMAADCB/wdFi5aGCDL2\nJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7ae440f210>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "result = dummy_task.map(delay_times)\n",
    "visualize_tasks(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a nice and easy parallelization! We can see that we utilize all four engines quite well.\n",
    "\n",
    "But one short coming so far is that the tasks are not load balanced, so one engine might be idle while others still have more tasks to work on.\n",
    "\n",
    "However, the IPython parallel environment provides a number of alternative \"views\" of the engine cluster, and there is a view that provides load balancing as well (above we have used the \"direct view\", which is why we called it \"dview\").\n",
    "\n",
    "To obtain a load balanced view we simply use the `load_balanced_view` method in the engine cluster client instance `cli`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "lbview = cli.load_balanced_view()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "@lbview.parallel(block=True)\n",
    "def dummy_task_load_balanced(delay):\n",
    "    \"\"\" a dummy task that takes 'delay' seconds to finish \"\"\"\n",
    "    import os, time\n",
    "\n",
    "    t0 = time.time()\n",
    "    pid = os.getpid()\n",
    "    time.sleep(delay)\n",
    "    t1 = time.time()\n",
    "    \n",
    "    return [pid, t0, t1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAADSCAYAAAAYCPc3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHYtJREFUeJzt3X9wFPX9x/HXHQFiSoYEhYsNjnegCAmX3PHDjIPW5BsC\nUipiCygiSYX+sM7A6BRpkU5bOzSEkVJT0TrtUJLWGZCp0xil0l6mCRipxebHQccglnIFNMkYMCSB\npJHcff9wuGnIjwaSvb1Nno+ZzHh7++P9/tzmeLmb3bWFQqGQAAAAENXsZhcAAACA/43QBgAAYAGE\nNgAAAAsgtAEAAFgAoQ0AAMACYswuwGg2m83sEgAAAAasrxt7jIgjbaFQiJ8B/vzoRz8yvQYr/TBe\njBfjFT0/jBfjNRzGqz8jIrQBAABYHaENAADAAght6CYzM9PsEiyF8bo2jNe1YbyuDeN1bRivaxMN\n42UL/a8TqBZns9n+5zliAACAaNBfbuFIGwAAgAUQ2gAAACyA0AYAAGABhDYAAAALILQBAABYwLB/\njJUkPXngSTV3NBuy7oTYBEnqtn5/g1/pSenXvK7rXc7Kehu/4Wagn+tIGIuBMnMshsPnMJAeIvF9\nM5htjMTvw2js+VprGuzvTzSOQTQZEaGtuaNZzgSnIesONAckqdv6K09XXtf2rnc5K+tt/IabgX6u\nI2EsBsrMsRgOn8NAeojE981gtjESvw+jsedrrWmwvz/ROAbRhNOjAAAAFkBoAwAAsABCGwAAgAUQ\n2gAAACzAsNDW0dGhjIwMeTwepaSkaNOmTZKk8+fPKycnR9OmTdOCBQvU3Nwcnp6VlaX4+HitW7eu\n27p2794tt9ut9PR0LVq0SOfOnZMkFRUVaeLEifJ6vfJ6vfrNb35jVDsAAACmMiy0xcbGqry8XLW1\ntTp69KjKy8tVWVmpgoIC5eTk6MSJE8rOzlZBQUF4/i1btmj79u3d1tPZ2akNGzbo4MGD8vv9SktL\n086dOyV9/lDVlStXqqamRjU1NVqzZo1R7QAAAJjK0NOjcXFxkj4PXl1dXUpMTFRpaany8vIkSXl5\neSopKQnPO2/ePI0dO7bbOmJiYpSYmKi2tjaFQiFduHBBycnJkqRQKKRQKGRkCwAAAFHB0NAWDAbl\n8XjkcDiUlZWl1NRUNTY2yuFwSJIcDocaGxu7LWOz2boXaLersLBQM2fOVHJyso4fPx4+omaz2fTa\na68pLS1Ny5cv19mzZ41sBwAAwDSGhja73a7a2lqdPXtWhw4dUnl5ebf3bTZbj5B2tZaWFq1fv15+\nv18ff/yx3G63tm7dKkm6//779e9//1tHjx5VTk5O+Aje1Wr31qqiqEIVRRUK1AaGpDcAAIDBCtQG\nwhmloqii33kj8kSE8ePHa/HixaqqqpLD4VBDQ4OSkpJUX1+vSZMm9btsXV2dXC6XXC6XJGn58uXa\ntm2bJGnChAnh+dauXauNGzf2ug7Pwx7usAwAAKKO0+OU0+MMvz5YfLDPeQ070tbU1BS+MrS9vV0+\nn09er1dLlixRcXGxJKm4uFhLly7tttzVf6M2ZcoUHT9+XE1NTZIkn8+nlJQUSVJ9fX14vtLS0vB0\nAACA4cawI2319fXKy8tTMBhUMBjU6tWrlZ2dLa/XqxUrVmjXrl1yOp3at29feBmn06nW1lZ1dnaq\npKREPp9P06dPV35+vrKysmS32+V0OlVUVCRJeuGFF1RaWqqYmBjdeOON4ekAAADDjWGhze12q7q6\nusf0CRMmqKysrNdlAoFAr9Nzc3OVm5vbY3p+fr7y8/MHVScAAIAV8EQEAAAACyC0AQAAWAChDQAA\nwAIIbQAAABZAaAMAALAAQhsAAIAFROSJCGZLiE1QoDlg2LoldVt//Jj469re9S5nZb2N33Az0M91\nJIzFQJk5FsPhcxhID5H4vhnMNkbi92E09nytNQ329ycaxyCa2EJXP4JgmLHZbD2esgAAABCN+sst\nnB4FAACwAEIbAACABRDaAAAALIDQBgAAYAGENgAAAAsgtAEAAFgAoQ0AAMACCG0AAAAWQGgDAACw\nAEIbAACABRDaAAAALIDQBgAAYAGENgAAAAsgtAEAAFgAoQ0AAMACCG0AAAAWQGgDAACwAEIbAACA\nBRDaAAAALCDG7AIi4ckDT6q5o9nsMqKSv8Gv9KR0s8swRUJsgiRZct8YyOc2Ej/bSPQ80P3GyvvX\n9RrqnqNlH45UHVe2Ey19G8Xf4JekqO4xITZBz9/3vNll9DAiQltzR7OcCU6zy4hKlacrR+zYBJoD\nkmTJ/gfyuY3EzzYSPQ90v7Hy/nW9hrrnaNmHI1XHle1ES99GqTxdKSm6fzeu7MvRhtOjAAAAFkBo\nAwAAsABCGwAAgAUQ2gAAACzAsNDW0dGhjIwMeTwepaSkaNOmTZKk8+fPKycnR9OmTdOCBQvU3Nwc\nnp6VlaX4+HitW7eu27p2794tt9ut9PR0LVq0SOfOnZMkvfzyy0pLS5PX69Vdd90lv99vVDsAAACm\nMiy0xcbGqry8XLW1tTp69KjKy8tVWVmpgoIC5eTk6MSJE8rOzlZBQUF4/i1btmj79u3d1tPZ2akN\nGzbo4MGD8vv9SktL086dOyVJq1at0tGjR1VTU6NnnnlG3/3ud41qBwAAwFSGnh6Ni4uT9Hnw6urq\nUmJiokpLS5WXlydJysvLU0lJSXjeefPmaezYsd3WERMTo8TERLW1tSkUCqmlpUXJycmSpPj4+PB8\nbW1tuummm4xsBwAAwDSG3qctGAxq1qxZOnnypL7zne8oNTVVjY2NcjgckiSHw6HGxsZuy9hstm6v\n7Xa7CgsLNXPmTI0bN07Tpk3Tiy++GH7/pZde0o4dO3Tx4kUdPnzYyHYAAABMY2hos9vtqq2t1YUL\nF7Rw4UKVl5d3e99ms/UIaVdraWnR+vXr5ff75XK5tG7dOm3dulWbN2+WJD3xxBN64okntGfPHq1Z\ns6bHNiSpdm+tArEBSZLT45TT4xyS/gAAAAajoqJCFRUVA5o3Ik9EGD9+vBYvXqyqqio5HA41NDQo\nKSlJ9fX1mjRpUr/L1tXVyeVyyeVySZKWL1+ubdu29ZjvoYce0uOPP97rOjwPe6L6zssAAGBkyszM\nVGZmZvj1s88+2+e8hv1NW1NTU/jK0Pb2dvl8Pnm9Xi1ZskTFxcWSpOLiYi1durTbcqFQqNvrKVOm\n6Pjx42pqapIk+Xw+paSkSJI+/PDD8Hz79+9XWlqaUe0AAACYyrAjbfX19crLy1MwGFQwGNTq1auV\nnZ0tr9erFStWaNeuXXI6ndq3b194GafTqdbWVnV2dqqkpEQ+n0/Tp09Xfn6+srKyZLfb5XQ6VVRU\nJEl68cUXVVZWptGjR2vixInavXu3Ue0AAACYyrDQ5na7VV1d3WP6hAkTVFZW1usygUCg1+m5ubnK\nzc3tMf35558fVI0AAABWwRMRAAAALIDQBgAAYAGENgAAAAsgtAEAAFgAoQ0AAMACCG0AAAAWEJEn\nIpgtITZBgeaA2WVEpfgx8SN2bBJiEyTJkv0P5HMbiZ9tJHoe6H5j5f3reg11z9GyD0eqjivbiZa+\njRI/Jl5SdP9uXNmXo40tdPUjCIYZm83W4ykLAAAA0ai/3MLpUQAAAAsgtAEAAFhAv6GtqKhIs2bN\nUlxcnOLi4jRnzpzww94BAAAQOX1eiFBcXKzCwkLt2LFDXq9XoVBINTU1evrpp2Wz2Xp9FigAAACM\n0eeFCBkZGdq7d69cLle36YFAQA899JD+9re/RaTAweJCBAAAYBXXdSFCa2trj8AmSU6nU62trUNX\nHQAAAP6nPkNbbGxsnwv19x4AAACGXp+nR2+44QbddtttvS508uRJXbp0ydDChgqnRwEAgFX0l1v6\nvBChrq7OsIIAAABwbXgiAgAAQJS4riNt48aNk81m63OFLS0tQ1MdAAAA/ieOtAEAAESJ6zrS1t7e\nrpdfflknT56U2+3W2rVrFRPT5+wAAAAwUJ9H2lasWKExY8bo7rvv1ltvvSWn06nCwsJI1zdoHGkD\nAABW0V9u6TO0ud1uHTt2TJJ0+fJlzZ07VzU1NcZVaRBCGwAAsIrreiLCf58K5bQoAACAufo80jZq\n1CjFxcWFX7e3t+uGG274fCELXT3KkTYAAGAV13UhQldXl2EFAQAA4Nr0eXoUAAAA0WNE/LHakwee\nVHNHs+Hb8Tf4lZ6Ubvh2IsWq/fRWt1V7McJQjkU0jGtCbIIkDcnveDT0c8VgajG6j8Gufyg/MyNE\n035gtJHUa28SYhP0/H3Pm13GgI2I0Nbc0SxngtPw7VSerozIdiLFqv30VrdVezHCUI5FNIxroDkg\nSUNSRzT0c8VgajG6j8Gufyg/MyNE035gtJHUa2+u7ItWwelRAAAACyC0AQAAWAChDQAAwAIIbQAA\nABZgWGjr6OhQRkaGPB6PUlJStGnTJknS+fPnlZOTo2nTpmnBggVqbm4OT8/KylJ8fLzWrVvXbV27\nd++W2+1Wenq6Fi1apHPnzkmSduzYodTUVKWnp2v+/Pk6ffq0Ue0AAACYyrDQFhsbq/LyctXW1uro\n0aMqLy9XZWWlCgoKlJOToxMnTig7O1sFBQXh+bds2aLt27d3W09nZ6c2bNiggwcPyu/3Ky0tTTt3\n7pQkzZo1S1VVVfL7/Vq2bJk2btxoVDsAAACmMvT06JXHYHV2dqqrq0uJiYkqLS1VXl6eJCkvL08l\nJSXheefNm6exY8d2W0dMTIwSExPV1tamUCiklpYWJScnS5IyMzMVGxsrScrIyNDZs2eNbAcAAMA0\nhoa2YDAoj8cjh8OhrKwspaamqrGxUQ6HQ5LkcDjU2NjYbRmbzda9QLtdhYWFmjlzppKTk1VXV6c1\na9b02NauXbv05S9/2bhmAAAATGTozXXtdrtqa2t14cIFLVy4UOXl5d3et9lsPULa1VpaWrR+/Xr5\n/X65XC6tW7dOW7du1ebNm8PzvPLKK6qurtbPf/7zXtdRu7dWgdiAJMnpccrpcQ6qLwAAgKFQUVGh\nioqKAc0bkScijB8/XosXL1ZVVZUcDocaGhqUlJSk+vp6TZo0qd9l6+rq5HK55HK5JEnLly/Xtm3b\nwu+XlZUpPz9fhw4d0ujRo3tdh+dhz4i+4zMAAIhOmZmZyszMDL9+9tln+5zXsNOjTU1N4StD29vb\n5fP55PV6tWTJEhUXF0uSiouLtXTp0m7LhUKhbq+nTJmi48ePq6mpSZLk8/mUkpIiSaqpqdHjjz+u\nN954QzfddJNRrQAAAJjOsCNt9fX1ysvLUzAYVDAY1OrVq5WdnS2v16sVK1Zo165dcjqd2rdvX3gZ\np9Op1tZWdXZ2qqSkRD6fT9OnT1d+fr6ysrJkt9vldDpVVFQkSdq4caMuXryoZcuWSZJuvfXW8IUN\nAAAAw4lhoc3tdqu6urrH9AkTJqisrKzXZQKBQK/Tc3NzlZub22O6z+cbVI0AAABWwRMRAAAALIDQ\nBgAAYAGENgAAAAsgtAEAAFgAoQ0AAMACInJzXbMlxCYo0BwwfDvxY+Ijsp1IsWo/vdVt1V6MMJRj\nEQ3jmhCbIElDUkc09HPFYGoxuo/Brn8oPzMjRNN+YLSR1GtvruyLVmELXX0322HGZrP1uGEvAABA\nNOovt3B6FAAAwAIIbQAAABZAaAMAALAAQhsAAIAFENoAAAAsgNAGAABgAYQ2AAAACyC0AQAAWACh\nDQAAwAIIbQAAABZAaAMAALAAQhsAAIAFENoAAAAsgNAGAABgAYQ2AAAACyC0AQAAWAChDQAAwAII\nbQAAABZAaAMAALCAGLMLiISvl3zd7BIswd/gV3pSuuW2F+m6I2kgvQ3X/q+3L6uOx1DWnRCbIElq\n7mg2ZP1G8zf4Janfenvr0QiR2o5VmLEfWWnf7UtCbIKev+/5Qa9nRIQ2Z4LT7BIsofJ0ZUTHaqi2\nF+m6I2kgvQ3X/q+3L6uOx1DWHWgOSOr+3Welcak8XSmp/+/u3no0QqS2YxVm7EdW2nf7cmU/GixO\njwIAAFgAoQ0AAMACCG0AAAAWQGgDAACwAEIbAACABRgW2jo6OpSRkSGPx6OUlBRt2rRJknT+/Hnl\n5ORo2rRpWrBggZqbm8PTs7KyFB8fr3Xr1nVb1+7du+V2u5Wenq5Fixbp3LlzkqRDhw5p1qxZGj16\ntF577TWjWgEAADCdYaEtNjZW5eXlqq2t1dGjR1VeXq7KykoVFBQoJydHJ06cUHZ2tgoKCsLzb9my\nRdu3b++2ns7OTm3YsEEHDx6U3+9XWlqadu7cKUm69dZbVVxcrEceecSoNgAAAKKCoadH4+LiJH0e\nvLq6upSYmKjS0lLl5eVJkvLy8lRSUhKed968eRo7dmy3dcTExCgxMVFtbW0KhUJqaWlRcnKypM9D\nm9vtlt3OWV4AADC8GZp2gsGgPB6PHA6HsrKylJqaqsbGRjkcDkmSw+FQY2Njt2VsNlv3Au12FRYW\naubMmUpOTlZdXZ3WrFljZNkAAABRx9AnItjtdtXW1urChQtauHChysvLu71vs9l6hLSrtbS0aP36\n9fL7/XK5XFq3bp22bt2qzZs3D7iOiqKK8H87PU45Pc5raQMAAMAQFRUVqqioGNC8EXmM1fjx47V4\n8WJVVVXJ4XCooaFBSUlJqq+v16RJk/pdtq6uTi6XSy6XS5K0fPlybdu2rcd8/YW/zK9nDqp+AAAA\nI2RmZiozMzP8+tlnn+1zXsNOjzY1NYWvDG1vb5fP55PX69WSJUtUXFwsSSouLtbSpUu7LRcKhbq9\nnjJlio4fP66mpiZJks/nU0pKSo9lrl4OAABgODHsSFt9fb3y8vIUDAYVDAa1evVqZWdny+v1asWK\nFdq1a5ecTqf27dsXXsbpdKq1tVWdnZ0qKSmRz+fT9OnTlZ+fr6ysLNntdjmdThUVFUmS3nvvPX31\nq1/Vp59+qjfffFM//vGPdezYMaNaAgAAMI1hoc3tdqu6urrH9AkTJqisrKzXZQKBQK/Tc3NzlZub\n22P63LlzdebMmUHVCQAAYAXcKwMAAMACCG0AAAAWQGgDAACwAEIbAACABRDaAAAALCAiN9c1W6A5\nYHYJlhA/Jj6iYzVU24t03ZE0kN6Ga//X25dVx2Mo606ITZDU/bvPSuMSPyZeUv/f3b31aIRIbccq\nzNiPrLTv9uXKfjRYttAwvyutzWbjxrsAAMAS+sstnB4FAACwAEIbAACABRDaAAAALIDQBgAAYAGE\nNgAAAAsgtKGbiooKs0uwFMbr2jBe14bxujaM17VhvK5NNIwXoQ3dRMNOaSWM17VhvK4N43VtGK9r\nw3hdm2gYL0IbAACABRDaAAAALGBEPBEBAADAKvqKZsP+2aPDPJMCAIARgtOjAAAAFkBoAwAAsABC\nGwAAgAUM69B24MABTZ8+Xbfffru2bdtmdjlRbc2aNXI4HHK73WaXYglnzpxRVlaWUlNTNXPmTP3i\nF78wu6So1tHRoYyMDHk8HqWkpGjTpk1ml2QJXV1d8nq9uv/++80uJeo5nU6lpaXJ6/XqzjvvNLuc\nqNfc3Kxly5ZpxowZSklJ0bvvvmt2SVHrgw8+kNfrDf+MHz/etO/8YXv1aFdXl+644w6VlZUpOTlZ\nc+fO1Z49ezRjxgyzS4tKb7/9tsaNG6fc3FwdO3bM7HKiXkNDgxoaGuTxeNTW1qbZs2erpKSE/asf\nly5dUlxcnC5fvqy7775b27dv19133212WVFtx44dqqqqUmtrq0pLS80uJ6q5XC5VVVVpwoQJZpdi\nCXl5ebr33nu1Zs0aXb58WRcvXtT48ePNLivqBYNBJScn68iRI7rlllsivv1he6TtyJEjuu222+R0\nOjV69Gg9/PDDev31180uK2rdc889SkxMNLsMy0hKSpLH45EkjRs3TjNmzNDHH39sclXRLS4uTpLU\n2dmprq4u/nH9H86ePas//vGP+sY3vsFV8APEOA3MhQsX9Pbbb2vNmjWSpJiYGALbAJWVlWnq1Kmm\nBDZpGIe2jz76qNugTp48WR999JGJFWG4CgQCqqmpUUZGhtmlRLVgMCiPxyOHw6GsrCylpKSYXVJU\ne+qpp/Tcc8/Jbh+2X9NDymazaf78+ZozZ45+/etfm11OVDt16pQmTpyoxx57TLNmzdI3v/lNXbp0\nyeyyLGHv3r165JFHTNv+sP024Ka6iIS2tjYtW7ZMhYWFGjdunNnlRDW73a7a2lqdPXtWhw4diorn\n+EWrN998U5MmTZLX6+Xo0QC98847qqmp0VtvvaUXX3xRb7/9ttklRa3Lly+rurpaTzzxhKqrq/WF\nL3xBBQUFZpcV9To7O/XGG29o+fLlptUwbENbcnKyzpw5E3595swZTZ482cSKMNx89tln+trXvqZH\nH31US5cuNbscyxg/frwWL16sv//972aXErUOHz6s0tJSuVwurVy5Un/5y1+Um5trdllR7eabb5Yk\nTZw4UQ8++KCOHDlickXRa/LkyZo8ebLmzp0rSVq2bJmqq6tNrir6vfXWW5o9e7YmTpxoWg3DNrTN\nmTNHH374oQKBgDo7O/Xqq69qyZIlZpeFYSIUCmnt2rVKSUnRk08+aXY5Ua+pqUnNzc2SpPb2dvl8\nPnm9XpOril75+fk6c+aMTp06pb179+r//u//9Nvf/tbssqLWpUuX1NraKkm6ePGi/vznP3MlfD+S\nkpJ0yy236MSJE5I+/zut1NRUk6uKfnv27NHKlStNrWHYPsYqJiZGO3fu1MKFC9XV1aW1a9dyZV8/\nVq5cqYMHD+rcuXO65ZZb9JOf/ESPPfaY2WVFrXfeeUevvPJK+BYDkrR161bdd999JlcWnerr65WX\nl6dgMKhgMKjVq1crOzvb7LIsgz/36F9jY6MefPBBSZ+f+lu1apUWLFhgclXR7YUXXtCqVavU2dmp\nqVOnavfu3WaXFNUuXryosrIy0/9ectje8gMAAGA4GbanRwEAAIYTQhsAAIAFENoAAAAsgNAGAABg\nAYQ2AAAACyC0AQAAWAChDQAipKKiQvfff7/ZZQCwKEIbAACABRDaAIwYFy9e1OLFi+XxeOR2u7Vv\n3z5VVVUpMzNTc+bM0X333aeGhgZJ0j//+U/Nnz9fHo9Hs2fP1qlTpyRJTz/9tNxut9LS0rRv3z5J\nnx9By8zM1PLlyzVjxgw9+uij4W0eOHBAM2bM0OzZs/WHP/whPP3gwYPyer3yer2aNWuW2traIjgS\nAKxo2D7GCgCuduDAASUnJ2v//v2SpJaWFi1atEilpaW68cYb9eqrr2rz5s3atWuXVq1apWeeeUYP\nPPCAOjs71dXVpddee01+v19Hjx7VJ598orlz5+pLX/qSJKm2tlbvv/++br75Zs2bN0+HDx/WrFmz\n9K1vfUvl5eWaOnWqHnroofAjqX72s5/ppZde0l133aVLly5p7Nixpo0LAGvgSBuAESMtLU0+n0/f\n//73VVlZqdOnT+sf//iH5s+fL6/Xq5/+9Kf66KOP1NbWpo8//lgPPPCAJGnMmDG64YYb9M477+iR\nRx6RzWbTpEmTdO+99+q9996TzWbTnXfeqS9+8Yuy2WzyeDw6deqUjh8/LpfLpalTp0qSHn30UV15\ncuC8efP01FNP6YUXXtCnn36qUaNGmTYuAKyBI20ARozbb79dNTU12r9/v37wgx8oKytLqampOnz4\ncLf5Wltb+1zH1Y9rvnLk7L+PlI0aNUqXL1/u8aD3/172e9/7nr7yla9o//79mjdvnv70pz/pjjvu\nuO7eAAx/HGkDMGLU19crNjZWq1at0oYNG3TkyBE1NTXp3XfflSR99tlnev/99xUfH6/Jkyfr9ddf\nlyT95z//UXt7u+655x69+uqrCgaD+uSTT3To0CHdeeedPYKc9HmYmz59ugKBgP71r39Jkvbs2RN+\n/+TJk0pNTdXGjRs1d+5cffDBBxEYAQBWxpE2ACPGsWPH9PTTT8tut2vMmDH65S9/qVGjRmn9+vW6\ncOGCLl++rKeeekopKSn63e9+p29/+9v64Q9/qNGjR+v3v/+9HnzwQf31r39Venq6bDabnnvuOU2a\nNEl1dXU9jqpJnx99+9WvfqXFixcrLi5O99xzjy5evChJKiwsVHl5uex2u2bOnKlFixZFejgAWIwt\n1Nv/IgIAACCqcHoUAADAAghtAAAAFkBoAwAAsABCGwAAgAUQ2gAAACyA0AYAAGABhDYAAAAL+H/I\n32FnzmeQowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7ae43ec990>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "result = dummy_task_load_balanced.map(delay_times)\n",
    "visualize_tasks(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example above we can see that the engine cluster is a bit more efficiently used, and the time to completion is shorter than in the previous example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many other ways to use the IPython parallel environment. The official documentation has a nice guide:\n",
    "\n",
    "* http://ipython.org/ipython-doc/dev/parallel/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When more communication between processes is required, sophisticated solutions such as MPI and OpenMP are often needed. MPI is process based parallel processing library/protocol, and can be used in Python programs through the `mpi4py` package:\n",
    "\n",
    "http://mpi4py.scipy.org/\n",
    "\n",
    "To use the `mpi4py` package we include `MPI` from `mpi4py`:\n",
    "\n",
    "    from mpi4py import MPI\n",
    "\n",
    "A MPI python program must be started using the `mpirun -n N` command, where `N` is the number of processes that should be included in the process group.\n",
    "\n",
    "Note that the IPython parallel enviroment also has support for MPI, but to begin with we will use `mpi4py` and the `mpirun` in the follow examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting mpitest.py\n"
     ]
    }
   ],
   "source": [
    "%%file mpitest.py\n",
    "\n",
    "from mpi4py import MPI\n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "rank = comm.Get_rank()\n",
    "\n",
    "if rank == 0:\n",
    "   data = [1.0, 2.0, 3.0, 4.0]\n",
    "   comm.send(data, dest=1, tag=11)\n",
    "elif rank == 1:\n",
    "   data = comm.recv(source=0, tag=11)\n",
    "    \n",
    "print(\"rank =\", rank, \", data =\", data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank = 0 , data = [1.0, 2.0, 3.0, 4.0]\r\n",
      "rank = 1 , data = [1.0, 2.0, 3.0, 4.0]\r\n"
     ]
    }
   ],
   "source": [
    "!mpirun -n 2 python mpitest.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Send a numpy array from one process to another:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting mpi-numpy-array.py\n"
     ]
    }
   ],
   "source": [
    "%%file mpi-numpy-array.py\n",
    "\n",
    "from mpi4py import MPI\n",
    "import numpy\n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "rank = comm.Get_rank()\n",
    "\n",
    "if rank == 0:\n",
    "   data = numpy.random.rand(10)\n",
    "   comm.Send(data, dest=1, tag=13)\n",
    "elif rank == 1:\n",
    "   data = numpy.empty(10, dtype=numpy.float64)\n",
    "   comm.Recv(data, source=0, tag=13)\n",
    "    \n",
    "print(\"rank =\", rank, \", data =\", data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank = 0 , data = [ 0.71397658  0.37182268  0.25863587  0.08007216  0.50832534  0.80038331\r\n",
      "  0.90613024  0.99535428  0.11717776  0.48353805]\r\n",
      "rank = 1 , data = [ 0.71397658  0.37182268  0.25863587  0.08007216  0.50832534  0.80038331\r\n",
      "  0.90613024  0.99535428  0.11717776  0.48353805]\r\n"
     ]
    }
   ],
   "source": [
    "!mpirun -n 2 python mpi-numpy-array.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 3: Matrix-vector multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare some random data\n",
    "N = 16\n",
    "A = numpy.random.rand(N, N)\n",
    "numpy.save(\"random-matrix.npy\", A)\n",
    "x = numpy.random.rand(N)\n",
    "numpy.save(\"random-vector.npy\", x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting mpi-matrix-vector.py\n"
     ]
    }
   ],
   "source": [
    "%%file mpi-matrix-vector.py\n",
    "\n",
    "from mpi4py import MPI\n",
    "import numpy\n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "rank = comm.Get_rank()\n",
    "p = comm.Get_size()\n",
    "\n",
    "def matvec(comm, A, x):\n",
    "    m = A.shape[0] / p\n",
    "    y_part = numpy.dot(A[rank * m:(rank+1)*m], x)\n",
    "    y = numpy.zeros_like(x)\n",
    "    comm.Allgather([y_part,  MPI.DOUBLE], [y, MPI.DOUBLE])\n",
    "    return y\n",
    "\n",
    "A = numpy.load(\"random-matrix.npy\")\n",
    "x = numpy.load(\"random-vector.npy\")\n",
    "y_mpi = matvec(comm, A, x)\n",
    "\n",
    "if rank == 0:\n",
    "    y = numpy.dot(A, x)\n",
    "    print(y_mpi)\n",
    "    print(\"sum(y - y_mpi) =\", (y - y_mpi).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 6.40342716  3.62421625  3.42334637  3.99854639  4.95852419  6.13378754\r\n",
      "  5.33319708  5.42803442  5.12403754  4.87891654  2.38660728  6.72030412\r\n",
      "  4.05218475  3.37415974  3.90903001  5.82330226]\r\n",
      "sum(y - y_mpi) = 0.0\r\n"
     ]
    }
   ],
   "source": [
    "!mpirun -n 4 python mpi-matrix-vector.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 4: Sum of the elements in a vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare some random data\n",
    "N = 128\n",
    "a = numpy.random.rand(N)\n",
    "numpy.save(\"random-vector.npy\", a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting mpi-psum.py\n"
     ]
    }
   ],
   "source": [
    "%%file mpi-psum.py\n",
    "\n",
    "from mpi4py import MPI\n",
    "import numpy as np\n",
    "\n",
    "def psum(a):\n",
    "    r = MPI.COMM_WORLD.Get_rank()\n",
    "    size = MPI.COMM_WORLD.Get_size()\n",
    "    m = len(a) / size\n",
    "    locsum = np.sum(a[r*m:(r+1)*m])\n",
    "    rcvBuf = np.array(0.0, 'd')\n",
    "    MPI.COMM_WORLD.Allreduce([locsum, MPI.DOUBLE], [rcvBuf, MPI.DOUBLE], op=MPI.SUM)\n",
    "    return rcvBuf\n",
    "\n",
    "a = np.load(\"random-vector.npy\")\n",
    "s = psum(a)\n",
    "\n",
    "if MPI.COMM_WORLD.Get_rank() == 0:\n",
    "    print(\"sum =\", s, \", numpy sum =\", a.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sum = 64.948311241 , numpy sum = 64.948311241\r\n"
     ]
    }
   ],
   "source": [
    "!mpirun -n 4 python mpi-psum.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* http://mpi4py.scipy.org\n",
    "\n",
    "* http://mpi4py.scipy.org/docs/usrman/tutorial.html\n",
    "\n",
    "* https://computing.llnl.gov/tutorials/mpi/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OpenMP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about OpenMP? OpenMP is a standard and widely used thread-based parallel API that unfortunaltely is **not** useful directly in Python. The reason is that the CPython implementation use a global interpreter lock, making it impossible to simultaneously run several Python threads. Threads are therefore not useful for parallel computing in Python, unless it is only used to wrap compiled code that do the OpenMP parallelization (Numpy can do something like that). \n",
    "\n",
    "This is clearly a limitation in the Python interpreter, and as a consequence all parallelization in Python must use processes (not threads).\n",
    "\n",
    "However, there is a way around this that is not that painful. When calling out to compiled code the GIL is released, and it is possible to write Python-like code in Cython where we can selectively release the GIL and do OpenMP computations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This system has 12 cores\n"
     ]
    }
   ],
   "source": [
    "N_core = multiprocessing.cpu_count()\n",
    "\n",
    "print((\"This system has %d cores\" % N_core))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a simple example that shows how OpenMP can be used via cython:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext cythonmagic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -f -c-fopenmp --link-args=-fopenmp -c-g\n",
    "\n",
    "cimport cython\n",
    "cimport numpy\n",
    "from cython.parallel import prange, parallel\n",
    "cimport openmp\n",
    "\n",
    "def cy_openmp_test():\n",
    "\n",
    "    cdef int n, N\n",
    "\n",
    "    # release GIL so that we can use OpenMP\n",
    "    with nogil, parallel():\n",
    "        N = openmp.omp_get_num_threads()\n",
    "        n = openmp.omp_get_thread_num()\n",
    "        with gil:\n",
    "            print(\"Number of threads %d: thread number %d\" % (N, n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of threads 12: thread number 0\n",
      "Number of threads 12: thread number 10\n",
      "Number of threads 12: thread number 8\n",
      "Number of threads 12: thread number 4\n",
      "Number of threads 12: thread number 7\n",
      "Number of threads 12: thread number 3\n",
      "Number of threads 12: thread number 2\n",
      "Number of threads 12: thread number 1\n",
      "Number of threads 12: thread number 11\n",
      "Number of threads 12: thread number 9\n",
      "Number of threads 12: thread number 5\n",
      "Number of threads 12: thread number 6\n"
     ]
    }
   ],
   "source": [
    "cy_openmp_test()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: matrix vector multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare some random data\n",
    "N = 4 * N_core\n",
    "\n",
    "M = numpy.random.rand(N, N)\n",
    "x = numpy.random.rand(N)\n",
    "y = numpy.zeros_like(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first look at a simple implementation of matrix-vector multiplication in Cython:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "cimport cython\n",
    "cimport numpy\n",
    "import numpy\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.wraparound(False)\n",
    "def cy_matvec(numpy.ndarray[numpy.float64_t, ndim=2] M, \n",
    "              numpy.ndarray[numpy.float64_t, ndim=1] x, \n",
    "              numpy.ndarray[numpy.float64_t, ndim=1] y):\n",
    "\n",
    "    cdef int i, j, n = len(x)\n",
    "\n",
    "    for i from 0 <= i < n:\n",
    "        for j from 0 <= j < n:\n",
    "            y[i] += M[i, j] * x[j]\n",
    "            \n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check that we get the same results\n",
    "y = numpy.zeros_like(x)\n",
    "cy_matvec(M, x, y)\n",
    "numpy.dot(M, x) - y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000 loops, best of 3: 2.93 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit numpy.dot(M, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000 loops, best of 3: 5.4 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit cy_matvec(M, x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Cython implementation here is a bit slower than numpy.dot, but not by much, so if we can use multiple cores with OpenMP it should be possible to beat the performance of numpy.dot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -f -c-fopenmp --link-args=-fopenmp -c-g\n",
    "\n",
    "cimport cython\n",
    "cimport numpy\n",
    "from cython.parallel import parallel\n",
    "cimport openmp\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.wraparound(False)\n",
    "def cy_matvec_omp(numpy.ndarray[numpy.float64_t, ndim=2] M, \n",
    "                  numpy.ndarray[numpy.float64_t, ndim=1] x, \n",
    "                  numpy.ndarray[numpy.float64_t, ndim=1] y):\n",
    "\n",
    "    cdef int i, j, n = len(x), N, r, m\n",
    "\n",
    "    # release GIL, so that we can use OpenMP\n",
    "    with nogil, parallel():\n",
    "        N = openmp.omp_get_num_threads()\n",
    "        r = openmp.omp_get_thread_num()\n",
    "        m = n / N\n",
    "        \n",
    "        for i from 0 <= i < m:\n",
    "            for j from 0 <= j < n:\n",
    "                y[r * m + i] += M[r * m + i, j] * x[j]\n",
    "\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check that we get the same results\n",
    "y = numpy.zeros_like(x)\n",
    "cy_matvec_omp(M, x, y)\n",
    "numpy.dot(M, x) - y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000 loops, best of 3: 2.95 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit numpy.dot(M, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 209 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit cy_matvec_omp(M, x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, this implementation is much slower than numpy.dot for this problem size, because of overhead associated with OpenMP and threading, etc. But let's look at the how the different implementations compare with larger matrix sizes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_vec  = numpy.arange(25, 2000, 25) * N_core"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "duration_ref = numpy.zeros(len(N_vec))\n",
    "duration_cy = numpy.zeros(len(N_vec))\n",
    "duration_cy_omp = numpy.zeros(len(N_vec))\n",
    "\n",
    "for idx, N in enumerate(N_vec):\n",
    "    \n",
    "    M = numpy.random.rand(N, N)\n",
    "    x = numpy.random.rand(N)\n",
    "    y = numpy.zeros_like(x)\n",
    "    \n",
    "    t0 = time.time()\n",
    "    numpy.dot(M, x)\n",
    "    duration_ref[idx] = time.time() - t0\n",
    "    \n",
    "    t0 = time.time()\n",
    "    cy_matvec(M, x, y)\n",
    "    duration_cy[idx] = time.time() - t0\n",
    "    \n",
    "    t0 = time.time()\n",
    "    cy_matvec_omp(M, x, y)\n",
    "    duration_cy_omp[idx] = time.time() - t0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAGGCAYAAACnlj1bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8zWcbx/HPCUKoEVvs2kXVVpo2aiuqURpb0FZtumiN\nRFGKVrW0RYvSpmYRI7Q0qFF77xFSodQmIus8f/yepE2NRM755Zwk3/frlZecc36/+75PnqTPlTvX\nfV0Wq9VqRURERERETOPi6AWIiIiIiKR1CrpFREREREymoFtERERExGQKukVERERETKagW0RERETE\nZAq6RURERERMpqBbRERERMRkCrpFREREREzmtEH3mTNn6NmzJ23btnX0UkREREREbOK0QXfJkiWZ\nOXOmo5chIiIiImKzFA26u3fvToECBahcuXKC54OCgihfvjxlypRh/PjxKbkkERERERHTpWjQ7evr\nS1BQUILnYmJi6Nu3L0FBQRw+fJiAgACOHDmSkssSERERETFVigbdnp6euLu7J3hu+/btlC5dmhIl\nSpApUyZ8fHxYtmwZV69epVevXuzdu1e73yIiIiKSqmV09ALOnz9P0aJF4x8XKVKEP/74g9y5c/P1\n118/8l6LxWL28kREREREALBarcm+1+EHKW0NnK1Wqz7+/zFy5EiHr8FZ1mbGfPYa05ZxknPv49xj\n1rVp/cOZvxZp4WfPXuM688/e41zvzN9vKf3hzF8L/ezZZxxn+tmzlcOD7sKFCxMaGhr/ODQ0lCJF\nijhwRamXl5eXo5fwUCm9NjPms9eYtoyTnHsf5x5n/h5yZs78dUsLP3v2GteZf/aSO0d658xfM/3s\n2WectPSzZ7HaI3R/DCEhIbRs2ZIDBw4AEB0dTbly5Vi3bh0eHh7UqlWLgIAAKlSokOhYFouFkSNH\n4uXl5dQ/eCJpkZ+fH35+fo5ehki6o589kZQVHBxMcHAw/v7+Nu14p2jQ3b59ezZs2MCVK1fInz8/\no0aNwtfXl9WrVzNw4EBiYmLo0aMHQ4cOTdJ4FovFLtv9IvL4goOD9cuuiAPoZ0/EMWyNO1N8p9ue\nFHSLiIiISEqwNe50ePUSM6iqidhCv8iJiIiIvaX6oNvPz++BOd0KnCQ59AubiIiI/FtcTret0mR6\nidJOJLn0vSMiIiIPYmuM4PCSgSIiIiIiaZ2CbhERERERk6X6oNvPz88ueTYiIiIiIv8VHBxsl9r4\nyukW+Rd974iIiMiDKKdbRERERMTJKehOYSVKlGDSpElUqVKFXLly4ePjw71795g9ezaenp4JrnVx\nceH06dMAdOvWjd69e9O8eXOyZ8+Op6cnFy9eZMCAAbi7u1OhQgX27t2bYJ5x48ZRsWJFcufOTffu\n3bl37x4AlSpVYsWKFfHXRkVFkTdvXvbt25cCXwERERGR9CfVB92pLafbYrGwcOFC1qxZw5kzZ9i/\nfz+zZ89OUn3ohQsXMmbMGP7++29cXV2pU6cONWvW5OrVq7z66qsMHjw4wfU//vgja9eu5dSpUxw/\nfpzRo0cD0LVrV+bNmxd/3apVqyhcuDBVqlSx75sVERERSeXsldOdJoLu/zbGcXb9+/enYMGCuLu7\n07JlywQ71A9jsVjw9vamatWqZM6cmVdeeYVs2bLRqVMnLBYL7dq1Y8+ePQmu79u3L4ULF8bd3Z0P\nP/yQgIAAADp27MjKlSu5ffs2AHPnzqVz587mvFkRERGRVMzLy0tBty0sFvt8JEfBggXjP8+aNWt8\n8JuY/Pnzx3+eJUuWBI/d3NzuG6do0aLxnxcrVoywsDAAPDw8qFevHosWLeL69esEBQXRsWPHZL0X\nEREREUlcqm8Dn1zOVqAiW7ZshIeHxz++ePGizWOeO3cuweceHh7xj7t27cq3335LVFQUdevWpVCh\nQjbPJyIiIiIPlm53up1FXOmZKlWqcOjQIfbt20dERMR9f8Z43BI1VquVadOmcf78ea5evcqYMWPw\n8fGJf/2VV15h9+7dTJkyhS5dutj8PkRERETk4VJ90J3aDlL+l8ViwWKxUKZMGUaMGEHDhg0pV64c\nnp6eCQ5Xxl33sMdxz/378w4dOtC4cWNKlSpFmTJlGDZsWPzrWbJkwdvbm5CQELy9vU18hyIiIiKp\nl5rjoOY4j1KyZEm+/fZbXnzxxYde89FHH3HixAm+//77FFyZc9P3joiIiDyIrTFCus3pTu+uXr3K\nd999x9y5cx29FBEREZE0L9Wnl8jjmzFjBsWKFaNZs2Y899xzjl6OiIiISJqn9BKRf9H3joiIiDyI\nrTGCdrpFREREREymoFtERERExGSpPuhO7SUDRURERMR5qWQgyukW+9P3joiIiDyIcrpFRERERJyc\ngu5Uys/Pj86dOzt6GSIiIiKSBAq6U4Hg4GCKFi2a4Ln/toAXEREREeeloDuVUt6xiIiISOqhoNsB\nQkND8fb2Jn/+/OTNm5c+ffqQJ08eDh48GH/NpUuXyJYtG+fOnaNZs2aEhYWRPXt2cuTIwYULF7BY\nLERGRtK1a1dy5MhBpUqV2LVrV/z9R44cwcvLC3d3dypVqkRgYGD8a926daNPnz60aNGCHDlyUKdO\nHU6fPp2iXwMRERGR9ERBdwqLiYmhRYsWlCxZkrNnzxIWFkaHDh3w8fFh3rx58dcFBATQsGFDihUr\nRlBQEB4eHty6dYubN29SqFAhrFYry5cvp3379ty4cYNWrVrRt29fAKKiomjZsiVNmzbl8uXLfPHF\nF3Ts2JHjx4/Hjz9//nz8/Py4du0apUuX5sMPP0zxr4WIiIhIepHqg+7UVqd7+/btXLhwgQkTJuDm\n5oarqyv16tWjS5cuBAQExF83d+7c+IOSD0sl8fT0pGnTplgsFjp16sS+ffsA2LZtG3fu3GHIkCFk\nzJiR+vXr06JFiwTje3t7U6NGDTJkyEDHjh3Zu3evie9aREREJHWyV53ujLYvxbGS+0Ww+NvnIKJ1\n5OPlVoeGhlK8eHFcXBL+vlO7dm3c3NwIDg6mYMGCnDp1ilatWj1yrAIFCsR/njVrViIiIoiNjSUs\nLOy+g5fFixcnLCwMMA5h/vteNzc3bt++/VjvQ0RERCQ98PLywsvLC39/f5vGSfVBd3I9brBsL0WL\nFuXcuXPExMSQIUOGBK917dqVefPmUaBAAdq2bYurqyvw4Eolj6pe4uHhQWhoKFarNf66s2fPUr58\neTu+ExERERFJqlSfXpLa1K5dm0KFCjFkyBDCw8OJiIhgy5YtAHTq1IklS5bwww8/0KVLl/h7ChQo\nwJUrV7h582b8c4+qXlK7dm2yZs3KJ598QlRUFMHBwaxYsQIfH59E7xURERER+1PQncJcXFwIDAzk\n5MmTFCtWjKJFi7JgwQLA2AWvVq0aLi4uPPfcc/H3lC9fnvbt2/Pkk0+SO3fu+Ool/93tjnvs6upK\nYGAgq1evJl++fPTt25e5c+dStmzZ+Osedq+IiIiI2J/Fmoq3PS0WywN3bR/2fGrQo0cPChcuzKhR\noxy9lHQpNX/viIiIiHlsjRHSbU63MwoJCWHJkiWqJCIiIiKSxii9xEkMHz6cypUr895771G8eHFH\nL0dERERE7EjpJSL/ou8dEREReRBbYwTtdIuIiIiImExBt4iIiIiIyRR0i4iIiIiYLNUH3X5+fgQH\nBzt6GSIiIiKSBgUHB+Pn52fzODpIKfIv+t4RERGRB9FBynTKz8+Pzp07O3oZIiIiIpIECrpTgeDg\nYIoWLZrgObVtFxEREUk9FHSnUimVAtGtWzfmzJmTInOJiIiIpFUKuh0gNDQUb29v8ufPT968eenT\npw958uTh4MGD8ddcunSJbNmyce7cOZo1a0ZYWBjZs2cnR44cXLhwAYvFQmRkJF27diVHjhxUqlSJ\nXbt2xd9/5MgRvLy8cHd3p1KlSgQGBsa/1q1bN/r06UOLFi3IkSMHderU4fTp0w9c66N21Lds2ULN\nmjXJlSsXtWrVYuvWrfGveXl5MXToUGrXrk3OnDlp3bo1165di39927Zt1K1bF3d3d5555hk2bNiQ\n4N4RI0bw3HPPkSNHDpo0acKVK1cACAkJwcXFhdmzZ1OsWDHy5MnD119/zY4dO3j66adxd3enX79+\n8WPNnj2bevXq0a9fP3LlykWFChVYv359Uv5nEhEREbEbBd0pLCYmhhYtWlCyZEnOnj1LWFgYHTp0\nwMfHh3nz5sVfFxAQQMOGDSlWrBhBQUF4eHhw69Ytbt68SaFChbBarSxfvpz27dtz48YNWrVqRd++\nfQGIioqiZcuWNG3alMuXL/PFF1/QsWNHjh8/Hj/+/Pnz8fPz49q1a5QuXZoPP/zwoWt+UOB99epV\nXnrpJQYOHMjVq1cZPHgwL730UoLAeu7cucyaNYsLFy6QMWNG+vfvD8D58+dp0aIFI0aM4Nq1a0yc\nOJE2bdrEB9Zx73/27NlcunSJyMhIJk6cmGD+7du3c/LkSX766ScGDBjA2LFjWb9+PYcOHWLBggVs\n3LgxwbWlS5fmypUr+Pv74+3tnWCdIiIiImZT0J3Ctm/fzoULF5gwYQJubm64urpSr149unTpQkBA\nQPx1c+fOjT8o+bBUEk9PT5o2bYrFYqFTp07s27cPMHaR79y5w5AhQ8iYMSP169enRYsWCcb39vam\nRo0aZMiQgY4dO7J3794HzmG1Wh84/8qVKylXrhwdO3bExcUFHx8fypcvz/LlywEjUO/SpQtPPfUU\nWbNm5aOPPmLBggXExsYyb948mjdvTtOmTQFo2LAhNWrUYOXKlfH3+vr6Urp0abJkyUK7du3uW9/w\n4cNxdXWlUaNGZM+enQ4dOpA3b148PDzw9PRkz5498dfmz5+fAQMGkCFDBtq1a0e5cuXi5xIRERFJ\nCek36LZY7PPxmEJDQylevDguLgm/9LVr18bNzY3g4GCOHj3KqVOnaNWq1SPHKlCgQPznWbNmJSIi\ngtjYWMLCwu47eFm8eHHCwsL+/9YtCe51c3Pj9u3b8Y/j0jTc3d0JCAigd+/e8Y/jdtPDwsIoVqzY\nQ+cAEqyhWLFiREVF8ffff3P27FkWLlwYP6a7uzubN2/m4sWL8dcXLFjwoev773t3c3O77/GdO3fi\nHxcuXPiR6xQRERExW0ZHL8BhHFSLuWjRopw7d46YmBgyZMiQ4LWuXbsyb948ChQoQNu2bXF1dQUe\nnN7xqFxrDw8PQkNDsVqt8dedPXuW8uXLJ2mN+/fvj//c19eX+vXr06VLlwTXFC5cmCVLliR47uzZ\nszRr1iz+8blz5xJ8nilTJvLly0exYsXo3Lkz06dPT9J6bHX+/Pn71vnyyy+nyNwiIiIikJ53uh2k\ndu3aFCpUiCFDhhAeHk5ERARbtmwBoFOnTixZsoQffvghQZBboEABrly5ws2bN+Ofe1T1ktq1a5M1\na1Y++eQToqKiCA4OZsWKFfj4+CR674M86PrmzZtz/PhxAgICiI6OZv78+Rw9epQWLVrE3zNv3jyO\nHDlCeHg4I0aMoG3btvGpMIGBgaxdu5aYmBgiIiIIDg5OEBzbWp3l3/dfunSJKVOmEBUVxcKFCzl2\n7BjNmze3aXwRERGRx6GgO4W5uLgQGBjIyZMnKVasGEWLFmXBggWAsQterVo1XFxceO655+LvKV++\nPO3bt+fJJ58kd+7c8dVL/rvbHffY1dWVwMBAVq9eTb58+ejbty9z586lbNmy8dc97N4HedBruXPn\nZsWKFUyaNIm8efMyceJEVqxYQe7cuePv6dy5M926daNQoUJERkYyZcoUAIoUKcKyZcsYO3Ys+fPn\np1ixYkyaNClBoPzvOf+73qTUKP/3NbVr1+bEiRPky5eP4cOHs2jRItzd3RMdQ0RERMRe1AbeyfTo\n0YPChQszatQoRy/FJvXr16dz5850797doeuYPXs23377LZs2bUrS9an5e0dERETMY2uMkH5zup1Q\nSEgIS5YseWglkdRGwauIiIiIQeklTmL48OFUrlyZ9957j+LFizt6OXbhDK3qH5RKIyIiIpLSnDa9\n5M6dO/Tu3ZvMmTPj5eVFhw4d7rsmLaaXiGPpe0dEJH2IjIlk+q7pvFXjLTK4ZEj8Bkn3bI0RnHan\ne8mSJbRr147p06fHN1wRERERsdWpq6d4cc6LrDuzjnsx9xy9HEknUjTo7t69OwUKFKBy5coJng8K\nCqJ8+fKUKVOG8ePHA0Zt5bjmKv+tZy0iIiLyuO5G3eXdte9Se2ZtXi73MovbLSZrpqyOXpakEyka\ndPv6+hIUFJTguZiYGPr27UtQUBCHDx8mICCAI0eOUKRIEUJDQwGIjY1NyWWKiIhIGmK1Wll8eDGV\nvqpE6M1QDvc5zLv13sXF4rR/8Jc0KEWrl3h6ehISEpLgue3bt1O6dGlKlCgBgI+PD8uWLaN///70\n7duXlStXJtoOXUREROS/Lty6wO4Lu5m7fy4HLh3gmxbf0PDJho5elqRTDi8Z+O80EjAap/zxxx9k\nzZqV7777LtH7/fz84j/38vLCy8sLcI7KGSIiIuI47Ra1w4KFukXrMuvlWbhlcnP0kiQVCQ4OJjg4\n2G7jOTzotjU4/nfQHUfVJ0RERNK3S3cuceCvA/z1zl9kzpjZ0cuRVOjfm7kA/v7+No3n8GSmwoUL\nx+duA4SGhlKkSBEHrkhERERSu8BjgTQp3UQBtzgNhwfdNWrU4MSJE4SEhBAZGcn8+fMfK4fbz8/P\nrlv/IiIikvotO7aM1uVaO3oZkgYEBwc/MLPicaVoc5z27duzYcMGrly5Qv78+Rk1ahS+vr6sXr2a\ngQMHEhMTQ48ePRg6dGiSxlMjExERkfQtKiaKzaGbOXn1JJ2f7syfN/9k7KaxLDu2jJP9T5IrSy5H\nL1HSCFvjTqftSJkUCrpFRETSp+jYaGbsmsHHv39Mvmz5yOOWh5NXT3Lz3k361+7PWzXeIl+2fI5e\npqQhtsadDj9IKSIiIvI4tp/fTq8VvciVJRcL2y6kdpHaAASHBFPKvRRFcxZNZASRlJfqg24/P7/7\nTpeKiIhI2nPpziXe/eVdfjn1C+MajqPz050TVEHzKuHluMVJmmWv0oFKLxERERGnFh4VzoxdMxj7\n+1i6PN2FES+MIHvm7I5elqQzSi8RERGRNGvT2U20X9yemoVrsrbTWqoUrOLoJYkkS5KC7vPnzxMS\nEkJMTAxWqxWLxcLzzz9v9tqSROklIiIiaU9EdASfbP6EqTum8n3r72lSuomjlyTpVIqll7z//vvM\nnz+fp556igwZMsQ/HxgYaPPktlJ6iYiISNoSFRPFDwd+4KONH1GlQBUmN51MsZzFHL0sEfNLBpYt\nW5YDBw6QObPzdXRS0C0iIpJ27L24lw6LO1AoeyGGeQ6jfsn6No23dClUrgylStlpgZKu2Rp3JtqR\nslSpUkRGRiZ7AhEREZFHuRd9j9EbR9NobiM+9PyQdV3W2RxwX74MnTvD1Kn3v2a1ws8/Q+PGEB5u\n0zQiSZZoTrebmxvPPPMMDRo0iN/ttlgsTJkyxfTFiYiISNp1J/IO3+/7nk+3fUrFfBXZ+fpOiucq\nbpexx42DGjVg+XKYNAksFoiMhBMnYOhQCAmB0aPBzc0u04kkKtGgu1WrVrRq1Sq+DmbcQUpnoYOU\nIiIiqc/qE6t5c8WbVCtUjZktZ/JCiRfsMm5MDEycCPPmwd69ULs27NsHCxfCN9+Auzu0aQOLFoGr\nq12mlDQuRet037t3j+PHjwNQvnx5MmXKZPPE9qCcbhERkdQl9EYo7/36HltDtzKz1UwaPtnQ5jGj\noyEwEFq3hhEjYN06CAiA4sWhTx+YPx/q1YPPPoMnn7TDm5B0yfQ63cHBwXTt2pXixY0/95w7d445\nc+bwwgv2+Y1URERE0oeTV09SZ2YdetfszcyWM8nmms0u427ZAt7e0K0brFhh7Gx7eBiv9ehhpJCM\nGwcZ1Z1EHCjRne5q1aoREBBAuXLlADh+/Dg+Pj7s3r07RRb4KNrpFhERST26Lu1KKfdSjHhhhF3H\nHTUKTp+GHTtg8GAj0BaxN9N3uqOjo+MDbjBKCEZHRyd7QhEREUl/jv19jFUnVnGy30m7j71hgxFs\nz5plHJgUcUaJBt3Vq1enZ8+edOrUCavVyg8//ECNGjVSYm1JooOUIiIizuv4lePsDNvJnH1zGFh7\nIDmz5LTr+Pfuwfbt8NxzCrjFHCl2kDIiIoKpU6eyefNmADw9Pendu7dTNMtReomIiIjzOnPtDHW+\nrUP9EvUpnL0wfl5+ZM+c3a5z/P47DBgAu3bZdViR+5jekdKZKegWERFxTnej7lLvu3p0rdKVAXUG\n2H18qxU2bYIxY6BSJaMWt4iZTMvpbtu2LQsXLqRSpUr31eW2WCzs378/2ZOKiIhI2hR0MojZe2ez\n7c9t1C1al/61+9t1/JgYmDEDJkyALFmMiiVvvGHXKURM8dCd7rCwMDw8PDh79ux9Ub3FYokvIehI\n2ukWERFxvJjYGHaG7WTi1onsubCHoc8NpU6ROlTIVwEXi4vN4//9N+TJA3/+Ce3aQYYMxs52rVrK\n45aUY2vc+dCfBI//F7icNm0aJUqUSPAxbdq0ZE8oIiIiaceiw4vIPzE/PQN7UqVAFQ68dYAe1XpQ\nMX9FuwTcR45AkSLw9NNQpw688gps3Gh0mlTALalJojndVatWZc+ePQmeq1y5MgcOHDB1YUmhnW4R\nERHHWXd6He0Xt2d1x9VU96hu9/FjY+H556F9e6hQwWjb/txzdp9GJElMy+n+6quvmDZtGqdOnaJy\n5crxz9+6dYt69eole0J7U8lAERGRlBd0MoguP3dhYduFpgTc0dEwdKjx+VtvgYvtm+YiyWJ6ycAb\nN25w7do1hgwZwvjx4+Mj++zZs5MnTx6bJ7YH7XSLiIikrKt3r/L22rf57cxvzHp5FvVL1rfr+JGR\nsHo1jB8P2bLB3LlQsKBdpxBJlhQrGXjp0iUiIiLiHxcrVizZk9qLgm4REZGUERMbw4zdMxgZPJLX\nKr7G2AZjecL1iUTvW7QInnkGSpdOfI5bt8DLC9zcoGdP6NJFO9ziPEwPupcvX87bb79NWFgY+fPn\n5+zZs1SoUIFDhw4le1J7UdAtIiJijqN/H6V4zuK4ZXIDoOfynhy+fJivXvqKKgWrJGmMK1egeHGo\nVs1o1f6og4/R0fDyy+DhAdOn65CkOB/TqpfEGTZsGFu3bqVs2bKcOXOGdevWUbt27WRPKCIiIs5r\nzck1eM32os7MOjT4vgFXwq/w3Z7v2BK6hbWd1yY54Ab4/HN47TUIDzfSRP7NaoUzZ+DOHaP2dvfu\nxsHJadMUcEva9NCDlHEyZcpE3rx5iY2NJSYmhvr16zNggP07S4mIiIjjXLpziT6r+rD/r/2MfGEk\nrz71KsPWD+PZb5/lWsQ1NnbbmKR0kjg3bhgB9LZtcO0atGoFNWoYVUgmToRx44zmNvfuQcWKRu3t\nFSsgUyYT36SIAyUadLu7u3Pr1i08PT3p2LEj+fPn54knkv5DJyIiIs7LarUScDCAwWsG07VKV75v\n/X18SsknjT6hlHspCucoTIV8FR5jTBgxApo0+SeXe8wYeOEFI+i+exe2b4dSpeD0aVi50tjpzprV\njHco4hwSzem+c+cOWbJkITY2lh9++IGbN2/SsWNHp6hgopxuERGR5Dt/8zxvrHiDczfO8W2rb6lV\nuJbNY8bEQN++RlAdFAT58v3z2tmzsGwZvP66cVhSJDUx9SBldHQ0jRo14rfffkv2BGayWCyMHDlS\ndbpFREQeU6w11sjbLtkA//r+uGZwtWm8W7dg6VKjPXvevLBkCeTIYafFijhQXJ1uf39/c6uXNGjQ\ngMWLF5MrV65kT2IW7XSLiIgkz8zdM5m1dxabfDfZ1K49MBDefRdCQ6FePRg0CJo21WFISXtM60gZ\nJ1u2bFSuXJnGjRuT9f/JVhaLhSlTpiR7UhEREXGcq3evMmz9MII6BdkUcC9dCm++CQEBRrv2jIlG\nFSLpV6I/Ht7e3nh7eyd4zqJfX0VERFKtYeuH0aZCG54p+Eyy7j9wAL791gi2V62C6vbvAi+S5iS5\nI6UzUnqJiIhI0h26dIh3fnmHM9fOsKXHFnK75U7yvadPw0cfwcaNRpm/bt2MA5HFi5u3XhFnYnp6\nScmSJR846enTp5M9qYiIiKSciOgIPlz3IXP3z2XY88PoVaPXYx+cnDbN6BoZGAjly6s9u8jjSjTo\n3rFjR/znERERLFq0iCtXrpi6KBEREbGP3Rd20/nnzlTMV5HDfQ6TN2veZI2zZYtRa/upp+y8QJF0\nIlnpJdWqVWP37t1mrOexKL1ERETEsPncZpYeXUqrcq14rthzXA6/zMjfRrL4yGImN51M+0rtk30m\nKyIC8uSBv/4C9ceT9Mr09JJdu3bF/5DGxsayc+dOYmJikj2hiIiI2M+fN//kvV/eY9O5TXSo1IFe\nK3sRHhXOrXu36Px0Z472PfpYudsPsnu3kVKigFsk+RINut9+++34oDtjxoyUKFGCBQsWmL4wERER\nebjQG6FM2DKBefvn0btmb2a0nEE212yMaziOfX/tI0fmHDzp/qRd5tqyxajBnaYcPgx9+hhJ6vpt\nQlJAokF3cHBwCixDREREksov2I8vtn9Bj6o9ONT7EIWyF4p/zWKxJLsU4MNs3gyvvWbXIR1r715o\n3BgmTlTALSnmoTndkyZNMi54SP7X4MGDzVtVEimnW0RE0pvN5zbTblE79r65l3zZ8pk+n9UKBQvC\njh1QrJjp09nOaoWQENiwAUaPBj8/6NTpn9ejo6F2bejXz6h7KJJEtsadDy34c+vWLW7fvs3OnTv5\n6quvOH/+PH/++Sdff/21UxyijOPn56fdeBERSReiYqLotbIXnzb+1PSA++JFiIw06nNnygRFi5o6\nnX1YrUYf+jp1YNEi8PGBGTP+eS0gAJo1g1y5oGtXx65VUo3g4GD8/PxsHifR6iWenp6sWrWK7Nmz\nA0Yw3rx5czZt2mTz5LbSTreIiKQnk7ZMYs2pNazptMa07tCbNxtZF+vWGbFpnTrG8059nOvcORgx\nAm7fhhOx6Rk/AAAgAElEQVQnjA4+OXMavzUULmw8HjfOOBE6dCi0aAE5cjh61ZLKmLbTHefSpUtk\nypQp/nGmTJm4dOlSsicUERGRR4u1xjL/4HwqTqtIpWmVGLNxDL+f+52Pf/+Yqc2nmhJwR0XBO+9A\n+/bQqBFcuACLF0NsLLz6qt2ns5+//jIWXKAAvPACrFljBNwArq7Qrp2RThIZCdu2QYcOCrjFIRI9\nSNmlSxdq1aqFt7c3VquVpUuX0lV/khEREbErq9XKqWun2HxuM59u+5TMGTLzaeNPyeaajYADAbwy\n/xX61+5PmTxlkjV+dLQRc27ebKSMjB1r1N4GuHwZWrc2drb37Pnn+Zo1jSwNp/bee9CyJYwf/+DX\nBw+GypXhzTfBpL8OiCRFkprj7Nq1i02bNmGxWHj++eepWrVqSqwtUUovERGR1C48KpzXA18n6GQQ\n2TJlo2bhmnR5ugutyrVKsKMd9/93ydnljo01qo8cPgwNG8K1a3D+PAQFGcF4/frg6WnEramqvfv1\n61CihJFSks/8Q6WSvtkadyarI6WzUNAtIiKpWWRMJC//9DJ53PIwodGEBKX/7OnDD41iHuvWQebM\nEBNjbA6XLGlkZ2TODPPmpcKN4K++gt9+c/KEc0krFHSn3uWLiEg6Fh0bjc8iH2KsMSxsu5CMLolm\nfCbL99+Dv7+RWvLvzeAbN4xDkvnywS+/GIF3qhIbC9WrG9vzjRs7ejWSDpjeBl5ERETsKzo2mtcD\nX+fmvZsEtg80LeCeNw/efdfYDP5v9kXOnEanycyZU2HADfDRR5AtGzRo4OiViCSJgm4REZEU9Ovp\nXxkYNJAiOYrw82s/kzmjORHvpEnw+edGwP3UUw++xt3dlKnNt3y5UX97507IkMHRqxFJkkSD7sWL\nFzNkyBD++uuvBIc4bt68afriRERE0orT104zeM1gDlw6wMRGE2ldvrVptbZHjjSqjmzenEqa2jyO\nw4ehRw8IDDRaZYqkEomeUX7vvfdYvnw5N2/e5NatW9y6dUsBt4iISBJFx0YzcctEas2oRe3CtTnU\n+xCvVHjFtID7m2/gxx+NHe40FXBfuQJTp0KTJkb3nriuPSKpRKI73QULFqRChQopsRYREZE0Zf9f\n++mxvAc5Mufgj55/UCp3KVPnW7EC/Pxg0ybIn9/UqVLGqVNGGsnatXDyJDRtavxG4enp6JWJPLZE\nq5cMGDCAixcv0rp1a1xdXY2bLBa8vb1TZIGPouolIiLijK7dvYb/Bn9+PPAj4xqOw/cZX9N2tuP8\n9ptRizsw0GjAmOrt2gUvvQRduxpt22vXNjpMijiI6dVLbty4gZubG2vXrk3wvDME3SIiIs4kJjaG\nGbtnMDJ4JK+Uf4VDvQ+RL5v5TVt++gkGDDDKVaeJgHvtWujUCaZPN1pliqQBqtMtIiJio6iYKH46\n+BPjNo8jf7b8TG4ymSoFq5g+759/wuzZRh73qlVGt/NUKybG6ODzyy8wa5ZxEvS55xy9KpF4tsad\niR6kDA0N5ZVXXiFfvnzky5ePNm3a8OeffyZ7wqQ6c+YMPXv2pG3btqbPJSIikhzhUeF8uf1LynxR\nhll7ZzG5yWTWd1lvasD955/Qu7fR/bxqVdizx6hS4vQB940bD3/NajXeVL9+EBVlFBBXwC1pTKJB\nt6+vL61atSIsLIywsDBatmyJr6+v6QsrWbIkM2fONH0eERGR5Nh8bjNPf/U0v5z+hfmvzmd91/U0\nKtXItNzt69dhyBCoUgVy5IA1a+DSJVi8GIoVM2VK+7lzx/gt4aefHvz6Z5/B1q1G28yJE+HJJ1N0\neSIpIdGg+/Lly/j6+pIpUyYyZcpEt27duHTpUpIn6N69OwUKFKDyf34FDwoKonz58pQpU4bx48c/\n/spFREQc4F70PYb+OpRXF77KhEYTWOazjNpFzE2kPnDA2NW+fBn274dx46BcOTD5bGbyXLgAI0ZA\no0ZGD3owgu0iRWDgQONNxPnzT+P059SpxgnQ7Nkds2aRFJBo0J0nTx7mzp1LTEwM0dHRzJs3j7x5\n8yZ5Al9fX4KCghI8FxMTQ9++fQkKCuLw4cMEBARw5MgR5s6dy6BBgwgLC3v8dyIiImKyPRf2UGtm\nLY78fYR9vfbxSoVXTJ9z5Uqj0/mYMfDtt1C4sOlT2ubtt+H4cejc2fj8r7+MA5HjxkHHjsYByS++\ngO7d4emnoWxZOHgQihd39MpFTJXoQcqQkBD69evHtm3bAKhbty5ffPEFxR7jb1khISG0bNmSAwcO\nALB161b8/f3jg/Fx48YBMGTIkPh7rl69ygcffMC6devo2bMn77///v2L10FKERFJAdfuXmP4b8NZ\ndHgR4xuOp0uVLqalkWzebKQ0X7kCYWHw669GCsmzz5oynX2dOmWUTzl92siBeecd2L4dQkLgzBm4\ndw8+/hiuXTN2vl9/HfLkcfSqRZLE9JKBJUqUIDAwMNkTPMj58+cp+q82WUWKFOGPP/5IcE3u3Ln5\n+uuvEx3Lz88v/nMvLy+8vLzstUwREUnnYq2xzNk7hw/Wf0Drcq053Ocwud1ymzbfV1/BqFHQoYMR\ni5YsCWPHGvFpqjBxIvTqZQTcAMOHGzvZffpAhgyQNSt89JFj1yiSRMHBwQQHB9ttvIcG3ePHj+f9\n99+nX79+971msViYMmVKsie15+7Av4NuERERe9lzYQ99VvUhOjaawPaB1PCoYdpcsbEwdCgsXQq/\n/w6lzG1caV/nz8PcueDuDvPnw9Gj/7yWM6fRHjPV/NYg8o//bub6+/vbNN5Dg+6nnnoKgOrVqycI\nkq1Wq81Bc+HChQkNDY1/HBoaShH9QIqIiBO4HnGd4euHs+DwAkbXH02Paj1wsSR6BOqxRUfDoUNG\n9sWSJXD7tpFWkqqyLSIjwdvb2JLPksXI2/5v//myZR2zNhEn89Cgu2XLlgBkzZqVdu3aJXhtwYIF\nNk1ao0YNTpw4QUhICB4eHsyfP5+AgIBkjeXn56e0EhERsVmsNZa5++YyZN0QXi73Mod7HyZPVnMi\n4DVrjKIdBQsaKdAvvwzduhlxa6ry/vvGmwgIcNJSKiK2s1eaSaIHKatWrcqePXsSfe5h2rdvz4YN\nG7hy5Qr58+dn1KhR+Pr6snr1agYOHEhMTAw9evRg6NChj794HaQUERE72BW2i36r+xEdG83U5lOp\nWbimaXNt2WIE2T//nEr7v/z6K/j5GYXD79yBXbsgt3l57iLOwta486FB9+rVq1m1ahXz58/Hx8cn\nfpJbt25x+PBhtm/fnuxJ7UVBt4iI2OLv8L/5cN2HLDu2jDEvjsG3qq8pqSRx9u83ylfPmQNNm5o2\njXmOHQNPT6Pk31NPGWX+4g5NiqRxplUv8fDwoHr16ixbtozq1avHT5IjRw4+++yzZE8oIiLiaDGx\nMXyz6xv8gv1oX6k9R/seJVeWXKbOefIkNGsGn3+eygLuKVOMOoaNG8OECUbB8Ndec/SqRFKdRNNL\nIiMjcXV1Tan1PBaLxcLIkSOV0y0iIkm2+dxm+q7uS87MOfmy+ZdUyl/J9DlXrjR6wYwZAz17mj6d\n/axdC76+MHKkkVZSqZLRbVIkHYnL6fb39zcnvSTO8ePH+eCDDzh8+DB37941brJYOH36dLIntRel\nl4iISFJdvH2R9399n3Wn1zGx8UReq/iaaQ1u4kRHw7Bh8MMPxlnDVJXDHRYG1avDjz9C/fqOXo2I\nw9kadyaauObr60uvXr3ImDEjwcHBdO3alY4dOyZ7QhERkZRktVqZvms6lb+qTMFsBTna9yg+lXxM\nD7j/+MOIVffuhd27U0HA/d9g4o03jEY3CrhF7CLRoPvu3bs0bNgQq9VK8eLF8fPzY+XKlSmxtiTx\n8/Oza7cgERFJO0Kuh9B4XmNm7p5JcNdgxjcazxOuT5g2X2wsrFgBL7wAPj5GZ8lVqyBfPtOmtI85\nc+CZZ4w27QA7dsC+fUbHHpF0Ljg42C7NGBNNL6lbty6bNm3i1VdfpUGDBnh4eDB06FCOHTtm8+S2\nUnqJiIg8SKw1lm92fsPw34bzbt13ebvu22R0eWjtgGSJijK6mwcGGpvEVqvRK6ZkSXjvPXj1Vcho\n3ynNcfSoUZGkTBlo0QI++ABat4YGDeABXalF0ivTqpfEmTx5MuHh4UyZMoXhw4dz8+ZN5syZk+wJ\nRUREzLT/r/0MDBrInag7bPTdyFP5nrL7HOHhENc3bts2cHU1esO4uBi72qmmT0xEhFGJZMwYo5Zh\nzZrw9NNGbkwym9aJyIMlutPtzLTTLSIicY79fQy/DX78duY3hj43lD61+ti8uz1sGJw6ZXSLbNgQ\nMmSAq1eNDeEyZWDmTMiUyT7rT3F79xpb8rlywfz5xm8Ko0eDv78RhL/3nqNXKOJUTD9I2ahRI65f\nvx7/+Nq1azRp0iTZE4qIiNjT2etn6b6sO/W+q0elfJU42f8kA+oMsDng3rsXZsyAOnXgww+hRAkj\n88LTE+rVg1mzUmHAHRsL69aBt7dRNPyll+D77//Zmn/nHXj9dXjrLceuUyQNSjTovnz5Mrly/dMw\nwN3dnb/++svURT0OHaQUEUmfwqPCGfLrEKpNr4ZHdg9O9DvBh89/aJeDklYrDBxobPoOGAA7dxoH\nIiMjjTzuCROMVJJUZd06Y3t+8GB48UVjC3/AAMiS5Z9rsmSBadMge3bHrVPEyaTYQcrq1auzZMkS\nihcvDkBISAje3t7s3r3b5sltpfQSEZH0ae2ptby18i1qFa7FZ00+o+ATBe06/sKFRqbF7t1GSkmq\n99dfRnWSmTOhefNUlHQu4jxMP0g5ZswYPD09ef755wHYuHEj06dPT/aEIiIiyXXpziUGrxnM5tDN\nTGs+jWZlmtl9jrt34d13YfbsNBJwW61GO8zu3Y10EhFxiCQdpLx8+TLbtm3DYrFQp04d8ubNmxJr\nS5R2ukVE0odYayyz985m6LqhdHm6C35efmRzzWbKXB99BPv3G7vdacK0aUYC+pYtqTAJXcR52Bp3\nPjToPnLkCBUqVGDXrl0JJonr4FWtWrVkT2ovCrpFRNK+jWc3MnjNYDK6ZOSrl76iaqGqps0VGgpV\nqxo53CVKmDZNyjh0CMaOhV9/hU2boGxZR69IJFUzLb3k008/ZcaMGbz99tsPbJX722+/JXtSe/Lz\n88PLywsvLy9HL0VEROzo1NVTvPfre+wM28m4BuNMb90eGwtvvw29e6eBgHvkSPj6axg0CL76CnLk\ncPSKRFKt4OBguxTtUJ1uERFxKtcjrjN642hm7Z3F28++zaA6g3DL5GbqnOfPQ5cuRpfJ1ashmzmZ\nKynjxx+NAuPbt4OTpIOKpAWm7XQvXrz4kTsK3t7eyZ5URETkv2KtsczaM4sP139Iy7ItOdT7kN2r\nkjzIzz9Dr15Gx/OhQ5348OSoUVCwILzxxsOv2b3bKAO4bp0CbhEn89CgOzAwUEG3iIikiJ1hO+mz\nqg8uFhdWdVxFtULmnxv6+2+j6c2vv8KyZUYTHKdltRrl/m7ehNq1oUqVf147d844+XniBHz+uZFO\n8vTTjluriDyQ0ktERMRhroRf4YN1H7Ds2DI+bvAxXZ/piovF/l1nrFbYvBm2bYMdO4yPv/+Gtm3h\ns89SQcrzkSPQpInRnv3jj42Tnm5uRoA9YgTUrGk0vnnhBWjTxtGrFUmTTKteEufvv//G39+f33//\nHYvFgqenJyNGjCBPnjzJntReFHSLiKROMbExzNw9kxHBI3it4muMqj+KXFlyJX5jMkRHG13Ng4Oh\naVMjPq1ZE8qVS0VdJSdPNqqRzJgBnToZnSNdXIwygMuWQalSjl6hSJpnenMcHx8fXnjhBZYsWYLV\nauXHH3/ktdde49dff032pCIikn6tObmGd355h9xuuVnbaS1VClZJ/KaHiImBe/cga9YHv373LnTo\nALdvG+nOqba7+Zo10KOH8fnUqUZdw2eega1bU/GbEklfEt3prlSpEgcPHkzwXOXKlTlw4ICpC0sK\ni8XCyJEjVTJQRCQVOHjpIO+sfYdT104xodEEXi73ss0lAEeOhE8/NQ5CDhoEHh7/vHb9Orz8MhQq\nBN9/D66uNr4BR4mIgHz5jNxtd3fjubt3jd1utXMXMV1cyUB/f3+bdroT/cNa48aNCQgIIDY2ltjY\nWObPn0/jxo2TPaG9xdXpFhER53Tx9kXeCHyDF+e8SLPSzTjU+xCty7e2OeD+80/48ksICjJK/VWq\nBD17wrFjcOGCkd5cpYpRQS/VBtxgNLapXPmfgBuMfG4F3CIpwsvLCz8/P5vHSXSn+4knniA8PByX\n/ye+xcbGku3/BUwtFgs3b960eRHJpZxuERHnFRkTyaQtk5i0dRJdq3Rl2PPDcHdzT/zGJOrWzdjZ\nHjvWeHzlihGET51qNLoZPNgoAZjqY9N33jFSSEaOdPRKRNI10w9SOjMF3SIizmnzuc28ueJNiucq\nzpSmUyiV274H/fbsgWbN4Pjx+yuP3LljPF/VvG7x5rlyxahhWLv2P20xK1c2DlA6dU1DkbTPtKD7\nyJEjVKhQgd27dz/wxmrVzK+hmhgF3SIizuV6xHWG/DqE5ceWM7npZNo+1dburdutVmjYEF591ahK\nkibMm2fU4d6zxwiu9+wxygHWqWPU3L50yYm79oikD6ZVL/n000+ZMWMGgwcPfuB/MH/77bdkTyoi\nImmL1Wpl0eFFDFwzkJZlW3K4z2HTSgCuWgVhYfD666YMn/KOHDFyYb791vhtws3NKCTu4wO5c0OD\nBgq4RdKARNNLIiIiyJIlS6LPOYJ2ukVEHO/s9bP0WdWHM9fPML3FdOoVq2faXNHRxsbvJ59Aixam\nTZOy2rc33tTQoQmfv3HDaOnepg20bOmYtYlIPFvjzkSrl9StWzdJz4mISPoSER3BmI1jqDa9Gs8W\neZY9b+4xNeAGIwOjUCF46SVTp0k5Bw/C+vXQt+/9r+XMCbNnK+AWSSMeml5y4cIFwsLCCA8PZ/fu\n3Vit1vhqJeHh4Sm5xkeKKxmosoEiIinDarWy4vgKBq0ZROUCldnx+g6edH/S9Hlv3gR/fyO9JNVX\nJInj7/9PdRIRcUpxdbpt9dD0ktmzZzN79mx27dpFjRo14p/Pnj073bp1w9vb2+bJbaX0EhGRlHX8\nynEGBA3gzLUzTGk2hcalUq5vwwcfGLncs2en2JSJW7rUqKP9/PPGh/tjlETctw+aNIFTp+D/pXhF\nxHmZXjJw8eLFtGnTJtkTmElBt4hIyrh17xajN47m2z3fMvS5ofSr3Q/XDOZ3nImONuLaL76AM2dg\nyxYoUsT0aZNm924jaO7d22jHvnUrlC5ttL+sXDnx+195xQjUBw0yf60iYjPTgu5JkybFD/7v6iVx\njwcPHpzsSe1FQbeIiLmsVis/HfyJd395l4ZPNuTjBh9TKHuhJN8/bRoULGjEl4+TEnL5slGa+quv\noHhx6NcPvL0hU6ZkvAkzXL8O1asbnXlee814LjISJkyAbdsgMPDR92/aZFQnOXnSqFYiIk7PtJKB\nt27demCpwP8G4SIikjadvX6Wt1a+xZ83/2Rh24U8W/TZx7p/82YYPRry54fx4+Hjj+HFFx99z44d\nRqC+dKkRqC9f7oRNbqxW8PU1uvPEBdxg9Jp/+21jt3vnTvhXamYC4eHQvbvxRhVwi6Qb6kgpIiIJ\nxMTGMG3HNPw3+DOoziDeq/cemTI83hZzZKQRLPv5GRXv5s+H4cOhVCljc7h69X+uDQ83Xp82Df7+\nG3r1gh49IG9e+74vu5k0yVjwpk2QOfP9r3/5JaxZ8/Dd7kGDjGY3P/xg7jpFxK5Mz+n29fV94KTf\nffddsie1FwXdIiL2dfjyYXou70lGl4xMbzmd8nnLJ2uc0aP/ybKI++NoZKTR/+Wjj8DT0wiuAwON\nFOg6dYzU6CZNnLgPzK1bsGgRvP8+bN/+T5v2/4qIMHa7ly69f7f799+hXTs4cADy5DF9ySJiP6al\nl8R56aWX4tNJ7t69y88//4yHh0eyJxQREecTGRPJx5s+5ssdX/JR/Y94o/obuFgSbeXwQMePw+TJ\nsGtXwjxuV1ejbXuXLjBlCvTpA61aGSklJUva6Y3Y27VrRvC8ZAls2AD16sHixQ8PuAGyZDEC81Gj\njPyYOHFpJVOnKuAWSYceO70kNjaWevXqsXXrVrPWlGTa6RYRsd3W0K28Hvg6T7o/ybSXplEkR/LL\ng1itRtfyli1TeVGOiAjjN4MJE4wKI23aGB15cuZM+v2lShlB9zPPGMXFJ0yAwoUhIMDctYuIKUzf\n6f6v48ePc/ny5WRPKCIizuHi7YsMWz+MVSdWMbnpZNo+1dbmg/Jz5hhNbPr1s9MiU1psLPz0k1EU\nvGpVIx2kXLnHHydLFhgyxDhweeuWcZq0T5+EBy9FJF1JNOh+4okn4v8jbLFYKFCgAOPHjzd9YSIi\nYo570feYvG0yE7ZMwPcZX470OULOLEncwX2Ey5eNrIrVqyHjY2/pOIHYWGOb/s4dI9H8+edtG+/1\n1+Gvv+Dll6FmTfusUURSLVUvERFJJ6xWK0uPLuWdX96hUv5KTGw0kTJ5ytht/M6doUABmDjRbkOm\nrPXrYeBA2LsXXJKXzy4iaVeKpJfs37+fkJAQoqOj459zhjbwAH5+fnh5eeHl5eXopYiIOK39f+1n\nYNBALt25xDctvqHhkw3tOv4vvxgV9A4dsuuwKWv6dHjjDQXcIpJAcHAwwcHBNo+TpJKBBw4coGLF\nirj86z9Es2bNsnlyW2mnW0Tk0S7ducTw9cNZemwpI18YyRvV3yCji31zP8LDja7nX3wBzZvbdWj7\n2bbNqK/9008Prkl4+TKUKQMhIZArV4ovT0Scn+k73X/88QeHDh1SF0oRkVQkKiaKL7d/ydjfx9Kp\ncieO9jmKu5u7KXP5+xspy04bcJ88abS3zJ4dvvvOyLX+rzlzoHVrBdwiYppEg+6aNWty+PBhKlas\nmBLrERERG609tZYBQQMonrM4m3w3JbvBTWKuXDFSoLdsMVq+O6XLl4127X5+xm8GL70EPj5GAB7H\najVSS2bPdtQqRSQdSDTo9vX15dlnn6VgwYJk/n+7W4vFwv79+01fnIiIJN3pa6cZvGYwBy4dYHKT\nybQo28K0v1IuWQJ9+xrNFffvh2zZTJnGNnfvGt132raFN980nmvSBMaNgzFj/rluwwajc8+zzzpm\nnSKSLiSa012qVCk+++wzKlWqlCCnu8SjunGlEOV0i4jAncg7fPz7x3y18yveefYdBj07iCwZs5gy\n16VLRrC9b5+RqVGvninT2C4mxgi23dxg7tx/DkeePw9PPw27d0Px4sZzHToYfej793fcekXE6dka\ndyYadD/77LNO0X3yQRR0i0h6ZrVamX9oPu/+8i7PF3+e8Q3H29RN8tFzGY0UBw2Cbt2MbA03N1Om\nejxjx8KMGfDkk8ZByLJljX+DguDwYePf//+VNp6fH5w4AT/8AH//DaVLw+nTkDu3Q96CiKQOpgfd\nvXv35vr167Rs2RJXV9f4SZ2hZKCCbhFJr/Zd3Ee/1f24FXmLKU2n4Fnc07S5wsLgrbfg1CmYNcuJ\n+ryMG2ccgAwIMLbgT5yA48eNf8EIqt0fcHj0zh2jy+TixUYy+t69RjMcEZFHML16SXh4OJkzZ2bt\n2rUJnneGoFtEJL25dvcaw9YPY9GRRYzyGkXPaj3J4PKAEnh2snGj0bn89ddhwYL7N40d5rPP4Ntv\njXxsDw/jucaNk3ZvtmwwejQMHgxXrxo75SIiJlNHShGRVMBqtTJv/zze+/U9WpdrzZgGY8jtZl46\nhNUK06bBqFEwbx40amTaVI9v2jSYMMEIuIsVS94YsbFQowZERBgdfVQWV0QSkSIdKUVExHEOXTpE\n71W9uR15m2U+y6hVuJap8927B336wB9/GOUAS5Uydbr7RUcbDWweFAh/+y2MHw/BwckPuME4WPn9\n90bdQwXcIpICFHSLiDip25G3GbVhFLP2zsLvBT961ehlaioJwIUL0KYNFCoEW7fCE0+YOt2D1apl\n5GaXLm18lClj/HvjBnz6Kfz2G5Qsafs8lSrZPoaISBK5POrF2NhYFixYkFJrERERjFSSJUeW8NTU\np7h4+yIH3zpIn1p9TA+4t237p3/MokUOCrhPnzZOboaGGrvabdsaOdibNsEvv8DatUYQLiKSyiSa\n0129enV27dqVUut5LMrpFpG05tTVU/Rb3Y+zN84ytflUvEp4pci8330HQ4YYcW7Lliky5YN99plR\n6k+HG0XEydgadz5ypxugUaNGTJw4kdDQUK5evRr/ISIi9hMRHYF/sD+1ZtbCq4QXe97ckyIB9507\n0K+fkSa9caODA26AZcvg5ZcdvAgREftLdKe7RIkS97URtlgsnD592tSFASxbtoyVK1dy8+ZNevTo\nQaP/HJ/XTreIpAWrT6ymf1B/KuevzOSmkymW04YDgkkUEQHffGOUun7xRZg6FXLlMn3aR7tyxWhy\nc/Gik3TeERH5h+nNcZzB9evXeeedd5g5c2aC5xV0i0hqFnI9hEFrBnHw0kGmNJ1CszLNTJ8zMtJI\nJRkzBqpXB39/qFLF9GmT5vvv4eefjQ8RESdjenpJZGQkn3/+OW3atOHVV1/liy++ICoq6rEm6d69\nOwUKFKBy5coJng8KCqJ8+fKUKVOG8ePHP/T+0aNH07dv38eaU0TEWUVERzB642iqT69O9ULVOfDW\nAdMD7uhoo5tkuXJGBseSJbB0qRMF3KDUEhFJ0xLd6e7RowfR0dF07doVq9XK3LlzyZgx4327zo+y\nadMmnnjiCbp06cKBAwcAiImJoVy5cvz6668ULlyYmjVrEhAQwM6dO9m9ezfvvvsuhQoVYsiQITRu\n3JgGDRrcv3jtdItIKhOXSlIpfyU+a/IZJXKVMHW+mBj46SdjR9vDw2jE+Nxzpk6ZPBERUKCA0Ws+\nb0nqhFAAACAASURBVF5Hr0ZE5D6mN8fZsWMH+/fvj3/coEEDnn766ceaxNPTk5CQkATPbd++ndKl\nS1OiRAkAfHx8WLZsGUOGDKFz584ATJkyhXXr1nHz5k1OnjzJm2+++Vjziog4i7hUkgN/HWBKsyk0\nL9Pc9DlXrYJ334WcOeGrr4zcbaftA7N+vbHtroBbRNKoRIPujBkzcvLkSUqXLg3AqVOnyJjR9p46\n58+fp2jRovGPixQpwh9//JHgmv79+9O/f/9HjuPn5xf/uZeXF15eXjavTUTEXiKiI5i4ZSKfbfuM\nQXUGEdAmgCwZs5g659mzMHAgHDxoVOB76SUHBtsnT0KRIpAlkfes1BIRcTLBwcEEBwfbbbxEo+cJ\nEybw4osvUvL/3b9CQkKYNWuWzRP/tyJKcv076BYRcSZxqSQV81Vk1xu7TE8liYyESZNg4kQj6A4I\nSDzWNdW1a1CnDnTsCJ9//vDrYmNh+XKjZqGIiJP472auv7+/TeMlGnTXq1eP48ePc+zYMQDKlStn\n04RxChcuTGhoaPzj0NBQihQpYpexRUQcyRGpJOvWQZ8+Rrf0HTuMynsON3IkNG4MCxYYgXetWg++\nbscOcHdXp0kRSdMSrV5St25dsmTJQpUqVahSpQpZsmShbt26Nk9co0YNTpw4QUhICJGRkcyfP59W\nrVo99jh+fn523foXEUmue9H3GLNxDNWnV6dawWoc7H3Q9IA7LAzat4cePYwGN4GBThJwHzxonOCc\nMsXYfn/9dXhY5SulloiIEwsODrZLZsVDq5dcuHCBsLAwOnbsyI8//ojVasVisXDz5k169erF0aNH\nkzxJ+/bt2bBhA1euXCF//vyMGjUKX19fVq9ezcCBA4mJiaFHjx4MHTr08Rav6iUi4iTWnlpL31V9\nqZCvApObTKake0lT54uOhi++MOptv/kmfPghZM1q6pRJZ7VCo0ZGIN2vn/G4WTOoXx/ef//+6ytW\nNPrP16mT8msVEUki05rjzJkzh9mzZ7Nz505q1KgR/3z27Nnp1q0b3t7eyZ7UXhR0i4ijhd4IZdCa\nQey+sJspzabQomwL0+fcuBH69oWCBeHLL6FsWZMn7NPHCIx7907a9T//DMOHw969EHfw/swZqFkT\n/vgDSpX659qTJ8HTE86fB5dE//gqIuIwpnekXLx4MW3atEn2BGZS0C0ijhIZE8lnWz9jwpYJ9K3V\nl/frvY9bJnNblx85AkOGwP798Mkn8OqrKVCVZPduo/xJ5szG5L16Pfr6iAh46imYMQP+219h4kQI\nCoJffvln4ZMmwbFjMH26OesXEbET0ztS7ty5k+vXr8c/vnbtGsOGDUv2hPamnG4RSWnrz6ynytdV\n2HhuI3/0/AM/Lz9TA+6LF41Y94UX4Pnn4ehRaNs2hcoAjhgBH3xg1NEeO9ZIA3mUSZPgmWfuD7jB\nKKly5QrMnfvPc8rnFhEnZ3pOd5xnnnmGvXv3JniuatWq7Nmzx+bJbaWdbhFJSedvnuedX95ha+hW\nPm/6Oa3KtbJb+dMHuX3biGGnTAFfXyP2zZ3btOnut20btGsHx48btQdPnDDysseMga5d77/+zz+N\ngHvHDij5kJz2nf9r787jbK7bP46/UMiNtCBbd8q+JMKdZRghWQZJ1ky2lH0pil8YW3cRkjWSrXvE\nLQwxZmwzspN93xpZKm4l0zAzzsz398cnJGnMzPnO95wz7+fj4ZE5c873c83ky3U+c32ua6fZOT9w\nwLxreOop+Oknh3sbiogkzfad7sTERGJjY29+fO3aNeLj41O8oIiIt7mecJ2xm8dSblo5nnroKQ51\nP0STEk1sS7hdLlNtUayYyXd37jSVGWmacIPZ5X7vvVsJcdGisGYNDBwIwcF3Pv+dd8yW/N0SboCK\nFU37wLfeghUroE4dJdwiki4k2ae7bdu21K5dm44dO2JZFrNmzSIwMDAtYrsnQUFBmkQpIraJjIqk\n+8ru5M+Rn82dNlPsEftOLVoWfP21yV3z5jXzYv5wjj1tbdhgDjl26HD74yVKmJrsOnXg/vtNnQvA\npk0QGQmffpr0tYcPhzJlYMsWk9SLiHgwd02mTLK8BCA0NJS1a9cCULduXerVq5fqhd1B5SUiYpcf\nf/uRt8PfZsPpDYyvN55mJZvZWkqyYwf07w8XL5pDkg0aODi63bLA398k3O3b//Vz9u6FevVg6lRT\nk12pEvTrZ3ax78XKleZ1P/wAjz7qrshFRGxje/cSMKPfjx8/Tt26dbl69SoJCQnkyJEjxYu6i5Ju\nEXG3hMQEpuyYwvANw+lUvhPv1XiP7Jmz27bexYumKUhoqNkAbt/+Vpc9x6xZY9oEHjz498Hs2mX6\nbzdoYOpgNm5M3juFH36AfPlSH6+ISBpIbd6Z5F/t06dPZ8aMGfz888+cPHmSs2fP0rVr15s73yIi\nvmL7ue10XdGVHJlzENk+klK5S9m2VkKCqcQICoJ27UxHkpw5bVvu3lmW6bEdFJR09l+hgqmHadLE\njMJM7ta8Em4RSUeSTLonT57M9u3bee73SWHFihXjwoULtgcmIpJWfrn2C4PWDmLp0aWMrjOaV59+\n1dZSki1bzEZyzpymE1+ZMrYtlXyhoRAdDS1b3tvzK1Uyg20cq4UREfEOSXYvyZIlC1myZLn5scvl\nsvUfo+RSn24RSSnLspi3dx6lppgd7UPdDtGuXDvb/o67cAE6djRDbd5+G9av97CE+8Yu97BhyZsO\n6UH/JoiIuFua9enu378/uXLlYu7cuUyaNIkpU6ZQqlQpRo0alerFU0s13SKSUocuHqLbim5Ex0cz\nteFUKheobNtaCQkwbZrJZdu1g6FDPaSU5M+WLDGF5d9+q5HsIiJ/YvtBysTERD777DPCw8MBqFev\nHp07d/aI3W4l3SKSXDHxMYzcMJLPdn/G0JpD6VqxK5kyZrJtvT+Wkkya5GE723+UmAjlysG//w2N\nGjkdjYiIx7E96V68eDENGza8rcTEUyjpFpHkWHZ0Gb1Ce1G1UFXGvjCWfDnsO8h34YLptx0ebgbb\ntGrl4VUYCxbA+PHmXYJHByoi4gzbJ1IuW7aMokWL0q5dO77++mtcLleKFxMRccLpy6dp8mUT+q/u\nz8zGMwl+Odi2hNvlurWj/cgjcPgwtG7t4Xmsy2VqXoYP9/BARUS8V5JJ9+zZszlx4gTNmzdn/vz5\nPPnkk3Tq1CktYrsnOkgpIncTnxDPBxs/4Nnpz1IpfyX2vbmP2k/Wtm29yEh49llYvNgckvzoIw+t\n3f6z4GDInRvq1nU6EhERj5NmBylviI+PJywsjM8//5wNGzZw6dKlVC+eWiovEZG7iYyKpNvKbjyR\n6wkm1p/Ikw89adtaZ8+aaZKbN5tEu3lzL9owjouD0qVh5kyoWdPpaEREPJbt5SUrV66kffv2FC1a\nlEWLFvH666/z008/pXhBERE7XYi5wGtLX6PdknaMqDWCr1t/bVvCHRsL778PzzwDRYuaUpJXXvGi\nhNuyoGdPKF9eCbeIiM2SHI4zb948WrZsybRp08iaNWtaxCQikiJfHfqK7iu78+rTr3Ko+yHbxrdb\nlhnE2LcvlC0L27fDk/ZtpNtn2jSzPb9li9ORiIj4vHsuLwFTP+2OmhZ3UXmJiABcjr1Mz9CebD27\nlTlN51C1UFXb1jp2DPr0gVOn4JNP4IUXbFvKXhs2mG35zZvhqaecjkZExOPZXl7yRyEhISleSETE\nDqtPrubpqU/zYJYH2fPGHtsS7uho0wKwalV4/nnYt8+LE+7vvzdj3r/4Qgm3iEgaSbK8RCTVzp2D\nHDm8pI2DeIuY+BgGrB7A8mPLmdl4JnWfsqfzhmWZ5h7vvAO1a8P+/ZDPvvbe9rt6FZo2NSc/1a1E\nRCTNJGune+fOnYDpZOIp1DLQCyxcCAULwosvmhrS8+edjki83OYzm3nm02eIjo9mX9d9tiXcu3eD\nnx+MG2f+GM+Z4+UJt2VB586mW0nfvk5HIyLiFdKsZWDNmjWZPXs2hQsXBmD79u107tyZffv2pXrx\n1FJNtxeJjoawMFi6FFauNK0emjQxO24lS3pRuwdxUpwrjqCIIGbtmcWUhlNoVrKZLev88AMMGwZL\nlsDIkdCxI2Syb1J82hkzxrx72LABHnjA6WhERLxKavPOJMtLBg0aRP369enZsyfnzp0jNDSU2bNn\np3hBSady5DDNi5s3h+vXzT/6ISFm9ztLFpN8N2kCVar4SHYj7rb3x70ELg2kcK7C7H1zL3mz53X7\nGpcuwejRMGMGdOgAR47AQw+5fRlnrFplxrxv26aEW0TEAffUvWT9+vXUrVuX3Llzs3v3bh577LG0\niC1J2un2AZYFe/aYHfClS+HHH6FRI5OE16mj5EBwJboYs2kM47aO46O6HxFYLpAMbv7JSHQ0fPwx\nTJhg3he+956piPIZx49D9erw1VfmvyIikmypzTuTTLpHjBjBggULmDFjBvv27WPcuHGMHTuWRo0a\npXhRd1HS7YO++87sgIeEwK5d5uRakyYmEX/kEaejkzR2/NJxApcGku3+bMxqMovHH3zcrdePjYWp\nU+HDD817vKAgKFLErUs478oVeO450+ewSxenoxER8Vq2J919+vTh3//+Nw/8vuN4+vRpOnfuzOrV\nq1O8qLso6fZxly7BihVmB3ztWqhQwSTgTZrA72cMxDclWolM2TGFoIgghtYcSvfK3cmYIVnnvv+W\nywWzZ8Pw4eaP1YgRZsiNz0lMhJdegvz5zbsLERFJMduTbk+mpDsduXYN1qwxO+DLl8Njj92qAy9f\nXgcxfciZX8/QcVlHrsRdYW7TuRR/tLjbrp2YaM4RDhliykfef99sAvuk7dtNn0PLgvBwyJzZ6YhE\nRLyabUl37969mTBhAgEBAX+56LJly1K8qLso6U6nEhJg69ZbdeBxcbc6odSoAfff73SEkgKWZTFv\n3zzeDn+b3v/qzTvV3+G+jO4ZJWBZ5ocm//d/kDWrSbZr13bLpT3PkSPmC922DYYONSdC79NIBhGR\n1LIt6f7222959tlniYyMvGOBDBkyULNmzRQv6i4ZMmRg6NCh+Pv74+/v73Q44gTLgsOHzQ740qXm\nwFiDBiYJf/FF0zVFPN7FmIu88fUbHP/5OHObzqV8vvJuu3ZkJAwaBL/+CqNGQePGPvqDkbNnTVH6\nsmVm8E2PHjqILCLiBhEREURERDBs2DD7yktcLheBgYEEBweneAE7aadb7nD+vEk6QkJg0ybTqaFp\nUwgI8PKpJr5r6ZGldF3RlcCnAxleazhZ7suS8otZlhkZ+dhj7Pw+D//3f3DihKndbtXKR7tRXroE\nH3wAn38Ob7wBAwZArlxORyUi4nNsr+muXr06a9euJUuWVPxDaBMl3fK3rlwxvYmXLoXQUChR4lYZ\nSokSTkeX7l2OvUyfVX3Y+P1G5jSdQ7XHq6X8YufPw3/+A3Pn4rpwiZ+u5uCF7FvoOfRhOnb00XLm\nmBjT43DcOHjlFRg82ByYFBERW9iedLdr144jR47QuHFjsmXLdnPRfv36pXhRd1HSLfcsPt7UGSxd\nanbB//GPWwcxn3sOMrqvM4YkbdWJVby+/HUaFW3EmBfGkD1z9uRf5OpV8/9z7lxzaLBZMzY+GcjL\n46uzvOQAnrV2kGlNuBm+5EtiY+Gzz0xheo0apvVK0aJORyUi4vNsT7pvzJr/8zCKoUOHpnhRd1HS\nLSliWfDtt7fqwC9eNOUnTZua03VZszodoc/6NfZX3gp/izWn1vBZ48+o82Sd5F0gMRG++QbmzDEz\n2p97DgIDSQxoQtDobMyeDYsWQeWKiWbKTfbs5rm+UMR99SpMn25GuT/7rKnfrlDB6ahERNIN25Pu\nhQsX0qJFiyQfc4KSbnGLkydvDeTZu9dMSWnaFBo29KEZ4M5bfXI1nZd35sWnXmTMC2PImSXnvb/4\n+HGzoz1vnjkc+9pr0LYt5MvH5cvmtzExph1gnjy/v+bqVfD3N4OVhgyx40tKG7/9BtOmwdixULWq\nGZdZ3n0HTUVE5N7YnnSXL1+e3bt3J/mYE5R0i9tdvHhrIM/69VCx4q2BPP/8p9PReaXouGj6r+7P\nyuMrmREwg3pF6t3bC3/5BRYsMMn2yZPQpo1JtsuVu7lzfeCAmf3SsKHZAL6jW+SPP5rd8JEj4dVX\n3fuF2S06GiZPhvHjzZuH997z0Qk+IiLeIbV5512bt4aGhrJy5UrOnTtHr169bi4SHR3N/eqDLL4q\nd25o3978unoVVq82O+AjRphpKjfqwP+Q+MndrftuHZ2WdeL5J55nf9f9PJj1wb9/wfXr5vDr3Lnm\ne1+vnuk5/cILd2TU//0vdOtmctK75tOPPQZffw3PPw+PP25qoD3d5cswcSJ88on5+iMioGRJp6MS\nEZFUuutO9969e9m9ezdDhgxhxIgRN5PunDlzUqtWLR7ygB+7a6db0kxCAmzefGsgT0LCrU4ofn4a\nPvInv8X/xjur32HZsWV82uhTGhRtcPcnWxbs3m0S7fnzoUgRs6P9yit/Wd7jcpk8fOFCWLz4Hist\nVq+Gdu1gwwYoVizlX5idfv7ZdCOZPNmUxAwa5LmxioikQ7aXl8THx5PZQ/ttKekWR1gWHDx46yDm\nqVOmvqFJE7MzmT0FnTh8SGRUJB1COlDjnzUYX288Dz1wlzfo8fFmN3fOHFO3HBhoEuMiRe567UuX\nTL9tMPn5o48mI7AZM2D0aNiyJZkvtEF8vPkztGvXrV8HD0KLFjBwIDz1lLPxiYjIHWxPuo8dO8ag\nQYM4dOgQ165du7noqVOnUryouyjpFo9w9qwZyLN0qRlPX6PGrYE8efM6HV2aiYmPYdDaQSw6vIhp\nDacRUDzg71/w9tuwY4cp3alePcm2jbt3Q7NmJi8dNSqFP1x4913YuBHWrEm7LjVXr8K+fbcn2EeO\nmMS6fHnTgaRCBXjmGciZjMOlIiKSpmxPuqtVq8awYcPo168fy5cvZ9asWSQkJDBixIgUL+ouSrrF\n41y+bAbxhISY2uTSpW+VofhwqcDG7zfSfml7qhSqwoQXJ/DwAw///QvCwqBzZ5NJ38Ou8xdfQN++\npvIiVY2TEhPNVvl995lhOnbU5btcZjrkxo0mwT51ytRk30iuK1QwByJ/n3sgIiLewfaku0KFCuza\ntYuyZcuyf//+2x5zWoYMGRg6dCj+/v74+/s7HY7I7eLizCG4pUvNTnjOnCb5fu01n5mIefX6Vd5b\n9x5fHviSKQ2n0LRE06Rf9NNPZof3P/+BWrX+9qnXr0P//uYs5JIlbmrece2aOVhZt66ZD+9O331n\n+hdmzWr+W6GCeePloSV6IiKStIiICCIiIhg2bJi9SXfVqlX55ptvaN68ObVr1yZ//vwMHDiQo0eP\npnhRd9FOt3iNxETYudMk4J99Zg4Nvvii01GlypYzW2gf0p4K+Sowsf5EHs12D3XSiYnQoIFpxThy\n5N8+9cIFs6udLZvJz916dvvCBahSxfTvfu0191zzxnb8wIHQp4+mnIqI+Bjbd7q3b99OyZIluXz5\nMoMHD+bKlSsMGDCA5557LsWLuouSbvFKmzeb5tIffWQODnqZa9evMWT9EL7Y/wWT6k/i5VIv3/uL\nx441IyM3bPiLptq37NgBL79s8uGgIMiUKfVx3+HwYdP/+ssvk9xx/1u//mp6F+7aZU53PvOM20IU\nERHPYXvS7cmUdIvXOnzY7HT36GEOFHpJz+9tZ7fRPqQ9ZfKUYUqDKeT+R+57f/G335qveft2KFz4\nrk/7/HNz3nH6dFONY6v1602Nd2Rkykp+Nm0yTcJffNG8oVCdtoiIz7I96d6xYwfvv/8+UVFRuFyu\nm4vu27cvxYu6i5Ju8WrnzplkrU4dk7B5cDlCnCuOoIggZu2ZxSf1P6FF6WSeZoyONvXNI0dCy5Z/\n+ZT4eFOVsW6dqd9Os3kws2ebDipbt5rhSPfC5TJfy7Rp8Omn5rCsiIj4NNuT7mLFivHRRx9RpkwZ\nMv4hKXjiiSdSvKi7KOkWr/fLLyZhK1DAJH9Zsjgd0R12nt/Ja0tfo/gjxZnacCp5s6egDWL79qZG\nZObMv/z0Dz9A8+Ym550714HOeYMHw9q15tcDD/z9c28clsye3fQYz5cvbWIUERFHpUnLwE2bNqV4\nATsp6RafEBtrkrhffzUjFj2kV3OcK44RG0YwY9cMPq73Ma3KtCJDSspggoNh2DBT8/yPf9zx6c2b\nzYHJN94wkyYd2fC3LPP/ICHB1GXfLQgdlhQRSbdsT7rDw8NZsGABderUuTmZMkOGDDRr1izFi7qL\nkm7xGQkJpr5761bT5/uxxxwNZ9cPu2i/tD2FHyrMtIbTyJcjhbu5p07Bv/4F4eF3zGu3LFOZMWQI\nzJplhno6KjbWlPrUqAHvv3/753RYUkQk3Utt3pnkTLc5c+Zw9OhRXC7XbeUlnpB0i/iMTJlgyhQz\narFaNTNYp2jRNA8jPiGe9795nyk7pjCu3jjalm2bst1tME22W7c229d/SrhdLujVy5xf3LTJkS/1\nTlmzmpaOVaqYaZGdOpnHbxyWrF/fHAbVYUkREUmBJJPunTt3cuTIkZT/wysi9yZDBnjvPbPLXaOG\nGahTqVKaLb/17Fbe/PpNCuYsyJ4395A/R/7UXXDIEDNtsnfv2x7+9VdTTpIxI2zZ4jHVNMajj8KK\nFeDnZ+rst241hyWnT4fGjZ2OTkREvFiSSXfVqlU5dOgQpUuXTot4RKRzZ8iTxwyRmTfP9iE6F2Iu\n8O6ad1l1YhWj645O3e72DWvXmhORu3ff1g7xu++gUSMzEHL8eDON3eMUKwb//a+ZWFmzpvkadFhS\nRERSKcma7hIlSnDy5EkKFy5Mlt87K6hloEgasHmIjivRxbSd0xgWOYzApwMZ6j+UnFncsO188aIp\nJ5k1yySuv9uyxQy8GTgQevZM/TK2+/57KFhQhyVFRARIg4OUUVFRf/m4WgaKpIEbQ3S6d4f+/d02\nRGfj9xvpvrI7jzzwCBPrT6R0Hjf9JMuyICAASpeGDz+8+fD8+abKZPZss4EvIiLibTSR0nvDF7k3\nbhyi80P0DwxYM4CIqAg+qvsRLUq3cO95jU8+MW31Nm6EzJmxLBg+3Gx6L18OZcu6bykREZG0lNq8\nUz83FfF0BQrAhg2mc0abNhAXl+xLXE+4zrgt4yg7tSwFchTgcPfDtCzT0r0J9549ZrLj/PmQOTOx\nsabpR2ioOY+ohFtERNIzj026jxw5QteuXWnRogUz7zLFTiTdeOgh0+v6+nVTn3Hlyj2/dP136yn/\naXnCToaxqeMmPqjzAdkzZ3dvfDExpj3g+PHw1FNcuAC1a5vWgOvXO952XERExHEeX16SmJhIq1at\nWLhw4R2fU3mJpDsJCeYU4pYtSQ7ROXvlLG+Hv82Ws1sYX288L5V4yb7Wn6+/bnbg587l0CHToaRt\nWzOIUucQRUTEF3h8eUnHjh3JmzcvZf/0s+VVq1ZRokQJihYtyod/OHD1R8uXL6dhw4a0atXK7jBF\nvEOmTDB5smkDUrUqHDt2x1PiE+L5cOOHPDPtGYo+XJTD3Q/TrGQz+xLuhQshIgImTyY8HPz9TbI9\nYoQSbhERkRts3+n+5ptvyJ49O4GBgezfvx+AhIQEihcvzpo1ayhQoACVKlVi/vz57Ny5k127dtG/\nf3/y5781mKNJkyaEhITcGbx2uiU9++wzGDz4tiE64SfD6Rnak6IPF+XjFz+myMNF7I0hKgoqV4aV\nK5m6oyLDh5sW19Wr27usiIhIWrN9DHxq+fn53dF2cPv27RQpUuRm28FWrVoREhLCu+++S7vf+xFH\nRkayePFiYmNjqVWrlt1hinifzp0hb15o0ICfpn5Et8Rl7P5hNxNenEBA8QD713e5oG1bEt/qT78v\nKhIWZpqWPPWU/UuLiIh4G0fmwZ07d45ChQrd/LhgwYJs27bttufUrFmTmjVrJnmtoKCgm7/39/fH\n39/fXWGKeLzY+nX5cthL1O/QgXY9m/Kf4YfIel/WtFl8+HBcWf7BSxve4lqcKTPPlSttlhYREbFb\nREQEERERbrueI0m3O2tL/5h0i6QnK46toPeq3pTNW5a6a8Jp2qIT5PrErUN07ioigoRpM3jh0d0U\n9cvIpElw//32LikiIpKW/ryZO2zYsFRdz5Gku0CBApw5c+bmx2fOnKFgwYJOhCLidU79coo+q/pw\n5H9HmNRgEi8WedF8YvNmM0Tn/HkYN86+U4yXLhHXsh2dEj6nUefH6NvX/hxfRETE2znSW6BixYoc\nP36cqKgo4uPjWbBgAY0bN07RtYKCgty69S/iqa5dv8bQ9UOpNKMSVQpWYX/X/bcSbjBDdL75Bnbv\nTvEQnSRZFufqd+bz6Ba0mFWffv2UcIuIiG+LiIhwS2WF7d1LWrduTWRkJJcuXSJPnjwMHz6cDh06\nEBoaSp8+fUhISKBTp04MHDgw2ddW9xJJDyzLIuRoCH3D+lIpfyXGvjCWQg8WuvsLYmNNk+zLl2HJ\nEsiZ001xQFjTqRQInUHixi2Uq5zFLdcVERHxBqnNOz1+OM7fUdItvu7YpWP0XtWb05dPM7H+RGo/\nWfveXvjHITorV0K+fKmKIy4ORrTYz1srnyd+3Sby+hVL1fVERES8jccPx7GbykvEF8XExzBo7SCq\nzqxK7cK12fPmnntPuOHWEJ3mzaFatb8conOv/vc/aPj8NbpEtCbbpDFKuEVEJF3xmvISO2mnW3yN\nZVksOrSIt8Lfwu+ffoypO4b8OfIn/cK/M3MmvPcehISYQTbJcPQoNGwI83J247nil8kQ/B8VcYuI\nSLrk8cNxROTeHL54mJ6hPbkQc4Evmn1BjX/WcM+FO3WCPHlM9jx3LtSvf08vW7cOWreG4FeWUGXl\nKpi2Wwm3iIhICnl9eYmIt4uOi6Z/eH9qzK5B4+KN2fXGLvcl3DcEBJid7g4dTOKdhM8+Mwn3lAhe\nAwAAEutJREFU0olnqP3fNyE4GB580L0xiYiIpCNev9MdFBSkSZTilSzLYv6B+QxYPYA6T9bhQNcD\n5M2e174Fq1aF9etNL+8ffoABA+7YuU5IgHffNfn5NxEJFHvzVejTB557zr64REREPJi7JlOqplvE\nAft/2k+P0B5Ex0UzqcEkqhaqmnaLnztnEu/atW8bovPzz9Cxo+k0uHgxPDx5BEREQHi4OZgpIiKS\njqX77iUi3uRy7GV6h/am9tzatCzdkh2v70jbhBtuH6LTujVWbBwLF0KZMlC4sMmxHz68yXQ/mTdP\nCbeIiIgbKOkWSQMJiQlM/3Y6JSaV4Or1qxzsdpBulbqRKaNDCW2uXBAWxtVoF/sKNWDs0CssXgzj\nx0PmmF/McJ0ZMyB/KjuniIiICOADNd0ini4yKpLeq3qTI0sOVrZdSYV8FZwOicREmD47K0O3L+Tr\nwj3ZEl+DjP8MBesx6NIFGjc2hy9FRETELbw+6dZBSvFUpy+fpv/q/mw9u5UxdcfQonQLMnhAy72j\nR+H118HlgvUbMlGq5GR4/30zRKdtWzNIZ948p8MUERHxCDpIiQ5SimeKiY/hw00fMnnHZHpV7kX/\nav3Jdn82p8Pi+nUYM8aUkAwZAt26/alce+ZMePtt2LwZSpZ0LE4RERFPlNq8U0m3iJvcaAH4zpp3\n8Hvcjw/rfEihBws5HRYAO3dC586mRHvaNHj88bs8MSFBBydFRET+giZSiniAned30ntVb2Jdscx/\neT7VH6/udEgAXL1qdrW/+ALGjoU2bZIYKqmEW0RExBZKukVS4cfffmTQ2kGEnghl1POjaP9MezJm\n8IymQGvXmjORVarA/v2QO7fTEYmIiKRfXp906yClOCHOFceEbRMYvWk0Hct35GiPo+TMktPpsAD4\n5Rd46y2TdE+dCg0aOB2RiIiI99JBSlTTLWnPsiyWHV3GW+FvUSp3Kca+MJaijxR1OiwALAu++gp6\n9YLmzWHUKMiRw+moREREfINqukXSyMELB+kT1odzV84xpeEUXnjqBadDuuncOeje3XT7W7QIqqbx\nkEsRERH5e55RfCriwX6+9jM9V/ak1pxaNC7WmL1v7vWYhDsxEaZPh/Ll4ZlnzGR3JdwiIiKeRzvd\nInfhSnTx6c5PGRY5jFdKvcKh7od4NNujTod107Fj5qBkbCysWwdlyjgdkYiIiNyNkm6Rv7D21Fp6\nr+pNnn/kYW3gWsrmLet0SDddv27a/330EQweDD16qNOfiIiIp1PSLfIHJ38+ydur32bvj3sZ+8JY\nmpZo6hGj22/YtQs6dYI8eczAmyeecDoiERERuRdeX9MdFBTkljYukr5Fx0UzcM1A/vXZv6icvzKH\nuh/ipZIveUzCffUqDBgA9etDv36wapUSbhERkbQQERFBUFBQqq+jloGSriVaiczdO5f/W/d/1Hmy\nDv+u/W/y58jvdFi3WbfO1G5Xrgwff2x2uUVERCRtqWWgSAptObOFXqt6kSlDJha3WMy/Cv7L6ZBu\n88sv0L8/hIfDlCnQqJHTEYmIiEhKeX15iUhynb1yllcXv0rz/zanV+VebO602eMS7q++Mt1IsmaF\nAweUcIuIiHg77XRLunHt+jXGbhnL+K3j6VqxK0d7HCV75uxOh3Wb8+dNN5LDh2HhQqhWzemIRERE\nxB200y0+z7IsFh1aRMnJJdnz4x52vr6Tkc+P9KiE27JgxgwoVw5KlzZDbpRwi4iI+A7tdItP2/Pj\nHvqs6sPP135mVpNZ1Cpcy+mQ7nDiBLz+OsTEwNq18PTTTkckIiIi7qadbvFJF2Mu8sbyN6j3RT1a\nlWnFrjd2eVzC7XLB6NHw3HPQuDFs2aKEW0RExFd5/U53UFAQ/v7++Pv7Ox2KeID4hHgmb5/M+xvf\np23ZthzpfoSHHnjI6bDusHu3GXLz6KOwYwcULux0RCIiIvJXIiIi3DITRn26xWeEHg+lb1hfnsj1\nBOPrjadk7pJOh3SH6GgYNQpmzTK73IGB4CHzd0RERORvqE+3pHtH/3eUfuH9OH7pOOPrjadB0QYe\nM0kSIC7OTJAMDjb/DQiAffsgb16nIxMREZG0op1u8VqXYy8zInIEc/bOYWD1gfT8V08yZ8rsdFgA\nJCTAhg0m0V682PTcbtMGmjeHRx5xOjoRERFJLu10S7qTkJjA57s/Z/D6wQQUC+Bgt4Pkze78trFl\nwa5dJtH+8kszrr1NG9izBwoVcjo6ERERcZKSbvEqG05voPeq3mTPnJ2VbVdSIV8Fp0Pi+HGTaAcH\nm44kbdrA6tVQqpTTkYmIiIinUNItXuH05dMMWDOArWe3MrrOaFqUbuFo3fb587BggUm0z5yBli1h\n7lyoXFkHI0VEROROqukWjxYTH8PoTaOZtGMSvSr3on+1/mS7P5sjsVy+DF99ZRLtXbugaVOzq12r\nFtynt68iIiI+TTXd4pMsy2L+gfm8u+Zdqj1ejd1v7ObxBx9P8ziuXYOvvzaJ9rp1UKcOdOsGDRrA\nAw+keTgiIiLipZR0i8fZeX4nvVf1JtYVS/DLwVR/vHqaru9ymQQ7OBhCQuDZZ6FtW9NbO1euNA1F\nREREfITKS8Rj/PjbjwxaO4jQE6GMen4Ur5V7jUwZM6XJ2pYF27aZRHvhQnj8cVM60rIl5MuXJiGI\niIiIB1N5iXi9OFccE7ZNYPSm0XQs35GjPY6SM0vONFn70KFbnUcyZzaJ9jffQNGiabK8iIiIpBNe\nn3QHBQXh7++Pv7+/06FIMlmWxfJjy+kX1o9SuUuxpdMWij5if7b7/femj3ZwMFy8CK1bw6JFUL68\nOo+IiIjI7SIiIoiIiEj1dVReIo44eOEgfcP6cvbKWcbXG0+9IvVsXe9//zOJdXAwHDwIL79sdrX9\n/CBT2lSwiIiIiBdLbd6ppFvS1M/Xfmbo+qF8efBLBtcYTNeKXbk/0/22rPXbb7BsmUm0v/kG6tc3\niXa9epAliy1LioiIiI9STbd4BVeii093fsqwyGE0L9Wcw90P82i2R92+Tnw8hIebRHvlSqha1STa\n8+dDjhxuX05ERETknijpFtttOL2Bbiu6kecfeVgTuIan8z7t1usnJsLGjSbRXrQISpQwifaECZA7\nt1uXEhEREUkRlZeI7UKOhJBgJfBSiZfcNrrdsmDvXpNoz58PDz1kEu1WreCJJ9yyhIiIiMhNqun2\n3vAlBU6eNEl2cDBcvWoS7datoWxZpyMTERERX6ak23vDl3v000+wYIFJtE+dghYtTLJdpYpa/ImI\niEjaUNLtveHL37hyBZYsMYn2tm0QEGBGsdeuDffb0+xERERE5K6UdHtv+PInsbEQGmoS7fBw8Pc3\nO9oBAZAtm9PRiYiISHqmpNt7wxcgIQEiIkyivWQJlCtnEu2XX4aHH3Y6OhERERFDSbf3hp9uWRbs\n3GkS7QULIF8+k2i3bAkFCzodnYiIiMidNBxHvMbRoybRDg42H7dpA+vWmb7aIiIiIr5MSbfYbvFi\nGDUKzp83fbSDg6FiRXUeERERkfRD5SViu02bzCFJf3/IlMnpaERERESSL7V5Z0Y3xuJ2MTExVKpU\niRUrVjgdiqRCtWqm1Z8SbhEREUmvPDrpHj16NC1btnQ6DBH5CxEREU6HIJIu6d4T8U62J90dO3Yk\nb968lP3TnO5Vq1ZRokQJihYtyocffnjH61avXk2pUqXInTu33SGKSAroH34RZ+jeE/FOtifdHTp0\nYNWqVbc9lpCQQI8ePVi1ahWHDh1i/vz5HD58mHnz5tG3b1/Onz9PZGQkW7duJTg4mBkzZqh2+x54\n8l/EaR2bHeu565qpuU5KXpuc13jynyFP5snfN1+499x1XU++91K6Rnrnyd8z3XvuuY4v3Xu2J91+\nfn489NBDtz22fft2ihQpwhNPPMH9999Pq1atCAkJoV27dowfP578+fMzcuRIxo8fT5s2bejSpQsZ\n1OoiSfrLx971nP6LJ6WvVdJtP0/+vvnCveeu63ryvZfSNdI7T/6e6d5zz3V86d5Lk+4lUVFRBAQE\nsH//fgAWLVpEWFgYM2bMAOCLL75g27ZtTJw4MVnXVSIuIiIiImnF64bjuCtZVsmJiIiIiHgDR7qX\nFChQgDNnztz8+MyZMxTU/G8RERER8VGOJN0VK1bk+PHjREVFER8fz4IFC2jcuLEToYiIiIiI2M72\npLt169ZUrVqVY8eOUahQIWbNmsV9993HpEmTqFevHqVKlaJly5aULFnS7lBERERERBzh1WPgRURE\nRES8gUdPpEyukJAQunTpQqtWrVi9erXT4YikG0eOHKFr1660aNGCmTNnOh2OSLoSExNDpUqVWLFi\nhdOhiKQrERER+Pn50bVrVyIjI5N8vk8l3U2aNGH69OlMmzaNBQsWOB2OSLpRokQJpk6dypdffklY\nWJjT4YikK6NHj6Zly5ZOhyGS7mTMmJEcOXIQFxd3Tw1BPD7pTskY+ZEjR9KjR4+0DFPE5yT33lu+\nfDkNGzakVatWaR2qiE9Jzr23evVqSpUqRe7cuZ0IVcTnJOf+8/PzY+XKlXzwwQcMHTo0yWt7fNKd\nnDHylmXxzjvvUL9+fZ555hmHIhbxDcm59wACAgIIDQ1lzpw5ToQr4jOSc+9FRkaydetWgoODmTFj\nhuZXiKRScu6/G3NncuXKRVxcXJLXdmQ4TnL4+fkRFRV122N/HCMP3Bwjv2bNGtauXcuVK1c4ceIE\nb7zxRtoHLOIjknPvXbhwgcWLFxMbG0utWrXSPlgRH5Kce2/kyJEAzJkzh9y5c2tSs0gqJef+O3Lk\nCGFhYVy+fJmePXsmeW2PT7r/yrlz5yhUqNDNjwsWLHhzjPy9fNEikjJ3u/dq1qxJzZo1HYxMxLfd\n7d674bXXXnMiLJF04W7337vvvstLL710z9fx+PKSv6J38iLO0L0n4gzdeyLOcdf955VJt8bIizhD\n956IM3TviTjHXfefVybdGiMv4gzdeyLO0L0n4hx33X8en3RrjLyIM3TviThD956Ic+y8/zQGXkRE\nRETEZh6/0y0iIiIi4u2UdIuIiIiI2ExJt4iIiIiIzZR0i4iIiIjYTEm3iIiIiIjNlHSLiIiIiNhM\nSbeIiIiIiM2UdIuI+IC9e/cSGhp6189/++239O7dO1VrnD9/nldeeSVV1xARSa80HEdExAfMnj2b\nb7/9lokTJ97xOZfLxX333edAVCIicoN2ukVEPEBUVBQlSpSgQ4cOFC9enLZt2xIeHk61atUoVqwY\nO3bsAGD79u1UrVqVChUqUK1aNY4dO0Z8fDxDhgxhwYIFlC9fnoULFxIUFES7du2oXr06gYGBREZG\nEhAQAECfPn0YMWIEAGFhYdSsWfOOeCIjIylfvjzly5enQoUKxMTEEBUVRdmyZQHo3Lnzzc/nyZPn\n5vXGjBlD5cqVKVeuHEFBQWnwnRMR8RKWiIg47rvvvrPuu+8+68CBA1ZiYqL17LPPWh07drQsy7JC\nQkKspk2bWpZlWVeuXLFcLpdlWZa1evVq6+WXX7Ysy7Jmz55t9ezZ8+b1hg4dalWsWNGKjY21LMuy\n1q9fbzVq1MiyLMu6evWqVbp0aWvdunVW8eLFrVOnTt0RT0BAgLV582bLsiwrJibGcrlc1nfffWeV\nKVPmtudFRUVZpUqVsr7//nsrLCzM6tKli2VZlpWQkGA1atTI2rBhg9u+RyIi3kw/bxQR8RCFCxem\ndOnSAJQuXZo6deoAUKZMGaKiogC4fPkygYGBnDhxggwZMuByuQCwLAvrD9WCGTJkoHHjxmTJkuWO\ndR544AFmzJiBn58fEyZMoHDhwnc8p1q1avTt25e2bdvSrFkzChQocMdzYmNjeeWVV5g4cSKFChVi\nwoQJhIeHU758eQBiYmI4ceIEfn5+qfvGiIj4ACXdIiIe4o8JcsaMGcmcOfPN399IrgcPHkzt2rVZ\nsmQJp0+fxt/f/67Xy5Yt210/t2/fPnLnzs25c+f+8vPvvPMOjRo1YsWKFVSrVo2wsLA7Evg333yT\n5s2b8/zzz998bODAgXTp0iXJr1VEJL1RTbeIiBe5cuUK+fPnB2DWrFk3H8+ZMyfR0dH3dI3Tp08z\nbtw4du/eTWhoKNu3b7/jOSdPnqR06dIMGDCASpUqcfTo0ds+P3nyZH777TcGDBhw87F69erx+eef\nExMTA8C5c+e4ePFisr9GERFfpKRbRMRDZMiQ4a4f3/j9gAEDGDhwIBUqVCAhIeHm47Vq1eLQoUM3\nD1L+1etvfNy5c2fGjh3LY489xsyZM+ncuTPx8fG3rT1hwgTKli1LuXLlyJw5M/Xr17/tmmPHjuXA\ngQM3D1NOnz6dunXr0qZNG6pUqcLTTz9NixYt+O2339z5LRIR8VpqGSgiIiIiYjPtdIuIiIiI2ExJ\nt4iIiIiIzZR0i4iIiIjYTEm3iIiIiIjNlHSLiIiIiNhMSbeIiIiIiM2UdIuIiIiI2ExJt4iIiIiI\nzf4fkxKyHOM9bYoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7ad90badd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 6))\n",
    "\n",
    "ax.loglog(N_vec, duration_ref, label='numpy')\n",
    "ax.loglog(N_vec, duration_cy, label='cython')\n",
    "ax.loglog(N_vec, duration_cy_omp, label='cython+openmp')\n",
    "\n",
    "ax.legend(loc=2)\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_ylabel(\"matrix-vector multiplication duration\")\n",
    "ax.set_xlabel(\"matrix size\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For large problem sizes the the cython+OpenMP implementation is faster than numpy.dot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this simple implementation, the speedup for large problem sizes is about:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0072232987815148"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((duration_ref / duration_cy_omp)[-10:]).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Obviously one could do a better job with more effort, since the theoretical limit of the speed-up is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N_core"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * http://openmp.org\n",
    " * http://docs.cython.org/src/userguide/parallelism.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OpenCL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OpenCL is an API for heterogenous computing, for example using GPUs for numerical computations. There is a python package called `pyopencl` that allows OpenCL code to be compiled, loaded and executed on the compute units completely from within Python. This is a nice way to work with OpenCL, because the time-consuming computations should be done on the compute units in compiled code, and in this Python only server as a control language."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting opencl-dense-mv.py\n"
     ]
    }
   ],
   "source": [
    "%%file opencl-dense-mv.py\n",
    "\n",
    "import pyopencl as cl\n",
    "import numpy\n",
    "import time\n",
    "\n",
    "# problem size\n",
    "n = 10000\n",
    "\n",
    "# platform\n",
    "platform_list = cl.get_platforms()\n",
    "platform = platform_list[0]\n",
    "\n",
    "# device\n",
    "device_list = platform.get_devices()\n",
    "device = device_list[0]\n",
    "\n",
    "if False:\n",
    "    print((\"Platform name:\" + platform.name))\n",
    "    print((\"Platform version:\" + platform.version))\n",
    "    print((\"Device name:\" + device.name))\n",
    "    print((\"Device type:\" + cl.device_type.to_string(device.type)))\n",
    "    print((\"Device memory: \" + str(device.global_mem_size//1024//1024) + ' MB'))\n",
    "    print((\"Device max clock speed:\" + str(device.max_clock_frequency) + ' MHz'))\n",
    "    print((\"Device compute units:\" + str(device.max_compute_units)))\n",
    "\n",
    "# context\n",
    "ctx = cl.Context([device]) # or we can use cl.create_some_context()\n",
    "\n",
    "# command queue\n",
    "queue = cl.CommandQueue(ctx)\n",
    "\n",
    "# kernel\n",
    "KERNEL_CODE = \"\"\"\n",
    "//\n",
    "// Matrix-vector multiplication: r = m * v\n",
    "//\n",
    "#define N %(mat_size)d\n",
    "__kernel\n",
    "void dmv_cl(__global float *m, __global float *v, __global float *r)\n",
    "{\n",
    "    int i, gid = get_global_id(0);\n",
    "    \n",
    "    r[gid] = 0;\n",
    "    for (i = 0; i < N; i++)\n",
    "    {\n",
    "        r[gid] += m[gid * N + i] * v[i];\n",
    "    }\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "kernel_params = {\"mat_size\": n}\n",
    "program = cl.Program(ctx, KERNEL_CODE % kernel_params).build()\n",
    "\n",
    "# data\n",
    "A = numpy.random.rand(n, n)\n",
    "x = numpy.random.rand(n, 1)\n",
    "\n",
    "# host buffers\n",
    "h_y = numpy.empty(numpy.shape(x)).astype(numpy.float32)\n",
    "h_A = numpy.real(A).astype(numpy.float32)\n",
    "h_x = numpy.real(x).astype(numpy.float32)\n",
    "\n",
    "# device buffers\n",
    "mf = cl.mem_flags\n",
    "d_A_buf = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=h_A)\n",
    "d_x_buf = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=h_x)\n",
    "d_y_buf = cl.Buffer(ctx, mf.WRITE_ONLY, size=h_y.nbytes)\n",
    "\n",
    "# execute OpenCL code\n",
    "t0 = time.time()\n",
    "event = program.dmv_cl(queue, h_y.shape, None, d_A_buf, d_x_buf, d_y_buf)\n",
    "event.wait()\n",
    "cl.enqueue_copy(queue, h_y, d_y_buf)\n",
    "t1 = time.time()\n",
    "\n",
    "print(\"opencl elapsed time =\", (t1-t0))\n",
    "\n",
    "# Same calculation with numpy\n",
    "t0 = time.time()\n",
    "y = numpy.dot(h_A, h_x)\n",
    "t1 = time.time()\n",
    "\n",
    "print(\"numpy elapsed time =\", (t1-t0))\n",
    "\n",
    "# see if the results are the same\n",
    "print(\"max deviation =\", numpy.abs(y-h_y).max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/pyopencl-2012.1-py2.7-linux-x86_64.egg/pyopencl/__init__.py:36: CompilerWarning: Non-empty compiler output encountered. Set the environment variable PYOPENCL_COMPILER_OUTPUT=1 to see more.\n",
      "  \"to see more.\", CompilerWarning)\n",
      "opencl elapsed time = 0.0188570022583\n",
      "numpy elapsed time = 0.0755031108856\n",
      "max deviation = 0.0136719\n"
     ]
    }
   ],
   "source": [
    "!python opencl-dense-mv.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* http://mathema.tician.de/software/pyopencl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/json": {
       "Software versions": [
        {
         "module": "Python",
         "version": "3.3.2+ (default, Oct  9 2013, 14:50:09) [GCC 4.8.1]"
        },
        {
         "module": "IPython",
         "version": "2.0.0-b1"
        },
        {
         "module": "OS",
         "version": "posix [linux]"
        },
        {
         "module": "numpy",
         "version": "1.9.0.dev-d4c7c3a"
        },
        {
         "module": "mpi4py",
         "version": "1.3.1"
        },
        {
         "module": "Cython",
         "version": "0.20.post0"
        }
       ]
      },
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>3.3.2+ (default, Oct  9 2013, 14:50:09) [GCC 4.8.1]</td></tr><tr><td>IPython</td><td>2.0.0-b1</td></tr><tr><td>OS</td><td>posix [linux]</td></tr><tr><td>numpy</td><td>1.9.0.dev-d4c7c3a</td></tr><tr><td>mpi4py</td><td>1.3.1</td></tr><tr><td>Cython</td><td>0.20.post0</td></tr><tr><td colspan='2'>Mon Mar 17 17:32:10 2014 JST</td></tr></table>"
      ],
      "text/latex": [
       "\\begin{tabular}{|l|l|}\\hline\n",
       "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
       "Python & 3.3.2+ (default, Oct  9 2013, 14:50:09) [GCC 4.8.1] \\\\ \\hline\n",
       "IPython & 2.0.0-b1 \\\\ \\hline\n",
       "OS & posix [linux] \\\\ \\hline\n",
       "numpy & 1.9.0.dev-d4c7c3a \\\\ \\hline\n",
       "mpi4py & 1.3.1 \\\\ \\hline\n",
       "Cython & 0.20.post0 \\\\ \\hline\n",
       "\\hline \\multicolumn{2}{|l|}{Mon Mar 17 17:32:10 2014 JST} \\\\ \\hline\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "Software versions\n",
       "Python 3.3.2+ (default, Oct  9 2013, 14:50:09) [GCC 4.8.1]\n",
       "IPython 2.0.0-b1\n",
       "OS posix [linux]\n",
       "numpy 1.9.0.dev-d4c7c3a\n",
       "mpi4py 1.3.1\n",
       "Cython 0.20.post0\n",
       "\n",
       "Mon Mar 17 17:32:10 2014 JST"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%load_ext version_information\n",
    "\n",
    "%version_information numpy, mpi4py, Cython"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}