02-array-programming.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Array programming in NumPy\n",
    "\n",
    "While many of the array operations that are familiar to you from Python lists will also work with NumPy arrays, an important difference between NumPy arrays and Python lists is how we _should_ operate on them.  Specifically, to use NumPy arrays idiomatically, you sometimes need to think in terms of programming with arrays.  This means that some operations will have results that you don't expect if you're thinking of how Python lists work.\n",
    "\n",
    "For example, consider the `max` function in Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So `max` returns its maximum argument.  If its arguments are lists, it will return the list with the largest first element (if there is a tie, it will return the list from the tie with the largest second element, and so on).  Try it out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max([1,1,5,9],[1,2,4,6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The similarly-named NumPy `maximum` function has a slightly different interface and subtly-different behavior.  While the built-in `max` function will take the maximum from any number of arguments, NumPy's `maximum` takes two.  Furthermore, its result is slightly different, as you can see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.maximum([1,1,5,9],[1,2,4,6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Notice that we invoked `np.maximum` on two Python lists, but it would behave similarly on two NumPy arrays.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.maximum(np.array([1,1,5,9]),np.array([1,2,4,6]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operations that work on entire arrays and can change an element at a time (rather than that work on entire arrays but only  as atomic objects) are extremely important in NumPy.\n",
    "\n",
    "These elementwise operations are useful and they're also far more efficient than code you'd write in Python.  To see why, write a Python function to take the elementwise maximum of two arrays and compare its performance to `np.maximum`.   First, let's set up some test data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0xDA7ABA5E)\n",
    "ls1 = list(np.random.randint(256, size=5000))\n",
    "ls2 = list(np.random.randint(256, size=5000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "def nmax(la, lb):\n",
    "    # FIXME: replace this function's body with a function that computes \n",
    "    # and returns the elementwise maximum of a and b.  Do not use NumPy.\n",
    "    return max(la,lb)\n",
    "\n",
    "for i in range(10000):\n",
    "    nmax(ls1, ls2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get stuck on the previous cell, [see the hint](hints/hint02.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "def pmax(a, b):\n",
    "    return np.maximum(a, b)\n",
    "\n",
    "for i in range(10000):\n",
    "    pmax(ls1, ls2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see another example of array programming in the next two cells.  Assuming `ls` is a Python list and `na` is a NumPy array, how is `ls * 3`  different from `na * 3`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[1,2,5,9] * 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array([1,2,5,9]) * 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we multiply a list by a list or an array by an array?  Try the next two cells to find out; we'll switch things up and do NumPy first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Array-array multiplication\n",
    "np.array([1,2,5,9]) * np.array([1,2,5,9])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[1,2,5,9] * [1,2,5,9]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write some Python code in the next cell that returns the same result as the NumPy code in \n",
    "the cell beginning with `# Array-array multiplication`.\n",
    "\n",
    "(If you need a hint, the [hint for elementwise maximum](hints/hint02.ipynb) may be a good place to start!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ewprod(l1, l2):\n",
    "    # replace this function with one that computes \n",
    "    # the elementwise product of l1 and l2\n",
    "    \n",
    "    return l1\n",
    "\n",
    "ewprod([1,2,5,9], [1,2,5,9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll get a quick timing to see the performance difference between these approaches.  First, though, we'll make NumPy arrays from the random Python lists we used to test elementwise maximum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "na1 = np.array(ls1)\n",
    "na2 = np.array(ls2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "acc = None\n",
    "\n",
    "for i in range(10000):\n",
    "    acc = na1 * na2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "acc = None\n",
    "\n",
    "for i in range(10000):\n",
    "    acc = ewprod(ls1, ls2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As Python programmers, our instinct is to perform elementwise transformations on arrays with a `for` loop (in place) or with a list comprehension (making a copy).  The array programming features in NumPy enable us to perform these kinds of transformations _concisely_ (with less code) and _efficiently_ (with less time).  They are faster for several reasons: \n",
    "\n",
    "1. many of these library routines are implemented in native code (a `for` loop in C is often faster than a `for` loop in the Python interpreter),\n",
    "2. some of these routines can use optimized primitives or hardware acceleration (depending on how you've installed NumPy), and\n",
    "3. many of these routines can modify an argument in-place -- which concise pure-Python solutions cannot -- thus saving the cost of allocating a new array to store the result (and reducing memory pressure, memory fragmentation, and the workload for the garbage collector).\n",
    "\n",
    "We'll see how to take advantage of modifying arrays in place in  the next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1,2,3])\n",
    "b = np.array([3,4,5])\n",
    "c = np.array([7,8,9])\n",
    "np.maximum(a, b, c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that we said `maximum` takes the maximum of up to two arrays.  But in the previous cell we've given it three arguments and gotten a confusing result -- it seems to have ignored `c` and just taken the maximum of `a` and `b`.  What's going on?  Let's inspect `a`, `b`, and `c` and see if we can figure it out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(a)\n",
    "print(b)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the three-argument version of `maximum`, NumPy places the output in a preexisting array.  We can see that the contents of `c` are changed but its identity is not:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1,2,3])\n",
    "b = np.array([3,4,5])\n",
    "c = np.array([7,8,9])\n",
    "\n",
    "before = id(c)\n",
    "\n",
    "np.maximum(a, b, c)\n",
    "\n",
    "after = id(c)\n",
    "\n",
    "print(c)\n",
    "print(before == after)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's okay to reuse one of the input arguments as the output argument, which is helpful if you're keeping a running aggregate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1,2,3,7])\n",
    "b = np.array([3,1,5,2])\n",
    "\n",
    "before = id(a)\n",
    "\n",
    "np.maximum(a, b, a)\n",
    "\n",
    "after = id(a)\n",
    "\n",
    "print(a)\n",
    "print(before == after)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might be wondering how we can use elementwise multiplication to store a result in place, since there's no place for a third argument in the `*` operator. We'll need to use the `multiply` function in this case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1,1,2,3,5])\n",
    "b = np.array([2,3,5,7,11])\n",
    "\n",
    "np.multiply(a, b, a)\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selecting and filtering arrays\n",
    "\n",
    "In conventional Python, you might filter an array by adding a condition to a list comprehension.  So if you had some Fibonacci numbers, you could find the even ones like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fibs = [1,1,2,3,5,8,13,21,34,55,89]\n",
    "[f for f in fibs if f % 2 == 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In NumPy, we have a different way to filter arrays by applying a predicate to every element of an array and generating a *boolean mask*, or an array of truth values indicating whether or not the predicate is satisfied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1,1,2,3,5,8,13,21,34,55,89])\n",
    "a % 2 == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy provides functions that use these masks to filter and transform entire arrays at once.  For example, the `extract` function takes a mask (which coincidentally looks like a conditional expression) and returns an array consisting of only the elements that satisfy the expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.extract(a % 2 == 0, a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also count the number of elements satisfying a conditional expression by summing the `True` values in the boolean mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(a % 2 == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that Python's built-in `sum` function works in this example, but don't use it for this application in general — it has some subtle differences to `numpy.sum`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(a % 2 == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several other functions that work similarly, like `where`, which takes a boolean  mask and two arrays of the same length.  `where` returns an array consisting of the elements from either the second or third arguments:  for elements in the condition that are true, the result will have the corresponding element from the second argument; for other elements, the result will have the corresponding element from the third argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.where(a % 2 == 0, a, np.zeros(len(a), int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's a little clunky to have that `np.zeros(len(a), int)`, isn't it?  I bet we can do better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.where(a % 2 == 0, a, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait a minute — `where` is supposed to take three arrays, but here we're giving it two arrays and an `int`.  What's going on?  Can we always substitute numbers for arrays?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array([2,4,6,8]) * 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like we can!  What about substituting a shorter array for an array of the appropriate length — does that work?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array([2,4,6,8]) * np.array([2, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This doesn't work — we get an error:\n",
    "\n",
    "`ValueError: operands could not be broadcast together with shapes (4,) (2,)`\n",
    "\n",
    "To understand what this means, we'll have to talk a bit more about the words \"broadcast\" and \"shape.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Broadcasting\n",
    "\n",
    "The *shape* of a NumPy array is a representation of its dimensions.  We've only talked about one-dimensional arrays so far, so the shape of one of the arrays we've seen is simply a 1-tuple containing the number of elements it has."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array([1,2,3,4]).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array([1,2,3,4,5,6]).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Broadcasting* is how NumPy deals with function calls that expect arguments of the same shape but receive arguments of varying shapes.  We'll look more at broadcasting in the next notebook; for now, know that if you pass a one-dimensional array and a scalar to a NumPy function, NumPy will treat the scalar as if it were an array with the same shape as the other argument.\n",
    "\n",
    "So:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1,2,3,4])\n",
    "b = np.array([3,3,3,3])\n",
    "\n",
    "# == is elementwise in NumPy, so np.all(x == y) will return true \n",
    "# if each corresponding element in x and y are identical\n",
    "np.all(np.multiply(a, 3) == np.multiply(a, b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case study:  a bit vector class\n",
    "\n",
    "A _bit vector_ is an efficient way to store sets of values from a finite domain.  (By \"finite domain,\" we mean that the set has a fixed number of possible elements and each can be consistently numbered, like letters of the alphabet.)  Think of it as an array of truth values represented as an array of integers.  So if our finite domain is `{red, green, blue}`, we might have a bit vector `001` to represent the set `{blue}` or `101` to represent the set `{red, blue}`.\n",
    "\n",
    "We can use NumPy's bitwise logical operations to efficiently implement intersection and union for bit vectors.  The class below shows how to put this functionality together by using these operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BitVector(object):\n",
    "    def __init__(self, size):\n",
    "        self._size = size\n",
    "        ct = size % 64 == 0 and (size / 64) or (size / 64 + 1)\n",
    "        self._entries = np.zeros(int(ct), np.uint64)\n",
    "    \n",
    "    def __len__(self):\n",
    "        return self._size\n",
    "    \n",
    "    def __getitem__(self, key):\n",
    "        k = int(key)\n",
    "        return (self._entries[int(k / 64)] & np.uint64(1 << (k % 64))) > 0\n",
    "    \n",
    "    def __setitem__(self, key, value):\n",
    "        k = int(key)\n",
    "        if value:\n",
    "            update = np.uint64(1 << key % 64)\n",
    "            self._entries[int(k / 64)] = self._entries[int(k / 64)] | update\n",
    "        else:\n",
    "            update = np.uint64(1 << key % 64)\n",
    "            self._entries[int(k / 64)] = self._entries[int(k / 64)] ^ update\n",
    "        \n",
    "    def intersect(self, other, out=None):\n",
    "        if out is None:\n",
    "            result = BitVector(self._size)\n",
    "        else: \n",
    "            result = out\n",
    "        \n",
    "        np.bitwise_and(self._entries, other._entries, result._entries)\n",
    "        return result\n",
    "    \n",
    "    def union(self, other):\n",
    "        if out is None:\n",
    "            result = BitVector(self._size)\n",
    "        else: \n",
    "            result = out\n",
    "        \n",
    "        np.bitwise_or(self._entries, other._entries, result._entries)\n",
    "        return result"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}