course2/01_resource_usage_measure_and_profiling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook 1 - Measuring resource usage and profiling of python code\n",
    "----------------------------------------------\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Content <a id='toc'></a>\n",
    "\n",
    "0. **[Introduction - meet the code](#0)**  \n",
    "   <br>\n",
    "   \n",
    "1. **[Timing](#2)**  \n",
    "    1.1 [Timing a single object](#1.1)  \n",
    "    1.2 [Timing a set of lines](#1.2)  \n",
    "    1.3 [Profiling](#1.3)  \n",
    "    <br>\n",
    "    \n",
    "2. **[Measuring RAM usage](#3)**  \n",
    "    2.1 [Line-by-line memory](#2.1)  \n",
    "    2.2 [Time-based memory usage](#2.2)  \n",
    "    <br>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "## Introduction - meet the code <a id='0'></a>\n",
    "---------------------------------------\n",
    "\n",
    "The first step of any code optimization process should be measuring what your code is doing, in order to pinpoint where your effort should be focused."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the IPython `autoreload` extension.\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "In this notebook we will mostly focus on a simple function which computes pairwise distances between a set of vectors (rows of a 2-dimensional matrix), a very classical operation present in many data analysis methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pairwise_distance(x):\n",
    "    \"\"\"Compute pairwise Euclidean distances between a rows of the input matrix `x`.\n",
    "    \n",
    "    Arguments:\n",
    "        x: a 2-dimensional numpy array (matrix) of numbers, or any nested sequence\n",
    "           of sequences that all have the same length.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Create a square matrix whose size is the number of vectors (rows) in the\n",
    "    # input matrix.\n",
    "    # This matrix will be used to store the Euclidean distances between each pair\n",
    "    # of vectors (rows of the matrix).\n",
    "    num_vectors = len(x)          # Number of rows of the input matrix.\n",
    "    num_measurements = len(x[0])  # Number of columns of the input matrix.\n",
    "    distance_matrix = [[0]*num_vectors for x in range(num_vectors)]\n",
    "    \n",
    "    # Loop over all possible combinations of vectors (rows of the input matrix).\n",
    "    for i in range(num_vectors):\n",
    "        for j in range(num_vectors):\n",
    "            \n",
    "            # Compute the squared distances between all elements of vectors (rows) \"i\" and \"j\".\n",
    "            d = []\n",
    "            for k in range(num_measurements):\n",
    "                d.append((x[i][k] - x[j][k]) ** 2)\n",
    "            \n",
    "            # Euclidean distance between vectors (rows) \"i\" and \"j\". This is\n",
    "            # computed as the square root of the sum of squared distances between\n",
    "            # individual elements of the vectors.\n",
    "            distance_matrix[i][j] = sum(d) ** 0.5\n",
    "\n",
    "    return(distance_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "As a test dataset, let's generate a 200 x 100 matrix filled with random data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "num_vector = 200    # Number of rows.\n",
    "num_measures = 100  # Number of columns.\n",
    "\n",
    "data = np.random.uniform(size=(num_vector, num_measures))\n",
    "print(data)\n",
    "print(data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Let's also generate some random nucleotide data that we will need for another example, later on.\n",
    "We store these nucleotide sequences in files that we create on-disk: `data/large_file.fas` and `data/medium_file.fas`.  \n",
    "Notes:\n",
    "* The cell can take 1-2 minutes to run.\n",
    "* **`%%time`** is an\n",
    "  **[IPython \"magic\" command](https://ipython.readthedocs.io/en/stable/interactive/magics.html#)**\n",
    "  that measures and displays the time needed to run a call.  \n",
    "  More specifically, it will give the following values:\n",
    "  * CPU time: time during which the CPU was in use.\n",
    "  * Wall time: actual time it took to run the cell (includes things such as time spent reading/writing to disk)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "# Create a \"data\" directory, if not already present.\n",
    "if not os.path.isdir(\"data\"):\n",
    "    os.mkdir(\"data\")\n",
    "\n",
    "# 1'000'000 random sequences of 500 nucleotides.\n",
    "with open(\"data/large_file.fas\", \"w\") as OUT:\n",
    "    for i in range(1000000):\n",
    "        print(\">seq{}\".format(i), file=OUT)\n",
    "        s = \"\".join(np.random.choice(list(\"ATGC\"), size=500))\n",
    "        print(s, file=OUT)\n",
    "\n",
    "# 500 random sequences of 500 nucleotides.\n",
    "with open(\"data/medium_file.fas\", \"w\") as OUT :\n",
    "    for i in range(500):\n",
    "        print(\">seq{}\".format(i), file=OUT)\n",
    "        s = \"\".join(np.random.choice(list(\"ATGC\"), size=500))\n",
    "        print(s, file=OUT)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "# 1. Measuring time usage <a id='2'></a>\n",
    "\n",
    "\n",
    "## 1.1 timing a single object <a id='1.1'></a>\n",
    "\n",
    "On you terminal, you may measure up the time taken by a python script execution using:\n",
    "\n",
    "* On **Linux** and **Mac OS**.\n",
    "```sh\n",
    "time python my_script.py\n",
    "```\n",
    "\n",
    "* On **Windows**, the following should work when using the Windows PowerShell (but not the old windows shell).\n",
    "```sh\n",
    "Measure-Command {start-process command-to-benchmark -Wait}\n",
    "```\n",
    "\n",
    "<br>\n",
    "\n",
    "Or, by using **[IPython \"magic\" commands](https://ipython.readthedocs.io/en/stable/interactive/magics.html#)**:\n",
    "* **`%time`**: measures the time to run a line of code.\n",
    "* **`%%time`**: measures the time to run an entire cell.\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "> **Warning:** IPython magic commands only work in Jupyter-Notebook (or other IPython compatible shells).\n",
    "  They will not work if you try to run them in a classic python interpreter.\n",
    "\n",
    "<div>\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Examples:**\n",
    "* Measuring execution time of **a single line**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time D = pairwise_distance(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Measuring execution time of **an entire Jupyter cell**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Applies on the whole cell.\n",
    "\n",
    "D = pairwise_distance(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "That is all nice, but **there is always a bit a variability between runs**, which becomes very apparent on lines (or cells) that have a smaller execution time.  \n",
    "* To illustrate this, let's try to run the `pairwise_distance()` on a small 100 x 10 matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "small_data = np.random.uniform(size=(100, 10))\n",
    "\n",
    "%time D = pairwise_distance(small_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "To solve this problem, we have the **`timeit`** module:\n",
    "\n",
    "```sh\n",
    "python -m timeit myScript.py\n",
    "```\n",
    "\n",
    "Or the equivalent magic command:\n",
    "* **`%timeit`**: benchmarks the time needed to execute a line of code. `%timeit` carries-out multiple\n",
    "  replicates of the line to benchmark and averages the results.\n",
    "* By default, `timeit` will perform a number of replicates that it deems *optimal*, but\n",
    "  the number of **repeats** and **loops** can also be manually specified:\n",
    "  * **`-n`:** number of loops, i.e. the number of times the code is run in a given repeat/run.\n",
    "    The benchmark time is then averaged over all `n` loops.\n",
    "  * **`-r`:** number of runs/repeats. The benchmark is repeated `r` times and only the best run/repeat is kept.\n",
    "    This allows to get rid of runs that might have been slow for some reason external to the code (e.g.\n",
    "    the machine was busy doing something else at the same time).\n",
    "  * Example: `%timeit -n 2 -r 10`, 2 loops and 10 repeats.  \n",
    "\n",
    "<br>\n",
    "\n",
    "* **`%%timeit`** is essentially the same as `%timeit`, but applies to an entire cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit D = pairwise_distance(small_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Using the **`-r`** and **`-n`** options, the number of runs and loops to perform can be manually specified.\n",
    "* `%timeit` will perform `-r` **repeats** of `-n` loops each. The repeat with the best execution time is kept."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit -n 2 -r 10 D = pairwise_distance(small_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "    \n",
    "**Question:** why does `%timeit` take the \"best\" out of `r` loops (rather than the average) ?\n",
    "\n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Ok, already this is nice, we will definitely be using this to compare together different implementations.  \n",
    "For example, I have rewritten the function using **`numpy`**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here is the old version for comparison:\n",
    "def pairwise_distance(X):\n",
    "\n",
    "    num_vectors = len(X)\n",
    "    num_measurements = len(X[0])\n",
    "    D = [[0]*num_vectors for x in range(num_vectors)]\n",
    "    \n",
    "    for i in range(num_vectors):\n",
    "        for j in range(num_vectors):\n",
    "            d = []\n",
    "            for k in range(num_measurements):\n",
    "                d.append((X[i][k] - X[j][k]) ** 2)\n",
    "            \n",
    "            D[i][j] = sum(d) ** 0.5\n",
    "    return(D)\n",
    "\n",
    "def pairwise_distance_numpy(X):\n",
    "\n",
    "    num_vectors = X.shape[0]\n",
    "    num_measurements = X.shape[1] \n",
    "    D = np.empty((num_vectors, num_vectors), dtype=np.float64)\n",
    "    \n",
    "    for i in range(num_vectors):\n",
    "        for j in range(num_vectors):\n",
    "            d = np.square(np.subtract(X[i], X[j]))\n",
    "            D[i, j] = np.sqrt(np.sum(d))\n",
    "    return(D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "**Let's run a benchmark** on these 2 implementations using a **100 x 100 matrix of random values**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.random.uniform(size=(100, 100))\n",
    "\n",
    "print(\"native python:\")\n",
    "%timeit -n 5 -r 3 D = pairwise_distance(data)\n",
    "\n",
    "print(\"numpy:\")\n",
    "%timeit -n 5 -r 3 D = pairwise_distance_numpy(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nice! \n",
    "\n",
    "**Trick**: adding the option **`-o`** to the `%timeit` command allows to save the outputs to a variable - `timeit_res` in the example below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeit_res = %timeit -n 5 -r 3 -o D = pairwise_distance_numpy(data)\n",
    "print(\"average:\",timeit_res.average , \"standard-dev\", timeit_res.stdev )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is a neat trick to help us investigate how execution time evolves with the data size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "times = []\n",
    "N = []\n",
    "\n",
    "num_measures = 10\n",
    "num_vector_list = range(10, 200, 10)\n",
    "\n",
    "for num_vector in num_vector_list:\n",
    "    data = np.random.uniform(size=(num_vector,num_measures))\n",
    "    \n",
    "    # Here we also use the \"-q\" option to suppress the text output.\n",
    "    timeit_res = %timeit -n 2 -r 7 -o -q D = pairwise_distance_numpy(data)\n",
    "    times.extend(timeit_res.timings)\n",
    "    N.extend( [num_vector]*len(timeit_res.timings) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "fig = plt.figure(figsize = (8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "sns.lineplot(x = N, y = times, ax = ax)\n",
    "\n",
    "ax.set_aspect(1.0/ax.get_data_ratio(), adjustable=\"box\")\n",
    "\n",
    "plt.xlabel(\"number of vectors\", size=18)\n",
    "plt.ylabel(\"time(s)\", size=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Question:** What do you think about the shape of this curve? Was this expected?\n",
    "\n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## 1.2 timing a set of lines<a id='1.2'></a>\n",
    "\n",
    "\n",
    "Sometimes you will want to measure the execution time of a particular step in you code, but it may not be easy to isolate it for a usage with `time` or `timeit`.\n",
    "\n",
    "For instance, consider the following code which reads a [FASTA file](https://en.wikipedia.org/wiki/FASTA_format) and computes the [GC content](https://en.wikipedia.org/wiki/GC-content) of each sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "GC=[]\n",
    "\n",
    "with open(\"data/large_file.fas\", \"r\") as IN :\n",
    "    for l in IN:\n",
    "        if not l.startswith(\">\"):\n",
    "            GC.append((l.count(\"C\") + l.count(\"G\")) * 100 / len(l.strip()))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "The codes takes about 4.5s to complete, but how much of this is file reading, and how much is GC content computation ?\n",
    "* Here we do not really want to re-write the code to have neatly separated steps to apply `%time` on.\n",
    "  Additionally, we may not be able store all the sequences in memory.\n",
    "* In these situation, the **`time()`** function from module **`time`** (so, **`time.time()`**) is quite useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It returns the time (in second) since the Epoch, which is 00:00:00 on 1 January 1970.\n",
    "\n",
    "> **Epoch**, as defined in [wikipedia](https://en.wikipedia.org/wiki/Epoch_(computing)):  \n",
    "> In computing, an epoch is a date and time from which a computer measures system time. Most computer systems determine time as a number representing the seconds removed from particular arbitrary date and time. For instance, Unix and POSIX measure time as the number of seconds that have passed since Thursday 1 January 1970 00:00:00 UT, a point in time known as the Unix epoch.  \n",
    "Windows NT systems, up to and including Windows 11 and Windows Server 2022, measure time as the number of 100-nanosecond intervals that have passed since 1 January 1601 00:00:00 UTC, making that point in time the epoch for those systems.\n",
    "\n",
    "Applied to our code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GC_time = 0\n",
    "\n",
    "start = time.time()  # Record time at the start of the block to benchmark.\n",
    "with open(\"data/large_file.fas\", \"r\") as IN:\n",
    "    GC = []\n",
    "    for l in IN:\n",
    "        if not l.startswith(\">\"):\n",
    "            t1 = time.time()   # Record time at the start of the block to benchmark.\n",
    "            GC.append((l.count(\"C\") + l.count(\"G\")) * 100 / len(l.strip()))\n",
    "            t2 = time.time()   # Record time at the end of the block to benchmark.\n",
    "            GC_time += t2 - t1 # Compute elapsed time.\n",
    "\n",
    "stop = time.time()         # Record time at the end of the block to benchmark.\n",
    "total_time = stop - start  # Compute elapsed time.\n",
    "\n",
    "print(\"Total: {:.3f}s\".format(total_time))\n",
    "print(\"  GC%: {:.3f}s\".format(GC_time))\n",
    "print(\" read: {:.3f}s\".format(total_time - GC_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br> \n",
    "\n",
    "[back to ToC](#toc)\n",
    "\n",
    "## 1.3 profiling<a id='1.3'></a>\n",
    "\n",
    "\n",
    "Most of the time your code is more complex than a single function, and before optimizing you first want to see which function you should optimize. That is when **profiling** comes in handy.\n",
    "\n",
    "For instance, consider the following code, which reads sequences from a FASTA file, sort them by GC content, then computes a matrix of pairwise distance between all sequences and finally writes this matrix to a file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def read_fasta(filename):\n",
    "    \"\"\"Reads a FASTA file and returns its sequences as a dictionary.\"\"\"\n",
    "    Dseq={}\n",
    "    curseq = \"\"\n",
    "    cur = \"\"\n",
    "    with open(filename, \"r\") as IN:\n",
    "        for l in IN:\n",
    "            if l.startswith(\">\"):\n",
    "                if cur != \"\":\n",
    "                    Dseq[cur] = curseq\n",
    "                cur = l[1:].strip()\n",
    "                curseq = \"\"\n",
    "            else:\n",
    "                curseq += l.strip()\n",
    "                \n",
    "        if cur != \"\":\n",
    "            Dseq[cur] = curseq\n",
    "            \n",
    "    return Dseq\n",
    "            \n",
    "\n",
    "def computeGC( seq ):\n",
    "    \"\"\"Takes a sequence of nucleotides (str) and compute its\n",
    "    GC percentage (float).\n",
    "    \"\"\"\n",
    "    gc = 0\n",
    "    for i in range(len(seq)):\n",
    "        if seq[i] in \"GC\":\n",
    "            gc += 1\n",
    "    return 100 * gc / len(seq)\n",
    "\n",
    "def computeGC_dict(Dseq):\n",
    "    \"\"\"Takes a dictionary containing sequences as values \n",
    "    and compute a dictionary containing their GC%.\n",
    "    \"\"\"\n",
    "    Dgc = {}\n",
    "    for k in Dseq:\n",
    "        Dgc[k] = computeGC(Dseq[k])\n",
    "    return Dgc\n",
    "\n",
    "def compute_sequence_similarity(seqA  ,seqB):\n",
    "    \"\"\"Compute similarity between 2 sequence as the fraction\n",
    "    of position where they have the same value.\n",
    "    \"\"\"\n",
    "    l = len(seqA)\n",
    "    similar = 0\n",
    "    for i in range(l):\n",
    "        if seqA[i] == seqB[i]:\n",
    "            similar += 1\n",
    "    return similar / l\n",
    "\n",
    "\n",
    "\n",
    "def main_script(input_filename, output_filename):\n",
    "    \n",
    "    # Step 1: read fasta\n",
    "    Dseq = read_fasta(input_filename)\n",
    "    \n",
    "    # Step 2: compute GC%\n",
    "    Dgc = computeGC_dict( Dseq )\n",
    "    \n",
    "    # Step 3: sort by GC%.\n",
    "    ordered_seq = sorted(Dgc.keys(), key = lambda x: Dgc[x])\n",
    "    \n",
    "    # Step 4: compute pairwise distance matrix.\n",
    "    sim = np.zeros((len(Dseq), len(Dseq)))\n",
    "    for i, s1 in enumerate(ordered_seq):\n",
    "        for j, s2 in enumerate(ordered_seq):\n",
    "            sim[i, j] = compute_sequence_similarity(Dseq[s1], Dseq[s2])\n",
    "\n",
    "    # Step 5: write the matrix.\n",
    "    with open(output_filename, \"w\") as OUT:\n",
    "        print(\",\".join(ordered_seq), file=OUT)\n",
    "        for i in range(len(ordered_seq)):\n",
    "            print( *(sim[i]), sep=\",\", file=OUT)\n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if you have the eye for it, it looks like most of these function could be rewritten to be faster.\n",
    "\n",
    "For instance, the function to compute the GC%:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def computeGC_better( seq ):\n",
    "    \"\"\"takes a sequence (str) and compute it GC% (float).\"\"\"\n",
    "    return 100 * (seq.count(\"G\") + seq.count(\"C\")) / len(seq)\n",
    "\n",
    "seq = \"ATGC\" * 5000\n",
    "%timeit -n 100 -r 10 computeGC(seq)\n",
    "%timeit -n 100 -r 10 computeGC_better(seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is a commendable speedup!  \n",
    "But, considering that coding time is a finite resource, where should we start? Where is out effort better spent?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We recommend using **[cProfile](https://docs.python.org/3/library/profile.html)**.  \n",
    "In the terminal :\n",
    "\n",
    "```sh\n",
    "python -m cProfile -o profile.log -s cumtime myScript.py\n",
    "```\n",
    "\n",
    "will execute the script, and profile time usage of every functions\n",
    "* `-o` : output file for the profiling log.\n",
    "* `-s cumtime` : to sort by cumulative time spent in a single function.\n",
    "\n",
    "Then to interpret the output log from cProfile, we recommend the\n",
    "**[snakeviz](https://jiffyclub.github.io/snakeviz)** library.\n",
    "\n",
    "Otherwise, in Jupyter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# IMPORTANT: keep single quotes around the file names/path arguments, otherwise\n",
    "#            it won't run on Windows.\n",
    "\n",
    "%prun -l 30 -s cumtime  main_script('data/medium_file.fas', 'test.out')\n",
    "# The %prun magic command activate profiling\n",
    "#  -l 30: limits the report to 30 lines\n",
    "#  -s cumtime: sort by decreasing cumtime\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The columns correspond to:\n",
    " * `ncalls`: for the **number of calls**.\n",
    " * `tottime`: for the **total time** spent in the given function (excluding time spent in calls to sub-functions).\n",
    " * `percall`: is the quotient of tottime divided by ncalls.\n",
    " * `cumtime`: is the **cumulative time spent in this and all subfunctions** (from invocation till exit).\n",
    "   This figure is accurate even for recursive functions.\n",
    " * `percall`: is the quotient of cumtime divided by primitive calls.\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "### Micro-Exercise 1\n",
    "\n",
    "* Look at the profile. Where should optimization efforts go first?\n",
    "  What would be the effect of using our better implementation of the GC% computing function?\n",
    "\n",
    "<div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "We now have the tools to help us diagnose which part of our code takes the most time. But, before we move on to optimization, let's see what would happen if we launched our code on a larger dataset:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time main_script('data/large_file.fas', 'test.out')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the problem is not the execution time, but the RAM (random access memory) usage. \n",
    "\n",
    "While time is a somewhat flexible constraint (it is always possible to wait a bit longer), memory is a hard limit: you either make your code less memory-hungry, or you move to another computer...  \n",
    "Let's focus on memory now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "\n",
    "# 2. Measuring RAM usage <a id='3'></a>\n",
    "\n",
    "## 2.1 Line-by-line memory <a id='2.1'></a>\n",
    "\n",
    "To measure the memory imprint of your code - a nice tool is **[memory-profiler](https://pypi.org/project/memory-profiler)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Install it with :\n",
    "\n",
    "```python\n",
    "!pip install --user memory_profiler\n",
    "```\n",
    "\n",
    "Basically, in your code you add a **decorator** to the function of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from memory_profiler import profile\n",
    "\n",
    "@profile\n",
    "def my_func():\n",
    "    a = [1] * (10 ** 6)\n",
    "    b = [2] * (2 * 10 ** 7)\n",
    "    time.sleep(1)\n",
    "    del b\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, you may either do it command-line style  :\n",
    "\n",
    "```\n",
    "python -m memory_profiler example.py\n",
    "```\n",
    "\n",
    "or prefer jupyter-magic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%memit _= my_func()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It gives us a result: `peak memory: 210.19 MiB, increment: 152.58 MiB`\n",
    "            \n",
    "But it does not like that our code is in the same notebook. Let's have it in another script :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import my_func\n",
    "\n",
    "_ = my_func()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, that is quite nice: we have a run down of the RAM usage, line-by-line.\n",
    "\n",
    "Let's see how that works for our `pairwise_distance` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile tmp.py\n",
    "\n",
    "import numpy as np\n",
    "from memory_profiler import profile\n",
    "\n",
    "# The memory increments are fairly small here, so we set the precision higher.\n",
    "@profile(precision=3)\n",
    "def pairwise_distance_profile(X):\n",
    "\n",
    "    num_vectors = len(X)\n",
    "    num_measurements = len(X[0])\n",
    "    D = [[0]*num_vectors for x in range(num_vectors)]\n",
    "    \n",
    "    for i in range(num_vectors):\n",
    "        for j in range(num_vectors):\n",
    "            d = []\n",
    "            for k in range(num_measurements):\n",
    "                d.append((X[i][k] - X[j][k]) ** 2)\n",
    "            \n",
    "            D[i][j] = sum(d) ** 0.5\n",
    "    return(D)\n",
    "\n",
    "\n",
    "# The precision parameter does not work in jupyter notebooks :-( \n",
    "# so I integrate the test to the script.\n",
    "num_vector = 100\n",
    "num_measures = 10\n",
    "\n",
    "data = np.random.uniform(size=(num_vector,num_measures))\n",
    "_ = pairwise_distance_profile(data) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "!python tmp.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> NB : the precision parameter does not work in jupyter notebooks, so here the tests are directly integrated to the script.\n",
    "\n",
    "That is super neat, but we have to note that this also took a while : tracking all of this memory creates a lot of overhead.\n",
    "\n",
    "Compare with the version without `@profile`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_vector = 100\n",
    "num_measures = 10\n",
    "data = np.random.uniform(size=(num_vector,num_measures))\n",
    "\n",
    "%time _=pairwise_distance(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is a big slow down: x100. \n",
    "For scripts with longer execution times it can get fairly prohibitive to profile the memory in such a fine way.\n",
    "\n",
    "Let's explore an alternative with less overhead."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## 2.2 time-based memory usage<a id='2.2'></a>\n",
    "\n",
    "[**mprof**](https://pypi.org/project/memory-profiler) is an executable which let's you monitor any script memory usage over-time (you can install mprof with `pip` or `conda`).\n",
    "\n",
    "It comes packaged with `memory_profiler` and allows some nice integration: \n",
    "it uses `@profile` to annotate its report and plot.\n",
    "\n",
    "**However,** this only works if you **don't import memory_profiler in the script**... otherwise it defaults back the line-by-line profiling.\n",
    "\n",
    "On the previous script:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile tmp.py\n",
    "\n",
    "import numpy as np\n",
    "# from memory_profiler import profile --> do not import this\n",
    "\n",
    "@profile(precision=3)\n",
    "def pairwise_distance_profile(X):\n",
    "\n",
    "    num_vectors = len(X)\n",
    "    num_measurements = len(X[0])\n",
    "    D = [[0]*num_vectors for x in range(num_vectors)]\n",
    "    \n",
    "    for i in range(num_vectors):\n",
    "        for j in range(num_vectors):\n",
    "            d = []\n",
    "            for k in range(num_measurements):\n",
    "                d.append( ( X[i][k] - X[j][k] )**2 )\n",
    "            \n",
    "            D[i][j] = sum(d) **0.5\n",
    "    return(D)\n",
    "\n",
    "\n",
    "num_vector = 500\n",
    "num_measures = 10\n",
    "\n",
    "data = np.random.uniform(size=(num_vector,num_measures))\n",
    "\n",
    "_=pairwise_distance_profile(data)\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we run this script with mprof:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!time mprof run tmp.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This generated a file with a title looking like `mprofile_20220711084326.dat`.\n",
    "\n",
    "We can generate a plot of this profile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mprof plot -o tmp.png"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we load this image in the jupyter notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"tmp.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we clearly see the initial memory loading (linked to the importation of libraries and generation of the data), and the function of interest is clearly marked. \n",
    "We can see it leads to an increase of RAM usage of about 8MiB.\n",
    "\n",
    "Also note that the overhead is much lower:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_vector = 500\n",
    "num_measures = 10\n",
    "data = np.random.uniform(size=(num_vector,num_measures))\n",
    "\n",
    "%time _=pairwise_distance(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1.5s, vs. 2.3s with `mprof`: much more reasonable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* `mprof` is very useful to explore and pinpoint memory spike, especially since it\n",
    "  **works with all executables** and not only python scripts.\n",
    "* You can increase the granularity of the report using the `--interval` parameter (default: 0.1s).\n",
    "* `mprof` also has a mode designed to monitor executables using multiprocessing.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## 2.3 getting the size of a single object <a id='2.3'></a>\n",
    "\n",
    "Last but not least, when you know your code enough you can often point to the precise object who represents the majority of RAM usage in your code.\n",
    "\n",
    "In the case of the script for distance computation between sequences in a FASTA file, the error message pinpointed the problematic line:\n",
    "\n",
    "```python\n",
    "---------------------------------------------------------------------------\n",
    "MemoryError                               Traceback (most recent call last)\n",
    "File <timed eval>:1, in <module>\n",
    "\n",
    "Input In [59], in main_script(input_filename, output_filename)\n",
    "     57 ordered_seq = sorted(Dgc.keys() , key=lambda x: Dgc[x])\n",
    "     59 # step 4 : compute pairwise distance matrix\n",
    "---> 60 sim = np.zeros((len(Dseq), len(Dseq)))\n",
    "     61 for i, s1 in enumerate(ordered_seq):\n",
    "     62     for j, s2 in enumerate(ordered_seq):\n",
    "\n",
    "MemoryError: Unable to allocate 7.28 TiB for an array with shape (1000000, 1000000) and data type float64\n",
    "```\n",
    "\n",
    "So here, no need to use advanced tools : the problem is this square matrix.\n",
    "\n",
    "We can investigate the memory needed by an object in memory with `sys.getsizeof`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "a = 0.5 # a simple float\n",
    "print(\"size of a float:\" , sys.getsizeof(a))\n",
    "\n",
    "b =\"abcdef\" # a simple string\n",
    "print(\"size of a str:\" , sys.getsizeof(b))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reported size is in bytes. To get kib: divide by 1024. \n",
    "\n",
    "To get MiB, divide by 1024*1024"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "c = np.zeros((N,N)) #NxN matrix \n",
    "print(\"size of a {}x{} matrix: {:.2f} MiB\".format(N,N, sys.getsizeof(c)/(1024*1024)) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A little *caveat* though. Consider the following code: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "a = 0.5 # a simple float\n",
    "print(\"size of a float:\" , sys.getsizeof(a))\n",
    "\n",
    "b = [np.random.random() for i in range(10)]\n",
    "print(\"size of a list of 10 floats:\" , sys.getsizeof(b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See anything strange? \n",
    "\n",
    "If a single `float` is 24bytes, then how can a list of 10 floats be less than 10*24=240 bytes ?\n",
    "\n",
    "<br>\n",
    "\n",
    "This is because `getsizeof` only account for direct memory does not go follow references to objects. \n",
    "In practice, that means it struggles with containers.\n",
    "\n",
    "The official documentation point to this function if you want to get the total size of an object, including everything it contains or refers to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sys import getsizeof, stderr\n",
    "from itertools import chain\n",
    "from collections import deque\n",
    "try:\n",
    "    from reprlib import repr\n",
    "except ImportError:\n",
    "    pass\n",
    "\n",
    "def total_size(o, handlers={}, verbose=False):\n",
    "    \"\"\"Returns the approximate memory footprint an object and all of its contents.\n",
    "\n",
    "    Automatically finds the contents of the following builtin containers and\n",
    "    their subclasses:  tuple, list, deque, dict, set and frozenset.\n",
    "    To search other containers, add handlers to iterate over their contents:\n",
    "\n",
    "        handlers = {SomeContainerClass: iter,\n",
    "                    OtherContainerClass: OtherContainerClass.get_elements}\n",
    "    \"\"\"\n",
    "    dict_handler = lambda d: chain.from_iterable(d.items())\n",
    "    all_handlers = {tuple: iter,\n",
    "                    list: iter,\n",
    "                    deque: iter,\n",
    "                    dict: dict_handler,\n",
    "                    set: iter,\n",
    "                    frozenset: iter,\n",
    "                   }\n",
    "    all_handlers.update(handlers)     # user handlers take precedence\n",
    "    seen = set()                      # track which object id's have already been seen\n",
    "    default_size = getsizeof(0)       # estimate sizeof object without __sizeof__\n",
    "\n",
    "    def sizeof(o):\n",
    "        if id(o) in seen:       # do not double count the same object\n",
    "            return 0\n",
    "        seen.add(id(o))\n",
    "        s = getsizeof(o, default_size)\n",
    "\n",
    "        if verbose:\n",
    "            print(s, type(o), repr(o), file=stderr)\n",
    "\n",
    "        for typ, handler in all_handlers.items():\n",
    "            if isinstance(o, typ):\n",
    "                s += sum(map(sizeof, handler(o)))\n",
    "                break\n",
    "        return s\n",
    "\n",
    "    return sizeof(o)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"total size of a list of 10 floats:\" , total_size(b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, `getsizeof` returned 184, which is the size of the list, add to this the 10 float : 10*24, and you get\n",
    "$184+24*10=424$.\n",
    "It works!\n",
    "\n",
    "Note that for numpy array this does not change anything (because numpy arrays do not access their data by reference)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"array getsizeof :\", sys.getsizeof(c))\n",
    "print(\"array total_size:\", total_size(c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "### Micro-Exercise 2\n",
    "\n",
    "* Find out which is the largest square matrix your RAM can reasonably accommodate.\n",
    "    (NB: you can do this by testing it, or deducing it from the size of smaller objects)\n",
    "* **Additional task (if you have time):** how could we modify the `main_script` to make it less memory hungry?\n",
    "\n",
    "<div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "# Additional material\n",
    "------------------------------\n",
    "\n",
    "This section will not be covered in the course, but you can read it on your own if you are interested.\n",
    "\n",
    "## Annex 1 - Kmeans implementation profiling\n",
    "\n",
    "Imagine you have a script, implementing a Kmeans algorithm. \n",
    "Here are the functions which look like the best candidates for optimization:\n",
    "* `computeDistanceToCentroid`: compute the distance between a point and a centroid.\n",
    "* `computeNearestCentroid`: compute which centroid is the closest to each point (actually\n",
    "  calls `computeDistanceToCentroid`, but also possess some other potentially costly computations).\n",
    "* `computeCentroids`: computes the position of the centroid of a points with a given cluster assignment.\n",
    "\n",
    "Are they really the best candidate ? which one should we go for first ?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generating some random data.\n",
    "\n",
    "def generateCluster(n, means, sds):\n",
    "    P = np.random.randn(len(means), n)\n",
    "    for i in range(len(means)):\n",
    "        P[i,] = P[i,] * sds[i] + means[i]\n",
    "    \n",
    "    return P\n",
    "\n",
    "\n",
    "clusterSizes = [4000, 2000, 4000, 4000, 2000]\n",
    "clusterMeans = [ [0, -2], [3, 3], [-1, 3], [-5, 0], [5,-1]]\n",
    "clusterSDs = [[0.5,1], [1,0.5], [0.5,0.5], [2,1], [1,1]]\n",
    "C = []\n",
    "A = []\n",
    "for i in range(len(clusterSizes)):\n",
    "    C.append(generateCluster(clusterSizes[i], clusterMeans[i], clusterSDs[i]))\n",
    "    A += [i]*clusterSizes[i]\n",
    "Points = np.concatenate(C, axis=1)\n",
    "realAssignment = np.array(A)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Run the memory profiling with `%pruns` IPython magic command:\n",
    "* `-l 30`: limits the report to 30 lines.\n",
    "* `-s cumtime`: sort output lines by decreasing cumtime (cumulative time)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Kmeans import Kmeans\n",
    "\n",
    "# Performing Kmeans.\n",
    "k = 5\n",
    "%prun -l 30 -s cumtime  kmeanAssignment = Kmeans(Points, k, maxNbRounds=1000)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Environment (conda_py311)",
   "language": "python",
   "name": "conda_py311"
  },
  "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}