LMB_Python_Numeric_2023withAnswers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numeric and scientific Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NumPy arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though mathematical operations can be performed in regular Python the NumPy module is often faster and more convenient. Working with NumPy involves using n-dimensional array objects which store regular arrays of data of a specified type (usually, but not limited to, numeric types). Such arrays can be used in a similar way to lists: they contain an ordered sequence of values and can be used in loops etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%autosave 0\n",
    "%matplotlib inline\n",
    "import numpy as np  # Load the NumPy module, assign it the name \"np\" for convienence\n",
    "\n",
    "l = [1,4,9,16,25] # A list\n",
    "a = np.array(l)   # An array built from a list\n",
    "\n",
    "print(l)  #  Printed with commas\n",
    "print(a)\n",
    "print(type(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, an important idea with NumPy arrays is that operations can be performed on the array as a whole, rather than looping though all the component elements. This means that a fast internal implementation can be used. Also, it results in syntax similar to matrix algebra, where each variable is an entire array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plain Python : loops\n",
    "a = [1,2,3,4,5]\n",
    "b = []\n",
    "for x in a:\n",
    "    b.append(x*3)\n",
    "    \n",
    "# List comprehension    \n",
    "c = [x*x for x in a]\n",
    "\n",
    "print(a)\n",
    "print(b)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numpy whole-array operations\n",
    "a = np.array([1,2,3,4,5])\n",
    "b = 3 * a   # All elements multiplied by three\n",
    "c = a * a   # All elements squared\n",
    "\n",
    "print(a)\n",
    "print(b)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An array will contain elements of the same data type. This data type is determined from the contents of the array when it is constructed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array([2, 5, 1, 8, 0])\n",
    "z = np.array([3.14159, 2.71828, 1.41421])\n",
    "\n",
    "print(y.dtype)  # int - whole numbers\n",
    "print(z.dtype)  # float - fixed precision real numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data type of an array may be specifically stated (i.e. forced), irrespective of its initial elements. Also, the data type may be converted (making a new array in the process) using `astype()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array([2, 5, 1, 8, 0], float)  # Force float dtype\n",
    "x = z.astype(int)                     # Convert to ints\n",
    "\n",
    "print(y, y.dtype)\n",
    "print(z)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many common operations are applied in an element-wise manner, i.e. to each value individually, and operations can work between two arrays if they have compatible sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([1.0, 2.0, 3.0])\n",
    "y = np.array([3.0, 4.0, 5.0])\n",
    "\n",
    "print(x + 2)    # Add 2 to all values\n",
    "print(y / 2)    # Divide all values by 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Element-wise operations bewtween different arrays with compatible sizes\n",
    "print(x + y)   \n",
    "print(x * y)   \n",
    "print(x - y)\n",
    "print(x / y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 1: Fahrenheit to Celsius</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>Convert an array of temperatures in degrees Fahrenheit to degrees Celsius: subtract 32 and divide by 1.8.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code\n",
    "temps = [-30,0,32,70,98,212]\n",
    "f = np.array(temps)\n",
    "c = (f-32.0)/1.8\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-dimensional arrays\n",
    "\n",
    "An array can have a number of different dimensions/axes, i.e. so that it can represent vectors, matrices tensors etc. For example a 2D array, with 'row' and 'column' axes could be created as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[1,2,3],\n",
    "              [4,5,6]])  # Make 2D array from list of lists\n",
    "\n",
    "print(x)\n",
    "print(x.shape)           # (2,3) - rows x columns\n",
    "print(x.size)            # 6 - elements in total\n",
    "print(x.ndim)            # 2 - two axes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And similarly for a 3D array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array([[[0,1], [2,3]],\n",
    "              [[4,5], [6,7]],\n",
    "              [[8,9], [10,11]]])\n",
    "print(y)\n",
    "print(y.shape)\n",
    "print(y.ndim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of array axes can be forced to be larger than what is automatically suggested by the input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.array([0,1,4,9], ndmin=2)  # Forces a 2D array with one row\n",
    "\n",
    "print(z)\n",
    "print(z.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operations are generally applied in an element-wise manner even when an array has multiple axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(x)\n",
    "print(x * 0.01)  # Element-wise operation on 2D array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Array creation from differently sized sub-collections is not allowed in more recent versions of NumPy. (In older versions it works, but gives a 1D array of Python objects, which may not be the expected result...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[0,1], [2,3,4], [5,6,7,8]]) # Fails! Inconsistent 'inhomogeneous' shape."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arrays are easily reshaped to change the number of rows, columns axes etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(1,101)     # The NumPy equivalent of range()\n",
    "y = x.reshape(10,10)     # Same data as 10 rows x 10 columns (makes a new array)\n",
    "\n",
    "print(x.shape, y.shape)\n",
    "print(x.size, y.size)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filled arrays\n",
    "\n",
    "There are various functions to create particular kinds of filled arrays (e.g. all zeros), with a specified size/shape and a specified data type for the elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.zeros((2,3))     # 2 x 3 array full of 0.0\n",
    "b = np.ones((3,2), int) # 3 x 2 array full of 1\n",
    "c = np.full(4, 7.0)     # array, length four, full of 7.0\n",
    "d = np.identity(3)      # 3 x 3 identity matrix\n",
    "\n",
    "print(a)\n",
    "print(b)\n",
    "print(c)\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, arrays can be created containing regular ranges of numbers, using start, end and step specifications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(1.2, 4.4, 0.2)  # From 1.2 to <4.4 in steps of 0.2\n",
    "\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.linspace(1.2, 4.4, 5)  # From 1.2 to 4.4 in five even steps\n",
    "z = np.logspace(2, 7, 5, base=10) # From 100 (10^2) to 10 million (10^7) in 5 even log_10 steps\n",
    "\n",
    "print(y)\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 2: Checkerboard</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>Using modular arithmetic (e.g. `%` operator) create a two-dimensional, integer array of size 5x5 with alternating ones and zeros.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code\n",
    "a = np.arange(25).reshape(5,5)\n",
    "b = a % 2\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Broadcasting\n",
    "\n",
    "The shape of two arrays determines whether they are compatible for operations that combine their elements. As shown above, arrays are compatible if they are the same shape/size. However, they may also be compatible if their last axes are the same size. Here a 1D array of size 2 (`y`) can be added to a 2D array of size 3 x 2 (`x`) because the 1D size matches the number of 2D columns, so the operation can be applied separately (\"broadcast\") to each row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[2,3], [4,5], [6,7]])\n",
    "y = np.array([10,100])\n",
    "\n",
    "print(x.shape, y.shape)\n",
    "print(x+y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Naturally, if an operation is applied to arrays with incompatible shapes, an error is triggered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([2, 4, 6, 8, 10])\n",
    "y = np.array([1, 3, 5])\n",
    "\n",
    "print(x*y) # Incompatible shapes : last axis different sizes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, it is sometimes possible to reshape arrays so they become compatible. In the next example an extra last axis, with size 1, is added to an otherwise incompatible array. When this is multiplied by the 1D array `x`, each column of `y` (only a single value) can be multiplied separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([2, 4, 6, 8, 10])\n",
    "y = np.array([1, 3, 5])\n",
    "\n",
    "y = y[:,None]      # Add a new axis of size 1 : .reshape() could also be used here\n",
    "\n",
    "print(x.shape)\n",
    "print(y.shape)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(x*y) # Shapes now compatible : operation is spread along last axis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Array sub-selections\n",
    "\n",
    "There is an index and range (slicing) syntax of the form `start:limit:step`, i.e. similar to the system used with lists and strings etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(10) ** 2   # Sequential ints, squared\n",
    "\n",
    "print(x)\n",
    "print(x[2])         # number at index 2\n",
    "print(x[1:4])       # slice range (makes new array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ranges may also have a third argument to specify the increment (step) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(x[2:9:2])   # Start at index 2, increment 2\n",
    "print(x[::2])     # Every other element (start:end is implicit)\n",
    "print(x[::-1])    # Negative increments means backwards (end:start implicit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The indexing and slicing syntax for arrays with more than one axis uses a comma, such as `data[i,j]` for a 2-dimensional array. This differs from regular Python where separate brackets are needed for each sub-list, e.g. using `data[i][j]`. Accessing with multiple brackets will still work with NumPy arrays, but will be slower, as it makes an intermediate array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.arange(25).reshape(5,5)  # Two dimensional array : 5 rows x 5 columns\n",
    "print(y)\n",
    "print(y[2,3])     # One element specified with [Row, Column]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For multi-dimensional arrays, indices and slices/ranges may be specified for each axis to select sub-arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(y[1:4,1:4]) # A range of rows and columns (makes a new array)\n",
    "print(y[1,0:5])   # One entire row: all columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entire axis ranges may be selected using `:`, but this is implicit for the last axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(y[-1])      # The entire last row: columns implicit\n",
    "print(y[-1,:])    # The entire last row, again \n",
    "print(y[:,2])     # One column (all rows)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tuples of indices may be specified to make selections from an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = (0,2,4)      # Row/column indices\n",
    "print(y[idx,:])    # Select specific rows, all columns\n",
    "print(y[idx, idx]) # Select row, column pairs : (0,0) (2,2) (4,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ranges and explicit indices can be used to modify the contents of an array, as well as for selecting sub-parts shown above. For example the below 5x5 array has its central 3x3 elements modified, firstly to a single value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.arange(25).reshape(5,5)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y[1:4, 1:4] = 0 # Set the central ranges of the 2D array to be zero \n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next an internal row is set in its entirity to a corresponsing sequence of numbers (of matching size):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y[2] = range(-5,0) # Set all of row with index 2 to new values.\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 3: Power columns</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>Create an array containing the first five powers (i.e. $x$, $x^2$ ... $x^5$ as columns) for integers in range [1, 10]. Extract the third and last columns (for all rows) as a separate array.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code : Hint use .reshape() to broadcast two arrays.\n",
    "x = np.arange(11)\n",
    "p = np.arange(1,6)\n",
    "x = x[:,None]\n",
    "\n",
    "y = x ** p\n",
    "print(y)\n",
    "z = y[:,(2,-1)]\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mathematical operations\n",
    "\n",
    "NumPy provides various mathematical functions, many of which are named like those in the standard Python `math` module. These operate on all the elements of an array, though they also work on single numbers. The functions will also accept regular iterable Python objects (lists or tuples etc.) as input, but a NumPy array is created where appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = range(1,5)\n",
    "\n",
    "print(np.sqrt(x))   # Square root\n",
    "print(np.exp(x))    # E to the power\n",
    "print(np.log(x))    # Natural log\n",
    "print(np.mod(x, 2)) # Modular arithmetic : remainder in base 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a number of functions to round numbers up, down or to a specified number if decimal places."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array([3.14159, 2.71828, 1.41421, -0.70717])\n",
    "\n",
    "print(np.abs(y))      # absolute/positive value\n",
    "print(np.round(y, 3)) # round to 3 d.p.\n",
    "print(np.ceil(y))     # whole number above (as float type)\n",
    "print(np.floor(y))    # whole number below (as float type)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is notable that integer rounding of half integer values uses *convergent rounding* where these edge values get rounded to the nearest even number to avoid statistical bias (this is actually the same as plain Python)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.round([0.5, 1.5, 2.5, 3.5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy has various trigonometric functions, though it is notable that these are sometimes named differently to the `math` module equivalent. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.radians([30.0, 60.0, 90.0, 135.0]) # Convert degrees to radians\n",
    "print(z)\n",
    "c = np.cos(z)        # Cosine \n",
    "print(c)\n",
    "print(np.arccos(c))  # Inverse cosine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inbuilt array methods\n",
    "\n",
    "There are various methods inbuilt into an array object, for example to calculate sums, extrema, mean and standard deviation. Without extra arguments these consider all elements of the array, irrespective of array shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[3,9], [3,1], [5,7]])\n",
    "\n",
    "print(x.min())          # minimum value (of all elements)\n",
    "print(x.max())          # maximum value\n",
    "print(x.sum())          # summation\n",
    "print(x.mean())         # mean\n",
    "print(x.std())          # standard devaiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These methods accept an `axis` specification so that the operation is performed multiple times along that array direction, thus generating a new array. It might be semantically confusing as to which axis number should be used to operate over rows or columns etc. For example, in a 2D array to get the summation for all rows you do `arr.sum(axis=1)`. The key idea here is that the sums are done over the column axis (`=1`) to give a value for each row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(x.max(axis=0))    # max along row axis : a value for each column\n",
    "print(x.sum(axis=1))    # sum along column axis : a value for each row\n",
    "print(x.mean(axis=1))   # mean along column axis : a value for each row"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 4: Mean angle</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>Calculate the mean angle of the below list of angles in degrees, noting that you cannot take a simple average; convert to radians and find the inverse tangent of the average sine and cosine using `np.arctan2()`</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code\n",
    "angles = [6, 45, 90, 262, 355]\n",
    "rads = np.radians(angles)\n",
    "s = np.sin(rads)\n",
    "c = np.cos(rads)\n",
    "m = np.arctan2(s.mean(), c.mean())\n",
    "m = np.degrees(m)\n",
    "print(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boolean arrays and selection\n",
    "\n",
    "Comparison operations between arrays generate Boolean arrays giving `True` or `False` for each elemental comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([1, 8, 9, 3])\n",
    "y = np.array([5, 0, 3, 7])\n",
    "\n",
    "z = x > y\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Element selections from arrays can be made by passing such arrays of `True` and `False` values (for each axis). This technique of using a 'mask' extends the index (row, column) based slection already shown."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(x)\n",
    "print(z)\n",
    "print(x[z])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Boolean masks can also be used to set a subset of array elements. For example here the negative values are set to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([3,-5,1,-2,-7,8,0])\n",
    "y = x.copy() # A copy of the original x\n",
    "z = x < 0    # Z has True for negative elements of x\n",
    "\n",
    "print(x)\n",
    "print(z)     # Boolean mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[z] = 0     # Negative elements in x set to zero\n",
    "\n",
    "print(y)     # Original\n",
    "print(x)     # Modified"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 5: Distance mask</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>For the given 2D array of coordinates (e.g. $(x,y)$) create a Boolean mask to select all the points at least 1.0 unit radius away from the orgin (0.0, 0.0). Make a new array containing only these points.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code\n",
    "a = np.array([[-0.49,-0.14], [ 0.64,-0.47], [1.17, 0.55], [-0.56, 0.82], [1.22, 0.61],\n",
    "              [ 0.84, 0.81], [-0.89,-0.50], [0.89, 0.44], [ 0.57,-0.82], [1.31, 0.61],\n",
    "              [ 0.97, 1.40], [ 1.37, 0.55], [1.69, 0.10], [ 0.59, 1.76], [0.01,-0.97]])\n",
    "\n",
    "sq = a * a\n",
    "dists = np.sqrt(sq.sum(axis=1))\n",
    "selection = dists >= 1.0\n",
    "print(selection)\n",
    "\n",
    "b = a[selection]\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Concatenating arrays\n",
    "\n",
    "Sometimes there is a need to combine the data in NumPy arrays, to join them together into a larger array. There are several NumPy functions to do this, depending on which axes of a multidimensional array should be linked. The `hstack()`, `vstack()` and `dstack()` join arrays along the horizontal (first), vertical (second) and depth (third) axes respectively. It is notable that the arrays to be joined are passed to these functions in a collection (usually a list or tuple) and make a new array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([1,2,3,4])\n",
    "y = np.array([7,8,9,0])\n",
    "\n",
    "h = np.hstack([x, y]) # Join along column axis\n",
    "\n",
    "print(h, h.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = np.vstack([x, y]) # Join along row axis\n",
    "\n",
    "print(v, v.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are also similarly named `hsplit()`, `vsplit()` and `dsplit()` which perform the opposite functions. Though it is notable that these function require a second argument to say how many equal-sizd arrays to split into.s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[1,11],\n",
    "              [2,12],\n",
    "              [3,13],\n",
    "              [4,14]])\n",
    "\n",
    "a, b = np.hsplit(x, 2) # Separate columns\n",
    "\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c, d = np.vsplit(x, 2) # Separate rows (two segments)\n",
    "\n",
    "print(c)\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `concatenate()` function is more general and a specific axis number, to join on, can be stated.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 x 2 input arrays\n",
    "x = np.array([[1,2],\n",
    "              [3,4]])\n",
    "y = np.array([[7,8],\n",
    "              [9,0]])\n",
    "\n",
    "# Join along rows (same as hstack)\n",
    "c = np.concatenate([x, y], axis=0) \n",
    "\n",
    "print(c)\n",
    "\n",
    "# Join columns\n",
    "d = np.concatenate([x, y], axis=1) \n",
    "\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Index selection methods\n",
    "\n",
    "Sometimes it is handy to convert from a Boolean array to indices, e.g. to know which (row, column) values are true. For this `nonzero()` can be used get a tuple representing the positions of the true (non-zero) elements. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array(range(10))\n",
    "y = x % 2 == 1           # True if odd\n",
    "\n",
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = y.nonzero()  # Positions where y is true\n",
    "\n",
    "print(idx)         # A tuple of arrays (one for each axis)\n",
    "print(x[idx])      # Values in y where z was true"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a similar manner `argsort()` can be used to get the indices of an array in value order, which can then be used for sorting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array([9,2,7,1,5,3,6,1])\n",
    "idx = y.argsort()  # An array of element positions, from smallest to largest value\n",
    "\n",
    "print(idx)\n",
    "print(y[idx])      # Values from y at sorted positions : makes a sorted array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here it is notable that NumPy also has a more simple `sort()` function that gives a sorted copy of an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.sort(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also `argmax` and `argmin` can be used to get the index positions of extrema."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array([-8,2,7,-9,5,3,6,11])\n",
    "idx = y.argmin()   #  Position of min value\n",
    "\n",
    "print(idx)\n",
    "print(y[idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[1,7,2],\n",
    "              [0,1,9]])\n",
    "idx = x.argmax(axis=1)  # Column pos. of max. for each row\n",
    "\n",
    "print(x)\n",
    "print(idx) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 6: Filtered indices</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>For 20 evenly spaced angle values in the range 0 to $2\\pi$, print the indices where the absolute value of the cosine is more than 0.5. Also, print the corresponding angles for these indices in degrees.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code\n",
    "\n",
    "angles = np.linspace(0.0, 2*np.pi, 20)\n",
    "c = np.cos(angles)\n",
    "idx = (c>0.5).nonzero()\n",
    "\n",
    "print(c)\n",
    "print(idx)\n",
    "print(np.degrees(c[idx]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear algebra\n",
    "\n",
    "There are many functions that relate to matrix operations and linear algebra, for example to calculate various products (inner, outer, cross), to transpose and find inverse matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array(((1,1),(1,0)))\n",
    "y = np.array(((0,1),(1,1)))\n",
    "\n",
    "print(np.dot([1,2,3], [4,5,6]))   # scalar product\n",
    "print(np.dot(x, y))               # matrix multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[1, 2, 3], [4, 5, 6]])\n",
    "\n",
    "print(x)\n",
    "print(x.transpose())   # swap rows with cols : x.T also used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array(((1,1),(1,0)))\n",
    "\n",
    "print(np.linalg.inv(y))  # inverse transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NumPy sub-modules\n",
    "\n",
    "NumPy has a number of handy submodules such as the very useful `random` module for pseudo-random number generation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import uniform, normal # Import from sub-module\n",
    "\n",
    "# Sample uniform distribution\n",
    "a = uniform(0.0, 2.0, (5,5)) # From [0, 2), as 5 by 5 array\n",
    "\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample normal distribution\n",
    "b = normal(0.0, 2.5, 10)  # Mean 0.0, SD 2.5, 10 values\n",
    "\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also the `fft` module for Fourier transforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import pi\n",
    "\n",
    "l = 0.04  # Decay\n",
    "w = 0.1   # Frequency\n",
    "\n",
    "t = np.arange(0.0, 30.0, 1.0)         # Time values\n",
    "x = np.exp(2j*pi*w*t) * np.exp(-l*t)  # Wave equation using complex numbers\n",
    "y = np.fft.fft(x)                     # Fast Fourier transform array of wave\n",
    "\n",
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting with matplotlib\n",
    "\n",
    "NumPy installations include the `matplotlib` module, which provides the ability to make graphs. For example, plotting the above using the `plot()` function makes a line graph. The way this operates is that datasets are added individually before the completed plot is shown at the end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (16,8) # Set plot size\n",
    "\n",
    "plt.plot(x.real)  #  Real values from complex number\n",
    "plt.plot(x.imag)  #  Imaginary values\n",
    "plt.show()  #  Display on-screen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is notable that invoking the `show()` funtion clears all of the plotted data,."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.show()  # No data plotted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many options for controlling the presentation and style of the plots. A few examples are illustrated below to add an axis label, a figure legend and line style. For full options see the [matplotlib documentation](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.plot.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot same data as before\n",
    "plt.plot(x.real, color='red', label='Real')\n",
    "plt.plot(x.imag, color='#0080FF', linewidth=5,\n",
    "         alpha=0.3, label='Imaginary')\n",
    "\n",
    "plt.xlabel('Time (s)')\n",
    "plt.legend()           # Add legend : made with labels\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As well as on-screen presentation, all the charts may be saved to a file using `savefig()` and specifying the save location. There are many file formats available including pixel formats such as PNG and JPEG and vector graphics such as PDF or SVG. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x.real, color='#FF4000')\n",
    "plt.plot(x.imag, color='#0080FF')\n",
    "\n",
    "# file format is specified in the file name extension\n",
    "plt.savefig(\"TestGraph.png\", dpi=72, bbox_inches='tight')  # .png : PNG pixmap\n",
    "plt.savefig(\"TestGraph.pdf\", dpi=100, bbox_inches='tight') # .pdf : PDF "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Several types of chart and graph are available. For example the `scatter()` function makes a scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_vals = 2000\n",
    "x_vals = np.arange(num_vals)\n",
    "norm_vals = np.random.normal(0.0, 1.0, num_vals) # Mean, STD, count\n",
    "\n",
    "plt.scatter(x_vals, norm_vals, s=4, alpha=0.5)  # size 4 spots\n",
    "plt.scatter(x_vals, np.sort(norm_vals), s=4, alpha=0.5)  # Height order\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another handy plot uses `hist()` to make histogram of the data (this has many options from `numpy.histogram`). Here, we show a histogram of the same random, normally distributed data used above: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Histogram with set number of bins in a set range\n",
    "plt.hist(norm_vals, bins=50, range=(-3.5, 3.5))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Different graph types can be mixed together in the same plot if required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(norm_vals, color='red')\n",
    "\n",
    "x = np.arange(-3.5, 3.6, 0.1)\n",
    "y = num_vals/4 * np.exp(-x*x/2)\n",
    "\n",
    "plt.plot(x, y, color='black')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 7: 2D Gaussians</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>Make a scatter plot for two random two-dimensional datasets that sample normal distributions: make one with $\\mu$=0.0, $\\sigma$=1.0 (for both axes) and the other with $\\mu$=1.0, $\\sigma$=0.5 (for both axes). Plot 1000 2D points for each and use `alpha=0.5` to add transparency to the `scatter()`.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code - Hint: use np.random.normal(mean, stdev, shape)\n",
    "s = (2, 1000)\n",
    "d1 = np.random.normal(0.0, 1.0, s)\n",
    "d2 = np.random.normal(1.0, 0.5, s)\n",
    "x1, y1 = d1\n",
    "x2, y2 = d2\n",
    "\n",
    "plt.scatter(x1, y1, alpha=0.5)\n",
    "plt.scatter(x2, y2, alpha=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sub-plots\n",
    "\n",
    "So far we have generated plots with a single chart. However, Matplotlib has the ability to embed several sub-plots within a larger figure. There are a few ways to do this, but one of the easiest is via the `pyplot.subplots()` function. This takes a number of rows and columns as input and outputs an array of sub-plots; here `sp`. The sub-plot objects can be used in a similar manner to that already shown for making the graphs etc. Note that the `fig` variable is a reference to the larger figure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, sp = plt.subplots(1, 3) # Rows, columns : Gives figure and sub-plot array\n",
    "\n",
    "# Random values for two dims\n",
    "n_points = 2000\n",
    "x_vals = np.random.gamma(5.0, 1.0, n_points)\n",
    "y_vals = np.random.normal(0.0, 1.0, n_points)\n",
    "\n",
    "sp[0].scatter(x_vals, y_vals, alpha=0.5, s=3)    # Simple scatter\n",
    "sp[0].set_title('Scatter')\n",
    "\n",
    "sp[1].hist2d(x_vals, y_vals, bins=(20,20))       # 2D histogram\n",
    "sp[1].set_title('Hist2D')\n",
    "\n",
    "sp[2].hexbin(x_vals, y_vals, gridsize=20)        # Hex-bin histogram\n",
    "sp[2].set_title('Hexbin')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 8: Sines and cosines </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>Using the below arrays of `x` and `y` values, create a plot with two sub-plots. In the first sub-plot make a line plot of `x` versus `y`. In the second sub-plot create a 2D array of all `x` multiplied by all `y` (you can use broadcasting or `np.outer()`) and display this using the `.matshow()` function</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code\n",
    "a = np.arange(-2*np.pi,2*np.pi,0.05)\n",
    "x = np.sin(a)\n",
    "y = np.cos(a)\n",
    "fig, (ax1, ax2) = plt.subplots(2)\n",
    "\n",
    "ax1.plot(x, y)\n",
    "\n",
    "m = x[:,None] * y[None,:] # or np.outer(x,y)\n",
    "ax2.matshow(m)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example : Mandelbrot fractal\n",
    "\n",
    "The following example illustrates the use of several previously array creation and manipulation techniques and applies them inside a Python function. Here the objective is to generate a Mandelbrot fractal; constructed from the number of iterations it takes for the iterative process *z = z<sup>2</sup> + c* to grow large, for different complex values of *c*. The function has mandatory input arguments `x0`, `x1`, `y0` and `y1` to specify a 2D numeric range, an optional `step` to specify granularity and lastly `maxiter` and `limit` to control the fractal generation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mandelbrot(x0, x1, y0, y1, step=5e-3, maxiter=64, limit=4.0):\n",
    "    \n",
    "    x_vals = np.arange(x0, x1, step)\n",
    "    y_vals = np.arange(y0, y1, step) * 1j  # Y is complex\n",
    "    \n",
    "    z = x_vals + y_vals[:,None]            # Complex number grid\n",
    "    c = z.copy()\n",
    "\n",
    "    out = np.zeros(z.shape) # Iterations before Z becomes large\n",
    "    \n",
    "    for i in range(maxiter):\n",
    "        keep = abs(z) < limit          # Bool mask for small values\n",
    "        z[keep] = z[keep]**2 + c[keep] # Update subset\n",
    "        out[keep] = i                  # Record last iteration\n",
    "   \n",
    "    return out "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inside the function the first lines are set-up for two grids of complex numbers: `c` will stay constant and `z` is updated in an iterative loop. The iteration is then done in a `for` loop where small values of `z` are updated and the number if iterations reached is recorded in `out`. The output 2D array is then readily visualised in a graphical form using `pyplot.matshow()`. This uses a named colour scheme to map numeric values (e.g. intensity) to colours. For a full list of schemes see the [matplotlib colormap documentation](https://matplotlib.org/users/colormaps.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = mandelbrot(-2.2, 1.0, -1.5, 1.5)\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.matshow(m, cmap='RdYlBu')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 9: Julia set</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>Based on the `mandelbrot()` example above make a new function that creates a Julia set fractal. Do this by changing the complex array `c` to a single, specified complex number. Plot the fractal for different values of `c`.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code\n",
    "\n",
    "def julia(c, x0=-2, x1=2, y0=-1.5, y1=1.5, step=5e-3, maxiter=64, limit=4.0):\n",
    "\n",
    "    x_vals = np.arange(x0, x1, step)\n",
    "    y_vals = np.arange(y0, y1, step) * 1j  # Y is complex\n",
    "    \n",
    "    z = x_vals + y_vals[:,None]            # Complex number grid\n",
    " \n",
    "    out = np.zeros(z.shape) # Iterations before Z becomes large\n",
    "    \n",
    "    for i in range(maxiter):\n",
    "        keep = abs(z) < limit          # Bool mask for small values\n",
    "        z[keep] = z[keep]**2 + c       # Update subset\n",
    "        out[keep] = i                  # Record last iteration\n",
    "\n",
    "   \n",
    "    return out\n",
    "    \n",
    "j1 = julia(0.8j)\n",
    "j2 = julia(-0.4-0.6j)\n",
    "\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.matshow(j2, cmap='RdYlBu')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SciPy\n",
    "\n",
    "SciPy is an extensive library that builds upon the NumPy arrays and their functions. It provides specialized scientific functionality for areas including further linear algebra, optimization, integration, interpolation, signal processing, image processing and differential equations. These are the modules listed at [scipy.org](https://docs.scipy.org):\n",
    "\n",
    "Subpackage |\tDescription\n",
    "-----------|----------\n",
    "cluster \t|Clustering algorithms\n",
    "constants \t|Physical and mathematical constants\n",
    "fftpack \t|Fast Fourier Transform routines\n",
    "integrate \t|Integration and ordinary differential equation solvers\n",
    "interpolate \t|Interpolation and smoothing splines\n",
    "io \t|Input and Output\n",
    "linalg \t|Linear algebra\n",
    "ndimage \t|N-dimensional image processing\n",
    "odr \t|Orthogonal distance regression\n",
    "optimize \t|Optimization and root-finding routines\n",
    "signal \t|Signal processing\n",
    "sparse \t|Sparse matrices and associated routines\n",
    "spatial \t|Spatial data structures and algorithms\n",
    "special \t|Special functions\n",
    "stats \t|Statistical distributions and functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistical distributions\n",
    "\n",
    "SciPy provides objects representing random variables for a variety of different statistical distributions. These are created by specifying the particular parameters for the distribution (e.g. mean and standard deviation for normal/Gaussian). Here normal and poisson distributions are illustrated. A random variable object is created by specifying the parameters for each distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm, poisson\n",
    "\n",
    "# Make \"random variable\" objects with params\n",
    "norm_rv = norm(1.0, 0.5)   # (Mean, STD)\n",
    "pois_rv = poisson(1.5)     # (Event_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From these objects we access the probability density function `.pdf()` for the (continuous) normal distribution and the probability mass function `.pmf()` for the Poisson distribution. These generate the probabilities of obtaining the input values (here `x`), according to the probability distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(0, 5, 0.1)  # Value to plot probs for\n",
    "\n",
    "# Get probabilites from the distributions\n",
    "y_norm = norm_rv.pdf(x)\n",
    "y_pois = pois_rv.pmf(x)  # Only works at discrete, integer values\n",
    "\n",
    "# Make line plots\n",
    "plt.plot(x, y_norm, label='Normal')\n",
    "plt.plot(x, y_pois, label='Poisson') \n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The random variable objects have many more associated methods. These are listed in detail in the [SciPy .stats documentation](https://docs.scipy.org/doc/scipy/reference/stats.html). Taking the normal distribution as an example we can calculate the cumulative density function with `.cdf()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot cumulative probability \n",
    "plt.plot(x, norm_rv.cdf(x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also notable that the random variable's module often has a `.fit()` method to estimate statistical parameters from data. Here, for example, we use `norm` (without parenthesis and specific parameters) to estimate mean and standard deviation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "y = [-0.516, -0.486, 0.495, -0.07, 0.974,\n",
    "     1.27, 1.179, 0.458, 0.815, 1.163]\n",
    "\n",
    "mean, std = norm.fit(y)\n",
    "print(mean, std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next example uses the random variable objects in a different way, wrapping them in a function to perform a tailed-test, as would be done to estimate a p-value: the probability of obtaining a value (or more extreme) from a given random distribution with stated parameters. \n",
    "\n",
    "The function `normal_tail_test()` performs the test on an input array of numbers for a normal distribution with given mean value, `mv` and standard deviation `std`. There is an option to state if we want to do a one- or two-tailed test, i.e. consider only values on the same side of the mean or both sides."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normal_tail_test(values, mv, std, one_sided=True):\n",
    "  \n",
    "    norm_rv = norm(mv, std)        # Normal distrib. object\n",
    "    diffs = abs(values-mv)         # Pos. separations from mean\n",
    "    integ = norm_rv.cdf(mv-diffs)  # CDF : integral of initial tail region\n",
    "  \n",
    "    if not one_sided:   # Two-tailed test\n",
    "        integ *= 2      # Symmetric: double area\n",
    "\n",
    "    return integ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function can be tested given some parameters and test values. In this case the values could represent the heights of people."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean    = 1.76\n",
    "std_dev = 0.075\n",
    "values  = np.array([1.8, 1.9, 2.0])\n",
    "\n",
    "result = normal_tail_test(values, mean, std_dev, one_sided=True)\n",
    "print('Normal one-tail p-vals:', result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we plot the curve integral that was considered for computing the p-value at the `1.9` point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rv = norm(mean, std_dev)\n",
    "x_vals = np.arange(1.5,2.2,0.01)\n",
    "plt.plot(x_vals, rv.pdf(x_vals))\n",
    "\n",
    "x_vals2 = np.arange(1.9,2.2,0.01)\n",
    "plt.fill_between(x_vals2, rv.pdf(x_vals2), alpha=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting\n",
    "\n",
    "While there are modules to do specialised fitting for statistical distributions, as shown above, SciPy has more general numerical fitting/optimization routines that can be used with arbitrary functions. We can illustrate this using an example that tries to fit a scaled bimodal Gaussian distribution. Out test data is constructed from randomly sampling two different normal distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = np.concatenate([np.random.normal(0.0, 1.0, 10000),\n",
    "                       np.random.normal(1.5, 0.5, 5000)])\n",
    "\n",
    "y_vals, x_edges = np.histogram(test, bins=100, range=(-4.0, 4.0))\n",
    "half_width = 4.0/100.0\n",
    "x_vals = x_edges[:-1] + half_width\n",
    "\n",
    "plt.plot(x_vals, y_vals)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To fit this we first create a function `bimodal_norm()` which maps data to the probability density curve for a bimodal distribution using specified parameters. Here the parameters will be the amplitudes, means and standard deviations for the two underlying, component curves. Internally this function simply uses `norm.pdf()` twice and adds the scaled result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "# Bimodal scaled normal distributions\n",
    "def bimodal_norm(data,mu1,sigma1,amp1,\n",
    "                      mu2,sigma2,amp2):\n",
    "\n",
    "    v1 = amp1 * norm.pdf(data, mu1, sigma1)\n",
    "    v2 = amp2 * norm.pdf(data, mu2, sigma2)\n",
    "    return v1 + v2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then let SciPy fit our test data using this function via `.curve_fit()`, and vire the best-fit parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "params, cov = curve_fit(bimodal_norm, x_vals, y_vals)\n",
    "print(params[:3])\n",
    "print(params[3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the fitted parameters to visualise the optimised bimodal and component curves in relation to the initial input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1,sigma1,amp1,mu2,sigma2,amp2 = params\n",
    "\n",
    "# Optimized bimodal\n",
    "d0 = bimodal_norm(x_vals, mu1,sigma1,amp1,mu2,sigma2,amp2)\n",
    "# Separate normal distribs\n",
    "d1 = amp1 * norm.pdf(x_vals, mu1, sigma1)\n",
    "d2 = amp2 * norm.pdf(x_vals, mu2, sigma2)\n",
    "\n",
    "plt.plot(x_vals, d0)\n",
    "plt.plot(x_vals, d1)\n",
    "plt.plot(x_vals, d2)\n",
    "plt.scatter(x_vals, y_vals, alpha=0.5) # Orginal data\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Image data\n",
    "\n",
    "SciPy has various handy functions to manipulate image data sorted as a NumPY ND array.\n",
    "To illustrate some of these an example image file will be downloaded from Wikipedia using Python's inbuilt `urllib` module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from urllib.request import urlopen\n",
    "plt.rcParams[\"figure.figsize\"] = (16,8)\n",
    "\n",
    "img_file = 'example_img.jpg'\n",
    "url = 'https://upload.wikimedia.org/wikipedia/commons/0/09/FluorescentCells.jpg'\n",
    "\n",
    "with open(img_file, 'wb') as file_obj:\n",
    "    file_obj.write(urlopen(url).read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The saved image data can the be read, creating a 3D NumPy array, and displayed using Matplotlib:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixmap = plt.imread(img_file)     # Read image data as array\n",
    "\n",
    "plt.imshow(pixmap)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The the array of pixel data `pixmap` is an array of unsigned 8-bit integers (in this instance), in the range 0 to 255. The size of the array is `width` x `height` for each colour channel, which forms the last axis. Note that for greyscale images there would be no channel axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "height, width, channels = pixmap.shape\n",
    "\n",
    "print(pixmap.dtype)\n",
    "print(pixmap.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For many operations is generally easier to work with data stored as floating-point numbers in the interval \\[0.0,1.0\\], rather than integers: this avoids many rounding and scale issues. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixmap = pixmap.astype(float)/255.0  # Convert from 0 .. 255 to 0.0 .. 1.0\n",
    "\n",
    "print(pixmap.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The colour channels may be separately accessed by taking the array slice and selecting a index for the last axis. Here, the histogram of intensity values is presented, after converting the 3D pixmap into a 1D array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red = pixmap[:,:,0]   # Color channels are last axis\n",
    "grn = pixmap[:,:,1]\n",
    "blu = pixmap[:,:,2]\n",
    "\n",
    "grn_values = grn.flatten()\n",
    "plt.hist(grn_values, bins=256, range=(0.0,1.0), log=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The individual colour channels may be displayed separately, and their intensities can be encoded as a color map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(grn, cmap='hot')  # Green channel with a recoloured intensity scheme \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Colours can be transformed by matrix multiplication: here the orginal, separate channels are mixed to make new colours."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = np.array([[0.5,0.0,0.5],    # Red maps to magenta : half red, half blue\n",
    "                [0.5,0.5,0.0],    # Green maps to yellow : half red : half green\n",
    "                [0.0,0.0,1.0]])   # Blue is unchanged\n",
    "\n",
    "pixmap2 = np.dot(pixmap, mat)      # Apply color trsnform matrix\n",
    "pixmap2 /= pixmap2.max(axis=(0,1)) # Normalize channels separately\n",
    "\n",
    "plt.imshow(pixmap2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SciPy module `ndimage` has various functions thjat are commonly used to manipulate images, and which would otherwise be tricky to achieve with just NumPy array manipulations. For example images can be rotated and mapped back to a rectangular array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import ndimage\n",
    "\n",
    "rot = ndimage.rotate(pixmap, 45)\n",
    "rot = np.clip(rot, 0.0, 1.0)    # Interpolation can give values just outside [0,1]\n",
    "\n",
    "plt.imshow(rot) \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=purple>Exercise 10: Greyscale</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>Convert the cell image used in above exmples to greyscale by finding the average of the red, geen and blue values for each pixel.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code\n",
    "\n",
    "pixmap = plt.imread(img_file)     # Read image data as array\n",
    "grey = pixmap.mean(axis=2)\n",
    "grey -= grey.min()\n",
    "grey *= 255/grey.max()\n",
    "grey = grey.astype(int)\n",
    "pixmap2 = np.dstack([grey, grey, grey])\n",
    "\n",
    "plt.imshow(pixmap2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ther are also various useful functions to apply convolutional filters, e.g. for blurring or edge detection, and other anlytic functions, like for computing the distance transform. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (sp1, sp2, sp3, sp4) = plt.subplots(1, 4)\n",
    "\n",
    "gauss = ndimage.gaussian_filter(grn, 5.0)  # Sigma is 5.0\n",
    "sp1.imshow(gauss)\n",
    "sp1.set_title('Gauss blur')\n",
    "\n",
    "sobel = ndimage.sobel(grn)\n",
    "sp2.imshow(sobel, cmap='Greys')\n",
    "sp2.set_title('Sobel edge')\n",
    "\n",
    "mask = (grn > grn.mean()).astype(int) # Makes an array of 0 or 1\n",
    "sp3.imshow(mask, cmap='Greys')\n",
    "sp3.set_title('Binary mask')\n",
    "\n",
    "dtrans = ndimage.distance_transform_edt(mask)\n",
    "sp4.imshow(dtrans, cmap='hot')\n",
    "sp4.set_title('Distance transform')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font color=purple>Exercise 11: High-pass filter</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=purple>Create a high-pass Gaussian filter for the greyscale image created in the previous exercise. Do this by creating a Gaussian blurred version of the image (try $\\sigma=3$) and then subtract this from the original greyscale array. Display the original and high-pass images side-by-side for comparison. </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise code\n",
    "fig, (sp1, sp2) = plt.subplots(1, 2)\n",
    "\n",
    "g1 = ndimage.gaussian_filter(pixmap2, 3.0)\n",
    "d1 = (pixmap2-g1).astype(float)\n",
    "d1 -= d1.min()\n",
    "d1 /= d1.max()\n",
    "\n",
    "sp1.imshow(pixmap2)\n",
    "sp2.imshow(d1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "### Example : PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below function illustrates the use of the NumPy module and performs a principle component analysis: treating input data as vectors it finds the orthogonal directions in the data of maximal variance (the Eigenvectors of the covariance matrix). This is often used on high-dimensionality data to create simpler representations that still preserve the most important features. The function takes two arguments, the input data, which assumed to be equivalent to a list of vectors, and the number of principle components to extract. There is a small complication in this function as `linalg.eig()` outputs a matrix (`p_comp_mat` below) where each Eigenvector is a column, rather than a row. This orientation is useful for applying the matrix as a transformation, as we demonstrate below. Though, in the code it means the matrix is sorted and selected on its last axis (using `[:,:n]` etc.)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "def pca(data, n_comp=2):\n",
    "\n",
    "    data = np.array(data)               # Convert input to array\n",
    "    mean = data.mean(axis=0)            # Mean vector\n",
    "    centred_data = (data - mean).T      # Centre all vectors and transpose\n",
    "\n",
    "    covar = np.cov(centred_data)        # Get covariance matrix\n",
    "    evals, evecs = np.linalg.eig(covar) # Get Eigenvalues and Eigenvectors\n",
    "\n",
    "    idx = evals.argsort()[::-1]         # E. value indices by decreasing size\n",
    "    evecs = evecs[:,idx]                # Sort Eigenvecs according to Eigenvals\n",
    "  \n",
    "    p_comp_mat = evecs[:,:n_comp]       # Select required principle components\n",
    "\n",
    "    return p_comp_mat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function is tested using some random 3D data. Here the random module from NumPy is used to create three clusters of vector points. Initially each has the same mean (0.0) and standard deviation (0.5), but the last two clusters are transposed by adding an offset vector. The clusters are concatenated together (along the long axis) to create the final test dataset with three ‘blobs’."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = (100, 3)   # 100 points times 3 dimensions\n",
    "d1 = np.random.normal(0.0, 0.5, size) \n",
    "d2 = np.random.normal(0.0, 0.8, size) + np.array([4.0, 1.0, 2.0])\n",
    "d3 = np.random.normal(0.0, 0.7, size) + np.array([2.0, 0.0, -1.0])\n",
    "test_data = np.concatenate([d1, d2, d3], axis=0)\n",
    "x, y, z = test_data.T\n",
    "\n",
    "fig, (sp1, sp2, sp3) = plt.subplots(1,3)\n",
    "sp1.scatter(x, y, alpha=0.3)\n",
    "sp2.scatter(y, z, alpha=0.3)\n",
    "sp3.scatter(z, x, alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The principle components are computed and given back as an array, albeit in transposed (data_dim, pc) form. These can be plotted on the original data axes to show the direction of the principle components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcomps = pca(test_data, 2)   # Run PCA\n",
    "plt.scatter(x, y, alpha=0.5, color='#FFB000')\n",
    "\n",
    "for pc in (0,1):\n",
    "    x_val = (0, 3 * pcomps[0,pc])  # X dim, pc\n",
    "    y_val = (0, 3 * pcomps[1,pc])\n",
    "    plt.plot(x_val, y_val)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The principle component matrix can be used to transform the test data. Here the first two principle components are used as new X and Y axis directions, illustrating that the transformation gives a better separated 2D view of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transformed = np.dot(test_data, pcomps)\n",
    "\n",
    "xt, yt = transformed.T\n",
    "plt.scatter(xt, yt, alpha=0.5, color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}