04_numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numpy\n",
    "\n",
    "\n",
    "This section introduces you to [`numpy`](http://www.numpy.org/), Python's\n",
    "numerical computing library.\n",
    "\n",
    "\n",
    "Numpy is not actually part of the standard Python library. But it is a\n",
    "fundamental part of the Python ecosystem - it forms the basis for many\n",
    "important Python libraries, and it (along with its partners\n",
    "[`scipy`](https://www.scipy.org/), [`matplotlib`](https://matplotlib.org/) and\n",
    "[`pandas`](https://pandas.pydata.org/)) is what makes Python an attractive\n",
    "alternative to Matlab as a scientific computing platform.\n",
    "\n",
    "\n",
    "The `fslpython` environment currently includes [Numpy\n",
    "1.11.1](https://docs.scipy.org/doc/numpy-1.11.0/index.html), which is a little\n",
    "out of date, but we will update it for the next release of FSL.\n",
    "\n",
    "\n",
    "## Contents\n",
    "\n",
    "\n",
    "* [The Python list versus the Numpy array](#the-python-list-versus-the-numpy-array)\n",
    "* [Numpy basics](#numpy-basics)\n",
    " * [Creating arrays](#creating-arrays)\n",
    " * [Loading text files](#loading-text-files)\n",
    " * [Array properties](#array-properties)\n",
    " * [Descriptive statistics](#descriptive-statistics)\n",
    " * [Reshaping and rearranging arrays](#reshaping-and-rearranging-arrays)\n",
    "* [Operating on arrays](#operating-on-arrays)\n",
    " * [Scalar operations](#scalar-operations)\n",
    " * [Multi-variate operations](#multi-variate-operations)\n",
    " * [Matrix multplication](#matrix-multiplication)\n",
    " * [Broadcasting](#broadcasting)\n",
    " * [Linear algebra](#linear-algebra)\n",
    "* [Array indexing](#array-indexing)\n",
    " * [Indexing multi-dimensional arrays](#indexing-multi-dimensional-arrays)\n",
    " * [Boolean indexing](#boolean-indexing)\n",
    " * [Coordinate array indexing](#coordinate-array-indexing)\n",
    "* [Exercises](#exercises)\n",
    " * [Load an array from a file and do stuff with it](#load-an-array-from-a-file-and-do-stuff-with-it)\n",
    " * [Concatenate affine transforms](#concatenate-affine-transforms)\n",
    "\n",
    "* [Appendix A: Generating random numbers](#appendix-generating-random-numbers)\n",
    "* [Appendix B: Importing Numpy](#appendix-importing-numpy)\n",
    "* [Appendix C: Vectors in Numpy](#appendix-vectors-in-numpy)\n",
    "\n",
    "* [Useful references](#useful-references)\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"the-python-list-versus-the-numpy-array\"></a>\n",
    "## The Python list versus the Numpy array\n",
    "\n",
    "\n",
    "Numpy adds a new data type to the Python language - the `array` (more\n",
    "specifically, the `ndarray`). A Numpy `array` is a N-dimensional array of\n",
    "homogeneously-typed numerical data.\n",
    "\n",
    "\n",
    "You have already been introduced to the Python `list`, which you can easily\n",
    "use to store a handful of numbers (or anything else):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = [10, 8, 12, 14, 7, 6, 11]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You could also emulate a 2D or ND matrix by using lists of lists, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyz_coords = [[-11.4,   1.0,  22.6],\n",
    "              [ 22.7, -32.8,  19.1],\n",
    "              [ 62.8, -18.2, -34.5]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simple tasks, you could stick with processing your data using python\n",
    "lists, and the built-in\n",
    "[`math`](https://docs.python.org/3.5/library/math.html) library. And this\n",
    "might be tempting, because it does look quite a lot like what you might type\n",
    "into Matlab.\n",
    "\n",
    "\n",
    "But __BEWARE!__ A Python list is a terrible data structure for scientific\n",
    "computing!\n",
    "\n",
    "\n",
    "This is a major source of confusion for people who are learning Python, and\n",
    "are trying to write efficient code. It is _crucial_ to be able to distinguish\n",
    "between a Python list and a Numpy array.\n",
    "\n",
    "\n",
    "___Python list == Matlab cell array:___ A list in Python is akin to a cell\n",
    "array in Matlab - they can store anything, but are extremely inefficient, and\n",
    "unwieldy when you have more than a couple of dimensions.\n",
    "\n",
    "\n",
    "___Numy array == Matlab matrix:___ These are in contrast to the Numpy array\n",
    "and Matlab matrix, which are both thin wrappers around a contiguous chunk of\n",
    "memory, and which provide blazing-fast performance (because behind the scenes\n",
    "in both Numpy and Matlab, it's C, C++ and FORTRAN all the way down).\n",
    "\n",
    "\n",
    "So you should strongly consider turning those lists into Numpy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "data = np.array([10, 8, 12, 14, 7, 6, 11])\n",
    "\n",
    "xyz_coords = np.array([[-11.4,   1.0,  22.6],\n",
    "                       [ 22.7, -32.8,  19.1],\n",
    "                       [ 62.8, -18.2, -34.5]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you look carefully at the code above, you will notice that we are still\n",
    "actually using Python lists. We have declared our data sets in exactly the\n",
    "same way as we did earlier, by denoting them with square brackets `[` and `]`.\n",
    "\n",
    "\n",
    "The key difference here is that these lists immediately get converted into\n",
    "Numpy arrays, by passing them to the `np.array` function.  To clarify this\n",
    "point, we could rewrite this code in the following equivalent manner:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Define our data sets as python lists\n",
    "data       = [10, 8, 12, 14, 7, 6, 11]\n",
    "xyz_coords = [[-11.4,   1.0,  22.6],\n",
    "              [ 22.7, -32.8,  19.1],\n",
    "              [ 62.8, -18.2, -34.5]]\n",
    "\n",
    "# Convert them to numpy arrays\n",
    "data       = np.array(data)\n",
    "xyz_coords = np.array(xyz_coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, in practice, we would never create a Numpy array in this way - we\n",
    "will be loading our data from text or binary files directly into a Numpy\n",
    "array, thus completely bypassing the use of Python lists and the costly\n",
    "list-to-array conversion.  I'm emphasising this to help you understand the\n",
    "difference between Python lists and Numpy arrays. Apologies if you've already\n",
    "got it, forgiveness please.\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"numpy-basics\"></a>\n",
    "## Numpy basics\n",
    "\n",
    "\n",
    "Let's get started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"creating-arrays\"></a>\n",
    "### Creating arrays\n",
    "\n",
    "\n",
    "Numpy has quite a few functions which behave similarly to their equivalents in\n",
    "Matlab:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('np.zeros gives us zeros:                       ', np.zeros(5))\n",
    "print('np.ones gives us ones:                         ', np.ones(5))\n",
    "print('np.arange gives us a range:                    ', np.arange(5))\n",
    "print('np.linspace gives us N linearly spaced numbers:', np.linspace(0, 1, 5))\n",
    "print('np.random.random gives us random numbers [0-1]:', np.random.random(5))\n",
    "print('np.random.randint gives us random integers:    ', np.random.randint(1, 10, 5))\n",
    "print('np.eye gives us an identity matrix:')\n",
    "print(np.eye(4))\n",
    "print('np.diag gives us a diagonal matrix:')\n",
    "print(np.diag([1, 2, 3, 4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> See the [appendix](#appendix-generating-random-numbers) for more information\n",
    "> on generating random numbers in Numpy.\n",
    "\n",
    "\n",
    "The `zeros` and `ones` functions can also be used to generate N-dimensional\n",
    "arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.zeros((3, 4))\n",
    "o = np.ones((2, 10))\n",
    "print(z)\n",
    "print(o)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that, in a 2D Numpy array, the first axis corresponds to rows, and the\n",
    "> second to columns - just like in Matlab.\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"loading-text-files\"></a>\n",
    "### Loading text files\n",
    "\n",
    "\n",
    "The `numpy.loadtxt` function is capable of loading numerical data from\n",
    "plain-text files. By default it expects space-separated data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt('04_numpy/space_separated.txt')\n",
    "print('data in 04_numpy/space_separated.txt:')\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But you can also specify the delimiter to expect<sup>1</sup>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt('04_numpy/comma_separated.txt', delimiter=',')\n",
    "print('data in 04_numpy/comma_separated.txt:')\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> <sup>1</sup> And many other things such as file headers, footers, comments,\n",
    "> and newline characters - see the\n",
    "> [docs](https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html)\n",
    "> for more information.\n",
    "\n",
    "\n",
    "  Of course you can also save data out to a text file just as easily, with\n",
    "[`numpy.savetxt`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.random.randint(1, 10, (10, 10))\n",
    "np.savetxt('mydata.txt', data, delimiter=',', fmt='%i')\n",
    "\n",
    "with open('mydata.txt', 'rt') as f:\n",
    "    for line in f:\n",
    "        print(line.strip())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"array-properties\"></a>\n",
    "### Array properties\n",
    "\n",
    "\n",
    "Numpy is a bit different than Matlab in the way that you interact with\n",
    "arrays. In Matlab, you would typically pass an array to a built-in function,\n",
    "e.g. `size(M)`, `ndims(M)`, etc. In contrast, a Numpy array is a Python\n",
    "object which has _attributes_ that contain basic information about the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.zeros((2, 3, 4))\n",
    "print(z)\n",
    "print('Shape:                     ', z.shape)\n",
    "print('Number of dimensions:      ', z.ndim)\n",
    "print('Number of elements:        ', z.size)\n",
    "print('Data type:                 ', z.dtype)\n",
    "print('Number of bytes:           ', z.nbytes)\n",
    "print('Length of first dimension: ', len(z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> As depicted above, passing a Numpy array to the built-in `len` function will\n",
    "> only give you the length of the first dimension, so you will typically want\n",
    "> to avoid using it - use the `size` attribute instead.\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"descriptive-statistics\"></a>\n",
    "### Descriptive statistics\n",
    "\n",
    "\n",
    "Similarly, a Numpy array has a set of methods<sup>2</sup> which allow you to\n",
    "calculate basic descriptive statistics on an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.random(10)\n",
    "print('a: ', a)\n",
    "print('min:          ', a.min())\n",
    "print('max:          ', a.max())\n",
    "print('index of min: ', a.argmin())  # remember that in Python, list indices\n",
    "print('index of max: ', a.argmax())  # start from zero - Numpy is the same!\n",
    "print('mean:         ', a.mean())\n",
    "print('variance:     ', a.var())\n",
    "print('stddev:       ', a.std())\n",
    "print('sum:          ', a.sum())\n",
    "print('prod:         ', a.prod())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> <sup>2</sup> Python, being an object-oriented language, distinguishes\n",
    "> between _functions_ and _methods_. Hopefully we all know what a function is\n",
    "> - a _method_ is simply the term used to refer to a function that is\n",
    "> associated with a specific object. Similarly, the term _attribute_ is used\n",
    "> to refer to some piece of information that is attached to an object, such as\n",
    "> `z.shape`, or `z.dtype`.\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"reshaping-and-rearranging-arrays\"></a>\n",
    "### Reshaping and rearranging arrays\n",
    "\n",
    "\n",
    "A numpy array can be reshaped very easily, using the `reshape` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(16).reshape(4, 4)\n",
    "b = a.reshape((2, 8))\n",
    "print('a:')\n",
    "print(a)\n",
    "print('b:')\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this does not modify the underlying data in any way - the `reshape`\n",
    "method returns a _view_ of the same array, just indexed differently:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[3, 3] = 12345\n",
    "b[0, 7] = 54321\n",
    "print('a:')\n",
    "print(a)\n",
    "print('b:')\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you need to create a reshaped copy of an array, use the `np.array`\n",
    "function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(16).reshape((4, 4))\n",
    "b = np.array(a.reshape(2, 8))\n",
    "a[3, 3] = 12345\n",
    "b[0, 7] = 54321\n",
    "print('a:')\n",
    "print(a)\n",
    "print('b:')\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `T` attribute is a shortcut to obtain the transpose of an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(16).reshape((4, 4))\n",
    "print(a)\n",
    "print(a.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `transpose` method allows you to obtain more complicated rearrangements\n",
    "of an N-dimensional array's axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(24).reshape((2, 3, 4))\n",
    "b = a.transpose((2, 0, 1))\n",
    "print('a: ', a.shape)\n",
    "print(a)\n",
    "print('b:', b.shape)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note again that the `T` attribute and `transpose` method return _views_ of\n",
    "> your array.\n",
    "\n",
    "\n",
    "Numpy has some useful functions which allow you to concatenate or stack\n",
    "multiple arrays into one. The `concatenate` function does what it says on the\n",
    "tin:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.zeros(3)\n",
    "b = np.ones(3)\n",
    "\n",
    "print('1D concatenation:', np.concatenate((a, b)))\n",
    "\n",
    "a = np.zeros((3, 3))\n",
    "b = np.ones((3, 3))\n",
    "\n",
    "print('2D column-wise concatenation:')\n",
    "print(np.concatenate((a, b), axis=1))\n",
    "\n",
    "print('2D row-wise concatenation:')\n",
    "\n",
    "# The axis parameter defaults to 0,\n",
    "# so it is not strictly necessary here.\n",
    "print(np.concatenate((a, b), axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `hstack`, `vstack` and `dstack` functions allow you to concatenate vectors\n",
    "or arrays along the first, second, or third dimension respectively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.zeros(3)\n",
    "b = np.ones(3)\n",
    "\n",
    "print('a: ', a)\n",
    "print('b: ', b)\n",
    "\n",
    "hstacked = np.hstack((a, b))\n",
    "vstacked = np.vstack((a, b))\n",
    "dstacked = np.dstack((a, b))\n",
    "\n",
    "print('hstacked: (shape {}):'.format(hstacked.shape))\n",
    "print( hstacked)\n",
    "print('vstacked: (shape {}):'.format(vstacked.shape))\n",
    "print( vstacked)\n",
    "print('dstacked: (shape {}):'.format(dstacked.shape))\n",
    "print( dstacked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"operating-on-arrays\"></a>\n",
    "## Operating on arrays\n",
    "\n",
    "\n",
    "If you are coming from Matlab, you should read this section as, while many\n",
    "Numpy operations behave similarly to Matlab, there are a few important\n",
    "behaviours which differ from what you might expect.\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"scalar-operations\"></a>\n",
    "### Scalar operations\n",
    "\n",
    "\n",
    "All of the mathematical operators you know and love can be applied to Numpy\n",
    "arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(1, 10).reshape((3, 3))\n",
    "print('a:')\n",
    "print(a)\n",
    "print('a + 2:')\n",
    "print( a + 2)\n",
    "print('a * 3:')\n",
    "print( a * 3)\n",
    "print('a % 2:')\n",
    "print( a % 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"multi-variate-operations\"></a>\n",
    "### Multi-variate operations\n",
    "\n",
    "\n",
    "Many operations in Numpy operate on an element-wise basis. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.ones(5)\n",
    "b = np.random.randint(1, 10, 5)\n",
    "\n",
    "print('a:     ', a)\n",
    "print('b:     ', b)\n",
    "print('a + b: ', a + b)\n",
    "print('a * b: ', a * b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This also extends to higher dimensional arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.ones((4, 4))\n",
    "b = np.arange(16).reshape((4, 4))\n",
    "\n",
    "print('a:')\n",
    "print(a)\n",
    "print('b:')\n",
    "print(b)\n",
    "\n",
    "print('a + b')\n",
    "print(a + b)\n",
    "print('a * b')\n",
    "print(a * b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait ... what's that you say? Oh, I couldn't understand because of all the\n",
    "froth coming out of your mouth. I guess you're angry that `a * b` didn't give\n",
    "you the matrix product, like it would have in Matlab.  Well all I can say is\n",
    "that Numpy is not Matlab. Matlab operations are typically consistent with\n",
    "linear algebra notation. This is not the case in Numpy. Get over it. Take a\n",
    "calmative.\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"matrix-multiplication\"></a>\n",
    "### Matrix multiplication\n",
    "\n",
    "\n",
    "When your heart rate has returned to its normal caffeine-induced state, you\n",
    "can use the `@` operator or the `dot` method, to perform matrix\n",
    "multiplication:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(1, 5).reshape((2, 2))\n",
    "b = a.T\n",
    "\n",
    "print('a:')\n",
    "print(a)\n",
    "print('b:')\n",
    "print(b)\n",
    "\n",
    "print('a @ b')\n",
    "print(a @ b)\n",
    "\n",
    "print('a.dot(b)')\n",
    "print(a.dot(b))\n",
    "\n",
    "print('b.dot(a)')\n",
    "print(b.dot(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The `@` matrix multiplication operator is a relatively recent addition to\n",
    "> Python and Numpy, so you might not see it all that often in existing\n",
    "> code. But it's here to stay, so if you don't need to worry about\n",
    "> backwards-compatibility, go ahead and use it!\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"broadcasting\"></a>\n",
    "### Broadcasting\n",
    "\n",
    "\n",
    "One of the coolest features of Numpy is _broadcasting_<sup>3</sup>.\n",
    "Broadcasting allows you to perform element-wise operations on arrays which\n",
    "have a different shape. For each axis in the two arrays, Numpy will implicitly\n",
    "expand the shape of the smaller axis to match the shape of the larger one. You\n",
    "never need to use `repmat` ever again!\n",
    "\n",
    "\n",
    "> <sup>3</sup>Mathworks have shamelessly stolen Numpy's broadcasting behaviour\n",
    "> and included it in Matlab versions from 2016b onwards, referring to it as\n",
    "> _implicit expansion_.\n",
    "\n",
    "\n",
    "Broadcasting allows you to, for example, add the elements of a 1D vector to\n",
    "all of the rows or columns of a 2D array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(9).reshape((3, 3))\n",
    "b = np.arange(1, 4)\n",
    "print('a:')\n",
    "print(a)\n",
    "print('b: ', b)\n",
    "print('a * b (row-wise broadcasting):')\n",
    "print(a * b)\n",
    "print('a * b.T (column-wise broadcasting):')\n",
    "print(a * b.reshape(-1, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Here we used a handy feature of the `reshape` method - if you pass `-1` for\n",
    "> the size of one dimension, it will automatically determine the size to use\n",
    "> for that dimension. Take a look at [the\n",
    "> appendix](#appendix-vectors-in-numpy) for a discussion on vectors in Numpy.\n",
    "\n",
    "\n",
    "Here is a more useful example, where we use broadcasting to de-mean the rows\n",
    "or columns of an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(9).reshape((3, 3))\n",
    "print('a:')\n",
    "print(a)\n",
    "\n",
    "print('a (cols demeaned):')\n",
    "print(a - a.mean(axis=0))\n",
    "\n",
    "print('a (rows demeaned):')\n",
    "print(a - a.mean(axis=1).reshape(-1, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> As demonstrated above, many functions in Numpy accept an `axis` parameter,\n",
    "> allowing you to apply the function along a specific axis. Omitting the\n",
    "> `axis` parameter will apply the function to the whole array.\n",
    "\n",
    "\n",
    "Broadcasting can sometimes be confusing, but the rules which Numpy follows to\n",
    "align arrays of different sizes, and hence determine how the broadcasting\n",
    "should be applied, are pretty straightforward. If something is not working,\n",
    "and you can't figure out why refer to the [official\n",
    "documentation](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html).\n",
    "\n",
    "\n",
    "<a class=\"anchor\" id=\"linear-algebra\"></a>\n",
    "### Linear algebra\n",
    "\n",
    "\n",
    "Numpy is first and foremost a library for general-purpose numerical computing.\n",
    "But it does have a range of linear algebra functionality, hidden away in the\n",
    "[`numpy.linalg`](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html)\n",
    "module. Here are a couple of quick examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.linalg as npla\n",
    "\n",
    "a = np.array([[1, 2, 3,  4],\n",
    "              [0, 5, 6,  7],\n",
    "              [0, 0, 8,  9],\n",
    "              [0, 0, 0, 10]])\n",
    "\n",
    "print('a:')\n",
    "print(a)\n",
    "\n",
    "print('inv(a)')\n",
    "print(npla.inv(a))\n",
    "\n",
    "eigvals, eigvecs = npla.eig(a)\n",
    "\n",
    "print('eigenvalues and vectors of a:')\n",
    "for val, vec in zip(eigvals, eigvecs):\n",
    "    print('{:2.0f} - {}'.format(val, vec))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"array-indexing\"></a>\n",
    "## Array indexing\n",
    "\n",
    "\n",
    "Just like in Matlab, slicing up your arrays is a breeze in Numpy.  If you are\n",
    "after some light reading, you might want to check out the [comprehensive Numpy\n",
    "Indexing\n",
    "reference](https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html).\n",
    "\n",
    "\n",
    "> As with indexing regular Python lists, array indices start from 0, and end\n",
    "> indices (if specified) are exclusive.\n",
    "\n",
    "\n",
    "Let's whet our appetites with some basic 1D array slicing. Numpy supports the\n",
    "standard Python\n",
    "[__slice__](https://www.pythoncentral.io/how-to-slice-listsarrays-and-tuples-in-python/)\n",
    "notation for indexing, where you can specify the start and end indices, and\n",
    "the step size, via the `start:stop:step` syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(10)\n",
    "\n",
    "print('a:                              ', a)\n",
    "print('first element:                  ', a[0])\n",
    "print('first two elements:             ', a[:2])\n",
    "print('last element:                   ', a[a.shape[0] - 1])\n",
    "print('last element again:             ', a[-1])\n",
    "print('last two elements:              ', a[-2:])\n",
    "print('middle four elements:           ', a[3:7])\n",
    "print('Every second element:           ', a[1::2])\n",
    "print('Every second element, reversed: ', a[-1::-2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that slicing an array in this way returns a _view_, not a copy, into that\n",
    "array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(10)\n",
    "print('a:', a)\n",
    "every2nd = a[::2]\n",
    "print('every 2nd:', every2nd)\n",
    "every2nd += 10\n",
    "print('a:', a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"indexing-multi-dimensional-arrays\"></a>\n",
    "### Indexing multi-dimensional arrays\n",
    "\n",
    "\n",
    "Multi-dimensional array indexing works in much the same way as one-dimensional\n",
    "indexing but with, well, more dimensions. Use commas within the square\n",
    "brackets to separate the slices for each dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(25).reshape((5, 5))\n",
    "print('a:')\n",
    "print(a)\n",
    "print(' First row:     ', a[  0, :])\n",
    "print(' Last row:      ', a[ -1, :])\n",
    "print(' second column: ', a[  :, 1])\n",
    "print(' Centre:')\n",
    "print(a[1:4, 1:4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For arrays with more than two dimensions, the ellipsis (`...`) is a handy\n",
    "feature - it allows you to specify a slice comprising all elements along\n",
    "more than one dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(27).reshape((3, 3, 3))\n",
    "print('a:')\n",
    "print(a)\n",
    "print('All elements at x=0:')\n",
    "print(a[0, ...])\n",
    "print('All elements at z=2:')\n",
    "print(a[..., 2])\n",
    "print('All elements at x=0, z=2:')\n",
    "print(a[0, ..., 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"boolean-indexing\"></a>\n",
    "### Boolean indexing\n",
    "\n",
    "\n",
    "A numpy array can be indexed with a boolean array of the same shape. For\n",
    "example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(10)\n",
    "\n",
    "print('a:                          ', a)\n",
    "print('a > 5:                      ', a > 4)\n",
    "print('elements in a that are > 5: ', a[a > 5])"
   ]
  },