4-Memory-Layout-and-Views.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Memory Layout and Views\r\n",
    "\r\n",
    "This last section deals with a more theoretical topic that you tend to not encounter directly when approaching numpy as a user. However, it is important to be aware of it, because it can have significant performance implications.\r\n",
    "\r\n",
    "## Memory Layout\r\n",
    "\r\n",
    "The first topic is a peak under the hood of modern HPC array libraries; numpy in particular. How are arrays stored in your memory, and why does this make them so fast? Internally, every array is a C-style byte array with some metadata associated with it, i.e., - if you speak some C - `malloc(some_size)`. This means that, at the lowest level, every array is 1 dimensional and has no notion of data type, shape, or number of dimensions. \r\n",
    "\r\n",
    "To bring structure into this buffer, we introduce the notion of an `itemsize`, which is simply a measure of how many bytes make up a single item in the array, e.g., a float32 is made up of 32 bit (hence the name) or 4 bytes, so it has `itemsize=4`. Next we introduce a `stride`, which you can think of as a step size. It tells numpy how many bytes it has to move from its current position to find the next element (along an axis) in the array. \r\n",
    "\r\n",
    "Notice, that `stride` and `itemsize` are different variables, which is very useful for optimization. Setting `stride = itemsize` gives us a standard 1d array that we are familiar with. To get the location of the next element from our block of bytes, numpy moves `strides`(`=itemsize`) many bytes forward. If we decide that we actually only need every second element in our array, we can achieve this very easily by setting `strides=2*itemsize` and without copying any data (!). This is huge compared to the alternative of having to copy the array, because memory allocation and copy operations are dreadfully slow compared to the speeds at which a CPU operates. In general, we want to avoid copying data as much as possible, as less copies means more speed (usually).\r\n",
    "\r\n",
    "Next, we can add `offset` and `size` as metadata to allow trunciating of the array. This follows the same motivation to avoid copying data as `strides` does: `offset` allows choosing a different starting point of the array, and `size` allows choosing a different end point. Combining this with `strides`, we can now remodel the array according to any slice (`3:7:2`) by simply setting `offset`, `size`, and `stride` appropriately; no copies needed. We only modify 3 numbers, regardless of size, which obviously beats making a copy for large arrays. We can even flip the array by using a negative `stride` (and setting `offset` appropriately).\r\n",
    "\r\n",
    "To get support for multi-dimensional arrays we simply allow multiple `strides` and add a `shape` parameter. Recall that a `stride` indicates how many bytes we have to move to find the next element along the given axis. If our array has 2 dimensions (a matrix) and a row contains 5 elements at every `row_stride` byte, then we can get to the next column by making 5 `row_stride` steps, i.e. `column_stride = 5 * row_stride`. We can play this game in any number of dimensions, and even switch between them without copying, as long as we can assign a constant `stride` to each axis in the array as well as a valid `shape` so that we never exceed the underlying buffer. This is why `ndarray.reshape` or `ndarray.transpose` are so fast. They don't need to create a new array; all that is needed is an update of the array's metadata.\r\n",
    "\r\n",
    "Going into numpy itself, the parameters (`strides`, `itemsize`, `shape`, etc.) are available as read-only values should you need them.\r\n",
    "\r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\r\n",
    "\r\n",
    "some_array = np.arange(100, dtype=np.float32)\r\n",
    "\r\n",
    "print(f\"The size of a single element: {some_array.itemsize}\")\r\n",
    "print(f\"The strides of each dimension: {some_array.strides}\")\r\n",
    "print(f\"Number of elements along each axes: {some_array.shape}\")\r\n",
    "print(f\"Total number of elements: {some_array.size}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the take-away from all of this? The main message here is that, if you are interested in speed, you should avoid copying data where possible, and that numpy has an array of tricks (see what I did there?) up its sleve to help you avoid those copies.\r\n",
    "\r\n",
    "One more advanced usage of this, which is - unfortunately - out of scope here, is to use a technique called stride tricks. When doing this, you manually change the `strides` and `shape` of your array to avoid a copy and speed up your code. If you recall the rolling average filter that you were asked to implement in the previous section, a faster alternative to the implemtation shown there would be to use stride tricks. Instead of constructing `vector_batch` using fancy indexing, we can modify the `strides` and `shape` of `vectorA`. We can let the first dimension advance by `stride=itemsize` (just like the original array), but we let it stop at the last element for which we can compute the filter. To this we add a new dimension that also advances by `sride=itemsize` and has a maximum length equal to the filter length. We obtain the same `vector_batch` as we did in the exercise, but we avoid copying data. This has little effect for the exercise itself; however when we increase the amount of data points into the billions such optimizations can have a massive impact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Views\r\n",
    "\r\n",
    "The final topic of this introduction to numpy are `view`s. Views are a direct manifestation of all the theory presented above and the idea that arrays should be able to share data. You can get a view explicitly by calling `ndarray.view()` on any array, or indirectly as the result of tuple-indexing (with slices) or from `ndarray.reshape()`. Essentially, a view is a copy of the metadata of the original array to which modications may have been applied, but that shares the same data as the original array. Where possible, numpy will return such a view (it is faster) and - as a consequence of sharing data - writing to the view will modify the original array and vice versa. This is the cause of a lot of confusion among beginners."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "foo = np.arange(9)\r\n",
    "\r\n",
    "# bar will share data with foo\r\n",
    "bar = foo.reshape(3,3)\r\n",
    "\r\n",
    "print(f\"bar before modification of foo:\\n {bar}\")\r\n",
    "foo[4] = 42\r\n",
    "print(f\"bar after modification of foo:\\n {bar}\")\r\n",
    "\r\n",
    "print(f\"foo before modification of bar:\\n {foo}\")\r\n",
    "bar[2, 2] = 43\r\n",
    "print(f\"foo after modification of bar:\\n {foo}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Exercise\r\n",
    "\r\n",
    "What is the value of `result` in the following program?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "foo = np.ones(9)\r\n",
    "bar = foo.reshape(3,3)\r\n",
    "bar += np.array([1, 2, 3])[:, None]\r\n",
    "baz = bar.ravel()\r\n",
    "baz[7] = 2\r\n",
    "\r\n",
    "result = np.sum(foo[-3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\r\n",
    "<summary>Answer (click me to reveal)</summary>\r\n",
    "\r\n",
    "The answer is 10. `foo.reshape` produces a view into `foo`, and all modifications to bar will also affect foo. Similarly, baz is a view into bar, which again shares data with foo. After the dust has settled from the modifications, foo looks like this `[2., 2., 2., 3., 3., 3., 4., 2., 4.]` and summing the last three elements is 4 + 2 + 4 = 10.\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
}