monte-carlo-pi.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of active threads running on client init 5\n"
     ]
    }
   ],
   "source": [
    "import ipyparallel as ipp\n",
    "rc = ipp.Client()\n",
    "engines = rc[:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Approximate π with random $x,y$ samples on a unit square:\n",
    "\n",
    "\n",
    "$$\n",
    "x \\in [0, 1],\\\\\n",
    "y \\in [0, 1]\n",
    "$$\n",
    "\n",
    "where the probability that any given point $x,y$ is inside the unit circle is\n",
    "the relative area of the circle within the square:\n",
    "\n",
    "$$\n",
    "P(x^2 + y^2 \\leq 1) = \\frac{\\pi}{4}\n",
    "$$\n",
    "\n",
    "so we can approximate π by collecting random samples on the unit square and counting how many are inside the unit circle:\n",
    "\n",
    "$$\n",
    "\\pi \\approx \\frac{4}{N}  \\sum_{i=1}^N \\left( x_i^2 + y_i^2 \\right) \\leq 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mc_pi(n):\n",
    "    \"\"\"Monte Carlo approximation of π\n",
    "    \n",
    "    Throw darts uniformly distributed on a square,\n",
    "    count how many land in circumscribed circle.\n",
    "\n",
    "    Fraction inside circle approaches π / 4\n",
    "    \"\"\"\n",
    "    import random\n",
    "    samples = []\n",
    "    for i in range(n):\n",
    "        x = random.random()\n",
    "        y = random.random()\n",
    "        in_circle = (x * x) + (y * y) <= 1\n",
    "        samples.append(in_circle)\n",
    "    return 4 * sum(samples) / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run our tests for a series of sample sizes in serial and parallel,\n",
    "comparing the times of each to see how the parallel performance\n",
    "improves over serial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using 2 processes\n",
      "\n",
      "Monte Carlo sampling of π: 100 samples (50 per engine)\n",
      "serial:\n",
      "33.3 µs ± 2.79 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "parallel:\n",
      "7.99 ms ± 731 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "speedup: 0.00x\n",
      "\n",
      "Monte Carlo sampling of π: 1000 samples (500 per engine)\n",
      "serial:\n",
      "320 µs ± 7.52 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "parallel:\n",
      "8.37 ms ± 941 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "speedup: 0.04x\n",
      "\n",
      "Monte Carlo sampling of π: 10000 samples (5000 per engine)\n",
      "serial:\n",
      "3.29 ms ± 185 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "parallel:\n",
      "15.2 ms ± 4.34 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "speedup: 0.22x\n",
      "\n",
      "Monte Carlo sampling of π: 100000 samples (50000 per engine)\n",
      "serial:\n",
      "35.4 ms ± 2.07 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "parallel:\n",
      "24.3 ms ± 362 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "speedup: 1.45x\n",
      "\n",
      "Monte Carlo sampling of π: 1000000 samples (500000 per engine)\n",
      "serial:\n",
      "326 ms ± 3.78 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "parallel:\n",
      "193 ms ± 7.08 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "speedup: 1.69x\n"
     ]
    }
   ],
   "source": [
    "n_engines = len(engines)\n",
    "\n",
    "print(f\"Using {n_engines} processes\")\n",
    "\n",
    "for n_samples in [100, 1000, 10_000, 100_000, 1_000_000]:\n",
    "    if n_samples % len(engines):\n",
    "        n_samples += len(engines)\n",
    "    samples_per_engine = n_samples // n_engines\n",
    "    print(f\"\\nMonte Carlo sampling of π: {n_samples} samples ({samples_per_engine} per engine)\")\n",
    "    print(\"serial:\")\n",
    "    tr_serial = %timeit -o mc_pi(n_samples)\n",
    "    print(\"parallel:\")\n",
    "    tr_parallel = %timeit -o sum(engines.apply(mc_pi, samples_per_engine)) / len(engines)\n",
    "    print(f\"speedup: {tr_serial.average / tr_parallel.average:.2f}x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that in this example with four local engines that the parallel implementation\n",
    "is slower up to about 100,000 samples, where the serial case takes about 30 milliseconds.\n",
    "\n",
    "The overhead of scheduling and waiting for a single parallel task\n",
    "means that parallelism on tasks quicker than several milliseconds will overwhelm the benefit of the the concurrent workload.\n",
    "\n",
    "Tasks must take long enough for the overhead of managing parallel and distributed tasks to be worth it.\n",
    "In the case of IPython Parallel, that's on the order of 10s of milliseconds per task.\n",
    "If there are many tasks that can be queued concurrently, the overhead can be pipelined and shared,\n",
    "reducing the effective per-task overhead."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}