03NumpyAndPlotting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/weymouth/NumericalPython/blob/main/03NumpyAndPlotting.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "# Numerical Python and Plotting\n",
    "\n",
    "We have established the main features of Python, but up to this point we have only created simple functions and applied them to create lists of numbers. We need to extend the base-language of Python for more advanced numerical work because:\n",
    " - There are no built-in data structures for arrays, matrices, and tables (unlike, say, `Matlab` or `Julia`).\n",
    " - Using lists of `float` numbers is generally very slow and lacks useful built-in features like matrix multiplication.\n",
    " - There is no built-in method to visualize data in plots.\n",
    "\n",
    "This notebook will introduce the `NumPy` and `PyPlot` libraries to address these issues."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ivbCdFpECLl8"
   },
   "source": [
    "# NumPy\n",
    "\n",
    "The numerical python, or [NumPy](https://numpy.org/), library enables fast and simple numerical methods in Python. To starting using this library (or any other) we need to use a new python keyword `import`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NMoD8fv_B_J3"
   },
   "outputs": [],
   "source": [
    "import numpy as np # importing numpy\n",
    "\n",
    "np.set_printoptions(precision=3)   # This sets numpy printing precision \n",
    "np.set_printoptions(suppress=True) # Don't use scientific notation by default"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us access to all the methods and functions in `NumPy` using the short name `np`. \n",
    "\n",
    "In the second and third line I've used the method `set_printoptions` to keep the notebook output pretty. Note that the format is `np` and then `.` and then the function name - just like when we used built-in methods for strings and lists in the previous notebooks. \n",
    "\n",
    "There are [far too many](https://numpy.org/doc/stable/reference/routines.html) new methods available to go through in this introduction, but most can be grouped into a few basic categories\n",
    "\n",
    "| Category       | Sub module   | Description                                                 |\n",
    "|----------------|--------------|-------------------------------------------------------------|\n",
    "| math           | numpy        | Scientific operations like $\\sqrt{a},\\log(a),\\sin(a)$, etc |\n",
    "| arrays         | numpy        | Array and matrix creation, and array operations like multiplication |\n",
    "| linear algebra | numpy.linalg | Matrix decomposition and solving linear systems |\n",
    "| fft            | numpy.fft    | Discrete Fourier Transform (of many types) and their inverse |\n",
    "| random sampling| numpy.random | Sample from different random variable distributions |\n",
    "\n",
    "<span style=\"display:none\"></span>\n",
    "\n",
    "Just as an example let see what is in the `numpy.random` submodule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"it contains methods such as...\",dir(np.random)[30:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get help on one of those. Notice we use the dot notation twice to access a submodule method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "?np.random.randint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All NumPy functions are well documented like this; with all the argument explained and useful examples. Using this function we can generate a sample of what might happen if we, say, roll a 20 sided-dice 4 times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = np.random.randint(1,21,size=4)\n",
    "samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every time you run this code, the result will be different. Try it!\n",
    "\n",
    "## Arrays and matrices\n",
    "\n",
    "Note that `samples` is a new data type called `array`. This is the building block for everything in NumPy. The easiest way to make arrays from scratch is to pass a list to `np.array`. Arrays can also have two dimensions, defining _a matrix_, or more. These can be created by passing a list of lists to `np.array`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_array = np.array(range(4))\n",
    "r_matrix = np.array([[1,2],[3,4],[5,6]])\n",
    "print(r_array)\n",
    "print('')\n",
    "print(r_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looping and slicing works on arrays just like it did on lists. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in r_array:\n",
    "    print(x)\n",
    "print('')\n",
    "print(r_array[-2:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It also works on matrices (and higher dimensional arrays), but there are more options since there are more indices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for row in r_matrix: # row by row\n",
    "    print(row)\n",
    "    for element in row: # element by element\n",
    "        print(element)\n",
    "print('')\n",
    "print(r_matrix[0,0]) # first element\n",
    "print(r_matrix[0])   # first row\n",
    "print(r_matrix[:,0]) # first column"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are also a few handy built-in functions to create arrays and matrices. See if you can guess what these functions will create before running the code block below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.linspace(0,1,6))\n",
    "print('')\n",
    "print(np.ones([3,2]))\n",
    "print('')\n",
    "print(np.eye(4))\n",
    "print('')\n",
    "print(np.diag([1,2,3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `numpy.linspace` function is especially useful as it makes an array equally spaced floats between the values you specify.\n",
    "\n",
    "Get some practice with this by defining your own matrix and slicing it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1. Define a 3x3 diagonal matrix called `my_matrix` with 2 along the diagonal in two different ways.\n",
    "#2. Slice this matrix to get the array [0 0 2] in two different ways."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Array operations and broadcasting\n",
    "\n",
    "There are a number of built-in NumPy operations that only make sense to apply to arrays vectors and matrices. Things like the transpose, matrix multiplication, and vector products. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1,2,3])\n",
    "b = np.array([4,5,6])\n",
    "print(r_matrix.T)      # transpose\n",
    "print(r_matrix.T @ a)  # matrix multiplication\n",
    "print(np.inner(a,b))   # inner product\n",
    "print(np.cross(a,b))   # cross product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are certainly useful in many situations. However, the _really_ nice thing about NumPy arrays is that they automatically and efficiently apply simple operations to each element of the array. Contrast this to lists, which required us to use list-comprehensions to get the same results!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_a,list_b = [1,2,3],[4,5,6] # lists, not arrays\n",
    "\n",
    "# list operators give unwanted \"string-like\" result\n",
    "print(2*list_a+list_b)        \n",
    "\n",
    "# list-comprehensions give desired result, but add complexity\n",
    "print([2*a+b for a,b in zip(list_a,list_b)])\n",
    "\n",
    "# Numpy arrays give desired result automatically!\n",
    "print(2*a+b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is called _operator broadcasting_ (or vectorizing) and is almost always the behavior we want in numerical programming. \n",
    "\n",
    "And this isn't limited to simple operations like `*,+`: NumPy has a slew of built-in scientific functions that will broadcast. Guess what these statements will print before running:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.sqrt(9))\n",
    "print(np.sqrt(r_array))\n",
    "print(np.log(np.e))\n",
    "print(np.log(samples))\n",
    "print(np.sin(np.pi/2))\n",
    "print(np.sin(r_matrix*np.pi/4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that NumPy defines some important constants like $e,\\pi$ that we can use as well. \n",
    "\n",
    "Get some more practice by completing the exercises below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#3. Create an array `c` which is my_matrix times the array `a`.\n",
    "#4. Look up the help for np.linalg.solve and use it to solve the system `my_matrix * x = c`\n",
    "#5. Use assert to make sure `c` is `2a` and `x` is `a`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Array Functions\n",
    "\n",
    "Not surprisingly, you can write functions for arrays as well, building on any of the original variable operators or NumPy operators. As a first trivial example, let's copy-paste two of our old functions from previous notebooks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_odd(n):\n",
    "    return n%2==1\n",
    "def double_then_add(a,b):\n",
    "    return 2*a+b\n",
    "\n",
    "print(is_odd(a))\n",
    "print(double_then_add(a,b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that operator broadcasting has applied automatically. _We didn't have to change anything!_\n",
    "\n",
    "As a second example let's create a few functions that wouldn't make sense on individual numbers. Let's rotate a point $p$ in 2D space by an angle $q$ around the origin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotation_matrix(q):\n",
    "    return np.array([[np.cos(q),-np.sin(q)],\n",
    "                     [np.sin(q), np.cos(q)]])\n",
    "\n",
    "def rotate_point(a,q):\n",
    "    return rotation_matrix(q) @ a\n",
    "\n",
    "a = np.array([1,1])\n",
    "for q in np.linspace(0,np.pi,9):\n",
    "    print(\"q={:.3g} rad, new p={}\".format(q,rotate_point(a,q)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new `rotate_point` function seems to have worked but it is a little hard to tell... \n",
    "\n",
    "# PyPlot\n",
    "\n",
    "Even for the previous simple example, it would be much easier to check the results if we could plot them. [Matplotlib](https://matplotlib.org/) is the most developed plotting library in python. `Matplotlib.PyPlot` are a collection of functions to make plot copied after the popular Matplot interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt   # import plotting library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note this `from-import` syntax lets you pick specific submodules or even specific functions to import. This can help with code readability. \n",
    "\n",
    "As a first example, let's just plot $\\sqrt{x}$ from $x=0...1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0,1)    # array of x values\n",
    "plt.plot(x,np.sqrt(x))  # plot x vs it's sqrt\n",
    "plt.xlabel('x')         # label x-axis\n",
    "plt.ylabel('$\\sqrt{x}$',rotation=0) # label y-axis (don't rotate it)\n",
    "plt.show()              # show the plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example shows that you can control pretty much everything in PyPlot, the trick is to just build the figure up one element at a time. \n",
    "\n",
    "Lets use PyPlot to test the `rotate_point` function above. We'll exactly the same loop as before, but instead of printing, we will add each rotated point to a scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for q in np.linspace(0,np.pi,9):\n",
    "    x,y = rotate_point(a,q)\n",
    "    plt.scatter(x,y,label='q={:.3g}'.format(q)) # plot point and assign label\n",
    "plt.legend()               # show legend\n",
    "plt.xlabel('x')            # label x-axis\n",
    "plt.ylabel('y',rotation=0) # label y-axis (don't rotate it)\n",
    "plt.axis('equal')          # scale the x,y axis equally \n",
    "plt.show()                 # show the plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The points form a half-circle starting from our original point $a=[1,1]$, visually confirming that the `rotate_point` function is working. \n",
    "\n",
    "Let's develop a function to create a polar plot of a given function $r(q)$ as a more advanced example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def polar_plot(func,q=np.linspace(0,2*np.pi)):\n",
    "    r = func(q)                # evaluate the function on q array\n",
    "    plt.plot(r*np.cos(q),r*np.sin(q)) # compute and plot x,y arrays\n",
    "    plt.xlabel('x')            # label x-axis\n",
    "    plt.ylabel('y',rotation=0) # label y-axis (don't rotate it)\n",
    "    plt.axis('equal')          # scale the x,y axis equally \n",
    "    plt.show()                 # show the plot\n",
    "\n",
    "def cardiod(q): return 2*(1-np.cos(q))\n",
    "polar_plot(cardiod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the function `polar_plot` doesn't `return` anything - it just `show`s a plot. The function `polar_plot` is also interesting because its first argument is a *function*, not a number or an array. This is common enough in Python that there is a [lambda](https://realpython.com/python-lambda/) syntax to create functions _in-place_ instead of giving the function a name and then using it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polar_plot(lambda q: q/np.pi) # anonymous function to make a spiral"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `lambda` syntax is very handy when doing more advanced analysis like function optimization, root-finding, integration, etc. Try making a more interesting polar plot, like $2+\\sin(10q)$ or any other function you like.\n",
    "\n",
    "\n",
    "# Failure Modeling Example\n",
    "\n",
    "Let's finish with an engineering simulation. Let's model the potential failure of a fuel oil pump using a simplistic model. \n",
    "\n",
    "First we'll model the pressure upstream of the pump as a random walk: at each time step the pressure can go up, down or stay constant with equal probability. We could do this one random step at a time, but it is faster to generate all the steps at once, and then add them together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pressure_walk(n_steps):\n",
    "    pressure_steps = np.random.randint(-1,2,size=n_steps)\n",
    "    return np.cumsum(pressure_steps) # cumulative summation of the steps\n",
    "\n",
    "for i in range(5):\n",
    "    plt.plot(pressure_walk(200))\n",
    "plt.xlabel('time'); plt.ylabel('pressure')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the cumulative summation function to get the signal from the steps. The plot shows 5 random pressure histories and every time you run the block above the histories will change.\n",
    "\n",
    "Now let's think about the FO pump. A pump only has a finite range of operation - if the pressure is too low, the pump will immediately stop working. So lets write a function to check if the pump fails, given a pressure signal. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pump_failure(pressure,lower=-15):\n",
    "    for i,p in enumerate(pressure):\n",
    "        if p<=lower: return i\n",
    "    return len(pressure)\n",
    "\n",
    "pressure_test = np.zeros(10) # zero array for testing\n",
    "i = np.random.randint(1,11)  # pick a failure time\n",
    "pressure_test[i:] = -20      # replace test values after this time \n",
    "print(pressure_test)\n",
    "print(\"Failure time = {}\".format(pump_failure(pressure_test)))\n",
    "assert(i==pump_failure(pressure_test))\n",
    "print(\"Test passed!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `enumerate` function iterates through the index and values in a list or array. This lets us test if the pressure limit has been crossed and return the index when that happens.\n",
    "\n",
    "The code below the function is a test to see if the function is working properly. [Writing little tests like this lets you catch errors before integrating functions together!](https://docs.python-guide.org/writing/tests/)\n",
    "\n",
    "Finally, we can use these two function to get a feeling for how quickly our pump is likely to fail. Let's run our code a few thousand times and plot the results as a histogram. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_runs = 10000\n",
    "n_steps = 1000\n",
    "failure_data = [pump_failure(pressure_walk(n_steps)) for i in range(n_runs)]\n",
    "plt.hist(failure_data,bins=range(1,n_steps,10))\n",
    "plt.xlabel('time of pump failure'); plt.ylabel('number of failures')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.count_nonzero(np.array(failure_data)==n_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that most failures happen between 100 and 200 time steps, but there is a \"long tail\" of pumps that last much much longer. \n",
    "\n",
    "## Additional exercises\n",
    "\n",
    "1. Our random walk function isn't very realistic. Find a function in the `numpy.random` submodule to generate steps from a Gaussian distribution instead. How does this change the signals and the failure results?\n",
    "1. A pump can also fail if it is exposed to excessively high pressure for too long and becomes damaged. Write a function to return the accumulated time the pump spends above a threshold value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyPPRJif4V4H6mucN5PMLqW1",
   "include_colab_link": true,
   "name": "03NumpyAndPlotting.ipynb",
   "provenance": []
  },
  "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}