Python/Parallel Programming.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parallel Programming in Python\n",
    "\n",
    "Before we get into it, let's first discuss two common parallel programming paradigms: \n",
    "\n",
    "* shared memory parallel programming (sometimes referred to as threading or vectorization)\n",
    "* distributed memory parallel programming\n",
    "\n",
    "The first of these, is a paradigm where there are a number of threads being executed in parallel, but all the threads have access to the same memory. For example, openMP is a shared memory programming paradigm. Shared memory is both a blessing and a curse. While sharing memory doesn't have very much overhead in terms of communication, it also means that the programmer must be very careful when changing addresses in memory, since depending on the order in which different threads access memory and execute their programs, the overall program may yield different results (this is called a race condition). \n",
    "\n",
    "Distributed memory parallel programming avoids this issue by doing exactly what its name says, distributing all data in memory to all processes. In practice, this means that when we run a distributed memory program, each process runs its own separate version of the code, and the programmer has to tell the processes exactly how memory needs to be distributed.\n",
    "\n",
    "We will address both of these paradigms in this lesson, focusing more on distributed memory programming for reasons that will become clear in a moment. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Threading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "import urllib3\n",
    "import shutil\n",
    "import re\n",
    "import os\n",
    "\n",
    "urllib3.disable_warnings()\n",
    "\n",
    " \n",
    "def download_xkcd(start, stop):\n",
    "    \n",
    "    http = urllib3.PoolManager()\n",
    "    \n",
    "    path = \"xkcd_\"\n",
    "    content = True\n",
    "    for i in range(start, stop):\n",
    "            url = \"http://www.xkcd.com/\"+str(i)+\"/\"\n",
    "            rd = http.request('GET', url)\n",
    "            data = rd.data\n",
    "            res = re.search(b\"/comics/[a-z0-9_()]*.(jpg|png)\", data)\n",
    "            if res:\n",
    "                imgurl = \"http://imgs.xkcd.com\"+res.group()\n",
    "                with http.request('GET', imgurl, preload_content=False) as r, open(path+str(i)+imgurl[-4:], 'wb') as out_file:       \n",
    "                    shutil.copyfileobj(r, out_file)                \n",
    "            else:\n",
    "                if re.search(b\"Not Found\", data) and i != 404:\n",
    "                        content = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "195 ms ± 6.23 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "download_xkcd(10,20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from queue import Queue\n",
    "from threading import Thread\n",
    "import threading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def download_worker(q):\n",
    "    \n",
    "    http = urllib3.PoolManager()\n",
    "\n",
    "    path = \"xkcd_\"\n",
    "    \n",
    "    while not q.empty():\n",
    "        i = q.get()\n",
    "        q.task_done()\n",
    "        url = \"http://www.xkcd.com/\"+str(i)+\"/\"\n",
    "        rd = http.request('GET', url)\n",
    "        data = rd.read()\n",
    "        res = re.search(b'/comics/[a-z0-9_()]*.(jpg|png)', data)\n",
    "        if res:\n",
    "            imgurl = \"http://imgs.xkcd.com\"+res.group()\n",
    "            with http.request('GET', imgurl, preload_content=False) as r, open(path+str(i)+imgurl[-4:], 'wb') as out_file:       \n",
    "                shutil.copyfileobj(r, out_file)               \n",
    "        else:\n",
    "            if re.search(b\"Not Found\", data) and i != 404:\n",
    "                print('Failed to download comic {}'.format(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DownloadThread(threading.Thread):\n",
    "    def __init__(self, function, args):\n",
    "        self.running = False\n",
    "        self.function = function\n",
    "        self.args = args\n",
    "        super(DownloadThread, self).__init__()\n",
    "\n",
    "    def start(self):\n",
    "        self.running = True\n",
    "        super(DownloadThread, self).start()\n",
    "\n",
    "    def run(self):\n",
    "        while self.running:\n",
    "            self.function(*self.args)\n",
    "\n",
    "    def stop(self):\n",
    "        self.running = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def download_xkcd_threaded(nthreads, start, stop):\n",
    "    threads = []\n",
    "\n",
    "    q = Queue(maxsize=0)\n",
    "    \n",
    "    for i in range(start, stop):\n",
    "        q.put(i)\n",
    "        \n",
    "    for i in range(nthreads):\n",
    "        t = DownloadThread(download_worker, (q,))\n",
    "        t.start()\n",
    "        threads.append(t)  \n",
    "        \n",
    "    q.join()\n",
    "    \n",
    "    while len(threads)>0:\n",
    "        del threads[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Automatic pdb calling has been turned OFF\n"
     ]
    }
   ],
   "source": [
    "%pdb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "176 ms ± 29.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "download_xkcd_threaded(5, 10, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "345 ms ± 53.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "download_xkcd_threaded(10, 10, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "774 ms ± 189 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "download_xkcd_threaded(20, 10, 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1\n",
    "Why does increasing the number of threads more not yield faster results?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The GIL\n",
    "Unfortunately, we can't use threading to do all of our parallel programming in python. The reason for this is something called the Global Interpreter Lock (GIL). The GIL is the bane of threading in python. It is a mutex (mutual exclusion object) that protects access to python objects, preventing multiple threads from executing python bytecodes at once. Note that the GIL only exists in Cython (the most common implementation of the python language in C). If many threads are running then the GIL is passed between the different threads. Sometimes the GIL is passed voluntarily, but if one thread holds the GIL for a long time (1000 bytecodes in python 2) then it will have the GIL forcibly revoked.\n",
    "\n",
    "**Note that the GIL does not prevent all race conditions**\n",
    "\n",
    "\n",
    "Let's take a quick look at why:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def threaded_foo():\n",
    "\n",
    "    def foo():\n",
    "        global n\n",
    "        n += 1\n",
    "                \n",
    "    threads = []\n",
    "    for i in range(100):\n",
    "        t = threading.Thread(target=foo)\n",
    "        threads.append(t)\n",
    "\n",
    "    for t in threads:\n",
    "        t.start()\n",
    "\n",
    "    for t in threads:\n",
    "        t.join()\n",
    "\n",
    "    print(n)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100\n",
      "200\n",
      "300\n",
      "400\n",
      "500\n",
      "600\n",
      "700\n",
      "800\n",
      "900\n",
      "1000\n",
      "1100\n",
      "1200\n",
      "1300\n",
      "1400\n",
      "1500\n",
      "1600\n",
      "1700\n",
      "1800\n",
      "1900\n",
      "2000\n",
      "2100\n",
      "2200\n",
      "2300\n",
      "2400\n",
      "2500\n",
      "2600\n",
      "2700\n",
      "2800\n",
      "2900\n",
      "3000\n",
      "3100\n",
      "3200\n",
      "3300\n",
      "3400\n",
      "3500\n",
      "3600\n",
      "3700\n",
      "3800\n",
      "3900\n",
      "4000\n",
      "4100\n",
      "4200\n",
      "4300\n",
      "4400\n",
      "4500\n",
      "4600\n",
      "4700\n",
      "4800\n",
      "4900\n",
      "5000\n",
      "5100\n",
      "5200\n",
      "5300\n",
      "5400\n",
      "5500\n",
      "5600\n",
      "5700\n",
      "5800\n",
      "5900\n",
      "6000\n",
      "6100\n",
      "6200\n",
      "6300\n",
      "6400\n",
      "6500\n",
      "6600\n",
      "6700\n",
      "6800\n",
      "6900\n",
      "7000\n",
      "7100\n",
      "7200\n",
      "7300\n",
      "7400\n",
      "7500\n",
      "7600\n",
      "7700\n",
      "7800\n",
      "7900\n",
      "8000\n",
      "8100\n",
      "8200\n",
      "8300\n",
      "8400\n",
      "8500\n",
      "8600\n",
      "8700\n",
      "8800\n",
      "8900\n",
      "9000\n",
      "9100\n",
      "9200\n",
      "9300\n",
      "9400\n",
      "9500\n",
      "9600\n",
      "9700\n",
      "9800\n",
      "9900\n",
      "10000\n"
     ]
    }
   ],
   "source": [
    "n=0\n",
    "for i in range(100):\n",
    "    threaded_foo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "Why doesn't n increment by 100 every time?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Practically, what this means is that only latency bound programs, programs that spend most of their waiting for other things to happen such as reading or writing a file, will experience speed ups from threading.\n",
    "\n",
    "The way to get around the issue of the GIL is to use multiple python interpreters in order to parallelize code. One common implementaion of such parallelism is the multiprocessing package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most basic way to use the multiprocessing module is to use the ``Process`` class. The syntax is very similar to ``Threading``, except now each process is given its own python interpreter, and thus each one has its own GIL, and doesn't need to pass it around."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['ZVby6', 'mgGYI', '1eUxg', 'dxdX4']\n"
     ]
    }
   ],
   "source": [
    "import multiprocessing as mp\n",
    "import random\n",
    "import string\n",
    "\n",
    "random.seed(123)\n",
    "\n",
    "# Define an output queue\n",
    "output = mp.Queue()\n",
    "\n",
    "# define a example function\n",
    "def rand_string(length, output):\n",
    "    \"\"\" Generates a random string of numbers, lower- and uppercase chars. \"\"\"\n",
    "    rand_str = ''.join(random.choice(\n",
    "                        string.ascii_lowercase\n",
    "                        + string.ascii_uppercase\n",
    "                        + string.digits)\n",
    "                   for i in range(length))\n",
    "    output.put(rand_str)\n",
    "\n",
    "# Setup a list of processes that we want to run\n",
    "processes = [mp.Process(target=rand_string, args=(5, output)) for x in range(4)]\n",
    "\n",
    "# Run processes\n",
    "for p in processes:\n",
    "    p.start()\n",
    "\n",
    "# Exit the completed processes\n",
    "for p in processes:\n",
    "    p.join()\n",
    "\n",
    "# Get process results from the output queue\n",
    "results = [output.get() for p in processes]\n",
    "\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The ``Pool`` class\n",
    "A more convenient approach to simple parallel processing is using ``Pool``. Let's set up a simple Monte Carlo integration. We'll go over ``apply``, and ``apply_async``, but there are also ``map`` and ``map_asyn`` functions that are similar to they builtin python map functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def counts_in_a_circle(n):\n",
    "    np.random.seed()\n",
    "    x = np.random.uniform(-1, 1, size=n)\n",
    "    y = np.random.uniform(-1, 1, size=n)\n",
    "    \n",
    "    n_in_circle = np.sum(np.sqrt(x**2 + y**2)<1)\n",
    "    \n",
    "    return n_in_circle, n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def monte_carlo_area_mp(n_proc, n_samples):\n",
    "\n",
    "    pool = mp.Pool(processes=n_proc)\n",
    "    results = np.array([pool.apply(counts_in_a_circle, args=(n_samples//n_proc,)) for x in range(n_proc)])    \n",
    "    \n",
    "    area = 4*np.sum(results[:,0])/np.sum(results[:,1])\n",
    "    \n",
    "    print('Monte carlo area, fractional error: {}, {}'.format(area, np.abs(area-np.pi)/np.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Monte carlo area, fractional error: 3.14186344, 8.619399141301981e-05\n",
      "Monte carlo area, fractional error: 3.1415376, 1.7524101901078094e-05\n",
      "Monte carlo area, fractional error: 3.14180432, 6.737551094184428e-05\n",
      "Monte carlo area, fractional error: 3.14147244, 3.8265174084851955e-05\n",
      "Monte carlo area, fractional error: 3.14145104, 4.50770056490629e-05\n",
      "Monte carlo area, fractional error: 3.14133904, 8.072771290168311e-05\n",
      "Monte carlo area, fractional error: 3.14142656, 5.286923166283999e-05\n",
      "Monte carlo area, fractional error: 3.14111948, 0.00015061583151223393\n",
      "7.71 s ± 1.18 s per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "monte_carlo_area_mp(2, 100000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Monte carlo area, fractional error: 3.14158552, 2.2706921551221597e-06\n",
      "Monte carlo area, fractional error: 3.1414928, 3.178438480209791e-05\n",
      "Monte carlo area, fractional error: 3.1418304, 7.567703277353008e-05\n",
      "Monte carlo area, fractional error: 3.1415322, 1.924297528648465e-05\n",
      "Monte carlo area, fractional error: 3.14151052, 2.6143933618850822e-05\n",
      "Monte carlo area, fractional error: 3.14190312, 9.88245276967883e-05\n",
      "Monte carlo area, fractional error: 3.1415612, 1.0011988587136655e-05\n",
      "Monte carlo area, fractional error: 3.14154516, 1.511767916145237e-05\n",
      "7.4 s ± 708 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "monte_carlo_area_mp(4, 100000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def monte_carlo_area_async(n_proc, n_samples):\n",
    "\n",
    "    pool = mp.Pool(processes=n_proc)\n",
    "    results = [pool.apply_async(counts_in_a_circle, args=(n_samples//n_proc,)) for x in range(n_proc)]\n",
    "    print('Function called')\n",
    "\n",
    "    output  = np.array([p.get() for p in results])\n",
    "    print('Results returned')\n",
    "    \n",
    "    area = 4*np.sum(output[:,0])/np.sum(output[:,1])\n",
    "    \n",
    "    print('Monte carlo area, fractional error: {}, {}'.format(area, np.abs(area-np.pi)/np.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Function called\n",
      "Results returned\n",
      "Monte carlo area, fractional error: 3.14151568, 2.450145460622952e-05\n"
     ]
    }
   ],
   "source": [
    "monte_carlo_area_async(8, 100000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPI\n",
    "MPI (Message Passing Interface) is probably the most useful parallel and widely used programming paradigm for scientific computing. It is similar to multiprocessing, but is much more explicit about communication, is more flexible and can scale to many machines.\n",
    "\n",
    "The idea behind MPI is very simple. You write your program as if each process is running it independently. Every process will allocate different memory, and run all computations independently. The programming must explicitly outline exceptions to this. We will see how this works below using the same monte carlo program as above. \n",
    "\n",
    "Th\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I will first demonstrate how to use the simplest functionality in MPI, sending and receiving buffers between individual tasks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing monte_carlo_area_sendrcv.py\n"
     ]
    }
   ],
   "source": [
    "%%file monte_carlo_area_sendrcv.py\n",
    "\n",
    "from __future__ import print_function, division\n",
    "from mpi4py import MPI\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "def monte_carlo_area_sendrcv(n_samples):\n",
    "    \n",
    "    comm = MPI.COMM_WORLD\n",
    "    \n",
    "    rank = comm.Get_rank()\n",
    "    size = comm.Get_size()\n",
    "    \n",
    "    print('I am rank {} of {}'.format(rank, size))\n",
    "    \n",
    "    n_samples_per_task = n_samples // size\n",
    "    \n",
    "    x = np.random.uniform(-1, 1, size=n_samples_per_task)\n",
    "    y = np.random.uniform(-1, 1, size=n_samples_per_task)\n",
    "    \n",
    "    n_in_circle = np.sum(np.sqrt(x**2 + y**2)<1)\n",
    "\n",
    "    n_in_circle_sum = np.array([0])\n",
    "    draws_sum = n_samples_per_task * size\n",
    "    \n",
    "    if rank==0:\n",
    "        n_in_circle_sum = n_in_circle\n",
    "        \n",
    "        for i in range(size-1):\n",
    "            temp = comm.recv()\n",
    "            print(temp)\n",
    "            n_in_circle_sum += temp[0]\n",
    "            \n",
    "    else:\n",
    "        comm.send([n_in_circle], 0)\n",
    "\n",
    "    area = 4 * n_in_circle_sum/draws_sum\n",
    "    \n",
    "    print('[rank {}]: Monte carlo area, fractional error: {}, {}'.format(rank, area, np.abs(area-np.pi)/np.pi))\n",
    "\n",
    "if __name__=='__main__':\n",
    "    \n",
    "    n_samples = int(sys.argv[1])\n",
    "    \n",
    "    monte_carlo_area_sendrcv(n_samples)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first few lines in this code are generic to almost any MPI code. If you have used MPI in C or some other language before, you'll notice that I don't have to call MPI_Init. This is done automatically in the MPI import statement.\n",
    "\n",
    "Note the Barrier command. This forces all tasks to wait at that command until every one of them has reached that point in the program."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another commonly used set of functions is scatter and gather. These are inverses in some sense, with scatter dividing up an array on the root process and distributing those pieces to the other processes. Gather then takes many arrays on different processes and combines them in the root process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing monte_carlo_area_scattergather.py\n"
     ]
    }
   ],
   "source": [
    "%%file monte_carlo_area_scattergather.py\n",
    "\n",
    "from __future__ import print_function, division\n",
    "from mpi4py import MPI\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "def monte_carlo_area_scattergather(n_samples):\n",
    "    \n",
    "    comm = MPI.COMM_WORLD\n",
    "    \n",
    "    rank = comm.Get_rank()\n",
    "    size = comm.Get_size()\n",
    "    \n",
    "    print('I am rank {} of {}'.format(rank, size))\n",
    "    \n",
    "    comm.Barrier()\n",
    "    \n",
    "    if rank == 0:\n",
    "        n_samples_per_task = [n_samples // size for i in range(size)]\n",
    "    else:\n",
    "        n_samples_per_task = None\n",
    "    \n",
    "    n_samples_per_task = comm.scatter(n_samples_per_task, root=0)\n",
    "    \n",
    "    x = np.random.uniform(-1, 1, size=n_samples_per_task)\n",
    "    y = np.random.uniform(-1, 1, size=n_samples_per_task)\n",
    "    \n",
    "    n_in_circle = np.sum(np.sqrt(x**2 + y**2)<1)\n",
    "\n",
    "    n_in_circle_sum = np.array([0])\n",
    "    draws_sum = n_samples_per_task * size\n",
    "\n",
    "    results = comm.gather(n_in_circle, root=0)\n",
    "    \n",
    "    if rank==0:\n",
    "        n_in_circle_sum = np.sum(results)\n",
    "\n",
    "    area = 4 * n_in_circle_sum/draws_sum\n",
    "    \n",
    "    print('[rank {}]: Monte carlo area, fractional error: {}, {}'.format(rank, area, np.abs(area-np.pi)/np.pi))\n",
    "\n",
    "if __name__=='__main__':\n",
    "    \n",
    "    n_samples = int(sys.argv[1])\n",
    "    \n",
    "    monte_carlo_area_scattergather(n_samples)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll notice that in the last example I gathered everything and then had to sum the elements manually. This is a very common operation that you will have to perform, and as such there is a command which does this for you, ``Reduce``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing monte_carlo_area_reduce.py\n"
     ]
    }
   ],
   "source": [
    "%%file monte_carlo_area_reduce.py\n",
    "\n",
    "from __future__ import print_function, division\n",
    "from mpi4py import MPI\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "def monte_carlo_area_reduce(n_samples):\n",
    "    \n",
    "    comm = MPI.COMM_WORLD\n",
    "    \n",
    "    rank = comm.Get_rank()\n",
    "    size = comm.Get_size()\n",
    "    \n",
    "    print('I am rank {} of {}'.format(rank, size))\n",
    "    \n",
    "    comm.Barrier()\n",
    "    \n",
    "    n_samples_per_task = n_samples // size\n",
    "    \n",
    "    x = np.random.uniform(-1, 1, size=n_samples_per_task)\n",
    "    y = np.random.uniform(-1, 1, size=n_samples_per_task)\n",
    "    \n",
    "    n_in_circle = np.sum(np.sqrt(x**2 + y**2)<1)\n",
    "\n",
    "    n_in_circle_sum = np.array([0])\n",
    "    draws_sum = n_samples_per_task * size\n",
    "    \n",
    "    n_in_circle_sum = comm.reduce(n_in_circle, root=0, op=MPI.SUM)\n",
    "    \n",
    "    if n_in_circle_sum is None:\n",
    "        n_in_circle_sum = 0\n",
    "        \n",
    "        \n",
    "    area = 4 * n_in_circle_sum/draws_sum\n",
    "    \n",
    "    print('[rank {}]: Monte carlo area, fractional error: {}, {}'.format(rank, area, np.abs(area-np.pi)/np.pi))\n",
    "\n",
    "if __name__=='__main__':\n",
    "    \n",
    "    n_samples = int(sys.argv[1])\n",
    "    \n",
    "    monte_carlo_area_reduce(n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "Look into the bcast function. Write a new monte carlo function which makes use of it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many variations of these functions that will be useful in various applications. Some common variations are \n",
    "\n",
    "* all{gather, reduce}\n",
    "* i{send,recv,...}\n",
    "* Gather, Reduce, ...."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}