Initial Commit

Tutorial/0-numpy.ipynb

@@ -0,0 +1,490 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# NumPy - Numeric Python\n",
+    "NumPy is the core library for scientifc computing Python. The module provides high-performance implementations for dealing with n-dimensional arrays and tools for working with these arrays."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## How to work with NumPy\n",
+    "The fundamental data structure in NumPy is called _array_. An array is a n-dimensional container. It can be used for \n",
+    "\n",
+    "* vectors\n",
+    "* matrices\n",
+    "* n-dimensional data\n",
+    "\n",
+    "Internally a NumPy array is stored row major (C order) by default. Note that this is different to Matlab where arrays are stored column major."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Import NumPy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib notebook\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Creating NumPy arrays"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# integer vector, 1D array\n",
+    "v = np.array([1, 2, 3])\n",
+    "v"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# matrix, 2D array\n",
+    "M = np.array([[1, 2], \n",
+    "              [3, 4]], dtype=float)\n",
+    "M"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# volume, 3D array\n",
+    "V = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8 , 9], [10, 11, 12]]])\n",
+    "V"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Initial placeholders"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create array of zeros with 3 rows and 4 columns\n",
+    "np.zeros((3, 4))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create array of ones \n",
+    "np.ones((2, 3, 4), dtype=np.int16)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create array of evenly spaced values\n",
+    "# args: start, stop, step: inclusive start, exclusive stop\n",
+    "np.arange(3, 11, 2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create array with evenly spaced values using number of samples\n",
+    "# start, stop, num_samples\n",
+    "x = np.linspace(-1, 4, 20)\n",
+    "\n",
+    "_, ax = plt.subplots(figsize=(5, 3.5))\n",
+    "ax.plot(x, np.sin(x))\n",
+    "\n",
+    "_, ax = plt.subplots(figsize=(5, 3.5))\n",
+    "ax.plot(np.sin(x))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create a constant array \n",
+    "np.full((2, 3), 7)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create a n x n identity matrix\n",
+    "np.eye(3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create an array with random values\n",
+    "np.random.random((2, 3))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create an empty array (attention: memory is uninitialized!)\n",
+    "np.empty((3, 2))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "base = np.array([1, 2, 3])\n",
+    "\n",
+    "a = []\n",
+    "for i in range(10):\n",
+    "    a.append(base.copy())\n",
+    "a = np.array(a)\n",
+    "    \n",
+    "b = np.stack([base.copy()] * 10)\n",
+    "c = np.repeat(base[None], 10, axis=0)\n",
+    "d = np.tile(base, (10, 1))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(a)\n",
+    "assert np.allclose(a, b) and np.allclose(b, c) and np.allclose(c, d)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Inspecting array"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# array dimensions\n",
+    "V.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# number of dimensions\n",
+    "V.ndim"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# number of elements\n",
+    "V.size"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# data type of array elements\n",
+    "V.dtype"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Array operations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# sum\n",
+    "V.sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# min\n",
+    "V.min()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# max\n",
+    "V.max()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# mean\n",
+    "V.mean()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# standard deviation\n",
+    "V.std()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Arithmetic Operations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "a = np.arange(3)\n",
+    "a"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# addition\n",
+    "a + a"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "a + 1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# exp\n",
+    "np.exp(a)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# * is elementwise multiplication\n",
+    "# @ is __matmul__\n",
+    "# a is one-dimensional\n",
+    "for inner in [(a * a).sum(), a @ a, a.dot(a), a.T.dot(a)]:\n",
+    "    print(inner)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for outer in [np.outer(a, a), a[None] * a[:, None]]:\n",
+    "    print(outer)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Slicing, Indexing\n",
+    "In python indexing starts at 0!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# get first element\n",
+    "a\n",
+    "a[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "M"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# get 1st row\n",
+    "M[0] # = M[0,:]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# get last column\n",
+    "M[:,-1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "v = np.arange(10)\n",
+    "v"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# get every second element\n",
+    "# [start:stop:step]\n",
+    "v[1::2]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# boolean indexing\n",
+    "v[v > 3]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "v > 3"
+   ]
+  }
+ ],
+ "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.10.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

