3-Broadcasting-and-Fancy-Indexing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Broadcasting and Fancy Indexing\r\n",
    "\r\n",
    "In the previous two sections we have looked at features that are fairly common for numeric libraries. You will find similar patterns in other libraries from diffent languages, e.g. in Eigen. Now, we will take a look at two advanced features of numpy: broadcasting and fancy indexing.\r\n",
    "\r\n",
    "\r\n",
    "## Broadcasting\r\n",
    "\r\n",
    "Array broadcasting allows numpy to perform operations on arrays that would otherwise not be compatible with each other. It allows numpy to reapeat (broadcast) an array along dimensions with length 1 in order to match the shape of another array. This works in tandem with the tuple-indexing variable `np.newaxis` (`None`) we have encountered in the first section. The syntax `some_vector[None, :]` allows us to insert a new dimension along which the array can be broadcasted/repeated. This may sound confusing, but will become quite intuitive after having seen a few examples, which you can find below. Once you start working with it for a bit, it will become second nature.\r\n",
    "\r\n",
    "**Note**: Two (or more) arrays that can have their shapes matched via broadcasting are called broadcastable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\r\n",
    "\r\n",
    "# Examples of boradcasting\r\n",
    "foo =  np.arange(4)\r\n",
    "bar = np.zeros((4, 4))\r\n",
    "\r\n",
    "broadcasted = bar + foo[None, :]\r\n",
    "print(f\"foo was repeated and added to each row of bar:\\n {broadcasted}\")\r\n",
    "\r\n",
    "# np.newaxis is also valid, but less common\r\n",
    "broadcasted = foo[:, np.newaxis] + bar\r\n",
    "print(f\"foo was repeated and added to each column of bar:\\n {broadcasted}\")\r\n",
    "\r\n",
    "bar = np.array([[2, 4, 5], [8, 1, 3]])\r\n",
    "foo = np.array([3, 2, 1])\r\n",
    "\r\n",
    "# same as foo[None, :], but more verbose\r\n",
    "broadcasted = np.maximum(bar, foo.reshape(-1, *foo.shape))\r\n",
    "print(f\"Element-wise maximum between foo and bar along each row of bar:\\n {broadcasted}\")\r\n",
    "\r\n",
    "bar = np.arange(25).reshape(5,5)\r\n",
    "bar_mean = np.mean(bar, axis=1)\r\n",
    "bar = bar - bar_mean[:, None]\r\n",
    "print(f\"Bar with row-mean removed:\\n {bar}\")\r\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\r\n",
    "\r\n",
    "1. Compute the values of $f(x,y) = \\sin(x)^{10} + \\cos(10+xy)\\cos(x)$ for each point of an evenly spaced grid. The grid should range from (0,0) to (5, 5) and have 50 points along each axis.\r\n",
    "2. `point_batch` contains a batch of points in 3d space stacked along the first axis. Retrieve the point that is closest to `estimated_point`.\r\n",
    "3. Compute the kronecker product of `vectorA` and `vectorB`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point_batch = np.array([\r\n",
    "    [0.42790781, 0.22270672, 0.38540413],\r\n",
    "    [0.92987344, 0.30558244, 0.33538531],\r\n",
    "    [0.17709567, 0.58556512, 0.03374166],\r\n",
    "    [0.32993173, 0.1081541 , 0.9863263 ],\r\n",
    "    [0.59581806, 0.14109558, 0.93316973],\r\n",
    "    [0.83190349, 0.94947723, 0.16788122],\r\n",
    "    [0.54719748, 0.4821661 , 0.61150104],\r\n",
    "    [0.08215392, 0.26003697, 0.1960123 ],\r\n",
    "    [0.01890119, 0.98642855, 0.21429429],\r\n",
    "    [0.36139499, 0.87299214, 0.06041583]\r\n",
    "])\r\n",
    "estimated_point = np.array([0.73566544, 0.6376392 , 0.81113481])\r\n",
    "\r\n",
    "vectorA = np.arange(3)\r\n",
    "vectorB = np.arange(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\r\n",
    "<summary>Answers (click me to reveal)</summary>\r\n",
    "\r\n",
    "Exercise 1.1\r\n",
    "```\r\n",
    "x = np.linspace(0, 5, 50)\r\n",
    "y = np.linspace(0, 5, 50)[:, None]\r\n",
    "z = np.sin(x) ** 10 + np.cos(10 + y * x) * np.cos(x)\r\n",
    "```\r\n",
    "\r\n",
    "Exercise 1.2\r\n",
    "```\r\n",
    "distances = np.sqrt(np.sum(np.squared(point_batch - estimated_point[None, :]), axis=1))\r\n",
    "closest_idx = np.min(distances)\r\n",
    "nearest_point = point_batch[closest_idx]\r\n",
    "```\r\n",
    "\r\n",
    "Exercise 1.3\r\n",
    "```\r\n",
    "kron_product = vectorA[:, None] * vectorB[None, :]\r\n",
    "```\r\n",
    "\r\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fancy Indexing\r\n",
    "\r\n",
    "Fancy indexing, nowadays known as advanced indexing, expands the traditional tuple indexing mechanism by allowing you to use integer and boolean arrays as indexes. \r\n",
    "\r\n",
    "Let's start with integer arrays. In its simplest form, you can think of this type of fancy indexing as a generalization of slicing. The slice `3:9:2` will select the elements `[3, 5, 7]` along the axis where it is used. Instead of using the slice `3:9:2`, e.g. `some_array[3:9:2]`, you can use the equivalent array instead `some_array[[3, 5, 7]]` and the result will be the same. The array generalizes the slice, because it does not have to be ordered and nothing stops the array from having repeated elements. Hence, `some_array[[3, 5, 7, 5, 3]` is also valid, and will select elements along the given axis in the order they were listed.\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "some_array = 2 * np.arange(10)\r\n",
    "\r\n",
    "print(f\"Some Array: {some_array}\")\r\n",
    "print(f\"Traditional Slice: {some_array[3:9:2]}\")\r\n",
    "print(f\"Advanced/Fancy Indexing: {some_array[[3, 5, 7]]}\")\r\n",
    "print(f\"Repeated Elements: {some_array[[3, 5, 5, 7]]}\")\r\n",
    "print(f\"Reordering Elements: {some_array[[3, 7, 5, 8, 1, 2]]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The array generalizes the slice further by allowing higher-dimensional arrays to be used. This uses the same selection mechanism along the indexed dimension as the 1d case, but changes the shape of the resulting array by replacing the indexed dimension of the original array with the dimensions of the indexing array. For example, if `index_array.shape == (2,2)` and `some_array.shape == (3, 3)`, we can index `some_array` via `some_array[index_array]` and then `some_array[index_array].shape == (2, 2, 3)`, i.e., `some_array`'s first dimension was replaced by `index_array.shape`.\r\n",
    "\r\n",
    "This being said, multi-dimensional fancy indexing (apart from the broadcasting case explained below) is rare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "some_array = some_array = 2 * np.arange(10)\r\n",
    "index_array = np.array([\r\n",
    "    [1, 2],\r\n",
    "    [2, 1]\r\n",
    "])\r\n",
    "\r\n",
    "print(f\"Some Array: {some_array}\")\r\n",
    "print(f\"Indexing using a 2d-array:\\n {some_array[index_array]}\")\r\n",
    "\r\n",
    "some_array = np.arange(9).reshape(3, 3)\r\n",
    "index_array = np.array([\r\n",
    "    [1, 2],\r\n",
    "    [2, 1]\r\n",
    "])\r\n",
    "print(\"-------\")\r\n",
    "print(f\"some_array.shape: {some_array.shape}\")\r\n",
    "print(f\"index_array.shape: {index_array.shape}\")\r\n",
    "print(f\"result.shape: {some_array[index_array].shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If `some_array` has more than one dimension, we can use fancy indexing along more than one dimension of it at the same time. In this case, the arrays used to index each dimension must match in shape or must be broadcastable. The result will be an array with elements coming from `some_array` at the positions described by the (broadcasted) arrays. An example will make this clear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "some_array = np.arange(9).reshape(3, 3)\r\n",
    "\r\n",
    "print(f\"The original array:\\n {some_array}\")\r\n",
    "print(f\"Traditional Slicing:\\n {some_array[:2, :2]}\")\r\n",
    "print(f\"Advanced/Fancy indexing:\\n {some_array[[[0], [1]], [0, 1]]}\")\r\n",
    "\r\n",
    "index = np.array([0, 2])\r\n",
    "print(f\"Corners via broadcasting:\\n {some_array[index[:, None], index[None, :]]}\")\r\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fancy indexing can also be combined with traditional tuple-indexing. This scenario is a bit more involved, and I haven't encountered a real-world use-case for it yet. Nontheless, it is a very cool feature, and you can read more about it in the official documentation: [https://numpy.org/doc/stable/reference/arrays.indexing.html#combining-advanced-and-basic-indexing](https://numpy.org/doc/stable/reference/arrays.indexing.html#combining-advanced-and-basic-indexing).\r\n",
    "\r\n",
    "Apart from integer arrays, fancy indexing also allows you to use boolean arrays to index an array. The best way to think of boolean array indexing is to think of it as a mask that selects individual elements along the chosen axis. This allows us to use logical expressions to index arrays, such as `some_array[some_array < 5]` to retrieve all elements of the array with value less than 5. Such expressions are, indeed, their most common use-case.\r\n",
    "\r\n",
    "Their integer array equivalent is `np.nonzero(bool_array)` and the same rules as above apply. One quirk of boolean array indexing occurs when using boolean indexing in combination with fancy indexing along other axes: the resulting array is _allways_ one dimensional and broadcasting to higher dimensions will fail. There is, again, non-trivial interaction with traditional tuple indexing, but this is rare and out of scope here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "some_array = np.arange(5)\r\n",
    "\r\n",
    "print(f\"The original array:\\n {some_array}\")\r\n",
    "print(f\"Boolean indexing:\\n {some_array[[True, False, True, False, True]]}\")\r\n",
    "print(f\"Indexing via expressions:\\n {some_array[some_array <= 3]}\")\r\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\r\n",
    "\r\n",
    "1. Sum all elements of `vectorA` except those with the highest value.\r\n",
    "2. Create a new matrix from `matrixA` that dublicates those rows of `matrixA` that only contain uneven elements. Remove all rows that contain even elements.\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorA = np.array([1, 3, 4, 5, 3, 3, 4, 5, 5, 3, 2, 1])\r\n",
    "matrixA = np.array([\r\n",
    "    [3, 5, 6],\r\n",
    "    [1, 3, 5],\r\n",
    "    [2, 6, 2],\r\n",
    "    [7, 11, 43],\r\n",
    "    [19, 16, 3]\r\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\r\n",
    "<summary>Answers (click me to reveal)</summary>\r\n",
    "\r\n",
    "Exercise 2.1\r\n",
    "```\r\n",
    "np.sum(vectorA[vectorA != np.max(vectorA)])  # == 24\r\n",
    "```\r\n",
    "\r\n",
    "Exercise 2.2\r\n",
    "```\r\n",
    "matrixA[[1, 1, 3, 3]]\r\n",
    "```\r\n",
    "\r\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Exercise\r\n",
    "\r\n",
    "`vectorA` contains position measurements of a person walking down a (1D) corridor. Samples are taking at a rate of 100Hz and the vector stores 10 seconds worth data. Unfortunatley, the data is quite noisy (accuracy ~ 10 cm). To improve accuracy you would like to use your a priori knowledge that people tend to not teleport and apply a moving average filter. You decide that 0.1s is a reasonable window size for this scenario. \r\n",
    "\r\n",
    "In order to easily apply this filter in a vectorized fashion, you need to reorganize the data stored in `vectorA`. Create a batch of vectors (`vector_batch`, a matrix) that, in each row, stores 0.1 seconds worth of position data, and where the column corresponds to progressing time. That is, the first row of the matrix contains the measurements from 0s to 0.1s, the second row contains the measurements from 0.01s to 0.11s, the third row contains the measurements from 0.02s to 0.62s, and so on until the last row, which contains measurements from 9.9s to 10s. (Hint: Use broadcasting to create a matrix containing the required indices and then use fancy indexing to retrieve the actual measurements.)\r\n",
    "\r\n",
    "After you have created `vector_batch`, apply the rolling average filter by computing the average value of each row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorA = np.load(\"position_data.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\r\n",
    "<summary>Answer (click me to reveal)</summary>\r\n",
    "\r\n",
    "```\r\n",
    "# at 100Hz 50 samples make up 0.5s\r\n",
    "window = np.arange(10)\r\n",
    "\r\n",
    "# a new window starts at every sample that has at least len(window) samples following it\r\n",
    "row_idx = np.arange(len(vectorA) - len(window))  \r\n",
    "\r\n",
    "# matrix of sample indices via broadcasting\r\n",
    "batch_idxs = row_idx[:, None] + window[None, :]\r\n",
    "\r\n",
    "# fancy-indexing to retrieve the values from the indices\r\n",
    "vector_batch = vectorA[batch_idxs]\r\n",
    "\r\n",
    "# filtering\r\n",
    "rolling_average = np.mean(vector_batch, axis=-1)\r\n",
    "\r\n",
    "# Bonus: new accuracy for this scenario\r\n",
    "vectorA_true = np.linspace(0, 5, 5*100)\r\n",
    "absolute_error = np.abs(vectorA_true[len(window):] - rolling_average)\r\n",
    "print(f\"The accuracy of the new sequence is ~ {np.mean(absolute_error)*100:.2f} cm\")\r\n",
    "```\r\n",
    "\r\n",
    "</details>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.8 64-bit",
   "name": "python388jvsc74a57bd05d69d85edf19066ea65ec45e9305decd722722c415050988d8ca48a84314af7c"
  },
  "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.8"
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}