numpy_ndarrays.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4bf3d82c-4562-4c75-9196-0c8e352261f6",
   "metadata": {},
   "source": [
    "- Archer numpy example for the same problem along with profiler output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "608dc987-845d-4132-87c7-66136590917b",
   "metadata": {},
   "source": [
    "# The `numpy` library for efficient numeric computing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cbe6b42-7e9e-449d-98ab-6a6758ddd7dd",
   "metadata": {},
   "source": [
    "## Performance of lists for aggregate numeric data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf05c57b-34c1-40b8-be7f-6b0a78c53b02",
   "metadata": {},
   "source": [
    "Often it is essential to hold a set of related numeric values in a variable e.g. a dataset of voltage, current and temperature measurements from a scientific experiment. Although, python provides built-in compound data structures like lists (for mutable data), tuples (for immutable data) and sets (for unique, non-repeating data), these are general purpose containers that have not been specifically optimised for numeric computing.\n",
    "\n",
    "A list can contain elements of different types of data. For example, an entity of any data type can be a member a list. To obtain this flexibility of being a general purpose container, the underlying implementation trades-off various aspects of performance (e.g. list elements may not always be allocated contiguous memory chunks). A python library called `numpy` provides a specialised data container known as `ndarray` that is capable of handling the performance demands of high performance computing, and forms the basis of the modern scientific computing stack."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3e95a66-0fbc-466c-995d-78e01d373448",
   "metadata": {},
   "source": [
    "## Performance of loops on aggregate data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1136802-6554-490b-9e0e-8eb343605d67",
   "metadata": {},
   "source": [
    "In scientific computing, it is a common requirement to handle large data sets and perform mathematical operations on them. The traditional method of accessing elements in a dataset by looping through them individually (to perform element-by-element computations) does not scale well in terms of performance for large data vectors.\n",
    "\n",
    "In numeric computing, we care about performance of (nested) loops reading/writing to/from arrays of floats/doubles.\n",
    "\n",
    "- Often performing same operation on all elements (individually, pairwise, or across >2 arrays)\n",
    "- Explicit loops in Python have very significant overhead\n",
    "- Type checking is done for each element (even though all elements of an array of same type)\n",
    "- Array operations are not vectorised by CPython (the reference implementation)\n",
    "\n",
    "### SIMD vectorisation\n",
    "One way to facilitate such large-scale computations is to take advantage of a bank of floating point units on modern processors known as [SIMD units](https://en.wikipedia.org/wiki/Single_instruction,_multiple_data) as illustrated below.\n",
    "\n",
    "![simd image](./images/simd_addition.png)\n",
    "\n",
    "\n",
    "While such a construct is not available natively in python for numeric computations, the `numpy` library provides one."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce7d8ff5-3ec1-4a39-af0f-1c02e6c33aa8",
   "metadata": {},
   "source": [
    "## Using the `numpy` library for numeric computing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "401e6525-341d-4eb4-a2f5-1228fb0e43e9",
   "metadata": {},
   "source": [
    "[`numpy`](https://numpy.org/) is a third party library usable within python projects that implements efficient data containers for handling numeric computing requirements. Numpy defines a fast array data type that is scalable to arbitrary dimensions. It provides a huge amount of operations on this array type including all common linear algebra operations. For instance, the ability to '_broadcast_' array operations is particularly useful to perform numeric operations on non-identical (but still meldable to be compatible) array shapes.\n",
    "\n",
    "![numpy_broadcasting image](./images/numpy_broadcasting.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4f46385-324a-48dd-b3c4-b0f7364b6f9e",
   "metadata": {},
   "source": [
    "NumPy (\"Numeric Python\") provides:\n",
    "\n",
    " - Statically typed, fixed-size multidimensional arrays (matrices) - type `numpy.ndarray`\n",
    " - Like C arrays wrapped in Python objects\n",
    " - Linear algebra operations, random number generation\n",
    " - Fourier transforms, polynomials etc\n",
    " - Fast precompiled functions for numerical routines\n",
    " - Efficient storage (contiguously in memory) and execution\n",
    "\n",
    "https://docs.scipy.org/doc/numpy/user/whatisnumpy.html\n",
    "\n",
    "Numpy can perform vectorised operations on all elements of arrays with a single function call"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e419e126-b0e9-408d-9156-24d80290b64b",
   "metadata": {},
   "source": [
    "The `requirements.txt` provided as part of setting up an isolated python environment for this workshop already includes the `numpy` library, and therefore it can be readily imported for use in our projects. The recommended method of importing `numpy` functionality within our project's python codes is as follows:\n",
    "\n",
    "```python\n",
    "import numpy as np\n",
    "```\n",
    "\n",
    "This allows us to access the methods of the `numpy` library using the `np` namespace. This prevents clashing of several popular function names implemented within numpy from similarly named entities from the standard library or other third party libraries. There are several features provided by numpy such as slicing and provision of shared access data views held in numpy's `ndarrays`. For a basic overview of using `numpy`, please refer to the following resources:\n",
    "\n",
    "- [Numpy quickstart tutorial](https://numpy.org/devdocs/user/quickstart.html)\n",
    "- [Numpy chapter on scientific python lectures](https://lectures.scientific-python.org/intro/numpy/index.html)\n",
    "- [Stanford's CS231n lecture notes chapter on numpy](https://cs231n.github.io/python-numpy-tutorial/)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac2d563e-9991-4dad-8f76-628fa54994ed",
   "metadata": {},
   "source": [
    "## Exercise: `numpy` vs native looping over lists for adding two vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "136bceb8-8c59-45d4-8204-432beb58033d",
   "metadata": {},
   "source": [
    "First import numpy for use in our codes using the recommended syntax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2ce3dd45-9b35-4184-b9ce-8efd7d271c76",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8f0b5904-42d3-4f67-a61e-67b5965feb09",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000000  # choose a fairly large, realistic array size in scientific problems (say a million elements) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ee50f31-d3c4-4e4f-98a2-8af62872a182",
   "metadata": {},
   "source": [
    "Let us add two vectors of size N."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "204986f5-94d5-4bd0-bc19-51264e3509a0",
   "metadata": {},
   "source": [
    "### First using native python lists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c1ba24ab-ba27-431a-b8d1-a73d0e0fb929",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_native = range(N)\n",
    "Y_native = range(N)\n",
    "Z_native = [0]*N    # pre-allocate a result list of the same size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29c22e89-8ed6-4b5d-a081-2e91ea7f065c",
   "metadata": {},
   "source": [
    "### and loops ...."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "711411bb-69f3-4174-acfb-b5bebe3acd4e",
   "metadata": {},
   "source": [
    "%%time is a jupyter notebook magic directive that instructs that the computation of the cell must be timed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d40447fe-b16e-4246-a88a-c77f0dd75de6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 416 ms, sys: 19.3 ms, total: 436 ms\n",
      "Wall time: 435 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for i in range(N):\n",
    "    Z_native[i] = X_native[i] + Y_native[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0426cab-523b-4d2c-a64b-fab8c3ae61c5",
   "metadata": {},
   "source": [
    "### List comprehension vs loops"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9097dab-1d36-47e3-a07d-43a3d3982c70",
   "metadata": {},
   "source": [
    "Using python's list comprehension facilities can be ever so slightly faster than the conventional loop syntax (using an explicit loop statement and an indented block of body of code to be repeated). However, this still executes the loop to construct the result list, but has just a slightly more pythonic syntax that is sometimes preferred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8ed2e28c-aa05-47fc-9579-aa153c87cf88",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 278 ms, sys: 39.8 ms, total: 318 ms\n",
      "Wall time: 316 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "Z = [X_native[i] + Y_native[i] for i in range(N)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfc32198-58bc-479b-b872-0a7a6e9bb599",
   "metadata": {},
   "source": [
    "### Now, switch to using Numpy arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "778b67f7-f713-4487-9f7e-f5fecb054cad",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_numpy = np.arange(N)\n",
    "Y_numpy = np.arange(N)\n",
    "Z_numpy = np.zeros(N)    # pre-allocate a result vector initialised with zeros"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19d95a3d-ce03-478c-bd81-8279c7093cbd",
   "metadata": {},
   "source": [
    "### and performing vectorised addition instead of loops ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ecf151c1-f89f-4e76-b3b5-840aab720c3d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.08 ms, sys: 6.66 ms, total: 9.75 ms\n",
      "Wall time: 6.98 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "Z_numpy = X_numpy + Y_numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7876597-c39d-42f9-a994-8311a8a4c237",
   "metadata": {},
   "source": [
    "We achieved a few orders of magnitude speedup!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d54bbdb4-30d7-45f3-988e-f6a4c8441025",
   "metadata": {},
   "source": [
    "### Exercise: Accelerating the performance of CFD code with numpy arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29f8316d-ef4d-4166-8c54-87b1ff34acf6",
   "metadata": {},
   "source": [
    "In the `cfd_numpy` directory, we have provided a refactored version of the CFD code that uses Numpy arrays rather than Python lists, which also replaces loops with vectorised operations. Inspect and understand the `numpy` version of the code. In particular, notice how the absence of loops results in much reduced code clutter, which improves readability by conveying the program's intent better.\n",
    "\n",
    "Navigate to the `cfd_numpy` directory, and run the program with the same problem size as before, i.e. a 128 x 128 grid (scaling factor of 4) and 10000 Jacobi iterations:\n",
    "\n",
    "```bash\n",
    "prompt:/path/to/cfd_numpy> python cfd.py 4 10000\n",
    "```\n",
    "\n",
    "The timing output of the program should reveal a good improvement in the speed e.g.\n",
    "\n",
    "```\n",
    "2D CFD Simulation\n",
    "=================\n",
    "Scale factor = 4\n",
    "Iterations   = 10000\n",
    "\n",
    "Initialisation took 0.00024s\n",
    "\n",
    "Grid size = 128 x 128\n",
    "\n",
    "Starting main Jacobi loop...\n",
    "\n",
    "...finished\n",
    "\n",
    "Calculation took 1.52996s\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfe860ba-79db-437c-a90c-8602f0797721",
   "metadata": {},
   "source": [
    "### Need to go faster!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db4f82a3-29ad-48a8-b1ff-a8368be4976b",
   "metadata": {},
   "source": [
    "We achieved a decent speedup of the CFD code by replacing lists and loops with `numpy` arrays and vectorised operations, but let us investigate if further speed-up is possible."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}