Tutorial/1-images.ipynb

@@ -0,0 +1,339 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Working with images"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib notebook\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import imageio.v2 as imageio"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Image loading and display"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "image = imageio.imread('taj_mahal_gray.jpg')\n",
+    "\n",
+    "start, step = 300, 30\n",
+    "print(f'The image \\n{image}\\nhas {image.shape[0]} rows and {image.shape[1]} columns.')\n",
+    "print(f'The {step} pixels in the zero-th row, starting at column {start} contain the values \\n{image[0, start:start + step]}.')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "_, ax = plt.subplots(1, 2, figsize=(9, 3))\n",
+    "ax[0].imshow(image)\n",
+    "ax[1].imshow(image, cmap='gray')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Deep and shallow copy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "image = imageio.imread('taj_mahal_gray.jpg')\n",
+    "\n",
+    "# Assign another name to the same underlying storage (\"Reference\")\n",
+    "# copy = image\n",
+    "\n",
+    "# Copy the underlying data to a new storage location (\"Copy\")\n",
+    "copy = image.copy()\n",
+    "\n",
+    "copy[200:350, 300:450] = 0\n",
+    "\n",
+    "fig, ax = plt.subplots(1, 2, figsize=(9, 3))\n",
+    "ax[0].imshow(image)\n",
+    "ax[1].imshow(copy, cmap='gray')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Image statistics & manipulations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "im_integral = imageio.imread('taj_mahal_gray.jpg')\n",
+    "im_float = im_integral / 255.  # normalize grayvalues to [0,1]\n",
+    "\n",
+    "fig, axs = plt.subplots(1, 2, figsize=(9, 3))\n",
+    "for ax, im in zip(axs, [im_integral, im_float]):\n",
+    "    artist = ax.imshow(im, cmap='gray')\n",
+    "    ax.set_title(f'min={im.min()}, max={im.max()}')\n",
+    "    fig.colorbar(artist, ax=ax)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "im = imageio.imread('taj_mahal_small.jpg')\n",
+    "print(im.shape)\n",
+    "_, ax = plt.subplots()\n",
+    "ax.imshow(im)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Conversion to grayscale image"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "im = im / 255.0\n",
+    "\n",
+    "coeffs = np.array([0.299, 0.587, 0.114])\n",
+    "gray_explicit = im[:,:,0] * 0.299 + im[:,:,1] * 0.587 + im[:,:,2] * 0.114\n",
+    "gray_broadcasting = (im * coeffs).sum(axis=2)\n",
+    "assert np.allclose(gray_broadcasting, gray_explicit)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gray_mean = im.mean(-1)\n",
+    "\n",
+    "_, ax = plt.subplots(1, 2, figsize=(9, 3))\n",
+    "ax[0].imshow(gray_broadcasting, cmap='gray')\n",
+    "ax[1].imshow(gray_mean, cmap='gray')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": true
+   },
+   "source": [
+    "## Image filtering"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import scipy.ndimage as sn\n",
+    "from skimage.transform import rescale"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# load image\n",
+    "im = imageio.imread('taj_mahal_gray.jpg') / 255.0\n",
+    "im_aa = rescale(im, 0.2, anti_aliasing=True)\n",
+    "im_noaa = rescale(im, 0.2, anti_aliasing=False)\n",
+    "\n",
+    "_, ax = plt.subplots(1, 2, figsize=(9,3))\n",
+    "ax[0].imshow(im_aa, cmap='gray')\n",
+    "ax[0].set_title('AA')\n",
+    "ax[1].imshow(im_noaa, cmap='gray')\n",
+    "ax[1].set_title('No AA')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Gaussian filter"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sigma = 2.3\n",
+    "kernel_size = 15\n",
+    "\n",
+    "space = np.arange(-kernel_size // 2 + 1, kernel_size // 2 + 1)\n",
+    "x, y = np.meshgrid(space, space)\n",
+    "\n",
+    "kernel = np.exp(-(x ** 2 + y ** 2) / (2 * sigma ** 2))\n",
+    "kernel = kernel / np.sum(kernel)\n",
+    "_, ax = plt.subplots(subplot_kw={'projection': '3d'})\n",
+    "ax.plot_surface(x, y, kernel, cmap='viridis')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Apply gaussian filter to image"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# convolution with the gaussian filter\n",
+    "filtered = sn.convolve(im, kernel)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "_, ax = plt.subplots(1, 2, figsize=(9,3));\n",
+    "ax[0].imshow(im, cmap='gray')\n",
+    "ax[0].set_title('Original')\n",
+    "ax[1].imshow(filtered, cmap='gray')\n",
+    "ax[1].set_title('2D filter')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Edge filter"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dx = np.array([\n",
+    "    [-1, 0, 1],\n",
+    "    [-1, 0, 1],\n",
+    "    [-1, 0, 1]\n",
+    "])\n",
+    "dy = dx.T\n",
+    "\n",
+    "_, ax = plt.subplots(1, 2, figsize=(5,2))\n",
+    "ax[0].imshow(dx, cmap='gray')\n",
+    "ax[1].imshow(dy, cmap='gray')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nabla_x = sn.convolve(im, dx)\n",
+    "nabla_y = sn.convolve(im, dy)\n",
+    "grad_magnitude = np.sqrt(nabla_x ** 2 + nabla_y ** 2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fi, ax = plt.subplots(1, 3, figsize=(9,3))\n",
+    "ax[0].imshow(nabla_x, cmap='gray')\n",
+    "ax[1].imshow(nabla_y, cmap='gray')\n",
+    "ax[2].imshow(grad_magnitude, cmap='gray')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nabla_x = sn.convolve(filtered, dx)\n",
+    "nabla_y = sn.convolve(filtered, dy)\n",
+    "grad_magnitude = np.sqrt(nabla_x ** 2 + nabla_x ** 2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fi, ax = plt.subplots(1, 3, figsize=(9,3))\n",
+    "ax[0].imshow(nabla_x, cmap='gray')\n",
+    "ax[1].imshow(nabla_y, cmap='gray')\n",
+    "ax[2].imshow(grad_magnitude, cmap='gray')"
+   ]
+  }
+ ],
+ "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.10.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

