numpy_advanced.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end",
     "slide_type": "subslide"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# NumPy\n",
    "\n",
    "## A structure for efficient numerical computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end",
     "slide_type": "subslide"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"images/overview.png\" width=\"80%\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Slicing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[9, 7, 5, 3, 1]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = list(range(10))\n",
    "x[2]\n",
    "x[-2]\n",
    "x[1:5]\n",
    "x[3:8]\n",
    "x[::2]\n",
    "x[::-2]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end",
     "slide_type": "subslide"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Num-what?\n",
    "\n",
    "This talk discusses some of the more advanced NumPy features. If you’ve\n",
    "never seen NumPy before, you may have more fun doing the tutorial at\n",
    "\n",
    "http://mentat.za.net/numpy/intro/intro.html\n",
    "\n",
    "or studying the course notes at\n",
    "\n",
    "http://scipy-lectures.github.io\n",
    "\n",
    "You can always catch up by reading:\n",
    "\n",
    ">  **The NumPy array: a structure for efficient numerical computation**. Stéfan\n",
    ">  van der Walt, S. Chris Colbert and Gaël Varoquaux. In IEEE Computing in\n",
    ">  Science Engineering, March/April 2011.\n",
    "\n",
    "Nicolas Rougier also has a beautiful tutorial at\n",
    "\n",
    "https://github.com/rougier/numpy-tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np  # We always use this convention,\n",
    "                    # also in the documentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The NumPy ndarray\n",
    "\n",
    "### Revision: Structure of an array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2],\n",
       "       [3, 4, 5],\n",
       "       [6, 7, 8]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.array([[ 0,  1,  2],\n",
    "              [ 3,  4,  5],\n",
    "              [ 6,  7,  8]])\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end",
     "slide_type": "subslide"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"images/ndarray_struct.svg.png\" width=\"50%\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "<img src=\"images/array_memory_presentation.png\" width=\"80%\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## A homogeneous container"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<memory at 0x7fbbef344210>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Data is just a pointer to bytes in memory:\n",
    "\n",
    "<img src=\"images/array_dtype.png\" width=\"60%\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of elements: 9\n",
      "number of dimensions: 2\n",
      "shape: (3, 3)\n",
      "data type: int64\n",
      "strides: (24, 8)\n",
      "flags:\n",
      "  C_CONTIGUOUS : True\n",
      "  F_CONTIGUOUS : False\n",
      "  OWNDATA : True\n",
      "  WRITEABLE : True\n",
      "  ALIGNED : True\n",
      "  WRITEBACKIFCOPY : False\n",
      "  UPDATEIFCOPY : False\n",
      "\n"
     ]
    }
   ],
   "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(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "<img src=\"images/array_flat_extended.png\" width=\"80%\"/>\n",
    "<br/>\n",
    "<div style=\"text-align: center;\">\n",
    "<img style=\"display: inline;\" src=\"images/array_reshape.png\" width=\"40%\"/>\n",
    "<img style=\"display: inline;\" src=\"images/array_reshape_nd.png\" width=\"40%\"/>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 2, 3, 3)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.random.random((3, 2, 3, 3))\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.ndim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('float64')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Common types in include ``int``, ``float``, ``bool``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1,  0,  1])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array ([-1, 0, 1], dtype=int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.,  0.,  1.])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([-1, 0, 1], dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True, False,  True])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([-1 , 0, 1], dtype=bool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Each item in the array has to have the same type (occupy a fixed nr of bytes\n",
    "in memory), but that does not mean a type has to consist of a single item:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([(0,  True), (1, False)], dtype=[('value', '<i8'), ('status', '?')])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dt = np.dtype([('value', np.int), ('status', np.bool)])\n",
    "x = np.array([(0, True), (1, False)], dtype=dt)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "This is called a structured array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Strides"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([(0,  True), (1, False)], dtype=[('value', '<i8'), ('status', '?')])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3],\n",
       "       [ 4,  5,  6,  7],\n",
       "       [ 8,  9, 10, 11]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.arange(12).reshape((3, 4))\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('int64')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.dtype.itemsize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(32, 8)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.strides"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We see ``(4 * itemsize, itemsize)`` or ``(skip_bytes_row, skip_bytes_col)``."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "<img src=\"images/array_memory_dtype.png\" width=\"100%\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  2,\n",
       "        0,  0,  0,  0,  0,  0,  0,  3,  0,  0,  0,  0,  0,  0,  0,  4,  0,\n",
       "        0,  0,  0,  0,  0,  0,  5,  0,  0,  0,  0,  0,  0,  0,  6,  0,  0,\n",
       "        0,  0,  0,  0,  0,  7,  0,  0,  0,  0,  0,  0,  0,  8,  0,  0,  0,\n",
       "        0,  0,  0,  0,  9,  0,  0,  0,  0,  0,  0,  0, 10,  0,  0,  0,  0,\n",
       "        0,  0,  0, 11,  0,  0,  0,  0,  0,  0,  0], dtype=uint8)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xbytes = x.ravel().view(dtype=np.uint8)\n",
    "xbytes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0, 0, 0, 0, 0, 0, 0, 0], dtype=uint8),\n",
       " array([1, 0, 0, 0, 0, 0, 0, 0], dtype=uint8),\n",
       " array([2, 0, 0, 0, 0, 0, 0, 0], dtype=uint8))"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xbytes[:8], xbytes[8:16], xbytes[16:24]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end",
     "slide_type": "subslide"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example: Taking the transpose\n",
    "\n",
    "We can take the transpose of the array *without* copying any data by simply manipulating ``shape`` and ``strides``.  What should they be?\n",
    "\n",
    "<img src=\"images/array_transpose_question.png\" width=\"100%\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of elements: 12\n",
      "number of dimensions: 2\n",
      "shape: (4, 3)\n",
      "data type: int64\n",
      "strides: (8, 32)\n",
      "flags:\n",
      "  C_CONTIGUOUS : False\n",
      "  F_CONTIGUOUS : True\n",
      "  OWNDATA : False\n",
      "  WRITEABLE : True\n",
      "  ALIGNED : True\n",
      "  WRITEBACKIFCOPY : False\n",
      "  UPDATEIFCOPY : False\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print_info(x.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Flags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of elements: 12\n",
      "number of dimensions: 2\n",
      "shape: (3, 4)\n",
      "data type: int64\n",
      "strides: (32, 8)\n",
      "flags:\n",
      "  C_CONTIGUOUS : True\n",
      "  F_CONTIGUOUS : False\n",
      "  OWNDATA : False\n",
      "  WRITEABLE : True\n",
      "  ALIGNED : True\n",
      "  WRITEBACKIFCOPY : False\n",
      "  UPDATEIFCOPY : False\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print_info(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  C_CONTIGUOUS : False\n",
       "  F_CONTIGUOUS : False\n",
       "  OWNDATA : False\n",
       "  WRITEABLE : True\n",
       "  ALIGNED : True\n",
       "  WRITEBACKIFCOPY : False\n",
       "  UPDATEIFCOPY : False"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = x[::2]\n",
    "\n",
    "z.flags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Copies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 5  1  2  3]\n",
      " [ 4  5  6  7]\n",
      " [ 8  9 10 11]]\n"
     ]
    }
   ],
   "source": [
    "y = x\n",
    "\n",
    "y[0, 0] = 5\n",
    "print(y)\n",
    "y[0, 0] = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[100,   1,   2,   3],\n",
       "       [  4,   5,   6,   7],\n",
       "       [  8,   9,  10,  11]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.flags.owndata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## x.func vs. np.func(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced operations: axis-wise evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "expr = np.load('expr.npz')['expr']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[     0,      0,      0, ...,      0,     18,      0],\n",
       "       [    12,      1,      0, ...,      2,      3,      4],\n",
       "       [ 64134,  52319,  48016, ...,  80938, 105589,  92008],\n",
       "       ...,\n",
       "       [   401,   2518,   1511, ...,   2260,   3574,   3936],\n",
       "       [   458,   2164,    540, ...,    977,   1302,    681],\n",
       "       [     0,      0,      0, ...,     10,      0,      0]],\n",
       "      dtype=uint32)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of elements: 7687500\n",
      "number of dimensions: 2\n",
      "shape: (20500, 375)\n",
      "data type: uint32\n",
      "strides: (4, 82000)\n",
      "flags:\n",
      "  C_CONTIGUOUS : False\n",
      "  F_CONTIGUOUS : True\n",
      "  OWNDATA : False\n",
      "  WRITEABLE : True\n",
      "  ALIGNED : True\n",
      "  WRITEBACKIFCOPY : False\n",
      "  UPDATEIFCOPY : False\n",
      "\n"
     ]
    }
   ],
   "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": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(375,)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lib_size = np.sum(expr, axis=0)\n",
    "lib_size.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20500,)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(expr, axis=1).shape"
   ]
  },
  {
   "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": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Structured arrays\n",
    "\n",
    "Like we saw earlier, each item in an array has the same type, but that does not\n",
    "mean a type has to consist of a single item:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dt = np.dtype([('value', np.int), ('status', np.bool)])\n",
    "np.array([(0, True), (1, False)], dtype=dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "This is called a structured array, and is accessed like a dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = np.array([(0, True), (1, False)], dtype=dt)\n",
    "x['value']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x['status']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Reading structured data from file\n",
    "\n",
    "<img src=\"images/radar_struct.png\" width=\"100%\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Reading this kind of data can be troublesome:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "```matlab\n",
    "while ((count > 0) && (n <= NumPoints))\n",
    "  % get time - I8 [ ms ]\n",
    "  [lw, count] = fread(fid, 1, 'uint32');\n",
    "  if (count > 0) % then carry on\n",
    "  uw = fread(fid, 1, 'int32');\n",
    "  t(1, n) = (lw + uw * 2^32) / 1000;\n",
    "\n",
    "  % get number of bytes of data\n",
    "  numbytes = fread(fid, 1, 'uint32');\n",
    "\n",
    "  % read sMEASUREMENTPOSITIONINFO (11 bytes)\n",
    "  m(1, n) = fread(fid, 1, 'float32'); % az [ rad ]\n",
    "  m(2, n) = fread(fid, 1, 'float32'); % el [ rad ]\n",
    "  m(3, n) = fread(fid, 1, 'uint8');   % region type\n",
    "  m(4, n) = fread(fid, 1, 'uint16');  % region ID\n",
    "  g(1, n) = fread(fid, 1, 'uint8');\n",
    "\n",
    "  numsamples = ( numbytes -12)/2; % 2 byte int\n",
    "  ...\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using structured arrays:\n",
    "\n",
    "```python\n",
    "dt = np.dtype([('time', np.uint64),\n",
    "               ('size', np.uint32),\n",
    "               ('position', [('az', np.float32),\n",
    "                             ('el', np.float32),\n",
    "                             ('region_type', np.uint8),\n",
    "                             ('region_ID', np.uint16)]),\n",
    "               ('gain', np.uint8),\n",
    "               ('samples', (np.int16, 2048))])\n",
    "\n",
    "data = np.fromfile(f, dtype=dt)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We can then access this structured array as before:\n",
    "\n",
    "```python\n",
    "data['position']['az']\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Broadcasting\n",
    "\n",
    "Combining of differently shaped arrays without creating large intermediate\n",
    "arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = np.arange(4)\n",
    "print(x)\n",
    "print(x + 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "![Broadcasting in 1D](images/broadcast_1D.png)\n",
    "\n",
    "See docstring of ``np.doc.broadcasting`` for more detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Broadcasting (2D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = np.arange(12).reshape((3, 4))\n",
    "b = np.array([1, 2, 3])[:, np.newaxis]\n",
    "\n",
    "print(a)\n",
    "print()\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"images/broadcast_2D.png\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Broadcasting (3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = np.zeros((3, 5))\n",
    "y = np.zeros(8)\n",
    "z = x[..., np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "```python\n",
    "x = np.zeros((3, 5))\n",
    "y = np.zeros(8)\n",
    "z = x[..., np.newaxis]\n",
    "```\n",
    "\n",
    "```\n",
    "X shape = (3, 5, 1)\n",
    "Y shape = (      8)\n",
    "Z shape = (3, 5, 8)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"images/broadcast_3D.png\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Broadcasting rules\n",
    "\n",
    "The broadcasting rules are straightforward: compare dimensions, starting from\n",
    "the last. Match when either dimension is one or ``None``, or if dimensions are\n",
    "equal:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "```\n",
    "Scalar    2D           3D           Bad\n",
    "\n",
    "( ,)     (3, 4)     (3, 5, 1)    (3, 5, 2)\n",
    "(3,)     (3, 1)     (      8)    (      8)\n",
    "----     ------     ---------    ---------\n",
    "(3,)     (3, 4)     (3, 5, 8)       XXX\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Explicit broadcasting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = np.zeros((3, 5, 1))\n",
    "y = np.zeros((3, 5, 8))\n",
    "\n",
    "xx, yy = np.broadcast_arrays(x, y)\n",
    "\n",
    "xx.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "```python\n",
    "np.broadcast_arrays([1, 2, 3],\n",
    "                    [[1], [2], [3]])\n",
    "```\n",
    "\n",
    "```\n",
    "[[1 2 3]   [[1 1 1]\n",
    " [1 2 3]    [2 2 2]\n",
    " [1 2 3]]   [3 3 3]]\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "np.broadcast_to(x, (3,5,8)).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fancy indexing\n",
    "\n",
    "An ndarray can be indexed in two ways:\n",
    "\n",
    " - Using slices and scalars\n",
    " - Using other ndarrays (\"fancy indexing\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x = np.arange(9).reshape((3, 3))\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x[:, [1, 1, 2]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Which is equivalent to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "np.array((x[:, 1], x[:, 1], x[:, 2])).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Output shape of an indexing op\n",
    "\n",
    "1. Broadcast all index arrays against one another.\n",
    "2. Use the dimensions of slices as-is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start with the following input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "What would happen if we do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "idx0 = np.array([[0, 1],\n",
    "                 [1, 2]])  # row indices\n",
    "\n",
    "idx1 = np.array([[0, 1]])  # column indices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "```python\n",
    "x[idx0, idx1]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Output shape of indexing op (cont'd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "idx0.shape, idx1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a, b = np.broadcast_arrays(idx0, idx1)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x[idx0, idx1]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Output shape of an indexing op (cont'd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = np.random.random((15, 12, 16, 3))\n",
    "\n",
    "index_one = np.array([[0, 1], [2, 3], [4, 5]])\n",
    "index_two = np.array([[0, 1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(index_one.shape)\n",
    "print(index_two.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Predict the output shape of"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "```python\n",
    "x[5:10, index_one, :, index_two]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Warning! When mixing slicing and fancy indexing, the\n",
    "*order* of the output dimensions is less easy to predict.\n",
    "Play it safe and **don't mix the two!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Output shape of an indexing op (cont'd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = np.random.random((15, 12, 16, 3))\n",
    "index_one = np.array([[0, 1], [2, 3], [4, 5]])\n",
    "index_two = np.array([[0, 1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(index_one.shape)\n",
    "print(index_two.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Broadcast ``index1`` against ``index2``:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "```\n",
    "(3, 2)  # shape of index1\n",
    "(1, 2)  # shape of index2\n",
    "------\n",
    "(3, 2)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The shape of ``x[5:10, index_one, :, index_two]`` is\n",
    "\n",
    "``(3, 2, 5, 16)``"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Jack's dilemma\n",
    "\n",
    "Indexing and broadcasting are intertwined, as we’ll see in the following\n",
    "example. One of my favourites from the NumPy mailing list:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "```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",
    "\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end",
     "slide_type": "subslide"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"images/jack.png\" width=\"90%\" alt=\"Jack's cube\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Test setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ni, nj, nk = (10, 15, 20)  # dimensions\n",
    "half_width = 3  # take 7 elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Make a fake data block such that ``block[i, j, k] == k`` for all i, j, k."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "block = np.empty((ni, nj, nk) , dtype=int)\n",
    "block[:] = np.arange(nk)[np.newaxis, np.newaxis, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Pick out a random fake horizon in ``k``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.RandomState(0)\n",
    "\n",
    "k = rng.randint(5, 15, size=(ni, nj))\n",
    "k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solution\n",
    "\n",
    "These two indices assure that we take a slice at each (i, j) position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "idx_i = np.arange(ni)[:, np.newaxis, np.newaxis]\n",
    "idx_j = np.arange(nj)[np.newaxis, :, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "This is the substantitive part that picks out the window:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "idx_k = k[:, :, np.newaxis] + np.arange(-half_width, half_width + 1)\n",
    "slices = block[idx_i, idx_j, idx_k]  # Slice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Apply the broadcasting rules:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "subslide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "```\n",
    "(ni,  1, 1                 )  # idx_i\n",
    "(1,  nj, 1                 )  # idx_j\n",
    "(ni, nj, 2 * half_width + 1)  # idx_k\n",
    "-\n",
    "(ni, nj, 7)  <-- this is what we wanted!\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solution verification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "slices = block[idx_i, idx_j, idx_k]\n",
    "slices.shape\n",
    "(10, 15, 7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Now verify that our window is centered on ``k`` everywhere:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "slices[:, :, 3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "np.all(slices[:, :, 3] == k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The \\_\\_array\\_interface\\_\\_\n",
    "\n",
    "Any object that exposes a suitable dictionary named\n",
    "``__array_interface__`` may be converted to a NumPy array. This is\n",
    "handy for exchanging data with external libraries. The array interface\n",
    "has the following important keys (see\n",
    "http://docs.scipy.org/doc/numpy/reference/arrays.interface):\n",
    "\n",
    " - **shape**\n",
    " - **typestr**\n",
    " - **data**: (20495857, True); 2-tuple—pointer to data and boolean to\n",
    "indicate whether memory is read-only\n",
    " - **strides**\n",
    " - **version**: 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Stride Tricks: Repeating sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "data = np.array([1, 2, 3, 4])\n",
    "print(data.shape, data.strides)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N = 400000\n",
    "repeat = np.lib.stride_tricks.as_strided(data, shape=(N, 4), strides=(0, 8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "repeat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Did you know?\n",
    "\n",
    "NumPy implements a Diophantine (integer) [equation solver](https://github.com/numpy/numpy/blob/master/numpy/core/src/private/mem_overlap.c) for `np.shares_memory` to calculate whether two arrays possibly overlap in memory.\n",
    "\n",
    "Q: Can you imagine why this is a tricky problem to solve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "sum(a[i] * x[i] for i in range(n)) == b\n",
    "\n",
    "a[i] > 0, 0 <= x[i] <= ub[i]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You can even access this solver from Python.  For example, given that a year has 365 days, how do we fit in months of 28, 30, and 31 days?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.core.multiarray_tests import solve_diophantine\n",
    "\n",
    "month_days = (28, 30, 31)\n",
    "limits = (20, 20, 20)\n",
    "\n",
    "solve_diophantine(month_days, limits, 365)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(2 * 28) + (1 * 30) + (9 * 31) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solve_diophantine(month_days, (1, 20, 20), 365)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(1 * 28) + (4 * 30) + (7 * 31)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Questions, discussion and exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "---\n",
    "<div style=\"height: 200px;\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%reload_ext load_style\n",
    "%load_style css/notebook.css"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}