05_numerical_data_in_numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical data in `numpy`\n",
    "\n",
    "[`numpy`](https://docs.scipy.org/doc/numpy/) is a very powerful library for working with numerical data in Python. It introduces the __Array__ data structure, which can contain multi-dimensional numerical data. The `numpy` library is __not__ part of the Python standard library. However, it comes bundles with [Anaconda](01_anaconda.ipynb) so you should already have it installed. The usual way to import `numpy` is as follows:\n",
    "\n",
    "    import numpy as np\n",
    "    \n",
    "This gives us access to all the `numpy`-functions using the prefix `np`. This is a convention, and you should do the same in your own code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.set_printoptions(precision=3, linewidth=75, edgeitems=1)\n",
    "np.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple way to instantiate an array is by giving it a list or another iterable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_in_list = [1, 2, 5, 3]\n",
    "data_as_np = np.array(data_in_list)\n",
    "data_as_np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1 / data_as_np, data_as_np**3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading data with `numpy`\n",
    "\n",
    "The `numpy` library comes with a few functions that can read numerical data from text files. The `np.loadtxt`-function is a fast reader with some basic functionality for skipping header lines, converting values to floats, etc. For more sophisticated files the `np.genfromtxt`-function can be used. It also supports handling missing values, but is slower.\n",
    "\n",
    "Note that these functions take a filename as input, so that you do not need to use the built-in `open`-function to open a file beforehand. Assume we have a text-file containing data of the following form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat data/numpy_simple.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These data can be loaded with a very simple `np.loadtxt`-command as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt('data/numpy_simple.txt')\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arrays in `numpy` have a shape specifying how big the dataset is. In this case we have 6 rows and 8 columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing and vectorization\n",
    "\n",
    "Similarly to lists and other sequences in Python, `numpy` arrays can be indexed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[0]    # First row (= row 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[2:5]    # Rows 2, 3 and 4 (3rd, 4th and 5th)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, with multi-dimensional arrays (in this case 2-dimensional), we can specify each dimension in the index separated by commas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[4, 5]    # Element in row 4, column 5 (5th row, 6th column)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[:, :2]    # All rows, first two columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `numpy`-arrays most operations are vectorized. That means that we do not need to explicitly loop over the elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[:, 4] + data[:, 6]    # Add columns 4 and 6 together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.exp(data[:, -1])    # Exponentiate the last (-1) column"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`numpy` also comes with summary functions like `sum`, `mean`, `std`, `var`, `median`, etc that can operate on a whole array or a given dimension (`axis`) of the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(data)    # Calculate the mean of the whole array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(data, axis=0)    # Calculate the mean of each column (along the rows, the 0th axis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More advanced reading of data\n",
    "\n",
    "Most datafiles are not as clean as the simple datafile we have been working with above. Let us instead try to load the following file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat data/numpy_header.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A naive use of `np.loadtxt` will fail because `numpy` tries to interpret the header as data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2 = np.loadtxt('data/numpy_header.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead we must give the `np.loadtxt` some more information. To get some help about a function and which parameters it takes, you can write a question mark after its name,\n",
    "\n",
    "    np.loadtxt?\n",
    "    \n",
    "or press `<shift>` and `<tab>` inside the paranthesis. Pressing `<shift>` and `<tab>` twice will give even more information.\n",
    "\n",
    "In this case, we notice that there is an argument called `skiprows` that can be used to ignore the header."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2 = np.loadtxt('data/numpy_header.txt', skiprows=13)\n",
    "np.allclose(data, data2)    # Test if data and data2 contains the same elements (within a tolerance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking more closely at the data, we also notice that there is one datapoint with the value of -999 that probably designates a missing data point. We can convert this to a `nan`-value to handle it properly as we are reading the data. However, to do so we need to use the more sophisticated `np.genfromtxt`-function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data3 = np.genfromtxt('data/numpy_header.txt', skip_header=13,\n",
    "                      missing_values='-999', usemask=True).filled(np.nan)\n",
    "data3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually handling the missing data is a two-step process. First we create a masked array, where the mask denotes which data are missing. This can be seen if we look directly at the output from `np.genfromtxt` (note the single `True`-value in the mask)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.genfromtxt('data/numpy_simple.txt', missing_values='-999', usemask=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the [`numpy` docs](https://docs.scipy.org/doc/numpy-dev/user/basics.io.genfromtxt.html) for more information about reading data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Back to our real world example\n",
    "\n",
    "Let us pick up the example we started in the previous lecture. We are able to read in data in a dirty format using the following code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from collections import OrderedDict\n",
    "\n",
    "lines = list()\n",
    "data = OrderedDict()\n",
    "\n",
    "def strip_field(field):\n",
    "    key, plus, value = field.partition('+')\n",
    "    return value.lstrip('0')\n",
    "\n",
    "with open('data/dirty_data.txt', mode='r') as fid:\n",
    "    next(fid)   # Throw away header line\n",
    "    for line in fid:\n",
    "        fields = [strip_field(f) for f in line.strip().split()]\n",
    "        lines.append(fields)\n",
    "        \n",
    "for line in lines:\n",
    "    name = line[0]   # First column is the identifier\n",
    "    values = [float(v) for v in line[1:-1]]   # Convert to float, ignoring the first and last column\n",
    "    data.setdefault(name, list()).append(values)   # setdefault creates the list if it does not already exist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converting list of lists to `numpy.array`s\n",
    "\n",
    "In order to make these data easier to work with, we'll now convert each list of lists to a `numpy.array`. We do this by passing each list of lists to a new `OrderedDict` as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data_np = OrderedDict((k, np.array(v)) for k, v in data.items())\n",
    "data_np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating the mean of the first column of the `CP4`-data is now trivial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_np['CP4'][:, 0].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The real world again ...\n",
    "\n",
    "Of course, we are not interested in a simple mean of these data ...\n",
    "\n",
    "For each identifier we have 10 measurements. Let us focus on the first column. We notice that every second measurement is approximately 200 away from the previous one. This is by design of the experiment, and an effect we wish to remove. Also, the first two measurements are part of the calibration and should not be included in the final mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "means = OrderedDict()\n",
    "\n",
    "for name, values in data_np.items():\n",
    "    angles = values[2:, 0].copy()  # First column, skip first two measurements\n",
    "    angles[1::2] += 200            # Add 200 to every second measurement\n",
    "    angles = np.mod(angles, 400)   # Modulo 400 to keep in interval (0, 400)\n",
    "    means[name] = angles.mean()\n",
    "    \n",
    "means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Views and its scary consequences ...\n",
    "\n",
    "One of the great advantages with `numpy` (in addition to making operations on big data sets convenient) is its speed. When doing operations on whole arrays at once, the heavy lifting is done inside C-code giving quite fast calculations.\n",
    "\n",
    "One important thing to be aware of is that `numpy`-arrays and views of `numpy`-arrays (sliced sub-arrays) are passed around as references. This means that the code is efficient both in terms of speed and memory, but can also bring some surprises.\n",
    "\n",
    "Note the `.copy()` done in the line `angles = values[2:, 0].copy()` above. This is not necessary, and actually gives us a small speed penalty. However, lets explore what happens without the copy ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = data_np['CP5']\n",
    "values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angles = values[2:, 0]\n",
    "angles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, so good ... Now we'll add 200 to every second element of `angles` ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angles[1::2] += 200\n",
    "angles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case (without the `.copy()`) `angles` is a view to the data in `values`. That means that it references the exact same data. So now that we have added 200 to `angles` we have done the same to `values`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And values was a reference to the data inside the `data_np`-dictionary. So those data are also changed (they are the __same__ data!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data_np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our case this would not have had any immediate consequences. We are after all not using the original data again. However, if we were to run the calculations again, or do some new calculations on the original data we would have gotten weird head-scratching bugs to deal with ...\n",
    "\n",
    "Note that this is a __feature__ of Python and `numpy` (so no reason to rush over to [bugs.python.org/](https://bugs.python.org/)). But it is important to be aware of, and to be conscious about changing data referenced by several names."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with matrices\n",
    "\n",
    "`numpy` data can be thought of both as (possibly multidimensional) arrays of numbers, and as vectors and matrices. So far we have worked with `numpy` data as arrays of numbers, where `numpy` provides an efficient way of handling many numbers at once.\n",
    "\n",
    "However, if we think of our data as vectors and matrices, we'll have additional needs such as finding lengths of vectors or doing dot products and matrix multiplication.\n",
    "\n",
    "---\n",
    "\n",
    "If you are have experience working with Matlab, `numpy` provides all the functionality you are used to. However there are some subtle differences to be aware of that might trip you up. These include:\n",
    "\n",
    "+ In Matlab `*` means matrix multiplication, while the special notation `.*` is used for elementwise multiplication. In Python it is essentially the opposite: `*` is used for elementwise multiplication, while the special notation `@` is used for matrix multiplication (the `@` notation was introduced in Python 3.5, in earlier versions you must use the function `.dot`).\n",
    "\n",
    "+ In Matlab, 1-dimensional vectors have an orientation. That is, they are either row or column vectors. In Python, that is not the case. 1-dimensional vectors are purely 1-dimensional (or look a lot like row vectors). The `@` still mainly does the right thing. However stacking many vectors together or doing an outer product of vectors is a bit cumbersome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([1, 2, 3])\n",
    "x, x.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.T @ x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above calculates the square of the norm of the vector, $\\| \\vec x \\|^2 = x^T x$. This is also available in the `np.linalg`-package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.norm(x)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now say we want to calculate the outer product $x x^T$ instead of the dot product $x^T x$. The naive way will not work (and might lead to some head scratching):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x @ x.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead we need to use the `np.outer` function to calculate $x x^T$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.outer(x, x)\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix multiplication is easily done with the `@` operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A @ A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While matrix inverse and other matrix operations are in the `numpy.linalg` package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = np.random.randint(1, 7, (2, 2))\n",
    "B_inv = np.linalg.inv(B)\n",
    "print(B)\n",
    "print(B_inv)\n",
    "B @ B_inv"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}