{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code Parallelization with mpi4py\n",
    "================\n",
    "\n",
    "In our previous lesson on the topic of parallelization, we used the idea of _code vectorization_ to speed up computation. By thinking about how we could use better data structures (numpy arrays, instead of Python lists), as well as how to do operations on many elements of a vector concurrently, we were able to speed up our simulation of a random walk. In this lesson, we build on those ideas and introduce another way to gain speedups.\n",
    "\n",
    "To do so, we take ideas from high-performance computing (HPC) to take advantage of the fact that modern computers (and supercomputers/computing clusters) often have more than 1 CPU. As such, it makes sense that, if we have a problem where the compute work can be split apart, that we take advantage of those additional CPUs. To do so, we have to develop and understanding of how our computer's hardware works, as well as what kinds of tools we need.\n",
    "\n",
    "General Ideas of Multi-Purpose Parallelization\n",
    "----------------------------\n",
    "\n",
    "We can think of multi-purpose parallelization as the more general problem of how to do parallelization. For example, some problems which fall under this category could include:\n",
    "\n",
    "* Monte Carlo simulation (such as the random walk simulator from the previous lesson)\n",
    "* Solving systems of differential equations\n",
    "* Performing analytics on \"big data\"\n",
    "* Machine learning algorithms\n",
    "\n",
    "Each problem above has its own challenges to being parallelized. Before we discuss one particular set of problems for which parallelization is straightforward, we need to learn a little bit about some jargon associated with compute hardware.\n",
    "\n",
    "**Cores, Nodes, Threads, Processes, Oh My!**\n",
    "\n",
    "These days -- particularly in the realm of building supercomputers -- a big question is \"How much work can you do, and how fast can you do it?\". To answer such questions, people may discuss:\n",
    "\n",
    "* How many **CPU**s are on the machine: The CPU is often called the \"brain\" of the computer, as it is the piece of hardware responsible for controlling the execution of instructions and code. Most modern computers have several CPUs.\n",
    "* The number of **cores**: Each CPU on your machine may have several cores. Each core is capable of receiving instructions from the CPU and performing calculations. (You may have heard the phrase \"dual-core\" or \"quad-core\" when looking at the specs for your CPU.)\n",
    "* How many **processes** may be executed by the machine: A process is a set of instructions with memory specifically reserved for it. Your machine runs many processes while you are surfing the web, writing a document, etc.\n",
    "* How many **threads** each core can support: A thread is also a set of instructions, but threads do not necessarily have their own specific memory allocation.\n",
    "\n",
    "Different CPU architectures support different numbers of threads. Your operating system typically controls the number of processes you as a user may spawn.\n",
    "\n",
    "In the supercomputing world, the idea of a **node** is also important. A node consists of CPUs, memory, and connections to other nodes. (For example, the [Tianhe-2](https://en.wikipedia.org/wiki/Tianhe-2) supercomputer has 16000 nodes, and each has its own set of CPUs and co-processors. More close to home, UNM's Center for Advanced Research Computing [will be installing](http://www.carc.unm.edu/news/2017/03/lanl-donation-adding-to-unm-supercomputing-power.html) a new machine which has 500 nodes.) If we are doing development on our local laptop/desktop, the idea of a node is basically irrelevant.\n",
    "\n",
    "(For more information and perspective on parallelized computation, see [this page](https://computing.llnl.gov/tutorials/parallel_comp/) from Lawrence Livermore National Laboratory.)\n",
    "\n",
    "**Communicating Between Processes - the MPI Standard**\n",
    "\n",
    "When executing parallelized code, it is often the case that we send some part of the problem to one process, another part to a different process, and we expect that those two processes will communicate with one another as the computation proceeds. [MPI](https://en.wikipedia.org/wiki/Message_Passing_Interface) is a standard which specifies a protocol for how that communication should occur. Although understanding the exact specification is \n",
    "beyond the scope of this lesson, we do need to understand a few things about it:\n",
    "\n",
    "* Communication between processes takes place using a ``Communicator``. The communicator \"hooks up\" processes to talk to one another.\n",
    "\n",
    "* To keep track of which process is doing what work, MPI assigns each process a value for its ``rank``. The rank is a unique identifier for that process.\n",
    "\n",
    "* There are communicators which enable _global_ communication, allowing processes to talk to all other processes. Such communication is useful when, for example, a single process is suppose to be \"orchestrating\" the computation, and sending out work to other processes.\n",
    "\n",
    "* ``mpi4py`` is a particular _implementation_ of the MPI standard\n",
    "\n",
    "Evaluating Multiple Function Inputs - the \"Embarrasing/Pleasingly Parallel\" Case\n",
    "-------------------------\n",
    "\n",
    "Because general-purpose parallelization is a broad topic, this lesson is going to focus on a particular sub-topic which often comes up in physics simulations, known as \"embarassingly\" or \"pleasingly\" parallel problems. These problems have the property that the work which needs to be done by any given process can be computed _completely independently of all other processes_. ([Appararently](https://en.wikipedia.org/wiki/Embarrassingly_parallel), the phrase \"embarassingly\" has the same sentiment as in the phrase \"an embarassment of riches\".)\n",
    "\n",
    "Such problems generally take the form of evaluating some function $f(\\mathbf{x})$ on many inputs $\\mathbf{x}$, where $f(\\mathbf{x})$ is independent of the value of some other input $\\mathbf{x}'$. These kinds of problems naturally arise in, say, creating some kind of phase diagram - we are interested in the behavior of some order parameter or other quantity of interest as we sweep the parameters of our simulation.\n",
    "\n",
    "In this lesson, we'll figure out how to parallelize the computation of the mean displacement of a random walker over various values of $p$, the probability the walker moves to the left or right."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Saying Hello to Your Computer Hardware\n",
    "===================\n",
    "\n",
    "As mentioned above, modern laptops typically have more than 1 CPU. Find out how many CPUs your laptop has:\n",
    "\n",
    "* Linux/MacOS: This information should be available through your \"About this Mac/About this Computer\" prompt.\n",
    "* Windows: This information is probably available by looking up the specs for your device.\n",
    "\n",
    "To give an example, information about my Linux laptop is below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAFGCAIAAABYHilOAAAAA3NCSVQICAjb4U/gAAAAGXRFWHRT\nb2Z0d2FyZQBnbm9tZS1zY3JlZW5zaG907wO/PgAAIABJREFUeJzs3XdYFEcbAPB3r+wV7ui9KVWk\nCKLYuyIq9q7R5LMmtsTYNfbee4wlsRM19l4TYxc1ir2gFAWltzuu3+73xyEi3h0Hd3oq7+/Jk+fY\nnZmdXeHdudnZGaIgPw8QQgiZA8PcFUAIocoLQzBCCJkNhmCEEDIbDMEIIWQ2GIIRQshsMAQjhJDZ\nYAhGCCGzwRCMEEJmgyEYIYTMBkMwQgiZDYZghBAyGwzBCCFkNhiCEULIbDAEI4SQ2bDMXQFUOShk\nqlfPVElP1Bmv1JkpVPZrSpRPi3NplYqWigkGAzh8wkLI4PIJroAQWjMd3JmOHiwPX6ZHNYaNo7lr\nj9DHQuB8wegjoSUixf2rikexqsc3VK/igVJXrBymrRPLP5ysHsEOa8p0rmLaSiJkXhiCkYlReVny\nm6cV104qn9yi1SrTFs50rsKpG8WpH83yCjJtyQiZBYZgZCI0pbj9r+z8X4rb/5o88n6I5e7HbdGT\n27QLIbD+2MdC6OPBEIyMRatV8stHJIc3qlOff+JDE1w+N7IPP3oQ9hejLxSGYFRxtEohO79PcmQj\nlZFixmoQJIfbvCe/w2CGg5sZq4FQBWAIRhWkvHdJ9MdMdVqyuStShCA5/I7f8zp9T5Acc9cFIUNh\nCEblps5+U7htnjz2lLkrogXTyVMwYDpZs5m5K4KQQTAEo/KgaenpHYW7ltIyibmrog+nbhvBkDkM\noY25K4JQGTAEI0NRBTmiX8cr4i6YuyIGYdg4Wg5fzK7RyNwVQUgfDMHIIKr4OwXLR6pz0s1dkfJg\nMC16/sTvPAwIwtxVQUg7DMGobLILB8SbptJKhbkrUhGcOq2Fo1bgMzr0ecIQjMog2bu6cP8aoGlz\nV6TiWL6hVhM3MSxtzV0RhErDEIx0o2nx1jnSU9vNXY8SGExW1UCmozvQtColXp36wsB8TDcf66k7\nGLb4Bgf6vGAIRjrQtGjjVNk/e8xdj7cIBq9Nf37nYQxr++JtqsSH4q2zlU/+M6QAhqO79azdTFvn\nj1ZFhMoNQzDSTrxtjvTENnPXogjBZAl/WsWpG/XhLlqlLNw6V3o2xpBymC5e1rP3YI8E+nzglO1I\nC8nh9Z9P/AUAfu+xWuMvABAstmDQTHaNxoaUo36TmD/vf6CQmbR2CFUchmBUmvzy4cJdy8xdi3cY\nlra8dt/pS0EQgv6TDCxNlfQof9VPFZ68GCHTwhCM3qNKfCjaMOWzGv/AadSRYJH607A8qzHdfA0s\nUHHr78Jdy42uF0ImgCEYvUOL8/KXDacVcnNX5D0sDz+Dkrl5G16m5OhG+Y3TFa0RQiaDa8ehd0Qb\nf6EyUz/BgQgWyXB0Z1o7EJa2BI9PMNhUYT6tkFHZb9QZKbRU/F5ikmtQoRx+OWpA06L1U1g+NZh2\nLuWpOEImhiEYFZGd3yeP/YgNQ4InIEMbs0MasP3CmO5+BFPn7x6VmapMuK/47x/ZhQMAoMpIMeTN\nNirjZbnqQxfmi38dbzVtB76+jMwIB6UhAAAqJyNnbBtaUmD6ohlMMrwZt1kPMqwJwS6jS/c9ClnW\ngHBapWB5h9gsOKg/LS3Oy/6hQQVeoRYMns2L7FveXAiZCvYFIwAA8dZZHyP+sjz8bFedsxq/gRPR\nqnzxFwBILju0EQCoEu4rHsXqTys58nvFprAo/HPJFzb3EPq6YAhGoHhw9SN1QRi5pga/4/eaXgLR\nr+Op3AxdyRT3r0iP/V6xQ9ASkWTX0grWDyGjYQiu9GiqcNtcE5ZHCKzJWi2KylYqCncuqnBR7IBa\n3ObdAYDKep07qZPi4fXSKWhKejamYNEQY9Zsll06pHp+t8LZETIG9gVXdvJLhwrWjjNVaexq4ZY/\nrmQIbbLHRFFZrwEACMJ61h52tfAKlqiQ5S0YqHx0Q/MTGVSPHdaU6ehO0xT1Kl529bj6TaLx1SZD\nGlpN/YzeBkSVB4bgSo1Wq3LHtlG/STJJaZzGnYU/zNe8RiG/frJgxSjNdpZfTZs5f1V44AGtkIs3\nT5ed32+SSupiPXM3u3rtj3oIhD6EHRGVmuLGGVPFX36noZYjlhS/xsap15YMaaj5rIq/I7t+osIl\nEyRH+MMi65m7yZrNCBb7vV0Ca4JbnuHAuhUe/NUk5SBULtgKrtRyJ3dRJdw3vhxeu+8E300rtVGd\n+iJ3QntapQQAhoOb7Yqz2gdFKGS0Uk7whUCU3SCgJSJV0mOqIJtgcxj2rixPf2V8XN6cb00w8w5B\n2Cw9wXI36E08hEwFQ3DlpUp8mDupk/HlcBp1shy5VGs/Q+GfiyWHN2o+W3wzkd9xCNC0MvmxMu6C\n8tkd9at4ddbr4hlzGEIbhnNVtlcQO7g+J6wJcHgGVkBx6+/8ZcONn3mH17qfYNBMIwtBqFwwBFde\n4q1zpCeNfQbF8gqymb0HdLxDTMskuT9HagbeEnwhv8sPsvP71a8TyiyW4PI5DaJ50QMNbJZKDv5W\nuNvY2d0IvtB+wzVd54LQx8CcPNnQWf7Q14RWq8TrJ9FyqTGFEHyhzYwYQvcM6ASLzbBxLBp0rFQo\n71+lRbkGFa1SqhIfyc7+qX6dyPILY/AE+pOzA2qrEh8aOzpCqWC5+7A8A4wqBKHywMdxlZQy7iKV\nn21kIRbfTGA4uOlPw/INrfgTM5qWXzmaO7ad/PLhMlIShHDwXIJvWcEDvSW7etzIEhAqFwzBlZTs\n6jEjS2AH1uG17K0/jfLe5dxJnWmZxJgD0ZKCgjVjxdvn6e/tZdg6WvSbYMyBAEB57zItKzSyEIQM\nhyG4UqLUiriLRpVAEIL+k/UP9ZXfOJO/aKippp6QHt9S8Ot4/VGY16Iny9PfmKPQSoXi/lVjSkCo\nXDAEV0bK+DhabNRjWE6d1izvEH2HeHRDtGo0rarI1Dm6yC8fEW+erS8FweD3GG3kURS3/zWyBIQM\nhyG4MlI+vmFkCbyuI/TsVeekFawcZdr4qyE9G6N/sWRORKThKxhppXxy05jsCJULhuDKSPn0tjHZ\n2dVqsasG6txN0+JfJxj/rE+Xwu3z9Q1rIwheZB9jyle/SaQKcowpASHDYQiujIycGIzTtIuevfIr\nRxQPPmJ3Kq2Qi/+YqScBt2mXcs9N/N4BaJO8MYiQITAEVzpUXqZRrTyCwanVStdOWq0S7zL2FYky\nKR5cVdz5V9degm/JCqxrTPmq5MfGZEfIcBiCKx3Vq2fGZGf7hTGs7XXtVVw5WjRH5UcmPbJJz15O\nhM6bhCFUL426RAgZDpfvrHTU6QYtc0lYWJE1GjLsXUGlVL98pnxyUzMtuirpoeL2eTK8udZc0nO7\nTVlX3RSPYtVpyUznKlr3koH1jClcnf7KmOwIGQ5DcKVDZZQVX0iu4JsJ3Fa9i2eeBAAqJ0Mcs0h+\n+TCtkOevGGUzdy+rSvXSJedlKp/dMXmFdZFfP8HvPEzrLqarF8ET0FJxxUqmMlOMqBdC5YAdEZUO\nlZupbzfJtZn5J6/NtyXjLwAwbB0tRy3j9/oZAEAhE62bCDRVKqsi7tKHGz8e5b0rOvcRDJaX7jEb\nZaHyMo1ZCQkhw2EIrnSoAn3DxYRD5rB8aujaa9F1BCciEgBUSY/kt/4ptVeV+MAkNTSQIj5Oz8ty\nTAd3YwqnxfnGZEfIQBiCKx1Kd3BhefhxG3fWn92iz1jNhw+nzlElPzKybuWjkKl1P/pj2DobU7aR\nbw8iZCDsC658FDonqOTUbVvmCm9MN1+Wp7/q5TNVfFypXVSOznXmPxLxhl8IgRUAAKWmpe+m16Fk\nEjr7jTEl02pjJ4BHyBAYgtE7DHcfQ5IxXbxUL5+ps98ATZcM2Xra1x/JR3wHRI7zpaFPATsiKh+K\n1rWn1CM4nRjMt5/eL4rWWfIXh9Z9lRAyIWwFVz4MnV0N6tcGrTqhTksGAIa1Q6kFNwmBpammpjQQ\nyzuEweVTJWb4pWUSUCkBgJaKKQNX6NCGKGudDoRMAkNw5cPi6Nojv3mW32mo/txUZqoq6REAsP3C\nSu1i2jhSGZ90RK3lTyt1vZ0h2bu6cN/qihddVp84QiaBHRGVDiGw1rVLFX9HfvOs/uyFe1dpOhw4\n9duV2sXyqGZ89QxHkBw9yyaps416T5phYWVMdoQMhCG40mHoXm0TAMQbf9EzFaT07J+yCwcAgOni\nRdZtU2ov0zvIJDU0EMs7mGDq/BpHZaQaUzgh1HmjQsiEMARXOgwbBz17qYKc3Kk95JePlHq2RksK\nxDvmi/+YAQDAYAq/X0Cw2KXyckIbf8rv72RQfT17lUbMdsYQ2hg13SVCBsO+4EqHWdaax3RhfsGa\nMcy9q8haLZjOVUClUr64p7j9b9GjNgZTOHQeu3rtDzMy7N1YXkGqhE/0jhxZN0rXLnVasjHvVjCc\nPCucF6FywRBc6TAdPQxJpk5Llh7f8uF2tn8Yt3l3Xbm4rXqLN06teOUMxvKr+eE8QcWUj28ZVbhx\nLzcjZDjsiKh0WB5GrTGsen5PTwOT27iL/r5mU+F3GKRnr+K/c8YUzqzySZ8rosoMQ3Clw7B3JfiW\nFc5Oq5SK2+d17SVIDr/HjxUu3EAs31BOHZ29ECCXKu5eNKr8Tzu0A1VmGIIrI5avzrnQDCG7eEjP\nXl6rPvrXtzcSwWQJB83U89xPdvkorZAbcwiWz0esP0IlYQiujNjVwo3JrnhwVZUSr3M3g2k5cgmQ\nXGMOoQev20j9IV7/KvdlYji6M2wcjSkBIcNhCK6M2Matbgk0rX/pNqabr+WoZSWmkjAZTt0oi64j\n9CRQ3L+iSnxozCHIwDrGZEeoXDAEV0Zs/3CCyzemBNmlw/oawgCcOlHCYQtNG4XZoY2FI5frG3pM\n04W7Vxh5FLJGYyNLQMhwGIIrI4LFJms0MqoISl24a6n+JNwmXYQ/LjdVjwS3QXurcesJUucEFwAg\nv35S9bz0LMblQjBZZCiGYPTpYAiupDgNOxhZguLW3/LrJ/Wn4daPtpm7j+BaGHMggk0K+k0W/rhC\nf/ylC/PFW+cYcyAAYAXV0zOHBkImhyG4kiLDWxgzNE1DvHkmVZCjPw2V/YaWVXz6czKwrvXCw7wO\ng8p89Vm0dS6Vp3dlUgNw60cbWQJC5YIhuJIiSA73g6nOyovKzxat/FHPGpq0UiHePq/oB5LLbfOt\noW1MBpOs2cxq6jarGTEsd78yk0v/3i2/eNCgkvUguZx6peceQuijwheUKy9Oky7Sv3cbWYji4XVx\nzCJB/yla90qP/a5+k6T5bNF+EL/Xz4J+k5QPrsrjLqqe3VGlPgd5iYXsCAbD3oXtG8oOrMOpE8Ww\n1jedUEmq53GFm2cbcxYa3AbtCL7Q+HIQMhyG4MqLXS2c5e6nf2CDIaTHNjOENvzOw0ptp7JeFx5Y\np/nMsHbgdfoeAAg2SdZsRtZsBgBA01R+Fi2XgloFHB7D0q7s+cloSv0micrLIkgOw96VYe2gSonP\nXzSUVimMPAsA4EX2Nb4QhMoFQ3AlRhC86IGiDZONL6lw1zKCxeG1H1hyo3j7PFDINJ8teo3WMgyO\nIAxv6qpTn0uO/SGPPUMXvlsklOnmQxcWlNkfbQgysC7Lt/Q6IAh9bNgXXKlxmnRi2jqZpCjxjvni\nnQuKZxlW3rskjz2t+czy9Oc20zm5WtkodeHu5bnj28v+2Vsy/gKAOvWF8Y/gNHgdBpukHITKBVvB\nlRrBInldR4h/n26S0qRH/6DSkoU/LAQOX7R1bvF2Qb/JFX9Hg1IXrPyxOJozHNw4Ea2Yjp4EQSgT\nH8hjz9BSsfE1Z/nVJMObG18OQuVFFORXfGZr9BWgVYrcn6PUGa9MVSDD3pUMqS87v1/zIzu0sfUU\nLfMOG6gwZpHkyCYAYAhtLAbO4NZvV3LZZlomkR7ZWHjwNz2jMgxhNXUbGdLQmBIQqhjsiKjsCBZp\n0d8E3cHFqKzXxfEXGExBv4oXrkp8KDm2GQAYlrbWs/dwG7QvGX8BgODy+T1HWw5baMyCSZyISIy/\nyFwwBCPg1GlNBjf4GCWzPPwYVnYVzi459oemeSsYOp/p6q0rGadJF16rCg5mIEiOae9ACJULhmAE\nACAYOlf/678Vo0p+kjOisWjdhIpM3UCpNXPDs7yDORGt9Kfldx+pZzVlvRlHMXGlOGQ++DgOAQAw\nnTwteo4R71xg8pJppUJ24YDswgGGgxunVksypCHLL0x/05gS5aoTHsjvXaIlIgAga5cRfwGAYe3A\n9Aoub6BneQfz2uNACGROGIJREV77AYq4C4oHVz9S+VRmqvTUdump7QDAsHVkOldl2LowLG0IFgcA\ngKCpgmwqL1ud+qLUs0Gmq5ch5TOdPcsXgkmu5chlFWs7I2Qq+PuH3iIYwhFLcye2N8mbDvpRORlU\nToaBiQm1YaMdyjnuTfi/6Uw3n3JlQcjksC8YvcOwdRSOXv0xVrswhoGvUKtfPTO8TG6LHtyWPSta\nI4RMBkMweg8ZVE/wv6nmrsV75NdPAk3pT6NOS1YlPTKwQJZvmGDgLKPrhZAJYAhGpfGi+vMivzF3\nLd5Rv0mSnd+rLwVNi7fNLX43Wj+mcxWriRvLng8IoU8CQzDSQjBwOrdBe3PX4h3RljnKxzd17S3c\nv0YzfK1MTFtnq8mbGZa2pqsaQkbBF5SRdrRKKVr1k/zGGXNXpAjBJvndRvGi+pWc0pfKel0Ys1h2\n9ZghJTCs7KxnxDDdfD9aHREqNwzBSCdapRCvHS+7dtzcFSmB5HJCGzOcq4BaqU56onhyy8DZIRi2\njtZTd+AQCPS5wRCM9KLUok3TZf/sMXc9jMJ0qWo9ZSvD0d3cFUGoNAzBqGySw+sLdy0z8HnX54Zd\nvbbVuPW4LjL6PGEIRgZR/PdPwdpxtKTA3BUpH17rfhbfTSFYOP4BfaYwBCNDqdOSC1b/rHpxz9wV\nMQjBEwgHzeQ07mzuiiCkD4ZgVA60SinZu1p6dBOtVpm7Lvqwq9USjlzKdPQwd0UQKgOGYFRuqoQH\novWTVcmPzV0RLQiuhaDPWG5Uv1KTuyP0ecIQjCqCVqtkf++W/LWKEuWauy5vEQSnQXtBv8kMW0dz\nVwUhQ2EIRhVHSwokhzZKT2+nZRLz1oQMa2rR8yeWTw3zVgOh8sIQjIxFFeRIT2yRndtthhYxweDU\nieS1H8T2D//Uh0bIFDAEIxNRyGSXD0v//kv1/O4nOBrDyo7brBu3RS+mc5VPcDiEPhIMwcjEVC+f\nKa6fkF07oX6dYPLCCb6QrN2KW68tO7QxwWKbvHyEPjEMwehjUaclK+MuKB5eVz67Q+VlVrwgkssO\niCBD6pGB9VjewZ/bjPIIGQNDMPoUqMxUZdJjdcozdVoSlZGqzk2n8rM1q3OWRJAchsCGsLZnOrgx\nHD1Y7v6sKtWYHn74ehv6WmEIRuZD0++98czmEiTHfLVByAxw+U5kPgRBWFiZuxIImRO+QYQQQmaD\nIRghhMwGQzBCCJkNhmCEEDIbDMEIIWQ2GIIRQshsMAQjhJDZYAhGCCGzwRCMEEJmgyEYIYTMBkMw\nQgiZDYZghBAyGwzBCCFkNhiCEULIbDAEI4SQ2WAI/jIolcqrV69t3rLF3BVBCJkSrprxZXj27Fnt\niDpMJjM3J9vcdUEImQy2gpE+WVlZW7dtW7Fy5Ucqf9as2ZZW1ouXLNGf7MaNG5OnTGncpKmPr9+Q\nod+XWWx+fv6kyZODgkMcHJ2CQ2r8MnWqWCzWn2X79h2WVtaWVtYvX77UlUalUmnS3Llzp+T2H4YN\n02zX898Pw4ZpEiuVyk2//x7Vpo2HZxUbWzt3D8/GTZpOmDgpLi6uzPNCXx9cuAjpc+vWrR9//MnP\nz+/n0aNNW7JUKhs/Yfz27Tv0J0tJTR0+bPi/Fy4AAJfL9fKqamtroz9LdnZ2q8jWL168YLFYzk5O\nr1+/XrNm7ZkzZ8+dPWNlpX2dpDNnzvxkxAna2Ni4uroU/5ibmyeVSnk8no2Ndck0ACCRSDp17hIb\nGwsADAbDzs6usLDw7t27d+/ezc3N3bRxQ4XrgL5QGILRp0ZR1NFjx6ZPn5GYmKg/5ePHT9p36JCZ\nmdmwYYNxY8c2bdqUxSr7N3bs2HEvXryoERKy5689bq6ur1696tGz56NHjydMnLhh/foP07948WLA\nwEF8Pl8qlapUqgqc0YL58xfMn1/84w/Dhv35564ePbqvXbOmVMpVq1fHxsZyudwVy5d3796Nw+EA\nQGJi4okTJxs2bFCBQ6MvHXZEoE9t0KDB/ft/m5SU9O23/du2baMrmVgs7tO3b2Zm5vDhw04cP96y\nZUtD4m9CQsLBQ4cAYOPGjW6urgDg4eGx/rffAGDPnr9SUlJKpVcqlQMGDiosLNyy+Q+SJI06MQOc\nPHkKAAYNGvjNN3018RcAvLy8RowYHhYW9rGPjj5DGILN7NSpU5ZW1rVqR5Ta3qJlK0sr6z179pTa\nvmLlykaNGjs6Obu4ujVt1nzjxo2lGm4GFpifn//L1KkdOnQMCg5xcnaxsbWrUtWrU+cuZ86c+bCS\n8fHxJbs1Y2L+rEDNi7Xv0L5Hjx4XL15Yu2YNn8fXlWzdb78lJCTUq1dv/rx5BEHoSlbK8eMnaJoO\nDw8PDKxevDEsLKx69QCKoo4fP14q/Zy5c+Pi4n75ZUrr1q0NPIQxNF3SuvpDUCWEHRFfErVaPWPG\nTJIkXVxc8vLy7ty5c+fOnePHT/z1157iJpWB0tPT16xZCwAkSTrY2wNARmbm+fPnz58/v3LlioED\nBpRMzONxa9WqVfyjk5OjMWfRrWvXbl276k9DUdSGDRsBYPy4sQxGORoKcXfvAkDtErXVqFmz5uPH\nTzR7i129em316jUNGzYYO2aM4YcwRlBg4PPnzzdu3NSxQ8eSNwlUaWEr+EtCEMTSpUteJifdv3c3\nOSnxz5gYgUBw/t9/lyxZWrECmUxmRnra48ePHj9+9DI5adDAgQAwefIUkUhUMpm7u8eJ48eL/2vV\nqpUJTkav+/fvZ2ZmEgRx+vSZuvXqObu4url7tIpsvWPHTpqm9WTUjGewd7Avtd3R0REAkpOTi7co\nlcqfx/zM4ZDrf/utXFHeGDNmTBcIBFlZWY2bNBk9+uf7Dx58muOizxaG4C8Jg8EYOmQIn8/XfG7f\nPnrhggUAsGHjRrlcXuEyNR8sLCwWLVpobW0tlUqvXbtmqjpXzO07dwCApunNW7YwGcw6ERFOTk43\nbtwYMXJk8egurcRiEQBYWZb+pm9tbQ0AItG7oWm//rru8eMn48aNq1KliulPQAdfX9+TJ0+EhYUp\nlcrNW7Y0bNgosnXU/gMH9N9X0FcMQ/CXrVu3rgCQn59/9/2v2BVDkmTNsDAASElNNb40Y2RlZQFA\nRETE8/hnV69eOXLk8O3/bsXs3MlgMHbt2n3q1Ckjy8/IyFi4aJGHh8ePo0aZor7lEFqjxoV/zx87\ndrRXr148Hi82NnbAgIHtoqNfvXr1iWuCPgcYgr9sFhYWdnZ2YLqgqSktMzPTJKVVmFgk1lTG1ta2\neGOHDu07d+oEAAcOHtKVUSAQAkB+QX6p7Xl5eQAgFAo0Py5bvlwikUycMKG8fegmQRBEk8aNN23c\n8OJ5/PTp03g87pUrVzt36VrhrzLoy4WP48xM86xfrVZXuASSZBeXY3yBDCYTACg1VWZK42uuh4XA\nAgDyP3h7Piws7MDBg0lJSboyenp6xsbGZmVmldqekZEBAJo+B6lUtnXrNgA4eOjQ0WPHitNoIuCI\nESO5PN7ev3SO6DAhgUAwbuzY0Bo1unXvER8ff+r06U4dO36C46LPB4ZgM9O02irc6szLy0tLSweA\nqm87NI0ssBTNUFytcda0ByrFx8cHABISEmmaLjkiTa6QAwBXd9M1LDR07969t/77r9R2zSvFYaGh\nACCRFEqlUgD4+++/PyzhwsWLJjiB8oiMjPTw8Hj16tWzZ88+8aGR2WFHhJn5+voAgEgkunf/fgWy\nb9i4kaZpJyen0NBQkxRYimYEa0ZGBkWVbheb9kClNG7UiMlkpqena97lLXb27DkAqF27tq6M0dHt\nCIK4ffv248dPijfGxcU9fvyEwWBER0cDgJ2dXUF+3of/aZ5zPrh/7xPPXSWXy7OyMgHAgm/xKY+L\nPgcYgs3Mycmpbt26ADB82PCSjSBdj8iLQ6FUKlu1avWiRYsBYMKE8cUDG8pboH4BAQEkSYrF4vXr\ni6YvyM/P1/QDmPZApTg6Onbr1g0Axo0br+lDAICVK1fduHGDy+V+9923mi0nTpx0c/eo36BBTk6O\nZou3t3eXzp0BYOjQoW/epAFASmrqsOHDAaBXr57u7u7G180Y4ydM3LhxY8mxcfHx8f8bMEAqlTGZ\nzMjIjz7aD31usCPC/JYsXtS2XfS9+/drR9SxsrISCoUikSg/v/QDJQBQq9Uurm5OTk4EQbx+/VrT\nd/n90KFDBg+uWIFlEgqF//vuu42bNk2aPHnxkiUCgSA1NbV3716/rVtn2gN9aOGC+XFxcffu3w+p\nERoSEpKVlZWQkMBgMFYsX148jGzXrl0ikejhw0eXLl8u7kVdtmxp3N27d+/dCw4JcXJySktLU6lU\n1apVW7RwoUkqBgDde/Rks0v/7Xz37XeTJ0/SkyszM/P3339Xq9Xjxk/g8bhCoaVIVCCVygCAIIj5\n8+b5+fmZqoboS4Eh2PzCwsLO//P34sVLLl+5kpWVJZVKbWxsqlWrFhwUVDxvgJOT06hRI+/cuZOQ\nkPDmzRu1Wm1nZ1enTsTgwYObN2tWgQINt3DhAltb2z1//ZWamiqVSoODg2u+LcS0ByrF3t7+n7/P\nLV++4sjRo3fv3uVyuVFRUePGjtE0vTX69fvm4qVLnh4eTRo3Lt5oZ2d34d/z8xcsOHr0WEZGhouL\nS+fOnSZNnCgUCo2sUjGtPeB5ZXVfODg4nDl96sDBg7GxNxITE3Nycthstp+fX/369YYMHlzclYQq\nFZyyHSGEzAb7ghFCyGwwBCOEkNlzO3UuAAAgAElEQVRgCEYIIbPBEIwQQmaDIRghhMwGQzBCCJkN\nhmCEEDIbDMEIIWQ2GIIRQshsMAQjhJDZYAhGCCGzwRCMEEJmgyEYIYTMBkMwQgiZDYZghBAyGwzB\nCCFkNhiCEULIbDAEI4SQ2WAIRgghs8EQjBBCZoMhGCGEzAZDMEIImQ2GYIQQMhsMwQghZDYYghFC\nyGxY5q7A16OgQPTw0cOCApFMJpPJZHK5nKIoc1cKIfQREQTBKIHJZFpZWbm7ubm7uwmFQoNKKMjP\n+9i1/OoRBHHrv/8SE5OcnJzs7e2tra0ZDIbm34bJZGo+m7uOZkDT9IcnrnXjp6+GIckMrL+e8g0p\nwfCrpDVjBUqrcEbD61950DRNUZRaraYoSvM5Ly8vKysrLS3N3d2tQf36NE3rLwFbwSZQIBLFxz+P\njIwUCAQsFovFYgEAQRCaX83K/AuK0FdPE2RpmtZ8cHR09Pb2FovFZ8+erVKliquLi/7sGIKNRRBE\nQkKCn5+fUCjk8/nmrg5C6JMq1dJiMpkcDocgCD8/v5cvX7m5uupvCGMINoE3b9KCgoI0193cdUEI\nmR+Hw3F0dHz27GmZKXFEhLEIglCr1Rh/EULFCILgcDgqlbrMsICtYBNQq9UkSVbax24IoVIYDAZJ\nkiqVquyUn6A2XzmCoCiKJMnKGH/V2Vc3zV24L1FR6jNClRtBECRJqtVqKCssfMIQrEiIGTdg2Jrb\nBWUM0vjS0DQAsNns0iGYkr15eP16vMiUY4PlTzYMaNtj3tX8z2TAMZX/5MK/N5IlVKnPCFVuBEGw\n2WyKoqDig9Jo8X+rxi27Kdb8xLZ08aoe3iSqbZNq1hXsvCC49q6ubg6Cr6zvg6JpFoulpQmsSNg1\nfep/7TZs8ROSpjoYg+/oUcXT2ZJlcIOblqVe27djz+lrj98UUsCwcPYPb97zh++aOH1l/wwIfWYI\ngmCxWBRV8XHBlCxPDPatRw6NsKJkBZlJ9y6d+X322fOdJk7qUc2iAt+52a5RY+ZFlT/fF6F4FHCJ\nTfB2sIrpeijIKl3nbOhqcHK68NH2caN3PBOGduj9c6gHX5mb/vLJgzyKxTRRnYrOsOQgaDDl+SL0\nxTLw76CstpCFi39AdXsmANSs3zyyUczMeYdX7QxaOjSITwCAKu/e0a0xJ2+9KgSuc2jrbwZ1D+c+\nXDVqUXzTOcu/9SEBACR3V45Y/Cp64cLuLpmHJk04WnXKupGBHACgRE9OxcScuP4iV0XwXOr/b9Kw\nRnZMLQXavqsilXdz25rdN5PS8uU0AOkQ2LLPkN51HdgAAKqM82sW7LydLqOBZe3buMf3/Zu4cYhS\nWfie9dq1dEu7cvHWs0wZ8D3q9xg2pFVVDgHaz8XW+JZi+q6hbXcBAHgO2bappyDuj0Xrz95LzJYB\ncH37rlw9yCvnxOwJv119LaGBbRcYOWDciDZVuIU3Zvae/DDq152jAjgAAIU3Z/SclNDj983/c38T\nM2TAHt9le6aG8WhpwvG1y7aefZKrBq5dcPe5iwf4c947uOzptrk7ntl1XLR+VG2rt11OHXsVfVC9\n0XJoQnp3Qe8x58QlSgmctG9VpA1DlX1r1+r1+y8nioDvFtFp2Njv6muuvA7lTY9QpVSeIEMIAjv3\nrfX38msnn/YLqmkB8ud/zVl4mtPq23ED3agXZ7buXL5csGRG0wb+rBv34jJUPu4sAEXKrWdy2/q1\nnEr99SmSDsybfTA/pMO3YwLtoCCH8hUytRc4s73L20pSktSHTzPsu4z43l+gyn1yZueh1Sttls5q\n78ICYFr6Ne0xvLWtkFGYcH7bjk0rHXwXdHJjabI4dP1xeABf/upKzI59WwShHXsN62JHp1zcEbNt\ntXPAou4e7DIPXea10dIKJgDANmr63B6eJBAcO2c29ebplRvJtr2njIuwg0KZoxtJsKyD2g6a0sXO\nkiF+enzNr0tnuQRv/MYjoFUw+9LNG2mq6lXZALLkyw9l9i0auL174kcQhCp576wV5/g9Jy5v5s4q\neJPBcyz9QFDyZP+pdFbNX76tZc3Ucj/WemhPbsDQtb/3VgGAKuP80qkxb+p1DLNkErJHf/w88SC3\n44j5o6tSTw6t+m3aNMuta3s5a45IEO/OX/NBqj29BwZhVFmYqBVcqlCeR7AT/Jf2Mk9Vky+5u/dM\nhnf/Fd+1tGcCVPNkPIlbeuFWVlTjxtVY62Ljsjq5O7MUKbFxBdbh9dxIgBLDM2jxg73HU+zazhzT\n27e44UYX3NBWYGabDu8FQtKxeq3QQA5AsJf6/qg/rj0uaOdiywCC6x7ewB0AAPzc5LFXf733StrR\nTTNPBukQEBYcyIHqXlTc1ZV5NZo1iXBiQoiHLPb6+rj4gq4etmKt51L60BXAtvPw9vYu6gtWAgDw\n3CMa1q3Je5uAX7VB86oAABBURXrh7/k3X0j6egqDWgez5l2IzfimqhtLnnzheq5tg2ZVOEUlaKhF\n6SIQBNWqUyPAkgEBRRtlhVIlRQMAwWTnJyQVgktNP6H2Z66E1kNbce09ve0BFMl/Ld/z3KHLsomt\nHFl03tVtB19XG/7nqI5OTIAQb8a961NOXk7v1l37WdN5N7Wn7+OOXdAIlVTxvwhl1tNUJZ2zZXT/\nLe82CtNElFVwiwD2mitx2W3b2LyJvZljU7dp1VLPo5SZD5MVvIBwD44BBapBexxkWLi4CuBpVqEa\nbBmgyLx5ePvBy09eZ0uAz1UC7aP4YEgeQdo4WkBqloQCYALBs7UhISVfRpX70KajSLu8Y+32s/eT\nM8RgwVcCHSBX0cCwDm8fxp5z7lpG9+72r/69lGXfrI0v9/2cXP8e39S6tG7iN88bd+jcpXOrUEcO\nIYlb0GP8JRkAAHgPXxlFAxCMkvdiOv/K1IFLZYM2LG7nqNZ6aCAAgJY+2TFj/ZOqAzcMqSlkAMjT\n7icq6axVfSJXvSvLKjVPBdpfyFboTI8hGKH3lO8vgpa+up8OZJCnZlQEDezgwdO+K9EDybRwJAlW\ncGQIZ8WFW1lNg65cy7Jv1NT9g6+fNNBAfPikUHuBumpDMNgMoCiaBlCmHF+88rCyUf9h/wuwJ/P/\n27h4r9YsTDYL6LdTSBIEiwlU0Qvc5Tv0h3XR0RFBlHwcR7x7PqfZokzaM3l6jKL1qEmjaziTuVcW\nT9is2c2wqdWpDnfaySvp7WqeO5/uHNnWh0sQ75fA9e6xeG+D22cO/rV7+ZjdMW1mrx3fMGDw0pXd\nVAAAhIV7lex7XLj58KUEPK3eVo1QiHLzJTKaUCXrODQBtOTRltkxb2r8uKW3L6+oc4ECYIePW/1j\ncHHznWAJnbhEaomzJ0p81pUen9ShyuIjdETQ4ocHYm4rLJu2rcYHAHs/N9bplESlXTNv7vvHEga1\nqSuc//c/11Ou5Hm2a+r+QSRj2/u5sU4/i0tVBHqT72/UWmCZFG8evQaPb7u3CndgAkhzrFkgLjtX\nqfpU7NA6Mdh8EqR5Uj3jZGWv7rwE75EDOjVwZgIUZtixIF+zh7Cq2a2Z1ZgjR/9NPJfl1TO69PeI\nokPw3Gp3Glm7bYc/h/5v05/XBzaM9gwK9Xy7l3bu2MTq8pnNB56G/y+AX+q0dB6aljzaMndfbq3x\nKzq6vb11sp2CPVkHk54pHNq9X5CmY6TktMgUpS89Quh9ZYXgwtQnDx9Y04qCzMR7l89ejFd4d5rY\nL5BPAIAwtGtL+9mnly8he0SFuPDVea/TyZqREQ4sAJ5vm2b2F4/+kcYK+KG+44fHIISh3SMdZx5f\nvJToGRnsxJXn5vGCGwTqLrAspJO/Axw+fei8U1NvK7YsRaSC8vzlE3rOpcLYTsFe5N7zm3fX6hPE\nyclkhrWuVToJ1zXIGWIObD/uGl3NliVJzitRbYugbu2cT+1amsoOndRCS3+IMv3SkVh1VR9nCyrn\n/pM8IF2syFLNcEGt78e3iZu6fcSQFz26Nw9xt2RK39zOAuDrObQiae+SfWnOXUb4SJITEgCAIK3d\n3GzqftfR+acDU6eQg7pGePBV2a9ec+p3auLM4DsIIT/u7OXnVZp7lfjsoyM9Po1D6H16YgyDay2A\n+HPrFp0DALbQ2at60yEz2jX2tyrKQ/AC+s6cbLVj19mYlSeVAFz76m39WkY4AACQHi3beR/dnhbR\noZaNtodBBM+/97Qpwp1/nt6+4pgKWFb+7cfVDfLSXWAZ2B4dfh6U+8eBnYv/VQEAi29dta5VOf7e\n9Z1LGZhMpvbJ6Air+qPHt5/16445EyjC0r/bnKbh1qWr7d13zpjsZdvWTj6pBAC2hb1vM1t2USAl\nvTr2DNi1OqVx3wb2Wq6hKvvpxZ1712YqAAihZ51vZoypZ1X6tsOwbTh+028h27ftP/vb3D1KACCt\nXKs1qO/NI9hOWg+tyrj+TzIAHJz2w8GiQpz7/r5tiE+NYWsW2/y66dCvM/YpAfjOYd2D2jdxZjs2\n+773P/P3rdrWuP6c2iU+19GR3oALitDXgqZpJoup1jtTBK6aYSwWi3X8xMlWrVoZuE4JQqgyEIlE\n586di27XVv9kPThNjwlwOBy1Wl3mCiUIoUqCpmnNHLZlpsQxQiagmZVOy6AIhFBlpVKpSLLsQVXY\nCjYBkiRxsWSEUEmaOWzLTIYh2Fg0DRYW/KysLOyIQAhp0DSdlZVlYcEvMypgR4SxaJpycXG5e/eu\nn5+f5qaH3REIVVqapphKpUpOTgoNDaXpMr4f44gIE2Cz2Y8ePc7KznZwcHB0dLS0tGQymZp1jBgM\nBi5ohNDXiqZpiqIoitJ8UKvVBQUFGRkZmZmZ9nZ2gYHVlUql/hIwBJsAQRAsFlsml2VmZObk5IjE\nYs2/R/H/sY8Coa+S5iF8cXuLIAihQGBra+vg6MDlcFUqZZl/+xiCTeZtg/f9Ni82fxH66pUIszRd\n3DI26BE99gWbjOEXHSGENHBEBEIImQ2GYIQQMhsMwQghZDYYghFCyGwwBCOEkNlgCEYIIbPBEIwQ\nQmaDIfhLpkiIGTdg2JrbBfjyHUJfJrO9miG7v+KHhfcCflg2obFtifsALY7bMGbJZbdRv02vJ/wy\n3iyjCh6f3rP/7M3H6RKaZVU1vPU3AzoEWRafE61Iu7Jh3oanEbOWf+tdeuq6MvIWF5F++8SBoxfu\nxGdIaCD4Dj7B9dr1617HnuDau7q6OQhYAECL766fsfJKupwGACA4Nm6+oS269Iysbsn8yBcAIVRh\n5grBtCQrVwHKe3/uf1xrcFDxKruK5GPbL4sBCrMK1SD8Ml7dI5i0hOEVOSDanS9/ffPgzr1LlK6r\nx0ZYEkBJXt04sW/Pof/SabApZ95itPT5/nlzDiRaVG8RPTjQlafKz0x98UxEMRgADNeoMfOiitJR\n4vR0uWPbn4ZEWFJSUUbCjZOHts9NlC+e1cnty7iQCFVC5vrjVIkzRcCysRJdivm3w+x2ziwAACrn\nxu5T6QJ7UlyYXVj0qq8q797RrTEnb70qBK5zaOtvBnUPt2VReTe3rdl9MyktX04D8D3rtWvplnbl\n4q1nmTLge9TvMWxIq6ocAgCAKnh8csfO4zcS89Ucp6Dm3b7t1cCNQwBQBff+2vDn5Scvc+UAnKqd\nJnR7s3hZfNM5y7/1IQGAzr00Z+R69eA1M5rbltVXQ1gEdhsUqPkcUk3w9Pr8+PgcZYQlqc66tHnr\ndUGLYcMLdq27U768xSnkCfvXHki0bjVx7oAaxV8LWr29iq8PTZpwtOqUdSMDixZI4bv4VvO3YwAE\nhYZ5y55NOR6XIuvoJvgyvk8gVPmYKwRThdkScG4zoObVlYcPPW72fQifAEXSib8e8JqN6pCyNia7\nUA0AIH/+15yFpzmtvh030I16cWbrzuXLBUtmtneQpD58muHQ9cfhAXz5qysxO/ZtEYR27DWsix2d\ncnFHzLbVzgGLunuwQZF0YN6cg6KwLkPG+3Eybh2KWTczAxb93NCWSYkTbt1NsWo/fEgNG5AobNy9\n0/xYN+7dzVT5uLEAZCm3k2mXDv4f9AnoQ8uzn1449VDh3Li2MxsAmI6tp/8WRRCKhO27y5v3LdmL\nExeyWIHDu4aUp1uGVklzk2JPxWYxfdt6cjH+IvTZMlcIVhfmSoDn4NemR/iZX/dc6lo9ykF8e//f\n+X59uwQ5bmMr8vNlNHBFd/eeyfDuv+K7lvZMgGqejCdxSy/cymzTFgCAdAgICw7kQHUvKu7qyrwa\nzZpEODEhxEMWe319XHxBVw9byf29x1Pt2s78qbsvBwBC/CzSR6/bey6lTs8qBAAAx6VGRGhR+5G2\nbFyN+VtsXFZHN2eWIu3uc7lN7RAHgy+P9O6yYYvvKADY/j1mdPcrCnuGzROsPa/mMhW8SpWAQ3BV\nQ9uxiZtH9d/89gcyoPfo5s7YC4HQ58tMf5+ULF9CsywFXOvA7q3tphw9/qJh64f773EaTWlkx35p\nyYG0XJkalFlPU5V0zpbR/be8yypME6lBUKIsgrRxtIDULAkFwASCZ2tDQkq+jAJl5qNkBS8g3KPo\nSzpDGFDbDbY9fiWhq1iUrhFhFdIykL360p3sqLZ2uU8e5AhC6riXvfDTW9yA/82ZHZ358u7ZfXtn\nzKXnTe/saXBmfXlpGgAYjJIBmBbfWj5uk7zX/InNrT4oy7Xz5JF1rQi1QpyZeOvkrt0z5tLzpncs\nx4kghD4lc4VgaYEMSBcOA0jPqK4hpzfv2JmWll6l2zg/HkFxBByQFcgoAAAa2MGDp33n/24taKaF\nIwmZ75XGZLOAfjtPJEGwmEDRNF00ied747X0tCUJYXDrMO6yf2MzWzV8EpvCC+pbpewFqN/l5th5\n+Nh5+AQEuYp+mn3i2Iu2w6sbml1PXqbQxZGEe09fy9q7FDeEaWVhgVgq1zovJsfW1cPDjgEAVX2r\n+/ESR6775+qbtj2rsLUlRgiZm5nGBVPyAhmQfDYBwLCJ6NbEIuHSQyq8a2NHFgCD5LNBJpLSwLb3\nc2MpUxKVdq5uxZyt2QZ+KWc7BHiwpc/upCqKDlrw9NZrsK/mztNaACEIatvYOvXcPw/iriRxazT1\n5Vfo1AgGAWqlukIjdT/My/OLrCNQ3Nl7KkFW7hIphVgGwCJZ2BmM0OfKTK1gWi6WA4unaZtxvDv2\n66B85tQ+WEAAAIPkkyAvkKuBEIZ2bWk/+/TyJWSPqBAXvjrvdTpZMzLCwbCDEMIaPdu6/nJkxRp+\n30g/TsatQzsfMMJ+aO1BAqi0ZeB4R0W5ndu3caeaEz7Z4AgsTzp16B7Hu4qDBUOa/vTigX/yLRu3\n8NLZBKYLbq+dsvKh38iFo+pYK8vKS/CD+gxt8mj5gRm/JEe3buDvImDK0x/mAPC0ly55Hf/kcRol\nL8xNe3r95MmXLP/+9ZywNxihz5W5OiLkYgWQfFLTCGfZ1e0ztO7bfQTJ54BaKlLQwOcF9J052WrH\nrrMxK08qAbj21dv6tYxwMLTtzvHq/ssU7raYk+sW7qNJ+4Bm38/q08hW96sKLJdm3cIOrr4jbBnl\nZ2gEpuQFuUlXLpz4q0AFQNp4hXb+uU+nYAu9Lc+3y0kZkpdhHT50waxqBw+evrx7zTElAJBCJ+/w\ncA9eqavAsHCwJ+NPrZ53CgCALXTyCuw0vHPHBvg8DqHPF64d9z5l6v7Jv1xrOG9BFzfsPkUIfWzY\nQgIAAFqRlZQiYSoSz245Iq0/JtIV4y9C6BPAEAwAAMo3536dfuQNIfRq9N3k/9XAt8kQQp8EdkQg\nhJDZ4GSVCCFkNhiCEULIbDAEI4SQ2WAIRgghs8EQjBBCZoOD0j5/tCLt1uHdRy4+yZSzrX3qtOvV\nrXFVPg6bQ+hrgIPSPnfKN6eXLrzm0Xdwl9rufFXGvZPbNt+qMuqXHr7aJxv6CGiVRCSSqwGAYPO4\noGZaWLAJUMnECpYFH+cAQsgY2Ar+vFG517efFfQd7n5tzS+7FUzrkK4/DBspnrPlbKNZHV0/0T+e\nKv3KX3su3HggrlKzdpjgv6M5fRaPCRcdnzbtXpeVUxtYYgxGqOKwL/jzo877b/eaX0+kKADogoeX\nskPbhnDysoWdpi2e2Pj1vkMprpFNieux6Vone/sY2G6Rg3/o6u/esP8PPbxVDG7qsb8f3zx0leLJ\nsyXqT1UJhL5OH6EhRYv/WzVu2U1x0QEEjl5B9aI6d6jv+cm+OX+WaHlG/KNU0i+0qqCM+x4lenE9\n9l5QFA2gzHuZb+3rxCYeAAAQLDaLoIFlXdVFcTZbBZ9+ZWS1TGzZqIf/jTU77br18r/0SlpiEmNF\nQsyUOZerjFo0MhxbxggZ6GP8EVOyPDHYR478vp4NyETpCbEnD6795UHWwmkdv8rZx+iCK7OG/6r4\nbvWcSHvd82CCMvnIimX3ms5fWlVQvmWENFFOnrhv3s9bchy/me/DgXit07fTivTbx/cfvXAnPlNK\nA8vSzT+8Ze9+rb1N8eiOYFsIOExapWAJHfyihw4OEgbwz9xKUpaoCMG1d3V1cxCwANSZp38ZvePl\n+0XYRc9d1rfqB+dOSZJjT5z4J/be89eq0ImrRweJdeelCh4c/n3r8duvpWy76i36DulV34kEANC6\nXdu/iyJ597gpZ7wmr/s5mPvuop08eOzC7WfpEhq49n41m0Z3bV/bhaPtomm/wpGOj1fraHOo0o5O\nGbfLcvTvU2sXT35aeHPukJWib5bPi/5gImdamnrz7Im/r9179ipXTgNL4OwTHNEsqk1Df6sP/lAV\nWXeOx+w9eyO5gCYsXGu06DO4Wy1brX/Opa9wDZ6OK1Z8jicOHL1wJz5DQgPBd/AJrteuX/c69kyD\nrqcxaHnKxT+37r/4JEsBfLeabfoO7Bxm89V3lX60E7Rw9a9WzZ4JEBhaO9Q6c9QfZy+ktult++S9\npeOnzexZldS1zjwAUKInp2JiTlx/kasieC71/zdpWCM7ptZl7Wn5y/PbNu2/8iJPDaS1f7uxU3p4\nkVo3guEr20+b2fPDkPFJsa09rfKfpys9gePV/cch1B9LH6VKmxNJaaSn3fv/dLT0+f55cw4kWlRv\n0X5QoKuAFqUlPE4U0UzTNEh5IcMmhACA26gxAAD2NQCg7yTf9+rqGjVmXpTms039UXN9pEW9FFR+\n3M5VBwvCwxxK34Fp+cuzaxf8meTeuFX0oM6u9vbuHGCydeVVvTm9fNG+jJDuw79xF8Xu2blmPtgu\n6leNq9a+3YB1o2jJs71z5x1Ktgpu3XFINSeOLP3ptdP7V0641WHqtF7VSt26dF9h3W0OJ8OvMC1N\nPLVmccw9TmhU234d/d1sSEqcmfzg+tkds06daTd6fO9Q6xI3eFoSf2DrP/k12v3QxZWZdfvIrqPL\nV9munB7pUKoRoO0Kg84rSZQ4x+jBga48VX5m6otnIorxCTosaXHcpjkb4ry6DZ9ey1b65MQf25cs\n5Syd1cHlKw/Cn+L0CDafwwCVTAmll4531LfOPCiSDsybfTA/pMO3YwLtoCCH8hUydSxrH6U+vuqP\nq9zooVPrubBEGdlcOxaAMlXLRijHyvaO5W60U3k3t63ZfTMpLV9OA5AOgS37DOld913oyTo65X9H\nAQDcei9Z0MGFpfV2UvLSWQY3sT908n79QGcnC45j+HetH+67clN4Aer/9H4TSp6wf+2BRJvISXP+\nV7zcfYMWb2tVnltOWVXSSfX60KQJR6tOWTcykMOydPOy1GymC26vOZNgGz2zb/VS89jTsvi/lu0R\ntZ21tEOVEo1Oho688mcnjjzn1Bs/onOoBQHBTuKnk48eiOs0qWaa9u11oSzy5/vXHkp2iJ46u2/1\nonhbv2nziJ3T5h1duzd8ackVC/VdYboAQEebo6ehbUMqN3b9/JiEkO+Xfh8uvbhz87LfXuSqgCVw\nq9l90oTs35etWGY9f3r0uzlUCX7goCUrGaTmDhvsmBY3/d+4VEWkw3urqei4wvIE7VesHidh/9oD\nidatJs4dUOPtOUKrt4WVdTUfrx8x97KkxBbfYb/OamRVfFxVztMnCq8gZ5IAAHVe/ONCjyC3d+uF\nK1Ku3yt0bP1NxwgPNkCVft2vXN9456Wsg0vJtXq/Qh8tBNOUWq1SqmT5b+KvHdz5kLJsWs+NDXlQ\naul48X+61pn3lD/YezzFru3MMb19i/8U6IIbWpe1b+6fJQK+X3BodW8BAd6axMpCLRtpcTlWti83\nSpL68GmGfZcR3/sLVLlPzuw8tHqlzdJZ7d/eya2bjBrXzpUFBMfagaXjdtLevkSBDOu6/aMuL9z5\nuu/gLhYEg1Ojqce2zbcjfuz9fttA9uLEhSxW8MhuwcIPWr3luuXoqFLFWyLyF4e3Xqfrju3gU/qr\nPV1wZ9/fmbTDpaUj9+UoSJfAZj3+17OuE6kjryr78eMChlc9H02wZDuHhlgdvBb3Wur5VOt2RV1X\nAABKrVQpFG+fGqpKLswne376Qg479McOASXau4RFYKeewefWXjwd38u/xLdrfVe4tHdtDjAwBCsS\nju68SbaYNqQ+6+rSWVuTqvcYMbOmTcbRZeuevmEO7/NjjzvjD+1/1Hxk6LuaEkzybYuXluflSsGm\nSukv7TquMEPHlVSEUScuZLECh3cN0X2Ouq8n6d1z1sJoNQCosq5tWn44s2bkezddedKJDQv+dew/\n/ecoN8mtzbNW3vAcuXh0fZu3DWymtacdXLp9K7WTe1UOiBIfpBNubd1N0L/xmftoIfjljp8H7Cj6\nbOXXcuigftX5hLr0GGQ968y75Dx8b5cmvY5l7VlR7boE39qx6OfkiJaRUZGNAuxIAjheWjaWa2X7\niiEdq9cKDeQABHup74/649rjgnYutkVLNFm5uHt4aKIMXXBb6+2kTZv3SmO7tBozxebw7vWTN2fI\n2Da+ddv/PKlRVe57fyPqglepEnAM9Pyw27dctxxdd7g2Ff02SOXe/Ovv3Ko9uocKCQBQKyRyJU0D\nAMHkECk3n6v4fnU6RYe6CY1ISZQAACAASURBVFXJ/2zbtHqucvqSAdW4WvNS4iwx8H2FbxdvZVjY\nC6AwSyTTsV2p+fHVzrEDd75fKdK76KKlpEi1XDSCXyXIER68TilQB3PfRjk9V/jthdba5igw6DIp\nUq7ezHNuGeXHTj92/AGE/TSsU4QloSBtNa1elmOdJm57Ttx+rQr1/fCLmTr3Vszm/3iNx7bxKLVT\nrv0Kf0PouJIFOakScAiuqm+6bN3Xk+DYunnYAihSj/9xPNku6pfvG73XWcbx7T1lZMGsNbMXJVbP\nuvzQqe+MEfVsSnRwsJxaDe//aP6OKWPu16/Bi7+eHDxkavSnf+L8yX20M3RqO+aHenYcLt/Kzt6a\nq/MxlZ515mmggfjw24/2Ze1JdttJa8IfXjp9/Ngfc48dajJ61pDaNqT7hxstyrWyvXEYFi6uAnia\nVagG2w8703TdTtSll+YkSKfaPX6q3UPfoWgAAELLmZTrlqO7ShULwar0y8ceEDV+aqxZv072aN2o\nBTflAADg8e2cNllycIxo3qCGPQPAq9+gF//NvvJPUt9qARwteYvOseQZEmVsBwAAxzY/Da1j8/b3\nT/nm7JqN/0HJjNouGkEA0KV26LzCRbS3OQwLwcqcV3lsFx9bljrhZS7YNXQtFeiZPCs+yMUK6oOc\nqqxrv8/69U6VAdMHhVuW+h2jZTlar3D3ljquGE0DAINRcg8tvrV83CZ5r/kTm9sCgP7rCQC0/MXh\nlX8muPec2yfwg/UTWfb1B/30/PHMk7c5tUaPalel1PciWpqVlqm0j2hTx0N0P4HKi7twNSmii//X\n/iLoRwvBXMeqPj76RggAQNE686ef3UlVBHprHpQ9vfUa7KPceQTb3s+NdfpZ3NtdmvT2fm6s0ymJ\nSrtm3txS/zQExym41bfBTVoenTZ+95E7PcJb2DI+3Nhc9xFNfg2AYLAZQFG0rm407bcTSCvvcZgC\nF0cS7j1NldEupZsw5bzl6KhShajSr11I4dTo9fa7O+nTc9K0KDUAAMFzsnkdCyDKl6s1o9NZ1q5W\nIM0pVGvPCwyBgwAkWSIlDRwCAKjCrEKw8BJydWwvahHynHz8/d89wefcLj4XptDNmQNxT1Jl7ZxK\nXjRamvokE7g13QQlfnv1XWENrW0OgsVhASikchreRhJaLlYAsEhmiYBJ0xQNTAYBBM+aDw/SCkoN\nOFTmvsoBy6DSoxmV6RfWzdj0tNqgmSOau2j5N2KymdquMOi6kkLSkYR7T1/L2r87R1pZWCCWyouD\nv+7rCZq+5zWHMwL6L4r21FIfde7N7WtPib3rBedc37j+rMvoSHfyXYtLfHfbr+eYXZaMinZmQZvo\ntmfnT9y2/mSdRd3cv8ZxVO+Y+dUMzTrz2SdWrDl47d6D2+e2r9j0gBHWvbUHCYQwtHukY/bxxUt3\n/XPz/sP7ty5fephHCUO7trTP+3v5kphzN+4+fHD7ypmTNzNVoMq6depc7P34hISnD+++KADS0pIk\ntG7Uc8SPjM1lgVwke/ezvZ8bS5mSqLRzdSvmbF309ZCmynr8URLPL7KOQBG398QLaalsbIcAD7b0\n2Z1UhebnoltONa23HP1VKi9V9t3YNyz/Rn5vow+D7+pXvXpA9eoB1QOq2Fq4+tpBxoPnBZo/b3n6\niyyw8Szqz/wgL7DsqlezpBJjEzRnqEy/dz+PVTXMhadje9n/nlzf1g0sFbf3nEqQvbtotPTFib1x\nCst6kT7vdUPqvsJvS3Os6uPjVcXNqeR3PqbQo4oAXsU+zn+bicp7dD0FBFU8hCX+9NhWTgJVRkqB\nmnRrUMdWfHXnoXvZcoUkX0KDJCdH/Cb20Pksx3p1nEvGIloWv3fhpgdeA6dpj78ABMdF6xXm6Lpi\nmnO8s/e9C2I4Wha/77eT+cH/G9ZK26rd8ud75q+8at97xpRR42eOqvlq6+zfYvPetevVeYnJUtLJ\ns2j4HdsuoLoN5L7K++pf/jF7V4vudeYJnn/vaVOEO/88vX3FMRWwrPzbj6sb5KV1WXvLvIQbh089\nyVEAEBauNTqNGhgmINRpWjaWf2V7E2Hb+3uQJ67tOxrS3p9dkMWo3jgstGtL+9mnly8he0SFuPDV\nea/TyZqREQ4Mno0FiB9d+S/ZrW4VgUEVI/hBvQc1frTq0MxfXraLqufnaMGQ56bGP5PW6t+zRs+2\nrr8cWbGG3zfSj5Nx69DOB4ywH1p7kAAfvF9HCHVUqSK/JrQ0+e5rcG6q65Uc0q1FG58TMTs3HhZ0\nD+Mln/rjqrRq75aa1pO2vBzv6GjvC7s2/Ob5baSb6Prug2m2LYeEWRIcC+3bQVTWReNV7/l9y4dL\nDsyY+qpD2wYBTlx5+pMrJ47cSneOnNgn4P1q67nCPdx1H4PjG92u6qW/fp279k10HU++JDn2+MEH\nKu/e0T4lH3GQbhEBnH//uZzasqtPjzH9M5bvWPTTEQBgWQpUN5eNugmC4D6TOrzXsKRyrsScTPfo\n9F1VWWpiIgAAMPhOHo6quLVTVj70G7lwVB1rXVeYA9qvGEEE9Rna5NHyAzN+SY5u3cDfRcCUpz/M\ngdI9Y9opXp3YdDLToXW/qtLXL18BALAtnZytitu5nCqtvxvXxDfMnUsAt/7QWZb3RF5W7+5DLIeQ\nUNsDp7dtP9u/lb+l4lXs3pOZZPVebuYdFfoJfIQQTFg2nLmzofajuXZeurNzqY1Mq+BOPy7opDW9\nTXCnUfNL7WJZh3QaFVJqo2/P6Wt6ls6tbaO+I2qtnqkQwvABg1uu3nlgzWKaEHi1GV0vLEDr7cSB\ntK3Xp/21305tPlCr5s+hBv36AzBsIr5fMLPawUNnruxed1wJQHDtvEIahSio8txyCJ72KlXkhFU5\niRlq0t3NUtdNhOXSZuxk5ebtx1fO2EdzXGp2mzgoumjeC615Wa5tx0yU/7H10JqFcpZNtVYjx/at\nziN0bjegIUdYhg6YM7Pq3r0nzm1ddYQGIOyqNegycnR0PY8P7xu6r7C+Q7DdO0yZydu1++Tx32Nl\nADyngBYDh/du8f5LSoRFULcOnpP3rzsQOK1XQNuxq1rmZ2RK2DYONhxFdlou2DjZ80v9pSre3H9J\ng/zQ4mmHirdVH71harWiLl39V1jHlQSGdfjQBbOqHTx4+vLuNceUAEAKnbzDwz14ZX5dVmXfjU0F\ngDMrfjlTtMm+w/ylvT1L9CIG1yxOzrINCrd9rwCOb59fRjK37tuz+JIEgLT1qdN/8rfvPbD7OuFM\naQgBAJX977zRmzIb//jL/+o4kWZ4AqTKurJp9rorisCobu2bhFdzs+WyQCnOySes7S0++lc0ZD4Y\nghECAABV1vXtK9b/nay08ggMadp/UBttT5Q+Klr+5tbxQyfP33iWrfh/e/cZ0NTVxgH8uZkEQgiy\nl+wgEIYoqFW0DlSkitY90Lp4q3XrW1cd4Lbu1eGoA16tKIqoiFgtLhyADFEQBUEQ2RCQ7Nz3QxiJ\nJAGsGkjP74uS3Jw785ybm5vzb3jMNmj3+iHNfsmMaA5UghGkES6oyH2alpn73nLgUHemuk4+cRGP\nU11dJyRo6ekzaWhEZs2GSjCCIIjaaPzFbgRBkPYLlWAEQRC1QSUYQRBEbVAJRhAEURtUghEEQdQG\n3XDY/uGCd4lRZy7dzizlk5n2PsPGj/a10fThoxDkXwLdlNbeCYtid2xNsJo0a1R3S21RSVrMiWOJ\n1vNXj3X4d6ehIohm0KwLEXjti9j/hd0sFLQ8aQchqXxwMo4+aaplwv7VyxYv23q5tteceV2S/ogr\n+mIZ9giCfD6aVYLFnGdxV289q1YwQosgJ3zZ9Dn7kzkfMwzflyWuSjqz/+DVAgEAzsm4U+7h70at\nKtcNXLN9ue/bcxcLzP36YQ8eFnfcGtzxekpB0e2zl58qHzhRUpl4ctu2szmtWiOZ/atUmxr8THD+\n24TzF9I4KochQv6pL12CcUFx0sVf1i8Onjx5yqTJQbMWrd9z5lHZF6gmMhHr6oHzS148efK6tuUD\nWlLz6sHDtMI6HEBYlV/NdDCpH7AXI5FJGA4kpo2ZIL+845bgjtZTCgtj9x+JTS9TPnSthFuQlv6i\npHXD7Mrs30/T4GfDy39w7sCJ9rUzNM4XrUjqjMiWjVj/tHDOvZC5BwXT9m3wU5USIsy7tHtnWr/N\nO2zobRv+RXr883PPbVr8R4Xx5M32VMj+8D2haBkEeWeWrbpuu/LQYrYW4LVJe5ftfFwrfYrMMLN1\n9uo7xL+vE1PhEYALipOvRkbHP8kuqcMB0zayZ/ccNmWMj+HnPl7+cU+J1ySEfH/wRdMDzksPr+6m\nxXke++f5uMfPi+twkp6N1+DJ04e7SoN+JJynUUeOX0l+yyUbOA+YNHt8LxP5/YNXJ52+mN957M4+\nBkTA617fOh12+f7zEi4AkW5i39V/1oxBlu3sw6Ti3TfK5vH6Jafym09O77d2x+i6uPCIuEd5HBzT\nMXcfMHHW6G6dqLYjZvSP33TyynD3iTYaP26vunzJEsxXFZGtIFDdsurW/i1hycU8HEhMB9+x/wnq\na0HFQFKZePLg6Uc5RdUCAALd0mPAt1O+lRlgkJv865IZFZV8+Qx52Yh1AJDUZF4LD7/64FWlCKOZ\n9fpuxZze9MJbJw6fv/eqSgwUJmvY0lVjbdV+2JGZnfWqXxYLOwPVdsyC2ZKjO54Vcvtjr99ROhu0\nbddJeFW1YDh4XrC3noTHKX2dduf6kdC4W4HLV4x1+jBcviP3lDi/VkBg9Jmz2F8ad0fQNtECwIh4\nHcHWb3qApTb/7eMLYRE/C833LfVmYKKi2F3bzpW4jZk72bLm4Z9h+zdDp21TnGRCsUQlDy6lYN2W\n9DEmAc5J/CX0WJJu91Ezx9npEd6X5L4o19Mmthzw/iUp3X0kpk/wKsv3EgBxafxvRxIMxy0eZ08F\nIFAMDF9Hht6sdh/2/ShzYlnypdPRu/Z22rPWz4jmGDDM6mZU9PNA2fBm5FP6giVYGgOuLCJbQaA6\nkeHYb+zcwZ10Ce9zbp04dXiPkcOWQAuShFvwNLPEIHDebCcdcWX2ncjIfWuK5m1f2IspLRGYcfdR\noz1MSTWZ18Mi5TPkGwheR24KvVDtNnzqEhcD4FRIHHQlhZf3Hr2vFRD8U08zUk1JuVYbS1z9WlQ9\nPrH/zOPX76r5OMj3AQAAUBa96rtoAACLCT9vGW5GElWlRR8Pj0l88x60TD0GT545xquT7HwxBruv\n4cWY9F4upiY6VGOvaYMzzt17rBsPvRZ+zACGOmasLs6GRADo2qu/X5/w9Zui9oa57gh2lX17deye\nUlxXzcUY1va2tnIbSMdl9EwX6X/dnOhZDzZnZ1cIvRl4ztVLL6k9//vDSA8dDNgmtVkroyNTAlf0\nZDQsprgy5fZriuu4LnQMgP82OYtP9Vny/dj6eHufvtKpZC8KCQpjtq8Jq+y/MuRb8dEF27L7bdg1\n1Z4CAHjlnQ3zfhXP2j1DZmpRicKtBwobnMKmY9DCMaNy99m6GAAACN9k0oDIsHVhu9UnAuAzf95D\noBAxAAC28buUtX+nFAr8jGgkI29fy9OR8Tk8D3YrwwOQtvlyJVgaA646Ils2UB0AtCy9vpKGwjha\n8B/eP5j2hjvCQhcAACgmrt26ulABPLq66m1Y/EfEX4XdR1th0ld19/V2oQK4WItS5x+RzZCXwmuf\nRlwpMPBfv2SCQ+O8+FllNaDtyPZwtqNjYPeRKympK8zIKjEc9cN/WHRRZeb1sIvyfQCz7/xlw8xJ\ngFGZRiTgvzy7YWssddDUZTMsJK+uHw/btYv+8/pvDGUaJDB7BA25uzXs7aRZo3QwAtW9n9WJY8ne\nCyYoijSWiIUigaDhiqVIKFZ1cobRXUZO6vbXroSYrCmuXZvClDt4Tyl5X/ZeQoCaSg7TgEFtdt6O\n88uz4q9lCEx9u5uSQVT0/DmHYNvTXtoFkU093PQuJKS8FfRkNBwY/DephWA+0EoLAICob2MEd579\nlZDv8HXnD+NjpRu9IuHwllMFnnM3TmTTiZxejqRHaamlInsLEgCvIDkPNxvOYhBKm16gZOspbhAD\nZcdM46ZTvfuUwoiUhgtYOL+qkgv61tIUP5I+q4seL+VpiYBtrfYPhRrpC54Fq47I1mv+AkHp46iT\nF+5mvi2vA20tIeD2AgVfQRH13X0sIPxZAe9bK/l+mqBjqjBDXliaIZfrDgAAVNtho9iJp7YtzvMe\n6DfEr08Xg4+OTqAYO3fzcKECsG3F6fOPyvYBJD0zS6v6pFCckxxxvcQuaPe0gYZEAKfOhMyUHfGJ\npUOHyrVGNhu0ZJV+1JlfVx4r4ZH1HXp8s3hFHxuF7/83YUtnhH2wLCo6E4xmxTaBpHf5VaKuOg0H\nQkfvKSVioOsUn163OByA4dBv7IygAdb1G4ubunPO9icCADJr7LoxjloYCGrLakHbQbchn5SgY0iH\n92U1QoD6OYpqispFVDNjHQIAAMmk/4K5b3YfOboiOMrFt/8gvwHdbXSb3kLi6qSTRw8mW04LndPH\nkAQAuixfJ+IvD1PKRliYkgTvUl/y9bu7Gcl9YYCp2HrNG8Q5qYqPmeH1Nbg1u08lcWVi+LEkmu/S\noVbSz25EpoUexBdUicH641pEVPtyJZioa9ZyRLYMYcGV7XuihH2C5nzXxZBSnfT79gglLRMIGACu\nICleWYY8DjhgHz5IsfRfsd8r407slctHN16+2HdRyOzu/zC4iqBjprAPkBKWZRUK8Yo/FgX90fSg\n7rsa8YdpiRjFpPvYhd3Htjg/46ELg330G97hwqK4/b8ntXmhO3hPSWN/v+u370HCL3v5IOrY0aMb\necyd87oxMADQ6vLdhtCA0vzUuHMR6zbim9aONAUcAGQb+LBw4UKuEMi0xgxpiknvWVu8R79Kvht/\n668Dq88Z9/1+5YzeRhgAAPfRoQOg0/vHkMENyXCYnttAF/K+O0/Kh/gbVGY+raC7+VhSoKSVW695\ng8qPmYbTYNW7z0D1KPSisoQjIQefWE9fO9OLUb8vCCQaBYQ81Ql5yMf7gmfBNEc/H3ri3YhrOa6j\n7RWexckRFD17C1ZTxwzyMiICcCuYJKhVNB1em51aDAZ+5i03WY9s6GhBin2RUihwsZP7bIVRTdiD\nprL7Doxe898zl56M9RqgoHK2hbI+oHHZgcyetWYaq6nIEHWMKfDuI2dHM7FnsZruiKAmq/7kiHPf\npBcDxbWz7F0RGtJTEqiGrH5B/ylI+enGzVxeNw8aAGBUAyt7Ayv7Lq7mNQtDr15+5R+sZ0SHurIa\nIQ5UDKRXMUDHVrcpWxMj0cgg5AplFwOj6Dv0HO7Q0z/wzt7lvx4+4s5e4Q0AQLLt7ca5e+/wya6h\nwV9JL89iuuzBnlo7/35YOqh35sMCmuska9kOpYWtp6BBpcdMw//buPtkl6U4/tC6w1lOM9f/0N+s\n6ciRiHhCINEp7eymD83xBTcspu06Mbhvp4LIdat3nY5NSEpLT0l8mFGhbHKKCcsI3sRevJXyIic3\nN7egRu7UipsedTb2fnJqYtyp3b+nSVgBAyxbfaUK0/UY42dcfmX7jtM3H6dnpCfevZNRxS9LvHbj\nYXp2Tk5WRuorDlAYjE8c4kjWIgG/htf0t6GjBUlYkCs0MLdoZMqsP+HCJZ/5W3a8NiMyPFnA6OXv\npC37OM3Rz4cueBJxLadVN6VKe0r/MYO8WNadbRxtFd/k1thTOraxp6x7kdLsFxzSnnL59u0TLMpu\nX3pSpbq0KJodRsBALBTjJANnJ4Yk92EOFwcAEBanpVeRbDxlChBJ17QTiV9e+l7BTEgGrt3MQFRc\nwJEemmTTvnNDFvUj3D+w5Wxm/V2/GN3V35dZeOPm05R7r7Xc+zk0bGvp/lW99Zo3qPqYAWjz7muA\n87Ijth5+ajtjjVz9BQBxVWE16FmqLcVJ433R+4LbFJFNthq+eGbl0ciw7X+LAICkzbTpodd0ekLh\nZ105ElMuIDBsek5Z9d3gttwggNFYE9as0g37X+zJ3ZdFQNJjfbPMk5rzKOpaZoUAANMxdw+cP8Pz\nYy+nKUY2ZFlRriaci3b7hkXmlBGcfT09vh1oGBq762fK2CFuZtriqrfFlK5+3kYEmr4O1D67l5Rn\n0cOa/ikP/feFmRlPmbiAU5qbdjfudrbALnD5FBf5u42kPeWzXZHrVucFDP6KZUYn8oszKuDDyyP1\nKCYsI4iKvXjLpJ+dHplXUCOSLXrc9KizsVWuppTy1MvhaRLW1AGWFPn7B5ST9pTrr2zfgY3zY5to\n8SuraGwfk5e3UsQW1kY6kurMVxyg2Mv1lKKyB9G3Oaa25gwSryQr/kJMOf2rGXZawH997WIa1c7a\nSIfALc66HXmzmuE7wJYKVLuAALv407/90nmqn0XNgzMX3nUaONuTIdOilqW7GUSkF/L6G2iDIC/q\n1+gaG1dHKwNt4BY/u3EuB/QGsA3JUCmdmqTffcbKoNJVp3Yfc9w215tJAKDaDRlicePc72FiqtdK\nB20AkNm/Xiq3noIGdZUcM41Hfxt3Xz1Jxb3wmGKrwGk2vMLcXAAAIGibWJloE8RVr7KqqHYuRui7\nuM/kS/9YjEC37x+0rH9Q82fMR+4IGyn7AKZlPSA4dECwwnZoTuPWzHOhfvAoSb4RCmvWL2GKnyLp\nswPnbw6UfbHt2v3j2rAmbYXpek2fNXBfWOT+7ThGtx26qKdnly6T1q/UO3U6LnxPjBBAy9DZ33Gg\ntxGlU8+J3yT8cu1YZLeuiz0+0a1ABC0mHbJvHNp2AwDIuqa2zv1mrxvmy9JrfgR04J5Swq8qSLkS\nHVnOByAxrNwC5k8e7amLSTicytf34q+e5YgAKPq2HiMXTwxk62AAJHP/Jcv5R49f3L+VT9J3GjRv\n6SRnufGPiAaefTqHRf2dXefmQRODLqPu4fVT1yr4AEDWt3YLXDh5VBcaiCubNoiZ39zZT5ce/P0P\nb8eFPZgEIJl9Pdrzwr4nugOHOGoDABBl9u8ilVtPUYOKj5mP3X31BEXp+TjwL25fc7HxMedFv63x\nppYl3X6j5TbeHt2R9rl0wJHSPrh1FGmfNGY3SSrv/7z017oJO9YONv64jyTCwvMrVyf03rRllAW5\n5anbE96L4/8NTfJZvzPIAZ0FfyboIjuCqETQ7z5xuNmriJMPKpUPEqEILijLzcnPz4w/ue8St1eQ\nn3kHq78geH3laBzPa/I3dqj+fj4dcMh2UrNLFgjyOVE6fzNv+vvHOm38HbKw6MbBtZeKMF3bPtNW\nfuf+ab9a+AJwnGjkGfj9UJ9/eG8molIHvBCBIAiiKVD/hiAIojaoBCMIgqgNKsEIgiBqg0owgiCI\n2qASjCAIojaoBCMIgqgNKsFqhddm3/jz9K23QnUvCIIgaqERJbj14cTtjZiTcS36Roaiwb5ks4Q7\n7goiCKJSu/h1HC4oTr564fLt5BfFdThoGTp27Rfw7TfdzRoytJTE1jb+avKjw4k/JUld3sOrV28+\nTHv5VuSxfN8i96aBTSSc1LCQ7dererVlwATZLOF2sYIIgnx66i/BeN2LiI0bL+Yx2YNHzHYyofKK\nsxJiz+/5MXH4T2vGO2ljymNr2w+cnx93YMv/Xlv6DgqYOdLc0NCyqdDidS8iNu+Jr2nr541/nCWM\nIEj7p/YSzM8+d+BinnHAT6GTnOsHru3Vr7932JpN0QcivHZMY1EFymJr5cmEE2+YmLNxR/Pk2r0L\ndS8fUprp21KYsXI4L/vszj9r/EN2DLf+IPwWxBX3ftl1y2Dq0r5xm88reT036+yGBfk55QIiw7bn\niKnThjrSMQUjjX2YvtyqZUMQpF1T9xuZlx17u4LisXBEF5mBwzEdl8Dxbjf2x1/LHsdik1uMrQUA\nuXBiE5qJouRaOu+x0kzfloJpVcA5T879VYob3dkx71yFgGLm8vXY78b1MKEAgLDo+v4T+T0WbfI1\nvh+nvAUB1WnUzFGmlPIn0f8L21RG3LZa0bi68unLCIJoAjW/l8WcgkIuGLtY0eRLK6bd2dUE0osK\nOGK2ierY2gZy4cQKk2tJAEoyfS25LQTTqsIvePxSpO3oExjgYaEryrt54vC+jcK1P093IhXd+C2i\n+usVE13oWHnTKgvq+EIcBwCMSCUDANDcAscNdaECgIcjvXjxgcs3C76eaNPswoXsCiIIohk6xOmU\nktha5cOvKkmuLZebSCbT16TFYFrlcF5FGR+Mvft/5W5IALCdMvNVUui9m68nGJUdP1/ec8mPjloY\nNA00y3t2aP6Wx3wAALCauvMH+eXWdfA0hfTstzzcRi7QDUEQzaTmEkxkWJprQfLzAm6Aqa7MiTBe\n9+Z5CWh5WTIavnZTEFvbS0/pVYmWkmulZDN9VQfTqlwHMhGgppovlt7iR2Ka6wG34n1l6vV0biVs\nCv67ccq8zdNfjwpZsGLNEDEAAEYz1SPEyzcmXRisow0tiyDIx1H3WbCWw5A+zAc3/ryc7TKBpd14\nE1rm5bOpQqbfEIdmkVXS2NonxQUcETSmbH0YTtyQXBt64+ZTndzXWu7fOmg3S42sz/T1M9fCpLH2\nBblCg6/tWp3w2zAnqpmDAVx9+pIzysKAAMAvflUG+i76+j0W7GDx6wf5Flfe27f1qnHwuu+6Whrp\nNgXeit7KtSWufPqoEExGWHy4EM1XEEEQTaDuEgxarHHzhmdtig5Zle8/tLeTMZVfknk/Jia1wjJg\n1TiWFoCK2NoGzcOJuxqQmifXAoCSTF+M3FIwrQoUiwFD7a+Gh/0eRR/jScu7dvQ+12bCwM40bRKt\nca5iKpMCRG1Dc2Pd5pdPBMXPktNAV1yRfScy4oVW93kDLcnyHYbCFfyozY0gSLui/jcyptNlQsh2\n+8vnL98+dyiGB0A1sHMfueDHQB8LKgYAuNLY2qYmmocTG5CwZsm1UgozfTFai8G0ypHMhi5dKTx2\n8sqededwqlnX0ctnHOJDHgAABuRJREFUBpi3bsMSaJauTsaPruzfdhEAo1t4Dl8weXTzoBjFK9iq\nOSAI0p5pdHDRB8m1GpPpiyCIplD/WfCnhwvKXhfUEQW5cX9c4vZa0vGSaxEE+bfQxBLc0ZNrEQT5\n19DoCxEIgiDtm0YMVokgCNIxoRKMIAiiNqgEIwiCqA0qwQiCIGqDSjCCIIjaoBKMIAiiNqgEAwCA\npDLx5LZtZ3ME6l4QBEH+Vf4dJVg2jVghCbcgLf1FCU/Z8wiCIJ+DRvw6Dq9N2rts5+Na6V9khpmt\ns1ffIf59nZj1qyebRvzpCcqeXAmPiHuUx8ExHXP3ARNnje4mFzsnl6BMKo1dvehUvnwTBgEbd06y\naRyeWNyKaaA+eToyOv5JdkkdDpi2kT2757ApY3wMibWpv67bc69YOlYmRtW3cPAYMGqcnzOjXaWe\nIgiiESUYJLyqWjAcPC/YW0/C45S+Trtz/Uho3K3A5SvGOulgnzeNGK/Ljjx+s9p92PejzIllyZdO\nR+/a22nPWj8jYsPz8gnKRP1e8zfac+tzNCTVKWF7L3C8POUyQFozDc59eX7ThshcHecBAbNczGmi\n6tLCVy9qJAQCAEhqi4v5xv4LZ3szJNyakpxHMRdPbszlbw8JtNCMPY4gGkKD3pA6ZqwuzoZEAOja\nq79fn/D1m6L2hrnuCHbVxmTHSMP5+bdOHD5/71WVGChM1rClq8bayqVjCApjtq8Jq+y/MmQKuxXj\nS2DaLjN/3kOgEDEAALbxu5S1f6cUCvyMaAAKE5RJDAtbhvS/OCd5//WcTgHrJznryM2p5Wn4OecP\nROYyBy3fON29MW9kUP2/0usp2mYOTiwDAoCrh6cd78WqKykFvBEWaMgMBGlHNPRaMEZ3GTmpG4WT\nEJNVJ/+MsPDK3qP3Jb2Df9oQun7RlCHu8mOfiyoSDm85VeA5Z+XE1tRf6cyI9fUXAOdXVXJB31pf\n2mh9gvL8Ob7Givo6/quo4w/wHkHD7alKZ6VkGt6rq/FlJJfx37rptrCUuIhb8fL+tYdlRAcfFQHU\nCIKogwadBcvDaFZsE0h6l18l6qoj87jkfVkNaDuyPZzt6BjYyb1GXJ108ujBZMtpoXP6GH7ElhFX\nJoYfS6L5Lh1qRQYAkcIE5aYlqXx89q9Km7FjPJRXUWXTiDlvCuvAiG2jqpvIPTY/6FjDH5QuExb1\nN9XY3Y0gHdS/7j1JtR02ip14atviPO+BfkP8+nQxoDRUMe6jQwdAp/ePIYMtPmKEYVFZwpGQg0+s\np6+d6cUgAEgqHihKUG6avvju5aeY+0JfFXVR+TQ4DvXxo00P1SbuWnaYP37z8v6dAADAfOTKeT30\nMLGgtjQ3Meb0mXUb8U1rR1i2KpMUQZAvQkMvRADg3DfpxUAx7cz8oHpRLP1X7N+9crwz/8HRjQuX\n/Z5YKal/hmTbu6tB7b3DJxMqRM3aU01YHH/wp4NpdjPXLx5kQQEAwKtSr6dzK//eFDx58pRJUxcd\nzwNuwubpy87nCwEAQFScEF9AdR/EVnEhQfk0RF0zYwqUZr2VvY0OF77n1NbwG1YHqJ3Mraw6d7Zx\ncPH2C1owgyXOuXm/SNjGNUMQ5HPS0LNgvDYjMjxZwOjn76Td/FmMasIeNJXdd2D0mv+eufRkrNcA\nBgAA2bTv3KXdw9btObCF2SlkYhftVl44xXnZEVsPP7WdEfpDf7OGc0xMT0mCshkJAEBUnvqwiMQa\n76hiJqqmoTn6+dAT70Zcy3Edbd+KK7wSQS0PgERBiXMI0q5oUAl+X5iZ8ZSJCziluWl3425nC+wC\nl09x+bB8icoSb6SILayNdCTVma84QLFnUJomIel3n7EyqHTVqd3HHLfN9Wa25kOCpOJeeEyxVeA0\nG15hbi4AABC0TaxMtA3NlSYo49y81Ldg2q8zrXHeOCf5wKo9GY7zts73YRIUT9ME03adGNz32a7I\ndavzAgZ/xTKjE/nFGRUAMqmmdW+zM5+/k/DfV77LehATk09iBfU00aAdjiAaQDPekQQtJh2ybxza\ndgMAyLqmts79Zq8b5svSa7564qqcR1HXMisEAJiOuXvg/BmedEw2MZ5s5jd39tOlB3//w9txYY9W\nFGFBUXo+DvyL29dcbHzMedFva7x1lL9GVJFbIqZYWnzwUwkcb3maxnVmegVvCXG6cCH27pn9l4UA\nQNE1sfPysqIRAAg6RoaU7Gv7Nl0DACDrmti6BM4dOeIr9H0cgrQvKLgIQRBEbTT26zgEQZD2D5Vg\nBEEQtUElGEEQRG1QCUYQBFEbVIIRBEHUBpVgBEEQtUElGEEQRG1QCUYQBFEbVIIRBEHUBpVgBEEQ\ntUElGEEQRG1QCUYQBFGb/wNwktny3kXjJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image('system_info.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "My laptop has one Intel i7 CPU; according to Intel's documentation, each CPU has 4 cores and 8 threads. Thus, I could probably spwan up to 8 processes without any problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examining MPI - \"Hello World\"\n",
    "================\n",
    "\n",
    "Let's write a simple code to show how to write and execute MPI code. Create a new text file, ``mpi_hello.py``, and put the following code in it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "```python\n",
    "from mpi4py import MPI\n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "\n",
    "rank = comm.Get_rank()\n",
    "\n",
    "print('Hello from {0}!'.format(rank))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the code above is entirely valid Python, and Jupyter will happily execute it, there are some subtleties in how to actually exeucte parallelized code from within a Jupyter notebook. (Briefly, even if the code looks like its executing, at the hardware level, we might not actually be taking advantage of all the processes!)\n",
    "\n",
    "Once you've saved ``mpi_hello.py``, you can execute it from the command line as follows:\n",
    "\n",
    "```\n",
    "mpiexec -n [number of processes] python mpi_hello.py\n",
    "```\n",
    "\n",
    "You should get an output that looks similar to\n",
    "\n",
    "```\n",
    "Hello from 4!\n",
    "Hello from 9!\n",
    "Hello from 3!\n",
    "Hello from 7!\n",
    "Hello from 8!\n",
    "Hello from 2!\n",
    "Hello from 1!\n",
    "Hello from 6!\n",
    "Hello from 0!\n",
    "Hello from 5!\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(I called ``mpiexec`` with 10 processes.)\n",
    "\n",
    "**What happened?**\n",
    "\n",
    "What did the command line do? The function call ``mpiexec -n [number of processes]`` tells the MPI software on our machine to spawn new processes. The end of the function call ``python mpi_hello.py`` tells MPI to send to each process the instruction \"call ``python mpi_hello.py``\".\n",
    "\n",
    "Once each process makes its Python call, then a new Python interpreter is spun up for each process, and the code inside ``mpi_hello.py`` is executed. Let's take a look at that.\n",
    "\n",
    "The first line\n",
    "```python\n",
    "from mpi4py import MPI\n",
    "```\n",
    "\n",
    "simply tells Python to go get the MPI module from the ``mpi4py`` package. This module exposes the functions we need in order to be able to communicate between processes.\n",
    "\n",
    "The second line\n",
    "```python\n",
    "comm = MPI.COMM_WORLD\n",
    "```\n",
    "\n",
    "instantiates a communicator which allows communication with all the processes we created when we invoked ``mpiexec``.\n",
    "\n",
    "The third line\n",
    "```python\n",
    "rank = comm.Get_rank()\n",
    "```\n",
    "\n",
    "Uses the ``Get_rank()`` method of the communicator to get the rank assigned MPI assigned to the process. (Recall the rank is a unique identifer within the communicator.)\n",
    "\n",
    "Finally, we printed out the rank using \n",
    "```python\n",
    "print('Hello from {0}!'.format(rank))\n",
    "```\n",
    "Now that we've done a hello world, let's try something a bit more complicated. We have an eye towards parallelizing our random walk simulator over values of $p$. Let's figure out how to do some simple logic on each process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examining MPI - \"Hello Randomness\"\n",
    "===============\n",
    "\n",
    "Using randomness in a parallelized setting can be tricky. Recall that ``numpy`` provides access to several random number generators, as well as a method ``numpy.random.seed()`` to fix the seed of the random number generator. Fixing the seed allows us to have _reproducible randomness_, meaning that, every time we call our code with a fixed value of the seed, we get back the same (random) results.\n",
    "\n",
    "\n",
    "One way we could parallelize our code over $p$ goes as follows:\n",
    "\n",
    "* Spawn many processes.\n",
    "* Within each process, fix the random number seed.\n",
    "* Call ``numpy.random.uniform(0, 1)`` to generate a single random value for $p$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Randomness by Setting a Fixed Seed Across All Processes\n",
    "---------------\n",
    "\n",
    "Create a new file, ``mpi_hello_random-1.py``, and put the following code in it.\n",
    "\n",
    "```python\n",
    "from mpi4py import MPI\n",
    "import numpy as np\n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "\n",
    "rank = comm.Get_rank()\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "p = np.random.uniform(0, 1)\n",
    "\n",
    "print('Hello from {0}! I have value {1}'.format(rank, p))\n",
    "```\n",
    "\n",
    "From the command line, execute this file using ``mpiexec``:\n",
    "\n",
    "```\n",
    "mpiexec -n [number of processes] python mpi_hello_random-1.py\n",
    "```\n",
    "\n",
    "You should obtain an output which is similar to this\n",
    "\n",
    "```\n",
    "Hello from 2. My random value is 0.5488135039273248\n",
    "Hello from 0. My random value is 0.5488135039273248\n",
    "Hello from 1. My random value is 0.5488135039273248\n",
    "```\n",
    "\n",
    "though the exact random value may be different.\n",
    "\n",
    "At first glance, we might think \"_What's going on here_?\". After all, we don't want each process to use the _same_ value for $p$ in the random walk simulator - that's an unnecessary duplication of work!\n",
    "\n",
    "The mistake we made was in setting the _same random number seed for each process_ - numpy doesn't care what the rank of the process is, it simply set the seed to 0, and happily generated a single value for $p$.\n",
    "\n",
    "Is there some way we could get reproducible randomness, without each process returning the same value for $p$? **Yes - we can use a _process-dependent_ seed**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Randomness by Setting a Process-Dependent Seed\n",
    "-------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a new file, ``mpi_hello_random-2.py``, which has the same code as ``mpi_hello_random-1.py``, but where the value of the random number seed is the processes' rank:\n",
    "\n",
    "```python\n",
    "from mpi4py import MPI\n",
    "import numpy as np\n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "\n",
    "rank = comm.Get_rank()\n",
    "\n",
    "#The seed is no longer fixed to be zero!\n",
    "np.random.seed(rank)\n",
    "\n",
    "p = np.random.uniform(0, 1)\n",
    "\n",
    "print('Hello from {0}! I have value {1}'.format(rank, p))\n",
    "```\n",
    "\n",
    "From the command line, execute this file using ``mpiexec``:\n",
    "\n",
    "```\n",
    "mpiexec -n [number of processes] python mpi_hello_random-2.py\n",
    "```\n",
    "\n",
    "You should obtain an output which is similar to this\n",
    "\n",
    "```\n",
    "Hello from 0. My random value is 0.5488135039273248\n",
    "Hello from 1. My random value is 0.417022004702574\n",
    "Hello from 2. My random value is 0.43599490214200376\n",
    "```\n",
    "\n",
    "though your exact random values my be different. By keying the seed to the rank of the process, we're guaranteed that no two processes have the exact same seed (the rank is a _unique_ identifier), which ensures that no two processes will have the same value for $p$ (with high probability - it could be the case that even with two different values of the seed, ``numpy.random.uniform(0, 1)`` returns the same value).\n",
    "\n",
    "A Pitfall to This Approach\n",
    "------------\n",
    "\n",
    "This approach works really well if we are content with the values of $p$ we use in our simulation are random. However, suppose we are told \"Simulate the mean displacement of the random walk with $p=[0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1]$.\" It's incredibly improbable that spawning 11 process, and setting $p$ randomly, would give us those exact values. How can we parallelize over a given _list_ of values?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examining MPI - The \"Master/Slave\" Computing Paradigm\n",
    "=====================\n",
    "\n",
    "Our previous approaches to parallelization used MPI's communicator simply as a way to get the rank of the process. However, the communicator (as the name implies) allows us to communicate with other processes! Let's take advantage of that fact.\n",
    "\n",
    "Suppose we use one process to \"orchestrate\" the parallelization. In particular, given a list of values for $p$, that process should \"chunk\" the list into smaller lists, one for each process, and then pass a chunk to each of the remaining processes. In this way, we can fix the values of $p$ we want to simulate, but also gain a speedup using parallelization - instead of having a single process simulate every value in the list, we have many processes simulate values in a subset of that list.\n",
    "\n",
    "This approach to parallelization is often called the \"master/slave\" approach. As the name implies, one process is the \"master\" (or \"orchestrator\"), while the others are the \"slaves\" (or \"workers\"). For embarassingly parallel problems, this approach makes a lot of sense.\n",
    "\n",
    "To implement this approach, we need to use some _conditional logic_ in our code. We will designate the process with rank 0 as the \"master\" process, and use ``if/then`` statements to break up a designated list of values for $p$.\n",
    "\n",
    "Using the communicator's ``scatter`` method, we'll be able to assign various chunks to different processes, thereby parallelizing our work over $p$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Editing ``simulator.py``\n",
    "------------\n",
    "\n",
    "In your favorite text editor, open up the file ``simulator.py``. This file contains the code necessary to simulate the random walks. (If you look closely, you'll see it incorporates the ``covert_value`` and ``multiple_trajectories`` functions we wrote in the last lesson!)\n",
    "\n",
    "This file uses a little bit of [Python namespace knowledge](https://stackoverflow.com/questions/419163/what-does-if-name-main-do) to make it possible for us to both call this module from the command line (``python simulator.py``), as well as import it in other Python modules (``from simulator import X``, or ``import simulator``).\n",
    "\n",
    "For the purposes of this lesson, the two functions we are interested in are ``chunks`` and ``main``. ``chunks`` will be used to take the list of $p$ values and chunk it up into a number of pieces equal to the number of processes we spawned in our ``mpiexec`` call. It uses ideas based on [this StackOverflow post](https://stackoverflow.com/questions/312443/how-do-you-split-a-list-into-evenly-sized-chunks) to do the chunking.\n",
    "\n",
    "**Writing ``main()``**\n",
    "\n",
    "Navigate to the function ``main()``. This function will be executed by default if when we call ``simulator`` from the command line, and as such, it's where we need to implement parallelization.\n",
    "\n",
    "**Exposing the MPI Communicator**\n",
    "\n",
    "First, we need to instantiate an instance of MPI's communicator, using ``MPI.COMM_WORLD``. Then, we have access to its attributes ``Get_rank()`` and ``Get_size()``, which tell us the rank of the process, and the total number of processes, respectively.\n",
    "\n",
    "> Write 3 lines to create an MPI communicator, and assign the output of ``Get_X()`` to a variable called X.\n",
    "\n",
    "**Executing Conditional Logic Based on Processor Rank**\n",
    "\n",
    "We don't need every process to go and try to chunk up our list of values for $p$. Let's condition it on the rank of the process being zero. (That way, if we call the code from the command line using ``python simulator.py``, we will actually execute the simulator, albeit on a single process.)\n",
    "\n",
    "Conditioned on the rank of the process being zero, we need to do the following:\n",
    "\n",
    "* Check to see whether the directory ``simulated_data`` exists; if not, make it.\n",
    "* Generate a list of values for $p$.\n",
    "* Shuffle those values (for _load balancing_ purposes)\n",
    "* Chunk the data using the ``chunks`` function\n",
    "\n",
    "If the rank is _not_ zero, assign its values for $p$ to be ``None``.\n",
    "\n",
    "> Write code which implements the logic above.\n",
    "\n",
    "**Scatter the Chunks to Other Processors**\n",
    "\n",
    "Once we have the chunks, we need to send them to the other processes. We can do that using ``comm.scatter()``, where we set the _root_ of the scatter method to be the rank-0 process. (That way, what we scatter is the chunks, not the ``None`` type object.) By assigning the output of that scatter to some local variable, each process will get a different chunk.\n",
    "\n",
    "(You can check this by putting a ``print()`` statement after the scatter call.)\n",
    "\n",
    "> _Outside of the conditional logic statement_, assign the output ``comm.scatter()``, with ``root=0``, to a variable, where ``comm.scatter()`` takes as input your list of values for $p$.\n",
    "\n",
    "It is true the _number_ of values each process will have is not necessarily the same (particularly if the number of values for $p$ does not evenly divide the number of processes).\n",
    "\n",
    "**Generate Data**\n",
    "\n",
    "After the scatter command, we can now call ``data_gen`` on the local variable to generate the trajectories of the random walk.\n",
    "\n",
    "> Add a line to call ``data_gen()`` with the variable you assigned to ``comm.scatter()``."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Executing the Simulator\n",
    "==============\n",
    "\n",
    "Once you've finished writing ``simulator.py``, go ahead and call it from the command line using ``mpiexec``\n",
    "\n",
    "```\n",
    "mpiexec -n [number of processes] python simulator.py\n",
    "```\n",
    "\n",
    "Once that finishes, you can run ``python plotter.py`` to make a plot of the mean displacement versus time.\n",
    "\n",
    "**_Congratulations! You've successfully parallelized the computation of trajectories of a random walk!_**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extensions/Ideas\n",
    "===========\n",
    "\n",
    "* To see how this approach gives a speedup, go back to ``simulator`` and change the default options for ``multiple_trajectories`` to increase the total time $T$ and number of trajectories. (Set them to some really-high values, such as $10^{4}$ and $10^{3}$, respectively.) Now, compare the runtime of these two command line calls:\n",
    "\n",
    "``python simulator.py``\n",
    "\n",
    "and \n",
    "\n",
    "``mpiexec -n [number of processes] python simulator.py``\n",
    "\n",
    "The first runs the code in _serial_ mode, meaning it iterates over each value of $p$, while the second does the chunking described above. As the total time for simulating _each_ set of multiple trajectories increases, you'll gain an advantage when iterating over values of $p$ by using MPI.\n",
    "\n",
    "* Look up Python's ``sys.argv`` module, and modify ``simulator.py`` to accept command-line inputs. (For example, how many values of $p$ you want to use.)\n",
    "\n",
    "* Write a bash script to combine calling ``simulator`` and ``plotter``.\n",
    "\n",
    "* To get a sense of what's going on \"under the hood\", call `mpiexec` with the `-verbose` flag. Turns out there a lot which is hidden!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}