5-dask.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Dask & multiprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "```\n",
    "conda install dask distributed -c conda-forge\n",
    "```\n",
    "\n",
    "_(and numpy, matplotlib from before)_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Computers aren't getting any _faster._\n",
    "\n",
    "<center><img src=\"img/clock-rate.jpg\" width=\"60%\" /><div style=\"font-size: 14px\"><i>Computer Architecture: A Quantitative Approach,</i> David A. Patterson and John L. Hennessy</div></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "But Moore's Law is still in effect: the number of transistors per square inch continues to grow exponentially (for now). In the 21st century, however, those extra transitors are used to make more execution units, not to incrase the rate through smaller pipelines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "(The scaling of clock rates— Dennard's Law— ended because power dissipation scales with clock rate squared: anything faster than 3 GHz _cooks_ the chip!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Most programming languages, Python among them, cannot be transparently parallelized. You'll have to change your programs to use the extra processors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's take a moderately complex problem as an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "import time\n",
    "\n",
    "def prepare(height, width):\n",
    "    y, x = numpy.ogrid[-1:0:height*1j, -1.5:0:width*1j]\n",
    "    c = x + y*1j\n",
    "    fractal = numpy.zeros(c.shape, dtype=numpy.int32)\n",
    "    return c, fractal\n",
    "\n",
    "def run(c, fractal, maxiterations=20):\n",
    "    fractal *= 0                  # set fractal to maxiterations without replacing it\n",
    "    fractal += maxiterations\n",
    "    z = c\n",
    "    for i in range(maxiterations):\n",
    "        z = z**2 + c\n",
    "        diverge = numpy.absolute(z) > 2\n",
    "        divnow = diverge & (fractal == maxiterations)\n",
    "        fractal[divnow] = i\n",
    "        z[diverge] = 2\n",
    "    return fractal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "c, fractal = prepare(8000, 12000)\n",
    "\n",
    "starttime = time.time()\n",
    "fractal = run(c, fractal)\n",
    "time.time() - starttime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "ax.imshow(fractal)\n",
    "# ax.imshow(fractal[-200:, :300])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Python has built-in libraries for parallel processing:\n",
    "\n",
    "   * **threading:** lets you launch individual threads; you manage coordination.\n",
    "   * **multiprocessing:** same interface but it launches processes. Pro: can't make common mistakes due to shared memory. Con: memory isn't shared; have to ship data to and from workers.\n",
    "   * **concurrent.futures:** higher-level interface: Python manages workers; you send work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Here's an illustration of the threading interface. Since memory is shared, we don't have to send data to the workers or send results back— they can all see and modify the same array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import threading\n",
    "\n",
    "class Worker(threading.Thread):\n",
    "    def __init__(self, c, fractal, i):\n",
    "        super(Worker, self).__init__()\n",
    "        self.c, self.fractal, self.i = c, fractal, i\n",
    "    def run(self):\n",
    "        run(self.c[10*self.i : 10*(self.i + 1), :], self.fractal[10*self.i : 10*(self.i + 1), :])\n",
    "\n",
    "c, fractal = prepare(8000, 12000)\n",
    "workers = []\n",
    "for i in range(800):\n",
    "    workers.append(Worker(c, fractal, i))\n",
    "\n",
    "starttime = time.time()\n",
    "\n",
    "for worker in workers:\n",
    "    worker.start()\n",
    "for worker in workers:\n",
    "    worker.join()\n",
    "\n",
    "time.time() - starttime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now we have to check the result because it's easy to screw this up. (I did many times, preparing this talk.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "ax.imshow(fractal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "1 thread took 35 seconds to complete.\n",
    "\n",
    "8 threads took 12 seconds to complete."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "35 / 12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "3 ≠ 8.\n",
    "\n",
    "It's often difficult to get \"perfect scaling,\" N times more work from N threads, in real situations. Even though this problem is \"embarrassingly parallel\" (none of the workers need to know other workers' results), there can be scheduling overhead, contention for memory, or slow-downs due to Python's [Global Interpreter Lock](https://realpython.com/python-gil/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "One way to avoid the global interpreter lock is to send work to separate processes. Python interpreters in separate processes do not share memory and therefore do not need to coordinate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "However, that means that we can't send data by simply sharing variables. We have to send it through a `multiprocessing.Queue` (which serializes— pickles— the data so that it can go through a pipe)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "...usually. There's an exception to this: you can share arrays among processes if you declare them as shared memory before launching the subprocesses. Python has a special type for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiprocessing\n",
    "import ctypes\n",
    "\n",
    "sharedarray = multiprocessing.RawArray(ctypes.c_double, 100)\n",
    "sharedarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This is not a Numpy array, but it can be cast as a Numpy array (in the forked process) like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numpy.frombuffer(sharedarray, dtype=numpy.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "If the forked processes are not writing to different parts of the array, they can seriously garble the data if they write to the same element at the same time.\n",
    "\n",
    "\n",
    "It's not for the faint of heart, but it can be the fastest way to communicate between processes, and seperate processes are the only way to fully escape synchronization delays due to Python's global interpreter lock."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "By now, you may be wondering if there's a more \"high level\" approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Python 3 introduced an \"executor\" interface that manages workers for you. Instead of creating threads or processes with a `run` method, you create an executor and send work to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import concurrent.futures\n",
    "executor = concurrent.futures.ThreadPoolExecutor(max_workers=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c, fractal = prepare(8000, 12000)\n",
    "# fractal = numpy.asfortranarray(fractal)\n",
    "\n",
    "def piece(i):\n",
    "    ci =       c[10*i : 10*(i + 1), :]\n",
    "    fi = fractal[10*i : 10*(i + 1), :]\n",
    "    run(ci, fi)\n",
    "\n",
    "starttime = time.time()\n",
    "\n",
    "futures = executor.map(piece, range(800))\n",
    "for future in futures:         # iterating over them waits for the results\n",
    "    pass\n",
    "\n",
    "time.time() - starttime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Yay! A tiny bit better! What happens when we change to Fortran order? Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Always make sure we haven't screwed things up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "ax.imshow(fractal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Still, there needs to be a better way. Our array slices in `piece` are fragile: an indexing error can ruin the result. Can't the problem of scattering work be generalized?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import dask.array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c, fractal = prepare(8000, 12000)\n",
    "\n",
    "c = dask.array.from_array(c, chunks=(10, 12000))\n",
    "fractal = dask.array.from_array(fractal, chunks=(10, 12000))\n",
    "\n",
    "starttime = time.time()\n",
    "fractal = run(c, fractal)\n",
    "time.time() - starttime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "That was too fast: too good to be true."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fractal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This is not an array: it is a description of how to make an array. Dask has stepped through our procedure and built an execution graph, encoding all the dependencies so that it can correctly apply it to individual chunks. When we execute this graph, Dask will send a chunk to each processor in the computer and combine results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "starttime = time.time()\n",
    "fractal = fractal.compute()    # replace `fractal` the execution graph with `fractal` the array result\n",
    "time.time() - starttime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Now this check is a formality: Dask has managed the chunking, so we won't accidentally miss a slice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "ax.imshow(fractal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We seem to have paid for this simplicity: it took twice as long as the carefully sliced `pieces` in the executor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The reason is that our code is not as simple as it looks. It has masking and piecemeal assignments, which in principle could introduce complex dependencies. _We_ know that everything will be fine if you just chop up the array in independent sections— and thus we implemented our thread and executor-based solutions that way."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let me show you what Dask has to do for a 1×1 chunking of our problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "c, fractal = prepare(1, 1)                                 # try 2, 2\n",
    "c = dask.array.from_array(c, chunks=(1, 1))\n",
    "fractal = dask.array.from_array(fractal, chunks=(1, 1))\n",
    "fractal = run(c, fractal, maxiterations=1)                 # try more iterations\n",
    "fractal.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "If that were all, I'd probably stick to chopping up the grid by hand (when possible). However, _exactly the same interface_ that distributes work across cores in my laptop can distribute work around the world, just by pointing it to a remote scheduler.\n",
    "\n",
    "This is truly the ~~lazy~~ busy researcher approach!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Note to self: launch\n",
    "\n",
    "`dask-scheduler &`\n",
    "\n",
    "and\n",
    "\n",
    "`dask-worker --nthreads 8 127.0.0.1:8786 &`\n",
    "\n",
    "in a terminal now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import dask.distributed\n",
    "client = dask.distributed.Client(\"127.0.0.1:8786\")\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c, fractal = prepare(8000, 12000)\n",
    "\n",
    "c = dask.array.from_array(c, chunks=(100, 12000))\n",
    "fractal = dask.array.from_array(fractal, chunks=(100, 12000))\n",
    "fractal = run(c, fractal)\n",
    "\n",
    "starttime = time.time()\n",
    "fractal = client.compute(fractal, sync=True)\n",
    "time.time() - starttime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Well, that was exciting!\n",
    "\n",
    "In the end, this example took longer than the single-core version, but it illustrates how array operations _can be_ distributed in a simple way."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "I haven't shown very much of what Dask can do. It's a general toolkit for delayed and distributed evaluation. As such, it provides a nice way to work on Pandas-like DataFrames that are too large for memory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dask.dataframe\n",
    "\n",
    "df = dask.dataframe.read_csv(\"data/nasa-exoplanets.csv\")\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We don't see the data because they haven't been loaded. But we can get them if we need them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df[[\"pl_hostname\", \"pl_pnum\"]].compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Additionally, Dask isn't the only project filling this need. There's also:\n",
    "\n",
    "   * **Joblib:** annotate functions to execute remotely with decorators.\n",
    "   * **Parsl:** same, but work with conventional schedulers (Condor, Slurm, GRID); an academic project.\n",
    "   * **PySpark:** Spark is a big, scalable project, though its Python interface has performance issues.\n",
    "\n",
    "and many smaller projects.\n",
    "\n",
    "(Distributed computing hasn't been fully figured out yet.)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}