lectures/08_dynamic-programming-vectorization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"dsci512_header.png\" width=\"600\">\n",
    "\n",
    "# Lecture 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Outline:\n",
    "\n",
    "- Caching & memoization (15 min)\n",
    "- Dynamic programming intro (15 min)\n",
    "- Break (5 min)\n",
    "- Dynamic programming: lab 4 (10 min)\n",
    "- Dynamic programming: backtracking (10 min)\n",
    "- Code vectorization (20 min)\n",
    "- Course wrap-up (5 min)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning objectives\n",
    "\n",
    "- Identify cases of repeated/redundant computation and address them with caching/memoization.\n",
    "- Describe, at a high level, the basic idea behind dynamic programming.\n",
    "- Identify problems when dynamic programming is applicable.\n",
    "- Identify situations where vectorization is possible.\n",
    "- Convert un-vectorized code to vectorized code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import functools\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import joblib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn.metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Caching & Memoization (5 min)\n",
    "\n",
    "- Consider the Fibonacci sequence $0,1,1,2,3,5,8,13,\\ldots$.\n",
    "- The first two numbers are $F_0=0$ and $F_1=1$.\n",
    "- For the rest, $F_{n}=F_{n-1}+F_{n-2}$; that is, each number is the sum of the previous two numbers.\n",
    "- How to compute this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fib(n):\n",
    "    f = np.zeros(n+1, dtype=int)\n",
    "    f[1] = 1\n",
    "    \n",
    "    for i in range(2,n+1):\n",
    "        f[i] = f[i-1] + f[i-2]\n",
    "    \n",
    "    return f[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fib(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Easy. Well, I guess we could also do this recursively?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fib_rec(n):\n",
    "    if n == 0 or n == 1:\n",
    "        return n\n",
    "    \n",
    "    return fib_rec(n-1) + fib_rec(n-2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the code looks even cleaner! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fib_rec(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But wait..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30.3 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 fib(35)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.5 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 fib_rec(35)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Hmm... what's the problem?\n",
    "\n",
    "<br><br><br><br><br><br>\n",
    "\n",
    "- We can see A LOT of repeated computation. \n",
    "- For example, `fib(33)` is called by `fib(34)` and `fib(35)`, and it gets much worse deeper down the tree.\n",
    "  - Indeed, `fib(1)` is called millions of times..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "node A --- GIANT GRAPH --- node B --- node C\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a common theme! Other examples:\n",
    "\n",
    "- Lab 3 `largest_distance` had repeated computations of distances between vertices.\n",
    "- Factorio optional question from lab 2 recomputes the raw ingredients of an item over and over along different parts of the tree.\n",
    "- Recursive seam carving code in lab 4 has the same problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- How to solve this?\n",
    "- The general term **caching** refers to storing things for later.\n",
    "- The specific term **memoiziation** refers to storing the results of function calls for later. See [here](https://stackoverflow.com/questions/6469437/what-is-the-difference-between-caching-and-memoization).\n",
    "- Confusingly, memoization in Python is called cache; it's OK."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "@functools.lru_cache(maxsize=None)\n",
    "def fib_rec_cache(n):\n",
    "    if n == 0 or n == 1:\n",
    "        return n\n",
    "    \n",
    "    return fib_rec_cache(n-1) + fib_rec_cache(n-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14.4 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 fib_rec_cache(35)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Amazing! What happened here?\n",
    "- We have learned about [Python decorators](https://realpython.com/primer-on-python-decorators/) in DSCI 511. We use a particular type of decorator here to remember result of function calls:\n",
    "  - For example, if `fib_rec(30)` has already been computed, the result ($832040$) is saved and used if `fib_rec(30`) is ever called again.\n",
    "  - This, we only call `fib_rec(n)` once for each `n`.\n",
    "- This is a useful trick/concept for speeding things up. \n",
    "  - It's a tradeoff of more memory usage for faster running time. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other approaches:\n",
    "\n",
    "- `functools.lru_cache` stored the cache in memory.\n",
    "- We can use `joblib` to store the cache in a file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "memory = joblib.Memory(\"/tmp\", verbose=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@memory.cache\n",
    "def fib_rec_cache2(n):\n",
    "    if n == 0 or n == 1:\n",
    "        return n\n",
    "    \n",
    "    return fib_rec_cache2(n-1) + fib_rec_cache2(n-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "80.5 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 fib_rec_cache2(35)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dynamic programming: seam carving (25 min)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Looking at the lab \n",
    "\n",
    "For this part of the class we'll open up lab 4 and take a look at it.\n",
    "_In the lab, let's go through the intro and the recursive implementation._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "top of image\n",
    "  \n",
    "       x \n",
    "\n",
    "bottom of image\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Why is the code so slow?\n",
    "\n",
    "- This is the same issue as we saw earlier with the Fibonacci numbers.\n",
    "- We are repeating computations for subproblems.\n",
    "- We won't do it now, but I encourage you to draw out the recursion and see for yourself!\n",
    "- We could use the memoization trick again, but:\n",
    "  - We don't have a clear idea of how much memory it will take - will it be too much?\n",
    "  - We don't have a clear idea of how long it will run.\n",
    "  - It's a bit heavy handed and we cannot clearly see what the code is doing.\n",
    "  - It would be better to just directly write good code for this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Enter dynamic programming \n",
    "\n",
    "- Dynamic programming is an algorithm for solving certain optimization problems.\n",
    "  - It only works for some problems, but when it works it is FAST.\n",
    "  - This is just like linear programming (last lecture), except for a different set of problems.\n",
    "  - It is closely related to memoization; see [here](https://stackoverflow.com/questions/6184869/what-is-the-difference-between-memoization-and-dynamic-programming) for a discussion of memoization vs. dynamic programming. \n",
    "- Dynamic programming has a lot of applications:\n",
    "  - DNA sequence alignment\n",
    "  - Text hyphenation\n",
    "  - Running certain machine learning models (e.g. hidden Markov models)\n",
    "  - Finding the differences between two files (bonus - below!)\n",
    "  - Image resizing (this week's lab!)\n",
    "  - and more...\n",
    "- What does the name \"Dynamic Programming\" mean?\n",
    "  - \"Programming\" is often used for optimization, for historical reasons.\n",
    "  - There is no significance to \"Dynamic\", it just sounds good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Break (5 min)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### And now back to the lab\n",
    "\n",
    "- Let's look at the code for Exercise 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_vertical_seam(energy):\n",
    "    nrows, ncols = energy.shape\n",
    "\n",
    "    # Cumulative Minimum Energy\n",
    "    CME = np.zeros((nrows, ncols+2))\n",
    "    CME[:,0] = np.inf\n",
    "    CME[:,-1] = np.inf\n",
    "    CME[:,1:-1] = energy\n",
    "\n",
    "    for i in range(1,nrows):\n",
    "        for j in range(1,ncols+1):\n",
    "            CME[i,j] += min(CME[i-1,j-1],CME[i-1,j],CME[i-1,j+1])\n",
    "\n",
    "    # create the seam array\n",
    "    seam = np.zeros(nrows, dtype=int)\n",
    "    seam[-1] = np.argmin(CME[-1,:])\n",
    "\n",
    "    # track the path backwards\n",
    "    for i in range(nrows-2,-1,-1):\n",
    "        # min_index is 0, 1, or 2. I subtract 1 to give the offset from\n",
    "        # seam(i+1), namely -1, 0, or 1. Then I add this to the old value.\n",
    "        delta = np.argmin(CME[i, seam[i+1]-1:seam[i+1]+2]) - 1\n",
    "        seam[i] =  seam[i+1] + delta\n",
    "\n",
    "    return seam - 1 \n",
    "    # Above: -1 because the indices are all off by 1, due to padding of CME\n",
    "    #        This has nothing to do with the -1 in defining \"delta'\".\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9, 4, 9, 9, 1],\n",
       "       [6, 4, 6, 8, 7],\n",
       "       [1, 5, 8, 9, 2],\n",
       "       [7, 3, 3, 2, 4],\n",
       "       [6, 9, 2, 9, 8]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(3)\n",
    "nrows, ncols = 5,5\n",
    "energy = np.random.randint(1,10, size=(nrows,ncols))\n",
    "energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the lowest energy seam from the top to the bottom?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 0, 1, 2])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seam = find_vertical_seam(energy)\n",
    "seam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Note that the optimal seam does not include the \"1\" at the top!\n",
    "- This illustrates the point that a seam has to be globally optimal, not just locally!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- How does this work?\n",
    "- We need to compute the cumulative minimum energy from the top to a given location.\n",
    "- For the top row, it's just the energy values themselves. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[inf,  9.,  4.,  9.,  9.,  1., inf],\n",
       "       [inf,  6.,  4.,  6.,  8.,  7., inf],\n",
       "       [inf,  1.,  5.,  8.,  9.,  2., inf],\n",
       "       [inf,  7.,  3.,  3.,  2.,  4., inf],\n",
       "       [inf,  6.,  9.,  2.,  9.,  8., inf]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CME = np.zeros((nrows, ncols+2))\n",
    "CME[:,0] = np.inf\n",
    "CME[:,-1] = np.inf\n",
    "CME[:,1:-1] = energy\n",
    "CME"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Now, what is the least energy needed to get to position (1,0)? \n",
    "  - Well, we could do 9->6 or 4->6\n",
    "  - 4->6 is better. \n",
    "  - So the best we can do is 4+6=10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CME_1_1 = min(CME[0,0], CME[0,1], CME[0,2]) + CME[1,1]\n",
    "CME_1_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- About about for position (1,1)?\n",
    "  - 4->4 is best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CME_1_2 = min(CME[0,1], CME[0,2], CME[0,3]) + CME[1,2]\n",
    "CME_1_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Etc. For the whole row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[inf,  9.,  4.,  9.,  9.,  1., inf],\n",
       "       [inf, 10.,  8., 10.,  9.,  8., inf],\n",
       "       [inf,  1.,  5.,  8.,  9.,  2., inf],\n",
       "       [inf,  7.,  3.,  3.,  2.,  4., inf],\n",
       "       [inf,  6.,  9.,  2.,  9.,  8., inf]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for j in range(1,ncols+1):\n",
    "    CME[1,j] += min(CME[0,j-1],CME[0,j],CME[0,j+1])\n",
    "CME"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Ok! So now the 2nd row is done!\n",
    "- Now we can see why the padding with `np.inf` is useful along the sides:\n",
    "  - This way we can reuse the same code that takes the minimum of 3 values.\n",
    "  - We don't have to deal with the edges as special cases.\n",
    "  - After all, `np.inf` will never be the minimum (all the actual energy values are regular numbers).\n",
    "  - One could code it up without this but the code is a bit nicer this way.\n",
    "  - And you'll use this trick again in the last part of the lab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, let's do this for the rest of the rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[inf,  9.,  4.,  9.,  9.,  1., inf],\n",
       "       [inf, 10.,  8., 10.,  9.,  8., inf],\n",
       "       [inf,  9., 13., 16., 17., 10., inf],\n",
       "       [inf, 16., 12., 16., 12., 14., inf],\n",
       "       [inf, 18., 21., 14., 21., 20., inf]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i in range(2,nrows): # in the real code this loop starts from 1, but we already did row 1 manually\n",
    "    for j in range(1,ncols+1):\n",
    "        CME[i,j] += min(CME[i-1,j-1],CME[i-1,j],CME[i-1,j+1])\n",
    "CME"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Part 1 complete.\n",
    "- Let's take a breath for a moment before continuing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dynamic programming: backtracking (time permitting)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- So what we have so far as the minimum possible energy expenditure from the top to a particular location.\n",
    "- If we now look at the _bottom_ row, we'll see the minimum possible energy _total_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[inf,  9.,  4.,  9.,  9.,  1., inf],\n",
       "       [inf, 10.,  8., 10.,  9.,  8., inf],\n",
       "       [inf,  9., 13., 16., 17., 10., inf],\n",
       "       [inf, 16., 12., 16., 12., 14., inf],\n",
       "       [inf, 18., 21., 14., 21., 20., inf]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.0"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.min(CME[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Apparently this is 14.\n",
    "- OK, so now we have the optimal objective function value, 14.\n",
    "- But we still need the decision variables (seam) that yields this objective value.\n",
    "- We can retrace our steps here, which is called \"backtracking\":\n",
    "  - Start at 14\n",
    "  - Where could we have come from?\n",
    "  - Either the 12, the 16, or the other 12.\n",
    "  - The original energy value of the bottom point was 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9, 4, 9, 9, 1],\n",
       "       [6, 4, 6, 8, 7],\n",
       "       [1, 5, 8, 9, 2],\n",
       "       [7, 3, 3, 2, 4],\n",
       "       [6, 9, 2, 9, 8]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- So this 14 must have come from $12+2$. \n",
    "- In this (unusual) case there are two different 12s, typically it'd be unique.\n",
    "- What we're seeing here is two seams that are equally good!\n",
    "- Let's take the left branch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[inf,  9.,  4.,  9.,  9.,  1., inf],\n",
       "       [inf, 10.,  8., 10.,  9.,  8., inf],\n",
       "       [inf,  9., 13., 16., 17., 10., inf],\n",
       "       [inf, 16., 12., 16., 12., 14., inf],\n",
       "       [inf, 18., 21., 14., 21., 20., inf]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CME"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- So now we're at the 12.\n",
    "- Above it we see 9, 13, 16.\n",
    "- We could again look at the energy value there, which happens to be 3, meaning we came from the 9.\n",
    "- But actually there's an easier way:\n",
    "  - Where we came from must be the lowest of the three possible CME values (9, 13, 16)\n",
    "  - Because that is exacly how we computed the CME of 12 in the first place!!\n",
    "  - So we don't need to look and find that the energy there was 3, we just know it to be the case.\n",
    "- We continue this process all the way back up:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n"
     ]
    }
   ],
   "source": [
    "seam = np.argmin(CME[-1])\n",
    "print(seam)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.argmin(CME[-2, seam-1:seam+2]) - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means go left."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seam = seam - 1\n",
    "np.argmin(CME[-3, seam-1:seam+2]) - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means go left"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seam = seam - 1\n",
    "np.argmin(CME[-4, seam-1:seam+2]) - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means go right"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seam = seam + 1\n",
    "np.argmin(CME[-5, seam-1:seam+2]) - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- This means go up.\n",
    "- So from the bottom, we started at position 3 of the modified array.\n",
    "  - This is potition 2 in the original array.\n",
    "  - left, left, right, up:\n",
    "  - 2 -> 1 -> 0 -> 1 -> 1\n",
    "  - But this is backwards so it's [1, 1, 0, 1, 2]\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 0, 1, 2])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "find_vertical_seam(energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Good thing I don't ask you to write the code anymore, I think that was a bit unecessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### When can we use dynamic programming?\n",
    "\n",
    "- When can we use dynamic programming? \n",
    "  - Sequential decision-making: one can make an optimal decision by taking the optimal solution to a subproblem, and looking only at the next step.\n",
    "- Question: Can we use dynamic programming to solve our MDS TA assignment problem?\n",
    "\n",
    "\n",
    "<br><br><br>\n",
    "- Answer: No. If we have assigned 5 TAs optimally, we cannot just assign the 6th TA optimally without looking at the rest of the TA pool. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code Vectorization (20 min)\n",
    "\n",
    "Consider the code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10**7\n",
    "x = np.zeros(n)\n",
    "y = np.zeros(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.03 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -n1 -r1\n",
    "\n",
    "for i in range(n):\n",
    "    y[i] = x[i] + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.zeros(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16.3 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -n1 -r1\n",
    "\n",
    "y = x + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Wow, this this ran a lot faster! Why?!\n",
    "- This is called **vectorization**.\n",
    "- (Optional) This has to do with the Python language being _interpreted_ rather than _compiled_.\n",
    "  - You may have worked with other languages, like C or Java, that had a compilation step.\n",
    "  - Note: R is also interpreted.\n",
    "- Numpy has some super-optimized code under the hood, for common operations like vector addition\n",
    "  - Writing `x + 1` tells the interpreter to do vector addition which is fast.\n",
    "  - Writing the loop has Python execute each step one at a time.\n",
    "  - It doesn't \"know in advance\" that it will be a vector addition.\n",
    "- There are some ways around this, discussed below.\n",
    "  - There are also some languages that don't suffer from this issue, like [Julia](https://en.wikipedia.org/wiki/Julia_(programming_language)).\n",
    "- But in general, it's often good to vectorize your code when possible.\n",
    "- The speedups can be huge. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Case study: pairwise distances\n",
    "\n",
    "Consider two vectors (1D numpy arrays) $x$ and $y$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 100\n",
    "\n",
    "x = np.random.rand(d)\n",
    "y = np.random.rand(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The squared Euclidean distance between two vectors $x$ and $y$ is defined as\n",
    "\n",
    "$$\\sum_{i=1}^d (x_i - y_i)^2$$\n",
    "\n",
    "We definitely don't want to implement this with a loop - we can use `np.sum`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15.337953373218195"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum((x-y)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Now imagine we have a bunch of $x$ vectors and a bunch of $y$ vectors, and we want to compute the squared distance between all pairs.\n",
    "  - For example, think of KNN classification where you have a bunch of train vectors and a bunch of test vectors.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx = 5\n",
    "ny = 3\n",
    "\n",
    "X = np.random.rand(nx,d) # we have nx of them\n",
    "Y = np.random.rand(ny,d) # we have ny of them"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we have $5$ $x$-vectors and $3$ $y$-vectors, so we will have $5\\times 3 = 15$ distance pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[16.45044963, 16.57440775, 16.42308959],\n",
       "       [18.28721673, 15.00038461, 13.02764908],\n",
       "       [13.08522752, 17.96707571, 16.19017444],\n",
       "       [15.87051195, 16.95390091, 14.66619765],\n",
       "       [11.81458429, 18.42825745, 15.15205817]])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "squared_distances = np.zeros((nx,ny))\n",
    "\n",
    "for i in range(nx):\n",
    "    for j in range(ny):\n",
    "        squared_distances[i,j] = np.sum((X[i]-Y[j])**2)\n",
    "        \n",
    "squared_distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's wrap everything in a function for timing purposes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pairwise_squared_distances(X, Y):\n",
    "    nx = len(X)\n",
    "    ny = len(Y)\n",
    "    \n",
    "    squared_distances = np.zeros((nx,ny))\n",
    "    for i in range(nx):\n",
    "        for j in range(ny):\n",
    "            squared_distances[i,j] = np.sum((X[i]-Y[j])**2)\n",
    "\n",
    "    return squared_distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[18.38917922, 18.07747894, 17.96680562],\n",
       "       [16.50135296, 15.53137547, 16.33320306],\n",
       "       [15.40252805, 16.05500748, 17.45413749],\n",
       "       [16.24033663, 14.38279876, 19.69047128],\n",
       "       [16.59927659, 15.30817853, 15.75408885]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.random.rand(5, 100) \n",
    "Y = np.random.rand(3, 100)\n",
    "\n",
    "pairwise_squared_distances(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xbig = np.random.rand(1000, 100) \n",
    "Ybig = np.random.rand(1000, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.4 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 pairwise_squared_distances(Xbig, Ybig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question: what is the running time of `pairwise_squared_distances`?\n",
    "<br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The running time here is $O(dn_xn_y)$ because of the loops over $n_x$ and $n_y$, and the implicit loop over $d$ inside of `sum`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Broadcasting\n",
    "\n",
    "- You learned about broadcasting in DSCI 511. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 100)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 100)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 1, 100)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[:,np.newaxis].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 3, 100)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y[np.newaxis].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 3, 100)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(X[:,None] - Y[None]).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.90410728],\n",
       "       [0.76983655],\n",
       "       [0.12040406],\n",
       "       [0.61314058],\n",
       "       [0.76914491]])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xx = np.random.rand(5,1)\n",
    "xx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.47520232, 0.63405477, 0.09001631]])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yy = np.random.rand(1,3)\n",
    "yy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.42890497,  0.27005252,  0.81409098],\n",
       "       [ 0.29463423,  0.13578178,  0.67982024],\n",
       "       [-0.35479825, -0.5136507 ,  0.03038776],\n",
       "       [ 0.13793826, -0.02091419,  0.52312427],\n",
       "       [ 0.29394259,  0.13509014,  0.6791286 ]])"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xx-yy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 3, 100)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(X[:,None] - Y[None]).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.2762265262237462"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(X[:,None] - Y[None])[1,2,3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.2762265262237462"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[1,3] - Y[2,3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pairwise_squared_distances_broadcast(X, Y):    \n",
    "    return np.sum((X[:,None] - Y[None])**2,axis=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[18.38917922, 18.07747894, 17.96680562],\n",
       "       [16.50135296, 15.53137547, 16.33320306],\n",
       "       [15.40252805, 16.05500748, 17.45413749],\n",
       "       [16.24033663, 14.38279876, 19.69047128],\n",
       "       [16.59927659, 15.30817853, 15.75408885]])"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pairwise_squared_distances(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[18.38917922, 18.07747894, 17.96680562],\n",
       "       [16.50135296, 15.53137547, 16.33320306],\n",
       "       [15.40252805, 16.05500748, 17.45413749],\n",
       "       [16.24033663, 14.38279876, 19.69047128],\n",
       "       [16.59927659, 15.30817853, 15.75408885]])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pairwise_squared_distances_broadcast(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.42 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 pairwise_squared_distances(Xbig, Ybig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "781 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 pairwise_squared_distances_broadcast(Xbig, Ybig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- This gives the same result and is around 5x faster!!\n",
    "- BTW, here is the worst implementation, with no vectorization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pairwise_squared_distances_so_slow(X, Y):\n",
    "    nx = len(X)\n",
    "    ny = len(Y)\n",
    "    d = X.shape[1]\n",
    "    \n",
    "    squared_distances = np.zeros((nx,ny))\n",
    "    for i in range(nx):\n",
    "        for j in range(ny):\n",
    "            for k in range(d):\n",
    "                squared_distances[i,j] += (X[i,k]-Y[j,k])**2\n",
    "\n",
    "    return squared_distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[18.38917922, 18.07747894, 17.96680562],\n",
       "       [16.50135296, 15.53137547, 16.33320306],\n",
       "       [15.40252805, 16.05500748, 17.45413749],\n",
       "       [16.24033663, 14.38279876, 19.69047128],\n",
       "       [16.59927659, 15.30817853, 15.75408885]])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pairwise_squared_distances_so_slow(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %timeit -n1 -r1 pairwise_squared_distances_so_slow(Xbig, Ybig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The above takes around 3 minutes on my laptop, I don't want to run it live."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (optional) An even faster approach\n",
    "\n",
    "- Turns out there's a function that does this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[18.38917922, 18.07747894, 17.96680562],\n",
       "       [16.50135296, 15.53137547, 16.33320306],\n",
       "       [15.40252805, 16.05500748, 17.45413749],\n",
       "       [16.24033663, 14.38279876, 19.69047128],\n",
       "       [16.59927659, 15.30817853, 15.75408885]])"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sklearn.metrics.pairwise.euclidean_distances(X, Y, squared=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20.9 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 sklearn.metrics.pairwise.euclidean_distances(Xbig, Ybig, squared=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- How is it so fast?! That is almost 100x faster than broadcasting?!\n",
    "- This is a specific trick to this problem, using this identity:\n",
    "\n",
    "$$(x-y)^2=x^2+y^2-2xy$$\n",
    "\n",
    "If $x$ and $y$ are vectors, this becomes\n",
    "\n",
    "$$(x-y)^T(x-y)=x^Tx+y^Ty-2x^Ty$$\n",
    "\n",
    "- If we want to do this for all pairs for vectors, we can write it using matrix multiplication $XX^T + YY^T - 2XY^T$\n",
    "- Matrix multiplication is extremely fast, both because of the implementations and because of a fast algorithm called the [Strassen algorithm](https://en.wikipedia.org/wiki/Strassen_algorithm)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(end optional)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### When should I vectorize my code?\n",
    "\n",
    "- If you see a loop, ask yourself if the next iteration depends on the result from the previous one.\n",
    "- If not, try to vectorize it.\n",
    "- This is not so simple - when summing an array, does each iteration depend on the previous one?\n",
    "   - Yes, because you keep adding to the cumulative sum.\n",
    "   - But no, because you could, e.g., recursively sum pairs of numbers.\n",
    "   - This question is not always that simple, but often you will be able to tell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Faster algorithms vs. faster implementations\n",
    "\n",
    "- Sometimes we can speed up our code using a faster _algorithm_, like dynamic programming.\n",
    "  - This usually means that it's actually different in big O, or the number of operations.\n",
    "- Other times we speed it up using a faster _implementation_, like vectorization.\n",
    "  - More of these approaches coming below.\n",
    "- **These are two distinct approaches to speeding up your code.**\n",
    "  - But sometimes they can be combined, e.g. in the super-fast pairwise distances, or in your lab this week."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is all we have time for, but there are some bonus materials below on:\n",
    "\n",
    "- Profiling your code (to tell you what lines are slow)\n",
    "- Parallelizing code execution\n",
    "- Other tools to speed up Python\n",
    "- A really cool (but trickier) dynamic programming example: computing diffs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Course wrap-up (5 min)\n",
    "\n",
    "Course learning goals:\n",
    "\n",
    "- Analyze the scalability and trade-offs of various basic algorithms and data structures using Big-O notation.\n",
    "- Read and interpret recursive functions.\n",
    "- Select appropriate data structures, such as graphs, given a data set.\n",
    "- Map certain real-world problems to discrete optimization problems.\n",
    "- Diagnose rate-limiting aspects of slow Python code and enumerate various options for speeding it up.\n",
    "\n",
    "Hopefully this course has helped you more deeply understand some data structures you were already using, like sets, and also unlocked some new ways of thinking, like recursion or optimization.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------------------\n",
    "\n",
    "Everything below this line is optional / we don't have time for it.\n",
    "\n",
    "------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext snakeviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import njit, prange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsClassifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (optional) Profiling (5 min)\n",
    "\n",
    "- Profiling means measuring how long parts of your code takes.\n",
    "- We'll use [SnakeViz](https://jiffyclub.github.io/snakeviz/) for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "knn = KNeighborsClassifier(algorithm='brute')\n",
    "knn.fit(Xbig, np.zeros(Xbig.shape[0]));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file '/var/folders/sl/_n6p71ds1vb5tr2zmd73gnpr0000gp/T/tmp7_r80qv2'. \n",
      "Opening SnakeViz in a new tab...\n"
     ]
    }
   ],
   "source": [
    "%%snakeviz -t \n",
    "\n",
    "# Note: -t is to open in new tab which is necessary in JupyterLab (thanks to Joel for this tip)\n",
    "\n",
    "knn.predict(Ybig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Here, we can see that a lot of time was spent in the `pairwise_distances_chunked` function.\n",
    "- On the other hand, if we remove `algorithm='brute'` we'll see the code runs faster and we'll also see different profile results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Profiling is useful when your code is slow but you don't know which parts.\n",
    "- \"Premature optimization is the root of all evil\" -Donald Knuth\n",
    "  - So, make sure you know what part of your code needs speeding up.\n",
    "  - If something is only taking 1% of your code's total time, don't work on making it faster."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (optional) Other tools for speeding up Python code (5 min)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Numba\n",
    "\n",
    " \n",
    "- \"[Numba](http://numba.pydata.org/) is an open source JIT compiler that translates a subset of Python and NumPy code into fast machine code.\"\n",
    "  - \"JIT compiler\" = Just-in-time compiler.\n",
    "  - It sees a loop and figures out that it's actually vector addition, right before the loop is executed (hence \"just in time\").\n",
    "  - This is the easiest to deploy, because you don't need to change your code - just add a decorator.\n",
    "\n",
    "Recall that this code took about 7 seconds last time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10**7\n",
    "x = np.zeros(n)\n",
    "y = np.zeros(n)\n",
    "\n",
    "def run_code(x, y):\n",
    "    for i in range(len(x)):\n",
    "        y[i] = x[i] + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.59 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 run_code(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "@njit\n",
    "def run_code_jit(n, x, y):\n",
    "    for i in range(n):\n",
    "        y[i] = x[i] + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "397 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 run_code_jit(n, x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Bam! It's just that easy.\n",
    "- Note: with the `@` we're again using a decorator, like we did with the caching at the start of class.\n",
    "- This is rather confusing because `@` can also denote matrix multiplication..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is amazing - why isn't it the default?\n",
    "\n",
    "- This goes back to the idea that Python wasn't designed for numerical computation. \n",
    "- This is a 3rd party install\n",
    "- In this case we could have also vectorized our code instead."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cython\n",
    "\n",
    "- [Cython](https://cython.org/) allows the user to leverage the increased speed of C within Python.\n",
    "- You can call C routines.\n",
    "- You actually compile your Cython code - it is not interpreted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### PyPy\n",
    "\n",
    "- [PyPy](https://pypy.org/) is an implementation of Python in Python.\n",
    "  - The original Python, or CPython, is written in C.\n",
    "  - PyPy has several advantages including a just-in-time compiler which might make your code run faster.\n",
    "- You can install PyPy and run your code with `pypy test.py` instead of `python test.py`.\n",
    "  - You don't need to change your code, BUT some libraries may not work with PyPy. Things like numpy do work. For more info, see [here](https://pypy.org/compat.html).\n",
    "\n",
    "Note! This has nothing to do with [PyPI](https://pypi.org/), which is the package index, like [CRAN](https://cran.r-project.org/) for R. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (optional) Parallelization (5 min)\n",
    "\n",
    "- [Moore's law](https://en.wikipedia.org/wiki/Moore%27s_law) states that computing power (number of transistors on a chip) will double approximately every 2 years. \n",
    "  - Shockingly, this held true for ~40 years.\n",
    "  - Recently, progress has tapered off.\n",
    "  - However, recently progress has been in **parallel** and **distributed** computation.\n",
    "\n",
    "Consider our Fibonacci code from earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fib(n):\n",
    "    f = np.zeros(n+1, dtype=int)\n",
    "    f[1] = 1\n",
    "    \n",
    "    for i in range(2,n+1):\n",
    "        f[i] = f[i-1] + f[i-2]\n",
    "    \n",
    "    return f[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The running time here is $O(n)$.\n",
    "- Now consider our addition code:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_one(n, x):\n",
    "    y = np.zeros(n)\n",
    "    for i in range(n):\n",
    "        y[i] = x[i] + 1\n",
    "        \n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The running time here is also $O(n)$.\n",
    "- But this one can be done in parallel, because the $n$ computations _do not depend on each other_.\n",
    "- On the other hand, each Fibonacci iteration depends on the previous two values.\n",
    "- Both are $O(n)$, but they are fundamentally different.\n",
    "- The requirements for code to be parallelized and vectorized are similar, though the two are not the same thing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hardware for parallelization\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import multiprocessing\n",
    "\n",
    "multiprocessing.cpu_count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- My laptop has 8 cores. \n",
    "- I should theoretically be able to use them to run the code in parallel and achieve (close to) a 8x speedup.\n",
    "- Recently, [GPUs](https://en.wikipedia.org/wiki/Graphics_processing_unit) have become extremely popular for parallel computing, especially in deep learning.\n",
    "  - We'll talk about this in DSCI 572."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Software for parallelization\n",
    "\n",
    "- Numba can do some of this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 50_000_000\n",
    "\n",
    "x = np.zeros(n)\n",
    "y = np.zeros(n)\n",
    "\n",
    "@njit\n",
    "def run_code(x, y):\n",
    "    for i in range(len(x)):\n",
    "        y[i] = x[i] + np.random.rand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "951 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 run_code(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "@njit(parallel=True)\n",
    "def run_code_parallel(x, y):\n",
    "    for i in prange(len(x)): # <- note the prange instead of range\n",
    "        y[i] = x[i] + np.random.rand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "787 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 run_code_parallel(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- In this case the speedup is modest (I don't know the details), but it can be more.\n",
    "- Even a 2x speedup is great if your code takes 4 hours instead of 8 hours."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We saw [joblib](https://joblib.readthedocs.io/en/latest/) earlier, for on-disk memoization.\n",
    "- It can also be used for parallelization, but we won't go into it here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Parallel vs. distributed\n",
    "\n",
    "- Parallel computing can happen even on one machine, e.g. my 4 cores\n",
    "- Distributed computing means sharing the computation across machines.\n",
    "  - We'll talk about this in DSCI 525."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parallel computing is an entire field, so we don't do it much justice in 5 min. But this was hopefully better than nothing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus dynamic programming application: diffs between strings\n",
    "\n",
    "- When you look at a commit in GitHub, you just see the \"diff\" or difference between the two versions.\n",
    "\n",
    "<img src=\"diff.png\" width=\"600\">\n",
    "\n",
    "- Note the highlighted areas: git/GitHub knows about the added whitespace, changing \"see\" to \"SEE\", and the deleted exclamation mark.\n",
    "  - How does GitHub compute the diff?\n",
    "  - We want to highlight _characters that are not present in the other version_.\n",
    "  - (Optional note: I believe it's GitHub, rather than git, that's doing this. See [here](https://stackoverflow.com/questions/10398744/does-git-store-diff-information-in-commit-objects).)\n",
    "- Why is this an optimization problem?\n",
    "  - You can always say the diff is \"deleted the entire old file, added the entire new file\".\n",
    "  - But that is not a useful diff!\n",
    "  - What we are doing implicitly is _aligning_ the two sequences.\n",
    "  - We want an alignment with the _minimal_ changes that go from $x$ to $y$, i.e. we want to _minimize the number of highlighted characters._\n",
    "  - Now we have an optimization problem!\n",
    "- We have 4 highlighted characters in $x$ (red) and 6 highlighted characters in $y$ (green), for a total of 10."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Back to the example:\n",
    "\n",
    "<img src=\"diff_cropped.png\" width=\"600\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We're really _aligning_ the two strings\n",
    "\n",
    "- When we compare the strings and see a mismatch, how do we know if we should highlight in red or green?\n",
    "  - What we're really doing here is an _alignment_; it has to be the green.\n",
    "  - We have a binary decision to make everytime we hit a mismatch.\n",
    "- But we can't just do this one character at a time: the highlighting affects the alignment!\n",
    "  - We need to look globally at the whole of strings.\n",
    "- We want something that is optimal across _all possible alignments_.\n",
    "  - This brings us to the brute force approach: recursively check all possible alignments. \n",
    "- Optional note: this is a reformulation of the [longest common subsequence problem](https://en.wikipedia.org/wiki/Longest_common_subsequence_problem) and is also a form of [edit distance](https://en.wikipedia.org/wiki/Edit_distance)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "def num_diffs_rec(x, y):\n",
    "    \"\"\"\n",
    "    Find the number of characters in the diff between x and y.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x : str\n",
    "        The first string\n",
    "    y : str\n",
    "        The second string\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    int\n",
    "        The number of highlights\n",
    "        \n",
    "    Examples\n",
    "    --------\n",
    "    >>> num_diffs_rec(\"This is a!\", \"this  is a\")\n",
    "    4\n",
    "    >>> num_diffs_rec(\"xxHello\", \"Hellox\")\n",
    "    3\n",
    "    \"\"\"    \n",
    "    if len(x) == 0:\n",
    "        return len(y) # Highlight the rest of y\n",
    "    if len(y) == 0:\n",
    "        return len(x) # Highlight the rest of x\n",
    "    \n",
    "    if x[0] == y[0]:  # A match\n",
    "        return num_diffs_rec(x[1:], y[1:]) \n",
    "    else:\n",
    "        return 1 + min( num_diffs_rec(x[1:], y), num_diffs_rec(x, y[1:]) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_diffs_rec(\"This is a!\", \"this  is a\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_diffs_rec(\"xxHello\", \"Hellox\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: the code that actually returns the highlighted characters is slightly more complicated -- we'll defer that until later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How does this work?\n",
    "\n",
    "Back to the example:\n",
    "\n",
    "```\n",
    "This is a demonstration of diffs in git/GitHub. Let's see if it works!\n",
    "This    is a demonstration of diffs in git/GitHub. Let's SEE if it works\n",
    "```\n",
    "\n",
    "- The first mismatch is between `i` and ` `. \n",
    "- Remember, we're looking for a _minimum_ number of highlights.\n",
    "- The minimum total number of highlights is the lesser of two things:\n",
    "   1. The score if you highlight `i`, plus 1 for highlighting the `i`.\n",
    "   2. The score if you highlight ` `, plus 1 for highlighting the ` `.\n",
    "- So we do this recursively and that's it!\n",
    "- The base case is that if you reach the end of one string, you have to highlight the rest of the other string."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How slow is this code?\n",
    "\n",
    "- As usual, the brute force solution is unusably slow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30.9 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 num_diffs_rec(\"This'll be slow\", \"Yeah, right...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This implementation is pretty much useless!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Why is the code so slow?\n",
    "\n",
    "- This is the same issue as we saw earlier with the Fibonacci numbers.\n",
    "- We are repeating computations for subproblems.\n",
    "- We won't do it now, but I encourage you to draw out the recursion and see for yourself!\n",
    "- We could use the memoization trick again, but:\n",
    "  - We don't have a clear idea of how much memory it will take - will it be too much?\n",
    "  - We don't have a clear idea of how long it will run.\n",
    "  - It's a bit heavy handed and we cannot clearly see what the code is doing.\n",
    "  - It would be better to just directly write good code for this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Enter dynamic programming \n",
    "\n",
    "- Dynamic programming is an algorithm for solving certain optimization problems.\n",
    "  - It only works for some problems, but when it works it is FAST.\n",
    "  - This is just like linear programming (last lecture), except for a different set of problems.\n",
    "  - It is closely related to memoization; see [here](https://stackoverflow.com/questions/6184869/what-is-the-difference-between-memoization-and-dynamic-programming) for a discussion of memoization vs. dynamic programming. \n",
    "- Dynamic programming has a lot of applications:\n",
    "  - DNA sequence alignment\n",
    "  - Text hyphenation\n",
    "  - Running certain machine learning models (e.g. hidden Markov models)\n",
    "  - Finding the differences between two files (this!)\n",
    "  - Image resizing (this week's lab!)\n",
    "  - and more...\n",
    "- What does the name \"Dynamic Programming\" mean?\n",
    "  - \"Programming\" is often used for optimization, for historical reasons.\n",
    "  - There is no significance to \"Dynamic\", it just sounds good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dynamic programming implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "def num_diffs(x, y, return_table=False):\n",
    "    \"\"\"\n",
    "    Compute the number of highlighted characters in the\n",
    "    diff between x and y using dynamic programming.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x : str\n",
    "        The first string\n",
    "    y : str\n",
    "        The second string\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    numpy.ndarray\n",
    "        The dynamic programming table. \n",
    "        The last element is the result.\n",
    "        \n",
    "    Examples\n",
    "    --------\n",
    "    >>> num_diffs(\"This is a!\", \"this  is a\")[-1,-1]\n",
    "    4\n",
    "    >>> num_diffs(\"xxHello\", \"Hellox\")[-1,-1]\n",
    "    3\n",
    "    \"\"\"   \n",
    "    M = len(x)\n",
    "    N = len(y)\n",
    "    \n",
    "    opt = np.zeros((M+1, N+1), dtype=int)\n",
    "    opt[:,0] = np.arange(M+1)\n",
    "    opt[0,:] = np.arange(N+1)\n",
    "    \n",
    "    for i in range(1,M+1):\n",
    "        for j in range(1,N+1):\n",
    "            if x[i-1] == y[j-1]:\n",
    "                opt[i,j] = opt[i-1, j-1]\n",
    "            else:\n",
    "                opt[i,j] = 1 + min( opt[i-1,j], opt[i,j-1] )\n",
    "\n",
    "    return opt if return_table else opt[-1,-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_diffs(\"This is a!\", \"this  is a\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_diffs(\"xxHello\", \"Hellox\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How does this work?\n",
    "\n",
    "- The recursive implementation works \"top-down\":\n",
    "  - It starts from the big problem (the entire two strings) and computes the solutions for subproblems.\n",
    "- Dynamic programming works \"bottom-up\":\n",
    "  - It starts from the smallest subproblems and builds to the bigger ones.\n",
    "  - This way it doesn't recompute things multiple times.\n",
    "  - It is basically the recursive solution + memoization, but implemented deliberately.\n",
    "  \n",
    "To be concrete:\n",
    "\n",
    "- We defined a 2D array which we called `opt`. \n",
    "- `opt[i,j]` contains `num_diffs(x[:i], y[:j])`, that is the solution for the $(i,j)$ subproblem\n",
    "- The key logic is as follows:\n",
    "  - If `x[i]` equals `y[j]`, then `opt[i,j]` equals `opt[i-1,j-1]`\n",
    "    - Because we don't add any additional highlights by adding that next character!\n",
    "  - If they are not equal, then there are two possibilities: highlight in $x$ or highlight in $y$\n",
    "    - We want the **better** of those two possibilities.\n",
    "    - Good news is, we already know how good they are! This is the subproblem part.\n",
    "  - This is how we build up the solution iteratively instead of repeating computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3, 4, 5, 6],\n",
       "       [1, 2, 3, 4, 5, 6, 5],\n",
       "       [2, 3, 4, 5, 6, 7, 6],\n",
       "       [3, 2, 3, 4, 5, 6, 7],\n",
       "       [4, 3, 2, 3, 4, 5, 6],\n",
       "       [5, 4, 3, 2, 3, 4, 5],\n",
       "       [6, 5, 4, 3, 2, 3, 4],\n",
       "       [7, 6, 5, 4, 3, 2, 3]])"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, y = \"xxHello\", \"Hellox\"\n",
    "\n",
    "opt = num_diffs(x, y, return_table=True)\n",
    "opt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's recreate this table according to its definition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3, 4, 5, 6],\n",
       "       [1, 2, 3, 4, 5, 6, 5],\n",
       "       [2, 3, 4, 5, 6, 7, 6],\n",
       "       [3, 2, 3, 4, 5, 6, 7],\n",
       "       [4, 3, 2, 3, 4, 5, 6],\n",
       "       [5, 4, 3, 2, 3, 4, 5],\n",
       "       [6, 5, 4, 3, 2, 3, 4],\n",
       "       [7, 6, 5, 4, 3, 2, 3]])"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "same = 0*opt\n",
    "for i in range(len(x)+1):\n",
    "    for j in range(len(y)+1):\n",
    "        same[i,j] = num_diffs_rec(x[:i], y[:j])\n",
    "same"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- This illustrates how silly the recursive solution is!!\n",
    "- Each element of this table `same` was computed from scratch.\n",
    "- But actually we can get it in $O(1)$ time given the previous elements.\n",
    "  - That's how the dynamic programming code works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computational cost\n",
    "\n",
    "- Now we know how much it costs to compute this:\n",
    "  - Memory usage is $O(MN)$\n",
    "  - Runtime is also $O(MN)$\n",
    "  - (where $M$ is the length of $x$ and $N$ is the length of $y$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Recovering the highlights\n",
    "\n",
    "- This is called \"backtracking\" in dynamic programming.\n",
    "- We need access the table and make sense of it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3, 4, 5, 6],\n",
       "       [1, 2, 3, 4, 5, 6, 5],\n",
       "       [2, 3, 4, 5, 6, 7, 6],\n",
       "       [3, 2, 3, 4, 5, 6, 7],\n",
       "       [4, 3, 2, 3, 4, 5, 6],\n",
       "       [5, 4, 3, 2, 3, 4, 5],\n",
       "       [6, 5, 4, 3, 2, 3, 4],\n",
       "       [7, 6, 5, 4, 3, 2, 3]])"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Intuitively, we check where we might have come from at the previous step.\n",
    "- We start at the lower-right and \"backtrack\" to the upper left.\n",
    "- The function below has some extra HTML visualization code that you can ignore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import HTML, display\n",
    "\n",
    "dark_green  = '<span style=\"background-color: rgba(0,255,0,0.5)\">'\n",
    "dark_red    = '<span style=\"background-color: rgba(255,0,0,0.5)\">'\n",
    "light_green = '<span style=\"background-color: rgba(0,255,0,0.05)\">'\n",
    "light_red   = '<span style=\"background-color: rgba(255,0,0,0.05)\">'\n",
    "\n",
    "def show_diff(x, y, align=False):\n",
    "    opt = num_diffs(x, y, return_table=True)\n",
    "    \n",
    "    x_highlight = ''\n",
    "    y_highlight = ''\n",
    "    i = len(x) \n",
    "    j = len(y)\n",
    "    while i > 0 or j > 0:\n",
    "        if i > 0 and j > 0 and x[i-1] == y[j-1]:\n",
    "            x_highlight = x[i-1] + x_highlight\n",
    "            y_highlight = y[j-1] + y_highlight\n",
    "            i -= 1\n",
    "            j -= 1\n",
    "        elif j > 0 and opt[i, j] == opt[i, j-1] + 1:\n",
    "            y_highlight = dark_green + y[j-1] + '</span>' + y_highlight \n",
    "            if align:\n",
    "                x_highlight = \" \" + x_highlight\n",
    "            j -= 1\n",
    "        else:\n",
    "            x_highlight = dark_red   + x[i-1] + '</span>' + x_highlight \n",
    "            if align:\n",
    "                y_highlight = \" \" + y_highlight\n",
    "            i -= 1\n",
    "    \n",
    "    x_highlight = light_red   + x_highlight + \"</span>\"\n",
    "    y_highlight = light_green + y_highlight + \"</span>\"\n",
    "    \n",
    "    display(HTML('<code>' + x_highlight + '</code>' + '<br>' + '<code>' + y_highlight + '</code>'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<code><span style=\"background-color: rgba(255,0,0,0.05)\"><span style=\"background-color: rgba(255,0,0,0.5)\">x</span><span style=\"background-color: rgba(255,0,0,0.5)\">x</span>Hello</span></code><br><code><span style=\"background-color: rgba(0,255,0,0.05)\">Hello<span style=\"background-color: rgba(0,255,0,0.5)\">x</span></span></code>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show_diff(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note: the colours do not render properly on GitHub - you need to run the notebook locally to see the highlighting.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<code><span style=\"background-color: rgba(255,0,0,0.05)\">This is a demonstration of diffs in git/GitHub. Let's <span style=\"background-color: rgba(255,0,0,0.5)\">s</span><span style=\"background-color: rgba(255,0,0,0.5)\">e</span><span style=\"background-color: rgba(255,0,0,0.5)\">e</span> if it works<span style=\"background-color: rgba(255,0,0,0.5)\">!</span></span></code><br><code><span style=\"background-color: rgba(0,255,0,0.05)\">This<span style=\"background-color: rgba(0,255,0,0.5)\"> </span><span style=\"background-color: rgba(0,255,0,0.5)\"> </span><span style=\"background-color: rgba(0,255,0,0.5)\"> </span> is a demonstration of diffs in git/GitHub. Let's <span style=\"background-color: rgba(0,255,0,0.5)\">S</span><span style=\"background-color: rgba(0,255,0,0.5)\">E</span><span style=\"background-color: rgba(0,255,0,0.5)\">E</span> if it works</span></code>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "before = \"This is a demonstration of diffs in git/GitHub. Let's see if it works!\"\n",
    "after  = \"This    is a demonstration of diffs in git/GitHub. Let's SEE if it works\"\n",
    "\n",
    "show_diff(before, after)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br><br>vs.<br><br>\n",
    "\n",
    "<img src=\"diff_cropped.png\" width=\"600\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- They match.\n",
    "- Well, with a slight difference on highlighting of the spaces.\n",
    "- This is because the solution to the optimization problem is not unique - both are optimal.\n",
    "- Arbitrary choices made in the code will determine which one is returned by the code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Why could we use dynamic programaming to find the diffs between documents?\n",
    "  - This comes back to the discussion of how `opt` is computed.\n",
    "  - The optimal solution at a given step can be described in terms of optimal solutions of subproblems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3, 4, 5, 6],\n",
       "       [1, 2, 3, 4, 5, 6, 5],\n",
       "       [2, 3, 4, 5, 6, 7, 6],\n",
       "       [3, 2, 3, 4, 5, 6, 7],\n",
       "       [4, 3, 2, 3, 4, 5, 6],\n",
       "       [5, 4, 3, 2, 3, 4, 5],\n",
       "       [6, 5, 4, 3, 2, 3, 4],\n",
       "       [7, 6, 5, 4, 3, 2, 3]])"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, y = \"xxHello\", \"Hellox\"\n",
    "\n",
    "opt = num_diffs(x, y, return_table=True)\n",
    "opt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider `opt[3,1]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'xxH'"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'H'"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y[:1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, because the last letters match (both `H`), we can just grab the score from the diagonal entry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt[3,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt[2,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider `opt[5,2]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'xxHel'"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'He'"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y[:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know this is the right alignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<code><span style=\"background-color: rgba(255,0,0,0.05)\"><span style=\"background-color: rgba(255,0,0,0.5)\">x</span><span style=\"background-color: rgba(255,0,0,0.5)\">x</span>He<span style=\"background-color: rgba(255,0,0,0.5)\">l</span></span></code><br><code><span style=\"background-color: rgba(0,255,0,0.05)\">He</span></code>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show_diff(x[:5], y[:2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But how do we get there? \n",
    "\n",
    "- We need to look at two options: the `l` is highlighted or the `e` is highlighted.\n",
    "- Since they don't match, one of them has to be highlighted.\n",
    "- What's confusing here is that, in the final result, neither is highlighted.\n",
    "  - But for this subproblem, one of them has to be!\n",
    "  - This \"line of reasoning\" isn't going to end up being used in the final result at all.\n",
    "  - But we still have to compute it because we don't know in advance what will end up being optimal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<code><span style=\"background-color: rgba(255,0,0,0.05)\"><span style=\"background-color: rgba(255,0,0,0.5)\">x</span><span style=\"background-color: rgba(255,0,0,0.5)\">x</span>He</span></code><br><code><span style=\"background-color: rgba(0,255,0,0.05)\">He</span></code>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show_diff(x[:4], y[:2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt[4,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<code><span style=\"background-color: rgba(255,0,0,0.05)\"><span style=\"background-color: rgba(255,0,0,0.5)\">x</span><span style=\"background-color: rgba(255,0,0,0.5)\">x</span>H<span style=\"background-color: rgba(255,0,0,0.5)\">e</span><span style=\"background-color: rgba(255,0,0,0.5)\">l</span></span></code><br><code><span style=\"background-color: rgba(0,255,0,0.05)\">H</span></code>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show_diff(x[:5], y[:1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt[5,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt[5,2]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 ('dsci512')",
   "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.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "833530e805897aab5e7a30ee530b0bcb4cd19373b840f29eae5897b61b26002e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}