{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": true }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> All content here is under a Creative Commons Attribution [CC-BY 4.0](https://creativecommons.org/licenses/by/4.0/) and all source code is released under a [BSD-2 clause license](https://en.wikipedia.org/wiki/BSD_licenses). Parts of these materials were inspired by https://github.com/engineersCode/EngComp/ (CC-BY 4.0), L.A. Barba, N.C. Clementi.\n", ">\n", ">Please reuse, remix, revise, and [reshare this content](https://github.com/kgdunn/python-basic-notebooks) in any way, keeping this notice." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Module 5: Overview \n", "\n", "In the [prior module](https://yint.org/pybasic04) you created vectors, matrices and arrays, with numbers in them. Now you get to perform calculations.\n", "\n", "* addition, subtraction, multiplication, division\n", "* averages, minimum, maximum: over all entries, or just the rows and columns separately.\n", "* matrix multiplication, transposing and inverses.\n", "\n", "You will learn these by discovery: copy and paste the code, but look at the code first. Often there are other interesting ideas from prior modules, or even other new features we will learn about later.\n", "\n", "### Preparing for this module\n", "\n", "

Start a new version controlled repository for your work. Put each section in a new Python script file. Commit your work at the end of each section." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Accessing and changing entries in an array\n", "\n", "Learn by discovering...\n", "\n", "```python\n", "import numpy as np\n", "\n", "# Create an array with 3 rows and 4 columns\n", "x = np.array([[1, 2, 3, 4],\n", " [5, 6, 7, 8],\n", " [9, 10, 11, 12]])\n", "\n", "# Access individual values:\n", "print(x[0, 0]) # first\n", "print(x[2, 3]) # last\n", "\n", "# Get vector \"slices\"\n", "print(x[0, :]) # Row 1\n", "print(x[1, :]) # Row 2\n", "print(...) # Column 3\n", "\n", "# Get sub-arrays:\n", "print(x[1:3, 1:3]) # what section is this?\n", "print(x[1:, 1:3]) # and this?\n", "print(x[1:, 1:]) # and this?\n", "\n", "# Negative indexing is also allowed, like we saw with lists.\n", "print(x[-2:, :]) # The last 2 rows\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can also overwrite entries in the array:\n", "```python\n", "# Create an array with 3 rows and 4 columns\n", "x = np.array([[1.1, 2, 3, 4],\n", " [5, 6, 7, 8],\n", " [9, 10, 11, 12]])\n", "\n", "# Change individual values:\n", "x[0, 0] = 0\n", "x[1, 2] = 3 \n", "print(x)\n", "\n", "# Change an entire slice in row 2:\n", "x[1, :] = 2\n", "\n", "# Change the last column to be all -4\n", "x[:, -1] = -4\n", "\n", "# Reset the entire matrix to infinty!\n", "x[:, :] = np.inf\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simple elementwise functions and operations on a NumPy array\n", "\n", "Once we have created an array - [see the prior notebook](https://yint.org/pybasic04) - we are then ready to actually use them for calculations!\n", "\n", "Let us consider these calculations:\n", "1. Addition and subtraction\n", "2. Multiplication and division (element-by-element)\n", "3. Square roots and other powers\n", "4. Trigonometric and other functions\n", "\n", "

Commit your prior work, and start a new Python script for this section." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. Addition and subtraction\n", "\n", "Create 2 matrices, of the same shape, and let's get started:\n", "\n", "```python\n", "import numpy as np\n", "\n", "A = np.ones(shape=(5,5))\n", "B = np.ones(shape=(5,5))\n", "print('A = \\n{}\\n\\nB = \\n{}'.format(A, B))\n", "\n", "print(A + B)\n", "```\n", "\n", "The ``+`` operation on two arrays is actually just a convenience. The actual function in NumPy which is being called to do the work is the ``np.add(...)`` function. \n", "```python\n", "np.add(A, B) # does exactly the same\n", "\n", "```\n", "\n", "Similarly, we have the ``-`` and ``np.subtract()`` functions that serve the same purpose:\n", "\n", "```python\n", "print(A - B)\n", "print(np.subtract(A, B))\n", "print(np.add(A, -B)) \n", "```\n", "\n", "These are called ***element-by-element operations***. That means, NumPy performs the operation of addition on each corresponding element in the arrays `A` and `B` and then repeats that entry-by-entry. This is also called elementwise in NumPy's documentation.\n", "\n", "NumPy will also allow you take shortcuts. Imagine that you want to subtract the value of 3.5 from every entry in matrix `A`. You do not first need to create a matrix with the same shape as ``A`` with all entries containing the values of 3.5, and then go to subract that. \n", "\n", "**There is a shortcut:**\n", "\n", "```python\n", "# Also does element-by-element calculations\n", "print(A - 3.5) \n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Multiplication and division (element-by-element)\n", "\n", "Multiplication and division can also be done element-by-element using the ``*`` and ``/`` operators.\n", "\n", "Try these to learn something new:\n", "> Create a matrix with the numbers 1 to 25 inclusive, ``reshape``d into 5 rows and 5 columns. \n", ">\n", "> 1. Divide each entry by 10.\n", "> 2. Now multiply each entry by ``np.pi`` $\\approx 3.14159$\n", "> 3. Start from the same matrix, and divide each entry by ``np.e`` $\\approx 2.71828$\n", "> 4. Try dividing by zero, in regular Python (i.e. not with NumPy): ``5.6 / 0``\n", "> 5. Now try dividing a NumPy matrix by zero: why do you get a different result? Does the NumPy result make more sense, or less sense, than regular division by zero in Python?\n", "\n", "***Important note for MATLAB users***\n", "\n", "Matlab is unique among programming languages for overloading the ``*`` operator for matrix multiplication. In MATLAB that requires the matrices to the left and right of ``*`` to have a consistent shape. We will see the NumPy equivalent of regular matrix multiplication further down.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3. Square roots and other powers\n", "\n", "There are other elementwise operations that can be done on matrices. These often involve raising the individual matrix elements to a certain power, or taking a square root (which is the same as raising a number to the power of $0.5$), or calculating the logarithm.\n", "\n", "Let's try it out interactively.\n", "\n", "> 1. Use the ``**`` operation to raise a matrix to a power\n", "> 2. Use the ``.square()`` function \n", "> 3. Use the ``.power()`` function \n", "> 4. Use the ``.sqrt()`` function\n", "> 5. Verify that `**(0.5)` gives the same values as the `.sqrt()` function\n", "\n", "Extend your knowledge: try the above on a vector or matrix that contains values which are negative, zero and positive. What do you notice about the square root of a negative number?\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4. Trigonometric and other functions\n", "\n", "A wide variety of mathematical functions are possible to be calculated on NumPy arrays. See the full list in the [NumPy documentation](https://docs.scipy.org/doc/numpy/reference/routines.math.html).\n", "\n", "You will self-discover these function by running the code below.\n", "\n", "**Try the following on this matrix:**\n", "\n", "```python\n", "# A 4x4 matrix with positive and negative values:\n", "values = np.reshape(np.linspace(-2, +2, 16), (4, 4)) \n", "print(values)\n", "```\n", ">1. Try using the standard trigonometric functions ``np.sin(...)``, ``np.tan(...)``, etc on that matrix.\n", ">2. Rounding each value to the closest integer. Use the ``np.around()`` function. \n", ">>Which way do negative values round?\n", ">3. Round each value to a certain number of ``decimals``; try rounding to 1 decimal place. \n", ">>Are the results what you expect?\n", ">4. Similar to rounding: try the ``np.floor(...)`` and ``np.ceil(...)``: what is the difference between the floor and the ceiling? Hint: read the documentation for [`floor`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.floor.html) and [`ceil`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ceil.html).\n", ">> Watch for negative values: what do you expect ``np.floor(3.5)`` to do? And ``np.floor(-3.5)``\n", ">5. Logarithms and exponents are also part of the standard calculations we expect to do with matrices using the ``np.log(...)`` and ``np.exp(...)`` functions. Recall that $\\log(\\exp(x)) = x$. \n", ">\n", "> But this might be surprising to you:\n", "> ```python\n", "exponent = np.exp(values)\n", "print(exponent)\n", "recovered = np.log(exponent)\n", "print(recovered) \n", "print('-----')\n", "print(recovered - values)\n", "```\n", "> Does ``recovered`` match the original ``values`` matrix? It should: we first took the exponent, then the logarithm. This subtraction should be a matrix of all zeros.\n", ">\n", "> The last matrix in your printout above should be all zeros, but is not exactly equal to zero (it is very, very close to zero though).\n", ">\n", "> To test that we can use the ``np.isclose(...)`` function. It is another elementwise function that you can add to your toolbox. It tests if the entries in an array are close to another:\n", "```python\n", "np.isclose(recovered - values, 0)\n", "```\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## ➜ Challenge yourself: using NumPy and elementwise operations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "\n", "This problem uses all the skills you have learned so far: creating NumPy arrays from a list, and doing calculations on them. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back in [module 1](https://yint.org/pybasic01) we saw that $$n! \\approx \\sqrt{2 \\pi n} \\, \\cdot \\, n^n \\, \\cdot \\, e^{-n} $$ \n", "\n", "is what is known as the Stirling approximation for the factorial of an integer value $n$. \n", "\n", "1. Create a vector of **floating point values**, from 1.0 to 15.0 in NumPy, representing 15 values of $n$.\n", "2. Create another vector which calculates that Stirling approximation, elementwise, for every value of $n$ in your vector. The code for this should be a single line. See the template below.\n", "3. Create an empty list which will hold the true values of the factorial $n!$ and fill that vector in a for-loop. Append the true factorial values for $n = 1, 2, \\ldots N$ to that list. Remember we covered appending entries into a list in a prior module.\n", "4. Convert that Python list to a NumPy vector.\n", "5. Now you have 2 NumPy vectors: the true values, and the approximate values. Subtract the 2 vectors to calculate the percentage error. Does the percentage error get worse or better for larger values of $n$?\n", "\n", "\n", "\n", "```python\n", "N = 15\n", "\n", "# 1. Create a vector of values first: 1, 2, ... N\n", "n_values = np.linspace(...)\n", "\n", "# 2. Find Stirling's approximation for the factorial\n", "approximate = np.sqrt(...) * ....\n", "\n", "# 3. True value of n! (to compare with later)\n", "import math\n", "true_values = []\n", "for n in np.arange(N):\n", " # Use math.factorial(n) to calculate true value\n", " # Append that to the list\n", " ...\n", "\n", "# 4. Convert list to a NumPy array \n", "\n", "# 5. Finally, calculate and show the percentage relative error\n", "error = ...\n", "percentage_error = ...\n", "```\n", "\n", "***Observe:*** how the problem above (to discover the approximation error) was broken down into 5 sub-steps. Try doing that for your own programs. Create a new file and break the problem down:\n", ">```text\n", "1. Create a vector of values: 1, 2, ... N\n", "2. Find Stirling's approximation for the factorial\n", "3. Calculate true value of n!\n", "4. Convert list to NumPy array \n", "5. Calculate the relative error\n", "```\n", "\n", "then turn those lines into comments, and get started coding.\n", "\n", "

***Have you committed your script to your repo? Also push it to the remote server.***" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Operations on the entire matrix\n", "\n", "You have just obtained all the data in your matrix, and now you wish to find the largest, or smallest value in the entire matrix. Or perhaps the average of all numbers. These are global operations.\n", "\n", "```python\n", "import numpy as np\n", "rnd = np.array([[ 7, 3, 11, 12, 2], \n", " [10, 13, 8, 8, 2], \n", " [ 3, 13, 6, 2, 3], \n", " [ 5, 3, 9, 2, 6]])\n", "print('The matrix is:\\n{}'.format(rnd))\n", "\n", "max_value = np.amax(rnd)\n", "print('The maximum value is {}'.format(max_value))\n", "\n", "min_value = np.amin(rnd)\n", "print('The minimum value is {}'.format(min_value))\n", "\n", "# Note here how strings can be split across lines:\n", "print(('The mean value is {}, while the median '\n", " 'is {}').format(np.mean(rnd), np.median(rnd)))\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The ``np.amax(...)`` and ``np.amin(...)`` functions will find the maximum or minimum in the entire array: looking at every single entry in all dimensions.\n", "\n", "### Enrichment: (you can skip this for now)\n", "\n", "The NumPy library will internally unfold, or flatten the array into a single long vector. Take a look at what that looks like when you use the ``.flatten(...)`` method on the array: ``rnd.flatten()``. It works by going from column to column, down each row:\n", "\n", "```python\n", "print(rnd.flatten())\n", "```\n", "``[ 7 3 11 12 2 10 13 8 8 2 3 13 6 2 3 5 3 9 2 6]``\n", "\n", "This is actually how the data is stored internally in the computer's memory.\n", "\n", "The reason we point this ``array.flatten(...)`` function out is because sometimes knowing what the maximum value is is only half the work. The other half is knowing *where* that maximum value is. For that we have the ``np.argmax(...)`` function.\n", "\n", "Try this code:\n", "\n", "```python\n", "max_position = np.argmax(rnd)\n", "print(('The maximum value is in position {} of the '\n", " 'flattened array').format(max_position))\n", "```\n", "\n", "Verify that that is actually the case, using the space below to run the code. The `max_position` is the position in the ``flatten``ed array, which is why you should know about the ``array.flatten(...)`` function.\n", "\n", "Another tip: notice that the maximum value of ``13`` appears twice. To which one does ``max_position`` refer?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Other operations on the entire matrix, which reduce the entire matrix down to a single number. Try these out:\n", "```python\n", "# Related to maximum and minimim is `ptp`. Use \n", "# the help function to see what this does.\n", "rnd.ptp() \n", "\n", "rnd.mean()\n", "rnd.prod() # what does this do?\n", "rnd.std() \n", "rnd.sum() \n", "rnd.trace() # sum of the elements on the diagonal\n", "rnd.var()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Functions on rows or columns of a matrix\n", "\n", "Above you learned about elementwise operations. In other words, NumPy performed the mathematical calculation on every element (entry) in the array. You also saw functions which summarize the matrix by 1 single value (like the mean or standard deviation).\n", "\n", "Sometimes we need calculations that work on every row, or column, of an array. We will cover:\n", "1. Calculate the minimum value in every row (give back a column vector that has the minimum value of every row).\n", "2. Calculate the average value of every column (give back a row vector that has the average value of every column).\n", "3. Sort the columns (or rows), so every column (or row) goes from low to high.\n", "4. Show the (cumulative) sum going across each column (or row).\n", "\n", "We will talk about matrices in these examples, but the operations can also be applied to multi-dimensional arrays, with 3, 4, or more dimensions.\n", "\n", "> We also introduce 2 important terms: ***``axis``*** and ***``inplace``***, both of which you will regularly see in the NumPy documentation, but also later in Pandas.\n", "\n", "

Before you get started, commit prior work, and start this section in a new script." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. Minimum (or maximum) value in every row/column\n", "\n", "Think of a matrix containing the daily temperatures per city; one city per column, and every row is a day of the year. \n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
City 1City 2City 3City 4
Day 1791210
Day 21452
Day 3-31-2-3
Day 4-2-1-2-2
Day 5-3-1-2-4
\n", "\n", "\n", "* What is the max or min temperature for each city (i.e. per column)?\n", "* What is the max or min temperature each day for all cities (i.e. per row)?\n", "\n", "For this we also use the ``np.amax(array, axis=...)`` or ``np.amin(array, axis=...)`` function, which you saw above. But, now you must specify, as a second input, along which ***``axis``*** you want that extreme value to be calculated. \n", "* Axis 0 is the first axis, down each column, along the direction of the rows, going from top to bottom.\n", "* Axis 1 is the next axis, across each row, along the direction of the columns, going from left to right.\n", "* Axis 2, for arrays, is the next axis, the 3rd dimension.\n", "\n", "Try out the code below. Note that the array is created from a list-of-lists.\n", "```python\n", "\n", "import numpy as np\n", "temps = np.array([[7, 9, 12, 10], \n", " [1, 4, 5, 2], \n", " [-3, 1, -2, -3], \n", " [-2, -1, -2, -2], \n", " [-3, -1, -2, -4]])\n", "print(('The temperatures are given one column per '\n", " 'city, each row is a daily average '\n", " 'temperature:\\n{}').format(temps))\n", "\n", "max_value_0 = np.amax(temps, axis=0)\n", "print(('The maximum value along axis 0 (row-wise, '\n", " 'per city for all days) is {}').format(max_value_0))\n", "\n", "max_value_1 = np.amax(temps, axis=1)\n", "print(('The maximum value along axis 1 (column-wise, '\n", " 'per day for all cities) is {}').format(max_value_1))\n", "\n", "# Notice the above output is 'flatten'ed and returned as a row, \n", "# instead of a column, as you might hope for. We can use \n", "# the `keepdims` input though:\n", "max_value_1_col = np.amax(temps, axis=1, keepdims=True)\n", "print(('The maximum value along axis 1 (column-wise, '\n", " 'per day for all cities) is\\n{}').format(max_value_1_col))\n", "```\n", "\n", "You can visually verify that the maximum values returned are what you expected.\n", "\n", "Now try it for the minimum values of the matrix:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Average value in every row/column\n", "\n", "\n", "Many functions in NumPy take ***``axis``*** as in input argument, including the ``np.mean(...)``, ``np.std(...)`` and other functions you saw above. \n", "\n", "Complete this code:\n", "```python\n", "average_temperature_per_city = ...\n", "print(('The average temperature for each city '\n", " 'is: {}').format(average_temperature_per_city))\n", "\n", "average_temperature_per_day = ...\n", "print(('The average temperature per day for all '\n", " 'cities: {}').format(average_temperature_per_day))\n", "\n", "std_per_city = ...\n", "print(('The standard deviation of the temperature '\n", " 'is {}').format(std_per_city))\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3. Sorting every row/column\n", "\n", "Just like finding the [minimum or maximum](#2.--Average-value-in-every-row/column) value in every column, you might also be interested in sorting each column.\n", "\n", "We want to sort every column **independently** of the others. In other words every column will be sorted from low to high by the end. This is in contrast to sorting based on one column, and the other rows follow with. We will see that coming later.\n", "\n", "We have seen the ***``axis``*** input several times now, and here we will use it again to indicate which axis we would like to sort along.\n", "\n", "```python\n", "import numpy as np\n", "temps = np.array([[7, 9, 12, 10], \n", " [1, 4, 5, 2], \n", " [-3, 1, -2, -3], \n", " [-2, -1, -2, -2], \n", " [-3, -1, -2, -4]])\n", "\n", "sorted_columns = np.sort(temps, axis=0)\n", "print('The temperatures, sorted in each column, along axis 0: \\n{}'.format(sorted_columns))\n", "\n", "sorted_rows = np.sort(temps, axis=1)\n", "print('The temperatures, sorted in each row, along axis 1: \\n{}'.format(sorted_rows))\n", "\n", "print('Verify the original matrix is untouched:\\n{}'.format(temps))\n", "```\n", "\n", "The sort takes place and the result is set to a new variable. **This is not efficient**, especially if the input_array to be sorted is really large. It means that a copy of the data is made, using up memory, and computer time; and only then does the sort take place on the copy of the data.\n", "\n", "It is possible to simply sort the original array. This is called ***in-place sorting***. You will see that terminology in several places in NumPy's documentation: ***in-place***. It means that the original matrix is used, calculated on, and the result is left in the same variable as the original matrix.\n", "\n", "You perform an in-place sort as follows:\n", "```python\n", "input_array.sort(axis=...)\n", "```\n", "compared with a sort that is assigned to a new variable:\n", "```python\n", "output_array = np.sort(input_array, axis=...)\n", "```\n", "\n", "By definition an in-place operation means there is no need to assign the result to another variable as output. The second variation with ``np.sort(input_array, ...)`` follows good practice, where a function does not modify its inputs. So therefore the ``output_array`` is required.\n", "\n", "Let's see what the implication of in-place sorting is:\n", "\n", "```python\n", "import numpy as np\n", "temps = np.array([[7, 9, 12, 10], \n", " [1, 4, 5, 2], \n", " [-3, 1, -2, -3], \n", " [-2, -1, -2, -2], \n", " [-3, -1, -2, -4]])\n", "\n", "# In-place sort. We don't need to use `output=` but let's see what happens.\n", "# The in-place sort result is in the original variable \"temps\"\n", "output = temps.sort(axis=0) \n", "print('The sorted values along axis 0: \\n{}'.format(temps))\n", "print('Out of curiosity, the value of \"output\" is: {}'.format(output))\n", "\n", "# So you can simply say:\n", "temps.sort(axis=0)\n", "\n", "# and the result will be sorted in the original variable.\n", "```\n", "\n", "Note the above behaviour is consistent with what [we saw before in regular Python](https://yint.org/pybasic02)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4. Calculate the (cumulative) sum going across each row/column\n", "\n", "If the values in the rows or columns represent some property, such as weight in a container, then it might be interesting to calculate the cumulative value. \n", "\n", "We will create a matrix where the values in column 0 will all be ones. Values in column 1, 2, 3 and 4 will represent the weight (kilograms), added to 4 containers. Each row represents 1 minute. Column 0 is a counter. \n", "\n", "**The goal**: Find the point in time when the weight in the container just exceeds 100kg. You will see why we have a counter as column 1.\n", "\n", "```python\n", "import numpy as np\n", "n = 20\n", "k = 4\n", "\n", "counter = np.ones(shape=(n, 1)) \n", "weights = np.random.randint(low=4, high=9, size=(n, 4))\n", "\n", "# Stack the the objects side-by-side (h=horizontal)\n", "weight_matrix = np.hstack((counter, weights))\n", "\n", "print('The weight added to the {} containers every minute [ignore first column]:\\n{}'.format(k, weight_matrix))\n", "\n", "accumulation = np.cumsum(weight_matrix, axis=0)\n", "print('The cumulative weight over time is:\\n{}'.format(accumulation))\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5. Summary\n", "\n", "The methods to calculate something globally, or something on just a row or a column are the same methods: just add the ``axis`` keyword: \n", "* ``np.sum`` is over the entire array\n", "* ``np.sum(axis=...)`` is over the specified axis.\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
OperationDescription
np.sumSum of entries
np.prodProduct (multiplication) of all entries
np.meanArithmetic average
np.stdStandard deviation
np.varVariance
np.minMinimum value
np.maxMaximum value
np.argminIndex of the minimum
np.argmaxIndex of the maximum
np.medianMedian value
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## ➜ Challenge yourself: summarizing values from a matrix\n", "\n", " \n", "\n", "Use the skills just learned above regarding calculations, and in the [prior module](https://yint.org/pybasic04) regarding random numbers, to solve the following 2 problems.\n", "
\n", "\n", "

Complete this challenge in a new script. Have you pushed your prior work?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Tossing a coin can be seen as a random process that generates either a 0 for tails, or a 1 for heads. If the coin is fair, then the 0's or 1's are from a uniform distribution (equal chance of occurring).\n", "\n", " Simulate $n=[40, 400, 4000, 40000, 4\\times 10^6]$ coin tosses in Python, in a single NumPy vector. For each one of these 5 scenarios, calculate the average of the vector. \n", " \n", " * What should the theoretical average be for the vector of coin tosses, if it is a fair coin?\n", " * What average values do you actually achieve for each scenario above?\n", " * Repeat the simulations a few times, to get a feeling for repeatability.\n", " * How would you simulate a coin that gives 55% chance of landing on heads?\n", " \n", "2. The Dutch meteorological service (KNMI) makes historical weather data available on their website, [by the hour](https://projects.knmi.nl/klimatologie/uurgegevens/selectie.cgi). One particular data set available is the [homogenized temperature](http://projects.knmi.nl/klimatologie/onderzoeksgegevens/homogeen_260/index.html), which corrects the temperature readings for changes in the measuring station over time. This data set is [publicly available here](http://projects.knmi.nl/klimatologie/onderzoeksgegevens/homogeen_260/tg_hom_mnd260.txt). \n", "\n", " The first column is the station number, the next is the date, and the third column is the temperature, measured in units of °C.\n", "\n", " * Download the dataset and read in the data file by completing the following template:\n", " > ```python\n", " filename = 'KNMI-Homogenized-temperature-monthly-average-1901-2019.txt'\n", " full_filename = os.path.join(...)\n", " weather = np.loadtxt(full_filename, skiprows=27) \n", " ```\n", " * What is the lowest temperature recorded over that time period? And the highest?\n", " * What is the standard deviation of the temperature values?\n", " * What is the average of the first 700 rows in the data set? And the last 700 rows?\n", " * Repeat for the standard deviation as well. \n", " * What is the average temperature since the year you were born? And prior?\n", "\n", " In future modules we will plot and examine other statistics from this data file." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Matrix operations: multiplication, transpose, inverse\n", "\n", "Matrices can be considered like spreadsheets containing rows and observations. But the real power of matrices comes from some of the calculations we can do with them. \n", "\n", "* We have already seen [elementwise operations above](#Simple-elementwise-functions-and-operations-on-a-NumPy-array): addition and subtraction, as well as elementwise multiplication and division.\n", "* Here we consider ***matrix multiplication***, which is different from elementwise multiplication.\n", "* Matrix transponses (to flip a matrix around)\n", "* Matrix inverses (related to, but not the same as division)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can multiply a matrix, called $\\mathbf{A}$ that has 4 rows and 2 columns, with a second matrix, $\\mathbf{B}$, which contains 2 rows and 3 columns. The product of $\\mathbf{A}$ and $\\mathbf{B}$, or also written as $\\mathbf{AB}$ is another matrix, we have called it $\\mathbf{X} = \\mathbf{AB}$. That matrix $\\mathbf{X}$ has 4 rows and 3 columns.\n", "\n", "The value in first row and first column of $\\mathbf{A}$ is given by the variable $a_{1,1}$; the value in row 1 and column 2 is $a_{1,2}$, and so on. The entry in row $i$ and column $j$ is therefore given by $a_{i,j}$. We use a similar notation for matrix $\\mathbf{B}$ and $\\mathbf{X}$. You can see these values in the figure shown here. Figure credit: [Wikipedia](https://en.wikipedia.org/wiki/Matrix_multiplication#/media/File:Matrix_multiplication_diagram_2.svg).\n", "\n", "\n", "$$\\begin{eqnarray*}\\mathbf{A}\\mathbf{B} &=& \\mathbf{X}\\\\\n", "\\begin{bmatrix}\n", "{a_{1,1}} & {a_{1,2}} \\\\\n", "{a_{2,1}} & {a_{2,2}} \\\\\n", "{a_{3,1}} & {a_{3,2}} \\\\\n", "{a_{4,1}} & {a_{4,2}} \\\\\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "{b_{1,1}} & {b_{1,2}} & {b_{1,3}} \\\\\n", "{b_{2,1}} & {b_{2,2}} & {b_{2,3}} \\\\\n", "\\end{bmatrix}\n", "&=& \\begin{bmatrix}\n", "x_{1,1} & x_{1,2} & x_{1,3} \\\\\n", "x_{2,1} & x_{2,2} & x_{2,3} \\\\\n", "x_{3,1} & x_{3,2} & x_{3,3} \\\\\n", "x_{4,1} & x_{4,2} & x_{4,3} \\\\\n", "\\end{bmatrix}\\end{eqnarray*}\n", "$$\n", "\n", "The values at the intersection of $\\mathbf{X}$ are marked with circles:\n", "\n", "$\\begin{align*}\n", "x_{1,2} & = {a_{1,1}}{b_{1,2}} + {a_{1,2}}{b_{2,2}} \\\\\n", "x_{3,3} & = {a_{3,1}}{b_{1,3}} + {a_{3,2}}{b_{2,3}}\\\\\n", "\\end{align*}$\n", "\n", "Notice how the dimensions are consistent: a vector of 1 row and 2 columns $[{a_{1,1}} \\,\\,\\, {a_{1,2}}]$ is multiplied element-by-element with $[{b_{1,2}}\\,\\,\\, {b_{2,2}}]$, and then the multiplied values are added. That process is repeated over and over, until all the values in the new matrix $\\mathbf{X}$ is formed.\n", "\n", "There is also a pattern: the value that goes into row 1 and column 2 of $\\mathbf{X}$ comes from multiplying row 1 of $\\mathbf{A}$ and column 2 of $\\mathbf{B}$. This form of multiplying the vectors element-by-element and then adding up the result is called the ***dot product***. This explains why the function to do matrix multiplication in NumPy is called ``np.dot(...)``.\n", "\n", "Just a quick mention: you cannot just multiply any two matrices $\\mathbf{A}$ and $\\mathbf{B}$. We do require that:\n", "* all values be present (no missing values), and \n", "* that $\\mathbf{A}$ must have the same number of columns as the number of rows in matrix $\\mathbf{B}$.\n", "\n", "\n", "There is a lot of multiplying and adding happening here; 12 dot products. This is certainly not something you want to do by hand regularly. So ... let's try it with NumPy.\n", "\n", "```python\n", "import numpy as np\n", "A = np.array([[8, 7], [6, 5], [4, 3], [2, 1]])\n", "B = np.array([[1, 2, 3], [4, 5, 6]])\n", "\n", "print('The shape of matrix A is: {}\\n'.format(A.shape))\n", "print('The shape of matrix B is: {}\\n'.format(B.shape))\n", "\n", "X = np.dot(A, B)\n", "print('The matrix multiplication of AB = X \\n{}'.format(X))\n", "\n", "# The value of X can also be calculated from:\n", "X = A.dot(B)\n", "\n", "print(('The value in position (2,3) of X comes from row 2 '\n", " 'of A [{}] and column 3 of B [{}] which corresponds '\n", " 'to 6*3 + 5*6 = {}').format(A[1,:], B[:,2], 6*3 + 5*6))\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The `@` operator for matrix multiplication\n", "\n", "Recently ([in 2014](https://www.python.org/dev/peps/pep-0465/)) Python introduced a new operator: `@` specifically defined for matrix multiplication. Depending on if your Python version is up to date, you do not need to write `np.dot(A, B)` anymore, but rather, more compactly: `A @ B` for matrix multiplication. You will still see `np.dot(A, B)` in older NumPy code, but ``A @ B`` is much cleaner.\n", "\n", "```python\n", "# Requires a recent Python version \n", "import numpy as np\n", "A = np.array([[8, 7], [6, 5], [4, 3], [2, 1]])\n", "B = np.array([[1, 2, 3], [4, 5, 6]])\n", "X = A @ B\n", "print('Compact matrix multiplication: A@B = \\n{}'.format(X))\n", "\n", "# which is the same result as given above.\n", "```\n", "\n", "***Try this important exercise:***\n", "> Let $\\mathbf{A}$ be a matrix with 4 rows and 3 columns.\n", "> Let $\\mathbf{B}$ be a matrix with 3 rows and 2 columns.\n", ">```python\n", "import numpy as np\n", "A = np.array([[8, 7, 3], [6, 5, 4], [4, 3, 5], [2, 1, 6]])\n", "B = np.array([[1, 2], [4, 5], [7, 8]])\n", "```\n", ">\n", ">Verify if $\\mathbf{AB}$ is equal to $\\mathbf{BA}$. *Hint*: we already know from the shapes of $\\mathbf{A}$ and $\\mathbf{B}$ that this is not true. Remember the rule about the number of rows and columns?\n", ">\n", "> Now try it when the number of rows in $\\mathbf{B}$ do match the number of columns in $\\mathbf{A}$:\n", "> ```python\n", "A = np.array([[8, 7], [6, 5]])\n", "B = np.array([[1, 2], [3, 4]])\n", "```\n", "Is it true in general that $\\mathbf{AB} = \\mathbf{BA}$, or equivalently: $\\mathbf{AB} - \\mathbf{BA} = \\mathbf{0}$?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The transpose of a matrix\n", "\n", "In the above example you saw an error message when the dimensions of the matrix don't match. For example, you have two matrices:\n", "* $\\mathbf{A}$ with 4 rows and 3 columns\n", "* $\\mathbf{B}$ with 2 rows and 3 columns\n", "\n", "and you want to multiply them to get $\\mathbf{X}=\\mathbf{AB}$ with with 4 rows and 2 columns.\n", "\n", "Matrix $\\mathbf{B}$ does is not in the right order: it has 2 rows and 3 columns, but it should have 3 rows and 2 columns\n", "\n", "We first need to first ***transpose*** $\\mathbf{B}$, which means to flip the rows and columns: rows become columns and columns become rows. \n", "\n", "In this example once we have transposed matrix $\\mathbf{B}$ it will have 3 rows now and 2 columns, and be ready to multiply. We therefore write: $\\mathbf{X}=\\mathbf{AB}^T$, or you sometimes also see $\\mathbf{B}'$ to indicate the transpose.\n", "\n", "There above you see the notation we use to indicate a transpose: \n", "* $\\mathbf{B}$ has 2 rows and 3 columns\n", "* $\\mathbf{B}^T$ has 3 rows and 2 columns.\n", "\n", "Transposing happens so frequently that there is shortcut for it in NumPy: ``B.T`` will transpose matrix ``B``.\n", "\n", "***Try it:***\n", "```python\n", "import numpy as np\n", "A = np.array([[8, 7, 3], [6, 5, 4], [4, 3, 5], [2, 1, 6]])\n", "B = np.array([[1, 2, 3], [4, 5, 6]])\n", "print('Matrix B has this shape: {}'.format(B.shape))\n", "\n", "B_transpose = B.T\n", "print('B transposed is:\\n{}'.format(B_transpose))\n", "\n", "# Stack the operations in one line!\n", "print('B transposed has shape: {}'.format(B.T.shape))\n", "```\n", "\n", "We don't need to create an intermediate variable ``B_transpose`` to hold the transposed value. \n", "\n", "Finally, to calculate $\\mathbf{AB}^T$ we can simply write: ``np.dot(A, B.T)`` or even shorter ``A @ B.T``. Complete the code:\n", "```python\n", "# (out loud you would say: \"A times B transposed\")\n", "print('A times B transposed is {}'.format( ... ))\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The matrix inverse\n", "\n", "A square matrix (a matrix with the same number of rows and columns) can be ***inverted***. If the matrix is $\\mathbf{A}$, then the inverse is written as $\\mathbf{A}^{-1}$. Inversion is related to the concept of division. We won't go into the mathematical details of how an inverse is calculated though.\n", "\n", "The command for inverting is:\n", "```python\n", "A = np.random.random((4, 4))\n", "A_inv = np.linalg.inv(A)\n", "```\n", "\n", "Other linear algebra commands of interest in the ``linalg`` section of NumPy are:\n", "* ``np.linalg.eig`` for eigenvalues and eigenvectors\n", "* ``np.linalg.eigvals`` for the eigenvalues only\n", "* ``np.linalg.svd``for the singular value decomposition (related to PCA)\n", "* ``np.linalg.solve``for solving linear matrix equations\n", "* ``np.linalg.solve``for solving linear matrix equations\n", "\n", "Type ``help(np.linalg.info)`` for more details, and other linear algebra routines." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Check your knowledge\n", "\n", "You can apply most concepts learned above by solving this problem.\n", "\n", "If you build a least squares regression model with [more than one variable](https://learnche.org/pid/least-squares-modelling/multiple-linear-regression) you can assemble all you input data (predictors, columns used to make the prediction) in a single matrix $\\mathbf{X}$. For example: \n", "\n", "$$ y = b_0 + b_\\text{P} x_\\text{P} + b_\\text{T} x_\\text{T}$$\n", "\n", "has 2 prediction columns: $x_\\text{P}$ and $x_\\text{T}$; but there is also an intercept.\n", "\n", "Assemble matrix $\\mathbf{X}$ to have 3 columns: a column of ones, a column for $x_\\text{P}$ and one for the values of $x_\\text{T}$. The following data are from a designed experiment, where the 3 required columns are:\n", "\n", "$$ \\mathbf{X} = \n", "\\begin{bmatrix}\n", "1 & -1 & -1 \\\\\n", "1 & +1 & -1 \\\\\n", "1 & -1 & +1 \\\\\n", "1 & +1 & +1 \\\\\n", "1 & 0 & 0\\\\\n", "1 & 0 & 0 \n", "\\end{bmatrix}\n", "$$\n", "\n", "***Hint***: There are many way to create the above matrix, but consider using the ``np.hstack`` function to stack vectors next to each other (side-by-side).\n", "\n", "To build a regression model you also need your vector of $\\mathbf{y}$ values, the measured values:\n", "$$ \\mathbf{y} = \n", "\\begin{bmatrix}\n", "55\\\\90 \\\\ 190 \\\\ 170 \\\\135 \\\\ 142\n", "\\end{bmatrix}\n", "$$\n", "\n", "The regression coefficients, the 3 values of $\\mathbf{b} = \\left[ b_0, \\, b_\\text{P},\\, b_\\text{T}\\right]$ can be found from the following [linear algebra expression](https://learnche.org/pid/least-squares-modelling/multiple-linear-regression#estimating-the-model-parameters-via-optimization): $\\mathbf{b} = (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T \\mathbf{y}$\n", "* matrix multiplication\n", "* matrix transposes, and\n", "* matrix inversion.\n", "\n", "> 1. Calculate that vector $\\mathbf{b}$ with 3 values. A good check that you got the correct values is if the first entry is equal to ``np.average(y)``\n", "> 2. Once you have the coefficients, you can make predictions from the regression model: $\\hat{\\mathbf{y}} = \\mathbf{X} \\mathbf{b}$. Check the dimensions of $\\mathbf{X}$ and $\\mathbf{b}$ to ensure they are correct.\n", "> 3. Lastly calculate the residuals, the prediction error = actual $\\mathbf{y}$ minus predicted = $\\mathbf{e} = \\mathbf{y} - \\hat{\\mathbf{y}}$.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Wrap up this section by committing all your work. Have you used a good commit message? Push your work, to refer to later, but also as a backup." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">***Feedback and comments about this worksheet?***\n", "> Please provide any anonymous [comments, feedback and tips](https://docs.google.com/forms/d/1Fpo0q7uGLcM6xcLRyp4qw1mZ0_igSUEnJV6ZGbpG4C4/edit)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# IGNORE this. Execute this cell to load the notebook's style sheet.\n", "from IPython.core.display import HTML\n", "css_file = './images/style.css'\n", "HTML(open(css_file, \"r\").read())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "222px" }, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }