numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An introduction to improving Python performance with NumPy\n",
    "\n",
    "**The most important thing** when you write Python code using NumPy and vectorised libraries in general is to **think in term of vectors**. Even better for many operations is to **think in term of matrices**, i.e. 2D arrays. Well written vector and matrix operations can make good use of vector units, memory caches, and may even be multithreaded. \n",
    "\n",
    "Even the best vectorised functions can't perform vectorised operations if you pass in one value at a time. And BLAS operations on matrices may be even better optimised than the equivalent vector operation.\n",
    "\n",
    "**If you're unfamiliar with NumPy, there's a short NumPy tutorial at the bottom of this notebook**\n",
    "\n",
    "Let's compare scalar Python with vectorised NumPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from math import sin,pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting y = sin(x) using pure Python\n",
    "\n",
    "We want to plot $y = sin(x)$ for 10,000 $x$ values from -pi to pi .  First we'll do this using pure Python.\n",
    "\n",
    "The Python range function only works for integers so we have to explicitly loop over every integer and create the floating point values we need.\n",
    "\n",
    "Then, we need to loop over these x values and compute our y values.\n",
    "\n",
    "Here's how to do it using **list comprehensions**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "x = [-pi + n*pi/(N/2) for n in range(N+1)]\n",
    "y = [sin(xi) for xi in x]  \n",
    "plt.plot(x,y);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pure_python_speed = %timeit -o y = [sin(xi) for xi in x]  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll use NumPy. We can use the NumPy `linspace` function to generate the floating point values we need, and the NumPy `sin` function to generate the sines. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "x = np.linspace(-pi,pi,N)\n",
    "y = np.sin(x)\n",
    "plt.plot(x,y);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numpy_speed = %timeit -o y = np.sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('NumPy is {0:.1f} times faster than pure python for this operation'.format(pure_python_speed.best/numpy_speed.best))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elementwise addition of two vectors using pure python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following example, we create two vectors, $a$ and $b$ and add then together.  Here's one way of doing this using Python lists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -o\n",
    "# Creating and performing elementwise addition on two vectors - Using pure python\n",
    "N = 1000000\n",
    "a = [n*0.1 for n in range(N)]\n",
    "b = [n*0.2 for n in range(N)]\n",
    "c = [None] * N # Pre-allocate memory for the result\n",
    "for i in range(N):\n",
    "        c[i] = a[i] + b[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the result of the previous timeit cell. Must be run straight after the cell above\n",
    "pure_python_result = _"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -o\n",
    "# Creating and performing elementwise addition on two vectors - Using numpy\n",
    "N = 1000000\n",
    "numpy_a = np.arange(0,N/10,0.1)\n",
    "numpy_b = np.arange(0,N/5,0.2)\n",
    "numpy_c = numpy_a + numpy_b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numpy_result = _"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The NumPy version of vector addition is {0:.1f} times faster'.format(pure_python_result.best/numpy_result.best))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Where there are lists and loops - think vectors, matrices and Numpy!\n",
    "\n",
    "Whenever you see an algorithm written in Python using loops and lists, it is almost always better to switch to an array based approach using NumPy.  This array (or vectorised) based approach to scientific computing is extremely popular and forms the basis of other programming languages such as MATLAB and R where 'vectorising your code' has long been known as one of the primary ways to improve performance.\n",
    "\n",
    "In the Python world, Numpy forms the basis of a huge number of scientific and data science packages.  It is almost impossible to think of doing serious scientific computing without it.  As such, learning at least a little about it is essential."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "## Exercise 1: Use NumPy to accelerate the prime sieve\n",
    "\n",
    "Here is a pure Python implementation of the Prime Number Sieve, used in earlier training modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import  math\n",
    "def sieve_primes(n):\n",
    "    a = [True for x in range(n + 1)]\n",
    "    i = 2\n",
    "    while i <= math.sqrt(n):\n",
    "        if a[i]:\n",
    "            for j in range(i*i, n + 1, i):\n",
    "                a[j] = False\n",
    "        i += 1\n",
    "    return [i for i in range(2, len(a)) if a[i]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time how long it takes to find all primes below 10000000 and then write and test a numpy version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sieve_primes(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2: Use NumPy to accelerate option pricing\n",
    "\n",
    "The following function was ported from online MATLAB code online that claimed to price Asian call options.  I don't know much about option pricing but let's assume this claim is correct.  An initial port to Python used the standard **math** and **random** modules and looked like this.\n",
    "\n",
    "It was developed using the commonly used technique of copying and pasting the MATLAB code into Jupyter and editing it until it was valid Python code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math # This is the standard Pyton math module. Not the numpy one\n",
    "import random # This is the standard Python random module.  Not the numpy one\n",
    "\n",
    "def Asian(so,k,r,v,t,m,n):\n",
    "    \"\"\"\n",
    "    I have not identified what the arguments mean since the original MATLAB code didn't either. \n",
    "    This doc-string will be updated once I learn what they are!\n",
    "    \"\"\"\n",
    "    dt = t/m\n",
    "    AsianPayoffSum = 0\n",
    "    for i in range(1,n+1):\n",
    "        s = so\n",
    "        stSum = so\n",
    "        at = so\n",
    "        for j in range(1,m+1):\n",
    "            st = s * math.exp(((r-v**2/2)*dt) + (v*random.normalvariate(0,1)*math.sqrt(dt)))\n",
    "            stSum = stSum + st \n",
    "            at = stSum/(j+1)\n",
    "            s = st\n",
    "        AsianPayoff = max(at-k,0);\n",
    "        AsianPayoffSum = AsianPayoffSum + AsianPayoff;\n",
    "    AsianCall = math.exp(-r*t)*(AsianPayoffSum/n)\n",
    "    return(AsianCall)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is compatible with that shown by the author of the original MATLAB code.  It is common to discover Python code written like this -- typically written by those who have found something in MATLAB and who want to get it working in Python.\n",
    "\n",
    "It is slow, hence we'll use `%time` and not `%timeit`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time Asian(100,90,0.15,0.45,1,100,200000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use what you have learned so far to rewrite this function using NumPy and discover the speed-up you can reasonably expect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What to explore next\n",
    "\n",
    "Numpy is one of the fundamental libraries for data science and scientific computing. Now that you have learned how to use some of its functionality, you are ready to explore some of these other Python libraries.\n",
    "\n",
    "* [NAG Library for Python](https://www.nag.co.uk/nag-library-python) All 1900+ routines from the NAG C/Fortran library available from a Python interface.  Numpy arrays are the fundamental data type used.\n",
    "* [SciPy](https://scipy.github.io/devdocs/release.1.3.1.html) A large library of mathematical and statistical functions that operate on numpy arrays.\n",
    "* [Scikit-learn](https://scikit-learn.org/stable/) Traditional machine learning in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NumPy tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The convention for importing numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numpy is used so much that almost everyone imports it as np in order to save on typing\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating NumPy arrays manually"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From a list\n",
    "mydata = np.array([1.0,3.5,-2.657,8.2])\n",
    "mydata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ones, zeros and emptiness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ones and zeros\n",
    "np.ones(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.ones((5,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.zeros((3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Allocate memory but don't assign any values\n",
    "np.empty(12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accepted wisdom is that using `np.empty` is faster than `np.ones` because you don't have to pay the cost of initialising all those zeros.  For a 200x200 matrix, it certainly seems to be true on my machine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 200\n",
    "zero_time = %timeit -o some_zeros = np.zeros((N,N))\n",
    "empty_time = %timeit -o empty_matrix = np.empty((N,N))\n",
    "print('np.empty is {0} times faster than np.zeros'.format(zero_time.best/empty_time.best))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For N=10000, this is not the case for me.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "zero_time = %timeit -o some_zeros = np.zeros((N,N))\n",
    "empty_time = %timeit -o empty_matrix = np.empty((N,N))\n",
    "print('np.empty is {0} times faster than np.zeros'.format(zero_time.best/empty_time.best))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two approaches take the same amount of time.  When making code faster, it always pays to check your assumptions!\n",
    "Here, the reason why **empty** and **zeros** take the same amount of time is that the operating system is probably zeroing the allocated array for security reasons."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Common array constructors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A version of range that returns a numpy array and can handle floating point step sizes\n",
    "np.arange(0,1,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a linear spaced array from start to end with N points\n",
    "np.linspace(0,1,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A log-spaced array\n",
    "np.logspace(2,-1,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# An array filled with element and type of your choice\n",
    "np.full((3,3),7.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attributes of NumPy arrays\n",
    "\n",
    "A numpy array contains several attributes that tell us various things about it.  Here are a few of the most important"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.uniform(0,1,(3,3))\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.ndim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike Python lists where you can mix integers with doubles with strings, every element of a given numpy array must have the same type.  The type of a numpy array is given by its **dtype** attribute.  \n",
    "\n",
    "The default type for most operations is `float64` which is equivalent to double in C."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many numpy array creation functions have an optional **dtype** argument to select the type.  A common choice in machine learning is single precision"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.zeros((3,3),dtype='float32')\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.flags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of these attributes are not writeable. You can't change an array's size like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.size = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but you can change its shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.shape=(1,9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Construct the following 7 x 7 matrix using NumPy.\n",
    "\n",
    "$$\n",
    "\\begin{array}{rrrrrrr}\n",
    "  0 & 14 & 28 & 42 & 56 & 70 & 84 \\\\\n",
    "  2 & 16 & 30 & 44 & 58 & 72 & 86 \\\\\n",
    "  4 & 18 & 32 & 46 & 60 & 74 & 88 \\\\\n",
    "  6 & 20 & 34 & 48 & 62 & 76 & 90 \\\\\n",
    "  8 & 22 & 36 & 50 & 64 & 78 & 92 \\\\\n",
    "  10 & 24 & 38 & 52 & 66 & 80 & 94 \\\\\n",
    "  12 & 26 & 40 & 54 & 68 & 82 & 96 \\\\\n",
    "  \\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Broadcasting\n",
    "\n",
    "We can perform element wise operations mixing arrays and scalars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.randint(1,6,(1,10))\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "10*a - 2*a*a + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing in this way creates temporary matrices and it can sometimes pay to consider algebraic rearrangments of expressions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N  = 1000000\n",
    "a = np.random.random(N)\n",
    "%timeit 2*(3*a)*(5*a)\n",
    "%timeit 30*a*a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Be aware that changing the order of operations for floating point numbers may change the results by small amounts, due to rounding. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.all(2*(3*a)*(5*a) == 30*a*a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(np.abs(2*(3*a)*(5*a) - 30*a*a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We aren't limited to broadcasing scalars.  Consider the following"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.randint(1,6,(1,9))\n",
    "a.shape=(3,3)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([[101,102,103]]).T\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do you think the following would do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a+b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This behaviour is called broadcasting because the vector b was broadcasted across the matrix a. \n",
    "\n",
    "It's conceptually similar to the following but without the explicit creation of a tempory matrix like tiled_b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiled_b = np.tile(b,(1,3))\n",
    "tiled_b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a + tiled_b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Slicing\n",
    "\n",
    "If you have used MATLAB extensively, you are probably used to operations that look like slicing but with different syntax.  It's the standard way of indexing into NumPy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(1,11)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[::2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[2:8:1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[2:8:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Something that you really need to understand is that slices are just **views** into the original array. **They are not copies**. \n",
    "Perhaps the following helps illustrate the behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a slice called b that looks at all odd elements of the array\n",
    "b=a[::2]\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set every element of b to 10\n",
    "b[:] = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#What does a look like?\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Masking (Logical Indexing)\n",
    "\n",
    "If you have are used to MATLAB, you will have come across a technique called logical indexing.  In NumPy it is called masking."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.randint(1,10,(1,10))\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How might we sum all the elements that are bigger than, say, 4?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logical = a>4\n",
    "logical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[logical]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[logical].sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively done in one line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[a>4].sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### All and any\n",
    "\n",
    "Two commands that are very useful when dealing with boolean arrays are all and any"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "some = np.array([1,0,0,0,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.all(some)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.any(some)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Array-based mathematics\n",
    "\n",
    "There are several functions that prove to be useful when turning loops into vectorised code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(1,7,0.5)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.cumsum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.cumprod(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.argmin(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.std(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercises\n",
    "\n",
    "Create a 3x6 array of booleans containing all false values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the array \n",
    "\n",
    "`a = np.arange(20)`\n",
    "\n",
    "replace every element that's a multiple of 3 with 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct the following matrix\n",
    "\n",
    "$$\n",
    "\\begin{array}{rrrr}\n",
    "  2 & 0 & 0 & 0 \\\\\n",
    "  0 & 2 & 0 & 0  \\\\\n",
    "  0 & 0 & 2 & 0  \\\\\n",
    "  0 & 0 & 0 & 2  \\\\\n",
    "  \\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct the following matrix\n",
    "\n",
    "$$\n",
    "\\begin{array}{rrrr}\n",
    "  0 & 1 & 0 & 0 \\\\\n",
    "  4 & 0 & 1 & 0  \\\\\n",
    "  0 & 4 & 0 & 1  \\\\\n",
    "  0 & 0 & 4 & 0  \\\\\n",
    "  \\end{array}\n",
    "$$"
   ]
  },
  {
   "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.3"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}