pi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Profiling in Python\n",
    "\n",
    "There are several different ways to measure resource usage in Python applications, which help identify where an application should be parallelized and what kinds of design patterns from Ray core to use.\n",
    "to show how to determine where to parallelize.\n",
    "We'll introduce profiling methods to measure memory allocation, objects creation, deterministic profiling (call stack), sampling, etc., plus different means of visualization.\n",
    "\n",
    "The following example uses a Monte Carlo method for approximating the value of π, based on the tutorial [\"Ray Crash Course - Tasks\"](https://github.com/anyscale/academy/blob/main/ray-crash-course/01-Ray-Tasks.ipynb) by Dean Wampler.\n",
    "\n",
    "We'll compare/contrast the serial implementation versus use of *remote functions* to help parallelize this application.\n",
    "Then we'll analysze the per-CPU speedup, and also consider the overhead costs of using Ray core features -- to understand more about the trade-offs being made."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext watermark\n",
    "%watermark -v -m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Launch Ray\n",
    "\n",
    "First, start Ray and open its dashboard in another browser tab:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "import ray\n",
    "\n",
    "ray.init(\n",
    "    ignore_reinit_error=True,\n",
    "    logging_level=logging.ERROR,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application Code\n",
    "\n",
    "Define a simple function that uses a stochastic method of approximating π, repeated through some number of samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statistics\n",
    "import time\n",
    "\n",
    "def estimate_pi (num_samples):\n",
    "    xs = np.random.uniform(low=-1.0, high=1.0, size=num_samples)   # generate num_samples random samples for the x coordinate\n",
    "    ys = np.random.uniform(low=-1.0, high=1.0, size=num_samples)   # generate num_samples random samples for the y coordinate\n",
    "    xys = np.stack((xs, ys), axis=-1)                              # similar to Python's \"zip(a,b)\"; creates np.array([(x1,y1), (x2,y2), ...]).\n",
    "\n",
    "    inside = (xs**2.0 + ys**2.0) <= 1.0                            # create a predicate over all the array elements\n",
    "    xys_inside = xys[inside]                                       # select only those  elements inside the circle\n",
    "    in_circle = xys_inside.shape[0]                                # return the number of elements inside the circle\n",
    "    approx_pi = 4.0 * in_circle / num_samples                      # the Pi estimate\n",
    "\n",
    "    return approx_pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this approximation requires many samples, which are independent, a *task-parallel* pattern can be applied. We'll create a *remote function* version of the sampling function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@ray.remote\n",
    "def distrib_estimate_pi (num_samples):\n",
    "    return estimate_pi(num_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another function will collect measures for one epoch, i.e., some number of trials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_epoch (num_samples, num_trials, distrib=False):\n",
    "    start = time.time()\n",
    "    \n",
    "    if distrib:\n",
    "        refs = [distrib_estimate_pi.remote(num_samples) for _ in range(num_trials)]\n",
    "        pis = ray.get(refs)\n",
    "    else:\n",
    "        pis = [estimate_pi(num_samples) for _ in range(num_trials)]\n",
    "\n",
    "    # measure CPU time for the code section parallelized as a remote function\n",
    "    duration = time.time() - start\n",
    "\n",
    "    approx_pi = statistics.mean(pis)\n",
    "    stdev = statistics.stdev(pis)\n",
    "    error = 100.0 * abs(approx_pi - np.pi) / np.pi\n",
    "\n",
    "    return num_samples, duration, approx_pi, stdev, error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define class to manage the simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class Sim:\n",
    "    DF_COL_NAMES = [\"n\", \"duration\", \"approx_pi\", \"stdev\", \"error\"]\n",
    "    NUM_SAMPLES = 1000000\n",
    "    NUM_TRIALS = 20\n",
    "    STEP_SIZE = int(NUM_SAMPLES / 25)\n",
    "\n",
    "    def __init__ (self, distrib=False, num_samples=NUM_SAMPLES, num_trials=NUM_TRIALS, step_size=STEP_SIZE):\n",
    "        self.distrib = distrib\n",
    "        self.num_samples = num_samples\n",
    "        self.num_trials = num_trials\n",
    "        self.step_size = step_size\n",
    "        self.df = None\n",
    "\n",
    "\n",
    "    def run (self):\n",
    "        # use a minimum of 2 trials, to be able to calculate standard deviation\n",
    "        results = [\n",
    "            run_epoch(n_samp, self.num_trials, distrib=self.distrib)\n",
    "            for n_samp in range(2, self.num_samples, self.step_size)\n",
    "        ]\n",
    "    \n",
    "        self.df = pd.DataFrame(results, columns=self.DF_COL_NAMES)\n",
    "        return self\n",
    "\n",
    "\n",
    "    def plot (self):\n",
    "        plt.plot(\"n\", \"duration\", data=self.df, color=\"green\", linewidth=1, linestyle=\"dashed\")\n",
    "        plt.plot(\"n\", \"error\", data=self.df, color=\"red\", linewidth=0.5, linestyle=\"dashed\")\n",
    "        plt.plot(\"n\", \"stdev\", data=self.df, color=\"blue\", linewidth=2)\n",
    "\n",
    "        plt.yscale(\"log\")\n",
    "        plt.legend()\n",
    "        plt.show()\n",
    "        return self"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Profiling Tools\n",
    "\n",
    "Next, we'll set up to use the following tools for profiling in Python:\n",
    "\n",
    "  * [`objgraph`](https://mg.pov.lt/objgraph/)\n",
    "  * [`tracemalloc`](https://docs.python.org/3/library/tracemalloc.html)\n",
    "  * [`prun`](https://ipython.readthedocs.io/en/stable/interactive/magics.html#magic-prun)\n",
    "  * [`snakeviz`](https://jiffyclub.github.io/snakeviz/)\n",
    "  * [`pyinstrument`](https://github.com/joerick/pyinstrument/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import objgraph\n",
    "import tracemalloc\n",
    "import pyinstrument\n",
    "\n",
    "%load_ext snakeviz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Object Creation\n",
    "\n",
    "The first method of profiling shown here uses\n",
    "[`objgraph`](https://mg.pov.lt/objgraph/)\n",
    "to examine a before/after contrast of what objects are getting generated by the application.\n",
    "While other profiling methods can be repeated, this kind of analysis requires some special handling (isolation) and should be run first -- before other tools make the runtime environment confusing to analyze.\n",
    "The objectives here are to understand:\n",
    "\n",
    "  1. Which kinds of objects are growing (by count)\n",
    "  2. How the different objects refer to each other\n",
    "\n",
    "That can help identify if there are problems with allocating many objects, which might require the application code to be reworked.\n",
    "\n",
    "First we'll measure a baseline of the top 10 kinds of objects that have been created so far:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objgraph.show_growth(limit=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then run the simulation serially **once**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "sim_s = Sim(distrib=False, num_samples=Sim.NUM_SAMPLES)\n",
    "sim_s.run();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can measure the object counts again, and examine the growth (by contrast):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objgraph.show_growth()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objgraph.show_most_common_types() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The top few categories are function calls, dictionaries, and tuples – which is expected, given how our π approximations create lots of NumPy arrays.\n",
    "As expected, there's nothing much there to worry about – although in applications which large memory use we might need to parallelize so that each unit of work had more available memory.\n",
    "\n",
    "Next, let's examine a graph of how object make reference to other objects.\n",
    "This can be useful for tracing potential memory leaks..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roots = objgraph.get_leaking_objects()\n",
    "objgraph.show_refs(roots[:3], refcounts=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Memory Allocation\n",
    "\n",
    "The next profiling method uses\n",
    "[`tracemalloc`](https://docs.python.org/3/library/tracemalloc.html)\n",
    "to trace the memory blocks allocated by Python.\n",
    "This analysis provides the following information:\n",
    "\n",
    "  * traceback to where an object got allocated\n",
    "  * statistics about allocated memory blocks per filename, per line: total size, number and average size\n",
    "  * computing the differences between two snapshots to detect memory leaks\n",
    "  \n",
    "We'll run the application serially to capture a snapshot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tracemalloc.start()\n",
    "\n",
    "Sim(distrib=False, num_samples=Sim.NUM_SAMPLES).run()\n",
    "\n",
    "snapshot = tracemalloc.take_snapshot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "top_stats = snapshot.statistics('traceback')\n",
    "\n",
    "# pick the biggest memory block\n",
    "for stat in top_stats:\n",
    "    print(\"%s memory blocks: %.1f KiB\" % (stat.count, stat.size / 1024))\n",
    "\n",
    "    for line in stat.traceback.format():\n",
    "        print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There should not be any *memory leaks* in this application."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deterministic Profiling\n",
    "\n",
    "The\n",
    "[`%%prun` magic](https://ipython.readthedocs.io/en/stable/interactive/magics.html#magic-prun)\n",
    "invokes the\n",
    "[cProfile](https://docs.python.org/3/library/profile.html)\n",
    "*deterministic* profiler in Python to trace how often and for what duration the different functions get called.\n",
    "In other words, we'll track of the call stack statistics to understand more about CPU use.\n",
    "This configures the profiler to save data to the file `prof_cpu.txt` as text.\n",
    "\n",
    "The [`%%snakeviz` magic](https://jiffyclub.github.io/snakeviz/)\n",
    "then analyzes that profiling data and generates an interactive report.\n",
    "\n",
    "Note that *deterministic profiling* creates some overhead which can distort the measurements in relatively small code blocks or short-running code, although it's generally find for long-running programs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%prun -q -T prof_cpu.txt\n",
    "%%snakeviz\n",
    "\n",
    "Sim(distrib=False, num_samples=Sim.NUM_SAMPLES).run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *icicle* interactive chart illustrates how so very much of the overall CPU cost is in the `estimate_pi()` method.\n",
    "\n",
    "In some environments, there may be security checks that limit use of `SnakeViz` (pending updates) so here's a screenshot:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"snekviz.png\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sampling Profiler\n",
    "\n",
    "The\n",
    "[`pyinstrument`](https://github.com/joerick/pyinstrument/)\n",
    "library provides a *sampling profiler* – an alternative way to measure the call stack and CPU use.\n",
    "While not *exact*, this approach creates less overhead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "profiler = pyinstrument.Profiler()\n",
    "profiler.start()\n",
    "\n",
    "Sim(distrib=False, num_samples=Sim.NUM_SAMPLES).run()\n",
    "\n",
    "profiler.stop()\n",
    "display(HTML(profiler.output_html()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the `estimate_pi()` function is where most of the CPU cost occurs, although we can get clearer measures for how much overhead there is among the rest of code.\n",
    "In this case the cumulative overhead is relatively small, probably measured in milliseconds."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results of Serial Execution\n",
    "\n",
    "Now let's look at a visualization of the *number of samples* plotted versus *duration*, *error*, and *standard deviation*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_s.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the *standard deviation* drops quickly.\n",
    "We can calculate that measure without even knowing an \"exact\" value of π (the *error* measure), so we could potentially rework this application to have an *early termination* by setting a threshold on the *stdev* measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_s.df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results of Parallel Execution\n",
    "\n",
    "Now let's run again and this time parallelize the application using *remote functions*, by setting the `distrib=True` flag.\n",
    "We'll use the `%%time` magic to show the \"wall clock\" CPU time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "sim_d = Sim(distrib=True, num_samples=Sim.NUM_SAMPLES)\n",
    "sim_d.run();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_d.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_d.df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the *duration* of `estimate_pi()` CPU use scales pretty much the same as in serial execution, we'll show in the next section how the parallel run improves the application speed dramatically – depending on the number of available CPUs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing Aggregate Measures\n",
    "\n",
    "Now let's build a pipeline in `scikit-learn` to run a \n",
    "[*polynomial regression*](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html)\n",
    "of the *duration* measures.\n",
    "In other words, we'll perform some \"curve fitting\" to calculate the speedup of using *remote functions*.\n",
    "\n",
    "The `degree` parameter is set to the value `1`, which means we'll look at a linear regression.\n",
    "This pipeline allows you to evaluate fitting to higher order polynomials, by changing this parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "DEGREE = 1\n",
    "polyreg = make_pipeline(PolynomialFeatures(DEGREE), LinearRegression())\n",
    "\n",
    "X = sim_s.df.iloc[:, 0].values.reshape(-1, 1)\n",
    "Y_s = sim_s.df.iloc[:, 1].values.reshape(-1, 1)\n",
    "Y_d = sim_d.df.iloc[:, 1].values.reshape(-1, 1)\n",
    "\n",
    "polyreg.fit(X, Y_s)\n",
    "Y_s_pred = polyreg.predict(X)\n",
    "\n",
    "polyreg.fit(X, Y_d)\n",
    "Y_d_pred = polyreg.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X, Y_s, color=\"lightgray\")\n",
    "plt.scatter(X, Y_d, color=\"gray\")\n",
    "\n",
    "plt.plot(X, Y_s_pred, color=\"green\", linewidth=1, label=\"serial\")\n",
    "plt.plot(X, Y_d_pred, color=\"red\", linewidth=0.5, label=\"distrib\")\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distributed processing shows less time required for the tasks in aggregate.\n",
    "Let's calculate how much speedup, based on a ratio of the *slope* values for the two fitted regressions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "m = (Y_s_pred / Y_d_pred)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot these ratios, and take their median value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(m, color=\"red\", linewidth=1, label=\"speedup {:.2f}\".format(np.median(m)))\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depending on the number of available CPUs, you should see an asymptotically linear increase in performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ray.shutdown()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Try running the `Sim` simulation again in its parallelized mode, with a much larger number of samples.\n",
    "While that's running, check the Ray dashboard to see how processing for the *remote tasks* gets distributed across the available CPUs."
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}