Tutorial/2-implementation.ipynb

@@ -0,0 +1,361 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Explicit for-loop vs. built-ins"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def index_sum_for(arr, idx):\n",
+    "    sum_ = 0\n",
+    "    for a, i in zip(arr, idx):\n",
+    "        if i:\n",
+    "            sum_ += a\n",
+    "    return sum_\n",
+    "\n",
+    "\n",
+    "def index_sum_np(arr, idx):\n",
+    "    return arr[idx].sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_points = int(5e5)\n",
+    "array = np.random.randn((num_points))\n",
+    "indices = np.random.choice([False, True], size=num_points)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "44.7 ms ± 192 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit sum_for = index_sum_for(array, indices)\n",
+    "sum_for = index_sum_for(array, indices)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "2.17 ms ± 2.25 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit sum_np = index_sum_np(array, indices)\n",
+    "sum_np = index_sum_np(array, indices)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert np.allclose(sum_for, sum_np)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Pure allocation vs. allocation + initialization"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_points = int(1e4)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "6.5 µs ± 33.6 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit out = np.empty((num_points, num_points))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "72.1 ms ± 626 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit out = np.full((num_points, num_points), 1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Explicit for-loop vs. \"batched\" operations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rng = np.random.default_rng()\n",
+    "A = rng.standard_normal((3, 4))\n",
+    "batched_x = rng.standard_normal((1_000_000, 4))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def apply_matrix_for(A, xs):\n",
+    "    b = np.empty((xs.shape[0], A.shape[0]))\n",
+    "    for i, x in enumerate(xs):\n",
+    "        b[i] = A @ x\n",
+    "    return b"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def apply_matrix_batched(A, xs):\n",
+    "    return xs @ A.T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1.37 s ± 20.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit b_for = apply_matrix_for(A, batched_x)\n",
+    "b_for = apply_matrix_for(A, batched_x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "2.15 ms ± 270 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit b_batched = apply_matrix_batched(A, batched_x)\n",
+    "b_batched = apply_matrix_batched(A, batched_x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert np.allclose(b_for, b_batched)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Hints for Assignment 0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[ 6.,  2.,  3.,  4.],\n",
+       "       [ 5., 11.,  7.,  8.],\n",
+       "       [ 9., 10., 16., 12.],\n",
+       "       [13., 14., 15., 21.]])"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.set_printoptions(precision=1)\n",
+    "X = np.arange(1, 17, dtype=float).reshape(4, 4) + np.eye(4) * 5\n",
+    "kernel = np.array([\n",
+    "    [1, 2, 1],\n",
+    "    [2, 4, 2],\n",
+    "    [1, 2, 1],\n",
+    "]) / 16\n",
+    "X"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 49,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def variance_naive(X):\n",
+    "    M, N = X.shape\n",
+    "    var_im = np.zeros_like(X)\n",
+    "    for i in range(1, M - 1):\n",
+    "        for j in range(1, N - 1):\n",
+    "            x = X[i - 1:i + 2, j - 1:j + 2]\n",
+    "            mu_x = (x * kernel).sum()\n",
+    "            var_im[i, j] = (kernel * (x - mu_x) ** 2).sum()\n",
+    "\n",
+    "    return var_im"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 50,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0. , 0. , 0. , 0. ],\n",
+       "       [0. , 8.7, 9.4, 0. ],\n",
+       "       [0. , 7.9, 8.7, 0. ],\n",
+       "       [0. , 0. , 0. , 0. ]])"
+      ]
+     },
+     "execution_count": 50,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "naive = variance_naive(X)\n",
+    "naive"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 54,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import scipy.signal as ss\n",
+    "def variance(X):\n",
+    "    # in your code this corresponds to gaussian_filter calls\n",
+    "    mu_x = ss.convolve2d(X, kernel, mode='valid')\n",
+    "    mu_xx = ss.convolve2d(X * X, kernel, mode='valid')\n",
+    "    return np.pad(mu_xx - mu_x ** 2,  1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 55,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "smart = variance(X)\n",
+    "assert np.allclose(smart, naive)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.10.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

