Advanced NumPy Patterns.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# About me\n",
    "\n",
    "Juan Nunez-Iglesias  \n",
    "CZI Imaging Software Fellow  \n",
    "Monash University  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "http://elegant-scipy.org\n",
    "\n",
    "<a href=\"http://elegant-scipy.org\">\n",
    "<img src=\"https://github.com/elegant-scipy/elegant-scipy/blob/master/_images/cover.jpg?raw=true\"\n",
    " alt=\"Elegant SciPy Cover\" height=128 width=192 align=\"left\">\n",
    "</a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Quick example: gene expression, without numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Gene   | Cell type A | Cell type B | Cell type C | Cell type D |\n",
    "|--------|-------------|-------------|-------------|-------------|\n",
    "| Gene 0 | 100         | 200         | 50          | 400         |\n",
    "| Gene 1 | 50          | 0           | 0           | 100         |\n",
    "| Gene 2 | 350         | 100         | 50          | 200         |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene0 = [100, 200, 50, 400]\n",
    "gene1 = [50, 0, 0, 100]\n",
    "gene2 = [350, 100, 50, 200]\n",
    "expression_data = [gene0, gene1, gene2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is this a bad idea?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now with NumPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "a = np.array(expression_data)\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to:\n",
    "\n",
    "* Obtain an *RPKM* expression matrix\n",
    "* Quantile normalize the data\n",
    "\n",
    "using the awesome power of NumPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inside a numpy ndarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_info(a):\n",
    "    print('number of elements:', a.size)\n",
    "    print('number of dimensions:', a.ndim)\n",
    "    print('shape:', a.shape)\n",
    "    print('data type:', a.dtype)\n",
    "    print('strides:', a.strides)\n",
    "    print('flags:')\n",
    "    print(a.flags)\n",
    "    \n",
    "print_info(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(a.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "abytes = a.ravel().view(dtype=np.uint8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_info(abytes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(abytes[:24])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: take the transpose of `a`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_info(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_info(a.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: skipping rows and columns with *slicing*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_info(a.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_info(a.T[::2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_info(a.T[::2, ::2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting a copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[0, 0] = 5\n",
    "print(b)\n",
    "a[0, 0] = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced operations: axis-wise evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expr = np.load('expr.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_info(expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This has the raw read count data. However, each sample gets a different number of reads, so we want to normalize by the *library size*, which is the total number of reads across a column.\n",
    "\n",
    "The `np.sum` function returns the sum of all the elements of an array. With the `axis` argument, you can take the sum *along the given axis*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lib_size = np.sum(expr, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Generate a 10 x 3 array of random numbers. From each row, pick the column containing the number closest to 0.75.\n",
    "\n",
    "*Hint:* use of `np.abs` and `np.argmin` to find the column j that contains the closest element in each row i. The final result should be a vector of integers of shape (10,)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Some applications, such as clustering, are computationally expensive, and wouldn't work without first doing some form of *feature selection*, where we discard most of the data and keep only what we think will be most useful. One simple version is to keep only the genes with the most variance (as these will be more informative than genes that don't vary between patients).\n",
    "\n",
    "- Find the variance across patients of all the genes (rows) in the expression dataset.\n",
    "- Use `np.argsort` to find the location of the 1,500 most variable genes.\n",
    "- Use these indices to produce a shape (1500, 375) matrix containing only the most variable genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced operations: broadcasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to normalize every column by its corresponding library size, we have to *align* the two arrays' axes: each dimension must be either the same size, or one of the arrays must have size 1. Use `np.newaxis` to match the dimensions. But let's first do some simple examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a + 5  # simplest \"broadcasting\": scalar - array operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([1, 2, 3, 4])\n",
    "a + b  # broadcasting: coerce arrays to same shape by repeating as needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([1, 2, 3])\n",
    "a + b  # broadcasting: not just magic!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([[1], [2], [3]])\n",
    "a + b  # broadcasting: shape compatibility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, back to our expression data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(expr.shape)\n",
    "print(lib_size.shape)\n",
    "print(lib_size[np.newaxis, :].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, NumPy will automatically prepend singleton dimensions until the array shapes match or there is an error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.all(expr / lib_size ==\n",
    "       expr / lib_size[np.newaxis, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expr_lib = expr / lib_size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also multiply by $10^6$ in order to keep the numbers on a readable scale (reads per million reads)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expr_lib *= 1e6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, longer genes are more likely to produce reads. So we must normalize by the gene length (in kb) to produce a measure of expression called Reads Per Kilobase per Million reads (RPKM). We start by loading the gene lengths in *bases*. (1 kilobase = 1,000 bases.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_len = np.load('gene-lens.npy')\n",
    "print(gene_len.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: broadcast `expr`, `lib_size` and `gene_len` together to produce RPKM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reminder:\n",
    "\n",
    "$RPKM = \\frac{C}{N \\times 10^{-6} \\times L \\times 10^{-3}} = \\frac{10^9C}{NL}$\n",
    "\n",
    "where $C$ is the raw counts, $N$ is the library size (in reads) and $L$ is the gene length (in bases). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rpkm = ...  # FIX THIS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to admire our handywork, we will use a custom plotting function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "from scipy import stats\n",
    "\n",
    "def plot_col_density(data, xlim=None, ax=None, *args, **kwargs):\n",
    "    \"\"\"Use gaussian smoothing to estimate the density.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data : 2D array\n",
    "        The data to be plotted. Each column will result in a\n",
    "        different KDE line.\n",
    "    xlim : tuple of int, optional\n",
    "        The range of the x-axis. Useful when the output has a\n",
    "        few outlier values, squishing the bulk of the\n",
    "        distribution to a small range.\n",
    "    ax : matplotlib Axes object, optional\n",
    "        The axes on which to plot. If None, new axes will be\n",
    "        created.\n",
    "    *args, **kwargs\n",
    "        Extra arguments to ax.plot.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    ax : matplotlib Axes object\n",
    "        The axes on which the plot was made.\n",
    "    \"\"\"\n",
    "    density_per_col = [stats.kde.gaussian_kde(col) for col in data.T]\n",
    "    if xlim is not None:\n",
    "        m, M = xlim\n",
    "    else:\n",
    "        m, M = np.min(data), np.max(data)\n",
    "    x = np.linspace(m, M, 100)\n",
    "\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "    for density in density_per_col:\n",
    "        ax.plot(x, density(x), *args, **kwargs)\n",
    "    ax.set_xlabel('log-counts')\n",
    "    ax.set_ylabel('frequency')\n",
    "    if xlim is not None:\n",
    "        ax.set_xlim(xlim)\n",
    "    return ax\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_col_density(np.log(expr+1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_col_density(np.log(rpkm + 1), xlim=(0, 6));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should see that the most \"disparate\" data column is now a much better fit with the rest of the data. This should improve downstream scientific analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: 3D broadcasting\n",
    "\n",
    "Below, using broadcasting, produce the array containing the sum of every element in `x` with every element in `y`. That is, produce an array `z` such that `z[i, j, k]` contains either the sum of `x[i]` and `y[j, k]` OR the sum of `y[i, j]` and `x[k]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.random(size=(3, 5))\n",
    "y = np.random.randint(10, size=8)\n",
    "z = ... # FIX THIS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: explicit broadcasting and stride tricks\n",
    "\n",
    "Now, use `np.broadcast_arrays` to generate `xbroad` and `ybroad` that are the same shape as `z` (so that a simple element-wise addition will give `z`). Then use `print_info` on `xbroad` and `ybroad`. Notice anything weird?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stride tricks\n",
    "\n",
    "By manipulating the shape and strides of an array, we can produce a \"virtual\" array much bigger than its memory usage:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def repeat(arr, n):\n",
    "    return np.lib.stride_tricks.as_strided(arr,\n",
    "                                           shape=(n,) + arr.shape,\n",
    "                                           strides=(0,) + arr.strides)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "repeated_row = repeat(np.random.random(size=5), 4)\n",
    "repeated_row"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Be careful, though: some operations, such as `np.copy`, actually materialize the much bigger array!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_info(repeated_row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_info(np.copy(repeated_row))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.random((3, 2)).astype(np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try to answer these without looking at `x`. Then, try them out with the `print_info` function.\n",
    "\n",
    "- What is the shape of `x`?\n",
    "- What are the strides of `x`?\n",
    "- Is `x` C-contiguous, F-contiguous, or neither?\n",
    "\n",
    "Now let `y = repeat(x, 4)`. What is its shape? What are its strides? Is it contiguous?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: `np.lib.stride_tricks.as_strided`\n",
    "\n",
    "Use `as_strided` to produce a sliding-window view of a 1D array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sliding_window(arr, size=2):\n",
    "    \"\"\"Produce an array of sliding window views of `arr`\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    arr : 1D array, shape (N,)\n",
    "        The input array.\n",
    "    size : int, optional\n",
    "        The size of the sliding window.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    arr_slide : 2D array, shape (N - size + 1, size)\n",
    "        The sliding windows of size `size` of `arr`.\n",
    "        \n",
    "    Examples\n",
    "    --------\n",
    "    >>> a = np.array([0, 1, 2, 3])\n",
    "    >>> sliding_window(a, 2)\n",
    "    array([[0, 1],\n",
    "           [1, 2],\n",
    "           [2, 3]])\n",
    "    \"\"\"\n",
    "    return ...  # fix this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test your code here\n",
    "sliding_window(np.arange(8), 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: mean filtering\n",
    "\n",
    "Use `sliding_window` to implement mean filtering, in which every value in an array is replaced by the mean of it and its neighbours. This is a basic operation in signal processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_filter(signal, window_size=3):\n",
    "    \"\"\"Apply a mean filter to the input with the desired window size.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    signal : 1D array, shape (M,)\n",
    "        The input signal.\n",
    "    window_size : int, optional\n",
    "        The size of the window along which to compute the mean.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    filtered : 1D array, shape (M - window_size + 1,)\n",
    "        The filtered signal.\n",
    "    \"\"\"\n",
    "    return ...  # FIX THIS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test your function, we will use the example of a *difference filter*, which finds the location of changes in a signal using *convolution*. When the signal is perfectly noiseless, it works great:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = np.zeros(100, np.float)\n",
    "signal[30:60] = 1\n",
    "\n",
    "diff = np.array([1, 0, -1])\n",
    "from scipy import ndimage as ndi\n",
    "dsignal = ndi.convolve(signal, diff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2)\n",
    "ax[0].plot(signal)\n",
    "ax[0].set_title('signal')\n",
    "ax[1].plot(dsignal)\n",
    "ax[1].set_title('change')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, if the signal is corrupted by noise, a standard difference filter convolution doesn't work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "signal_noisy = signal + np.random.normal(0, 0.3, size=signal.shape)\n",
    "dsignal_noisy = ndi.convolve(signal_noisy, diff)\n",
    "\n",
    "fig, ax = plt.subplots(1, 2)\n",
    "ax[0].plot(signal_noisy)\n",
    "ax[0].set_title('signal')\n",
    "ax[1].plot(dsignal_noisy)\n",
    "ax[1].set_title('change')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Try mean filtering with different window sizes to see whether the change signal becomes more apparent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: padding\n",
    "\n",
    "What is the shape of your mean-filtered signal?\n",
    "\n",
    "...\n",
    "\n",
    "Oops! We've shortened the signal, which means that our indices have changed: `signal_filtered[i]` does not correspond to the signal around `signal[i]`.\n",
    "\n",
    "Use `np.pad` to add some \"fake\" data around `signal` before filtering, so that the filtered result has the same shape as the input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Gaussian filtering\n",
    "\n",
    "It turns out that mean filtering is not the \"optimal\" way to recover your signal, assuming certain properties of the noise. For that, we use *Gaussian* filtering, which uses a *weighted* mean of the sliding window elements. The weights are given by the famous Gaussian bell-shaped distribution. For example, here are the weights for a window size of 17 for a particular sigma:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weight = np.exp(-(np.arange(-8, 9) / (8/3))**2)\n",
    "weight /= np.sum(weight)  ## ensure overall intensity of signal doesn't change\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a function that uses sliding windows, broadcasting, and axis-wise operations to compute the Gaussian filter of a signal for a given window size. (You should also pad your input.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fancy indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can index arrays with slicing, but also with boolean arrays (including broadcasting!), integer arrays, and individual indices along multiple dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = np.array([0, 5, 99])\n",
    "selector = np.random.randint(0, 3, size=(3, 4))\n",
    "print(selector)\n",
    "print(values[selector])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relabeled = values[selector]\n",
    "has_large_cols = np.any(relabeled > 10, axis=1)\n",
    "print(relabeled[has_large_cols])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Use boolean indexing and broadcasting to select the columns of `relabeled` that do not contain 99."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: quantile normalization\n",
    "\n",
    "Quantile Normalization (https://en.wikipedia.org/wiki/Quantile_normalization) is a method to align distributions. Here we implement it using NumPy axis-wise operations and fancy indexing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qnorm(X):\n",
    "    \"\"\"Quantile normalize an input matrix.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X : 2D array of float, shape (M, N)\n",
    "        The input data, with each column being a\n",
    "        distribution to normalize.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    Xn : 2D array of float, shape (M, N)\n",
    "        The normalized data.\n",
    "    \"\"\"\n",
    "    ranks = \n",
    "    return Xn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logexpr = np.log(expr + 1)\n",
    "logrpkm = np.log(rpkm + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logexprn = qnorm(logexpr)\n",
    "logrpkmn = qnorm(logrpkm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_col_density(logexprn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_col_density(logrpkmn, xlim=(0, 0.25))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fancy indexing along multiple dimensions\n",
    "\n",
    "Combining fancy indexing and slicing selects entire rows/columns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relabeled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relabeled[[1, 1, 2], :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relabeled[:, [1, 3, 1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To select individual elements for a new array shape, we must use as many fancy indices as the array has dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selector_rows = [[0, 0],\n",
    "                 [1, 2]]\n",
    "selector_cols = [[0, 3],\n",
    "                 [1, 2]]\n",
    "\n",
    "arr = np.arange(12).reshape((3, 4))\n",
    "print(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(arr[selector_rows, selector_cols])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to think about this is:\n",
    "- make a \"coordinate array\", of the shape that you want plus one more axis, to hold the coordinates of each point (see below),\n",
    "- transpose that final axis to the front, and\n",
    "- convert to tuple\n",
    "\n",
    "For the above example, perhaps you find this \"notation\", with the individual coordinates in the final axis, more intuitive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selector_t = [[ [0, 0], [0, 3] ],\n",
    "              [ [1, 1], [2, 2] ]] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, we want element (0, 0) in the top left corner, (0, 3) in the top right, (1, 1) in the bottom left, and (2, 3) in the bottom right. However, *rows and columns* must be in the first dimension and presented as a tuple to index the original array. So, to use this notation, we can use np.transpose and cast the result to `tuple`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selector = tuple(np.transpose(selector_t, (2, 0, 1)))\n",
    "print(selector[0], selector[1], sep='\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For technical reasons that one might grasp for fleeting moments, the \"tuple of index arrays\" format is most consistent with other forms of multi-dimensional indexing in NumPy. It is a widespread convention (see e.g. `scipy.ndimage.map_coordinates`), so it's worth practicing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "What happens when you make `selector_col`:\n",
    "- a single number?\n",
    "- a 1D array with two elements?\n",
    "- a 2D array of shape (1, 2)?\n",
    "- a 1D array with three elements?\n",
    "\n",
    "Repeat similar experiments with `selector_row`.\n",
    "\n",
    "Does this remind you of any other NumPy feature we may have seen?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced exercise: Jack's dilemma\n",
    "\n",
    "(If time permits.)\n",
    "\n",
    "```email\n",
    "Date: Wed, 16 Jul 2008 16:45:37 -0500\n",
    "From: Jack Cook\n",
    "To: <numpy-discussion@scipy.org>\n",
    "Subject: Numpy Advanced Indexing Question\n",
    "\n",
    "Greetings,\n",
    "\n",
    "I have an I,J,K 3D volume of amplitude values at regularly sampled\n",
    "time intervals. I have an I,J 2D slice which contains a time (K)\n",
    "value at each I, J location. What I would like to do is extract a\n",
    "subvolume at a constant +/- K window around the slice. Is there an\n",
    "easy way to do this using advanced indexing or some other method?\n",
    "Thanks in advanced for your help.\n",
    "\n",
    "-- Jack\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In the command-line, run `python tools/generate-volume.py`\n",
    "\n",
    "data = np.load('geo.npz')  # npz file with multiple arrays, keyed like dict\n",
    "strata, density = data['strata'], data['density']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row = 256\n",
    "\n",
    "fig, (ax0, ax1) = plt.subplots(1, 2, sharey=True, figsize=(7.2, 3.6))\n",
    "ax0.imshow(strata[:, row, :], cmap='gray')\n",
    "ax0.set_ylabel('depth')\n",
    "ax0.set_xlabel('x')\n",
    "ax0.set_title(f'strata at row {row}')\n",
    "ax1.imshow(density[:, 256, :], cmap='magma')\n",
    "ax1.set_xlabel('x')\n",
    "ax1.set_title(f'fossil density at row {row}') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Exercise\n",
    "\n",
    "We want to quantify the apparent extinction event near a particular stratum. The strata, however, are not perfectly horizontal, and there is a large break in the rock, too. Therefore, we want to align the strata by using the solution to Jack's dilemma (which you must solve for him).\n",
    "\n",
    "Using a window size of 120 (so the half-width is K=60):\n",
    "- find the stratum of maximum intensity along the depth axis. This is a 2D slice of integers measuring the depth of the maximum intensity stratum at each (row, column) coordinate.\n",
    "- using broadcasting, create a volume of shape (2K, Nrow, Ncol) containing the depth coordinate at each (row, column) of your desired window.\n",
    "- create matching row and column index volumes to perform fancy indexing of the fossil density array.\n",
    "- extract the subvolume of fossil density around the stratum using fancy indexing\n",
    "- compute the mean fossil density at each depth.\n",
    "\n",
    "**-->Advice<--**: play around with a much tinier subset of the volume, say, a subset of size (10, 15, 20), and window size of 3. Evaluating the wrong expression can result in a giant dataset that blows up your memory and crashes your computer. (Yay Science!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear algebra with NumPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since version 3.5, Python supports the matrix multiplication operator, denoted by the `@` symbol. This makes it a new powerhouse for linear algebra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.array([[0, 1],\n",
    "              [1, 1],\n",
    "              [1, 0]])\n",
    "v = np.array([9, 2])\n",
    "\n",
    "print(M @ v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(M.T @ M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider the rotation matrix\n",
    "\n",
    "$\n",
    "R = \\begin{bmatrix}\n",
    "  \\cos \\theta &  -\\sin \\theta & 0 \\\\\n",
    "  \\sin \\theta & \\cos \\theta & 0\\\\\n",
    "  0 & 0 & 1\\\\\n",
    "\\end{bmatrix}\n",
    "$\n",
    "\n",
    "When $R$ is multiplied with a 3-dimensional column-vector $p =\n",
    "\\left[ x\\, y\\, z \\right]^T$, the resulting vector $R p$ is rotated\n",
    "by $\\theta$ degrees around the z-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_deg = 45\n",
    "theta = np.deg2rad(theta_deg)\n",
    "\n",
    "c = np.cos(theta)\n",
    "s = np.sin(theta)\n",
    "\n",
    "R = np.array([[c, -s, 0],\n",
    "              [s,  c, 0],\n",
    "              [0,  0, 1]]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_axis = np.array([1, 0, 0]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R @ x_axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rotated = R @ x_axis\n",
    "R @ rotated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: rotation matrices\n",
    "\n",
    "If you know some linear algebra, try to answer the following questions about $R$ in your head, before trying them out in Python. If you're totally stuck, just give them a go!\n",
    "\n",
    "- What does the matrix $S = RR$ do? (Remember, `S = R @ R` in code.)\n",
    "- What does $R^4 = RRRR$ do?\n",
    "- The inverse $R^{-1}$ of $R$ is defined as the matrix such that $R^{-1}R = I$, where $I$ is the identity matrix. What does $Q = R^{-1}$ do? (Note: NumPy provides matrix inverse computation in `np.linalg.inv`.)\n",
    "- What does $R$ do to the vector [0, 0, 1]?\n",
    "- What does that make the vector [0, 0, 1]?\n",
    "- Verify this with some function in `np.linalg`. (Look at the online documentation for this module.)\n",
    "- What does $R^8$ do? What does that make $R^8$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: revisiting Gaussian filtering\n",
    "\n",
    "Above, we performed Gaussian filtering with a broadcast product and a axis-wise sum. Use a matrix multiplication instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neural network with NumPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspired by https://iamtrask.github.io/2015/07/12/basic-python-network/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(128) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.random.random(size=(n, 3)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Make a vector y of length n that is 0 at position i if `X[i, 0] + X[i, 2] < 1`, and 1 otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Hint:* How do I get the 0th column of X?\n",
    "\n",
    "*Hint 2:* How can I convert an array of (True/False) to an array of float?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[2, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = (X[:, 0] + X[:, 2] > 1).astype(float) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X[:, 2], X[:, 0], color=plt.cm.viridis(y)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xb = np.column_stack((X, np.ones(X.shape[0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = 2 * np.random.random(size=Xb.shape[1]) - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    return 1 / (1 + np.exp(-x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = sigmoid(Xb @ W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(y, pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dsigmoid(x):\n",
    "    return x * (1 - x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train(X, y, W, n_iter=100):\n",
    "    W_history = np.empty((n_iter, W.size), dtype=float)\n",
    "    error_history = np.empty(n_iter, dtype=float)\n",
    "    W_out = np.copy(W)\n",
    "    for i in range(n_iter):\n",
    "        W_history[i] = W_out\n",
    "        predictions = sigmoid(X @ W_out)\n",
    "        derror = y - predictions\n",
    "        gradient = derror * dsigmoid(predictions)\n",
    "        update = X.T @ gradient\n",
    "        error_history[i] = np.sum(np.abs(error))\n",
    "        W_out += update\n",
    "    return W_out, W_history, error_history"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Wout, Whist, errhist = train(Xb, y, W, 10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.semilogx(errhist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mgrid[0:1:0.3, 0:1:0.3] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decision_space = np.mgrid[0:1:0.01, 0:1:1, 0:1:0.01, 1:2:1] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decision_space.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction_image = sigmoid(decision_space.T @ Wout).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction_image.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction_squeezed = np.squeeze(prediction_image) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(prediction_squeezed, cmap='magma', origin='lower',\n",
    "          extent=[0, 1, 0, 1])\n",
    "ax.scatter(X[:, 2], X[:, 0], color=plt.cm.viridis(y)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Why do some points appear misclassified?\n",
    "\n",
    "*Hint:* Look at `Wout`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References and sources\n",
    "\n",
    "In addition to original content, these notes use materials from the following sources:\n",
    "\n",
    "- [100 NumPy Exercises](http://www.labri.fr/perso/nrougier/teaching/numpy.100/), compiled by Nicolas Rougier, used here under the terms of the MIT License.\n",
    "- [SciPy Lecture Notes](http://www.scipy-lectures.org/), compiled by Gaël Varoquaux, Emmanuelle Gouillart, and Olav Vahtras, used under the terms of the CC-BY license.\n",
    "- [Elegant SciPy](http://elegant-scipy.org), written by Juan Nunez-Iglesias, Stéfan van der Walt, and Harriet Dashnow, used with permission from O'Reilly.\n",
    "- [iamtrask blog](https://iamtrask.github.io/2015/07/12/basic-python-network/), for a neural network from scratch that inspired this one.\n",
    "\n",
    "They are all freely available online and worth checking out in their entirety!"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}