day_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python Programming for Scientists - Day 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "In the third day we will learn the two most important libraries for data analysis in python:\n",
    "\n",
    "* [numpy](https://numpy.org/doc/stable/user/index.html) (for efficient numeric arrays)\n",
    "* [matplotlib](https://matplotlib.org/stable/plot_types/index.html) (for plotting)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [1] Numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far our only way to store and manipulate sets of numbers is with lists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [1,3,5,6,2]\n",
    "y = [0,2,9,1,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to sum these two lists, element-wise, to produce a new array $z = x + y$, in Python we would need to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = []\n",
    "for i in range(len(x)):\n",
    "    z.append(x[i]+y[i])\n",
    "    \n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Compute $z=x+y$ without a loop (using list comprehension instead)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works, but in general lists in Python are:\n",
    "\n",
    "* flexible - e.g. easy to make smaller or larger in size\n",
    "* heterogeneous - may contain mixes of different types\n",
    "\n",
    "However, flexibility comes at the cost of performance (slow) and storage (lots of memory), and lists are not good for numerical calculations.\n",
    "\n",
    "This is where **numpy** comes in: this is a Python library that defines a powerful **n-dimensional array** object ([ndarray](https://numpy.org/doc/stable/reference/arrays.ndarray.html)) which is: high performance (fast) and lean (no excess memory).\n",
    "\n",
    "> Note: numpy arrays can contain only numbers of exactly the same type, no mixing and matching!\n",
    ">\n",
    "> Note: no strings! Stick to plain Python.\n",
    "\n",
    "To start using numpy, we import the module. In all code you will ever see, it is renamed upon import to the shorthand \"np\", just for convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating numpy arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "One way to create an ndarray is from a Python list, using the `array` function to convert an input list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([11, 22, 33, 44, 55])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Numpy arrays have several attributes that give useful information about the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.ndim  # number of dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.size # total number of elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a.shape  # shape of the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a.dtype  # numerical type"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When you create a numpy array, you should always decide (and explicitly specify) its type!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([11, 22, 33, 44, 55], dtype='int32')\n",
    "a.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most common and useful data types are:\n",
    "\n",
    "* **int16** - integer numbers, from -32768 to 32767.\n",
    "* **int32** - integer numbers, from -2147483648 to 2147483647 (~two billion).\n",
    "* **int64** - integer numbers, from -9223372036854775808 to 9223372036854775807 (very big).\n",
    "* **float32** - \"single precision\" float number, with ~7 digits of precision (exponent can be anything up to ~100).\n",
    "* **float64** - \"double precision\" float number, with ~15 digits of precision (exponent can be anything up to ~1000).\n",
    "\n",
    "Only numbers within the minimum and maximum ranges can be stored by each type! Otherwise [overflow](https://en.wikipedia.org/wiki/Arithmetic_overflow) occurs, leading to big problems.\n",
    "\n",
    "> For example: An array of int64 takes 64 bits (8 bytes) per entry, so the total amount of memory (and space on disk) is -twice as much- as if you had used int32 (4 bytes per entry) instead.\n",
    "\n",
    "Also available: `int8` (integers up to 128), `complex64` (two single precision floats).\n",
    "\n",
    "> There are also \"unsigned\" versions of all integers, which allow only positive numbers and so double the maximum storable value, `uint16`, `uint32`, and `uint64`.\n",
    ">\n",
    "> Only use if you have a good reason to! Negative numbers can be silently converted to unexpected positive numbers, leading to hard to fix bugs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating numpy arrays (cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "There are several other ways, and helper functions, to create numpy arrays.\n",
    "\n",
    "The `arange()` function works just like the built-in Python `range()` function, but it can also produce floating point ranges:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arange(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "np.arange(3, 12, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "np.arange(1.2, 4.4, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The `linspace()` function creates linearly (evenly) spaced values between a minimum and maximum, with a given total number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linspace(1.0, 2.0, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The `logspace()` function is similar, but creates logarithmically spaced values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.logspace(1.0, 3.0, 11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The most general and best way to create an empty numpy array is with `zeros()`, which initializes all elements of the array to zero upon creation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.zeros(10, dtype='float32')\n",
    "y = np.zeros(1000, dtype='int64')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> You can also use `empty()`, which does not set the elements to zero (they can random random values, from whatever was left over in memory before!). A bit faster. Not worth the possible bugs, always use `zeros()` instead.\n",
    ">\n",
    "> You can also use `ones()`, to initialize all elements to one at creation. This is the same as `zeros()+1`.\n",
    "\n",
    "Every time you create a numpy array, there are **two critical choices: dtype and size**. You should assume that these **can never change** for a given array after it has been made."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Create a numpy array which contains the value 2 repeated 10 times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing and Slicing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "> \"An array can be indexed by a tuple of nonnegative integers, by booleans, by another array, or by integers.\"\n",
    "\n",
    "Similarly to lists, individual items in arrays can be accessed with an integer **index**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([9,8,7,6,5,1,1,1,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Arrays can be **sliced** by specifiying the start and end of the slice (where the last element is exclusive):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[0:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "When slicing, we can optionally specify a step size (or \"stride\"):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x[0:6:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "As for lists, the `start`, `end`, and `step` are all optional, and default to ``0``, ``len(array)``, and ``1`` respectively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[::2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: the common trick to reverse the order of an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of a single index, or a slice syntax, if you have a list (or tuple, or array) of integer **indices**, you can retrieve the corresponding elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inds = [0,1,6,7]\n",
    "x[inds]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Array subsets based on a condition (i.e. searching)\n",
    "\n",
    "One of the powerful features of numpy is that you can select elements from an array that fulfill a (logical) condition.\n",
    "\n",
    "There are two common ways to do this.\n",
    "\n",
    "1. First, you can use a boolean \"mask\": an array of bools, either True or False, with the same size as the original array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = (x > 1)\n",
    "mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You don't need to save this intermediate step, you can use it directly and more succintly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[x>1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Second, you can obtain and use an explicit list of indices. The **`np.where()`** function is essential, you will use it a lot!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.where(x > 1)\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[w]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In these cases, you can combine multiple conditions with the **logical operators** `&` (and), `|` (or):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.where((x>1) & (x<4))\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[(x>6) | (x<2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Create a numpy array from the list `[3,-np.pi,5,89,22,938,15,-43]`. Print the indices of all elements which are either negative, or even."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Special constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a few important constants and special \"values\":\n",
    "\n",
    "* `np.pi` - this is the usual $\\pi \\simeq 3.14$\n",
    "* `np.e` - this is the usual $e \\simeq 2.718$\n",
    "* `np.inf` - this is an official \"value\" in the standard for floating point numbers, which represents $+\\infty$.\n",
    "* `-np.inf` - represents $-\\infty$.\n",
    "* `np.nan` - this is also an official \"value\" for floating point numbers, it represents \"Not A Number\" (NaN), and should be used to mark e.g. missing data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-dimensional arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make 2D, 3D, and N-D arrays as well. The first argument `np.zeros()` is the desired shape, i.e. number of elements per dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2d = np.zeros((3,3), dtype='int32')\n",
    "\n",
    "print(x2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x3d = np.zeros((2,4,3), dtype='int32')\n",
    "\n",
    "print(x3d.size, x3d.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indexing and slicing multi-dimensional arrays works as in the 1D case. Here we assign new values, instead of extracting existing values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2d[0,0] = 42     # first index applies to first dimension, second index to second dimension\n",
    "x3d[0,3,1] = 11   # likewise for three dimensions\n",
    "x3d[1,1,2] = 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gt_zero = (x2d>0) # boolean mask array has the same dimensionality and shape as x2d\n",
    "gt_zero.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you noticed above, the return of `np.where()` was actually always a tuple with just one element, with the indices inside this first element. This is because we were searching a one-dimensional array. In general, the return is a tuple (i.e. list) with one entry per dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.where(x3d > 1) # two elements satisfy this condition, which are returned in a 3-tuple, one entry per dimension\n",
    "\n",
    "type(w), len(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can look a bit confusing, but this is exactly a format which can be used to index a multi-dimensional array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x3d[w]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How could we loop over a multi-dimensional array?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for i in range(x2d.shape[0]):\n",
    "    for j in range(x2d.shape[1]):\n",
    "        value = x2d[i,j]\n",
    "        print(f'{i = }, {j = }, {value = }')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Create a 2x2 numpy array with four unique entries (e.g. the numbers `1, 2, 3, 4`). Then write a function `transpose()` which takes such an array as input, and returns its transpose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "You can also slice multi-dimensional arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2d[0,:] # all elements along the second dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2d[0:2,0:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When fewer indices are provided than the number of axes, the missing indices are considered complete slices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2d[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expression within brackets `\"0\"` is treated as a `0` followed by as many instances of `:` as needed to represent the remaining axes. NumPy also allows you to write this using **multi-dimensional dot notation**: `x[i, ...]`.\n",
    "\n",
    "The dots (`...`) represent as many colons as needed to produce a complete indexing tuple.\n",
    "\n",
    "> For example, if `x` is an array with 5 axes, then\n",
    "> * `x[1, 2, ...]` is equivalent to `x[1, 2, :, :, :]`,\n",
    "> * `x[..., 3]` to `x[:, :, :, :, 3]` and\n",
    "> * `x[4, ..., 5, :]` to `x[4, :, :, 5, :]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x3d[1,...]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \"Row major\" vs \"column major\" (i,j vs j,i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As soon as an array is multidimensional, it becomes important to understand how it is actually represented on the computer.\n",
    "\n",
    "This is sometimes phrased as the question: \"is indexing row-major, or column-major\", in the sense of: does $x[i,j]$ take the element of $x$ corresponding to row $i$ and column $j$, or row $j$ and column $i$?\n",
    "\n",
    "* If $x$ represents a 2D \"matrix\" (or 3D \"tensor\"), then it can be useful to associate the indexing $x[i,j]$ with a mathematical notation like $x_{i,j}$ for matrices.\n",
    "\n",
    "* If $x$ represents a 2D \"image\", it can also be useful to associate the indexing $x[i,j]$ with a given pixel, ordered e.g. from the top left or bottom right corner.\n",
    "\n",
    "However, an important concept is that **there is no such thing as a multi-dimensional array in computer memory**! Memory is simply 1D, and any other structure is simply an -interpretation- of numbers arranged in 1D.\n",
    "\n",
    "We can see this by \"flattening\" a 2D array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2d.flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2d.flatten().shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, we can \"reshape\" an array of a given shape into any other shape, so long as the size (number of elements) does not change.\n",
    "\n",
    "> This reshaping does not change, in any way, the numbers stored in memory, nor their order. It only changes how they are interpreted by numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([9,8,7,6,5,1,1,1])\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = x.reshape((4,2))\n",
    "\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = x.reshape((2,4))\n",
    "\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most useful way of thinking about multi-dimensional array order is to never think in terms of rows, columns, or similar analogies.\n",
    "\n",
    "Instead, you should always think about which index (dimension) changes the \"fastest\" as you move through the array in memory.\n",
    "\n",
    "> NumPy uses C-order indexing. That means that **the last index is the most rapidly changing** (memory location), unlike Fortran or IDL, where the first index represents the most rapidly changing location in memory. This difference represents a great potential for confusion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Arithmetic (and broadcasting)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arithmetic operators on arrays apply **elementwise**. A new array is created and filled with the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([20, 30, 40, 50])\n",
    "b = np.arange(4) # [0,1,2,3]\n",
    "\n",
    "c = a - b\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "10 * np.sin(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: When operating with arrays of different types (e.g. int16 and int32), the type of the resulting array corresponds to the more general or precise one (\"upcasting\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a * b # element by element!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some operations, such as `+=` and `*=`, act in place to modify an existing array rather than create a new one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a += 1\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happened there? `a` is a one-dimensional array of size 4, but we added a single number (scalar) to it, i.e. a 0-dimensional array of size 1.\n",
    "\n",
    "Numpy has used **broadcasting** to make the computation. The term broadcasting describes how NumPy treats arrays with different shapes during arithmetic operations. Subject to certain constraints, the smaller array is “broadcast” across the larger array so that they have compatible shapes. This is quite powerful:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython import display\n",
    "display.Image(\"images/day3_broadcasting_1.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: The 'stretching' is only conceptual. Numpy does not actually create these intermediate arrays, saving memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([[ 0.0,  0.0,  0.0],\n",
    "              [10.0, 10.0, 10.0],\n",
    "              [20.0, 20.0, 20.0],\n",
    "              [30.0, 30.0, 30.0]])\n",
    "b = np.array([1,2,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = a + b\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The array `b` has been \"broadcast\" (extended) from a shape of (3,) to the shape of (4,3), so that it matches `a`. Then the two are summed element-wise.\n",
    "\n",
    "> To broadcast, when operating on two arrays, NumPy compares their *shapes*, element-wise.\n",
    ">\n",
    "> It starts with the trailing (i.e. rightmost) dimensions and works its way left. Two dimensions are compatible when\n",
    "> * they are equal, or\n",
    "> * one of them is 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython import display\n",
    "display.Image(\"images/day3_broadcasting_2.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> A one dimensional array added to a two dimensional array results in broadcasting if number of 1-d array elements matches the number of 2-d array elements along the last dimension."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Re-define `a` to have a shape of `(3,4)` instead of `(4,3)`, and try the same operation `c = a+b`. What happens?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython import display\n",
    "display.Image(\"images/day3_broadcasting_3.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> When the trailing dimensions of the arrays are unequal, broadcasting fails because it is impossible to align the values in the \"rows\" (first dimension) of the 1st array with the elements of the 2nd arrays for element-by-element addition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Create an array `x` with integers `0, 2, 4, ..., 10`, then compute a new array `dx` values where `dx[i] = x[i+1] - x[i]`. Do this without loops! Hint: Make `dx.size == x.size - 1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Built-in mathematical functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We saw above the use of `np.sin(x)` to compute the sin of all the elements of a numpy array `x`. These functions operate elementwise on an array, producing an array of the same shape as output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = np.arange(5)\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.cos(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.exp(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.std(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = np.linspace(0,3,12)\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.floor(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.ceil(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.round(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.isfinite(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consult the [mathematical functions documentation](https://numpy.org/devdocs/reference/routines.math.html) for a complete list."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many of these can be called \"reductions\", as they do not return a new array of the same size, but instead reduce it to e.g. a single number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many of these functions can be applied as methods to an ndarray. There is no functional difference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.cumsum(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.diff(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.clip(B,0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.max(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case of multi-dimensional arrays, you can perform a reduction along a specified axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([[ 0.0,  0.0,  0.0],\n",
    "              [10.0, 10.0, 10.0],\n",
    "              [20.0, 20.0, 20.0],\n",
    "              [30.0, 30.0, 30.0]])\n",
    "a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(a, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(a, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Compute the maximum of `a` along each axis. At what index, along each axis, does the maximum occur? (Hint: use [np.argmax()](https://numpy.org/doc/stable/reference/generated/numpy.argmax.html#numpy.argmax))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Create a new one-dimensional array with the following values: `[1.0, 4.2, np.nan, 7.9]`. The third value, \"not a number\", could represent e.g. missing or corrupted data. Compute the sum: what happens?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most numpy functions have \"nan-aware\" equivalents, e.g. `nansum()` instead of `sum()`, which ignore (skip) all NaN values on purpose.\n",
    "\n",
    "> Use when appropriate. Use only when needed, not everywhere. Seeing a NaN in your output can help you trace back to where such \"bad values\" started entering your computation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Random number generation is an important part of many algorithms and data analyses. In modern numpy, the way to do this is to first create a \"generator\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then sample from different random distributions, requesting $N=10$ numbers in each case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng.standard_normal(10) # Gaussian (mean = 0, stddev = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng.normal(loc=0.0, scale=2.0, size=10) # Gaussian (mean=0, stddev = 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng.random(10) # uniformly distributed floats in the half-open interval [0.0, 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng.integers(low=5, high=50, size=10) # random integers from low (inclusive) to high (exclusive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([3,5,1,2,4])\n",
    "rng.shuffle(x) # randomly permute\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: in addition to uniform and normal (Gaussian) distributions, there are many others you can draw from: Beta, binomial, Dirichlet, exponential, F, gamma, logistic, lognormal, poisson, power, rayleigh, ...\n",
    "\n",
    "It is often useful to control -which- random numbers are made: you may want random numbers, but the same random numbers every time you run your script. Otherwise, debugging a program with randomness can be a nightmare.\n",
    "\n",
    "This is solved by setting the **random seed** when initializing the generator. The seed is just a number, any number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(424242)\n",
    "rng.random(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Combining numpy arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Several arrays can be stacked together along different axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([[9, 7],\n",
    "              [5, 2]])\n",
    "b = np.array([[1, 9],\n",
    "              [5, 1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.hstack((a,b))\n",
    "print(c)\n",
    "print(c.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.vstack((a,b))\n",
    "print(c)\n",
    "print(c.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also use [np.concatenate()](https://numpy.org/doc/stable/reference/generated/numpy.concatenate.html#numpy.concatenate) for more complex cases.\n",
    "\n",
    "> Note: in general, **avoid** combining numpy arrays together, or changing the size of existing numpy arrays, if at all possible! In practice, you are always creating an entirely new array, and then filling it with some existing values (slow). There is probably a better way to do what you're trying."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Copies vs views"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When operating and manipulating arrays, their **data is sometimes copied (into a new array), and sometimes not (only a 'view' is made)**. This is often a source of confusion, and can make big problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([9,8,7,6,5,1,1,1])\n",
    "y = x # no copy, y is just a "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `y` is exactly the same as `x` - these are literally two \"views\" into the same data in memory (two names for the same ndarray). **Changing a view changes the original array as well!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y[2] = 0\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y is x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple assignments make no copy of objects or their data. Also, because all arguments to functions in python are passed [by reference](https://stackoverflow.com/questions/373419/whats-the-difference-between-passing-by-reference-vs-passing-by-value), no copy of numpy arrays are made when you pass them into a function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(my_array):\n",
    "    print(my_array[0]) # my_array is a view as well...\n",
    "    my_array[0] = -1\n",
    "    \n",
    "f(x)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While a view points to the same data, it can interpret it differently:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = x.reshape((4,2))\n",
    "print(z.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.shape # unchanged"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Important: **slicing an array returns a view** of it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = x[2:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z[:] = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Create the same view as above `zz = x[2:5]`, but then assign `zz = 20`. What happens to `zz`, and what happens to `x`? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you explicitly want to avoid a view, and make sure that you create an independent copy of the data, so that it can be modified without changing the original array, you can **make a deep copy**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = x[2:5].copy()\n",
    "arr[:] = 30\n",
    "print(arr)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr2 = x.copy()\n",
    "arr2[2] = 33\n",
    "print(arr2)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following questions are to check your understanding of this subtle issue - no need to write code, just figure out what the output will be in each case. After, you can check if you were correct.\n",
    "\n",
    "What will `x` be after the following?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [1, 2.0, [6,7,8], 'hello']\n",
    "y = x[0]\n",
    "y = 4.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What will `x` be after the following?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [1, 2.0, [6,7,8], 'hello']\n",
    "y = x[2]\n",
    "y.append(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What will `x` be after the following?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [1, 2.0, [6,7,8], 'hello']\n",
    "y = x[2]\n",
    "y = [4, 5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What will `x` be after the following?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [1, 2, 3, 4]\n",
    "y = x[::2]\n",
    "y[0] = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What will `x` be after the following?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([1, 2, 3, 4])\n",
    "y = x[1:]\n",
    "y[0] = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What will `x` be after the following?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([1, 2, 3, 4])\n",
    "y = x[::2]\n",
    "y[0] = 9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sorting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common task, which is the basis for many algorithms, is sorting (re-ordering an array such that elements are in order from smallest to largest, or vice versa)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([9,8,7,6,5,1,1,1,2])\n",
    "y = np.sort(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second way to sort is not to return the actual sorted elements, but instead to return the -indices- which, when used to index the original array, return it in sorted order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sort_inds = np.argsort(x)\n",
    "print(sort_inds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[sort_inds]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> This is useful if you have two (or more) arrays, and you want to sort them all, based on the elements of one (i.e. re-arrange them all in the same way)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "The [data/day3_munich_temps_bad_data.txt](data/day3_munich_temps_bad_data.txt) text file is the same as before, with the temperature (C) in Munich every day for several years, but it now contains some bad data points, e.g. caused by malfunctions of the measurements, or problems in the data collection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, read in the file using [np.loadtxt()](https://numpy.org/doc/stable/reference/generated/numpy.loadtxt.html) -- get used to reading the numpy documentation.\n",
    "\n",
    "Then, identify bad values by looking at the minimum and maximum values of the array. Use numpy masking/where selection to get rid of the bad temperature values. What fraction of data points remain? What is the true min and max across the data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [2] Matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several plotting libraries you can use in Python: there is no single \"way to plot\" in Python. Modern options include:\n",
    "\n",
    "* [matplotlib](https://matplotlib.org/) - most general\n",
    "* [plotly](https://plotly.com/python/) - emphasis on interactive plots\n",
    "* [bokeh](https://bokeh.org/) - also an emphasis on interactive plots\n",
    "* [seaborn](https://seaborn.pydata.org/) - emphasis on stats\n",
    "\n",
    "That said, **matplotlib** is by far the most widely used, definitely the most suggested, and we will describe how to use this package here.\n",
    "\n",
    "As with numpy, it is convention to import the library as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: In the past, you had to enter `%matplotlib inline` in a cell of a Jupyter notebook to get plots to appear. This is no longer needed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A basic figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A matplotlib figure looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython import display\n",
    "display.Image(\"images/day3_matplotlib.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Constructing a basic figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots() # a single axis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Remember! This syntax is just taking the return from `plt.subplots()`, which has two elements, and \"naming\" them into two separate variables.\n",
    ">\n",
    "> Exactly the same could be achieved with:\n",
    ">\n",
    "> ```\n",
    "> my_plot = plt.subplots()\n",
    "> fig = my_plot[0]\n",
    "> ax = my_plot[1]\n",
    "\n",
    "Let's make some random data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(424242)\n",
    "x = rng.random(50)\n",
    "y = rng.random(50) * 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(x,y) # in the Jupyter notebook, we can add a ';' after the last plotting command to suppress the object name/memory address from being printed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most basic command, [plot(x,y)](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html), plots $x$ versus $y$ as either connected lines, or individual markers.\n",
    "\n",
    "Matplotlib works naturally with numpy array as inputs (as done here)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Use the documentation of `plot()` to figure out how to change from lines to individual markers. Plot the same data with green circles (`o`), and then blue squares (`s`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving figures to a file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are writing a paper, or want to share a plot with a colleague, changes are you will want to save a figure either as a PDF (always the best option) or an image (e.g. PNG).\n",
    "\n",
    "* **PDF is a vector format**: all text, lines, symbols, etc are represented by equations/formulas. Can \"zoom in\" as much as you like, things never get blurry. (File size increases with amount of content).\n",
    "* PNG and JPG, **images are raster formats**: all information is represented by a grid of colored pixels. Things get blurry as you zoom in. (File size is ~constant, depends only on `N_px_height * N_px_width`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(x,y)\n",
    "\n",
    "# after all plotting commands are done, save to a file:\n",
    "fig.savefig('my_fig.png')\n",
    "plt.close(fig) # good habit (free memory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Styling markers, linestyles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, y, marker='o', color='orange');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, y, marker='D', linestyle='dotted', color='orange'); # solid, dotted, dashed, dashdot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, y, marker='o', color='orange', linewidth=0); # or \"lw=0\" for short"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consult the matplotlib [linestyles list](https://matplotlib.org/3.5.1/gallery/lines_bars_and_markers/linestyles.html) and the [markers reference](https://matplotlib.org/3.5.1/gallery/lines_bars_and_markers/marker_reference.html) for all the options.\n",
    "\n",
    "> You can also create custom linestyles and markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython import display\n",
    "display.Image(\"images/day3_matplotlib_markers.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Load the same `data/day2_munich_temps.txt` file we have used before, and make a line plot of temperature versus time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple datasets and labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can add as many datasets (e.g. point sets, lines) to an axis. If you provide a `label` for each, then the `ax.legend()` function can automatically create a legend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, y, marker='o', color='orange', lw=0, label='cats')\n",
    "ax.plot(x/1.5, y*2, marker='s', color='green', lw=0, label='dogs')\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Data points with errorbars\n",
    "\n",
    "The [errorbar()](https://matplotlib.org/3.5.1/api/_as_gen/matplotlib.axes.Axes.errorbar.html#matplotlib.axes.Axes.errorbar) function is very similar to `plot()`, except it can also draw lines representing errors or uncertainties, in either or both coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some data for \"errors\"\n",
    "xerr = np.sin(x*np.pi)/10\n",
    "yerr = y/10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(x, y, yerr, fmt='o');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(x, y, yerr, xerr, fmt='s');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scatter plots and colorbars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second plotting command after `ax.plot()` is `ax.scatter()`. Nearly identical, and the two can almost be used interchangeably. One difference is that [scatter()](https://matplotlib.org/3.5.1/api/_as_gen/matplotlib.pyplot.scatter.html) gives an easy way to assign different colors (or markers, or sizes) to each point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's generate some more data\n",
    "color_vals = np.sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(x, y, c=color_vals);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you use color in a figure, it is almost always a good idea to include a colorbar. Otherwise the mapping from color to value is unknown, which is not scientifically useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "s = ax.scatter(x, y, c=color_vals)\n",
    "fig.colorbar(s);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: We need to pass an 'object' to `fig.colorbar()` which has the needed information about the color mapping which has already been applied.\n",
    ">\n",
    "> Instead of calling `ax.plot()` or `ax.scatter()` and discarding the return value, we can instead save it (`s`) for this purpose.\n",
    "\n",
    "Optionally, give the colorbar a label:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "s = ax.scatter(x, y, c=color_vals)\n",
    "cbar = fig.colorbar(s)\n",
    "cbar.set_label('my color value');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: We use the exact same strategy here, saving the colorbar object `cbar` so that we can use its methods to later change its properties.\n",
    "\n",
    "Every matplotlib command, and object, has a huge number of options, methods, and so on. Always skim through the documentation to see to change or do anything.\n",
    "\n",
    "For example, we can set the [colormap](https://matplotlib.org/stable/tutorials/colors/colormaps.html) that is used with `cmap`, and we can use specified bounds for the color mapping `vmin,vmax` (instead of the automatic defaults). The colorbar automatically adapts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "s = ax.scatter(x, y, c=color_vals, cmap='inferno', vmin=0.1, vmax=0.6)\n",
    "cbar = fig.colorbar(s)\n",
    "cbar.set_label('my color value');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's generate some more data to use for symbol sizes\n",
    "size_vals = y * 600"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "s = ax.scatter(x, y, c=color_vals, s=size_vals, cmap='jet')\n",
    "cbar = fig.colorbar(s);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When markers or lines start to overlap and get crowded, the simplest solution is to make them partially transparent with **alpha**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "s = ax.scatter(x, y, c=color_vals, s=size_vals, cmap='jet', alpha=0.7)\n",
    "cbar = fig.colorbar(s);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Make a scatter plot of `x` versus `y`, but use our \"size_vals\" array to also specify color `c`. What happens? Try to use our \"color_vals\" array to specify *size*. What happens? Can you fix it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Customizing axes (labels, ranges, ticks), titles, grids, etc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every aspect of a matplotlib figure can be customized. By far the most important are the axes labels (by default: blank) and axes limits (by default: automatic).\n",
    "\n",
    "**Never make a figure without specifying axis labels and, if appropriate, physical units.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "s = ax.scatter(x, y, c=color_vals, s=size_vals, cmap='jet', alpha=0.7)\n",
    "ax.set_xlabel('x coordinate [cm]')\n",
    "ax.set_ylabel('y coordinate [$\\mu$m]')\n",
    "\n",
    "cbar = fig.colorbar(s)\n",
    "cbar.set_label('Density [g cm$^{-3}$]');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.set_title('Many Overlapping Circles')\n",
    "ax.set_xlabel('x coordinate [cm]')\n",
    "ax.set_ylabel('y coordinate [$\\mu$m]')\n",
    "\n",
    "ax.set_xlim([0,1])\n",
    "ax.set_ylim([-1,2.5])\n",
    "\n",
    "ax.set_xticks(np.linspace(0.1,0.9,5))\n",
    "\n",
    "s = ax.scatter(x, y, c=color_vals, s=size_vals, cmap='jet', alpha=0.7)\n",
    "\n",
    "cbar = fig.colorbar(s)\n",
    "cbar.set_label('Density [g cm$^{-3}$]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In science we often have large ranges of data values, and it is common to use **log scaling** for axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.set_ylim([1e-2,1e1]) # minimum value must be >0\n",
    "ax.set_yscale('log') # y-axis ticks and labels re-written to be in log\n",
    "ax.set_ylabel('y value') # label omits \"log\"\n",
    "\n",
    "s = ax.scatter(x, y, c=color_vals, s=size_vals, cmap='jet', alpha=0.7);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Histograms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "The `ax.hist(x)` computes and plots a (one-dimensional) histogram of the quantity `x`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.hist(y)\n",
    "ax.set_xlabel('N')\n",
    "ax.set_ylabel('y');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Use the `bins` argument to `ax.hist()` to increase the number of bins. What happens? Use the `range` argument to specify a 2-tuple of the `min,max` for the binning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "The `np.histogram(x)` function computes a histogram of the quantity `x`. Use it to compute the histogram of the same `x` array, then plot the result using matplotlib.\n",
    "\n",
    "> Consult the [np.histogram() documentation](https://numpy.org/doc/stable/reference/generated/numpy.histogram.html) if needed.\n",
    ">\n",
    "> Careful with the return of `histogram()`. Notice how it is a tuple containing two arrays. What does each mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple axes in one figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first line of all our figures so far has been the same `plt.subplots()`, which is a convenience function which (i) creates a new Figure, and (ii) creates one or more Axes within that Figure. By default, this is just one, but we can change this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,2) # nrows,ncols\n",
    "\n",
    "axes[0].hist(y)\n",
    "axes[0].set_xlabel('y')\n",
    "axes[0].set_ylabel('N')\n",
    "\n",
    "axes[1].scatter(x, y, c=color_vals, s=size_vals, cmap='jet', alpha=0.7);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default **figsize** is a bit small, and is roughly the golden ratio for a single panel. But if we have more than one panel (axis), this will need to be changed.\n",
    "\n",
    "> Note: The sizes of all \"Artists\" in a Figure (e.g. lines, markers, axis labels, line widths) are relative to figsize. For example, using a larger figsize will cause all text to become smaller (if you put the plot into the same final space, e.g. A4 paper width)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(12,6))\n",
    "\n",
    "axes[0].hist(y)\n",
    "axes[0].set_xlabel('y')\n",
    "axes[0].set_ylabel('N')\n",
    "\n",
    "axes[1].scatter(x, y, c=color_vals, s=size_vals, cmap='jet', alpha=0.7);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Re-make the same plot as above, but make the `figsize` twice as large, then half as large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make one last example of a multi-panel figure:\n",
    "\n",
    "> Note:  \"ul\" = upper left, \"lr\" = lower right, and so on.\n",
    "\n",
    "> Note: the `sharex` and `sharey` keywords provide one way to \"share\" axes between panels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ((ul, ur), (ll, lr)) = plt.subplots(2, 2, figsize=(14,8), sharex=True)\n",
    "\n",
    "ul.hist(y)\n",
    "ul.set_ylabel('N')\n",
    "\n",
    "ur.scatter(x, y, c=color_vals, s=size_vals, cmap='jet', alpha=0.7)\n",
    "ur.set_ylabel('y coordinate [$\\mu$m]');\n",
    "\n",
    "ll.plot(x, y, lw=3, color='#00ff00')\n",
    "ll.set_ylabel('Velocity [km/s]')\n",
    "ll.set_xlabel('x [cm]')\n",
    "\n",
    "lr.hist(x)\n",
    "lr.set_ylabel('N')\n",
    "lr.set_xlabel('x [$10^3$ km]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D figures (x, y, and a third value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often we want to have a third dimension in a plot represented by a colors in a colormap. In addition to `scatter()`, matplotlib has a number of plot types that are especially made for this.\n",
    "\n",
    "First, we will generate some data to plot. The [np.meshgrid()](https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html) function is a useful helper, which can \"replicate\" 1D arrays into 2D by repeating them along the second axis.\n",
    "\n",
    "This is particularly useful to define a set of $N \\times N$ points in 2D space, where we will compute a function $f(x,y)$ on each of those points, and visualize the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define coordinate ranges along both the x,y axes\n",
    "x_1d = np.linspace(-3, 3, 32)\n",
    "y_1d = np.linspace(-3, 3, 32)\n",
    "\n",
    "# replicate these x,y coordinates into two 2D arrays\n",
    "X, Y = np.meshgrid(x_1d, y_1d)\n",
    "\n",
    "# function value:\n",
    "Z = (1 - X/2 + X**5 + Y**3) * np.exp(-X**2 - Y**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_1d.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_1d.min(), x_1d.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.min(), X.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the range of actual \"values\" across this range of $x$ and $y$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z.min(), Z.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create a Figure with four subplots, and plot the same data with four different techniques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(10,6))\n",
    "\n",
    "pc = axes[0, 0].pcolormesh(X, Y, Z, vmin=-1, vmax=1, shading='auto', cmap='RdBu_r')\n",
    "fig.colorbar(pc, ax=axes[0, 0])\n",
    "axes[0, 0].set_title('pcolormesh()')\n",
    "\n",
    "co = axes[0, 1].contourf(X, Y, Z, levels=np.linspace(-1.25, 1.25, 11))\n",
    "fig.colorbar(co, ax=axes[0, 1])\n",
    "axes[0, 1].set_title('contourf()')\n",
    "\n",
    "pc = axes[1, 0].imshow(Z, cmap='plasma')\n",
    "fig.colorbar(pc, ax=axes[1, 0], extend='both')\n",
    "axes[1, 0].set_title('imshow()')\n",
    "\n",
    "pc = axes[1, 1].scatter(X.flatten(), Y.flatten(), c=Z.flatten(), s=10.0, marker='o', cmap='inferno')\n",
    "fig.colorbar(pc, ax=axes[1, 1], extend='both')\n",
    "axes[1, 1].set_title('scatter()');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Make a single panel contour plot of the same data, using `contour()` instead of `contourf()` (the \"f\" meaning \"filled\"). Try changing the levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D numpy arrays as images\n",
    "\n",
    "To visualize two-dimensional data, or plot a two dimensional function, we can \"plot\" the 2D numpy array as an image. **Careful** with the indexing conventions, and with the axes limits!\n",
    "\n",
    "Let's generate a set of `x` coordinates (from `-1` to `+1`) and `y` coordinates (from `0.0` to `5.0`), then use [np.meshgrid()](https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html) to expand these \"coordinate vectors\" into two 2D arrays `xx` and `yy` containing the positions, along the x and y directions separately, for every element (i.e. \"location\" in x,y space)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1,1,40)\n",
    "y = np.linspace(0,5,20)\n",
    "\n",
    "xx, yy = np.meshgrid(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, imagine we want to compute the function $f(x,y) = \\cos(8x) * \\exp(y/2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ff = np.cos(xx*8) * np.exp(yy/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.imshow(ff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two major problems with this figure:\n",
    "\n",
    "3. The axis extents are counting the number of \"pixels\", instead of representing the actual values of the $x$ and $y$ coordinates. The shape of the image also reflects the number of pixels (array shape). Fix with `extent`.\n",
    "4. Our function $f(x,y)$ is sinusoidal in the x-direction, which is good, but decays in the y-direction as we go upwards from the bottom: not good!\n",
    "   Mathematically, our function should increase in the positive y-direction. The current visualization, then, is flipped with respect to a normal coordinate system in a mathematical figure. Fix with `origin`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "im = ax.imshow(ff, origin='lower', extent=[x.min(),x.max(),y.min(),y.max()])\n",
    "\n",
    "cb = fig.colorbar(im)\n",
    "\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "cb.set_label('f(x,y)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note! On my personal computer (not the KIP Lab), the `origin` fix was not actually necessary. The default behavior of different matplotlib versions can change. Better be safe than sorry, always specify options when in doubt!\n",
    "\n",
    "The key is to match the output of `np.meshgrid()` to the input of `imshow()`:\n",
    "\n",
    "* With the default meshgrid ordering, you should use `origin='upper'` for imshow.\n",
    "* You can instead specify the argument `indexing='ij'` to meshgrid, and then use `origin='lower'` with imshow.\n",
    "\n",
    "Regardless, always carefully **check the consistency** between such a 2D visualization and the intended orientation/indexing of the underlying 2D array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vector fields in 2D: streamlines and quiver plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As opposed to a \"scalar function\" in 2D such as $f(x,y)$ with one value as a function of position, a \"vector function\" in 2D $\\vec{f}(x,y)$ has two values for each position. A common example is a velocity vector field, where $\\vec{f} = v_x \\hat{x} + v_y \\hat{y}$ has a component in both the $x$ and $y$ directions. Two common ways to visualize such functions are:\n",
    "\n",
    "* [streamlines](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.streamplot.html)\n",
    "* [quiver plots](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.quiver.html)\n",
    "\n",
    "First, let's define some coordinates, and a function $\\vec{f} = -y \\hat{x} + x \\hat{y}$ which looks a bit like a vortex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define 1d coordinates, 2d coordinates array\n",
    "x = np.arange(-10,10,1)\n",
    "y = np.arange(-10,10,1)\n",
    "\n",
    "xx,yy = np.meshgrid(x,y)\n",
    "\n",
    "# define a 2d vector function component-wise\n",
    "f_x = -yy\n",
    "f_y = xx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "plt.streamplot(x,y,f_x,f_y);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "plt.quiver(x,y,f_x,f_y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Create a streamline plot of the same data, but add a `density=2` parameter. Try `density=0.5`. Make the color of the lines change as a function of the x-coordinate (look at the documentation linked above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Use numpy to generate 10000 random values following a Gaussian (\"normal\") distribution, and plot a histogram. Choose a mean $\\mu$ and standard deviation $\\sigma$ which are interesting. Choose a number of bins so that you can visually recognize the function. Overplot the analytic Gaussian function on top of it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid #bbb; margin: 30px 0\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 3 Practice Problem - Monte Carlo Error Propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "You have likely encountered the concept of \"propagation of uncertainty\" before: given a mathematical expression for a quantity $f(x,y,z,...)$ which depends on one or more variables $x, y, z, ...$, then if these quantities represent measurements with uncertainties ($x \\pm \\sigma_x$, and so on), we can compute the uncertainty $\\sigma_f$ on the final value $f$.\n",
    "\n",
    "For example, let us consider the following equation:\n",
    "\n",
    "$$F = \\frac{G~M_1~M_2}{r^2}$$\n",
    "\n",
    "which gives the gravitational force between two masses $M_1$ and $M_2$ separated by a distance $r$.\n",
    "\n",
    "Let us now imagine that we have two masses: $M_1=40\\times10^3\\pm0.05\\times10^3\\rm{kg}$ and $M_2=30\\times10^3\\pm0.1\\times10^3\\rm{kg}$ at distance $r=3.8\\pm0.02~\\rm{m}$, where the uncertaintes are the standard deviations of Gaussian distributions which could be e.g. measurement errors. We also know the value of the gravitational constant $G = 6.67384\\times10^{-11}~\\rm{m}^3~\\rm{kg}^{-1}~\\rm{s}^{-2}$ (assume no uncertainty).\n",
    "\n",
    "The [standard error propagation rules](http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Example_formulas) give us that, for the equation above,\n",
    "\n",
    "$$\\sigma_F = F \\left[ \\left(\\frac{\\sigma_{M_1}}{M_1}\\right)^2 + \\left(\\frac{\\sigma_{M_1}}{M_1}\\right)^2 + 4 \\left(\\frac{\\sigma_{R}}{R}\\right)^2 \\right]^{1/2}$$\n",
    "\n",
    "## Task A\n",
    "\n",
    "Use the analytic equations to determine the resulting force, and the uncertainty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task B\n",
    "\n",
    "Now, we can try using a **Monte-Carlo** technique instead. The idea behind Monte-Carlo techniques is to generate many possible solutions using random numbers and using these to look at the overall results. In the above case, you can propagate uncertainties with a Monte-Carlo method by doing the following:\n",
    "\n",
    "* randomly sample values of $M_1$, $M_2$, and $r$, `N=1000000` times, using the means and standard deviations given above\n",
    "\n",
    "* compute the gravitational force for each set of values\n",
    "\n",
    "You should do this with numpy arrays, and without any loops.\n",
    "\n",
    "Make a plot of the normalized histogram of these values of the force, and then overplot a Gaussian function with the mean and standard deviation derived with the standard error propagation rules. Make sure the plot allows you to accurately compare the two results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task C\n",
    "\n",
    "Make the errors 100 times larger, and repeat the comparison of analytic and Monte Carlo uncertainty estimates. Why do you think the differ? Which is more accurate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid #bbb; margin: 30px 0\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 3 Challenge Problem - Mandelbrot Fractal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [Mandelbrot set](https://en.wikipedia.org/wiki/Mandelbrot_set) is often used as an introduction to the idea of fractals.\n",
    "\n",
    "Its formal definition is the set of all complex numbers $c$ for which the function $f(z) = c + z^2$ does not diverge to infinity, when starting from $z=0$. (This \"iterative solution\" can be computed using a recursive function, or just with a `for` loop for a given number of iterations).\n",
    "\n",
    "In practice, \"diverge to infinity\" can be thought of how large $z$ (its magnitude i.e. absolute value) has become after a number of iterations $N$.\n",
    "\n",
    "## Task A\n",
    "\n",
    "Write a function `mandel()` which takes as input `z`, `c`, and `N`, and returns the final value of $|f(z)|$.\n",
    "\n",
    "Test it out with a few different values of $c$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task B\n",
    "\n",
    "We will visualize the \"Mandelbrot fractal\" by computing the value of $f(z)$ for a grid of points in the complex plane. While doing so, we will record the \"time\" (iteration number) when $z$ diverges (becomes large). We will then treat this 2D array of \"divergence times\" values as an image, and plotting it with `imshow()`.\n",
    "\n",
    "Write a function `mandelbrot()` which accepts four input parameters: `num_x`, `num_y`, `max_val`, and `N_iter`.\n",
    "\n",
    "1. Generate `x` coordinates in the range `[-2.5,1.5]` and `y` coordinates in the range `[-1.5,1.5]`, according to the requested number.\n",
    "\n",
    "2. Use the same `meshgrid()` approach as above to create 2D arrays `xx` and `yy`, of the `x` and `y` coordinate values.\n",
    "\n",
    "3. Compute `c = x + iy`, which are all the values of $c$ we will evaluate (a 2D array).\n",
    "\n",
    "4. Create a 2D divergence time array to hold our final answers. \n",
    "\n",
    "5. Iterate `N_iter` times and for each, calculate the new value of $z$. For any values whose absolute value are greater than `max_val`, save the current iteration number as the divergence time.\n",
    "\n",
    "6. Return the divergence time array.\n",
    "\n",
    "Test out your function with a moderate number of \"pixels\" in both the `x` and `y` directions (e.g. a few hundred). It should take just a few seconds to run. Note: there should be **no loops**, except for the one for `N_iter`.\n",
    "\n",
    "Plot the result with `imshow()` - explore different choices for `max_val`, and create a big image!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}