6 - Algorithm Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ch. 6. Algorithm Analysis "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I. Intro to Parallel Algorithm Analysis\n",
    "This course assumes you have done some computer science, but does not assume any knowledge of algorithm design or formal algorithm analysis. However, I still want to do some basic analysis of some parallel algorithms and serial algorithms, so you have some understanding of what kinds of things will make your code faster, and what kind of code can even be made faster. We're going to be very informal about our algorithm analysis because it's more fun that way, and we will still get a few important points across. After this chapter, I want you to have some intuition that helps you understand why this graph looks the way it does:\n",
    "![parallel time graph](https://ars.els-cdn.com/content/image/1-s2.0-S0098300413001465-gr9.jpg)\n",
    "So, let's get started. We talked, as far back as Ch. 3, about what a parallel algorithm looks like and what makes an algorithm parallel. We talked about trivially parallelizable algorithms and monte carlo simulations. What we didn't talk about, specifically, is how these different algorithms change how workflows run. Essentially, if something is trivially parallelizable, if you go from 1 core to _n_ cores, without changing any parameters, you should see a roughly _n_ times speedup. This means that the run time of this program should be roughly $\\frac{1}{n}\\times (original\\ running\\ time)$. This means, as you add more cores to your trivially parallelizable code, you should see your code's runtime shrink as $\\frac{1}{n}$. This graph is what that looks like:\n",
    "![1/n](https://i.imgur.com/MbqqIAn.png)\n",
    "\n",
    "This graph assumes that your code initally took 10 seconds to run, and by the time you add five cores, you already only take ten seconds to run your code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 6.1 - Monte Carlo Analysis\n",
    "Given what we just talked about, Monte Carlo simulations should fit nicely into that class of algorithms. We should see a linear speedup and a resulting runtime proportional to $\\frac{1}{n}$. In this example, we're going to investigate that and try to show with empirical data that this holds for this case. This should hopefully complement our theoretical analysis from earlier and help convince you that this makes sense. We're going to use our parallel monte carlo pi from earlier, and time it as we scale it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining parallel monte carlo pi calculation \n",
    "\n",
    "import random\n",
    "import multiprocessing\n",
    "from multiprocessing import Pool\n",
    "\n",
    "\n",
    "#caculate the number of points in the unit circle\n",
    "#out of n points\n",
    "def monte_carlo_pi_part(n):\n",
    "    \n",
    "    count = 0\n",
    "    for i in range(n):\n",
    "        x=random.random()\n",
    "        y=random.random()\n",
    "        \n",
    "        # if it is within the unit circle\n",
    "        if x*x + y*y <= 1:\n",
    "            count=count+1\n",
    "        \n",
    "    #return\n",
    "    return count  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "# Note that this cell takes a while to run\n",
    "\n",
    "n = 100000000\n",
    "\n",
    "times = []\n",
    "# Call the process with different sized pools\n",
    "for np in range(1,32, 2):\n",
    "    start = time.time()\n",
    "    part_count=[n//np for i in range(np)]\n",
    "    pool = Pool(processes=np)   \n",
    "    # parallel map\n",
    "    count=pool.map(monte_carlo_pi_part, part_count)\n",
    "    pi = sum(count)/(n*1.0)*4\n",
    "    end = time.time()\n",
    "    tot = end - start\n",
    "    times.append(tot)\n",
    "    print(\"Esitmated value of Pi: {} ({} cores, time: {}s)\".format(pi, np, tot))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# Plot the times we collected\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "#h Helper function to plot an equation\n",
    "def graph(formula, x_range):  \n",
    "    x = np.array(x_range)  \n",
    "    y = formula(x) \n",
    "    plt.plot(x, y)  \n",
    "\n",
    "# Graph theoretical maximum (in blue)\n",
    "graph(lambda x: 32/x, range(1, 32))\n",
    "\n",
    "# Graph empirical data (in orange)\n",
    "plt.plot([i*2 for i in range(1,17)], times)\n",
    "\n",
    "plt.ylabel(\"Execution Time (sec)\")\n",
    "plt.xlabel(\"Number of cores\")\n",
    "plt.title(\"Execution Time vs Number of Cores for Monte Carlo Pi\")\n",
    "plt.show()\n",
    "# Hopefully, this looks a little like a nice 1/x curve. If not, something is probably broken"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II. Strictly Serial Algorithms\n",
    "We've seen how algorithms that respond very well to scale do, now what about the other end of the spectrum? If there are algorithms that can be parallelized really easily, it seems like there _should_ be some that can not. Why? Well, it just feels like it! While that's not really valid logic, the assumption is correct. There's a family of algorithms that can simply not be parallelized. This family is called _strictly serial algorithms_, and they are more common than you might think. One such common example is the Fibonacci sequence. The Fibonacci sequence is defined recursively, such that the _n_-th number in the sequence is the sum of the _n-1_-th and the _n-2_-th. That is to say that _fib(n)_ = _fib(n-1)+fib(n-2)_. Because of this, dependence on the previous values, computing the _n_-th Fibonacci number is a strictly serial task. (This is not actually true, as the Fibonacci sequence does indeed have a closed form solution, but that's not the point of this exercise. There are sequences which only have recursive forms, but they're more complicated). This image shows how Fibonacci number computations are always recursively broken down into calls of `fib(1)` or `fib(o)`:\n",
    "![fib call tree](https://i.stack.imgur.com/7iU1j.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 6.2 - Fibonacci Sequence Generator\n",
    "Because we mentioned how the Fibonacci sequence is a strictly serial algorithm, we're now going to use it as an example. We will write a recursive Fibonacci sequence generator and perform some analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute Nth fibonacci number\n",
    "def fib(n):\n",
    "    if n==0 or n==1:\n",
    "        return 1\n",
    "    return fib(n-1) + fib(n-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "# Feel free to play around with this\n",
    "times = []\n",
    "for i in range(0, 40, 2):\n",
    "    before = time.time()\n",
    "    res = fib(i)\n",
    "    after = time.time()\n",
    "    tot = after-before\n",
    "    print(\"fib({}) = {} ({}sec)\".format(i, res, tot))\n",
    "    times.append(tot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "# Plot the times we collected\n",
    "\n",
    "plt.plot([2*i for i in range(20)], times)\n",
    "plt.ylabel(\"Execution Time (sec)\")\n",
    "plt.xlabel(\"Number in Fibonacci sequence\")\n",
    "plt.title(\"Execution Time vs Number in Fibonacci Sequence\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully, you can see that attempting to parallelize any computation like this is hopeless. \n",
    "\n",
    "Let's examine the call tree for calculating the fifth fibonacci number:\n",
    "```\n",
    "fib(5) -> fib(4)+fib(3)\n",
    "            |      |\n",
    "      fib(3)+fib(2)|\n",
    "                fib(2)+fib(1)\n",
    "and so on\n",
    "```\n",
    "\n",
    "Now, the execution will still be sequential. This is because, even if you could fork a thread for each call in the tree, it would spawn so many new threads that it simply wouldn't be worth it to do that. It would slow the program down so much that it would end up slower than the serial algorithm you started with.  \n",
    "\n",
    "In general, you want to parallelize only \"embarrassingly parallel\" tasks. That is, tasks which are computationally expensive and can be computed independently. Many people forget about the first part. Threads are so expensive that it only makes sense to make one when you have a huge amount of work for them to do, and moreover, that you can devote an entire processor to the thread. If you have 8 processors then making 80 threads is going to force them to share the processor, slowing each one of them down tremendously. You'd do better to make only 8 threads and let each have 100% access to the processor when you have an embarrassingly parallel task to perform.\n",
    "\n",
    "In this case, it is theoretically possible to parallelize the algorithm (slightly), but practically not worth it. In some cases, such as non primitive recursive problems (like [the Ackermann function](https://en.wikipedia.org/wiki/Ackermann_function)), it's completely impossible. \n",
    "\n",
    "Also, note that this could be sped up by a lot through the use of [dynamic programming](https://en.wikipedia.org/wiki/Dynamic_programming), but that's not what this course is about."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III. Task Parallelism vs. Data Parallelism\n",
    "As you may gather from the title of this section, there are two flavors of parallelism which are important: task parallelism and data parallelism. Do you remember when we talked about Single Instruction, Multiple Data (SIMD)? Well, if you did (or if you didn't), you might be able to guess that SIMD is an example of data parallelism. Let's define data parallelism better. \n",
    "\n",
    "Data parallelism focuses on distributing the data across different nodes, which operate on the data in parallel. It can be applied on regular data structures like arrays and matrices by working on each element in parallel. It contrasts to task parallelism as another form of parallelism, where tasks are distributed rather than data (we'll get to this in a moment). \n",
    "\n",
    "As this definition suggests, data parallelism occurs when the same thing needs to happen to many different pieces of input data. This is a nice thing to have happen, because it is generally easy to make happen in parallel. Data parallelism is something you should look for when you're trying to speed up your code. It is easy to parallelize, and many different types of problems can be reduced to them.\n",
    "\n",
    "This image represents a data-parallel process:\n",
    "![data parallel job](https://upload.wikimedia.org/wikipedia/commons/a/a7/Sequential_vs._Data_Parallel_job_execution.png)\n",
    "The top describes a job happening in serial, while the bottom shows how one would divide up the top job into the four jobs represented on the bottom, but producing identical output.\n",
    "\n",
    "Now, let's get into task parallelism. What is task parallelism? Well, task parallelism is a form of parallelism where data is not distributed, but instructions are. Using our acronyms from before, it would be Multiple Instruction, Single Data (MISD). Task parallelism focuses on distributing tasks across different processors. In contrast to data parallelism which involves running the same task on different components of data, task parallelism is distinguished by running many different tasks at the same time on the same data.\n",
    "\n",
    "In a multiprocessor system, task parallelism is achieved when each processor executes a different thread (or process) on the same or different data. The threads may execute the same or different code. In the general case, different execution threads communicate with one another as they work, but is not a requirement. Communication usually takes place by passing data from one thread to the next as part of a workflow.\n",
    "\n",
    "As a simple example, if a system is running code on a 2-processor system (CPUs \"a\" & \"b\") in a parallel environment and we wish to do tasks \"A\" and \"B\", it is possible to tell CPU \"a\" to do task \"A\" and CPU \"b\" to do task \"B\" simultaneously, thereby reducing the run time of the execution.\n",
    "\n",
    "The image below represents a comparison between task and data parallelism, with task parallelism on the left and data parallelism is on the right.\n",
    "\n",
    "![parallelism comparison](http://lh3.googleusercontent.com/8hC1WoryRSmlqL7iQ3VtQwRnCLeO2Wa5zejBw3pmUlXxLpoqhRC2OhNPLZTPfz5OlD4tdmt7V6YqXqBqYmSgg_CQakRHs6IGbwtaliLuUqj7pM6_gCcSbZKSnwoUXaVZnZ8s2pWz4nVxTDnovHjjTxCdpu7p6jly4-Z7ddBH5G43k7ffP8Y-3tXL-o-4UJ-7LJj62Oh7oHEOybzfuw1_MfbDN6N_utlGwuWh1rvy97Niu9Oy9EUS4E7t4N6aJkYUlZptNY7bxC8i0MixcXnCxkH0PtRv7_P-eubpTl2VyKFBEqPyOTsq_8rpiSLsD0vVNjMB7erjSihoNpTI24891t-TOgkJIb2ShZjzUVeyJQQqT2f7tkLCGYRfr8E2V60evnzFw5jyvr06tmImsAvtL2jvo7mduggeusz7HGu6KpUL1APfqt9hFqeJIZOYxUWc-lk7tJ_-XVvhiU4ouwDsKJB2ZhSt_leZODmmE76JtFBSnikFv4ip4YHlt3TwVcsjYF5QNTDHV_rX7CupI2ODOV9gV7MHE6bEkd06edUP7kclqojT1H4C1OUS_yHqHu-3aEpOkuzF-65KGLZ9gRIQL-hyLAGcKJw=w300-h225-no)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 6.3 - Basic Task and Data Parallelism\n",
    "In this example, we will design one workflow that exhibits task parallelism and one that exhibits data parallelism. A nice real world example of something that is not only a nice simple example of a parallel algorithm, but is also an example of a real world computational task that comes up a lot is matrix multiplication. This image displays how data parallelism comes up in the problem of matrix multiplication:\n",
    "![matrix multiplication data parallelism](https://upload.wikimedia.org/wikipedia/commons/6/68/Data_Parallelism_in_matrix_multiplication.png)\n",
    "You can see that after some basic serial computation to set up the right hand side of the equation, each one of the cells in the 3x2 matrix on the right hand side can be computed concurrently or in parallel. This represents data parallelism because the same multiplication and addition instructions are used on different data. For multiplication, we can divide matrix A and B into blocks along rows and columns respectively. This allows us to calculate every element in matrix C individually thereby making the task parallel. \n",
    "\n",
    "For example: $A(m \\times n) \\cdot B (n \\times k)$ can be finished in $O(n)$ time instead of $O(m*n*k)$, when executed in parallel using $m*k$ processors. (This is the example from the image).\n",
    "\n",
    "As our example of task parallelism, we are going to use a very simple workflow - computing the square root of many numbers in parallel and computing the square of many numbers concurrently. This way, we have a single instruction (sqrt) which is run on some data, and a second instruction (square) which is run on other data. Based on the fact that all of the input data is in one data pool, a list, this is task parallelism because we are applying different tasks to data in different parts of the same data pool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data parallel matrix multiplication\n",
    "import random\n",
    "import multiprocessing\n",
    "from itertools import starmap, repeat\n",
    "from operator import mul\n",
    "\n",
    "def calc_row_of_product_matrix(a_row, b, izip=zip):\n",
    "    '''Calculate a row of the product matrix P = A * B\n",
    "    Arguments:\n",
    "      a_row is af A\n",
    "      b is the B matrix\n",
    "    returns the corresponding row of P matrix'''\n",
    "    return map(lambda col: sum(starmap(mul,zip(a_row,col))), zip(*b))\n",
    "\n",
    "def eval_func_tuple(f_args):\n",
    "    '''Takes a tuple of a function and args, evaluates and returns result'''\n",
    "    return f_args[0](*f_args[1:])\n",
    "\n",
    "class multimatrix(list):\n",
    "\n",
    "    def __mul__(self, b, izip=zip, repeat=repeat):\n",
    "        '''Concurrent matrix multiplication with multiprocessing.Pool. '''\n",
    "        pool = multiprocessing.Pool(multiprocessing.cpu_count())\n",
    "        return pool.map(eval_func_tuple, izip(repeat(calc_row_of_product_matrix), self, repeat(b))) \n",
    "\n",
    "class itermatrix(list):\n",
    "\n",
    "    @staticmethod\n",
    "    def sumprod(row, col, sum=sum, starmap=starmap, mul=mul):\n",
    "        '''Returns sumproduct of two vectors.'''\n",
    "        return sum(starmap(mul,zip(row,col)))\n",
    "\n",
    "    def __mul__(self, b, imap=map, izip=zip):\n",
    "        '''Matrix multiplication returning iterable of iterables'''\n",
    "        return map(lambda row: map(lambda col: itermatrix.sumprod(row, col), zip(*b)), self)\n",
    "\n",
    "def iterate_results(result):\n",
    "    '''Iterate over iterable of iterables,\n",
    "    and returns elements in list of lists.\n",
    "    Usage: if you want to run the whole calculation at once:\n",
    "    p = iterate_results(itermatrix([[1, 3], [-5, 6], [2, 4]]) * itermatrix([[1, 4], [8, 7]]))'''\n",
    "    return[[col for col in row] for row in result]\n",
    "\n",
    "def random_v(K=1000,min=-1000,max=1000):\n",
    "    '''Generates a random vector of dimension N;\n",
    "    Returns a list of integers.\n",
    "    The values are integers in the range [min,max].'''\n",
    "    return [random.randint(min,max) for k in range(K)]\n",
    "\n",
    "def random_m(N=1000, K=1000):\n",
    "    '''Generates random matrix. Returns list of list of integers.'''\n",
    "    return [random_v(K) for n in range(N)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[25, 25], [43, 22], [34, 36]] iter test OK\n",
      "CPU times: user 43 µs, sys: 135 µs, total: 178 µs\n",
      "Wall time: 146 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    a = [[1, 3], [-5, 6], [2, 4]]\n",
    "    b = [[1, 4], [8, 7]]\n",
    "    adotb = [[25, 25], [43, 22], [34, 36]]\n",
    "    #A = multimatrix(a)\n",
    "    #B = multimatrix(b)\n",
    "    #prod = A * B\n",
    "    #assert(adotb == prod)\n",
    "    #print(prod, \"multi test OK\")\n",
    "    A = itermatrix(a)\n",
    "    B = itermatrix(b)\n",
    "    iterprod = A * B\n",
    "    listprod = iterate_results(iterprod)\n",
    "    assert(adotb == listprod)\n",
    "    print(listprod, \"iter test OK\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0, 1.4142135623730951, 2.0, 2.449489742783178, 2.8284271247461903, 3.1622776601683795, 3.4641016151377544, 3.7416573867739413, 4.0]\n",
      "[1, 9, 25, 49, 81, 121, 169, 225, 289]\n",
      "[0.0, 1.4142135623730951, 2.0, 2.449489742783178, 2.8284271247461903, 3.1622776601683795, 3.4641016151377544, 3.7416573867739413, 4.0]\n",
      "[1, 9, 25, 49, 81, 121, 169, 225, 289]\n",
      "CPU times: user 984 ms, sys: 537 ms, total: 1.52 s\n",
      "Wall time: 1.38 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Task parallel suqare roots and squares\n",
    "import multiprocessing\n",
    "\n",
    "# I know there's a library function for this\n",
    "# I want it to be a bit slower so we can see the\n",
    "# performance benefits of the parallelism\n",
    "def sqrt(n):\n",
    "    return n**(0.5)\n",
    "\n",
    "def square(n):\n",
    "    return n**2\n",
    "\n",
    "run_serial = False # switch to True if you want to see it run slowly...\n",
    "if run_serial:\n",
    "    sqrts = []\n",
    "    squares = []\n",
    "    for i in range(1000000):\n",
    "        if i%2 == 0:\n",
    "            sqrts.append(sqrt(i))\n",
    "        else:\n",
    "            squares.append(square(i))\n",
    "    print(sqrts[0:9])\n",
    "    print(squares[0:9])\n",
    "                           \n",
    "\n",
    "pool = multiprocessing.Pool(32)\n",
    "sqrts = pool.map(sqrt, [i for i in range(1000000) if i%2==0])\n",
    "squares = pool.map(square, [i for i in range(1000000) if i%2==1])\n",
    "print(sqrts[0:9])\n",
    "print(squares[0:9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IV. Maximum Speedups\n",
    "By now, you should understand that different types of algorithms can be helped by parallelization different amounts. Let's go ahead and quantify that. First, let's do some mathematical analysis and figure out how well we could do. The first major breakthrough in terms of parallel algorithm analysis is called [Amdahl's Law](https://en.wikipedia.org/wiki/Amdahl%27s_law). Formulated by Gene Amdahl in 1967, Amdahl's Law allows us to compute theoretical maximum speedups based on how parallelizable our tasks are. In order to figure out your theoretical maximum, you need to know how many processors you have access to and what percentage of your program can be parallelized. For example, if a program needs 20 hours using a single processor core, and a particular part of the program which takes one hour to execute cannot be parallelized, while the remaining 19 hours (p = 0.95, where p is the fraction of the code that can be parallelized) of execution time can be parallelized, then regardless of how many processors are devoted to a parallelized execution of this program, the minimum execution time cannot be less than that critical one hour. Hence, the theoretical speedup is limited to at most 20 times (1/(1 − p) = 20). For this reason, parallel computing with many processors is useful only for highly parallelizable programs. This image sums up Amdahl's law for various values of p.\n",
    "![amdahl's law](https://upload.wikimedia.org/wikipedia/commons/thumb/e/ea/AmdahlsLaw.svg/640px-AmdahlsLaw.svg.png)\n",
    "\n",
    "Another parallel algorithm analysis law, the [Gustafson–Barsis Law](https://en.wikipedia.org/wiki/Gustafson%27s_law), gives the theoretical speedup in latency of the execution of a task at fixed execution time that can be expected of a system whose resources are improved. This law considers not how much faster your code can get, but instead how much more you can do with the same code in the same amount of time as you scale up your computer. Amdahl's law presupposes that the computing requirements will stay the same, given increased processing power. In other words, an analysis of the same data will take less time given more computing power.\n",
    "\n",
    "Gustafson, on the other hand, argues that more computing power will cause the data to be more carefully and fully analyzed: pixel by pixel or unit by unit, rather than on a larger scale. Where it would not have been possible or practical to simulate the impact of nuclear detonation on every building, car, and their contents (including furniture, structure strength, etc.) because such a calculation would have taken more time than was available to provide an answer, the increase in computing power will prompt researchers to add more data to more fully simulate more variables, giving a more accurate result.\n",
    "\n",
    "\n",
    "Amdahl's Law reveals a limitation in, for example, the ability of multiple cores to reduce the time it takes for a computer to boot to its operating system and be ready for use. Assuming the boot process was mostly parallel, quadrupling computing power on a system that took one minute to load might reduce the boot time to just over fifteen seconds. But greater and greater parallelization would eventually fail to make bootup go any faster, if any part of the boot process were inherently sequential.\n",
    "\n",
    "Gustafson's law argues that a fourfold increase in computing power would instead lead to a similar increase in expectations of what the system will be capable of. If the one-minute load time is acceptable to most users, then that is a starting point from which to increase the features and functions of the system. The time taken to boot to the operating system will be the same, i.e. one minute, but the new system would include more graphical or user-friendly features.\n",
    "\n",
    "This image sums up Gustafson's law for various values of p (same p as earlier):\n",
    "![Theoretical maximum speedup](https://cdn.comsol.com/wordpress/2014/03/Graph-depicting-how-the-size-of-the-job-increases-with-the-number-of-available-processes.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 6.4 - Timing Fibonacci and Monte Carlo\n",
    "Given all of that theory we just learned, what do you think the value of _p_ is for the Monte Carlo pi simulation? For the Fibonacci sequence generator? Think about it for a moment. Read ahead when you've got a guess. For the Monte Carlo Pi, the value of _p_ is very close to 1, while for the Fibonacci sequence generator, the value of _p_ is very close to 0. Hopefully the description and examples of those tasks make those values believable to you. In case they don't, we're going to time them right now, in this example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining parallel monte carlo pi calculation \n",
    "\n",
    "import random\n",
    "import multiprocessing\n",
    "from multiprocessing import Pool\n",
    "\n",
    "\n",
    "#caculate the number of points in the unit circle\n",
    "#out of n points\n",
    "def monte_carlo_pi_part(n):\n",
    "    \n",
    "    count = 0\n",
    "    for i in range(n):\n",
    "        x=random.random()\n",
    "        y=random.random()\n",
    "        \n",
    "        # if it is within the unit circle\n",
    "        if x*x + y*y <= 1:\n",
    "            count=count+1\n",
    "        \n",
    "    #return\n",
    "    return count  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Esitmated value of Pi: 3.1408732 \n",
      "CPU times: user 83.9 ms, sys: 674 ms, total: 758 ms\n",
      "Wall time: 908 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# Calling parallel monte carlo pi simulation\n",
    "\n",
    "# Nummber of points to use for the Pi estimation\n",
    "n = 10000000\n",
    "np = multiprocessing.cpu_count()\n",
    "# iterable with a list of points to generate in each worker\n",
    "# each worker process gets n/np number of points to calculate Pi from\n",
    "\n",
    "part_count=[n//np for i in range(np)]\n",
    "\n",
    "#Create the worker pool\n",
    "# http://docs.python.org/library/multiprocessing.html#module-multiprocessing.pool\n",
    "pool = Pool(processes=np)   \n",
    "\n",
    "# parallel map\n",
    "count=pool.map(monte_carlo_pi_part, part_count)\n",
    "\n",
    "print(\"Esitmated value of Pi: {} \".format(sum(count)/(n*1.0)*4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 573 ms, sys: 3.8 s, total: 4.38 s\n",
      "Wall time: 18.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import multiprocessing\n",
    "import time\n",
    "\n",
    "# Fibonacci sequence on 32 cores\n",
    "\n",
    "# Compute Nth fibonacci number\n",
    "def fib(n):\n",
    "    if n==0 or n==1:\n",
    "        return 1\n",
    "    return fib(n-1) + fib(n-2)\n",
    "\n",
    "def calc_single_fib(i):\n",
    "    before = time.time()\n",
    "    res = fib(i)\n",
    "    after = time.time()\n",
    "    tot = after-before\n",
    "\n",
    "pool = multiprocessing.Pool(32)\n",
    "tasks = pool.map(calc_single_fib, range(0, 40, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that even though we run each distinct fibonacci number in parallel, we can't parallelize the underlying function at all, so it's still very slow."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 6. Pi Over Many Nodes or Fibonacci Over Many Nodes\n",
    "Using Parsl, implement a fibonacci sequence generator and a mote carlo pi simulation. Run it on the entire cluster. Which one runs faster? Which one seems to scale better?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here"
   ]
  }
 ],
 "metadata": {
  "author": "mes",
  "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"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}