interpolation_error.py

@@ -0,0 +1,109 @@
+import numpy as np
+import imageio
+from skimage.transform import resize
+from scipy.ndimage import gaussian_filter
+
+
+def mssim(
+    x: np.ndarray,
+    y: np.ndarray,
+) -> float:
+    # Standard choice for the parameters
+    K1 = 0.01
+    K2 = 0.03
+    sigma = 1.5
+    truncate = 3.5
+    m = 1
+    C1 = (K1 * m) ** 2
+    C2 = (K2 * m) ** 2
+
+    # radius size of the local window (needed for
+    # normalizing the standard deviation)
+    r = int(truncate * sigma + 0.5)
+    win_size = 2 * r + 1
+    # use these arguments for the gaussian filtering
+    # e.g. filtered = gaussian_filter(x, **filter_args)
+    filter_args = {
+        'sigma': sigma,
+        'truncate': truncate
+    }
+
+    # Implement Eq. (9) from assignment sheet
+    # S should be an "image" of the SSIM evaluated for a window 
+    # centered around the corresponding pixel in the original input image
+    S = np.ones_like(x)
+
+    mu_x = gaussian_filter(x, **filter_args)
+    mu_y = gaussian_filter(y, **filter_args)
+
+    # Compute local variances and covariance
+    n = win_size ** 2 # number of pixels in the window
+    # n~ = n/(n-1)
+    # E[x^2] −E[x]^2). = (gaussian_filter(x * x, **filter_args) - mu_x * mu_x)
+    sigma_x_sq = (gaussian_filter(x * x, **filter_args) - mu_x * mu_x) * (n / (n - 1))
+    sigma_y_sq = (gaussian_filter(y * y, **filter_args) - mu_y * mu_y) * (n / (n - 1))
+    sigma_xy = (gaussian_filter(x * y, **filter_args) - mu_x * mu_y) * (n / (n - 1))
+    
+    # Compute SSIM
+    S = ((2 * mu_x * mu_y + C1) * (2 * sigma_xy + C2)) / \
+        ((mu_x ** 2 + mu_y ** 2 + C1) * (sigma_x_sq + sigma_y_sq + C2))
+    
+    # crop to remove boundary artifacts, return MSSIM
+    pad = (win_size - 1) // 2
+    return S[pad:-pad, pad:-pad].mean()
+
+
+def psnr(
+    x: np.ndarray,
+    y: np.ndarray,
+) -> float:
+    # Implement Eq. (2) without for loops
+    mse = np.mean((x - y) ** 2)
+    m = 1
+    psnr_value = 10 * np.log10(m ** 2 / mse)
+    return psnr_value
+
+
+def psnr_for(
+    x: np.ndarray,
+    y: np.ndarray,
+) -> float:
+    # Implement Eq. (2) using for loops
+    m, n = x.shape
+    mse = 0
+    # Eq. (1) from the assignment sheet for mse
+    for i in range(m):
+        for j in range(n):
+            mse += (x[i, j] - y[i, j]) ** 2
+    mse /= (m * n)
+    m = 1
+    psnr_value = 10 * np.log10(m ** 2 / mse)
+    return psnr_value
+
+
+def interpolation_error():
+    x = imageio.imread('./girl.png') / 255.
+    shape_lower = (x.shape[0] // 2, x.shape[1] // 2)
+    # downsample image to half the resolution
+    # and successively upsample to the original resolution
+    # using no nearest neighbor, linear and cubic interpolation
+    nearest, linear, cubic = [
+        resize(resize(
+            x, shape_lower, order=order, anti_aliasing=False
+        ), x.shape, order=order, anti_aliasing=False)
+        for order in [0, 1, 3]
+    ]
+
+    for label, rescaled in zip(
+        ['nearest', 'linear', 'cubic'],
+        [nearest, linear, cubic]
+    ):
+        print(label)
+        print(mssim(x, rescaled))
+        print(psnr(x, rescaled))
+        print(psnr_for(x, rescaled))
+        imageio.imwrite('girl_' + label + '.png', (rescaled * 255).astype(np.uint8))
+
+
+if __name__ == '__main__':
+    interpolation_error()

reference_output.txt

@@ -0,0 +1,12 @@
+nearest
+0.8031736958539067
+27.16729887950422
+27.167298879503985
+linear
+0.8890060634234201
+31.092116935553634
+31.092116935553662
+cubic
+0.9139147893771699
+31.97622728671044
+31.976227286710447

test_compare_output.py

@@ -0,0 +1,29 @@
+import unittest
+import subprocess
+
+class TestOutputComparison(unittest.TestCase):
+
+    def test_compare_output(self):
+        # Define the Python script to run and the expected output file
+        script_to_run = 'interpolation_error.py'  # Replace with your script's name
+        expected_file = 'reference_output.txt'
+
+        # Run the script and capture the output
+        result = subprocess.run(['python3', script_to_run], capture_output=True, text=True)
+        output_lines = result.stdout.splitlines()
+
+        # Open and read the expected output file
+        with open(expected_file, 'r') as expected:
+            expected_lines = expected.read().splitlines()
+
+            # Ensure the number of lines is the same
+            self.assertEqual(len(output_lines), len(expected_lines), 
+                             "The number of lines in the output does not match the expected output.")
+
+            # Compare each line
+            for i, (output_line, expected_line) in enumerate(zip(output_lines, expected_lines)):
+                self.assertEqual(output_line.strip(), expected_line.strip(), 
+                                 f"Line {i+1} does not match:\nOutput: {output_line}\nExpected: {expected_line}")
+
+if __name__ == '__main__':
+    unittest.main()