numpy/NumpyIntro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numpy\n",
    "\n",
    "[Numpy](https://www.numpy.org/) is the foundation of most of what you will do with scientific python. Modules like [scipy](https://www.scipy.org/) and [pandas](https://pandas.pydata.org/) are built on top of numpy, and it is the linga franca for most numerical work in python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you are first introduced to python, one of the big selling points is that it isn't statically typed. You can do things like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 1\n",
    "a = 'apple'\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And you don't get complaints from python about `a` being an integer. This is a _big_ advantage for python, and it works with collections as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myList = [1, \n",
    "          3., \n",
    "          'Ian', \n",
    "          [range(3)], \n",
    "          {'not': 'a good idea'}\n",
    "]\n",
    "myList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lists in python are about as general as they could be. This is very flexible and it lets you do some fancy things in Python, but it has a price. The Python interpreter can't make any assumptions about what will come next in a list; everything has to be treated as a generic object. Python does a good job of hiding the complexity of doing this, but as the lists get longer and more complex the overhead gets larger, and eventually perfomance becomes unacceptable.\n",
    "\n",
    "One solution to this is to use a statically typed language like C[<sup>1</sup>](#fn1 \"footnote and tooltip 1\"). There the burden of figuring out object types is left to the programmer, and the compiler can be much more efficient operating on them. A good example would be an array of `double`s. In memory, these will be allocated contiguously so when you need to jump to the 1402th double, you can do it with very simple arithmetic. Python has a much harder time because the memory allocated for your list could be a horrible mixture of all the different things you've stuffed in there. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The `ndarray`\n",
    "\n",
    "Numpy attempts let you keep the advantages of Python without sacrificing the speed of static typing by adding the concept of homogenous collections to python: `ndarray`s. The `ndarray` is the foundational concept in numpy, it is an array object which represents a multidimensional, homogeneous array of fixed-size items and most commonly these items will be numbers.\n",
    "\n",
    "```\n",
    "%%timeit\n",
    "for i in range(1000000):\n",
    "    i*i\n",
    "```\n",
    "```\n",
    "%timeit np.arange(1000000)**2\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "for i in range(1000000):\n",
    "    i*i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit np.arange(1000000)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`ndarray`s look like python lists, but they are fundamentally different, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = [1, 2, 3, 4]\n",
    "b = [5, 6, 7, 8]\n",
    "a+b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "na = np.array([1, 2, 3, 4])\n",
    "nb = np.array([5, 6, 7, 8])\n",
    "na + nb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`numpy` was able to assume that the things in the `ndarray` were the compatible types and vectorize the addition, if we want the same thing with python lists, we have to jump through some hoops"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[ i + j for i, j in zip(a, b) ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at what makes up an `ndarray`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "na = np.array([1,2,3,4,5])\n",
    "na"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`a` has the following attributes,\n",
    "\n",
    "  * **ndim**: The number of dimensions\n",
    "  * **shape**: A tuple giving the sizes of the dimensions\n",
    "  * **size**: The total number of elements (product of the shape items)\n",
    "  * **dtype**: The data type, int64, float64, complex128 etc.\n",
    "  * **itemsize**: The size of individual elements in memory\n",
    "  * **nbytes**: The total memory occupied by the ndarray\n",
    "  * **strides**: Tuple of bytes to step in each dimension when traversing an array\n",
    "  * **data**: The buffer storing the contiguous identically typed items\n",
    "  \n",
    "By knowing these attributes, numpy can use some of the same tricks that statically typed languges use because the Python interpreter can now infer the memory layout.\n",
    "\n",
    "**Exercise**: Check the dtype of `a`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(na)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "na.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "numpy has some convenience methods for creating new ndarrays. As you saw above, one way to create an ndarray is to pass a list with the values to `np.array`. Here are some others...\n",
    "\n",
    "Using `np.array` directly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 = np.array([0, 1, 1, 2, 3, 5, 8, 13])\n",
    "a1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`np.arange` will generate numbers between limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a2 = np.arange(0, 100)\n",
    "a2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arange?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a22 = np.arange(0, 10, dtype=np.float)\n",
    "a22"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`np.linspace` is extremely useful, you tell it where to start and stop and how many samples you need and linspace does the rest. Here we will create numbers 100 numbers between 0 (inclusive) and 1 (inclusive!), linearly spaced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "li1 = np.linspace(0, 1, 100)\n",
    "li1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `np.zeros` function will generate an `ndarray` full of zeros. The first argument is the shape which you can give as an integer (for 1d arrays) or a sequence (for n-dimensional arrays). You can also pass the `dtype=` argument to tell it what type of number you are looking for."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# z1 100 integer zeros\n",
    "z1 = np.zeros(100, dtype=np.int)\n",
    "z1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# z2 a 5,5 array of float64 zeros\n",
    "z2 = np.zeros((5, 5), dtype=np.float64)\n",
    "z2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `np.ones` function does something similar, but with unit value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# o1 100 integers ones\n",
    "o1 = np.ones(100, dtype='float64')\n",
    "o1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# o2 a 5x5 ndarray of complex128 one values. First argument is shape sequence\n",
    "o2 = np.ones((5, 5), dtype='complex128')\n",
    "o2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `np.eye` function genrates a 2D array with ones down the diagonal, zeros elsewhere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# e1 a 5x5 array with ones down the diagonal\n",
    "e1 = np.eye(5, dtype=np.float64)\n",
    "e1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll come back to the `numpy.random` module later on, but it has some useful convenience methods. `np.random.randint` returns random integers. Take a look at the help on the method then create a name `r1` with an array of `4x5` random numbers between `0` and `10`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# r1 a 4x5 array of random integers between 0 & 10, \n",
    "np.random.randint(0, 10, size=(4, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing and Slicing\n",
    "\n",
    "Now that we have some `ndarray`s to play with, lets look at using them. Of course, `ndarray`s are zero indexed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First element of a2\n",
    "a2[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The usual negative number notation works for couting from the end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2nd to last element of a2\n",
    "a2[-2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can update `ndarray` in place by index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a2[-1] = 0\n",
    "a2[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Slicing returns subarrays of the original `ndarray`. Crucially it does this inexpensively by returning a \"view\" on the data rather than copying it. This relies on numpy is using it's knowledge of the memory layout to return only the things you ask for and is a general tactic used by numpy, if it *can* return a view rather than a copy it *will*! If you really need a copy, try the `.copy()` method on `ndarrays`.\n",
    "\n",
    "Slicing uses the same notation as core python: `[start:stop:step]` where `start` is inclusive and `stop` is exclusive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Values between 3 and 19 of a2\n",
    "a2[3:19]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Every third value between 3 and 19\n",
    "a2[3:19:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using negative is allowed for all three parts of the slice, but for the step you have to think a bit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Values between -10 and -2\n",
    "a2[-10:-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Values between 10 and 0, reversed\n",
    "a2[10:2:-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the first argument of the slice is still inclusive and the second is not. If we omit a value when specifying the slice `start` defaults to 0, `end` defaults to the last element and `step` defaults to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a2[::3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For multi-dimensional arrays the indexing notation is similar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.arange(100)\n",
    "b.shape = (10, 10)\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How about row 0, column 3 (remember python is 0 indexed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b[0,3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b[:, -1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or the fifth column of the first two rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b[0:2, 4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b[0:2, 4:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b[::2, 9]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reshaping Arrays\n",
    "\n",
    "Sometimes it is convenient to reshape arrays. I did this above by setting the `.shape` attribute but numpy arrays also have a reshape method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.arange(27)\n",
    "c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`reshape` expects a sequence as the first argument (we gave a tuple) so"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = c.reshape((3, 9))\n",
    "d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reshaping isn't enough to provoke numpy to copy the data, all it needs to do is make a new view on the same data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d.base is c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d.base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another common reshaping task is to add dimension(s) to an existing array. numpy has a special `newaxis` object for this task. This is a powerful idea when combined with numpy's broadcasting rules (see later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e = np.arange(10)\n",
    "e[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e[np.newaxis, np.newaxis, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stacking & Splitting ndarrays\n",
    "\n",
    "You can combine general ndarrays with the `np.concatenate` and split them with `np.split`. There are also a number of convenience methods for commonly used shapes.\n",
    "\n",
    " * `hstack`\n",
    " * `vstack`\n",
    " * `hsplit`\n",
    " * `vsplit`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`hstack` is short for horizontal stack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3, 4])\n",
    "b = np.array([4, 3, 2, 1])\n",
    "\n",
    "np.hstack((a, b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For `2D` arrays this means we are joining columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.hstack((a[np.newaxis, :].T, b[np.newaxis, :].T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "vstack does the equivalent thing but vertically stacking (along axis 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.vstack((a, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.vstack((a[:, np.newaxis], b[:, np.newaxis]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In higher dimensions `concatenate` or `stack` should do what you need, but you need to manually tell it which axis to use for the stacking.\n",
    "\n",
    "`np.split` goes in the opposite direction. It will try to produce sub-arrays of equal size. Again there are `hsplit` and `vsplit` variants for common use cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.split(np.arange(64).reshape((8, 8)), 2, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.split(np.arange(64).reshape((8, 8)), 2, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Universal Functions (ufuncs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The real reason for using numpy is so you can do numerical operations, _quickly_. Python uses a concenpt called `ufuncs` or universal function. A ufunc is a function which operates on `ndarrays` element-by-element. In other language a ufunc is a vectorized wrapper around a function which can transform one `ndarray` element-by-element to produce another. This element-by-element behaviour is fundamentally different from the usual python behaviour.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "The key to writing fast numeric python code is: **Avoid for & while loops as far as you can, use numpy ufuncs as far as possible**\n",
    "\n",
    "\n",
    "Lets start with basic arithmetic operations. Numpy can use it's internal broadcasting to do these quickly and efficiently"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The usual operations are available\n",
    "\n",
    "  * +: addition\n",
    "  * -: subtraction\n",
    "  * *: multiplication\n",
    "  * /: division\n",
    "  * //: integer division\n",
    "  * **: power operator\n",
    "  * %: modulo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operations are element by element, and you can build up more complicated expressions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "la = np.linspace(0, 1, 100)\n",
    "\n",
    "(la ** 2 + la) / (la + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**: Try to calculate the terms of this sum as an `ndarray`\n",
    "$$\n",
    "\\sqrt{12}\\sum_{k=0}^{10}\\frac{(-3)^{-k}}{2k+1}\n",
    "$$\n",
    "\n",
    "*Hints*: Think term by term. `np.arange(10)` will give you an `ndarray` of k values, next raise `-3. * np.ones(k.shape)` to those powers. If you can calculate the terms you can use the `.sum()` method to sum them all up. How close to $\\pi$ did you get?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.arange(11)\n",
    "k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = (np.sqrt(12) * (-3.)**-k / (2*k + 1)).sum() - np.pi\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(12) * ((-3 * np.ones(len(k)))**-k / (2 * k + 1)).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The operators we were using `+,-,/,...` actually correspond to functions (`ufuncs`)\n",
    "\n",
    "|operator|function|description|\n",
    "|--------|--------|-----------|\n",
    "| + | np.add | Addition |\n",
    "| - | np.subtract | Subtraction |\n",
    "| - | np.negative | Unary negation |\n",
    "| * | np.multiply | Multiplication |\n",
    "| / | np.division | Ordinary floating point division |\n",
    "| // | np.floor_divide | floor (integer) division |\n",
    "| % | np.mod | Modulo/Remainder division |\n",
    "\n",
    "You can use either syntax, but in the function notation there are lots more functions to play with e.g.\n",
    "\n",
    "| function | description |\n",
    "|----------|-------------|\n",
    "| np.sin   | sin function |\n",
    "| np.cos   | cos function |\n",
    "| np.tan   | tan function |\n",
    "| np.abs   | absolute value |\n",
    "| np.exp   | exponential |\n",
    "| np.log   | natural log |\n",
    "| np.log2  | log base 2 |\n",
    "| np.log10 | log base 10 |\n",
    "|  ...     |    ...      |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = np.linspace(0, 2*np.pi, 25)\n",
    "p2 = np.sin(p1)\n",
    "p2\n",
    "\n",
    "# p2 sin of p1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1 = np.linspace(0, 1, 10) + np.linspace(1, 2, 10)*1j\n",
    "q1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.abs(q1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aggregate Functions\n",
    "\n",
    "Aggregate functions take an `ndarray` and reduce it along one (or more) axes. An example would be taking an array of numbers and calculating the mean value..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r1 = np.linspace(0, 10, 100)\n",
    "r1.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But there are lots of aggregate functions\n",
    "\n",
    "  * `min`: Minumum value\n",
    "  * `max`: Maximum value\n",
    "  * `sum`: Sum values\n",
    "  * `prod`: Product of values\n",
    "  * `mean`: Arthmetic mean\n",
    "  * `std`: Standard deviation\n",
    "  * `var`: Variance\n",
    "  * `argmin`: indices of the minimum value\n",
    "  * `argmax`: indices of the maximum value\n",
    "  * `all`: is a condition true in all elements\n",
    "  * `any`: is a condition true in any elements\n",
    "  * `allclose`: All the values are within a small tolerance **really useful!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default is to reduce along all axes, if you want to reduce along a specific axis you can pass that as an argument (the axes you specify are the ones which get squashed) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = np.arange(50)\n",
    "s2 = s1.reshape(5,10)\n",
    "s2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2.mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2.mean(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For binary operations (e.g. addition) you can also do reduction, so starting from\n",
    "\n",
    "[1, 3, 5, 7, 9]\n",
    "\n",
    "`np.add.reduce` will add 1 to 3, add that to 5 and so on, ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = np.arange(1, 10, 2)\n",
    "t1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.add.reduce(t1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "N.B. In this case we used a function from the module and passed in our ndarray."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing with Booleans and \"fancy\" indexing\n",
    "\n",
    "numpy has some more advanced tricks up it's sleeve when indexing. If we do a boolean test on an `ndarray` the result will be the same shape as the `ndarray` itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u1 = np.linspace(0, 1, 100)\n",
    "\n",
    "u1 < 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "numpy will let us use this to index arrays of the same shape (see also `np.where`)..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u1[u1 < 0.5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also do other things with the boolean values like `.count_nonzero` them, checking if `.all` values or `.any` values are true etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.count_nonzero(u1 < 0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.any(u1 > 0.75)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.any(u1[u1 < 0.5] > 0.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**N.B. Logical/Boolean indexing will generally force a copy of the ndarray rather than a view because, in general, the sequence of True/False values can be irregular**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fancy Indexing\n",
    "\n",
    "Fancy Indexing is the idea of using another array of indices, it is useful when the combinations you want to select become a bit more complicated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v1=np.arange(27)[::-1]\n",
    "v1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v1[[1, 4, 6]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v2 = v1.reshape(3,9)\n",
    "v2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In higher dimensions think of zipping together the arguements, e.g. 1st row, 4th column, ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v2[[0, 1, 1], [3, 7, 8]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*N.B. Fancy indexing usually creates copies of the `ndarray` because you usually can't reconstruct the selection with simple algebra*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Broadcasting\n",
    "\n",
    "`numpy` is at it's most efficient when it is operating element by element, but not all arrays are the same size. To work around this, `numpy` implements a set of rules under the name of `broadcasting` to make `ndarray`s conform whenever possible. This is great news, in that you don't have to worry about doing that yourself, but it is important to understand the rules so that you know how `numpy` will behave when combining differnt shaped `nbdarrays`. Conceptually, think of\n",
    "```python\n",
    "np.arange(10) * 5\n",
    "```\n",
    "`numpy` wants to operate element-by-element, but `5` isn't an `ndarray`, it's just a number. If we could prompte 5 to be a `1d` narray and put 5's in in every place `numpy` would be happy. This is the basic idea of broadcasting, in summary\n",
    "\n",
    "1. Given two arrays of different dimensions, prepend 1, to the shape of the smaller array\n",
    "1. Dimensions of size 1 are repeated (without copying)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(15)\n",
    "a = a.reshape(3, 5)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.arange(5)\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thinking through the rules manually, `b` has the smaller dimensions (1 vs. 2) so a dimension of length one will be prepended to `b`. `b` will then be repeated 3 times to conform with the shape of a."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "btmp = np.repeat(b[np.newaxis, :], 3, axis=0)\n",
    "btmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a * btmp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now again, but using the rules implicitly via broadcasting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a * b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Random\n",
    "\n",
    "Numpy has a few important submodules but `np.random` is probably the most important. As you might expect, it lets you work with random numbers but it goes far beyond simple a few random integers. You can make random samples from `ndarray`s, shuffle and permute sequences, draw samples from 30+ distributions and you can do it all quickly! Of course these are going to be pseudo-random numbers.\n",
    "\n",
    "`numpy.random` uses the concept of a `Generators` to implement sampling. The idea is you create a generator object then call methods on that generator to sample from the various distributions. The original `Generator` will normally get it's entropy from a (hopefully) reliable source then you can keep asking it for the `__next__` random elements distributed however you need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(47)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we've seeded the `Generator` with a specific value so the results are reporoducible but normally you would just call `np.random.default_rng()` to get a random value from the OS. Now we can sample from various distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng.normal(5.0, 1.0, (64, 64))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The binomial distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng.binomial(10, .5, size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5 random integers between 10 and 20, uniformly distributed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng.integers(10, 20, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`np.random` also has functions for shuffling arrays in-place (`np.shuffle`) and selecting elements at random from `ndarrays` (`np.choice`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(10)\n",
    "np.random.shuffle(a)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(10)\n",
    "np.random.choice(a, 3, replace=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also weight the selections in `np.choice` with probabilities. e.g. Here are the letter frequencies of ordinary english text."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "letter_freq = {' ': 0.19266420666588144,\n",
    " 'e': 0.09680968984305797,\n",
    " 't': 0.07140241019840815,\n",
    " 'a': 0.06361581392196947,\n",
    " 'o': 0.06183938572048604,\n",
    " 'i': 0.05349452953695084,\n",
    " 'n': 0.0521037201730283,\n",
    " 'h': 0.051232447485652234,\n",
    " 's': 0.049280151278754014,\n",
    " 'r': 0.04524648142979075,\n",
    " 'd': 0.03375374929612461,\n",
    " 'l': 0.03124157971419029,\n",
    " 'u': 0.02392934301198968,\n",
    " 'm': 0.021518821910249234,\n",
    " 'w': 0.020208685943305965,\n",
    " 'c': 0.02004895261728702,\n",
    " 'f': 0.016262143363080305,\n",
    " 'y': 0.016250849087503207,\n",
    " 'g': 0.013559584564274916,\n",
    " 'p': 0.012933559003715817,\n",
    " ',': 0.012259129404969158,\n",
    " '.': 0.01200420147051468,\n",
    " 'b': 0.011160357738111564,\n",
    " 'v': 0.0076220225466009876,\n",
    " 'k': 0.006392559976636984,\n",
    " 'x': 0.001353699601312072,\n",
    " 'j': 0.0008260955850676769,\n",
    " 'q': 0.0006663622590487316,\n",
    " 'z': 0.00031946665203789066\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use those relative frequency values as probabilities (they sum to ~1) and generate a random sample of letters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "letters = np.array(list(letter_freq.keys()))\n",
    "\n",
    "probabilities = np.array(list(letter_freq.values()))\n",
    "\n",
    "chosen = np.random.choice(letters, 1000, p=probabilities)\n",
    "''.join(chosen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other stuff?\n",
    "  * `_ix()`\n",
    "  * a.copy()\n",
    "  * np.linalg\n",
    "  * np.matrixlib\n",
    "  * np.fft\n",
    "  * Applications? Simple monte-carlo? Images?\n",
    "  * np.set_printoptions?\n",
    "  * memory layouts and strides?\n",
    "  * .flatten (copy) and ravel (view if possible)\n",
    "  * mgrid & ogrid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Footnotes\n",
    "[1]<span id=\"fn1\"> If you dig deep enough in some of the numpy/scipy code you will find that the actual instructions you are executing were compiled from fortran and C. In general you can pass things quite easily between existing libraries and python, but fortran and C use different storage orders for multi-dimensional arrays so you have to be a little bit careful_.</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra stuff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A manual gaussian blur"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = (16, 12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = plt.imread('blackbuck.png')\n",
    "\n",
    "plt.imshow(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We won't worry about the edges, we'll just throw away cells we can't compute, our stencil needs one value in each direction so (ignoring the 3rd dimension for now), think about this sumarray `[1:-1, 1:-1`. To calculate the average we need the values at those points along with the values offset 1 square up, one square down, one square left and one square right from those points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newimg = (\n",
    "    img[ :-2 , 1:-1 , :] +        # above\n",
    "    img[2:  ,  1:-1:, :] +        # below\n",
    "    img[1:-1,   :-2 , :] +        # left\n",
    "    img[1:-1,  2:   , :] +        # right\n",
    "    img[1:-1,  1:-1 , :] )/ 5     # center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(newimg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(10):\n",
    "    newimg = (\n",
    "        newimg[ :-2 , 1:-1 , :] +        # above\n",
    "        newimg[2:  ,  1:-1:, :] +        # below\n",
    "        newimg[1:-1,   :-2 , :] +        # left\n",
    "        newimg[1:-1,  2:   , :] +        # right\n",
    "        newimg[1:-1,  1:-1 , :] )/ 5     # center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(newimg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}