Lecture-6B-HPC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 6B - 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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import multiprocessing\n",
    "import os\n",
    "import time\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pool = multiprocessing.Pool(processes=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from IPython.parallel import Client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cli = Client()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the 'ids' attribute we can retrieve a list of ids for the IPython engines in the cluster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dview = cli[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "delay_times = numpy.random.rand(64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAADSCAYAAAAYCPc3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAHfdJREFUeJzt3X9wFPX9x/HXHQFiSpofChcN1jtQhIRL7vhhxoI1MQREKqLlh4DkKvSHdSaMjkiLdNra4RvCSCmpaB07FNLaAZk6jVEq7aVNwEgtmh8HHYNYyg2gScaAIQkEI7n7/sF4YwgJv3K32eT5mMmMt7f72ff7czfh5W521xIMBoMCAABAn2Y1ugAAAABcGqENAADABAhtAAAAJkBoAwAAMAFCGwAAgAlEGV1AuFksFqNLAAAAuGzd3dhjQBxpCwaDA/bn5z//ueE1GP0z0OeA/unf6Bron/7p//J/ejIgQhsAAIDZEdoAAABMgNDWz2VmZhpdguEG+hzQf6bRJRiK/jONLsFQ9J9pdAm9yhK81AlUk7NYLJc8RwwAANAX9JRbONIGAABgAoQ2AAAAEyC0AQAAmAChDQAAwAQIbQAAACbQ7x9jJUlP7HpCTWebIrY/X71P6Unpphk3kvpDDz0J92dv1PyZ+XMzc+2Xcrm9RWIOrmYfvnqfJBn++fAdMVZ8dLw23rvR6DJMYUCEtqazTbLH2yO2v4qjFWHZX7jGjaT+0ENPwv3ZGzV/Zv7czFz7pVxub5GYg6vZR8XRCkky/PPhO2Isf5Pf6BJMg9OjAAAAJkBoAwAAMAFCGwAAgAkQ2gAAAEwgbKHt7NmzysjIkMvlUkpKilatWiVJOnnypHJycjRmzBhNnz5dTU1NoeVZWVmKjY1VXl5ep7G2bNkip9Op9PR0zZw5UydOnJAkbd26VcOHD5fb7Zbb7dbvf//7cLUDAABgqLCFtujoaJWVlammpkb79+9XWVmZKioqVFBQoJycHB06dEjZ2dkqKCgIrb9mzRqtX7++0zjt7e1asWKFdu/eLZ/Pp7S0NG3atEnS+YeqLly4UNXV1aqurtbSpUvD1Q4AAIChwnp6NCYmRtL54NXR0aGEhASVlJTI4/FIkjwej4qLi0PrTpkyRUOHDu00RlRUlBISEtTa2qpgMKhTp04pOTlZkhQMBhUMBsPZAgAAQJ8Q1tAWCATkcrlks9mUlZWl1NRUNTQ0yGazSZJsNpsaGho6bWOxWDoXaLWqsLBQ48ePV3Jysg4ePBg6omaxWPTaa68pLS1N8+bN0/Hjx8PZDgAAgGHCenNdq9WqmpoanTp1SjNmzFBZWVmn9y0WS5eQdqHm5mYtX75cPp9PDodDeXl5Wrt2rVavXq37779fixYt0uDBg/Xyyy/L4/HoH//4R5cxarbXyB/tlyTZXXbZXfbeahEAAOCqlZeXq7y8/LLWjcgTEeLi4jRr1ixVVlbKZrOpvr5eSUlJqqur04gRI3rctra2Vg6HQw6HQ5I0b948rVu3TpKUmJgYWm/ZsmVauXLlRcdwPezq83eEBgAAA09mZqYyMzNDr5999tlu1w3b6dHGxsbQlaFtbW3yer1yu92aPXu2ioqKJElFRUWaM2dOp+0u/Bu1UaNG6eDBg2psbJQkeb1epaSkSJLq6upC65WUlISWAwAA9DdhO9JWV1cnj8ejQCCgQCCgJUuWKDs7W263W/Pnz9fmzZtlt9u1Y8eO0DZ2u10tLS1qb29XcXGxvF6vxo4dq/z8fGVlZclqtcput2vr1q2SpOeff14lJSWKiorS9ddfH1oOAADQ34QttDmdTlVVVXVZnpiYqNLS0otu4/f7L7o8NzdXubm5XZbn5+crPz//muoEAAAwA56IAAAAYAKENgAAABMgtAEAAJgAoQ0AAMAECG0AAAAmQGgDAAAwgYg8EcFo8dHx8jf5I7a/2CGxYdlfuMaNpP7QQ0/C/dkbNX9m/tzMXPulXG5vkZiDq9lH7JBYSTL88+E7Yqz46HijSzANS/DCRxD0MxaLpctTFgAAAPqinnILp0cBAABMgNAGAABgAoQ2AAAAEyC0AQAAmAChDQAAwAQIbQAAACZAaAMAADABQhsAAIAJENoAAABMgNAGAABgAoQ2AAAAEyC0AQAAmAChDQAAwAQIbQAAACZAaAMAADABQhsAAIAJENoAAABMgNAGAABgAoQ2AAAAE4gyuoBIeGLXE2o629Tt+756n9KT0iNYUfhcay9Xun1vz91Xx+tp7L74mYWjpqsd0+zfg6tRcbRCU78x1dAaelN3cxofHa+N9240oCIARhsQoa3pbJPs8fZu3684WtHj+2Zyrb1c6fa9PXdfHa+nsfviZxaOmq52TLN/D65GyYclhtfQm7qbU3+TP+K1AOgbOD0KAABgAoQ2AAAAEyC0AQAAmAChDQAAwATCFtrOnj2rjIwMuVwupaSkaNWqVZKkkydPKicnR2PGjNH06dPV1NQUWp6VlaXY2Fjl5eV1GmvLli1yOp1KT0/XzJkzdeLECUnSSy+9pLS0NLndbt15553y+XzhagcAAMBQYQtt0dHRKisrU01Njfbv36+ysjJVVFSooKBAOTk5OnTokLKzs1VQUBBaf82aNVq/fn2ncdrb27VixQrt3r1bPp9PaWlp2rRpkyRp8eLF2r9/v6qrq/XMM8/oqaeeClc7AAAAhgrr6dGYmBhJ54NXR0eHEhISVFJSIo/HI0nyeDwqLi4OrTtlyhQNHTq00xhRUVFKSEhQa2urgsGgmpublZycLEmKjY0Nrdfa2qobbrghnO0AAAAYJqz3aQsEApowYYIOHz6sH/3oR0pNTVVDQ4NsNpskyWazqaGhodM2Foul02ur1arCwkKNHz9ew4YN05gxY/TCCy+E3n/xxRe1YcMGnT59Wnv37g1nOwAAAIYJa2izWq2qqanRqVOnNGPGDJWVlXV632KxdAlpF2pubtby5cvl8/nkcDiUl5entWvXavXq1ZKkxx9/XI8//ri2bdumpUuXdtmHJNVsr5E/2i9JsrvssrvsvdIfAADAtSgvL1d5efllrRuRJyLExcVp1qxZqqyslM1mU319vZKSklRXV6cRI0b0uG1tba0cDoccDockad68eVq3bl2X9RYsWKDHHnvsomO4Hnb1qzulAwCA/iEzM1OZmZmh188++2y364btb9oaGxtDV4a2tbXJ6/XK7XZr9uzZKioqkiQVFRVpzpw5nbYLBoOdXo8aNUoHDx5UY2OjJMnr9SolJUWS9NFHH4XW27lzp9LS0sLVDgAAgKHCdqStrq5OHo9HgUBAgUBAS5YsUXZ2ttxut+bPn6/NmzfLbrdrx44doW3sdrtaWlrU3t6u4uJieb1ejR07Vvn5+crKypLVapXdbtfWrVslSS+88IJKS0s1ePBgDR8+XFu2bAlXOwAAAIYKW2hzOp2qqqrqsjwxMVGlpaUX3cbv9190eW5urnJzc7ss37hx4zXVCAAAYBY8EQEAAMAECG0AAAAmQGgDAAAwAUIbAACACRDaAAAATIDQBgAAYAIReSKC0eKj4+Vv8nf7fuyQ2B7fN5Nr7eVKt+/tufvqeD2N3Rc/s3DUdLVjmv17cDUSr0s0vIbe1N2cxkfHR74YAH2CJXjhIwj6GYvF0uUpCwAAAH1RT7mF06MAAAAmQGgDAAAwgR5D29atWzVhwgTFxMQoJiZGkyZNCj3sHQAAAJHT7YUIRUVFKiws1IYNG+R2uxUMBlVdXa2nn35aFovlos8CBQAAQHh0eyFCRkaGtm/fLofD0Wm53+/XggUL9O9//zsiBV4rLkQAAABmcVUXIrS0tHQJbJJkt9vV0tLSe9UBAADgkroNbdHR0d1u1NN7AAAA6H3dnh697rrrdOutt150o8OHD+vMmTNhLay3cHoUAACYRU+5pdsLEWpra8NWEAAAAK4MT0QAAADoI67qSNuwYcNksVi6HbC5ubl3qgMAAMAlcaQNAACgj7iqI21tbW166aWXdPjwYTmdTi1btkxRUd2uDgAAgDDq9kjb/PnzNWTIEE2dOlVvvfWW7Ha7CgsLI13fNeNIGwAAMIuecku3oc3pdOrAgQOSpHPnzmny5Mmqrq4OX5VhQmgDAABmcVVPRPjqqVBOiwIAABir2yNtgwYNUkxMTOh1W1ubrrvuuvMbmejqUY60AQAAs7iqCxE6OjrCVhAAAACuTLenRwEAANB3DIg/Vvtu8XeNLuGifPU+pSelX/M6ZtZTf5fq3YxzY8aau3OtvUR6LiKxvyvZR1+cv778/TSitr48H1+Kj47Xxns3Gl0GImRAhDZ7vN3oEi6q4mjFJWu7nHXMrKf+LtW7GefGjDV351p7ifRcRGJ/V7KPvjh/ffn7aURtfXk+vuRv8htdAiKI06MAAAAmQGgDAAAwAUIbAACACRDaAAAATCBsoe3s2bPKyMiQy+VSSkqKVq1aJUk6efKkcnJyNGbMGE2fPl1NTU2h5VlZWYqNjVVeXl6nsbZs2SKn06n09HTNnDlTJ06ckCRt2LBBqampSk9P17Rp03T06NFwtQMAAGCosIW26OholZWVqaamRvv371dZWZkqKipUUFCgnJwcHTp0SNnZ2SooKAitv2bNGq1fv77TOO3t7VqxYoV2794tn8+ntLQ0bdq0SZI0YcIEVVZWyufzae7cuVq5cmW42gEAADBUWE+PfvkYrPb2dnV0dCghIUElJSXyeDySJI/Ho+Li4tC6U6ZM0dChQzuNERUVpYSEBLW2tioYDKq5uVnJycmSpMzMTEVHR0uSMjIydPz48XC2AwAAYJiwhrZAICCXyyWbzaasrCylpqaqoaFBNptNkmSz2dTQ0NBpG4vF0rlAq1WFhYUaP368kpOTVVtbq6VLl3bZ1+bNm3XfffeFrxkAAAADhfXmularVTU1NTp16pRmzJihsrKyTu9bLJYuIe1Czc3NWr58uXw+nxwOh/Ly8rR27VqtXr06tM4rr7yiqqoq/frXv77oGOVby0P/bXfZZXfZr7onAACA3lJeXq7y8vLLWjciT0SIi4vTrFmzVFlZKZvNpvr6eiUlJamurk4jRozocdva2lo5HA45HA5J0rx587Ru3brQ+6WlpcrPz9eePXs0ePDgi46R+d3MXusFAACgt2RmZiozMzP0+tlnn+123bCdHm1sbAxdGdrW1iav1yu3263Zs2erqKhIklRUVKQ5c+Z02i4YDHZ6PWrUKB08eFCNjY2SJK/Xq5SUFElSdXW1HnvsMb3xxhu64YYbwtUKAACA4cJ2pK2urk4ej0eBQECBQEBLlixRdna23G635s+fr82bN8tut2vHjh2hbex2u1paWtTe3q7i4mJ5vV6NHTtW+fn5ysrKktVqld1u19atWyVJK1eu1OnTpzV37lxJ0i233BK6sAEAAKA/CVtoczqdqqqq6rI8MTFRpaWlF93G7/dfdHlubq5yc3O7LPd6vddUIwAAgFnwRAQAAAATILQBAACYAKENAADABAhtAAAAJkBoAwAAMIGI3FzXaP4mv9ElXFTskNhL1nY565hZT/1dqnczzo0Za+7OtfYS6bmIxP6uZB99cf768vfTiNr68nx8KT463ugSEEGW4IV3s+1nLBZLlxv2AgAA9EU95RZOjwIAAJgAoQ0AAMAECG0AAAAmQGgDAAAwAUIbAACACRDaAAAATIDQBgAAYAKENgAAABMgtAEAAJgAoQ0AAMAECG0AAAAmQGgDAAAwAUIbAACACRDaAAAATIDQBgAAYAKENgAAABMgtAEAAJgAoQ0AAMAECG0AAAAmEGV0AZHwxK4n1HS2yegyIs5X71N6UrrRZURMuPrtzXEvHMtMn9GV1mqm3i6lP/XypXD31Ffn7Frq6qs9XS6z1m/m35u9bUCEtqazTbLH240uI+IqjlYMqL7D1W9vjnvhWGb6jK60VjP1din9qZcvhbunvjpn11JXX+3pcpm1fjP/3uxtnB4FAAAwAUIbAACACRDaAAAATIDQBgAAYAKENgAAABMIW2g7e/asMjIy5HK5lJKSolWrVkmSTp48qZycHI0ZM0bTp09XU1NTaHlWVpZiY2OVl5fXaawtW7bI6XQqPT1dM2fO1IkTJyRJe/bs0YQJEzR48GC99tpr4WoFAADAcGELbdHR0SorK1NNTY3279+vsrIyVVRUqKCgQDk5OTp06JCys7NVUFAQWn/NmjVav359p3Ha29u1YsUK7d69Wz6fT2lpadq0aZMk6ZZbblFRUZEWLVoUrjYAAAD6hLCeHo2JiZF0Pnh1dHQoISFBJSUl8ng8kiSPx6Pi4uLQulOmTNHQoUM7jREVFaWEhAS1trYqGAyqublZycnJks6HNqfTKauVs7wAAKB/C2vaCQQCcrlcstlsysrKUmpqqhoaGmSz2SRJNptNDQ0NnbaxWCydC7RaVVhYqPHjxys5OVm1tbVaunRpOMsGAADoc8L6RASr1aqamhqdOnVKM2bMUFlZWaf3LRZLl5B2oebmZi1fvlw+n08Oh0N5eXlau3atVq9efdl11GyvkT/aL0myu+yyu+xX2goAAECv89f45a/xX9a6EXmMVVxcnGbNmqXKykrZbDbV19crKSlJdXV1GjFiRI/b1tbWyuFwyOFwSJLmzZundevWdVmvp/Dnetg1YB95AQAA+q4LDybtLtrd7bphOz3a2NgYujK0ra1NXq9Xbrdbs2fPVlFRkSSpqKhIc+bM6bRdMBjs9HrUqFE6ePCgGhsbJUler1cpKSldtrlwOwAAgP4kbEfa6urq5PF4FAgEFAgEtGTJEmVnZ8vtdmv+/PnavHmz7Ha7duzYEdrGbrerpaVF7e3tKi4ultfr1dixY5Wfn6+srCxZrVbZ7XZt3bpVkvTee+/poYce0meffaY333xTv/jFL3TgwIFwtQQAAGCYsIU2p9OpqqqqLssTExNVWlp60W38fv9Fl+fm5io3N7fL8smTJ+vYsWPXVCcAAIAZcK8MAAAAEyC0AQAAmAChDQAAwAQIbQAAACZAaAMAADCBiNxc12jx0fHyN/mNLiPiYofEDqi+w9Vvb4574Vhm+oyutFYz9XYp/amXL4W7p746Z9dSV1/t6XKZtX4z/97sbZZgP78rrcVi4ca7AADAFHrKLZweBQAAMAFCGwAAgAkQ2gAAAEyA0AYAAGAChDYAAAATILT1c+Xl5UaXYLiBPgf0X250CYai/3KjSzAU/ZcbXUKvIrT1c/3tC3s1Bvoc0H+50SUYiv7LjS7BUPRfbnQJvYrQBgAAYAKENgAAABMYEE9EAAAAMIvuolm/f/ZoP8+kAABggOD0KAAAgAkQ2gAAAEyA0AYAAGAC/Tq07dq1S2PHjtVtt92mdevWGV1ORC1dulQ2m01Op9PoUgxx7NgxZWVlKTU1VePHj9dvfvMbo0uKqLNnzyojI0Mul0spKSlatWqV0SUZoqOjQ263W/fff7/RpUSc3W5XWlqa3G637rjjDqPLibimpibNnTtX48aNU0pKit59912jS4qYDz/8UG63O/QTFxc34H4Hrl27VqmpqXI6nVq0aJE+//xzo0vqFf326tGOjg7dfvvtKi0tVXJysiZPnqxt27Zp3LhxRpcWEW+//baGDRum3NxcHThwwOhyIq6+vl719fVyuVxqbW3VxIkTVVxcPGA+f0k6c+aMYmJidO7cOU2dOlXr16/X1KlTjS4rojZs2KDKykq1tLSopKTE6HIiyuFwqLKyUomJiUaXYgiPx6O7775bS5cu1blz53T69GnFxcUZXVbEBQIBJScna9++fbr55puNLici/H6/7rnnHtXW1mro0KFasGCB7rvvPnk8HqNLu2b99kjbvn37dOutt8put2vw4MF6+OGH9frrrxtdVsTcddddSkhIMLoMwyQlJcnlckmShg0bpnHjxumTTz4xuKrIiomJkSS1t7ero6NjwP3jffz4cf31r3/V9773vQF7FflA7fvUqVN6++23tXTpUklSVFTUgAxsklRaWqrRo0cPmMAmSV//+tc1ePBgnTlzRufOndOZM2eUnJxsdFm9ot+Gto8//rjTl3TkyJH6+OOPDawIRvH7/aqurlZGRobRpURUIBCQy+WSzWZTVlaWUlJSjC4pop588kk999xzslr77a+5HlksFk2bNk2TJk3S7373O6PLiagjR45o+PDhevTRRzVhwgR9//vf15kzZ4wuyxDbt2/XokWLjC4johITE/XUU0/pG9/4hm666SbFx8dr2rRpRpfVK/rtbzNuqgtJam1t1dy5c1VYWKhhw4YZXU5EWa1W1dTU6Pjx49qzZ0+/ewZfT958802NGDFCbrd7wB5teuedd1RdXa233npLL7zwgt5++22jS4qYc+fOqaqqSo8//riqqqr0ta99TQUFBUaXFXHt7e164403NG/ePKNLiajDhw9r48aN8vv9+uSTT9Ta2qo//elPRpfVK/ptaEtOTtaxY8dCr48dO6aRI0caWBEi7YsvvtB3vvMdPfLII5ozZ47R5RgmLi5Os2bN0vvvv290KRGzd+9elZSUyOFwaOHChfrnP/+p3Nxco8uKqBtvvFGSNHz4cD344IPat2+fwRVFzsiRIzVy5EhNnjxZkjR37lxVVVUZXFXkvfXWW5o4caKGDx9udCkR9f777+ub3/ymrr/+ekVFRemhhx7S3r17jS6rV/Tb0DZp0iR99NFH8vv9am9v16uvvqrZs2cbXRYiJBgMatmyZUpJSdETTzxhdDkR19jYqKamJklSW1ubvF6v3G63wVVFTn5+vo4dO6YjR45o+/btuueee/SHP/zB6LIi5syZM2ppaZEknT59Wn//+98H1JXkSUlJuvnmm3Xo0CFJ5/+uKzU11eCqIm/btm1auHCh0WVE3NixY/Xuu++qra1NwWBQpaWl/ebPQ/rtY6yioqK0adMmzZgxQx0dHVq2bNmAunJw4cKF2r17t06cOKGbb75Zv/zlL/Xoo48aXVbEvPPOO3rllVdCtzyQzl8Cfu+99xpcWWTU1dXJ4/EoEAgoEAhoyZIlys7ONroswwy0P5doaGjQgw8+KOn8qcLFixdr+vTpBlcVWc8//7wWL16s9vZ2jR49Wlu2bDG6pIg6ffq0SktLB9zfM0pSenq6cnNzNWnSJFmtVk2YMEE/+MEPjC6rV/TbW34AAAD0J/329CgAAEB/QmgDAAAwAUIbAACACRDaAAAATIDQBgAAYAKENgAAABMgtAFAhJSXl+v+++83ugwAJkVoAwAAMAFCG4AB4/Tp05o1a5ZcLpecTqd27NihyspKZWZmatKkSbr33ntVX18vSfrvf/+radOmyeVyaeLEiTpy5Igk6emnn5bT6VRaWpp27Ngh6fwRtMzMTM2bN0/jxo3TI488Etrnrl27NG7cOE2cOFF/+ctfQst3794tt9stt9utCRMmqLW1NYIzAcCM+u1jrADgQrt27VJycrJ27twpSWpubtbMmTNVUlKi66+/Xq+++qpWr16tzZs3a/HixXrmmWf0wAMPqL29XR0dHXrttdfk8/m0f/9+ffrpp5o8ebK+9a1vSZJqamr0wQcf6MYbb9SUKVO0d+/e0ONzysrKNHr0aC1YsCD0SK1f/epXevHFF3XnnXfqzJkzGjp0qGHzAsAcONIGYMBIS0uT1+vVT37yE1VUVOjo0aP6z3/+o2nTpsntduv//u//9PHHH6u1tVWffPKJHnjgAUnSkCFDdN111+mdd97RokWLZLFYNGLECN1999167733ZLFYdMcdd+imm26SxWKRy+XSkSNHdPDgQTkcDo0ePVqS9Mgjj+jLJwdOmTJFTz75pJ5//nl99tlnGjRokGHzAsAcONIGYMC47bbbVF1drZ07d+qnP/2psrKylJqaqr1793Zar6WlpdsxLnxc85dHzr56pGzQoEE6d+5clwfVf3XbH//4x/r2t7+tnTt3asqUKfrb3/6m22+//ap7A9D/caQNwIBRV1en6OhoLV68WCtWrNC+ffvU2Niod999V5L0xRdf6IMPPlBsbKxGjhyp119/XZL0+eefq62tTXfddZdeffVVBQIBffrpp9qzZ4/uuOOOLkFOOh/mxo4dK7/fr//973+SpG3btoXeP3z4sFJTU7Vy5UpNnjxZH374YQRmAICZcaQNwIBx4MABPf3007JarRoyZIh++9vfatCgQVq+fLlOnTqlc+fO6cknn1RKSor++Mc/6oc//KF+9rOfafDgwfrzn/+sBx98UP/617+Unp4ui8Wi5557TiNGjFBtbW2Xo2rS+aNvL7/8smbNmqWYmBjdddddOn36tCSpsLBQZWVlslqtGj9+vGbOnBnp6QBgMpbgxf4XEQAAAH0Kp0cBAABMgNAGAABgAoQ2AAAAEyC0AQAAmAChDQAAwAQIbQAAACZAaAMAADCB/wdFi5aGCDL2JAAAAABJRU5ErkJggg==",
      "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lbview = cli.load_balanced_view()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAADSCAYAAAAYCPc3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAHYtJREFUeJzt3X9wFPX9x/HXHQFiSoYEhYsNjnegCAmX3PHDjIPW5BsCUipiCygiSYX+sM7A6BRpkU5bOzSEkVJT0TrtUJLWGZCp0xil0l6mCRipxebHQccglnIFNMkYMCSBpJHcff9wuGnIjwaSvb1Nno+ZzHh7++P9/tzmeLmb3bWFQqGQAAAAENXsZhcAAACA/43QBgAAYAGENgAAAAsgtAEAAFgAoQ0AAMACYswuwGg2m83sEgAAAAasrxt7jIgjbaFQiJ8B/vzoRz8yvQYr/TBejBfjFT0/jBfjNRzGqz8jIrQBAABYHaENAADAAght6CYzM9PsEiyF8bo2jNe1YbyuDeN1bRivaxMN42UL/a8TqBZns9n+5zliAACAaNBfbuFIGwAAgAUQ2gAAACyA0AYAAGABhDYAAAALILQBAABYwLB/jJUkPXngSTV3NBuy7oTYBEnqtn5/g1/pSenXvK7rXc7Kehu/4Wagn+tIGIuBMnMshsPnMJAeIvF9M5htjMTvw2js+VprGuzvTzSOQTQZEaGtuaNZzgSnIesONAckqdv6K09XXtf2rnc5K+tt/IabgX6uI2EsBsrMsRgOn8NAeojE981gtjESvw+jsedrrWmwvz/ROAbRhNOjAAAAFkBoAwAAsABCGwAAgAUQ2gAAACzAsNDW0dGhjIwMeTwepaSkaNOmTZKk8+fPKycnR9OmTdOCBQvU3Nwcnp6VlaX4+HitW7eu27p2794tt9ut9PR0LVq0SOfOnZMkFRUVaeLEifJ6vfJ6vfrNb35jVDsAAACmMiy0xcbGqry8XLW1tTp69KjKy8tVWVmpgoIC5eTk6MSJE8rOzlZBQUF4/i1btmj79u3d1tPZ2akNGzbo4MGD8vv9SktL086dOyV9/lDVlStXqqamRjU1NVqzZo1R7QAAAJjK0NOjcXFxkj4PXl1dXUpMTFRpaany8vIkSXl5eSopKQnPO2/ePI0dO7bbOmJiYpSYmKi2tjaFQiFduHBBycnJkqRQKKRQKGRkCwAAAFHB0NAWDAbl8XjkcDiUlZWl1NRUNTY2yuFwSJIcDocaGxu7LWOz2boXaLersLBQM2fOVHJyso4fPx4+omaz2fTaa68pLS1Ny5cv19mzZ41sBwAAwDSGhja73a7a2lqdPXtWhw4dUnl5ebf3bTZbj5B2tZaWFq1fv15+v18ff/yx3G63tm7dKkm6//779e9//1tHjx5VTk5O+Aje1Wr31qqiqEIVRRUK1AaGpDcAAIDBCtQGwhmloqii33kj8kSE8ePHa/HixaqqqpLD4VBDQ4OSkpJUX1+vSZMm9btsXV2dXC6XXC6XJGn58uXatm2bJGnChAnh+dauXauNGzf2ug7Pwx7usAwAAKKO0+OU0+MMvz5YfLDPeQ070tbU1BS+MrS9vV0+n09er1dLlixRcXGxJKm4uFhLly7tttzVf6M2ZcoUHT9+XE1NTZIkn8+nlJQUSVJ9fX14vtLS0vB0AACA4cawI2319fXKy8tTMBhUMBjU6tWrlZ2dLa/XqxUrVmjXrl1yOp3at29feBmn06nW1lZ1dnaqpKREPp9P06dPV35+vrKysmS32+V0OlVUVCRJeuGFF1RaWqqYmBjdeOON4ekAAADDjWGhze12q7q6usf0CRMmqKysrNdlAoFAr9Nzc3OVm5vbY3p+fr7y8/MHVScAAIAV8EQEAAAACyC0AQAAWAChDQAAwAIIbQAAABZAaAMAALAAQhsAAIAFROSJCGZLiE1QoDlg2LoldVt//Jj469re9S5nZb2N33Az0M91JIzFQJk5FsPhcxhID5H4vhnMNkbi92E09nytNQ329ycaxyCa2EJXP4JgmLHZbD2esgAAABCN+sstnB4FAACwAEIbAACABRDaAAAALIDQBgAAYAGENgAAAAsgtAEAAFgAoQ0AAMACCG0AAAAWQGgDAACwAEIbAACABRDaAAAALIDQBgAAYAGENgAAAAsgtAEAAFgAoQ0AAMACCG0AAAAWQGgDAACwAEIbAACABRDaAAAALCDG7AIi4ckDT6q5o9nsMqKSv8Gv9KR0s8swRUJsgiRZct8YyOc2Ej/bSPQ80P3GyvvX9RrqnqNlH45UHVe2Ey19G8Xf4JekqO4xITZBz9/3vNll9DAiQltzR7OcCU6zy4hKlacrR+zYBJoDkmTJ/gfyuY3EzzYSPQ90v7Hy/nW9hrrnaNmHI1XHle1ES99GqTxdKSm6fzeu7MvRhtOjAAAAFkBoAwAAsABCGwAAgAUQ2gAAACzAsNDW0dGhjIwMeTwepaSkaNOmTZKk8+fPKycnR9OmTdOCBQvU3Nwcnp6VlaX4+HitW7eu27p2794tt9ut9PR0LVq0SOfOnZMkvfzyy0pLS5PX69Vdd90lv99vVDsAAACmMiy0xcbGqry8XLW1tTp69KjKy8tVWVmpgoIC5eTk6MSJE8rOzlZBQUF4/i1btmj79u3d1tPZ2akNGzbo4MGD8vv9SktL086dOyVJq1at0tGjR1VTU6NnnnlG3/3ud41qBwAAwFSGnh6Ni4uT9Hnw6urqUmJiokpLS5WXlydJysvLU0lJSXjeefPmaezYsd3WERMTo8TERLW1tSkUCqmlpUXJycmSpPj4+PB8bW1tuummm4xsBwAAwDSG3qctGAxq1qxZOnnypL7zne8oNTVVjY2NcjgckiSHw6HGxsZuy9hstm6v7Xa7CgsLNXPmTI0bN07Tpk3Tiy++GH7/pZde0o4dO3Tx4kUdPnzYyHYAAABMY2hos9vtqq2t1YULF7Rw4UKVl5d3e99ms/UIaVdraWnR+vXr5ff75XK5tG7dOm3dulWbN2+WJD3xxBN64okntGfPHq1Zs6bHNiSpdm+tArEBSZLT45TT4xyS/gAAAAajoqJCFRUVA5o3Ik9EGD9+vBYvXqyqqio5HA41NDQoKSlJ9fX1mjRpUr/L1tXVyeVyyeVySZKWL1+ubdu29ZjvoYce0uOPP97rOjwPe6L6zssAAGBkyszMVGZmZvj1s88+2+e8hv1NW1NTU/jK0Pb2dvl8Pnm9Xi1ZskTFxcWSpOLiYi1durTbcqFQqNvrKVOm6Pjx42pqapIk+Xw+paSkSJI+/PDD8Hz79+9XWlqaUe0AAACYyrAjbfX19crLy1MwGFQwGNTq1auVnZ0tr9erFStWaNeuXXI6ndq3b194GafTqdbWVnV2dqqkpEQ+n0/Tp09Xfn6+srKyZLfb5XQ6VVRUJEl68cUXVVZWptGjR2vixInavXu3Ue0AAACYyrDQ5na7VV1d3WP6hAkTVFZW1usygUCg1+m5ubnKzc3tMf35558fVI0AAABWwRMRAAAALIDQBgAAYAGENgAAAAsgtAEAAFgAoQ0AAMACCG0AAAAWEJEnIpgtITZBgeaA2WVEpfgx8SN2bBJiEyTJkv0P5HMbiZ9tJHoe6H5j5f3reg11z9GyD0eqjivbiZa+jRI/Jl5SdP9uXNmXo40tdPUjCIYZm83W4ykLAAAA0ai/3MLpUQAAAAsgtAEAAFhAv6GtqKhIs2bNUlxcnOLi4jRnzpzww94BAAAQOX1eiFBcXKzCwkLt2LFDXq9XoVBINTU1evrpp2Wz2Xp9FigAAACM0eeFCBkZGdq7d69cLle36YFAQA899JD+9re/RaTAweJCBAAAYBXXdSFCa2trj8AmSU6nU62trUNXHQAAAP6nPkNbbGxsnwv19x4AAACGXp+nR2+44QbddtttvS508uRJXbp0ydDChgqnRwEAgFX0l1v6vBChrq7OsIIAAABwbXgiAgAAQJS4riNt48aNk81m63OFLS0tQ1MdAAAA/ieOtAEAAESJ6zrS1t7erpdfflknT56U2+3W2rVrFRPT5+wAAAAwUJ9H2lasWKExY8bo7rvv1ltvvSWn06nCwsJI1zdoHGkDAABW0V9u6TO0ud1uHTt2TJJ0+fJlzZ07VzU1NcZVaRBCGwAAsIrreiLCf58K5bQoAACAufo80jZq1CjFxcWFX7e3t+uGG274fCELXT3KkTYAAGAV13UhQldXl2EFAQAA4Nr0eXoUAAAA0WNE/LHakweeVHNHs+Hb8Tf4lZ6Ubvh2IsWq/fRWt1V7McJQjkU0jGtCbIIkDcnveDT0c8VgajG6j8Gufyg/MyNE035gtJHUa28SYhP0/H3Pm13GgI2I0Nbc0SxngtPw7VSerozIdiLFqv30VrdVezHCUI5FNIxroDkgSUNSRzT0c8VgajG6j8Gufyg/MyNE035gtJHUa2+u7ItWwelRAAAACyC0AQAAWAChDQAAwAIIbQAAABZgWGjr6OhQRkaGPB6PUlJStGnTJknS+fPnlZOTo2nTpmnBggVqbm4OT8/KylJ8fLzWrVvXbV27d++W2+1Wenq6Fi1apHPnzkmSduzYodTUVKWnp2v+/Pk6ffq0Ue0AAACYyrDQFhsbq/LyctXW1uro0aMqLy9XZWWlCgoKlJOToxMnTig7O1sFBQXh+bds2aLt27d3W09nZ6c2bNiggwcPyu/3Ky0tTTt37pQkzZo1S1VVVfL7/Vq2bJk2btxoVDsAAACmMvT06JXHYHV2dqqrq0uJiYkqLS1VXl6eJCkvL08lJSXheefNm6exY8d2W0dMTIwSExPV1tamUCiklpYWJScnS5IyMzMVGxsrScrIyNDZs2eNbAcAAMA0hoa2YDAoj8cjh8OhrKwspaamqrGxUQ6HQ5LkcDjU2NjYbRmbzda9QLtdhYWFmjlzppKTk1VXV6c1a9b02NauXbv05S9/2bhmAAAATGTozXXtdrtqa2t14cIFLVy4UOXl5d3et9lsPULa1VpaWrR+/Xr5/X65XC6tW7dOW7du1ebNm8PzvPLKK6qurtbPf/7zXtdRu7dWgdiAJMnpccrpcQ6qLwAAgKFQUVGhioqKAc0bkScijB8/XosXL1ZVVZUcDocaGhqUlJSk+vp6TZo0qd9l6+rq5HK55HK5JEnLly/Xtm3bwu+XlZUpPz9fhw4d0ujRo3tdh+dhz4i+4zMAAIhOmZmZyszMDL9+9tln+5zXsNOjTU1N4StD29vb5fP55PV6tWTJEhUXF0uSiouLtXTp0m7LhUKhbq+nTJmi48ePq6mpSZLk8/mUkpIiSaqpqdHjjz+uN954QzfddJNRrQAAAJjOsCNt9fX1ysvLUzAYVDAY1OrVq5WdnS2v16sVK1Zo165dcjqd2rdvX3gZp9Op1tZWdXZ2qqSkRD6fT9OnT1d+fr6ysrJkt9vldDpVVFQkSdq4caMuXryoZcuWSZJuvfXW8IUNAAAAw4lhoc3tdqu6urrH9AkTJqisrKzXZQKBQK/Tc3NzlZub22O6z+cbVI0AAABWwRMRAAAALIDQBgAAYAGENgAAAAsgtAEAAFgAoQ0AAMACInJzXbMlxCYo0BwwfDvxY+Ijsp1IsWo/vdVt1V6MMJRjEQ3jmhCbIElDUkc09HPFYGoxuo/Brn8oPzMjRNN+YLSR1GtvruyLVmELXX0322HGZrP1uGEvAABANOovt3B6FAAAwAIIbQAAABZAaAMAALAAQhsAAIAFENoAAAAsgNAGAABgAYQ2AAAACyC0AQAAWAChDQAAwAIIbQAAABZAaAMAALAAQhsAAIAFENoAAAAsgNAGAABgAYQ2AAAACyC0AQAAWAChDQAAwAIIbQAAABZAaAMAALCAGLMLiISvl3zd7BIswd/gV3pSuuW2F+m6I2kgvQ3X/q+3L6uOx1DWnRCbIElq7mg2ZP1G8zf4Janfenvr0QiR2o5VmLEfWWnf7UtCbIKev+/5Qa9nRIQ2Z4LT7BIsofJ0ZUTHaqi2F+m6I2kgvQ3X/q+3L6uOx1DWHWgOSOr+3Welcak8XSmp/+/u3no0QqS2YxVm7EdW2nf7cmU/GixOjwIAAFgAoQ0AAMACCG0AAAAWQGgDAACwAEIbAACABRgW2jo6OpSRkSGPx6OUlBRt2rRJknT+/Hnl5ORo2rRpWrBggZqbm8PTs7KyFB8fr3Xr1nVb1+7du+V2u5Wenq5Fixbp3LlzkqRDhw5p1qxZGj16tF577TWjWgEAADCdYaEtNjZW5eXlqq2t1dGjR1VeXq7KykoVFBQoJydHJ06cUHZ2tgoKCsLzb9myRdu3b++2ns7OTm3YsEEHDx6U3+9XWlqadu7cKUm69dZbVVxcrEceecSoNgAAAKKCoadH4+LiJH0evLq6upSYmKjS0lLl5eVJkvLy8lRSUhKed968eRo7dmy3dcTExCgxMVFtbW0KhUJqaWlRcnKypM9Dm9vtlt3OWV4AADC8GZp2gsGgPB6PHA6HsrKylJqaqsbGRjkcDkmSw+FQY2Njt2VsNlv3Au12FRYWaubMmUpOTlZdXZ3WrFljZNkAAABRx9AnItjtdtXW1urChQtauHChysvLu71vs9l6hLSrtbS0aP369fL7/XK5XFq3bp22bt2qzZs3D7iOiqKK8H87PU45Pc5raQMAAMAQFRUVqqioGNC8EXmM1fjx47V48WJVVVXJ4XCooaFBSUlJqq+v16RJk/pdtq6uTi6XSy6XS5K0fPlybdu2rcd8/YW/zK9nDqp+AAAAI2RmZiozMzP8+tlnn+1zXsNOjzY1NYWvDG1vb5fP55PX69WSJUtUXFwsSSouLtbSpUu7LRcKhbq9njJlio4fP66mpiZJks/nU0pKSo9lrl4OAABgODHsSFt9fb3y8vIUDAYVDAa1evVqZWdny+v1asWKFdq1a5ecTqf27dsXXsbpdKq1tVWdnZ0qKSmRz+fT9OnTlZ+fr6ysLNntdjmdThUVFUmS3nvvPX31q1/Vp59+qjfffFM//vGPdezYMaNaAgAAMI1hoc3tdqu6urrH9AkTJqisrKzXZQKBQK/Tc3NzlZub22P63LlzdebMmUHVCQAAYAXcKwMAAMACCG0AAAAWQGgDAACwAEIbAACABRDaAAAALCAiN9c1W6A5YHYJlhA/Jj6iYzVU24t03ZE0kN6Ga//X25dVx2Mo606ITZDU/bvPSuMSPyZeUv/f3b31aIRIbccqzNiPrLTv9uXKfjRYttAwvyutzWbjxrsAAMAS+sstnB4FAACwAEIbAACABRDaAAAALIDQBgAAYAGENgAAAAsgtKGbiooKs0uwFMbr2jBe14bxujaM17VhvK5NNIwXoQ3dRMNOaSWM17VhvK4N43VtGK9rw3hdm2gYL0IbAACABRDaAAAALGBEPBEBAADAKvqKZsP+2aPDPJMCAIARgtOjAAAAFkBoAwAAsABCGwAAgAUM69B24MABTZ8+Xbfffru2bdtmdjlRbc2aNXI4HHK73WaXYglnzpxRVlaWUlNTNXPmTP3iF78wu6So1tHRoYyMDHk8HqWkpGjTpk1ml2QJXV1d8nq9uv/++80uJeo5nU6lpaXJ6/XqzjvvNLucqNfc3Kxly5ZpxowZSklJ0bvvvmt2SVHrgw8+kNfrDf+MHz/etO/8YXv1aFdXl+644w6VlZUpOTlZc+fO1Z49ezRjxgyzS4tKb7/9tsaNG6fc3FwdO3bM7HKiXkNDgxoaGuTxeNTW1qbZs2erpKSE/asfly5dUlxcnC5fvqy7775b27dv19133212WVFtx44dqqqqUmtrq0pLS80uJ6q5XC5VVVVpwoQJZpdiCXl5ebr33nu1Zs0aXb58WRcvXtT48ePNLivqBYNBJScn68iRI7rlllsivv1he6TtyJEjuu222+R0OjV69Gg9/PDDev31180uK2rdc889SkxMNLsMy0hKSpLH45EkjRs3TjNmzNDHH39sclXRLS4uTpLU2dmprq4u/nH9H86ePas//vGP+sY3vsFV8APEOA3MhQsX9Pbbb2vNmjWSpJiYGALbAJWVlWnq1KmmBDZpGIe2jz76qNugTp48WR999JGJFWG4CgQCqqmpUUZGhtmlRLVgMCiPxyOHw6GsrCylpKSYXVJUe+qpp/Tcc8/Jbh+2X9NDymazaf78+ZozZ45+/etfm11OVDt16pQmTpyoxx57TLNmzdI3v/lNXbp0yeyyLGHv3r165JFHTNv+sP024Ka6iIS2tjYtW7ZMhYWFGjdunNnlRDW73a7a2lqdPXtWhw4diorn+EWrN998U5MmTZLX6+Xo0QC98847qqmp0VtvvaUXX3xRb7/9ttklRa3Lly+rurpaTzzxhKqrq/WFL3xBBQUFZpcV9To7O/XGG29o+fLlptUwbENbcnKyzpw5E3595swZTZ482cSKMNx89tln+trXvqZHH31US5cuNbscyxg/frwWL16sv//972aXErUOHz6s0tJSuVwurVy5Un/5y1+Um5trdllR7eabb5YkTZw4UQ8++KCOHDlickXRa/LkyZo8ebLmzp0rSVq2bJmqq6tNrir6vfXWW5o9e7YmTpxoWg3DNrTNmTNHH374oQKBgDo7O/Xqq69qyZIlZpeFYSIUCmnt2rVKSUnRk08+aXY5Ua+pqUnNzc2SpPb2dvl8Pnm9XpOril75+fk6c+aMTp06pb179+r//u//9Nvf/tbssqLWpUuX1NraKkm6ePGi/vznP3MlfD+SkpJ0yy236MSJE5I+/zut1NRUk6uKfnv27NHKlStNrWHYPsYqJiZGO3fu1MKFC9XV1aW1a9dyZV8/Vq5cqYMHD+rcuXO65ZZb9JOf/ESPPfaY2WVFrXfeeUevvPJK+BYDkrR161bdd999JlcWnerr65WXl6dgMKhgMKjVq1crOzvb7LIsgz/36F9jY6MefPBBSZ+f+lu1apUWLFhgclXR7YUXXtCqVavU2dmpqVOnavfu3WaXFNUuXryosrIy0/9ectje8gMAAGA4GbanRwEAAIYTQhsAAIAFENoAAAAsgNAGAABgAYQ2AAAACyC0AQAAWAChDQAipKKiQvfff7/ZZQCwKEIbAACABRDaAIwYFy9e1OLFi+XxeOR2u7Vv3z5VVVUpMzNTc+bM0X333aeGhgZJ0j//+U/Nnz9fHo9Hs2fP1qlTpyRJTz/9tNxut9LS0rRv3z5Jnx9By8zM1PLlyzVjxgw9+uij4W0eOHBAM2bM0OzZs/WHP/whPP3gwYPyer3yer2aNWuW2traIjgSAKxo2D7GCgCuduDAASUnJ2v//v2SpJaWFi1atEilpaW68cYb9eqrr2rz5s3atWuXVq1apWeeeUYPPPCAOjs71dXVpddee01+v19Hjx7VJ598orlz5+pLX/qSJKm2tlbvv/++br75Zs2bN0+HDx/WrFmz9K1vfUvl5eWaOnWqHnroofAjqX72s5/ppZde0l133aVLly5p7Nixpo0LAGvgSBuAESMtLU0+n0/f//73VVlZqdOnT+sf//iH5s+fL6/Xq5/+9Kf66KOP1NbWpo8//lgPPPCAJGnMmDG64YYb9M477+iRRx6RzWbTpEmTdO+99+q9996TzWbTnXfeqS9+8Yuy2WzyeDw6deqUjh8/LpfLpalTp0qSHn30UV15cuC8efP01FNP6YUXXtCnn36qUaNGmTYuAKyBI20ARozbb79dNTU12r9/v37wgx8oKytLqampOnz4cLf5Wltb+1zH1Y9rvnLk7L+PlI0aNUqXL1/u8aD3/172e9/7nr7yla9o//79mjdvnv70pz/pjjvuuO7eAAx/HGkDMGLU19crNjZWq1at0oYNG3TkyBE1NTXp3XfflSR99tlnev/99xUfH6/Jkyfr9ddflyT95z//UXt7u+655x69+uqrCgaD+uSTT3To0CHdeeedPYKc9HmYmz59ugKBgP71r39Jkvbs2RN+/+TJk0pNTdXGjRs1d+5cffDBBxEYAQBWxpE2ACPGsWPH9PTTT8tut2vMmDH65S9/qVGjRmn9+vW6cOGCLl++rKeeekopKSn63e9+p29/+9v64Q9/qNGjR+v3v/+9HnzwQf31r39Venq6bDabnnvuOU2aNEl1dXU9jqpJnx99+9WvfqXFixcrLi5O99xzjy5evChJKiwsVHl5uex2u2bOnKlFixZFejgAWIwt1Nv/IgIAACCqcHoUAADAAghtAAAAFkBoAwAAsABCGwAAgAUQ2gAAACyA0AYAAGABhDYAAAAL+H/I32FnzmeQowAAAABJRU5ErkJggg==",
      "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 environment 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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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 unfortunately 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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext cythonmagic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N_vec  = numpy.arange(25, 2000, 25) * N_core"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAGGCAYAAACnlj1bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8zWcbx/HPCUKoEVvs2kXVVpo2aiuqURpb0FZtumiNRFGKVrW0RYvSpmYRI7Q0qFF77xFSodQmIus8f/yepE2NRM755Zwk3/frlZecc36/+75PnqTPlTvXfV0Wq9VqRURERERETOPi6AWIiIiIiKR1CrpFREREREymoFtERERExGQKukVERERETKagW0RERETEZAq6RURERERMpqBbRERERMRkCrpFREREREzmtEH3mTNn6NmzJ23btnX0UkREREREbOK0QXfJkiWZOXOmo5chIiIiImKzFA26u3fvToECBahcuXKC54OCgihfvjxlypRh/PjxKbkkERERERHTpWjQ7evrS1BQUILnYmJi6Nu3L0FBQRw+fJiAgACOHDmSkssSERERETFVigbdnp6euLu7J3hu+/btlC5dmhIlSpApUyZ8fHxYtmwZV69epVevXuzdu1e73yIiIiKSqmV09ALOnz9P0aJF4x8XKVKEP/74g9y5c/P1118/8l6LxWL28kREREREALBarcm+1+EHKW0NnK1Wqz7+/zFy5EiHr8FZ1mbGfPYa05ZxknPv49xj1rVp/cOZvxZp4WfPXuM688/e41zvzN9vKf3hzF8L/ezZZxxn+tmzlcOD7sKFCxMaGhr/ODQ0lCJFijhwRamXl5eXo5fwUCm9NjPms9eYtoyTnHsf5x5n/h5yZs78dUsLP3v2GteZf/aSO0d658xfM/3s2WectPSzZ7HaI3R/DCEhIbRs2ZIDBw4AEB0dTbly5Vi3bh0eHh7UqlWLgIAAKlSokOhYFouFkSNH4uXl5dQ/eCJpkZ+fH35+fo5ehki6o589kZQVHBxMcHAw/v7+Nu14p2jQ3b59ezZs2MCVK1fInz8/o0aNwtfXl9WrVzNw4EBiYmLo0aMHQ4cOTdJ4FovFLtv9IvL4goOD9cuuiAPoZ0/EMWyNO1N8p9ueFHSLiIiISEqwNe50ePUSM6iqidhCv8iJiIiIvaX6oNvPz++BOd0KnCQ59AubiIiI/FtcTret0mR6idJOJLn0vSMiIiIPYmuM4PCSgSIiIiIiaZ2CbhERERERk6X6oNvPz88ueTYiIiIiIv8VHBxsl9r4yukW+Rd974iIiMiDKKdbRERERMTJKehOYSVKlGDSpElUqVKFXLly4ePjw71795g9ezaenp4JrnVxceH06dMAdOvWjd69e9O8eXOyZ8+Op6cnFy9eZMCAAbi7u1OhQgX27t2bYJ5x48ZRsWJFcufOTffu3bl37x4AlSpVYsWKFfHXRkVFkTdvXvbt25cCXwERERGR9CfVB92pLafbYrGwcOFC1qxZw5kzZ9i/fz+zZ89OUn3ohQsXMmbMGP7++29cXV2pU6cONWvW5OrVq7z66qsMHjw4wfU//vgja9eu5dSpUxw/fpzRo0cD0LVrV+bNmxd/3apVqyhcuDBVqlSx75sVERERSeXsldOdJoLu/zbGcXb9+/enYMGCuLu707JlywQ71A9jsVjw9vamatWqZM6cmVdeeYVs2bLRqVMnLBYL7dq1Y8+ePQmu79u3L4ULF8bd3Z0PP/yQgIAAADp27MjKlSu5ffs2AHPnzqVz587mvFkRERGRVMzLy0tBty0sFvt8JEfBggXjP8+aNWt88JuY/Pnzx3+eJUuWBI/d3NzuG6do0aLxnxcrVoywsDAAPDw8qFevHosWLeL69esEBQXRsWPHZL0XEREREUlcqm8Dn1zOVqAiW7ZshIeHxz++ePGizWOeO3cuweceHh7xj7t27cq3335LVFQUdevWpVChQjbPJyIiIiIPlm53up1FXOmZKlWqcOjQIfbt20dERMR9f8Z43BI1VquVadOmcf78ea5evcqYMWPw8fGJf/2VV15h9+7dTJkyhS5dutj8PkRERETk4VJ90J3aDlL+l8ViwWKxUKZMGUaMGEHDhg0pV64cnp6eCQ5Xxl33sMdxz/378w4dOtC4cWNKlSpFmTJlGDZsWPzrWbJkwdvbm5CQELy9vU18hyIiIiKpl5rjoOY4j1KyZEm+/fZbXnzxxYde89FHH3HixAm+//77FFyZc9P3joiIiDyIrTFCus3pTu+uXr3Kd999x9y5cx29FBEREZE0L9Wnl8jjmzFjBsWKFaNZs2Y899xzjl6OiIiISJqn9BKRf9H3joiIiDyIrTGCdrpFREREREymoFtERERExGSpPuhO7SUDRURERMR5qWQgyukW+9P3joiIiDyIcrpFRERERJycgu5Uys/Pj86dOzt6GSIiIiKSBAq6U4Hg4GCKFi2a4Ln/toAXEREREeeloDuVUt6xiIiISOqhoNsBQkND8fb2Jn/+/OTNm5c+ffqQJ08eDh48GH/NpUuXyJYtG+fOnaNZs2aEhYWRPXt2cuTIwYULF7BYLERGRtK1a1dy5MhBpUqV2LVrV/z9R44cwcvLC3d3dypVqkRgYGD8a926daNPnz60aNGCHDlyUKdOHU6fPp2iXwMRERGR9ERBdwqLiYmhRYsWlCxZkrNnzxIWFkaHDh3w8fFh3rx58dcFBATQsGFDihUrRlBQEB4eHty6dYubN29SqFAhrFYry5cvp3379ty4cYNWrVrRt29fAKKiomjZsiVNmzbl8uXLfPHFF3Ts2JHjx4/Hjz9//nz8/Py4du0apUuX5sMPP0zxr4WIiIhIepHqg+7UVqd7+/btXLhwgQkTJuDm5oarqyv16tWjS5cuBAQExF83d+7c+IOSD0sl8fT0pGnTplgsFjp16sS+ffsA2LZtG3fu3GHIkCFkzJiR+vXr06JFiwTje3t7U6NGDTJkyEDHjh3Zu3evie9aREREJHWyV53ujLYvxbGS+0Ww+NvnIKJ15OPlVoeGhlK8eHFcXBL+vlO7dm3c3NwIDg6mYMGCnDp1ilatWj1yrAIFCsR/njVrViIiIoiNjSUsLOy+g5fFixcnLCwMMA5h/vteNzc3bt++/VjvQ0RERCQ98PLywsvLC39/f5vGSfVBd3I9brBsL0WLFuXcuXPExMSQIUOGBK917dqVefPmUaBAAdq2bYurqyvw4Eolj6pe4uHhQWhoKFarNf66s2fPUr58eTu+ExERERFJqlSfXpLa1K5dm0KFCjFkyBDCw8OJiIhgy5YtAHTq1IklS5bwww8/0KVLl/h7ChQowJUrV7h582b8c4+qXlK7dm2yZs3KJ598QlRUFMHBwaxYsQIfH59E7xURERER+1PQncJcXFwIDAzk5MmTFCtWjKJFi7JgwQLA2AWvVq0aLi4uPPfcc/H3lC9fnvbt2/Pkk0+SO3fu+Ool/93tjnvs6upKYGAgq1evJl++fPTt25e5c+dStmzZ+Osedq+IiIiI2J/Fmoq3PS0WywN3bR/2fGrQo0cPChcuzKhRoxy9lHQpNX/viIiIiHlsjRHSbU63MwoJCWHJkiWqJCIiIiKSxii9xEkMHz6cypUr895771G8eHFHL0dERERE7EjpJSL/ou8dEREReRBbYwTtdIuIiIiImExBt4iIiIiIyRR0i4iIiIiYLNUH3X5+fgQHBzt6GSIiIiKSBgUHB+Pn52fzODpIKfIv+t4RERGRB9FBynTKz8+Pzp07O3oZIiIiIpIECrpTgeDgYIoWLZrgObVtFxEREUk9FHSnUimVAtGtWzfmzJmTInOJiIiIpFUKuh0gNDQUb29v8ufPT968eenTpw958uTh4MGD8ddcunSJbNmyce7cOZo1a0ZYWBjZs2cnR44cXLhwAYvFQmRkJF27diVHjhxUqlSJXbt2xd9/5MgRvLy8cHd3p1KlSgQGBsa/1q1bN/r06UOLFi3IkSMHderU4fTp0w9c66N21Lds2ULNmjXJlSsXtWrVYuvWrfGveXl5MXToUGrXrk3OnDlp3bo1165di39927Zt1K1bF3d3d5555hk2bNiQ4N4RI0bw3HPPkSNHDpo0acKVK1cACAkJwcXFhdmzZ1OsWDHy5MnD119/zY4dO3j66adxd3enX79+8WPNnj2bevXq0a9fP3LlykWFChVYv359Uv5nEhEREbEbBd0pLCYmhhYtWlCyZEnOnj1LWFgYHTp0wMfHh3nz5sVfFxAQQMOGDSlWrBhBQUF4eHhw69Ytbt68SaFChbBarSxfvpz27dtz48YNWrVqRd++fQGIioqiZcuWNG3alMuXL/PFF1/QsWNHjh8/Hj/+/Pnz8fPz49q1a5QuXZoPP/zwoWt+UOB99epVXnrpJQYOHMjVq1cZPHgwL730UoLAeu7cucyaNYsLFy6QMWNG+vfvD8D58+dp0aIFI0aM4Nq1a0ycOJE2bdrEB9Zx73/27NlcunSJyMhIJk6cmGD+7du3c/LkSX766ScGDBjA2LFjWb9+PYcOHWLBggVs3LgxwbWlS5fmypUr+Pv74+3tnWCdIiIiImZT0J3Ctm/fzoULF5gwYQJubm64urpSr149unTpQkBAQPx1c+fOjT8o+bBUEk9PT5o2bYrFYqFTp07s27cPMHaR79y5w5AhQ8iYMSP169enRYsWCcb39vamRo0aZMiQgY4dO7J3794HzmG1Wh84/8qVKylXrhwdO3bExcUFHx8fypcvz/LlywEjUO/SpQtPPfUUWbNm5aOPPmLBggXExsYyb948mjdvTtOmTQFo2LAhNWrUYOXKlfH3+vr6Urp0abJkyUK7du3uW9/w4cNxdXWlUaNGZM+enQ4dOpA3b148PDzw9PRkz5498dfmz5+fAQMGkCFDBtq1a0e5cuXi5xIRERFJCek36LZY7PPxmEJDQylevDguLgm/9LVr18bNzY3g4GCOHj3KqVOnaNWq1SPHKlCgQPznWbNmJSIigtjYWMLCwu47eFm8eHHCwsL+/9YtCe51c3Pj9u3b8Y/j0jTc3d0JCAigd+/e8Y/jdtPDwsIoVqzYQ+cAEqyhWLFiREVF8ffff3P27FkWLlwYP6a7uzubN2/m4sWL8dcXLFjwoev773t3c3O77/GdO3fiHxcuXPiR6xQRERExW0ZHL8BhHFSLuWjRopw7d46YmBgyZMiQ4LWuXbsyb948ChQoQNu2bXF1dQUenN7xqFxrDw8PQkNDsVqt8dedPXuW8uXLJ2mN+/fvj//c19eX+vXr06VLlwTXFC5cmCVLliR47uzZszRr1iz+8blz5xJ8nilTJvLly0exYsXo3Lkz06dPT9J6bHX+/Pn71vnyyy+nyNwiIiIikJ53uh2kdu3aFCpUiCFDhhAeHk5ERARbtmwBoFOnTixZsoQffvghQZBboEABrly5ws2bN+Ofe1T1ktq1a5M1a1Y++eQToqKiCA4OZsWKFfj4+CR674M86PrmzZtz/PhxAgICiI6OZv78+Rw9epQWLVrE3zNv3jyOHDlCeHg4I0aMoG3btvGpMIGBgaxdu5aYmBgiIiIIDg5OEBzbWp3l3/dfunSJKVOmEBUVxcKFCzl27BjNmze3aXwRERGRx6GgO4W5uLgQGBjIyZMnKVasGEWLFmXBggWAsQterVo1XFxceO655+LvKV++PO3bt+fJJ58kd+7c8dVL/rvbHffY1dWVwMBAVq9eTb58+ejbty9z586lbNmy8dc97N4HedBruXPnZsWKFUyaNIm8efMyceJEVqxYQe7cuePv6dy5M926daNQoUJERkYyZcoUAIoUKcKyZcsYO3Ys+fPnp1ixYkyaNClBoPzvOf+73qTUKP/3NbVr1+bEiRPky5eP4cOHs2jRItzd3RMdQ0RERMRe1AbeyfTo0YPChQszatQoRy/FJvXr16dz5850797doeuYPXs23377LZs2bUrS9an5e0dERETMY2uMkH5zup1QSEgIS5YseWglkdRGwauIiIiIQeklTmL48OFUrlyZ9957j+LFizt6OXbhDK3qH5RKIyIiIpLSnDa95M6dO/Tu3ZvMmTPj5eVFhw4d7rsmLaaXiGPpe0dEJH2IjIlk+q7pvFXjLTK4ZEj8Bkn3bI0RnHane8mSJbRr147p06fHN1wRERERsdWpq6d4cc6LrDuzjnsx9xy9HEknUjTo7t69OwUKFKBy5coJng8KCqJ8+fKUKVOG8ePHA0Zt5bjmKv+tZy0iIiLyuO5G3eXdte9Se2ZtXi73MovbLSZrpqyOXpakEykadPv6+hIUFJTguZiYGPr27UtQUBCHDx8mICCAI0eOUKRIEUJDQwGIjY1NyWWKiIhIGmK1Wll8eDGVvqpE6M1QDvc5zLv13sXF4rR/8Jc0KEWrl3h6ehISEpLgue3bt1O6dGlKlCgBgI+PD8uWLaN///707duXlStXJtoOXUREROS/Lty6wO4Lu5m7fy4HLh3gmxbf0PDJho5elqRTDi8Z+O80EjAap/zxxx9kzZqV7777LtH7/fz84j/38vLCy8sLcI7KGSIiIuI47Ra1w4KFukXrMuvlWbhlcnP0kiQVCQ4OJjg42G7jOTzotjU4/nfQHUfVJ0RERNK3S3cuceCvA/z1zl9kzpjZ0cuRVOjfm7kA/v7+No3n8GSmwoULx+duA4SGhlKkSBEHrkhERERSu8BjgTQp3UQBtzgNhwfdNWrU4MSJE4SEhBAZGcn8+fMfK4fbz8/Prlv/IiIikvotO7aM1uVaO3oZkgYEBwc/MLPicaVoc5z27duzYcMGrly5Qv78+Rk1ahS+vr6sXr2agQMHEhMTQ48ePRg6dGiSxlMjExERkfQtKiaKzaGbOXn1JJ2f7syfN/9k7KaxLDu2jJP9T5IrSy5HL1HSCFvjTqftSJkUCrpFRETSp+jYaGbsmsHHv39Mvmz5yOOWh5NXT3Lz3k361+7PWzXeIl+2fI5epqQhtsadDj9IKSIiIvI4tp/fTq8VvciVJRcL2y6kdpHaAASHBFPKvRRFcxZNZASRlJfqg24/P7/7TpeKiIhI2nPpziXe/eVdfjn1C+MajqPz050TVEHzKuHluMVJmmWv0oFKLxERERGnFh4VzoxdMxj7+1i6PN2FES+MIHvm7I5elqQzSi8RERGRNGvT2U20X9yemoVrsrbTWqoUrOLoJYkkS5KC7vPnzxMSEkJMTAxWqxWLxcLzzz9v9tqSROklIiIiaU9EdASfbP6EqTum8n3r72lSuomjlyTpVIqll7z//vvMnz+fp556igwZMsQ/HxgYaPPktlJ6iYiISNoSFRPFDwd+4KONH1GlQBUmN51MsZzFHL0sEfNLBpYtW5YDBw6QObPzdXRS0C0iIpJ27L24lw6LO1AoeyGGeQ6jfsn6No23dClUrgylStlpgZKu2Rp3JtqRslSpUkRGRiZ7AhEREZFHuRd9j9EbR9NobiM+9PyQdV3W2RxwX74MnTvD1Kn3v2a1ws8/Q+PGEB5u0zQiSZZoTrebmxvPPPMMDRo0iN/ttlgsTJkyxfTFiYiISNp1J/IO3+/7nk+3fUrFfBXZ+fpOiucqbpexx42DGjVg+XKYNAksFoiMhBMnYOhQCAmB0aPBzc0u04kkKtGgu1WrVrRq1Sq+DmbcQUpnoYOUIiIiqc/qE6t5c8WbVCtUjZktZ/JCiRfsMm5MDEycCPPmwd69ULs27NsHCxfCN9+Auzu0aQOLFoGrq12mlDQuRet037t3j+PHjwNQvnx5MmXKZPPE9qCcbhERkdQl9EYo7/36HltDtzKz1UwaPtnQ5jGjoyEwEFq3hhEjYN06CAiA4sWhTx+YPx/q1YPPPoMnn7TDm5B0yfQ63cHBwXTt2pXixY0/95w7d445c+bwwgv2+Y1URERE0oeTV09SZ2YdetfszcyWM8nmms0u427ZAt7e0K0brFhh7Gx7eBiv9ehhpJCMGwcZ1Z1EHCjRne5q1aoREBBAuXLlADh+/Dg+Pj7s3r07RRb4KNrpFhERST26Lu1KKfdSjHhhhF3HHTUKTp+GHTtg8GAj0BaxN9N3uqOjo+MDbjBKCEZHRyd7QhEREUl/jv19jFUnVnGy30m7j71hgxFsz5plHJgUcUaJBt3Vq1enZ8+edOrUCavVyg8//ECNGjVSYm1JooOUIiIizuv4lePsDNvJnH1zGFh7IDmz5LTr+Pfuwfbt8NxzCrjFHCl2kDIiIoKpU6eyefNmADw9Pendu7dTNMtReomIiIjzOnPtDHW+rUP9EvUpnL0wfl5+ZM+c3a5z/P47DBgAu3bZdViR+5jekdKZKegWERFxTnej7lLvu3p0rdKVAXUG2H18qxU2bYIxY6BSJaMWt4iZTMvpbtu2LQsXLqRSpUr31eW2WCzs378/2ZOKiIhI2hR0MojZe2ez7c9t1C1al/61+9t1/JgYmDEDJkyALFmMiiVvvGHXKURM8dCd7rCwMDw8PDh79ux9Ub3FYokvIehI2ukWERFxvJjYGHaG7WTi1onsubCHoc8NpU6ROlTIVwEXi4vN4//9N+TJA3/+Ce3aQYYMxs52rVrK45aUY2vc+dCfBI//F7icNm0aJUqUSPAxbdq0ZE8oIiIiaceiw4vIPzE/PQN7UqVAFQ68dYAe1XpQMX9FuwTcR45AkSLw9NNQpw688gps3Gh0mlTALalJojndVatWZc+ePQmeq1y5MgcOHDB1YUmhnW4RERHHWXd6He0Xt2d1x9VU96hu9/FjY+H556F9e6hQwWjb/txzdp9GJElMy+n+6quvmDZtGqdOnaJy5crxz9+6dYt69eole0J7U8lAERGRlBd0MoguP3dhYduFpgTc0dEwdKjx+VtvgYvtm+YiyWJ6ycAbN25w7do1hgwZwvjx4+Mj++zZs5MnTx6bJ7YH7XSLiIikrKt3r/L22rf57cxvzHp5FvVL1rfr+JGRsHo1jB8P2bLB3LlQsKBdpxBJlhQrGXjp0iUiIiLiHxcrVizZk9qLgm4REZGUERMbw4zdMxgZPJLXKr7G2AZjecL1iUTvW7QInnkGSpdOfI5bt8DLC9zcoGdP6NJFO9ziPEwPupcvX87bb79NWFgY+fPn5+zZs1SoUIFDhw4le1J7UdAtIiJijqN/H6V4zuK4ZXIDoOfynhy+fJivXvqKKgWrJGmMK1egeHGoVs1o1f6og4/R0fDyy+DhAdOn65CkOB/TqpfEGTZsGFu3bqVs2bKcOXOGdevWUbt27WRPKCIiIs5rzck1eM32os7MOjT4vgFXwq/w3Z7v2BK6hbWd1yY54Ab4/HN47TUIDzfSRP7NaoUzZ+DOHaP2dvfuxsHJadMUcEva9NCDlHEyZcpE3rx5iY2NJSYmhvr16zNggP07S4mIiIjjXLpziT6r+rD/r/2MfGEkrz71KsPWD+PZb5/lWsQ1NnbbmKR0kjg3bhgB9LZtcO0atGoFNWoYVUgmToRx44zmNvfuQcWKRu3tFSsgUyYT36SIAyUadLu7u3Pr1i08PT3p2LEj+fPn54knkv5DJyIiIs7LarUScDCAwWsG07VKV75v/X18SsknjT6hlHspCucoTIV8FR5jTBgxApo0+SeXe8wYeOEFI+i+exe2b4dSpeD0aVi50tjpzprVjHco4hwSzem+c+cOWbJkITY2lh9++IGbN2/SsWNHp6hgopxuERGR5Dt/8zxvrHiDczfO8W2rb6lVuJbNY8bEQN++RlAdFAT58v3z2tmzsGwZvP66cVhSJDUx9SBldHQ0jRo14rfffkv2BGayWCyMHDlSdbpFREQeU6w11sjbLtkA//r+uGZwtWm8W7dg6VKjPXvevLBkCeTIYafFijhQXJ1uf39/c6uXNGjQgMWLF5MrV65kT2IW7XSLiIgkz8zdM5m1dxabfDfZ1K49MBDefRdCQ6FePRg0CJo21WFISXtM60gZJ1u2bFSuXJnGjRuT9f/JVhaLhSlTpiR7UhEREXGcq3evMmz9MII6BdkUcC9dCm++CQEBRrv2jIlGFSLpV6I/Ht7e3nh7eyd4zqJfX0VERFKtYeuH0aZCG54p+Eyy7j9wAL791gi2V62C6vbvAi+S5iS5I6UzUnqJiIhI0h26dIh3fnmHM9fOsKXHFnK75U7yvadPw0cfwcaNRpm/bt2MA5HFi5u3XhFnYnp6ScmSJR846enTp5M9qYiIiKSciOgIPlz3IXP3z2XY88PoVaPXYx+cnDbN6BoZGAjly6s9u8jjSjTo3rFjR/znERERLFq0iCtXrpi6KBEREbGP3Rd20/nnzlTMV5HDfQ6TN2veZI2zZYtRa/upp+y8QJF0IlnpJdWqVWP37t1mrOexKL1ERETEsPncZpYeXUqrcq14rthzXA6/zMjfRrL4yGImN51M+0rtk30mKyIC8uSBv/4C9ceT9Mr09JJdu3bF/5DGxsayc+dOYmJikj2hiIiI2M+fN//kvV/eY9O5TXSo1IFeK3sRHhXOrXu36Px0Z472PfpYudsPsnu3kVKigFsk+RINut9+++34oDtjxoyUKFGCBQsWmL4wERERebjQG6FM2DKBefvn0btmb2a0nEE212yMaziOfX/tI0fmHDzp/qRd5tqyxajBnaYcPgx9+hhJ6vptQlJAokF3cHBwCixDREREksov2I8vtn9Bj6o9ONT7EIWyF4p/zWKxJLsU4MNs3gyvvWbXIR1r715o3BgmTlTALSnmoTndkyZNMi54SP7X4MGDzVtVEimnW0RE0pvN5zbTblE79r65l3zZ8pk+n9UKBQvCjh1QrJjp09nOaoWQENiwAUaPBj8/6NTpn9ejo6F2bejXz6h7KJJEtsadDy34c+vWLW7fvs3OnTv56quvOH/+PH/++Sdff/21UxyijOPn56fdeBERSReiYqLotbIXnzb+1PSA++JFiIw06nNnygRFi5o6nX1YrUYf+jp1YNEi8PGBGTP+eS0gAJo1g1y5oGtXx65VUo3g4GD8/PxsHifR6iWenp6sWrWK7NmzA0Yw3rx5czZt2mTz5LbSTreIiKQnk7ZMYs2pNazptMa07tCbNxtZF+vWGbFpnTrG8059nOvcORgxAm7fhhOx6Rk/AAAgAElEQVQnjA4+OXMavzUULmw8HjfOOBE6dCi0aAE5cjh61ZLKmLbTHefSpUtkypQp/nGmTJm4dOlSsicUERGRR4u1xjL/4HwqTqtIpWmVGLNxDL+f+52Pf/+Yqc2nmhJwR0XBO+9A+/bQqBFcuACLF0NsLLz6qt2ns5+//jIWXKAAvPACrFljBNwArq7Qrp2RThIZCdu2QYcOCrjFIRI9SNmlSxdq1aqFt7c3VquVpUuX0lV/khEREbErq9XKqWun2HxuM59u+5TMGTLzaeNPyeaajYADAbwy/xX61+5PmTxlkjV+dLQRc27ebKSMjB1r1N4GuHwZWrc2drb37Pnn+Zo1jSwNp/bee9CyJYwf/+DXBw+GypXhzTfBpL8OiCRFkprj7Nq1i02bNmGxWHj++eepWrVqSqwtUUovERGR1C48KpzXA18n6GQQ2TJlo2bhmnR5ugutyrVKsKMd9/93ydnljo01qo8cPgwNG8K1a3D+PAQFGcF4/frg6WnEramqvfv161CihJFSks/8Q6WSvtkadyarI6WzUNAtIiKpWWRMJC//9DJ53PIwodGEBKX/7OnDD41iHuvWQebMEBNjbA6XLGlkZ2TODPPmpcKN4K++gt9+c/KEc0krFHSn3uWLiEg6Fh0bjc8iH2KsMSxsu5CMLolmfCbL99+Dv7+RWvLvzeAbN4xDkvnywS+/GIF3qhIbC9WrG9vzjRs7ejWSDpjeBl5ERETsKzo2mtcDX+fmvZsEtg80LeCeNw/efdfYDP5v9kXOnEanycyZU2HADfDRR5AtGzRo4OiViCSJgm4REZEU9OvpXxkYNJAiOYrw82s/kzmjORHvpEnw+edGwP3UUw++xt3dlKnNt3y5UX97507IkMHRqxFJkkSD7sWLFzNkyBD++uuvBIc4bt68afriRERE0orT104zeM1gDlw6wMRGE2ldvrVptbZHjjSqjmzenEqa2jyOw4ehRw8IDDRaZYqkEomeUX7vvfdYvnw5N2/e5NatW9y6dUsBt4iISBJFx0YzcctEas2oRe3CtTnU+xCvVHjFtID7m2/gxx+NHe40FXBfuQJTp0KTJkb3nriuPSKpRKI73QULFqRChQopsRYREZE0Zf9f++mxvAc5Mufgj55/UCp3KVPnW7EC/Pxg0ybIn9/UqVLGqVNGGsnatXDyJDRtavxG4enp6JWJPLZEq5cMGDCAixcv0rp1a1xdXY2bLBa8vb1TZIGPouolIiLijK7dvYb/Bn9+PPAj4xqOw/cZX9N2tuP89ptRizsw0GjAmOrt2gUvvQRduxpt22vXNjpMijiI6dVLbty4gZubG2vXrk3wvDME3SIiIs4kJjaGGbtnMDJ4JK+Uf4VDvQ+RL5v5TVt++gkGDDDKVaeJgHvtWujUCaZPN1pliqQBqtMtIiJio6iYKH46+BPjNo8jf7b8TG4ymSoFq5g+759/wuzZRh73qlVGt/NUKybG6ODzyy8wa5ZxEvS55xy9KpF4tsadiR6kDA0N5ZVXXiFfvnzky5ePNm3a8OeffyZ7wqQ6c+YMPXv2pG3btqbPJSIikhzhUeF8uf1LynxRhll7ZzG5yWTWd1lvasD955/Qu7fR/bxqVdizx6hS4vQB940bD3/NajXeVL9+EBVlFBBXwC1pTKJBt6+vL61atSIsLIywsDBatmyJr6+v6QsrWbIkM2fONH0eERGR5Nh8bjNPf/U0v5z+hfmvzmd91/U0KtXItNzt69dhyBCoUgVy5IA1a+DSJVi8GIoVM2VK+7lzx/gt4aefHvz6Z5/B1q1G28yJE+HJJ1N0eSIpIdGg+/Lly/j6+pIpUyYyZcpEt27duHTpUpIn6N69OwUKFKDyf34FDwoKonz58pQpU4bx48c//spFREQc4F70PYb+OpRXF77KhEYTWOazjNpFzE2kPnDA2NW+fBn274dx46BcOTD5bGbyXLgAI0ZAo0ZGD3owgu0iRWDgQONNxPnzT+P059SpxgnQ7Nkds2aRFJBo0J0nTx7mzp1LTEwM0dHRzJs3j7x58yZ5Al9fX4KCghI8FxMTQ9++fQkKCuLw4cMEBARw5MgR5s6dy6BBgwgLC3v8dyIiImKyPRf2UGtmLY78fYR9vfbxSoVXTJ9z5Uqj0/mYMfDtt1C4sOlT2ubtt+H4cejc2fj8r7+MA5HjxkHHjsYByS++gO7d4emnoWxZOHgQihd39MpFTJXoQcqQkBD69evHtm3bAKhbty5ffPEFxR7jb1khISG0bNmSAwcOALB161b8/f3jg/Fx48YBMGTIkPh7rl69ygcffMC6devo2bMn77///v2L10FKERFJAdfuXmP4b8NZdHgR4xuOp0uVLqalkWzebKQ0X7kCYWHw669GCsmzz5oynX2dOmWUTzl92siBeecd2L4dQkLgzBm4dw8+/hiuXTN2vl9/HfLkcfSqRZLE9JKBJUqUIDAwMNkTPMj58+cp+q82WUWKFOGPP/5IcE3u3Ln5+uuvEx3Lz88v/nMvLy+8vLzstUwREUnnYq2xzNk7hw/Wf0Drcq053Ocwud1ymzbfV1/BqFHQoYMRi5YsCWPHGvFpqjBxIvTqZQTcAMOHGzvZffpAhgyQNSt89JFj1yiSRMHBwQQHB9ttvIcG3ePHj+f999+nX79+971msViYMmVKsie15+7Av4NuERERe9lzYQ99VvUhOjaawPaB1PCoYdpcsbEwdCgsXQq//w6lzG1caV/nz8PcueDuDvPnw9Gj/7yWM6fRHjPV/NYg8o//bub6+/vbNN5Dg+6nnnoKgOrVqycIkq1Wq81Bc+HChQkNDY1/HBoaShH9QIqIiBO4HnGd4euHs+DwAkbXH02Paj1wsSR6BOqxRUfDoUNG9sWSJXD7tpFWkqqyLSIjwdvb2JLPksXI2/5v//myZR2zNhEn89Cgu2XLlgBkzZqVdu3aJXhtwYIFNk1ao0YNTpw4QUhICB4eHsyfP5+AgIBkjeXn56e0EhERsVmsNZa5++YyZN0QXi73Mod7HyZPVnMi4DVrjKIdBQsaKdAvvwzduhlxa6ry/vvGmwgIcNJSKiK2s1eaSaIHKatWrcqePXsSfe5h2rdvz4YNG7hy5Qr58+dn1KhR+Pr6snr1agYOHEhMTAw9evRg6NChj794HaQUERE72BW2i36r+xEdG83U5lOpWbimaXNt2WIE2T//nEr7v/z6K/j5GYXD79yBXbsgt3l57iLOwta486FB9+rVq1m1ahXz58/Hx8cnfpJbt25x+PBhtm/fnuxJ7UVBt4iI2OLv8L/5cN2HLDu2jDEvjsG3qq8pqSRx9u83ylfPmQNNm5o2jXmOHQNPT6Pk31NPGWX+4g5NiqRxplUv8fDwoHr16ixbtozq1avHT5IjRw4+++yzZE8oIiLiaDGxMXyz6xv8gv1oX6k9R/seJVeWXKbOefIkNGsGn3+eygLuKVOMOoaNG8OECUbB8Ndec/SqRFKdRNNLIiMjcXV1Tan1PBaLxcLIkSOV0y0iIkm2+dxm+q7uS87MOfmy+ZdUyl/J9DlXrjR6wYwZAz17mj6d/axdC76+MHKkkVZSqZLRbVIkHYnL6fb39zcnvSTO8ePH+eCDDzh8+DB37941brJYOH36dLIntRell4iISFJdvH2R9399n3Wn1zGx8UReq/iaaQ1u4kRHw7Bh8MMPxlnDVJXDHRYG1avDjz9C/fqOXo2Iw9kadyaauObr60uvXr3ImDEjwcHBdO3alY4dOyZ7QhERkZRktVqZvms6lb+qTMFsBTna9yg+lXxMD7j/+MOIVffuhd27U0HA/d9g4o03jEY3CrhF7CLRoPvu3bs0bNgQq9VK8eLF8fPzY+XKlSmxtiTx8/Oza7cgERFJO0Kuh9B4XmNm7p5JcNdgxjcazxOuT5g2X2wsrFgBL7wAPj5GZ8lVqyBfPtOmtI85c+CZZ4w27QA7dsC+fUbHHpF0Ljg42C7NGBNNL6lbty6bNm3i1VdfpUGDBnh4eDB06FCOHTtm8+S2UnqJiIg8SKw1lm92fsPw34bzbt13ebvu22R0eWjtgGSJijK6mwcGGpvEVqvRK6ZkSXjvPXj1Vcho3ynNcfSoUZGkTBlo0QI++ABat4YGDeABXalF0ivTqpfEmTx5MuHh4UyZMoXhw4dz8+ZN5syZk+wJRUREzLT/r/0MDBrInag7bPTdyFP5nrL7HOHhENc3bts2cHU1esO4uBi72qmmT0xEhFGJZMwYo5ZhzZrw9NNGbkwym9aJyIMlutPtzLTTLSIicY79fQy/DX78duY3hj43lD61+ti8uz1sGJw6ZXSLbNgQMmSAq1eNDeEyZWDmTMiUyT7rT3F79xpb8rlywfz5xm8Ko0eDv78RhL/3nqNXKOJUTD9I2ahRI65fvx7/+Nq1azRp0iTZE4qIiNjT2etn6b6sO/W+q0elfJU42f8kA+oMsDng3rsXZsyAOnXgww+hRAkj88LTE+rVg1mzUmHAHRsL69aBt7dRNPyll+D77//Zmn/nHXj9dXjrLceuUyQNSjTovnz5Mrly/dMwwN3dnb/++svURT0OHaQUEUmfwqPCGfLrEKpNr4ZHdg9O9DvBh89/aJeDklYrDBxobPoOGAA7dxoHIiMjjTzuCROMVJJUZd06Y3t+8GB48UVjC3/AAMiS5Z9rsmSBadMge3bHrVPEyaTYQcrq1auzZMkSihcvDkBISAje3t7s3r3b5sltpfQSEZH0ae2ptby18i1qFa7FZ00+o+ATBe06/sKFRqbF7t1GSkmq99dfRnWSmTOhefNUlHQu4jxMP0g5ZswYPD09ef755wHYuHEj06dPT/aEIiIiyXXpziUGrxnM5tDNTGs+jWZlmtl9jrt34d13YfbsNBJwW61GO8zu3Y10EhFxiCQdpLx8+TLbtm3DYrFQp04d8ubNmxJrS5R2ukVE0odYayyz985m6LqhdHm6C35efmRzzWbKXB99BPv3G7vdacK0aUYC+pYtqTAJXcR52Bp3PjToPnLkCBUqVGDXrl0JJonr4FWtWrVkT2ovCrpFRNK+jWc3MnjNYDK6ZOSrl76iaqGqps0VGgpVqxo53CVKmDZNyjh0CMaOhV9/hU2boGxZR69IJFUzLb3k008/ZcaMGbz99tsPbJX722+/JXtSe/Lz88PLywsvLy9HL0VEROzo1NVTvPfre+wM28m4BuNMb90eGwtvvw29e6eBgHvkSPj6axg0CL76CnLkcPSKRFKt4OBguxTtUJ1uERFxKtcjrjN642hm7Z3F28++zaA6g3DL5GbqnOfPQ5cuRpfJ1ashmzmZKynjxx+NAuPbt4OTpIOKpAWm7XQvXrz4kTsK3t7eyZ5URETkv2KtsczaM4sP139Iy7ItOdT7kN2rkjzIzz9Dr15Gx/OhQ5348OSoUVCwILzxxsOv2b3bKAO4bp0CbhEn89CgOzAwUEG3iIikiJ1hO+mzqg8uFhdWdVxFtULmnxv6+2+j6c2vv8KyZUYTHKdltRrl/m7ehNq1oUqVf147d844+XniBHz+uZFO8vTTjluriDyQ0ktERMRhroRf4YN1H7Ds2DI+bvAxXZ/piovF/l1nrFbYvBm2bYMdO4yPv/+Gtm3hs89SQcrzkSPQpInRnv3jj42Tnm5uRoA9YgTUrGk0vnnhBWjTxtGrFUmTTKteEufvv//G39+f33//HYvFgqenJyNGjCBPnjzJntReFHSLiKROMbExzNw9kxHBI3it4muMqj+KXFlyJX5jMkRHG13Ng4OhaVMjPq1ZE8qVS0VdJSdPNqqRzJgBnToZnSNdXIwygMuWQalSjl6hSJpnenMcHx8fXnjhBZYsWYLVauXHH3/ktdde49dff032pCIikn6tObmGd355h9xuuVnbaS1VClZJ/KaHiImBe/cga9YHv373LnToALdvG+nOqba7+Zo10KOH8fnUqUZdw2eega1bU/GbEklfEt3prlSpEgcPHkzwXOXKlTlw4ICpC0sKi8XCyJEjVTJQRCQVOHjpIO+sfYdT104xodEEXi73ss0lAEeOhE8/NQ5CDhoEHh7/vHb9Orz8MhQqBN9/D66uNr4BR4mIgHz5jNxtd3fjubt3jd1utXMXMV1cyUB/f3+bdroT/cNa48aNCQgIIDY2ltjYWObPn0/jxo2TPaG9xdXpFhER53Tx9kXeCHyDF+e8SLPSzTjU+xCty7e2OeD+80/48ksICjJK/VWqBD17wrFjcOGCkd5cpYpRQS/VBtxgNLapXPmfgBuMfG4F3CIpwsvLCz8/P5vHSXSn+4knniA8PByX/ye+xcbGku3/BUwtFgs3b960eRHJpZxuERHnFRkTyaQtk5i0dRJdq3Rl2PPDcHdzT/zGJOrWzdjZHjvWeHzlihGET51qNLoZPNgoAZjqY9N33jFSSEaOdPRKRNI10w9SOjMF3SIizmnzuc28ueJNiucqzpSmUyiV274H/fbsgWbN4Pjx+yuP3LljPF/VvG7x5rlyxahhWLv2P20xK1c2DlA6dU1DkbTPtKD7yJEjVKhQgd27dz/wxmrVzK+hmhgF3SIizuV6xHWG/DqE5ceWM7npZNo+1dburdutVmjYEF591ahKkibMm2fU4d6zxwiu9+wxygHWqWPU3L50yYm79oikD6ZVL/n000+ZMWMGgwcPfuB/MH/77bdkTyoiImmL1Wpl0eFFDFwzkJZlW3K4z2HTSgCuWgVhYfD666YMn/KOHDFyYb791vhtws3NKCTu4wO5c0ODBgq4RdKARNNLIiIiyJIlS6LPOYJ2ukVEHO/s9bP0WdWHM9fPML3FdOoVq2faXNHRxsbvJ59AixamTZOy2rc33tTQoQmfv3HDaOnepg20bOmYtYlIPFvjzkSrl9StWzdJz4mISPoSER3BmI1jqDa9Gs8WeZY9b+4xNeAGIwOjUCF46SVTp0k5Bw/C+vXQt+/9r+XMCbNnK+AWSSMeml5y4cIFwsLCCA8PZ/fu3Vit1vhqJeHh4Sm5xkeKKxmosoEiIinDarWy4vgKBq0ZROUCldnx+g6edH/S9Hlv3gR/fyO9JNVXJInj7/9PdRIRcUpxdbpt9dD0ktmzZzN79mx27dpFjRo14p/Pnj073bp1w9vb2+bJbaX0EhGRlHX8ynEGBA3gzLUzTGk2hcalUq5vwwcfGLncs2en2JSJW7rUqKP9/PPGh/tjlETctw+aNIFTp+D/pXhFxHmZXjJw8eLFtGnTJtkTmElBt4hIyrh17xajN47m2z3fMvS5ofSr3Q/XDOZ3nImONuLaL76AM2dgyxYoUsT0aZNm924jaO7d22jHvnUrlC5ttL+sXDnx+195xQjUBw0yf60iYjPTgu5JkybFD/7v6iVxjwcPHpzsSe1FQbeIiLmsVis/HfyJd395l4ZPNuTjBh9TKHuhJN8/bRoULGjEl4+TEnL5slGa+quvoHhx6NcPvL0hU6ZkvAkzXL8O1asbnXlee814LjISJkyAbdsgMPDR92/aZFQnOXnSqFYiIk7PtJKBt27demCpwP8G4SIikjadvX6Wt1a+xZ83/2Rh24U8W/TZx7p/82YYPRry54fx4+Hjj+HFFx99z44dRqC+dKkRqC9f7oRNbqxW8PU1uvPEBdxg9Jp/+21jt3vnTvhXamYC4eHQvbvxRhVwi6Qb6kgpIiIJxMTGMG3HNPw3+DOoziDeq/cemTI83hZzZKQRLPv5GRXv5s+H4cOhVCljc7h69X+uDQ83Xp82Df7+G3r1gh49IG9e+74vu5k0yVjwpk2QOfP9r3/5JaxZ8/Dd7kGDjGY3P/xg7jpFxK5Mz+n29fV94KTfffddsie1FwXdIiL2dfjyYXou70lGl4xMbzmd8nnLJ2uc0aP/ybKI++NoZKTR/+Wjj8DT0wiuAwONFOg6dYzU6CZNnLgPzK1bsGgRvP8+bN/+T5v2/4qIMHa7ly69f7f799+hXTs4cADy5DF9ySJiP6all8R56aWX4tNJ7t69y88//4yHh0eyJxQREecTGRPJx5s+5ssdX/JR/Y94o/obuFgSbeXwQMePw+TJsGtXwjxuV1ejbXuXLjBlCvTpA61aGSklJUva6Y3Y27VrRvC8ZAls2AD16sHixQ8PuAGyZDEC81GjjPyYOHFpJVOnKuAWSYceO70kNjaWevXqsXXrVrPWlGTa6RYRsd3W0K28Hvg6T7o/ybSXplEkR/LLg1itRtfyli1TeVGOiAjjN4MJE4wKI23aGB15cuZM+v2lShlB9zPPGMXFJ0yAwoUhIMDctYuIKUzf6f6v48ePc/ny5WRPKCIizuHi7YsMWz+MVSdWMbnpZNo+1dbmg/Jz5hhNbPr1s9MiU1psLPz0k1EUvGpVIx2kXLnHHydLFhgyxDhweeuWcZq0T5+EBy9FJF1JNOh+4okn4v8jbLFYKFCgAOPHjzd9YSIiYo570feYvG0yE7ZMwPcZX470OULOLEncwX2Ey5eNrIrVqyHjY2/pOIHYWGOb/s4dI9H8+edtG+/11+Gvv+Dll6FmTfusUURSLVUvERFJJ6xWK0uPLuWdX96hUv5KTGw0kTJ5ytht/M6doUABmDjRbkOmrPXrYeBA2LsXXJKXzy4iaVeKpJfs37+fkJAQoqOj459zhjbwAH5+fnh5eeHl5eXopYiIOK39f+1nYNBALt25xDctvqHhkw3tOv4vvxgV9A4dsuuwKWv6dHjjDQXcIpJAcHAwwcHBNo+TpJKBBw4coGLFirj86z9Es2bNsnlyW2mnW0Tk0S7ducTw9cNZemwpI18YyRvV3yCji31zP8LDja7nX3wBzZvbdWj72bbNqK/9008Prkl4+TKUKQMhIZArV4ovT0Scn+k73X/88QeHDh1SF0oRkVQkKiaKL7d/ydjfx9KpcieO9jmKu5u7KXP5+xspy04bcJ88abS3zJ4dvvvOyLX+rzlzoHVrBdwiYppEg+6aNWty+PBhKlasmBLrERERG609tZYBQQMonrM4m3w3JbvBTWKuXDFSoLdsMVq+O6XLl4127X5+xm8GL70EPj5GAB7HajVSS2bPdtQqRSQdSDTo9vX15dlnn6VgwYJk/n+7W4vFwv79+01fnIiIJN3pa6cZvGYwBy4dYHKTybQo28K0v1IuWQJ9+xrNFffvh2zZTJnGNnfvGt132raFN980nmvSBMaNgzFj/rluwwajc8+zzzpmnSKSLiSa012qVCk+++wzKlWqlCCnu8SjunGlEOV0i4jAncg7fPz7x3y18yveefYdBj07iCwZs5gy16VLRrC9b5+RqVGvninT2C4mxgi23dxg7tx/DkeePw9PPw27d0Px4sZzHToYfej793fcekXE6dkadyYadD/77LNO0X3yQRR0i0h6ZrVamX9oPu/+8i7PF3+e8Q3H29RN8tFzGY0UBw2Cbt2MbA03N1Omejxjx8KMGfDkk8ZByLJljX+DguDwYePf//+VNp6fH5w4AT/8AH//DaVLw+nTkDu3Q96CiKQOpgfdvXv35vr167Rs2RJXV9f4SZ2hZKCCbhFJr/Zd3Ee/1f24FXmLKU2n4Fnc07S5wsLgrbfg1CmYNcuJ+ryMG2ccgAwIMLbgT5yA48eNf8EIqt0fcHj0zh2jy+TixUYy+t69RjMcEZFHML16SXh4OJkzZ2bt2rUJnneGoFtEJL25dvcaw9YPY9GRRYzyGkXPaj3J4PKAEnh2snGj0bn89ddhwYL7N40d5rPP4NtvjXxsDw/jucaNk3ZvtmwwejQMHgxXrxo75SIiJlNHShGRVMBqtTJv/zze+/U9WpdrzZgGY8jtZl46hNUK06bBqFEwbx40amTaVI9v2jSYMMEIuIsVS94YsbFQowZERBgdfVQWV0QSkSIdKUVExHEOXTpE71W9uR15m2U+y6hVuJap8927B336wB9/GOUAS5Uydbr7RUcbDWweFAh/+y2MHw/BwckPuME4WPn990bdQwXcIpICFHSLiDip25G3GbVhFLP2zsLvBT961ehlaioJwIUL0KYNFCoEW7fCE0+YOt2D1apl5GaXLm18lClj/HvjBnz6Kfz2G5Qsafs8lSrZPoaISBK5POrF2NhYFixYkFJrERERjFSSJUeW8NTUp7h4+yIH3zpIn1p9TA+4t237p3/MokUOCrhPnzZOboaGGrvabdsaOdibNsEvv8DatUYQLiKSyiSa0129enV27dqVUut5LMrpFpG05tTVU/Rb3Y+zN84ytflUvEp4pci8330HQ4YYcW7Lliky5YN99plR6k+HG0XEydgadz5ypxugUaNGTJw4kdDQUK5evRr/ISIi9hMRHYF/sD+1ZtbCq4QXe97ckyIB95070K+fkSa9caODA26AZcvg5ZcdvAgREftLdKe7RIkS97URtlgsnD592tSFASxbtoyVK1dy8+ZNevToQaP/HJ/XTreIpAWrT6ymf1B/KuevzOSmkymW04YDgkkUEQHffGOUun7xRZg6FXLlMn3aR7tyxWhyc/Gik3TeERH5h+nNcZzB9evXeeedd5g5c2aC5xV0i0hqFnI9hEFrBnHw0kGmNJ1CszLNTJ8zMtJIJRkzBqpXB39/qFLF9GmT5vvv4eefjQ8RESdjenpJZGQkn3/+OW3atOHVV1/liy++ICoq6rEm6d69OwUKFKBy5coJng8KCqJ8+fKUKVOG8ePHP/T+0aNH07dv38eaU0TEWUVERzB642iqT69O9ULVOfDWAdMD7uhoo5tkuXJGBseSJbB0qRMF3KDUEhFJ0xLd6e7RowfR0dF07doVq9XK3LlzyZgx4327zo+yadMmnnjiCbp06cKBAwcAiImJoVy5cvz6668ULlyYmjVrEhAQwM6dO9m9ezfvvvsuhQoVYsiQITRu3JgGDRrcv3jtdItIKhOXSlIpfyU+a/IZJXKVMHW+mBj46SdjR9vDw2jE+Nxzpk6ZPBERUKCA0Ws+b0nqhFAAACAASURBVF5Hr0ZE5D6mN8fZsWMH+/fvj3/coEEDnn766ceaxNPTk5CQkATPbd++ndKlS1OiRAkAfHx8WLZsGUOGDKFz584ATJkyhXXr1nHz5k1OnjzJm2+++Vjziog4i7hUkgN/HWBKsyk0L9Pc9DlXrYJ334WcOeGrr4zcbaftA7N+vbHtroBbRNKoRIPujBkzcvLkSUqXLg3AqVOnyJjR9p4658+fp2jRovGPixQpwh9//JHgmv79+9O/f/9HjuPn5xf/uZeXF15eXjavTUTEXiKiI5i4ZSKfbfuMQXUGEdAmgCwZs5g659mzMHAgHDxoVOB76SUHBtsnT0KRIpAlkfes1BIRcTLBwcEEBwfbbbxEo+cJEybw4osvUvL/3b9CQkKYNWuWzRP/tyJKcv076BYRcSZxqSQV81Vk1xu7TE8liYyESZNg4kQj6A4ISDzWNdW1a1CnDnTsCJ9//vDrYmNh+XKjZqGIiJP472auv7+/TeMlGnTXq1eP48ePc+zYMQDKlStn04RxChcuTGhoaPzj0NBQihQpYpexRUQcyRGpJOvWQZ8+Rrf0HTuMynsON3IkNG4MCxYYgXetWg++bscOcHdXp0kRSdMSrV5St25dsmTJQpUqVahSpQpZsmShbt26Nk9co0YNTpw4QUhICJGRkcyfP59WrVo99jh+fn523foXEUmue9H3GLNxDNWnV6dawWoc7H3Q9IA7LAzat4cePYwGN4GBThJwHzxonOCcMsXYfn/9dXhY5SulloiIEwsODrZLZsVDq5dcuHCBsLAwOnbsyI8//ojVasVisXDz5k169erF0aNHkzxJ+/bt2bBhA1euXCF//vyMGjUKX19fVq9ezcCBA4mJiaFHjx4MHTr08Rav6iUi4iTWnlpL31V9qZCvApObTKake0lT54uOhi++MOptv/kmfPghZM1q6pRJZ7VCo0ZGIN2vn/G4WTOoXx/ef//+6ytWNPrP16mT8msVEUki05rjzJkzh9mzZ7Nz505q1KgR/3z27Nnp1q0b3t7eyZ7UXhR0i4ijhd4IZdCaQey+sJspzabQomwL0+fcuBH69oWCBeHLL6FsWZMn7NPHCIx7907a9T//DMOHw969EHfw/swZqFkT/vgDSpX659qTJ8HTE86fB5dE//gqIuIwpnekXLx4MW3atEn2BGZS0C0ijhIZE8lnWz9jwpYJ9K3Vl/frvY9bJnNblx85AkOGwP798Mkn8OqrKVCVZPduo/xJ5szG5L16Pfr6iAh46imYMQP+219h4kQICoJffvln4ZMmwbFjMH26OesXEbET0ztS7ty5k+vXr8c/vnbtGsOGDUv2hPamnG4RSWnrz6ynytdV2HhuI3/0/AM/Lz9TA+6LF41Y94UX4Pnn4ehRaNs2hcoAjhgBH3xg1NEeO9ZIA3mUSZPgmWfuD7jBKKly5QrMnfvPc8rnFhEnZ3pOd5xnnnmGvXv3JniuatWq7Nmzx+bJbaWdbhFJSedvnuedX95ha+hWPm/6Oa3KtbJb+dMHuX3biGGnTAFfXyP2zZ3btOnut20btGsHx48btQdPnDDysseMga5d77/+zz+NgHvHDij5kJz2nf9r787jbK7bP46/UMiNtCBbd8q+JMKdZRghWQZJ1ky2lH0pil8YW3cRkjWSrXvELQwxZmwzspN93xpZKm4l0zAzzsz398cnJGnMzPnO95wz7+fj4ZE5c873c83ky3U+c32ua6fZOT9wwLxreOop+Oknh3sbiogkzfad7sTERGJjY29+fO3aNeLj41O8oIiIt7mecJ2xm8dSblo5nnroKQ51P0STEk1sS7hdLlNtUayYyXd37jSVGWmacIPZ5X7vvVsJcdGisGYNDBwIwcF3Pv+dd8yW/N0SboCKFU37wLfeghUroE4dJdwiki4k2ae7bdu21K5dm44dO2JZFrNmzSIwMDAtYrsnQUFBmkQpIraJjIqk+8ru5M+Rn82dNlPsEftOLVoWfP21yV3z5jXzYv5wjj1tbdhgDjl26HD74yVKmJrsOnXg/vtNnQvApk0QGQmffpr0tYcPhzJlYMsWk9SLiHgwd02mTLK8BCA0NJS1a9cCULduXerVq5fqhd1B5SUiYpcff/uRt8PfZsPpDYyvN55mJZvZWkqyYwf07w8XL5pDkg0aODi63bLA398k3O3b//Vz9u6FevVg6lRTk12pEvTrZ3ax78XKleZ1P/wAjz7qrshFRGxje/cSMKPfjx8/Tt26dbl69SoJCQnkyJEjxYu6i5JuEXG3hMQEpuyYwvANw+lUvhPv1XiP7Jmz27bexYumKUhoqNkAbt/+Vpc9x6xZY9oEHjz498Hs2mX6bzdoYOpgNm5M3juFH36AfPlSH6+ISBpIbd6Z5F/t06dPZ8aMGfz888+cPHmSs2fP0rVr15s73yIivmL7ue10XdGVHJlzENk+klK5S9m2VkKCqcQICoJ27UxHkpw5bVvu3lmW6bEdFJR09l+hgqmHadLEjMJM7ta8Em4RSUeSTLonT57M9u3bee73SWHFihXjwoULtgcmIpJWfrn2C4PWDmLp0aWMrjOaV59+1dZSki1bzEZyzpymE1+ZMrYtlXyhoRAdDS1b3tvzK1Uyg20cq4UREfEOSXYvyZIlC1myZLn5scvlsvUfo+RSn24RSSnLspi3dx6lppgd7UPdDtGuXDvb/o67cAE6djRDbd5+G9av97CE+8Yu97BhyZsO6UH/JoiIuFua9enu378/uXLlYu7cuUyaNIkpU6ZQqlQpRo0alerFU0s13SKSUocuHqLbim5Ex0czteFUKheobNtaCQkwbZrJZdu1g6FDPaSU5M+WLDGF5d9+q5HsIiJ/YvtBysTERD777DPCw8MBqFevHp07d/aI3W4l3SKSXDHxMYzcMJLPdn/G0JpD6VqxK5kyZrJtvT+Wkkya5GE723+UmAjlysG//w2NGjkdjYiIx7E96V68eDENGza8rcTEUyjpFpHkWHZ0Gb1Ce1G1UFXGvjCWfDnsO8h34YLptx0ebgbbtGrl4VUYCxbA+PHmXYJHByoi4gzbJ1IuW7aMokWL0q5dO77++mtcLleKFxMRccLpy6dp8mUT+q/uz8zGMwl+Odi2hNvlurWj/cgjcPgwtG7t4Xmsy2VqXoYP9/BARUS8V5JJ9+zZszlx4gTNmzdn/vz5PPnkk3Tq1CktYrsnOkgpIncTnxDPBxs/4Nnpz1IpfyX2vbmP2k/Wtm29yEh49llYvNgckvzoIw+t3f6z4GDInRvq1nU6EhERj5NmBylviI+PJywsjM8//5wNGzZw6dKlVC+eWiovEZG7iYyKpNvKbjyR6wkm1p/Ikw89adtaZ8+aaZKbN5tEu3lzL9owjouD0qVh5kyoWdPpaEREPJbt5SUrV66kffv2FC1alEWLFvH666/z008/pXhBERE7XYi5wGtLX6PdknaMqDWCr1t/bVvCHRsL778PzzwDRYuaUpJXXvGihNuyoGdPKF9eCbeIiM2SHI4zb948WrZsybRp08iaNWtaxCQikiJfHfqK7iu78+rTr3Ko+yHbxrdblhnE2LcvlC0L27fDk/ZtpNtn2jSzPb9li9ORiIj4vHsuLwFTP+2OmhZ3UXmJiABcjr1Mz9CebD27lTlN51C1UFXb1jp2DPr0gVOn4JNP4IUXbFvKXhs2mG35zZvhqaecjkZExOPZXl7yRyEhISleSETEDqtPrubpqU/zYJYH2fPGHtsS7uho0wKwalV4/nnYt8+LE+7vvzdj3r/4Qgm3iEgaSbK8RCTVzp2DHDm8pI2DeIuY+BgGrB7A8mPLmdl4JnWfsqfzhmWZ5h7vvAO1a8P+/ZDPvvbe9rt6FZo2NSc/1a1ERCTNJGune+fOnYDpZOIp1DLQCyxcCAULwosvmhrS8+edjki83OYzm3nm02eIjo9mX9d9tiXcu3eDnx+MG2f+GM+Z4+UJt2VB586mW0nfvk5HIyLiFdKsZWDNmjWZPXs2hQsXBmD79u107tyZffv2pXrx1FJNtxeJjoawMFi6FFauNK0emjQxO24lS3pRuwdxUpwrjqCIIGbtmcWUhlNoVrKZLev88AMMGwZLlsDIkdCxI2Syb1J82hkzxrx72LABHnjA6WhERLxKavPOJMtLBg0aRP369enZsyfnzp0jNDSU2bNnp3hBSady5DDNi5s3h+vXzT/6ISFm9ztLFpN8N2kCVar4SHYj7rb3x70ELg2kcK7C7H1zL3mz53X7GpcuwejRMGMGdOgAR47AQw+5fRlnrFplxrxv26aEW0TEAffUvWT9+vXUrVuX3Llzs3v3bh577LG0iC1J2un2AZYFe/aYHfClS+HHH6FRI5OE16mj5EBwJboYs2kM47aO46O6HxFYLpAMbv7JSHQ0fPwxTJhg3he+956piPIZx49D9erw1VfmvyIikmypzTuTTLpHjBjBggULmDFjBvv27WPcuHGMHTuWRo0apXhRd1HS7YO++87sgIeEwK5d5uRakyYmEX/kEaejkzR2/NJxApcGku3+bMxqMovHH3zcrdePjYWpU+HDD817vKAgKFLErUs478oVeO450+ewSxenoxER8Vq2J919+vTh3//+Nw/8vuN4+vRpOnfuzOrVq1O8qLso6fZxly7BihVmB3ztWqhQwSTgTZrA72cMxDclWolM2TGFoIgghtYcSvfK3cmYIVnnvv+WywWzZ8Pw4eaP1YgRZsiNz0lMhJdegvz5zbsLERFJMduTbk+mpDsduXYN1qwxO+DLl8Njj92qAy9fXgcxfciZX8/QcVlHrsRdYW7TuRR/tLjbrp2YaM4RDhliykfef99sAvuk7dtNn0PLgvBwyJzZ6YhERLyabUl37969mTBhAgEBAX+56LJly1K8qLso6U6nEhJg69ZbdeBxcbc6odSoAfff73SEkgKWZTFv3zzeDn+b3v/qzTvV3+G+jO4ZJWBZ5ocm//d/kDWrSbZr13bLpT3PkSPmC922DYYONSdC79NIBhGR1LIt6f7222959tlniYyMvGOBDBkyULNmzRQv6i4ZMmRg6NCh+Pv74+/v73Q44gTLgsOHzQ740qXmwFiDBiYJf/FF0zVFPN7FmIu88fUbHP/5OHObzqV8vvJuu3ZkJAwaBL/+CqNGQePGPvqDkbNnTVH6smVm8E2PHjqILCLiBhEREURERDBs2DD7yktcLheBgYEEBweneAE7aadb7nD+vEk6QkJg0ybTqaFpUwgI8PKpJr5r6ZGldF3RlcCnAxleazhZ7suS8otZlhkZ+dhj7Pw+D//3f3DihKndbtXKR7tRXroEH3wAn38Ob7wBAwZArlxORyUi4nNsr+muXr06a9euJUuWVPxDaBMl3fK3rlwxvYmXLoXQUChR4lYZSokSTkeX7l2OvUyfVX3Y+P1G5jSdQ7XHq6X8YufPw3/+A3Pn4rpwiZ+u5uCF7FvoOfRhOnb00XLmmBjT43DcOHjlFRg82ByYFBERW9iedLdr144jR47QuHFjsmXLdnPRfv36pXhRd1HSLfcsPt7UGSxdanbB//GPWwcxn3sOMrqvM4YkbdWJVby+/HUaFW3EmBfGkD1z9uRf5OpV8/9z7lxzaLBZMzY+GcjL46uzvOQAnrV2kGlNuBm+5EtiY+Gzz0xheo0apvVK0aJORyUi4vNsT7pvzJr/8zCKoUOHpnhRd1HSLSliWfDtt7fqwC9eNOUnTZua03VZszodoc/6NfZX3gp/izWn1vBZ48+o82Sd5F0gMRG++QbmzDEz2p97DgIDSQxoQtDobMyeDYsWQeWKiWbKTfbs5rm+UMR99SpMn25GuT/7rKnfrlDB6ahERNIN25PuhQsX0qJFiyQfc4KSbnGLkydvDeTZu9dMSWnaFBo29KEZ4M5bfXI1nZd35sWnXmTMC2PImSXnvb/4+HGzoz1vnjkc+9pr0LYt5MvH5cvmtzExph1gnjy/v+bqVfD3N4OVhgyx40tKG7/9BtOmwdixULWqGZdZ3n0HTUVE5N7YnnSXL1+e3bt3J/mYE5R0i9tdvHhrIM/69VCx4q2BPP/8p9PReaXouGj6r+7PyuMrmREwg3pF6t3bC3/5BRYsMMn2yZPQpo1JtsuVu7lzfeCAmf3SsKHZAL6jW+SPP5rd8JEj4dVX3fuF2S06GiZPhvHjzZuH997z0Qk+IiLeIbV5512bt4aGhrJy5UrOnTtHr169bi4SHR3N/eqDLL4qd25o3978unoVVq82O+AjRphpKjfqwP+Q+MndrftuHZ2WdeL5J55nf9f9PJj1wb9/wfXr5vDr3Lnme1+vnuk5/cILd2TU//0vdOtmctK75tOPPQZffw3PPw+PP25qoD3d5cswcSJ88on5+iMioGRJp6MSEZFUuutO9969e9m9ezdDhgxhxIgRN5PunDlzUqtWLR7ygB+7a6db0kxCAmzefGsgT0LCrU4ofn4aPvInv8X/xjur32HZsWV82uhTGhRtcPcnWxbs3m0S7fnzoUgRs6P9yit/Wd7jcpk8fOFCWLz4HistVq+Gdu1gwwYoVizlX5idfv7ZdCOZPNmUxAwa5LmxioikQ7aXl8THx5PZQ/ttKekWR1gWHDx46yDmqVOmvqFJE7MzmT0FnTh8SGRUJB1COlDjnzUYX288Dz1wlzfo8fFmN3fOHFO3HBhoEuMiRe567UuXTL9tMPn5o48mI7AZM2D0aNiyJZkvtEF8vPkztGvXrV8HD0KLFjBwIDz1lLPxiYjIHWxPuo8dO8agQYM4dOgQ165du7noqVOnUryouyjpFo9w9qwZyLN0qRlPX6PGrYE8efM6HV2aiYmPYdDaQSw6vIhpDacRUDzg71/w9tuwY4cp3alePcm2jbt3Q7NmJi8dNSqFP1x4913YuBHWrEm7LjVXr8K+fbcn2EeOmMS6fHnTgaRCBXjmGciZjMOlIiKSpmxPuqtVq8awYcPo168fy5cvZ9asWSQkJDBixIgUL+ouSrrF41y+bAbxhISY2uTSpW+VofhwqcDG7zfSfml7qhSqwoQXJ/DwAw///QvCwqBzZ5NJ38Ou8xdfQN++pvIiVY2TEhPNVvl995lhOnbU5btcZjrkxo0mwT51ytRk30iuK1QwByJ/n3sgIiLewfaku0KFCuzatYuyZcuyf//+2x5zWoYMGRg6dCj+/v74+/s7HY7I7eLizCG4pUvNTnjOnCb5fu01n5mIefX6Vd5b9x5fHviSKQ2n0LRE06Rf9NNPZof3P/+BWrX+9qnXr0P//uYs5JIlbmrece2aOVhZt66ZD+9O331n+hdmzWr+W6GCeePloSV6IiKStIiICCIiIhg2bJi9SXfVqlX55ptvaN68ObVr1yZ//vwMHDiQo0ePpnhRd9FOt3iNxETYudMk4J99Zg4Nvvii01GlypYzW2gf0p4K+Sowsf5EHs12D3XSiYnQoIFpxThy5N8+9cIFs6udLZvJz916dvvCBahSxfTvfu0191zzxnb8wIHQp4+mnIqI+Bjbd7q3b99OyZIluXz5MoMHD+bKlSsMGDCA5557LsWLuouSbvFKmzeb5tIffWQODnqZa9evMWT9EL7Y/wWT6k/i5VIv3/uLx441IyM3bPiLptq37NgBL79s8uGgIMiUKfVx3+HwYdP/+ssvk9xx/1u//mp6F+7aZU53PvOM20IUERHPYXvS7cmUdIvXOnzY7HT36GEOFHpJz+9tZ7fRPqQ9ZfKUYUqDKeT+R+57f/G335qveft2KFz4rk/7/HNz3nH6dFONY6v1602Nd2Rkykp+Nm0yTcJffNG8oVCdtoiIz7I96d6xYwfvv/8+UVFRuFyum4vu27cvxYu6i5Ju8WrnzplkrU4dk7B5cDlCnCuOoIggZu2ZxSf1P6FF6WSeZoyONvXNI0dCy5Z/+ZT4eFOVsW6dqd9Os3kws2ebDipbt5rhSPfC5TJfy7Rp8Omn5rCsiIj4NNuT7mLFivHRRx9RpkwZMv4hKXjiiSdSvKi7KOkWr/fLLyZhK1DAJH9Zsjgd0R12nt/Ja0tfo/gjxZnacCp5s6egDWL79qZGZObMv/z0Dz9A8+Ym550714HOeYMHw9q15tcDD/z9c28clsye3fQYz5cvbWIUERFHpUnLwE2bNqV4ATsp6RafEBtrkrhffzUjFj2kV3OcK44RG0YwY9cMPq73Ma3KtCJDSspggoNh2DBT8/yPf9zx6c2bzYHJN94wkyYd2fC3LPP/ICHB1GXfLQgdlhQRSbdsT7rDw8NZsGABderUuTmZMkOGDDRr1izFi7qLkm7xGQkJpr5761bT5/uxxxwNZ9cPu2i/tD2FHyrMtIbTyJcjhbu5p07Bv/4F4eF3zGu3LFOZMWQIzJplhno6KjbWlPrUqAHvv3/753RYUkQk3Utt3pnkTLc5c+Zw9OhRXC7XbeUlnpB0i/iMTJlgyhQzarFaNTNYp2jRNA8jPiGe9795nyk7pjCu3jjalm2bst1tME22W7c229d/SrhdLujVy5xf3LTJkS/1TlmzmpaOVaqYaZGdOpnHbxyWrF/fHAbVYUkREUmBJJPunTt3cuTIkZT/wysi9yZDBnjvPbPLXaOGGahTqVKaLb/17Fbe/PpNCuYsyJ4395A/R/7UXXDIEDNtsnfv2x7+9VdTTpIxI2zZ4jHVNMajj8KKFeDnZ+rst241hyWnT4fGjZ2OTkREvFiSSXfVqlU5dOgQpUuXTot4RKRzZ8iTxwyRmTfP9iE6F2Iu8O6ad1l1YhWj645O3e72DWvXmhORu3ff1g7xu++gUSMzEHL8eDON3eMUKwb//a+ZWFmzpvkadFhSRERSKcma7hIlSnDy5EkKFy5Mlt87K6hloEgasHmIjivRxbSd0xgWOYzApwMZ6j+UnFncsO188aIpJ5k1yySuv9uyxQy8GTgQevZM/TK2+/57KFhQhyVFRARIg4OUUVFRf/m4WgaKpIEbQ3S6d4f+/d02RGfj9xvpvrI7jzzwCBPrT6R0Hjf9JMuyICAASpeGDz+8+fD8+abKZPZss4EvIiLibTSR0nvDF7k3bhyi80P0DwxYM4CIqAg+qvsRLUq3cO95jU8+MW31Nm6EzJmxLBg+3Gx6L18OZcu6bykREZG0lNq8Uz83FfF0BQrAhg2mc0abNhAXl+xLXE+4zrgt4yg7tSwFchTgcPfDtCzT0r0J9549ZrLj/PmQOTOxsabpR2ioOY+ohFtERNIzj026jxw5QteuXWnRogUz7zLFTiTdeOgh0+v6+nVTn3Hlyj2/dP136yn/aXnCToaxqeMmPqjzAdkzZ3dvfDExpj3g+PHw1FNcuAC1a5vWgOvXO952XERExHEeX16SmJhIq1atWLhw4R2fU3mJpDsJCeYU4pYtSQ7ROXvlLG+Hv82Ws1sYX288L5V4yb7Wn6+/bnbg587l0CHToaRtWzOIUucQRUTEF3h8eUnHjh3JmzcvZf/0s+VVq1ZRokQJihYtyod/OHD1R8uXL6dhw4a0atXK7jBFvEOmTDB5smkDUrUqHDt2x1PiE+L5cOOHPDPtGYo+XJTD3Q/TrGQz+xLuhQshIgImTyY8HPz9TbI9YoQSbhERkRts3+n+5ptvyJ49O4GBgezfvx+AhIQEihcvzpo1ayhQoACVKlVi/vz57Ny5k127dtG/f3/y5781mKNJkyaEhITcGbx2uiU9++wzGDz4tiE64SfD6Rnak6IPF+XjFz+myMNF7I0hKgoqV4aVK5m6oyLDh5sW19Wr27usiIhIWrN9DHxq+fn53dF2cPv27RQpUuRm28FWrVoREhLCu+++S7vf+xFHRkayePFiYmNjqVWrlt1hinifzp0hb15o0ICfpn5Et8Rl7P5hNxNenEBA8QD713e5oG1bEt/qT78vKhIWZpqWPPWU/UuLiIh4G0fmwZ07d45ChQrd/LhgwYJs27bttufUrFmTmjVrJnmtoKCgm7/39/fH39/fXWGKeLzY+nX5cthL1O/QgXY9m/Kf4YfIel/WtFl8+HBcWf7BSxve4lqcKTPPlSttlhYREbFbREQEERERbrueI0m3O2tL/5h0i6QnK46toPeq3pTNW5a6a8Jp2qIT5PrErUN07ioigoRpM3jh0d0U9cvIpElw//32LikiIpKW/ryZO2zYsFRdz5Gku0CBApw5c+bmx2fOnKFgwYJOhCLidU79coo+q/pw5H9HmNRgEi8WedF8YvNmM0Tn/HkYN86+U4yXLhHXsh2dEj6nUefH6NvX/hxfRETE2znSW6BixYocP36cqKgo4uPjWbBgAY0bN07RtYKCgty69S/iqa5dv8bQ9UOpNKMSVQpWYX/X/bcSbjBDdL75BnbvTvEQnSRZFufqd+bz6Ba0mFWffv2UcIuIiG+LiIhwS2WF7d1LWrduTWRkJJcuXSJPnjwMHz6cDh06EBoaSp8+fUhISKBTp04MHDgw2ddW9xJJDyzLIuRoCH3D+lIpfyXGvjCWQg8WuvsLYmNNk+zLl2HJEsiZ001xQFjTqRQInUHixi2Uq5zFLdcVERHxBqnNOz1+OM7fUdItvu7YpWP0XtWb05dPM7H+RGo/WfveXvjHITorV0K+fKmKIy4ORrTYz1srnyd+3Sby+hVL1fVERES8jccPx7GbykvEF8XExzBo7SCqzqxK7cK12fPmnntPuOHWEJ3mzaFatb8conOv/vc/aPj8NbpEtCbbpDFKuEVEJF3xmvISO2mnW3yNZVksOrSIt8Lfwu+ffoypO4b8OfIn/cK/M3MmvPcehISYQTbJcPQoNGwI83J247nil8kQ/B8VcYuISLrk8cNxROTeHL54mJ6hPbkQc4Evmn1BjX/WcM+FO3WCPHlM9jx3LtSvf08vW7cOWreG4FeWUGXlKpi2Wwm3iIhICnl9eYmIt4uOi6Z/eH9qzK5B4+KN2fXGLvcl3DcEBJid7g4dTOKdhM8+Mwn3lAheAwAAEutJREFU0olnqP3fNyE4GB580L0xiYiIpCNev9MdFBSkSZTilSzLYv6B+QxYPYA6T9bhQNcD5M2e174Fq1aF9etNL+8ffoABA+7YuU5IgHffNfn5NxEJFHvzVejTB557zr64REREPJi7JlOqplvEAft/2k+P0B5Ex0UzqcEkqhaqmnaLnztnEu/atW8bovPzz9Cxo+k0uHgxPDx5BEREQHi4OZgpIiKSjqX77iUi3uRy7GV6h/am9tzatCzdkh2v70jbhBtuH6LTujVWbBwLF0KZMlC4sMmxHz68yXQ/mTdPCbeIiIgbKOkWSQMJiQlM/3Y6JSaV4Or1qxzsdpBulbqRKaNDCW2uXBAWxtVoF/sKNWDs0CssXgzjx0PmmF/McJ0ZMyB/KjuniIiICOADNd0ini4yKpLeq3qTI0sOVrZdSYV8FZwOicREmD47K0O3L+Trwj3ZEl+DjP8MBesx6NIFGjc2hy9FRETELbw+6dZBSvFUpy+fpv/q/mw9u5UxdcfQonQLMnhAy72jR+H118HlgvUbMlGq5GR4/30zRKdtWzNIZ948p8MUERHxCDpIiQ5SimeKiY/hw00fMnnHZHpV7kX/av3Jdn82p8Pi+nUYM8aUkAwZAt26/alce+ZMePtt2LwZSpZ0LE4RERFPlNq8U0m3iJvcaAH4zpp38Hvcjw/rfEihBws5HRYAO3dC586mRHvaNHj88bs8MSFBBydFRET+giZSiniAned30ntVb2Jdscx/eT7VH6/udEgAXL1qdrW/+ALGjoU2bZIYKqmEW0RExBZKukVS4cfffmTQ2kGEnghl1POjaP9MezJm8IymQGvXmjORVarA/v2QO7fTEYmIiKRfXp906yClOCHOFceEbRMYvWk0Hct35GiPo+TMktPpsAD45Rd46y2TdE+dCg0aOB2RiIiI99JBSlTTLWnPsiyWHV3GW+FvUSp3Kca+MJaijxR1OiwALAu++gp69YLmzWHUKMiRw+moREREfINqukXSyMELB+kT1odzV84xpeEUXnjqBadDuuncOeje3XT7W7QIqqbxkEsRERH5e55RfCriwX6+9jM9V/ak1pxaNC7WmL1v7vWYhDsxEaZPh/Ll4ZlnzGR3JdwiIiKeRzvdInfhSnTx6c5PGRY5jFdKvcKh7od4NNujTod107Fj5qBkbCysWwdlyjgdkYiIiNyNkm6Rv7D21Fp6r+pNnn/kYW3gWsrmLet0SDddv27a/330EQweDD16qNOfiIiIp1PSLfIHJ38+ydur32bvj3sZ+8JYmpZo6hGj22/YtQs6dYI8eczAmyeecDoiERERuRdeX9MdFBTkljYukr5Fx0UzcM1A/vXZv6icvzKHuh/ipZIveUzCffUqDBgA9etDv36wapUSbhERkbQQERFBUFBQqq+jloGSriVaiczdO5f/W/d/1HmyDv+u/W/y58jvdFi3WbfO1G5Xrgwff2x2uUVERCRtqWWgSAptObOFXqt6kSlDJha3WMy/Cv7L6ZBu88sv0L8/hIfDlCnQqJHTEYmIiEhKeX15iUhynb1yllcXv0rz/zanV+VebO602eMS7q++Mt1IsmaFAweUcIuIiHg77XRLunHt+jXGbhnL+K3j6VqxK0d7HCV75uxOh3Wb8+dNN5LDh2HhQqhWzemIRERExB200y0+z7IsFh1aRMnJJdnz4x52vr6Tkc+P9KiE27JgxgwoVw5KlzZDbpRwi4iI+A7tdItP2/PjHvqs6sPP135mVpNZ1Cpcy+mQ7nDiBLz+OsTEwNq18PTTTkckIiIi7qadbvFJF2Mu8sbyN6j3RT1alWnFrjd2eVzC7XLB6NHw3HPQuDFs2aKEW0RExFd5/U53UFAQ/v7++Pv7Ox2KeID4hHgmb5/M+xvfp23ZthzpfoSHHnjI6bDusHu3GXLz6KOwYwcULux0RCIiIvJXIiIi3DITRn26xWeEHg+lb1hfnsj1BOPrjadk7pJOh3SH6GgYNQpmzTK73IGB4CHzd0RERORvqE+3pHtH/3eUfuH9OH7pOOPrjadB0QYeM0kSIC7OTJAMDjb/DQiAffsgb16nIxMREZG0op1u8VqXYy8zInIEc/bOYWD1gfT8V08yZ8rsdFgAJCTAhg0m0V682PTcbtMGmjeHRx5xOjoRERFJLu10S7qTkJjA57s/Z/D6wQQUC+Bgt4Pkze78trFlwa5dJtH+8kszrr1NG9izBwoVcjo6ERERcZKSbvEqG05voPeq3mTPnJ2VbVdSIV8Fp0Pi+HGTaAcHm44kbdrA6tVQqpTTkYmIiIinUNItXuH05dMMWDOArWe3MrrOaFqUbuFo3fb587BggUm0z5yBli1h7lyoXFkHI0VEROROqukWjxYTH8PoTaOZtGMSvSr3on+1/mS7P5sjsVy+DF99ZRLtXbugaVOzq12rFtynt68iIiI+TTXd4pMsy2L+gfm8u+Zdqj1ejd1v7ObxBx9P8ziuXYOvvzaJ9rp1UKcOdOsGDRrAAw+keTgiIiLipZR0i8fZeX4nvVf1JtYVS/DLwVR/vHqaru9ymQQ7OBhCQuDZZ6FtW9NbO1euNA1FREREfITKS8Rj/PjbjwxaO4jQE6GMen4Ur5V7jUwZM6XJ2pYF27aZRHvhQnj8cVM60rIl5MuXJiGIiIiIB1N5iXi9OFccE7ZNYPSm0XQs35GjPY6SM0vONFn70KFbnUcyZzaJ9jffQNGiabK8iIiIpBNen3QHBQXh7++Pv7+/06FIMlmWxfJjy+kX1o9SuUuxpdMWij5if7b7/femj3ZwMFy8CK1bw6JFUL68Oo+IiIjI7SIiIoiIiEj1dVReIo44eOEgfcP6cvbKWcbXG0+9IvVsXe9//zOJdXAwHDwIL79sdrX9/CBT2lSwiIiIiBdLbd6ppFvS1M/Xfmbo+qF8efBLBtcYTNeKXbk/0/22rPXbb7BsmUm0v/kG6tc3iXa9epAliy1LioiIiI9STbd4BVeii093fsqwyGE0L9Wcw90P82i2R92+Tnw8hIebRHvlSqha1STa8+dDjhxuX05ERETknijpFtttOL2Bbiu6kecfeVgTuIan8z7t1usnJsLGjSbRXrQISpQwifaECZA7t1uXEhEREUkRlZeI7UKOhJBgJfBSiZfcNrrdsmDvXpNoz58PDz1kEu1WreCJJ9yyhIiIiMhNqun23vAlBU6eNEl2cDBcvWoS7datoWxZpyMTERERX6ak23vDl3v000+wYIFJtE+dghYtTLJdpYpa/ImIiEjaUNLtveHL37hyBZYsMYn2tm0QEGBGsdeuDffb0+xERERE5K6UdHtv+PInsbEQGmoS7fBw8Pc3O9oBAZAtm9PRiYiISHqmpNt7wxcgIQEiIkyivWQJlCtnEu2XX4aHH3Y6OhERERFDSbf3hp9uWRbs3GkS7QULIF8+k2i3bAkFCzodnYiIiMidNBxHvMbRoybRDg42H7dpA+vWmb7aIiIiIr5MSbfYbvFiGDUKzp83fbSDg6FiRXUeERERkfRD5SViu02bzCFJf3/IlMnpaERERESSL7V5Z0Y3xuJ2MTExVKpUiRUrVjgdiqRCtWqm1Z8SbhEREUmvPDrpHj16NC1btnQ6DBH5CxEREU6HIJIu6d4T8U62J90dO3Ykb968lP3TnO5Vq1ZRokQJihYtyocffnjH61avXk2pUqXInTu33SGKSAroH34RZ+jeE/FOtifdHTp0YNWqVbc9lpCQQI8ePVi1ahWHDh1i/vz5HD58mHnz5tG3b1/Onz9PZGQkW7duJTg4mBkzZqh2+x548l/EaR2bHeu565qpuU5KXpuc13jynyFP5snfN1+499x1XU++91K6Rnrnyd8z3XvuuY4v3Xu2J91+fn489NBDtz22fft2ihQpwhNPPMH9999Pq1atCAkJoV27dowfP578+fMzcuRIxo8fT5s2bejSpQsZ1OoiSfrLx971nP6LJ6WvVdJtP0/+vvnCveeu63ryvZfSNdI7T/6e6d5zz3V86d5Lk+4lUVFRBAQEsH//fgAWLVpEWFgYM2bMAOCLL75g27ZtTJw4MVnXVSIuIiIiImnF64bjuCtZVsmJiIiIiHgDR7qXFChQgDNnztz8+MyZMxTU/G8RERER8VGOJN0VK1bk+PHjREVFER8fz4IFC2jcuLEToYiIiIiI2M72pLt169ZUrVqVY8eOUahQIWbNmsV9993HpEmTqFevHqVKlaJly5aULFnS7lBERERERBzh1WPgRURERES8gUdPpEyukJAQunTpQqtWrVi9erXT4YikG0eOHKFr1660aNGCmTNnOh2OSLoSExNDpUqVWLFihdOhiKQrERER+Pn50bVrVyIjI5N8vk8l3U2aNGH69OlMmzaNBQsWOB2OSLpRokQJpk6dypdffklYWJjT4YikK6NHj6Zly5ZOhyGS7mTMmJEcOXIQFxd3Tw1BPD7pTskY+ZEjR9KjR4+0DFPE5yT33lu+fDkNGzakVatWaR2qiE9Jzr23evVqSpUqRe7cuZ0IVcTnJOf+8/PzY+XKlXzwwQcMHTo0yWt7fNKdnDHylmXxzjvvUL9+fZ555hmHIhbxDcm59wACAgIIDQ1lzpw5ToQr4jOSc+9FRkaydetWgoODmTFjhuZXiKRScu6/G3NncuXKRVxcXJLXdmQ4TnL4+fkRFRV122N/HCMP3Bwjv2bNGtauXcuVK1c4ceIEb7zxRtoHLOIjknPvXbhwgcWLFxMbG0utWrXSPlgRH5Kce2/kyJEAzJkzh9y5c2tSs0gqJef+O3LkCGFhYVy+fJmePXsmeW2PT7r/yrlz5yhUqNDNjwsWLHhzjPy9fNEikjJ3u/dq1qxJzZo1HYxMxLfd7d674bXXXnMiLJF04W7337vvvstLL710z9fx+PKSv6J38iLO0L0n4gzdeyLOcdf955VJt8bIizhD956IM3TviTjHXfefVybdGiMv4gzdeyLO0L0n4hx33X8en3RrjLyIM3TviThD956Ic+y8/zQGXkRERETEZh6/0y0iIiIi4u2UdIuIiIiI2ExJt4iIiIiIzZR0i4iIiIjYTEm3iIiIiIjNlHSLiIiIiNhMSbeIiIiIiM2UdIuI+IC9e/cSGhp6189/++239O7dO1VrnD9/nldeeSVV1xARSa80HEdExAfMnj2bb7/9lokTJ97xOZfLxX333edAVCIicoN2ukVEPEBUVBQlSpSgQ4cOFC9enLZt2xIeHk61atUoVqwYO3bsAGD79u1UrVqVChUqUK1aNY4dO0Z8fDxDhgxhwYIFlC9fnoULFxIUFES7du2oXr06gYGBREZGEhAQAECfPn0YMWIEAGFhYdSsWfOOeCIjIylfvjzly5enQoUKxMTEEBUVRdmyZQHo3Lnzzc/nyZPn5vXGjBlD5cqVKVeuHEFBQWnwnRMR8RKWiIg47rvvvrPuu+8+68CBA1ZiYqL17LPPWh07drQsy7JCQkKspk2bWpZlWVeuXLFcLpdlWZa1evVq6+WXX7Ysy7Jmz55t9ezZ8+b1hg4dalWsWNGKjY21LMuy1q9fbzVq1MiyLMu6evWqVbp0aWvdunVW8eLFrVOnTt0RT0BAgLV582bLsiwrJibGcrlc1nfffWeVKVPmtudFRUVZpUqVsr7//nsrLCzM6tKli2VZlpWQkGA1atTI2rBhg9u+RyIi3kw/bxQR8RCFCxemdOnSAJQuXZo6deoAUKZMGaKiogC4fPkygYGBnDhxggwZMuByuQCwLAvrD9WCGTJkoHHjxmTJkuWOdR544AFmzJiBn58fEyZMoHDhwnc8p1q1avTt25e2bdvSrFkzChQocMdzYmNjeeWVV5g4cSKFChViwoQJhIeHU758eQBiYmI4ceIEfn5+qfvGiIj4ACXdIiIe4o8JcsaMGcmcOfPN399IrgcPHkzt2rVZsmQJp0+fxt/f/67Xy5Yt210/t2/fPnLnzs25c+f+8vPvvPMOjRo1YsWKFVSrVo2wsLA7Evg333yT5s2b8/zzz998bODAgXTp0iXJr1VEJL1RTbeIiBe5cuUK+fPnB2DWrFk3H8+ZMyfR0dH3dI3Tp08zbtw4du/eTWhoKNu3b7/jOSdPnqR06dIMGDCASpUqcfTo0ds+P3nyZH777TcGDBhw87F69erx+eefExMTA8C5c+e4ePFisr9GERFfpKRbRMRDZMiQ4a4f3/j9gAEDGDhwIBUqVCAhIeHm47Vq1eLQoUM3D1L+1etvfNy5c2fGjh3LY489xsyZM+ncuTPx8fG3rT1hwgTKli1LuXLlyJw5M/Xr17/tmmPHjuXAgQM3D1NOnz6dunXr0qZNG6pUqcLTTz9NixYt+O2339z5LRIR8VpqGSgiIiIiYjPtdIuIiIiI2ExJt4iIiIiIzZR0i4iIiIjYTEm3iIiIiIjNlHSLiIiIiNhMSbeIiIiIiM2UdIuIiIiI2ExJt4iIiIiIzf4fkxKyHOM9bYoAAAAASUVORK5CYII=",
      "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}