From af74414449709fbcc9e93fd63afe844ff6818adf Mon Sep 17 00:00:00 2001
From: radoskov <radoslav.skoviera@cvut.cz>
Date: Tue, 25 Feb 2025 22:53:28 +0100
Subject: [PATCH] Lecture 2 updates and fixes
---
.../lecture_02/l2_matrices_processing.ipynb | 3144 +++++++++--------
.../lecture_02/l2_matrices_processing.qmd | 71 +-
2 files changed, 1803 insertions(+), 1412 deletions(-)
diff --git a/src/pge_lectures/lecture_02/l2_matrices_processing.ipynb b/src/pge_lectures/lecture_02/l2_matrices_processing.ipynb
index f4b80f9..c548ac2 100644
--- a/src/pge_lectures/lecture_02/l2_matrices_processing.ipynb
+++ b/src/pge_lectures/lecture_02/l2_matrices_processing.ipynb
@@ -1,1391 +1,1763 @@
{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# 2D arrays (matrices), array processing (search, cumulative sum)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "---\n",
- "title: Programming for Engineers\n",
- "subtitle: \"Lecture 2 - matrices and array processing\"\n",
- "author: Radoslav Ĺ koviera\n",
- "\n",
- "format:\n",
- " pdf:\n",
- " code-fold: false\n",
- " html:\n",
- " code-fold: false\n",
- "jupyter: python3\n",
- "toc: true\n",
- "---\n",
- "\n",
- "## Prelude: Visualizing (printing) arrays\n",
- "\n",
- "For debugging purposes, it is useful to see what is in an array. Smaller arrays are easy to print \"whole\":"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "import numpy as np\n",
- "a = np.random.permutation(10)\n",
- "print(a)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "You can use the `join` method of strings to make the output nicer. More on that later, when we get to strings."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def pretty_print_array(a):\n",
- " print('[' + ', '.join([str(x) for x in a]) + ']')\n",
- "\n",
- "pretty_print_array(a)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "For larger arrays, you need to be creative - the solution will depend on what you need. For example, it might be enough to print the first or last few items. Print the array as rows (Python kinda does that but \"uncontrollably\"):"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def print_as_rows(a, row_length=10):\n",
- " for i in range(0, len(a), row_length):\n",
- " print(a[i:i+row_length])\n",
- "\n",
- "a = np.random.permutation(100)\n",
- "print_as_rows(a)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "You can even improve it to have nicer output. Again, more on that next time, when we get to strings."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def pretty_print_as_rows(a, row_length=10):\n",
- " print('[')\n",
- " for i in range(0, len(a), row_length):\n",
- " print(', '.join([str(x) for x in a[i:i+row_length]]), end=',\\n')\n",
- " print(']')\n",
- "\n",
- "a = np.random.permutation(100).tolist()\n",
- "pretty_print_as_rows(a)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Printing matrices:\n",
- "\n",
- "(remember to define/run cell with the `pretty_print` function above)"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def print_matrix(m):\n",
- " print('[', end='')\n",
- " n_rows = len(m)\n",
- " for ri, row in enumerate(m):\n",
- " print(row, end=',\\n' if ri < n_rows - 1 else '')\n",
- " print(']')\n",
- "\n",
- "def formatted_print_matrix(m): # works only for integers\n",
- " n_rows = len(m)\n",
- " maximum_value = abs(max([max(row) for row in m])) + 1e-6\n",
- " if maximum_value > 0:\n",
- " max_digits = int(np.ceil(np.log10(maximum_value)))\n",
- " else: # all zeros\n",
- " max_digits = 1\n",
- " print('[', end='')\n",
- " for ri, row in enumerate(m):\n",
- " row_prefix = ('' if ri == 0 else ' ') + '['\n",
- " row_numbers = ', '.join([f\"{x:{max_digits}d}\" for x in row])\n",
- " row_end=',\\n' if ri < n_rows - 1 else ']'\n",
- " print(row_prefix + row_numbers, end=row_end)\n",
- " print(']')\n",
- "\n",
- "mat = np.random.randint(0, 20, size=(17, 11)).tolist()\n",
- "print(\"Row-by-row print:\")\n",
- "print_matrix(mat)\n",
- "print()\n",
- "print(\"Formatted print:\")\n",
- "formatted_print_matrix(mat)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Matrices\n",
- "\n",
- "Matrices are 2D arrays. They have two dimensions: rows and columns. The number of rows is often called the height and the number of columns is often called the width. This is similar to image, although, images have typically flipped dimensions: x-axis is width (columns) and y-axis is height (rows)."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "m = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]\n",
- "formatted_print_matrix(m)\n",
- "also_m = [\n",
- " [1, 2, 3],\n",
- " [4, 5, 6],\n",
- " [7, 8, 9]\n",
- "]\n",
- "formatted_print_matrix(also_m)\n",
- "m_as_well = [[j + i * 3 for j in range(1, 4)] for i in range(3)]\n",
- "formatted_print_matrix(m_as_well)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Creating a matrix in a loop (the following two functions are equivalent):"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def create_counting_matrix(height, width):\n",
- " return [[j + i * width for j in range(1, width + 1)] for i in range(height)]\n",
- "\n",
- "def create_counting_matrix_unpacked(height, width):\n",
- " mat = []\n",
- " for i in range(height):\n",
- " row = []\n",
- " for j in range(1, width + 1):\n",
- " row.append(j + i * width)\n",
- " mat.append(row)\n",
- " return mat\n",
- "\n",
- "formatted_print_matrix(create_counting_matrix(3, 4))"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Matrix pre-allocation:"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def create_empty(n_rows, n_cols):\n",
- " return [[0] * n_cols for _ in range(n_rows)]\n",
- "M, N = 3, 4\n",
- "mat = create_empty(M, N)\n",
- "print(f\"Empty {M} by {N} matrix:\")\n",
- "formatted_print_matrix(mat)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Filling matrices with values\n",
- "\n",
- "Filling matrix with a constant value. *Warning*, the following function might fail in case of \"jagged\" matrices - this is why it's best to only use 2D arrays (matrices) with equal column lengths (second dimension). There is a simple workaround - get `n_col` at each row-loop. Nonetheless, the best approach is to avoid jagged matrices in the first place."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def fill_matrix_with_value(matrix, value):\n",
- " n_row = len(matrix)\n",
- " if n_row == 0: # let's check if matrix is empty to avoid some errors\n",
- " print(\"Matrix is empty!\")\n",
- " return\n",
- " n_col = len(matrix[0])\n",
- " for r in range(n_row):\n",
- " for c in range(n_col):\n",
- " matrix[r][c] = value\n",
- "\n",
- "\n",
- "m = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]\n",
- "fill_matrix_with_value(m, 11)\n",
- "formatted_print_matrix(m)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Random fill:"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "import random\n",
- "\n",
- "def fill_matrix_with_randint(matrix, min_value=0, max_value=100):\n",
- " n_row = len(matrix)\n",
- " if n_row == 0: # let's check if matrix is empty to avoid some errors\n",
- " print(\"Matrix is empty!\")\n",
- " return\n",
- " n_col = len(matrix[0])\n",
- "\n",
- " for r in range(n_row):\n",
- " for c in range(n_col):\n",
- " matrix[r][c] = random.randint(min_value, max_value)\n",
- "\n",
- "mat = [[1] * 5 for _ in range(4)]\n",
- "fill_matrix_with_randint(mat)\n",
- "formatted_print_matrix(mat)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Shuffle matrix elements. There are few ways to do it. Here, we will create \"flattened\" indices, shuffle them and assign values to them from the original matrix."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def shuffle_matrix(matrix):\n",
- " n_row = len(matrix)\n",
- " m_cols = len(matrix[0])\n",
- " flat_indices = list(range(n_row * m_cols))\n",
- " random.shuffle(flat_indices)\n",
- " shuffled_matrix = create_empty(n_row, m_cols)\n",
- " for r in range(n_row):\n",
- " for c in range(m_cols):\n",
- " flat_index = r * m_cols + c\n",
- " shuffled_matrix[r][c] = matrix[int(flat_indices[flat_index] / m_cols)][flat_indices[flat_index] % m_cols]\n",
- "\n",
- " return shuffled_matrix\n",
- "\n",
- "m = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]\n",
- "sm = shuffle_matrix(m)\n",
- "print(\"Original matrix:\")\n",
- "formatted_print_matrix(m)\n",
- "print(\"Shuffled matrix:\")\n",
- "formatted_print_matrix(sm)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Matrix operations\n",
- "\n",
- "First, let's create some matrices:"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "mat_A = create_empty(3, 4)\n",
- "fill_matrix_with_randint(mat_A, 0, 20)\n",
- "mat_B = create_empty(3, 4)\n",
- "fill_matrix_with_randint(mat_B, 0, 20)\n",
- "print(\"Matrix A:\")\n",
- "formatted_print_matrix(mat_A)\n",
- "print(\"Matrix B:\")\n",
- "formatted_print_matrix(mat_B)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Addition\n",
- "\n",
- "Compute $C = A + B$. This one is easy, we just need to loop through elements of the two matrices and add them."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def add_two_matrices(mat_A, mat_B):\n",
- " n_row = len(mat_A)\n",
- " n_col = len(mat_A[0])\n",
- " mat_C = create_empty(n_row, n_col)\n",
- " for r in range(n_row):\n",
- " for c in range(n_col):\n",
- " mat_C[r][c] = mat_A[r][c] + mat_B[r][c]\n",
- " return mat_C\n",
- "\n",
- "mat_C = add_two_matrices(mat_A, mat_B)\n",
- "print(\"Result of addition of two matrices:\")\n",
- "formatted_print_matrix(mat_C)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Transposition\n",
- "\n",
- "We want to compute $A^T$. We just swap the indices for the output matrix."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def transpose(mat):\n",
- " n_row = len(mat) # get the matrix dimensions\n",
- " n_col = len(mat[0])\n",
- " mat_T = create_empty(n_col, n_row) # create a new matrix\n",
- " for r in range(n_row):\n",
- " for c in range(n_col):\n",
- " mat_T[c][r] = mat[r][c] # swap the row and column indices\n",
- " return mat_T\n",
- "\n",
- "mat_T = transpose(mat_A)\n",
- "print(\"Transposed matrix:\")\n",
- "formatted_print_matrix(mat_T)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Multiplication\n",
- "\n",
- "Compute $C = A \\times B$. It is not simple, we need to loop through the elements of the matrices and compute their product. The result will be a new matrix with dimensions $n \\times m$."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def multiply_two_matrices(mat_A, mat_B):\n",
- " n_row_A = len(mat_A)\n",
- " n_col_A = len(mat_A[0])\n",
- " n_row_B = len(mat_B)\n",
- " n_col_B = len(mat_B[0])\n",
- "\n",
- " if n_col_A != n_row_B: # the matrices need to have the correct shapes\n",
- " raise ValueError(f\"Matrices A and B cannot be multiplied. \"\n",
- " f\"Matrix B needs to have the same number of columns (has {n_col_B}) as matrix A has rows (has {n_row_A}).\")\n",
- " # this function will not be able to deal with broadcasting!\n",
- " mat_C = create_empty(n_row_A, n_col_B)\n",
- " for r in range(n_row_A):\n",
- " for c in range(n_col_B):\n",
- " mat_C[r][c] = 0\n",
- " for k in range(n_col_A):\n",
- " mat_C[r][c] += mat_A[r][k] * mat_B[k][c]\n",
- " return mat_C\n",
- "\n",
- "mat_C = multiply_two_matrices(mat_A, transpose(mat_B))\n",
- "print(\"Result of multiplication of two matrices:\")\n",
- "formatted_print_matrix(mat_C)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Miscellaneous operations on matrices\n",
- "\n",
- "#### Comparison\n",
- "\n",
- "Compare two matrices \"plain\" - just iterate over the rows & columns and compare corresponding elements in both matrices. Sometimes, we might want to compare elements in two matrices, regardless of their position. In that case, it is easier to \"flatten\" one of the matrices - it is easier to iterate over a 1D array. The 'unique' version of the second function requires us to remove duplicates from the flattened matrix first (you can remove them from any matrix but again, it is easier to do it for a 1D array). Otherwise, we would need to remove them at the end, when comparing the values."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def compare_two_matrices(mat_A, mat_B):\n",
- " n_row = len(mat_A)\n",
- " n_col = len(mat_A[0])\n",
- " if n_row != len(mat_B) or n_col != len(mat_B[0]):\n",
- " return False # matrices don't have the same shape => cannot be equal\n",
- " for r in range(n_row):\n",
- " for c in range(n_col):\n",
- " if mat_A[r][c] != mat_B[r][c]:\n",
- " return False # \"early\" termination for efficiency\n",
- " return True\n",
- "\n",
- "def remove_duplicates(arr):\n",
- " i = 0\n",
- " while i < len(arr) - 1:\n",
- " j = i + 1\n",
- " while j < len(arr):\n",
- " if arr[i] == arr[j]:\n",
- " del arr[j]\n",
- " j += 1\n",
- " i += 1\n",
- "\n",
- "def find_equal_values(mat_A, mat_B, unique=False):\n",
- " n_row_A = len(mat_A)\n",
- " n_col_A = len(mat_A[0])\n",
- " n_row_B = len(mat_B)\n",
- " n_col_B = len(mat_B[0])\n",
- " flat_values_B = []\n",
- " for r in range(n_row_B):\n",
- " flat_values_B.extend(mat_B[r])\n",
- "\n",
- " if unique: # remove duplicates, otherwise the `remove` method below might not be enough\n",
- " remove_duplicates(flat_values_B)\n",
- "\n",
- " equal_values = []\n",
- " for r in range(n_row_A):\n",
- " for c in range(n_col_A):\n",
- " item_A = mat_A[r][c]\n",
- " if item_A in flat_values_B:\n",
- " equal_values.append(item_A)\n",
- " if unique:\n",
- " flat_values_B.remove(item_A)\n",
- " return equal_values\n",
- "\n",
- "mat_A = create_empty(3, 4)\n",
- "fill_matrix_with_randint(mat_A, 0, 10)\n",
- "mat_B = create_empty(3, 4)\n",
- "fill_matrix_with_randint(mat_B, 0, 10)\n",
- "print(\"Matrix A:\")\n",
- "formatted_print_matrix(mat_A)\n",
- "print(\"Matrix B:\")\n",
- "formatted_print_matrix(mat_B)\n",
- "\n",
- "print(\"Matrices are equal:\", compare_two_matrices(mat_A, mat_B))\n",
- "print(\"Equal values:\", find_equal_values(mat_A, mat_B))\n",
- "print(\"Uniquely equal values:\", find_equal_values(mat_A, mat_B, unique=True))"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### All different\n",
- "\n",
- "Find whether the matrix contains only unique values. One option is to use the function above. Then, the number of uniquely equal values with itself must equal number of elements."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "# print(\"Matrix contains only unique values:\", len() == 0)\n",
- "mat_U = create_empty(3, 4)\n",
- "fill_matrix_with_randint(mat_U, 0, 50)\n",
- "formatted_print_matrix(mat_U)\n",
- "num_equal_values = len(find_equal_values(mat_U, mat_U, unique=True))\n",
- "if num_equal_values == len(mat_U) * len(mat_U[0]):\n",
- " print(\"Matrix contains only unique values.\")\n",
- "else:\n",
- " print(\"Matrix does not contain only unique values.\")"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Alternatively, a custom method with similar approach can be made. The easiest approach is to flatten the matrix and then basically compare \"two\" arrays."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def all_different(mat):\n",
- " # flatten the matrix\n",
- " flat_values = []\n",
- " for r in range(len(mat)):\n",
- " flat_values.extend(mat[r])\n",
- " n = len(flat_values)\n",
- " for i in range(n - 1):\n",
- " for j in range(i + 1, n):\n",
- " if flat_values[i] == flat_values[j]:\n",
- " return False\n",
- " return True\n",
- "\n",
- "mat_U = create_empty(3, 4)\n",
- "fill_matrix_with_randint(mat_U, 0, 50)\n",
- "formatted_print_matrix(mat_U)\n",
- "if all_different(mat_U):\n",
- " print(\"Matrix contains only unique values.\")\n",
- "else:\n",
- " print(\"Matrix does not contain only unique values.\")"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### \"Graphical\" matrices (images)\n",
- "\n",
- "Raw images are basically matrices, where each element represents an intensity value of the corresponding pixel. Color images are 3D matrices, where each element is a 3-tuple of intensity values of the color channels but we will not go there.\n",
- "\n",
- "However, we can have a little fun with gray-scale images."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "img = [\n",
- " [100, 100, 100, 100, 100, 100, 100, 100, 100, 100],\n",
- " [100, 0 , 0 , 0, 0 , 0 , 0, 0, 0, 100],\n",
- " [ 0, 0, 128, 128, 0, 0, 128, 128, 0, 0],\n",
- " [ 0, 0, 255, 255, 0, 0, 255, 255, 0, 0],\n",
- " [ 0, 0, 255, 255, 0, 0, 255, 255, 0, 0],\n",
- " [ 0, 0, 128, 128, 0, 0, 128, 128, 0, 0],\n",
- " [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
- " [ 0, 0, 0, 0, 0, 200, 0, 0, 0, 0],\n",
- " [ 0, 0, 0, 0, 0, 200, 0, 0, 0, 0],\n",
- " [ 0, 0, 0, 0, 200, 200, 0, 0, 0, 0],\n",
- " [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
- " [ 0, 255, 0, 0, 0, 0, 0, 0, 255, 0],\n",
- " [ 0, 255, 255, 255, 255, 255, 255, 255, 255, 0],\n",
- " [ 0, 0, 255, 128, 128, 128, 128, 255, 0, 0],\n",
- " [ 0, 0, 0, 255, 255, 255, 255, 0, 0, 0],\n",
- " [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
- " [ 0, 0, 0, 0, 100, 100, 0, 0, 0, 0],\n",
- "]\n",
- "\n",
- "print(img)\n",
- "\n",
- "from matplotlib import pyplot as plt\n",
- "\n",
- "plt.imshow(img, cmap='gray')\n",
- "plt.axis('off')\n",
- "plt.show()"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Array & matrix processing\n",
- "\n",
- "### Search\n",
- "\n",
- "#### Linear search\n",
- "\n",
- "Requires iterating over the array, asymptotic runtime is $O(n)$."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "import numpy as np\n",
- "N = 100\n",
- "a = np.random.permutation(N)\n",
- "searched_item = 50\n",
- "\n",
- "def linear_search(a, x):\n",
- " for i in range(len(a)):\n",
- " if a[i] == x:\n",
- " return i\n",
- " return None\n",
- "\n",
- "print(\"Linear search run time:\")\n",
- "%timeit linear_search(a, searched_item)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "##### Computational complexity side-quest\n",
- "\n",
- "Even though the 'asymptotic' time of linear search is $O(n)$, the 'actual' time depends on the 'properties' of the array - where the item is located in the array."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "N = 100\n",
- "a = np.random.permutation(N)\n",
- "searched_item = a[N // 2]\n",
- "print(\"Searching for an item in the middle of the array:\")\n",
- "%timeit linear_search(a, searched_item)\n",
- "\n",
- "searched_item = a[0]\n",
- "print(\"Searching for an item at the beginning of the array:\")\n",
- "%timeit linear_search(a, searched_item)\n",
- "\n",
- "searched_item = a[-1]\n",
- "print(\"Searching for an item at the end of the array:\")\n",
- "%timeit linear_search(a, searched_item)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Binary search\n",
- "\n",
- "Requires pre-sorted array but asymptotic runtime is $O(\\log n)$ - we only look which \"half\" of the array might contain the item."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def binary_search_loop(a, x):\n",
- " low = 0 # lower-bound index for search\n",
- " high = len(a) - 1 # upper-bound index for search\n",
- " loop_count = 1\n",
- " while low <= high: # while there are indices to search\n",
- " mid = (low + high) // 2 # compute the middle index\n",
- " if a[mid] < x: # if the middle point is less than the item\n",
- " low = mid + 1 # \"discard\" the left half (search from the middle + 1)\n",
- " elif a[mid] > x: # if the middle point is greater than the item\n",
- " high = mid - 1 # \"discard\" the right half (search from the middle - 1)\n",
- " else: # otherwise, the middle point is the item\n",
- " return mid, loop_count\n",
- " loop_count += 1\n",
- " return None, loop_count # we did not find the item\n",
- "\n",
- "def binary_search_recursion(a, x):\n",
- " if len(a) == 0: # we got an empty array\n",
- " return None # the item is for sure not in it\n",
- " mid_idx = len(a) // 2 # compute the middle point index\n",
- " if a[mid_idx] == x: # if the middle point is the item\n",
- " return mid_idx # we found it, job done\n",
- " elif a[mid_idx] < x: # if the middle point is less than the item\n",
- " # we look at the right half of the array\n",
- " return binary_search_recursion(a[mid_idx + 1:], x)\n",
- " else:\n",
- " # otherwise, we look at the left half of the array\n",
- " return binary_search_recursion(a[:mid_idx], x)\n",
- "\n",
- "# recursion is slower in this case, as we are copying the whole array, instead of just using the index.\n",
- "sa = np.sort(a)\n",
- "print(\"Binary search run time:\")\n",
- "%timeit binary_search_loop(sa, 5)\n",
- "%timeit binary_search_recursion(sa, 5)\n",
- "\n",
- "print(\"Buuuuuut...\")\n",
- "print(\"'Fair' binary search run time (includes sorting):\")\n",
- "%timeit binary_search_loop(np.sort(a), 5)\n",
- "%timeit binary_search_recursion(np.sort(a), 5)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The advantage of binary search is relatively constant search time. As you can see from the following, the search time is shorter than \"average\" search time for linear search: $O(\\log n) < O(n)$. This includes even sorting the array."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "print(\"Searching for an item in the middle of the array, using binary search:\")\n",
- "searched_item = sa[N // 2]\n",
- "%timeit binary_search_loop(np.sort(a), searched_item)\n",
- "\n",
- "print(\"Searching for an item at the beginning of the array, using binary search:\")\n",
- "searched_item = sa[0]\n",
- "%timeit binary_search_loop(np.sort(a), searched_item)\n",
- "\n",
- "print(\"Searching for an item at the end of the array, using binary search:\")\n",
- "searched_item = sa[-1]\n",
- "%timeit binary_search_loop(np.sort(a), searched_item)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Interpolation search\n",
- "\n",
- "Even more restrictive requirements than binary search - the array needs to be 'uniformly' (equal distribution of values) sorted. It uses similar approach as binary search but uses \"an educated guess\" of where the item might be (this is where the uniformity comes into play). The is $O(\\log n)$"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def interpolation_search(a, x):\n",
- " low = 0 # lower-bound index for search\n",
- " high = len(a) - 1 # upper-bound index for search\n",
- " loop_count = 1\n",
- " while low <= high and a[low] <= x <= a[high]: # while there are indices to search\n",
- " mid = low + int(((x - a[low]) * (high - low)) / (a[high] - a[low]))\n",
- " if a[mid] < x: # if the middle point is less than the item\n",
- " low = mid + 1 # \"discard\" the left half (search from the middle + 1)\n",
- " elif a[mid] > x: # if the middle point is greater than the item\n",
- " high = mid - 1 # \"discard\" the right half (search from the middle - 1)\n",
- " else: # otherwise, the middle point is the item\n",
- " return mid, loop_count\n",
- " loop_count += 1\n",
- " return None, loop_count # we did not find the item\n",
- "\n",
- "print(\"Searching for an item in the middle of the array, using interpolation search:\")\n",
- "searched_item = sa[N // 2]\n",
- "%timeit interpolation_search(np.sort(a), searched_item)\n",
- "\n",
- "print(\"Searching for an item at the beginning of the array, using interpolation search:\")\n",
- "searched_item = sa[0]\n",
- "%timeit interpolation_search(np.sort(a), searched_item)\n",
- "\n",
- "print(\"Searching for an item at the end of the array, using interpolation search:\")\n",
- "searched_item = sa[-1]\n",
- "%timeit interpolation_search(np.sort(a), searched_item)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The awesome speed of interpolation search:"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "a = np.random.permutation(N)\n",
- "sa = np.sort(a)\n",
- "searched_item = sa[N // 2 + 1]\n",
- "print(\"Interpolation search run time:\")\n",
- "%timeit interpolation_search(sa, searched_item)\n",
- "v, i_loop_count = interpolation_search(sa, searched_item)\n",
- "print(f\"Interpolation search 'step' count: {i_loop_count}\")\n",
- "\n",
- "print(\"Binary search run time:\")\n",
- "%timeit binary_search_loop(sa, searched_item)\n",
- "b, b_loop_count = binary_search_loop(sa, searched_item)\n",
- "print(f\"Binary search 'step' count: {b_loop_count}\")"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Unfortunately, the \"single step\" search only works if the array is completely uniform."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "a = np.random.choice(N**3, size=N, replace=False) / 100 # possibly non-uniform array\n",
- "sa = np.sort(a)\n",
- "searched_item = sa[N // 2 + 1]\n",
- "\n",
- "print(\"Interpolation search run time:\")\n",
- "%timeit interpolation_search(sa, searched_item)\n",
- "v, i_loop_count = interpolation_search(sa, searched_item)\n",
- "print(f\"Interpolation search 'step' count: {i_loop_count}\")\n",
- "print(\"Binary search run time:\")\n",
- "%timeit binary_search_loop(sa, searched_item)\n",
- "v, b_loop_count = binary_search_loop(sa, searched_item)\n",
- "print(f\"Binary search 'step' count: {b_loop_count}\")"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Find all\n",
- "\n",
- "The search finds the first occurrence of the searched item. Which one will depend on the search method. Sometimes, it is required to find all occurrences of the searched item in the array."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "import numpy as np\n",
- "def find_all(a, x):\n",
- " indices = []\n",
- " for i in range(len(a)):\n",
- " if a[i] == x:\n",
- " indices.append(i)\n",
- " return indices\n",
- "\n",
- "N = 100\n",
- "a = np.random.choice(N**3, size=N, replace=True)\n",
- "searched_item = a[N // 2]\n",
- "occurrences = find_all(a, searched_item)\n",
- "print(\"The value\", searched_item, \"occurs\", len(occurrences), \"times in the array at indices\", occurrences)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 'Statistical' computations on arrays\n",
- "\n",
- "#### Simple operations\n",
- "\n",
- "##### Minimum, maximum\n",
- "\n",
- "Python has built-in `min()` and `max()` functions. Here, we wil see simple implementation of these functions, simply to practice working with arrays. In practice, you can use the built-in functions."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "import numpy as np\n",
- "my_list = np.random.permutation(10).tolist()\n",
- "print(\"Input array:\", my_list)\n",
- "\n",
- "def minimum(a):\n",
- " min_value = a[0]\n",
- " for i in range(1, len(a)):\n",
- " if a[i] < min_value:\n",
- " min_value = a[i]\n",
- " # Alternatively:\n",
- " # min_value = min(min_value, a[i])\n",
- " return min_value\n",
- "\n",
- "def maximum(a):\n",
- " max_value = a[0]\n",
- " for i in range(1, len(a)):\n",
- " if a[i] > max_value:\n",
- " max_value = a[i]\n",
- " # Alternatively:\n",
- " # max_value = max(max_value, a[i])\n",
- " return max_value\n",
- "\n",
- "print(\"Minimum: \", minimum(my_list))\n",
- "print(\"Maximum: \", maximum(my_list))"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "##### Mean\n",
- "\n",
- "First, we need to compute the sum:"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "import numpy as np\n",
- "my_list = np.random.permutation(10).tolist()\n",
- "print(\"Input array:\", my_list)\n",
- "\n",
- "def sum(a):\n",
- " s = 0\n",
- " for i in range(len(a)):\n",
- " s += a[i]\n",
- " return s\n",
- "\n",
- "print(\"Sum: \", sum(my_list))"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Mean is then simply the sum divided by the number of items:"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def mean(a):\n",
- " return sum(a) / len(a)\n",
- "\n",
- "print(\"Mean: \", mean(my_list))"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Cumulative sum (prefix sum)\n",
- "\n",
- "Summing items between two indices is useful for many tasks. For example, computing average temperature between two specified dates. The following method computes sum between two indices in an array (inclusive of the start and the end indices)."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "import numpy as np\n",
- "temp_measurements = (np.random.rand(200) * 10).tolist()\n",
- "print(\"Input array: \" + ', '.join([f\"{x:.2f}\" for x in temp_measurements[:10]]) + \", ...\")\n",
- "\n",
- "def range_sum(a, start, end):\n",
- " s = 0\n",
- " for i in range(start, end + 1): # +1 to include the last element\n",
- " s += a[i]\n",
- " return s\n",
- "\n",
- "print(\"Sum between indices 5 and 15:\", range_sum(temp_measurements, 5, 15))"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "This is fine, if we need to to this once, but what if we need to do it many times?"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def generate_range_queries(a, n, query_len):\n",
- " numel_a = len(a)\n",
- " assert query_len <= numel_a, \"Query length must be less than the array length\"\n",
- " max_pos = numel_a - query_len\n",
- " queries = []\n",
- " for _ in range(n):\n",
- " start = np.random.randint(0, max_pos)\n",
- " end = start + query_len - 1\n",
- " queries.append((start, end))\n",
- " return queries\n",
- "\n",
- "print(\"Summing once run time:\")\n",
- "%timeit range_sum(temp_measurements, 50, 150)\n",
- "print(\"Summing 100 times run time:\")\n",
- "queries = generate_range_queries(temp_measurements, 100, 50)\n",
- "%timeit [range_sum(temp_measurements, start, end) for start, end in queries]"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Is there a better way? Why yes! We will use **cumulative sum** (also known as prefix sum). It is done bone adding all the previous values to the current value in the array."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def cumulative_sum(a):\n",
- " cs = [0] * len(a) # pre-allocate\n",
- " for i in range(1, len(cs)):\n",
- " cs[i] = a[i] + cs[i - 1] # add previous sum to the current element\n",
- " return cs\n",
- "my_list = list(range(10))\n",
- "print(\"Cumulative sum of ordered numbers from 0 to 9: \")\n",
- "print(\"Input array: [\" + ', '.join([f\"{x:2d}\" for x in my_list]) + \"]\")\n",
- "print(\"Cumulative sum: \" + str(cumulative_sum(my_list)))\n",
- "print()\n",
- "my_list = np.random.permutation(10).tolist()\n",
- "print(\"Cumulative sum of randomly shuffled numbers: \")\n",
- "print(\"Input array: [\" + ', '.join([f\"{x:2d}\" for x in my_list]) + \"]\")\n",
- "print(\"Cumulative sum: \" + str(cumulative_sum(my_list)))"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "| index | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |\n",
- "|--------|-----|-----|-------|---------|-----|-----|-----|-----|-----|-----|\n",
- "| value | 4 | 2 | 0 | 8 | 1 | 5 | 7 | 9 | 3 | 6 |\n",
- "| cumsum | 4 | 6 | 6 | 14 | 15 | 20 | 27 | 36 | 39 | 45 |\n",
- "| | 4 | 4+2 | 4+2+0 | 4+2+0+8 | ... | | | | | |\n",
- "\n",
- ": \"Visualizing\" cumulative sum\n",
- "\n",
- "If we have the cumulative sum precomputed, to get sum of elements between two indices $a$ and $b$, we can simply compute $cumsum[b] - cumsum[a-1]$. Therefore, instead of $b-a$ additions, we only compute one addition."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "my_list = np.random.permutation(10).tolist()\n",
- "\n",
- "def range_sum_cs(a_cumsum, start, end):\n",
- " return a_cumsum[end] - a_cumsum[start - 1]\n",
- "\n",
- "print(my_list)\n",
- "start, end = 3, 7\n",
- "print(f\"Sum between indices {start} and {end}:\", range_sum_cs(cumulative_sum(my_list), start, end))\n",
- "# sanity check with the \"simple\" method:\n",
- "print(\"This is the same as with the simple `range_sum` method:\", range_sum(my_list, start, end) == range_sum_cs(cumulative_sum(my_list), start, end))"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Now, let's say we have the cumulative sum precomputed and we want to compute the range sum for 100 different ranges. How long will it take?"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "print(\"Run time of cumulative sum computation:\")\n",
- "%timeit cumulative_sum(temp_measurements)\n",
- "\n",
- "cumsum_measurements = cumulative_sum(temp_measurements)\n",
- "print(\"Summing 100 times run time:\")\n",
- "%timeit [range_sum_cs(cumsum_measurements, start, end) for start, end in queries]"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Of course, the actual effectiveness of the cumulative / prefix sum depends on the number of queries and the length of the queries vs. the length of the array. For a few short queries in a long array, it might take longer to just compute the cumulative sum than to compute all the queries."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "short_queries = generate_range_queries(temp_measurements, 10, 5)\n",
- "print(\"Run time of simple range sum:\")\n",
- "%timeit [range_sum(temp_measurements, start, end) for start, end in short_queries]\n",
- "print(\"Run time of cumulative sum & range sum:\")\n",
- "%timeit cumsum_measurements = cumulative_sum(temp_measurements); [range_sum_cs(cumsum_measurements, start, end) for start, end in short_queries]"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Remember to always use the right tool for the job!\n",
- "\n",
- "Now, we can compute average values over varying ranges:"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "def range_average(csa, start, end):\n",
- " return range_sum_cs(csa, start, end) / (end - start + 1)\n",
- "\n",
- "a = list(range(10))\n",
- "csa = cumulative_sum(a)\n",
- "print(\"Input array:\", a)\n",
- "print(\"Average value of elements between indices 2 and 5 (inclusive):\", range_average(csa, 2, 5))"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Integral image (summed-area table, cumulative sum in 2D)\n",
- "\n",
- "In computer vision, a useful extension of the cumulative sum to 2 dimensions is used. It is called the **integral image** (or summed-area table)."
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "matrix = create_empty(11, 10)\n",
- "fill_matrix_with_randint(matrix, 0, 255)\n",
- "\n",
- "def integral_image(mat):\n",
- " ii = create_empty(len(mat), len(mat[0]))\n",
- " for i in range(len(mat)):\n",
- " for j in range(len(mat[0])):\n",
- " ii[i][j] = mat[i][j]\n",
- " if i > 0:\n",
- " ii[i][j] += ii[i - 1][j]\n",
- " if j > 0:\n",
- " ii[i][j] += ii[i][j - 1]\n",
- " if i > 0 and j > 0:\n",
- " ii[i][j] -= ii[i - 1][j - 1]\n",
- " return ii\n",
- "\n",
- "ii = integral_image(matrix)\n",
- "formatted_print_matrix(matrix)\n",
- "formatted_print_matrix(ii)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Now, to compute the sum of a sub-matrix, the same principle as in the cumulative \"range sum\" is applied. The \"only\" difference is that now we are in 2D. Therefore, the computation is slightly more complicated. To compute the sum over the \"area\" (sub-matrix) demarked by the indices $(a, b)$ and $(c, d)$, we need to compute: $ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]$. Let's visualize the integral image and this computation. Let us consider the problem to be defined as follows:"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "top_left = (5, 3) # a, b\n",
- "bottom_right = (7, 6) # c, d"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Let's break it down: We want to compute the sum of the gray area - area of interest (AOI) in @fig-integral_image_full. To do that, we need to take the \"total sum\" (delimited by red rectangle in @fig-integral_image_full) of the area from $[0, 0]$ to $[c, d]$. That is, the sum from the origin of the image to the bottom_right corner of the AOI. This sum has the value at $ii[c][d]$. Then, we subtract the two areas from origin to the top-right corner of the AOI (delimited by blue rectangle in @fig-integral_image_full) and the area from the origin to the bottom-left corner of the AOI (delimited by green rectangle in @fig-integral_image_full). The sums of these areas are located at $ii[a-1][d]$ and $ii[c][b-1]$. This way, we subtracted twice the sum of the area from the origin to just above the top-left corner of the AOI. We need to add it once back-in. Therefore, the last step is to add the sum of this area (delimited by orange rectangle in @fig-integral_image_full) that is located at $ii[a-1][b-1]$. To recap, the full equation is:\n",
- "\n",
- "$$area_sum[[a, b], [c, d]] = ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]$$\n",
- "\n",
- "Specifically, in our case:\n",
- "\n",
- "$$area_sum[[5, 3], [7, 6]] = ii[7][6] - ii[4][6] - ii[7][3] + ii[4][3]$$\n",
- "\n",
- "| general term | current case term | color in @fig-integral_image_full |\n",
- "|:-----------------|:------------------|:----------------------------------|\n",
- "| $+ ii[c][d]$ | $+ii[7][6]$ | red |\n",
- "| $- ii[a-1][d]$ | $-ii[4][6]$ | blue |\n",
- "| $- ii[c][b-1]$ | $-ii[7][3]$ | green |\n",
- "| $+ ii[a-1][b-1]$ | $+ii[4][3]$ | orange |"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "#| code-fold: true\n",
- "#| echo: false\n",
- "#| label: fig-integral_image_full\n",
- "#| fig-cap: Integral image with colored regions.\n",
- "import matplotlib as mpl\n",
- "import matplotlib.pyplot as plt\n",
- "\n",
- "# The following code is just for visualization. No need to understand it.\n",
- "def draw_matrix(matrix):\n",
- " rows = len(matrix)\n",
- " cols = len(matrix[0])\n",
- "\n",
- " fig, ax = plt.subplots(figsize=(cols, rows))\n",
- "\n",
- " # Create empty white image\n",
- " ax.imshow([[1]*cols for _ in range(rows)], cmap='gray', vmin=0, vmax=1)\n",
- "\n",
- " # Set ticks to be at the center of each cell\n",
- " ax.set_xticks(range(cols))\n",
- " ax.set_yticks(range(rows))\n",
- " ax.set_xticklabels(range(cols))\n",
- " ax.set_yticklabels(range(rows))\n",
- " ax.xaxis.tick_top()\n",
- "\n",
- " # Add the text from the matrix to each cell\n",
- " for i in range(rows):\n",
- " for j in range(cols):\n",
- " ax.text(j, i, matrix[i][j],\n",
- " ha='center', va='center', color='black', fontsize=12)\n",
- "\n",
- " # Draw grid lines by using minor ticks\n",
- " ax.set_xticks([x - 0.5 for x in range(cols + 1)], minor=True)\n",
- " ax.set_yticks([y - 0.5 for y in range(rows + 1)], minor=True)\n",
- " ax.grid(which='minor', color='black', linestyle='-', linewidth=1)\n",
- "\n",
- " ax.tick_params(which='major', top=False, left=False)\n",
- "\n",
- " return ax\n",
- "\n",
- "def draw_rectangle(ax, top_left, bottom_right, color, width=6, fill_only=False):\n",
- " start_row, start_col = top_left\n",
- " end_row, end_col = bottom_right\n",
- " rect = mpl.patches.Rectangle(\n",
- " (start_col - 0.5, start_row - 0.5), # lower-left corner\n",
- " (end_col - start_col + 1), # width\n",
- " (end_row - start_row + 1), # height\n",
- " edgecolor=color if not fill_only else 'none',\n",
- " facecolor=color if fill_only else 'none',\n",
- " linewidth=width, alpha=0.3 if fill_only else 1\n",
- " )\n",
- " ax.add_patch(rect)\n",
- "\n",
- "ax = draw_matrix(matrix)\n",
- "\n",
- "a, b = top_left\n",
- "c, d = bottom_right\n",
- "draw_rectangle(ax, top_left, bottom_right, 'k', fill_only=True)\n",
- "draw_rectangle(ax, [0, 0], bottom_right, 'r', width=16)\n",
- "draw_rectangle(ax, bottom_right, bottom_right, 'r', fill_only=True)\n",
- "draw_rectangle(ax, [0, 0], [a-1, d], 'b', width=12)\n",
- "draw_rectangle(ax, [a-1, d], [a-1, d], 'b', fill_only=True)\n",
- "draw_rectangle(ax, [0, 0], [c, b-1], 'g', width=8)\n",
- "draw_rectangle(ax, [c, b-1], [c, b-1], 'g', fill_only=True)\n",
- "draw_rectangle(ax, [0, 0], [a-1, b-1], 'orange', width=6)\n",
- "draw_rectangle(ax, [a-1, b-1], [a-1, b-1], 'orange', fill_only=True)\n",
- "\n",
- "\n",
- "plt.tight_layout()\n",
- "plt.show()"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Let's first compute the sum \"manually\":"
- ]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "#| code-fold: true\n",
- "#| echo: false\n",
- "sum_string_op = ''\n",
- "sum_string_value = ''\n",
- "sum_value = 0\n",
- "for i in range(a, c + 1):\n",
- " for j in range(b, d + 1):\n",
- " sum_value += matrix[i][j]\n",
- " sum_string_value += f' + {matrix[i][j]}'\n",
- " sum_string_op += f' + [{i}, {j}]'\n",
- "\n",
- "sum_string_op = sum_string_op[3:] # remove the first \" + \"\n",
- "sum_string_value = sum_string_value[3:]\n",
- "\n",
- "def print_str_nicely(string, max_line_length=80):\n",
- " for i in range(0, len(string), max_line_length):\n",
- " print(string[i:i+max_line_length])\n",
- "\n",
- "print(f'Sum (\"manual\" approach): {sum_value}')\n",
- "print(\"Summed elements:\")\n",
- "print_str_nicely(sum_string_op)\n",
- "print(\"Summed values:\")\n",
- "print_str_nicely(sum_string_value)"
- ],
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Now, let's compute the sum from the integral image:"
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 2D arrays (matrices), array processing (search, cumulative sum)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "---\n",
+ "title: Programming for Engineers\n",
+ "subtitle: \"Lecture 2 - matrices and array processing\"\n",
+ "author: Radoslav Ĺ koviera\n",
+ "\n",
+ "format:\n",
+ " pdf:\n",
+ " code-fold: false\n",
+ " html:\n",
+ " code-fold: false\n",
+ "jupyter: python3\n",
+ "toc: true\n",
+ "---\n",
+ "\n",
+ "## Prelude: Visualizing (printing) arrays\n",
+ "\n",
+ "For debugging purposes, it is useful to see what is in an array. Smaller arrays are easy to print \"whole\":"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[9 8 7 4 2 1 5 6 0 3]\n"
+ ]
+ }
+ ],
+ "source": [
+ "import numpy as np\n",
+ "a = np.random.permutation(10)\n",
+ "print(a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "You can use the `join` method of strings to make the output nicer. More on that later, when we get to strings."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[9, 8, 7, 4, 2, 1, 5, 6, 0, 3]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def pretty_print_array(a):\n",
+ " print('[' + ', '.join([str(x) for x in a]) + ']')\n",
+ "\n",
+ "pretty_print_array(a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For larger arrays, you need to be creative - the solution will depend on what you need. For example, it might be enough to print the first or last few items. Print the array as rows (Python kinda does that but \"uncontrollably\"):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[13 58 55 46 93 15 11 16 35 27]\n",
+ "[10 43 20 72 88 86 60 29 82 32]\n",
+ "[63 67 31 61 47 74 53 18 40 68]\n",
+ "[ 3 12 99 24 4 36 21 80 69 71]\n",
+ "[73 89 98 91 79 49 70 6 59 23]\n",
+ "[57 76 33 83 64 34 7 48 28 66]\n",
+ "[17 8 65 95 87 75 50 5 81 2]\n",
+ "[54 30 26 84 44 1 92 52 39 37]\n",
+ "[ 0 62 38 78 14 41 94 51 19 96]\n",
+ "[ 9 85 42 22 25 97 77 56 45 90]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def print_as_rows(a, row_length=10):\n",
+ " for i in range(0, len(a), row_length):\n",
+ " print(a[i:i+row_length])\n",
+ "\n",
+ "a = np.random.permutation(100)\n",
+ "print_as_rows(a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "You can even improve it to have nicer output. Again, more on that next time, when we get to strings."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[\n",
+ "74, 49, 38, 26, 58, 13, 66, 55, 59, 29,\n",
+ "93, 31, 5, 45, 8, 36, 57, 51, 65, 81,\n",
+ "97, 75, 22, 91, 47, 14, 84, 11, 0, 94,\n",
+ "23, 72, 98, 24, 17, 48, 30, 70, 87, 77,\n",
+ "16, 69, 39, 63, 6, 85, 42, 21, 27, 33,\n",
+ "64, 18, 61, 2, 88, 46, 78, 83, 1, 80,\n",
+ "37, 43, 25, 62, 12, 95, 68, 71, 89, 4,\n",
+ "76, 92, 79, 40, 44, 54, 53, 35, 99, 3,\n",
+ "73, 86, 60, 82, 28, 32, 19, 41, 90, 56,\n",
+ "15, 52, 50, 9, 10, 20, 67, 7, 96, 34,\n",
+ "]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def pretty_print_as_rows(a, row_length=10):\n",
+ " print('[')\n",
+ " for i in range(0, len(a), row_length):\n",
+ " print(', '.join([str(x) for x in a[i:i+row_length]]), end=',\\n')\n",
+ " print(']')\n",
+ "\n",
+ "a = np.random.permutation(100).tolist()\n",
+ "pretty_print_as_rows(a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Printing matrices:\n",
+ "\n",
+ "(remember to define/run cell with the `pretty_print` function above)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Row-by-row print:\n",
+ "[[17, 10, 6, 2, 3, 5, 19, 8, 18, 5, 3],\n",
+ "[11, 17, 15, 8, 5, 6, 3, 13, 7, 1, 19],\n",
+ "[7, 1, 18, 9, 18, 9, 6, 13, 14, 0, 13],\n",
+ "[2, 12, 1, 14, 5, 3, 1, 8, 19, 17, 16],\n",
+ "[19, 17, 0, 10, 19, 9, 4, 14, 15, 12, 7],\n",
+ "[6, 5, 16, 8, 1, 2, 4, 15, 9, 4, 0],\n",
+ "[19, 8, 15, 12, 11, 11, 4, 19, 18, 15, 17],\n",
+ "[18, 11, 5, 9, 6, 5, 19, 1, 2, 7, 0],\n",
+ "[5, 17, 10, 19, 10, 9, 18, 1, 10, 11, 13],\n",
+ "[6, 17, 7, 4, 7, 3, 16, 17, 19, 19, 15],\n",
+ "[9, 12, 5, 6, 15, 6, 1, 2, 16, 7, 19],\n",
+ "[2, 8, 11, 7, 13, 12, 5, 10, 8, 1, 7],\n",
+ "[3, 16, 3, 7, 17, 16, 18, 8, 0, 10, 0],\n",
+ "[8, 9, 0, 2, 5, 6, 5, 5, 18, 13, 3],\n",
+ "[15, 6, 6, 18, 10, 1, 3, 17, 15, 19, 7],\n",
+ "[5, 0, 10, 14, 16, 0, 8, 18, 6, 15, 7],\n",
+ "[14, 0, 6, 3, 18, 13, 5, 6, 6, 18, 0]]\n",
+ "\n",
+ "Formatted print:\n",
+ "[[17, 10, 6, 2, 3, 5, 19, 8, 18, 5, 3,\n",
+ " [11, 17, 15, 8, 5, 6, 3, 13, 7, 1, 19,\n",
+ " [ 7, 1, 18, 9, 18, 9, 6, 13, 14, 0, 13,\n",
+ " [ 2, 12, 1, 14, 5, 3, 1, 8, 19, 17, 16,\n",
+ " [19, 17, 0, 10, 19, 9, 4, 14, 15, 12, 7,\n",
+ " [ 6, 5, 16, 8, 1, 2, 4, 15, 9, 4, 0,\n",
+ " [19, 8, 15, 12, 11, 11, 4, 19, 18, 15, 17,\n",
+ " [18, 11, 5, 9, 6, 5, 19, 1, 2, 7, 0,\n",
+ " [ 5, 17, 10, 19, 10, 9, 18, 1, 10, 11, 13,\n",
+ " [ 6, 17, 7, 4, 7, 3, 16, 17, 19, 19, 15,\n",
+ " [ 9, 12, 5, 6, 15, 6, 1, 2, 16, 7, 19,\n",
+ " [ 2, 8, 11, 7, 13, 12, 5, 10, 8, 1, 7,\n",
+ " [ 3, 16, 3, 7, 17, 16, 18, 8, 0, 10, 0,\n",
+ " [ 8, 9, 0, 2, 5, 6, 5, 5, 18, 13, 3,\n",
+ " [15, 6, 6, 18, 10, 1, 3, 17, 15, 19, 7,\n",
+ " [ 5, 0, 10, 14, 16, 0, 8, 18, 6, 15, 7,\n",
+ " [14, 0, 6, 3, 18, 13, 5, 6, 6, 18, 0]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def print_matrix(m):\n",
+ " print('[', end='')\n",
+ " n_rows = len(m)\n",
+ " for ri, row in enumerate(m):\n",
+ " print(row, end=',\\n' if ri < n_rows - 1 else '')\n",
+ " print(']')\n",
+ "\n",
+ "def formatted_print_matrix(m): # works only for integers\n",
+ " n_rows = len(m)\n",
+ " maximum_value = abs(max([max(row) for row in m])) + 1e-6\n",
+ " if maximum_value > 0:\n",
+ " max_digits = int(np.ceil(np.log10(maximum_value)))\n",
+ " else: # all zeros\n",
+ " max_digits = 1\n",
+ " print('[', end='')\n",
+ " for ri, row in enumerate(m):\n",
+ " row_prefix = ('' if ri == 0 else ' ') + '['\n",
+ " row_numbers = ', '.join([f\"{x:{max_digits}d}\" for x in row])\n",
+ " row_end=',\\n' if ri < n_rows - 1 else ']'\n",
+ " print(row_prefix + row_numbers, end=row_end)\n",
+ " print(']')\n",
+ "\n",
+ "mat = np.random.randint(0, 20, size=(17, 11)).tolist()\n",
+ "print(\"Row-by-row print:\")\n",
+ "print_matrix(mat)\n",
+ "print()\n",
+ "print(\"Formatted print:\")\n",
+ "formatted_print_matrix(mat)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Matrices\n",
+ "\n",
+ "Matrices are 2D arrays. They have two dimensions: rows and columns. The number of rows is often called the height and the number of columns is often called the width. This is similar to image, although, images have typically flipped dimensions: x-axis is width (columns) and y-axis is height (rows)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[1, 2, 3,\n",
+ " [4, 5, 6,\n",
+ " [7, 8, 9]]\n",
+ "[[1, 2, 3,\n",
+ " [4, 5, 6,\n",
+ " [7, 8, 9]]\n",
+ "[[1, 2, 3,\n",
+ " [4, 5, 6,\n",
+ " [7, 8, 9]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "m = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]\n",
+ "formatted_print_matrix(m)\n",
+ "also_m = [\n",
+ " [1, 2, 3],\n",
+ " [4, 5, 6],\n",
+ " [7, 8, 9]\n",
+ "]\n",
+ "formatted_print_matrix(also_m)\n",
+ "m_as_well = [[j + i * 3 for j in range(1, 4)] for i in range(3)]\n",
+ "formatted_print_matrix(m_as_well)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Creating a matrix in a loop (the following two functions are equivalent):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[ 1, 2, 3, 4,\n",
+ " [ 5, 6, 7, 8,\n",
+ " [ 9, 10, 11, 12]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def create_counting_matrix(height, width):\n",
+ " return [[j + i * width for j in range(1, width + 1)] for i in range(height)]\n",
+ "\n",
+ "def create_counting_matrix_unpacked(height, width):\n",
+ " mat = []\n",
+ " for i in range(height):\n",
+ " row = []\n",
+ " for j in range(1, width + 1):\n",
+ " row.append(j + i * width)\n",
+ " mat.append(row)\n",
+ " return mat\n",
+ "\n",
+ "formatted_print_matrix(create_counting_matrix(3, 4))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Matrix pre-allocation:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Empty 3 by 4 matrix:\n",
+ "[[ 0, 0, 0, 0,\n",
+ " [ 0, 0, 0, 0,\n",
+ " [ 0, 0, 0, 0]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def create_empty(n_rows, n_cols):\n",
+ " return [[0] * n_cols for _ in range(n_rows)]\n",
+ "M, N = 3, 4\n",
+ "mat = create_empty(M, N)\n",
+ "print(f\"Empty {M} by {N} matrix:\")\n",
+ "formatted_print_matrix(mat)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Filling matrices with values\n",
+ "\n",
+ "Filling matrix with a constant value. *Warning*, the following function might fail in case of \"jagged\" matrices - this is why it's best to only use 2D arrays (matrices) with equal column lengths (second dimension). There is a simple workaround - get `n_col` at each row-loop. Nonetheless, the best approach is to avoid jagged matrices in the first place."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[11, 11, 11,\n",
+ " [11, 11, 11,\n",
+ " [11, 11, 11]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def fill_matrix_with_value(matrix, value):\n",
+ " n_row = len(matrix)\n",
+ " if n_row == 0: # let's check if matrix is empty to avoid some errors\n",
+ " print(\"Matrix is empty!\")\n",
+ " return\n",
+ " n_col = len(matrix[0])\n",
+ " for r in range(n_row):\n",
+ " for c in range(n_col):\n",
+ " matrix[r][c] = value\n",
+ "\n",
+ "\n",
+ "m = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]\n",
+ "fill_matrix_with_value(m, 11)\n",
+ "formatted_print_matrix(m)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Random fill:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[ 93, 17, 4, 37, 100,\n",
+ " [ 27, 73, 50, 82, 20,\n",
+ " [ 79, 77, 77, 43, 16,\n",
+ " [ 90, 84, 91, 10, 27]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "import random\n",
+ "\n",
+ "def fill_matrix_with_randint(matrix, min_value=0, max_value=100):\n",
+ " n_row = len(matrix)\n",
+ " if n_row == 0: # let's check if matrix is empty to avoid some errors\n",
+ " print(\"Matrix is empty!\")\n",
+ " return\n",
+ " n_col = len(matrix[0])\n",
+ "\n",
+ " for r in range(n_row):\n",
+ " for c in range(n_col):\n",
+ " matrix[r][c] = random.randint(min_value, max_value)\n",
+ "\n",
+ "mat = [[1] * 5 for _ in range(4)]\n",
+ "fill_matrix_with_randint(mat)\n",
+ "formatted_print_matrix(mat)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Shuffle matrix elements. There are few ways to do it. Here, we will create \"flattened\" indices, shuffle them and assign values to them from the original matrix."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Original matrix:\n",
+ "[[1, 2, 3,\n",
+ " [4, 5, 6,\n",
+ " [7, 8, 9]]\n",
+ "Shuffled matrix:\n",
+ "[[7, 8, 3,\n",
+ " [2, 1, 6,\n",
+ " [9, 5, 4]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def shuffle_matrix(matrix):\n",
+ " n_row = len(matrix)\n",
+ " m_cols = len(matrix[0])\n",
+ " flat_indices = list(range(n_row * m_cols))\n",
+ " random.shuffle(flat_indices)\n",
+ " shuffled_matrix = create_empty(n_row, m_cols)\n",
+ " for r in range(n_row):\n",
+ " for c in range(m_cols):\n",
+ " flat_index = r * m_cols + c\n",
+ " shuffled_matrix[r][c] = matrix[int(flat_indices[flat_index] / m_cols)][flat_indices[flat_index] % m_cols]\n",
+ "\n",
+ " return shuffled_matrix\n",
+ "\n",
+ "m = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]\n",
+ "sm = shuffle_matrix(m)\n",
+ "print(\"Original matrix:\")\n",
+ "formatted_print_matrix(m)\n",
+ "print(\"Shuffled matrix:\")\n",
+ "formatted_print_matrix(sm)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Matrix operations\n",
+ "\n",
+ "First, let's create some matrices:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Matrix A:\n",
+ "[[19, 2, 5, 10,\n",
+ " [13, 3, 0, 15,\n",
+ " [ 2, 2, 17, 10]]\n",
+ "Matrix B:\n",
+ "[[19, 3, 18, 8,\n",
+ " [ 0, 10, 5, 7,\n",
+ " [ 9, 6, 7, 7]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "mat_A = create_empty(3, 4)\n",
+ "fill_matrix_with_randint(mat_A, 0, 20)\n",
+ "mat_B = create_empty(3, 4)\n",
+ "fill_matrix_with_randint(mat_B, 0, 20)\n",
+ "print(\"Matrix A:\")\n",
+ "formatted_print_matrix(mat_A)\n",
+ "print(\"Matrix B:\")\n",
+ "formatted_print_matrix(mat_B)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Addition\n",
+ "\n",
+ "Compute $C = A + B$. This one is easy, we just need to loop through elements of the two matrices and add them."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Result of addition of two matrices:\n",
+ "[[38, 5, 23, 18,\n",
+ " [13, 13, 5, 22,\n",
+ " [11, 8, 24, 17]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def add_two_matrices(mat_A, mat_B):\n",
+ " n_row = len(mat_A)\n",
+ " n_col = len(mat_A[0])\n",
+ " mat_C = create_empty(n_row, n_col)\n",
+ " for r in range(n_row):\n",
+ " for c in range(n_col):\n",
+ " mat_C[r][c] = mat_A[r][c] + mat_B[r][c]\n",
+ " return mat_C\n",
+ "\n",
+ "mat_C = add_two_matrices(mat_A, mat_B)\n",
+ "print(\"Result of addition of two matrices:\")\n",
+ "formatted_print_matrix(mat_C)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Transposition\n",
+ "\n",
+ "We want to compute $A^T$. We just swap the indices for the output matrix."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Transposed matrix:\n",
+ "[[19, 13, 2,\n",
+ " [ 2, 3, 2,\n",
+ " [ 5, 0, 17,\n",
+ " [10, 15, 10]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def transpose(mat):\n",
+ " n_row = len(mat) # get the matrix dimensions\n",
+ " n_col = len(mat[0])\n",
+ " mat_T = create_empty(n_col, n_row) # create a new matrix\n",
+ " for r in range(n_row):\n",
+ " for c in range(n_col):\n",
+ " mat_T[c][r] = mat[r][c] # swap the row and column indices\n",
+ " return mat_T\n",
+ "\n",
+ "mat_T = transpose(mat_A)\n",
+ "print(\"Transposed matrix:\")\n",
+ "formatted_print_matrix(mat_T)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Multiplication\n",
+ "\n",
+ "Compute $C = A \\times B$. It is not simple, we need to loop through the elements of the matrices and compute their product. The result will be a new matrix with dimensions $n \\times m$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Result of multiplication of two matrices:\n",
+ "[[537, 115, 288,\n",
+ " [376, 135, 240,\n",
+ " [430, 175, 219]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def multiply_two_matrices(mat_A, mat_B):\n",
+ " n_row_A = len(mat_A)\n",
+ " n_col_A = len(mat_A[0])\n",
+ " n_row_B = len(mat_B)\n",
+ " n_col_B = len(mat_B[0])\n",
+ "\n",
+ " if n_col_A != n_row_B: # the matrices need to have the correct shapes\n",
+ " raise ValueError(f\"Matrices A and B cannot be multiplied. \"\n",
+ " f\"Matrix B needs to have the same number of columns (has {n_col_B}) as matrix A has rows (has {n_row_A}).\")\n",
+ " # this function will not be able to deal with broadcasting!\n",
+ " mat_C = create_empty(n_row_A, n_col_B)\n",
+ " for r in range(n_row_A):\n",
+ " for c in range(n_col_B):\n",
+ " mat_C[r][c] = 0\n",
+ " for k in range(n_col_A):\n",
+ " mat_C[r][c] += mat_A[r][k] * mat_B[k][c]\n",
+ " return mat_C\n",
+ "\n",
+ "mat_C = multiply_two_matrices(mat_A, transpose(mat_B))\n",
+ "print(\"Result of multiplication of two matrices:\")\n",
+ "formatted_print_matrix(mat_C)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Miscellaneous operations on matrices\n",
+ "\n",
+ "#### Comparison\n",
+ "\n",
+ "Compare two matrices \"plain\" - just iterate over the rows & columns and compare corresponding elements in both matrices. Sometimes, we might want to compare elements in two matrices, regardless of their position. In that case, it is easier to \"flatten\" one of the matrices - it is easier to iterate over a 1D array. The 'unique' version of the second function requires us to remove duplicates from the flattened matrix first (you can remove them from any matrix but again, it is easier to do it for a 1D array). Otherwise, we would need to remove them at the end, when comparing the values."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Matrix A:\n",
+ "[[ 6, 10, 4, 9,\n",
+ " [ 7, 9, 5, 8,\n",
+ " [ 6, 9, 7, 6]]\n",
+ "Matrix B:\n",
+ "[[ 3, 10, 4, 8,\n",
+ " [ 6, 5, 6, 5,\n",
+ " [ 6, 8, 7, 9]]\n",
+ "Matrices are equal: False\n",
+ "Equal values: [6, 10, 4, 9, 7, 9, 5, 8, 6, 9, 7, 6]\n",
+ "Uniquely equal values: [6, 10, 4, 9, 7, 5, 8]\n"
+ ]
+ }
+ ],
+ "source": [
+ "def compare_two_matrices(mat_A, mat_B):\n",
+ " n_row = len(mat_A)\n",
+ " n_col = len(mat_A[0])\n",
+ " if n_row != len(mat_B) or n_col != len(mat_B[0]):\n",
+ " return False # matrices don't have the same shape => cannot be equal\n",
+ " for r in range(n_row):\n",
+ " for c in range(n_col):\n",
+ " if mat_A[r][c] != mat_B[r][c]:\n",
+ " return False # \"early\" termination for efficiency\n",
+ " return True\n",
+ "\n",
+ "def remove_duplicates(arr):\n",
+ " i = 0\n",
+ " while i < len(arr) - 1:\n",
+ " j = i + 1\n",
+ " while j < len(arr):\n",
+ " if arr[i] == arr[j]:\n",
+ " del arr[j]\n",
+ " j += 1\n",
+ " i += 1\n",
+ "\n",
+ "def find_equal_values(mat_A, mat_B, unique=False):\n",
+ " n_row_A = len(mat_A)\n",
+ " n_col_A = len(mat_A[0])\n",
+ " n_row_B = len(mat_B)\n",
+ " n_col_B = len(mat_B[0])\n",
+ " flat_values_B = []\n",
+ " for r in range(n_row_B):\n",
+ " flat_values_B.extend(mat_B[r])\n",
+ "\n",
+ " if unique: # remove duplicates, otherwise the `remove` method below might not be enough\n",
+ " remove_duplicates(flat_values_B)\n",
+ "\n",
+ " equal_values = []\n",
+ " for r in range(n_row_A):\n",
+ " for c in range(n_col_A):\n",
+ " item_A = mat_A[r][c]\n",
+ " if item_A in flat_values_B:\n",
+ " equal_values.append(item_A)\n",
+ " if unique:\n",
+ " flat_values_B.remove(item_A)\n",
+ " return equal_values\n",
+ "\n",
+ "mat_A = create_empty(3, 4)\n",
+ "fill_matrix_with_randint(mat_A, 0, 10)\n",
+ "mat_B = create_empty(3, 4)\n",
+ "fill_matrix_with_randint(mat_B, 0, 10)\n",
+ "print(\"Matrix A:\")\n",
+ "formatted_print_matrix(mat_A)\n",
+ "print(\"Matrix B:\")\n",
+ "formatted_print_matrix(mat_B)\n",
+ "\n",
+ "print(\"Matrices are equal:\", compare_two_matrices(mat_A, mat_B))\n",
+ "print(\"Equal values:\", find_equal_values(mat_A, mat_B))\n",
+ "print(\"Uniquely equal values:\", find_equal_values(mat_A, mat_B, unique=True))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### All different\n",
+ "\n",
+ "Find whether the matrix contains only unique values. One option is to use the function above. Then, the number of uniquely equal values with itself must equal number of elements."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[50, 48, 26, 2,\n",
+ " [ 9, 44, 39, 41,\n",
+ " [ 4, 39, 19, 11]]\n",
+ "Matrix does not contain only unique values.\n"
+ ]
+ }
+ ],
+ "source": [
+ "# print(\"Matrix contains only unique values:\", len() == 0)\n",
+ "mat_U = create_empty(3, 4)\n",
+ "fill_matrix_with_randint(mat_U, 0, 50)\n",
+ "formatted_print_matrix(mat_U)\n",
+ "num_equal_values = len(find_equal_values(mat_U, mat_U, unique=True))\n",
+ "if num_equal_values == len(mat_U) * len(mat_U[0]):\n",
+ " print(\"Matrix contains only unique values.\")\n",
+ "else:\n",
+ " print(\"Matrix does not contain only unique values.\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Alternatively, a custom method with similar approach can be made. The easiest approach is to flatten the matrix and then basically compare \"two\" arrays."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[34, 6, 45, 24,\n",
+ " [ 3, 19, 16, 8,\n",
+ " [25, 28, 4, 44]]\n",
+ "Matrix contains only unique values.\n"
+ ]
+ }
+ ],
+ "source": [
+ "def all_different(mat):\n",
+ " # flatten the matrix\n",
+ " flat_values = []\n",
+ " for r in range(len(mat)):\n",
+ " flat_values.extend(mat[r])\n",
+ " n = len(flat_values)\n",
+ " for i in range(n - 1):\n",
+ " for j in range(i + 1, n):\n",
+ " if flat_values[i] == flat_values[j]:\n",
+ " return False\n",
+ " return True\n",
+ "\n",
+ "mat_U = create_empty(3, 4)\n",
+ "fill_matrix_with_randint(mat_U, 0, 50)\n",
+ "formatted_print_matrix(mat_U)\n",
+ "if all_different(mat_U):\n",
+ " print(\"Matrix contains only unique values.\")\n",
+ "else:\n",
+ " print(\"Matrix does not contain only unique values.\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### \"Graphical\" matrices (images)\n",
+ "\n",
+ "Raw images are basically matrices, where each element represents an intensity value of the corresponding pixel. Color images are 3D matrices, where each element is a 3-tuple of intensity values of the color channels but we will not go there.\n",
+ "\n",
+ "However, we can have a little fun with gray-scale images."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[100, 100, 100, 100, 100, 100, 100, 100, 100, 100], [100, 0, 0, 0, 0, 0, 0, 0, 0, 100], [0, 0, 128, 128, 0, 0, 128, 128, 0, 0], [0, 0, 255, 255, 0, 0, 255, 255, 0, 0], [0, 0, 255, 255, 0, 0, 255, 255, 0, 0], [0, 0, 128, 128, 0, 0, 128, 128, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 200, 0, 0, 0, 0], [0, 0, 0, 0, 0, 200, 0, 0, 0, 0], [0, 0, 0, 0, 200, 200, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 255, 0, 0, 0, 0, 0, 0, 255, 0], [0, 255, 255, 255, 255, 255, 255, 255, 255, 0], [0, 0, 255, 128, 128, 128, 128, 255, 0, 0], [0, 0, 0, 255, 255, 255, 255, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 100, 100, 0, 0, 0, 0]]\n"
+ ]
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI4AAADnCAYAAADIBm6aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAADFklEQVR4nO3d0Y3aUBBA0ThKURQEdYDrgILoijQQI/nKMQ/vOb/sYlhdjbQaPXt6vV6/YK3fn/4AfCfhkAiHRDgkwiH58+7Fy+XiX64f7H6/T0uvmTgkwiERDolwSIRDIhwS4ZAIh0Q4JMIhEQ6JcEjeLjmLx+Ox9VuygfP5vOn7mTgkwiERDolwSIRDIhwS4ZAIh0Q4JMIhEQ6JcEimd/fHmaZpl5Oct9vtENfY6zp7fZfX6+UkJ9sSDolwSIRDIhwS4ZAIh0Q4JMIhEQ6JcEiG2FXt8TyJaVpcu2zqYN/FroptCYdEOCTCIREOiXBIhEMiHBLhkAiHRDgkwiERDolwSIRDIhwS4ZAIh0Q4JMIhEQ6JcEiEQyIcEuGQbP5oxWKe509/hM0c6bu8Y+KQCIdEOCTCIREOiXBIhEMiHBLhkAiHRDgkwiEZ4na1jMntatmccEiEQyIcEuGQCIdEOCTCIREOiXBIhEMiHJIhTnKO6vl8rv6d0+n0Hz7JeEwcEuGQCIdEOCTCIREOiXBIhEMiHBLhkAiHRDgkwiERDolwSIRDIhwS4ZAIh0Q4JMIhEQ6JcEiEQ/JjDuSVw3UsM3FIhEMiHBLhkAiHRDgkwiERDolwSIRDIhwS4ZB4tCKLPFqRzQmHRDgkwiERDolwSIRDIhwS4ZAIh0Q4JMIh2fwk57ul6ZJpWtyl8Q8j/I1NHBLhkAiHRDgkwiERDolwSIRDIhwS4ZAIh0Q4JEPcrrYs7fgsE4dEOCTCIREOiXBIhEMiHBLhkAiHRDgkwiERDsnmS85yYnDtknOe59XXGNn1el318yOcfDVxSIRDIhwS4ZAIh0Q4JMIhEQ6JcEiEQyIcEuGQDHGSc+3S7mgnP0dYWq5l4pAIh0Q4JMIhEQ6JcEiEQyIcEuGQCIdEOCRD7KrW+sbdztGYOCTCIREOiXBIhEMiHBLhkAiHRDgkwiERDolwSL5yyVmcz+ddrvN4PHa5zqeZOCTCIREOiXBIhEMiHBLhkAiHRDgkwiERDolwSKaj3cGTfZg4JMIhEQ6JcEiEQyIckr+5ClQxu0ZGEwAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
]
- },
- {
- "cell_type": "code",
- "metadata": {},
- "source": [
- "sum_value = ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]\n",
- "print(f'Sum (integral image approach): {sum_value}')"
- ],
- "execution_count": null,
- "outputs": []
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "img = [\n",
+ " [100, 100, 100, 100, 100, 100, 100, 100, 100, 100],\n",
+ " [100, 0 , 0 , 0, 0 , 0 , 0, 0, 0, 100],\n",
+ " [ 0, 0, 128, 128, 0, 0, 128, 128, 0, 0],\n",
+ " [ 0, 0, 255, 255, 0, 0, 255, 255, 0, 0],\n",
+ " [ 0, 0, 255, 255, 0, 0, 255, 255, 0, 0],\n",
+ " [ 0, 0, 128, 128, 0, 0, 128, 128, 0, 0],\n",
+ " [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
+ " [ 0, 0, 0, 0, 0, 200, 0, 0, 0, 0],\n",
+ " [ 0, 0, 0, 0, 0, 200, 0, 0, 0, 0],\n",
+ " [ 0, 0, 0, 0, 200, 200, 0, 0, 0, 0],\n",
+ " [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
+ " [ 0, 255, 0, 0, 0, 0, 0, 0, 255, 0],\n",
+ " [ 0, 255, 255, 255, 255, 255, 255, 255, 255, 0],\n",
+ " [ 0, 0, 255, 128, 128, 128, 128, 255, 0, 0],\n",
+ " [ 0, 0, 0, 255, 255, 255, 255, 0, 0, 0],\n",
+ " [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
+ " [ 0, 0, 0, 0, 100, 100, 0, 0, 0, 0],\n",
+ "]\n",
+ "\n",
+ "print(img)\n",
+ "\n",
+ "from matplotlib import pyplot as plt\n",
+ "\n",
+ "plt.imshow(img, cmap='gray')\n",
+ "plt.axis('off')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Array & matrix processing\n",
+ "\n",
+ "### Search\n",
+ "\n",
+ "#### Linear search\n",
+ "\n",
+ "Requires iterating over the array, asymptotic runtime is $O(n)$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Linear search run time:\n",
+ "2.91 µs ± 188 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
+ ]
}
- ],
- "metadata": {
- "kernelspec": {
- "name": "python3",
- "language": "python",
- "display_name": "Python 3",
- "path": "/usr/share/jupyter/kernels/python3"
+ ],
+ "source": [
+ "import numpy as np\n",
+ "N = 100\n",
+ "a = np.random.permutation(N)\n",
+ "searched_item = 50\n",
+ "\n",
+ "def linear_search(a, x):\n",
+ " for i in range(len(a)):\n",
+ " if a[i] == x:\n",
+ " return i\n",
+ " return None\n",
+ "\n",
+ "print(\"Linear search run time:\")\n",
+ "%timeit linear_search(a, searched_item)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "##### Computational complexity side-quest\n",
+ "\n",
+ "Even though the 'asymptotic' time of linear search is $O(n)$, the 'actual' time depends on the 'properties' of the array - where the item is located in the array."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Searching for an item in the middle of the array:\n",
+ "4.77 µs ± 354 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
+ "Searching for an item at the beginning of the array:\n",
+ "317 ns ± 21.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n",
+ "Searching for an item at the end of the array:\n",
+ "8.87 µs ± 252 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "N = 100\n",
+ "a = np.random.permutation(N)\n",
+ "searched_item = a[N // 2]\n",
+ "print(\"Searching for an item in the middle of the array:\")\n",
+ "%timeit linear_search(a, searched_item)\n",
+ "\n",
+ "searched_item = a[0]\n",
+ "print(\"Searching for an item at the beginning of the array:\")\n",
+ "%timeit linear_search(a, searched_item)\n",
+ "\n",
+ "searched_item = a[-1]\n",
+ "print(\"Searching for an item at the end of the array:\")\n",
+ "%timeit linear_search(a, searched_item)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Binary search\n",
+ "\n",
+ "Requires pre-sorted array but asymptotic runtime is $O(\\log n)$ - we only look which \"half\" of the array might contain the item."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Binary search run time:\n",
+ "1.29 µs ± 21.3 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n",
+ "2.43 µs ± 93.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
+ "Buuuuuut...\n",
+ "'Fair' binary search run time (includes sorting):\n",
+ "2.81 µs ± 177 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
+ "4.09 µs ± 139 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
+ ]
}
+ ],
+ "source": [
+ "def binary_search_loop(a, x):\n",
+ " low = 0 # lower-bound index for search\n",
+ " high = len(a) - 1 # upper-bound index for search\n",
+ " loop_count = 1\n",
+ " while low <= high: # while there are indices to search\n",
+ " mid = (low + high) // 2 # compute the middle index\n",
+ " if a[mid] < x: # if the middle point is less than the item\n",
+ " low = mid + 1 # \"discard\" the left half (search from the middle + 1)\n",
+ " elif a[mid] > x: # if the middle point is greater than the item\n",
+ " high = mid - 1 # \"discard\" the right half (search from the middle - 1)\n",
+ " else: # otherwise, the middle point is the item\n",
+ " return mid, loop_count\n",
+ " loop_count += 1\n",
+ " return None, loop_count # we did not find the item\n",
+ "\n",
+ "def binary_search_recursion(a, x):\n",
+ " if len(a) == 0: # we got an empty array\n",
+ " return None # the item is for sure not in it\n",
+ " mid_idx = len(a) // 2 # compute the middle point index\n",
+ " if a[mid_idx] == x: # if the middle point is the item\n",
+ " return mid_idx # we found it, job done\n",
+ " elif a[mid_idx] < x: # if the middle point is less than the item\n",
+ " # we look at the right half of the array\n",
+ " return binary_search_recursion(a[mid_idx + 1:], x)\n",
+ " else:\n",
+ " # otherwise, we look at the left half of the array\n",
+ " return binary_search_recursion(a[:mid_idx], x)\n",
+ "\n",
+ "# recursion is slower in this case, as we are copying the whole array, instead of just using the index.\n",
+ "sa = np.sort(a)\n",
+ "print(\"Binary search run time:\")\n",
+ "%timeit binary_search_loop(sa, 5)\n",
+ "%timeit binary_search_recursion(sa, 5)\n",
+ "\n",
+ "print(\"Buuuuuut...\")\n",
+ "print(\"'Fair' binary search run time (includes sorting):\")\n",
+ "%timeit binary_search_loop(np.sort(a), 5)\n",
+ "%timeit binary_search_recursion(np.sort(a), 5)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The advantage of binary search is relatively constant search time. As you can see from the following, the search time is shorter than \"average\" search time for linear search: $O(\\log n) < O(n)$. This includes even sorting the array."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Searching for an item in the middle of the array, using binary search:\n",
+ "3.46 µs ± 255 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
+ "Searching for an item at the beginning of the array, using binary search:\n",
+ "3.52 µs ± 295 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
+ "Searching for an item at the end of the array, using binary search:\n",
+ "2.91 µs ± 113 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\"Searching for an item in the middle of the array, using binary search:\")\n",
+ "searched_item = sa[N // 2]\n",
+ "%timeit binary_search_loop(np.sort(a), searched_item)\n",
+ "\n",
+ "print(\"Searching for an item at the beginning of the array, using binary search:\")\n",
+ "searched_item = sa[0]\n",
+ "%timeit binary_search_loop(np.sort(a), searched_item)\n",
+ "\n",
+ "print(\"Searching for an item at the end of the array, using binary search:\")\n",
+ "searched_item = sa[-1]\n",
+ "%timeit binary_search_loop(np.sort(a), searched_item)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Interpolation search\n",
+ "\n",
+ "Even more restrictive requirements than binary search - the array needs to be 'uniformly' (equal distribution of values) sorted. It uses similar approach as binary search but uses \"an educated guess\" of where the item might be (this is where the uniformity comes into play). The is $O(\\log n)$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Searching for an item in the middle of the array, using interpolation search:\n",
+ "2.58 µs ± 88.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
+ "Searching for an item at the beginning of the array, using interpolation search:\n",
+ "2.62 µs ± 76.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
+ "Searching for an item at the end of the array, using interpolation search:\n",
+ "2.84 µs ± 284 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "def interpolation_search(a, x):\n",
+ " low = 0 # lower-bound index for search\n",
+ " high = len(a) - 1 # upper-bound index for search\n",
+ " loop_count = 1\n",
+ " while low <= high and a[low] <= x <= a[high]: # while there are indices to search\n",
+ " mid = low + int(((x - a[low]) * (high - low)) / (a[high] - a[low]))\n",
+ " if a[mid] < x: # if the middle point is less than the item\n",
+ " low = mid + 1 # \"discard\" the left half (search from the middle + 1)\n",
+ " elif a[mid] > x: # if the middle point is greater than the item\n",
+ " high = mid - 1 # \"discard\" the right half (search from the middle - 1)\n",
+ " else: # otherwise, the middle point is the item\n",
+ " return mid, loop_count\n",
+ " loop_count += 1\n",
+ " return None, loop_count # we did not find the item\n",
+ "\n",
+ "print(\"Searching for an item in the middle of the array, using interpolation search:\")\n",
+ "searched_item = sa[N // 2]\n",
+ "%timeit interpolation_search(np.sort(a), searched_item)\n",
+ "\n",
+ "print(\"Searching for an item at the beginning of the array, using interpolation search:\")\n",
+ "searched_item = sa[0]\n",
+ "%timeit interpolation_search(np.sort(a), searched_item)\n",
+ "\n",
+ "print(\"Searching for an item at the end of the array, using interpolation search:\")\n",
+ "searched_item = sa[-1]\n",
+ "%timeit interpolation_search(np.sort(a), searched_item)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The awesome speed of interpolation search:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Interpolation search run time:\n",
+ "1 µs ± 58.9 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n",
+ "Interpolation search 'step' count: 1\n",
+ "Binary search run time:\n",
+ "1.83 µs ± 94.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n",
+ "Binary search 'step' count: 7\n"
+ ]
+ }
+ ],
+ "source": [
+ "a = np.random.permutation(N)\n",
+ "sa = np.sort(a)\n",
+ "searched_item = sa[N // 2 + 1]\n",
+ "print(\"Interpolation search run time:\")\n",
+ "%timeit interpolation_search(sa, searched_item)\n",
+ "v, i_loop_count = interpolation_search(sa, searched_item)\n",
+ "print(f\"Interpolation search 'step' count: {i_loop_count}\")\n",
+ "\n",
+ "print(\"Binary search run time:\")\n",
+ "%timeit binary_search_loop(sa, searched_item)\n",
+ "b, b_loop_count = binary_search_loop(sa, searched_item)\n",
+ "print(f\"Binary search 'step' count: {b_loop_count}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Unfortunately, the \"single step\" search only works if the array is completely uniform."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Interpolation search run time:\n"
+ ]
+ }
+ ],
+ "source": [
+ "a = np.random.choice(N**3, size=N, replace=False) / 100 # possibly non-uniform array\n",
+ "sa = np.sort(a)\n",
+ "searched_item = sa[N // 2 + 1]\n",
+ "\n",
+ "print(\"Interpolation search run time:\")\n",
+ "%timeit interpolation_search(sa, searched_item)\n",
+ "v, i_loop_count = interpolation_search(sa, searched_item)\n",
+ "print(f\"Interpolation search 'step' count: {i_loop_count}\")\n",
+ "print(\"Binary search run time:\")\n",
+ "%timeit binary_search_loop(sa, searched_item)\n",
+ "v, b_loop_count = binary_search_loop(sa, searched_item)\n",
+ "print(f\"Binary search 'step' count: {b_loop_count}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Find all\n",
+ "\n",
+ "The search finds the first occurrence of the searched item. Which one will depend on the search method. Sometimes, it is required to find all occurrences of the searched item in the array."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "def find_all(a, x):\n",
+ " indices = []\n",
+ " for i in range(len(a)):\n",
+ " if a[i] == x:\n",
+ " indices.append(i)\n",
+ " return indices\n",
+ "\n",
+ "N = 100\n",
+ "a = np.random.choice(N**3, size=N, replace=True)\n",
+ "searched_item = a[N // 2]\n",
+ "occurrences = find_all(a, searched_item)\n",
+ "print(\"The value\", searched_item, \"occurs\", len(occurrences), \"times in the array at indices\", occurrences)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 'Statistical' computations on arrays\n",
+ "\n",
+ "#### Simple operations\n",
+ "\n",
+ "##### Minimum, maximum\n",
+ "\n",
+ "Python has built-in `min()` and `max()` functions. Here, we wil see simple implementation of these functions, simply to practice working with arrays. In practice, you can use the built-in functions."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "my_list = np.random.permutation(10).tolist()\n",
+ "print(\"Input array:\", my_list)\n",
+ "\n",
+ "def minimum(a):\n",
+ " min_value = a[0]\n",
+ " for i in range(1, len(a)):\n",
+ " if a[i] < min_value:\n",
+ " min_value = a[i]\n",
+ " # Alternatively:\n",
+ " # min_value = min(min_value, a[i])\n",
+ " return min_value\n",
+ "\n",
+ "def maximum(a):\n",
+ " max_value = a[0]\n",
+ " for i in range(1, len(a)):\n",
+ " if a[i] > max_value:\n",
+ " max_value = a[i]\n",
+ " # Alternatively:\n",
+ " # max_value = max(max_value, a[i])\n",
+ " return max_value\n",
+ "\n",
+ "print(\"Minimum: \", minimum(my_list))\n",
+ "print(\"Maximum: \", maximum(my_list))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "##### Mean\n",
+ "\n",
+ "First, we need to compute the sum:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "my_list = np.random.permutation(10).tolist()\n",
+ "print(\"Input array:\", my_list)\n",
+ "\n",
+ "def sum(a):\n",
+ " s = 0\n",
+ " for i in range(len(a)):\n",
+ " s += a[i]\n",
+ " return s\n",
+ "\n",
+ "print(\"Sum: \", sum(my_list))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Mean is then simply the sum divided by the number of items:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def mean(a):\n",
+ " return sum(a) / len(a)\n",
+ "\n",
+ "print(\"Mean: \", mean(my_list))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Cumulative sum (prefix sum)\n",
+ "\n",
+ "Summing items between two indices is useful for many tasks. For example, computing average temperature between two specified dates. The following method computes sum between two indices in an array (inclusive of the start and the end indices)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "temp_measurements = (np.random.rand(200) * 10).tolist()\n",
+ "print(\"Input array: \" + ', '.join([f\"{x:.2f}\" for x in temp_measurements[:10]]) + \", ...\")\n",
+ "\n",
+ "def range_sum(a, start, end):\n",
+ " s = 0\n",
+ " for i in range(start, end + 1): # +1 to include the last element\n",
+ " s += a[i]\n",
+ " return s\n",
+ "\n",
+ "print(\"Sum between indices 5 and 15:\", range_sum(temp_measurements, 5, 15))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This is fine, if we need to to this once, but what if we need to do it many times?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def generate_range_queries(a, n, query_len):\n",
+ " numel_a = len(a)\n",
+ " assert query_len <= numel_a, \"Query length must be less than the array length\"\n",
+ " max_pos = numel_a - query_len\n",
+ " queries = []\n",
+ " for _ in range(n):\n",
+ " start = np.random.randint(0, max_pos)\n",
+ " end = start + query_len - 1\n",
+ " queries.append((start, end))\n",
+ " return queries\n",
+ "\n",
+ "print(\"Summing once run time:\")\n",
+ "%timeit range_sum(temp_measurements, 50, 150)\n",
+ "print(\"Summing 100 times run time:\")\n",
+ "queries = generate_range_queries(temp_measurements, 100, 50)\n",
+ "%timeit [range_sum(temp_measurements, start, end) for start, end in queries]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Is there a better way? Why yes! We will use **cumulative sum** (also known as prefix sum). It is done bone adding all the previous values to the current value in the array."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def cumulative_sum(a):\n",
+ " cs = [0] * len(a) # pre-allocate\n",
+ " for i in range(1, len(cs)):\n",
+ " cs[i] = a[i] + cs[i - 1] # add previous sum to the current element\n",
+ " return cs\n",
+ "my_list = list(range(10))\n",
+ "print(\"Cumulative sum of ordered numbers from 0 to 9: \")\n",
+ "print(\"Input array: [\" + ', '.join([f\"{x:2d}\" for x in my_list]) + \"]\")\n",
+ "print(\"Cumulative sum: \" + str(cumulative_sum(my_list)))\n",
+ "print()\n",
+ "my_list = np.random.permutation(10).tolist()\n",
+ "print(\"Cumulative sum of randomly shuffled numbers: \")\n",
+ "print(\"Input array: [\" + ', '.join([f\"{x:2d}\" for x in my_list]) + \"]\")\n",
+ "print(\"Cumulative sum: \" + str(cumulative_sum(my_list)))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "| index | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |\n",
+ "|--------|-----|-----|-------|---------|-----|-----|-----|-----|-----|-----|\n",
+ "| value | 4 | 2 | 0 | 8 | 1 | 5 | 7 | 9 | 3 | 6 |\n",
+ "| cumsum | 4 | 6 | 6 | 14 | 15 | 20 | 27 | 36 | 39 | 45 |\n",
+ "| | 4 | 4+2 | 4+2+0 | 4+2+0+8 | ... | | | | | |\n",
+ "\n",
+ ": \"Visualizing\" cumulative sum\n",
+ "\n",
+ "If we have the cumulative sum precomputed, to get sum of elements between two indices $a$ and $b$, we can simply compute $cumsum[b] - cumsum[a-1]$. Therefore, instead of $b-a$ additions, we only compute one addition."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "my_list = np.random.permutation(10).tolist()\n",
+ "\n",
+ "def range_sum_cs(a_cumsum, start, end):\n",
+ " return a_cumsum[end] - a_cumsum[start - 1]\n",
+ "\n",
+ "print(my_list)\n",
+ "start, end = 3, 7\n",
+ "print(f\"Sum between indices {start} and {end}:\", range_sum_cs(cumulative_sum(my_list), start, end))\n",
+ "# sanity check with the \"simple\" method:\n",
+ "print(\"This is the same as with the simple `range_sum` method:\", range_sum(my_list, start, end) == range_sum_cs(cumulative_sum(my_list), start, end))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now, let's say we have the cumulative sum precomputed and we want to compute the range sum for 100 different ranges. How long will it take?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(\"Run time of cumulative sum computation:\")\n",
+ "%timeit cumulative_sum(temp_measurements)\n",
+ "\n",
+ "cumsum_measurements = cumulative_sum(temp_measurements)\n",
+ "print(\"Summing 100 times run time:\")\n",
+ "%timeit [range_sum_cs(cumsum_measurements, start, end) for start, end in queries]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Of course, the actual effectiveness of the cumulative / prefix sum depends on the number of queries and the length of the queries vs. the length of the array. For a few short queries in a long array, it might take longer to just compute the cumulative sum than to compute all the queries."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "short_queries = generate_range_queries(temp_measurements, 10, 5)\n",
+ "print(\"Run time of simple range sum:\")\n",
+ "%timeit [range_sum(temp_measurements, start, end) for start, end in short_queries]\n",
+ "print(\"Run time of cumulative sum & range sum:\")\n",
+ "%timeit cumsum_measurements = cumulative_sum(temp_measurements); [range_sum_cs(cumsum_measurements, start, end) for start, end in short_queries]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Remember to always use the right tool for the job!\n",
+ "\n",
+ "Now, we can compute average values over varying ranges:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def range_average(csa, start, end):\n",
+ " return range_sum_cs(csa, start, end) / (end - start + 1)\n",
+ "\n",
+ "a = list(range(10))\n",
+ "csa = cumulative_sum(a)\n",
+ "print(\"Input array:\", a)\n",
+ "print(\"Average value of elements between indices 2 and 5 (inclusive):\", range_average(csa, 2, 5))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Integral image (summed-area table, cumulative sum in 2D)\n",
+ "\n",
+ "In computer vision, a useful extension of the cumulative sum to 2 dimensions is used. It is called the **integral image** (or summed-area table)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "matrix = create_empty(11, 10)\n",
+ "fill_matrix_with_randint(matrix, 0, 255)\n",
+ "\n",
+ "def integral_image(mat):\n",
+ " ii = create_empty(len(mat), len(mat[0]))\n",
+ " for i in range(len(mat)):\n",
+ " for j in range(len(mat[0])):\n",
+ " ii[i][j] = mat[i][j]\n",
+ " if i > 0:\n",
+ " ii[i][j] += ii[i - 1][j]\n",
+ " if j > 0:\n",
+ " ii[i][j] += ii[i][j - 1]\n",
+ " if i > 0 and j > 0:\n",
+ " ii[i][j] -= ii[i - 1][j - 1]\n",
+ " return ii\n",
+ "\n",
+ "ii = integral_image(matrix)\n",
+ "formatted_print_matrix(matrix)\n",
+ "formatted_print_matrix(ii)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now, to compute the sum of a sub-matrix, the same principle as in the cumulative \"range sum\" is applied. The \"only\" difference is that now we are in 2D. Therefore, the computation is slightly more complicated. To compute the sum over the \"area\" (sub-matrix) demarked by the indices $(a, b)$ and $(c, d)$, we need to compute: $ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]$. Let's visualize the integral image and this computation. Let us consider the problem to be defined as follows:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "top_left = (5, 3) # a, b\n",
+ "bottom_right = (7, 6) # c, d"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's break it down: We want to compute the sum of the gray area - area of interest (AOI) in @fig-integral_image_full. To do that, we need to take the \"total sum\" (delimited by red rectangle in @fig-integral_image_full) of the area from $[0, 0]$ to $[c, d]$. That is, the sum from the origin of the image to the bottom_right corner of the AOI. This sum has the value at $ii[c][d]$. Then, we subtract the two areas from origin to the top-right corner of the AOI (delimited by blue rectangle in @fig-integral_image_full) and the area from the origin to the bottom-left corner of the AOI (delimited by green rectangle in @fig-integral_image_full). The sums of these areas are located at $ii[a-1][d]$ and $ii[c][b-1]$. This way, we subtracted twice the sum of the area from the origin to just above the top-left corner of the AOI. We need to add it once back-in. Therefore, the last step is to add the sum of this area (delimited by orange rectangle in @fig-integral_image_full) that is located at $ii[a-1][b-1]$. To recap, the full equation is:\n",
+ "\n",
+ "$$area_sum[[a, b], [c, d]] = ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]$$\n",
+ "\n",
+ "Specifically, in our case:\n",
+ "\n",
+ "$$area_sum[[5, 3], [7, 6]] = ii[7][6] - ii[4][6] - ii[7][3] + ii[4][3]$$\n",
+ "\n",
+ "| general term | current case term | color in @fig-integral_image_full |\n",
+ "|:-----------------|:------------------|:----------------------------------|\n",
+ "| $+ ii[c][d]$ | $+ii[7][6]$ | red |\n",
+ "| $- ii[a-1][d]$ | $-ii[4][6]$ | blue |\n",
+ "| $- ii[c][b-1]$ | $-ii[7][3]$ | green |\n",
+ "| $+ ii[a-1][b-1]$ | $+ii[4][3]$ | orange |"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#| code-fold: true\n",
+ "#| echo: false\n",
+ "#| label: fig-integral_image_full\n",
+ "#| fig-cap: Integral image with colored regions.\n",
+ "import matplotlib as mpl\n",
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "# The following code is just for visualization. No need to understand it.\n",
+ "def draw_matrix(matrix):\n",
+ " rows = len(matrix)\n",
+ " cols = len(matrix[0])\n",
+ "\n",
+ " fig, ax = plt.subplots(figsize=(cols, rows))\n",
+ "\n",
+ " # Create empty white image\n",
+ " ax.imshow([[1]*cols for _ in range(rows)], cmap='gray', vmin=0, vmax=1)\n",
+ "\n",
+ " # Set ticks to be at the center of each cell\n",
+ " ax.set_xticks(range(cols))\n",
+ " ax.set_yticks(range(rows))\n",
+ " ax.set_xticklabels(range(cols))\n",
+ " ax.set_yticklabels(range(rows))\n",
+ " ax.xaxis.tick_top()\n",
+ "\n",
+ " # Add the text from the matrix to each cell\n",
+ " for i in range(rows):\n",
+ " for j in range(cols):\n",
+ " ax.text(j, i, matrix[i][j],\n",
+ " ha='center', va='center', color='black', fontsize=12)\n",
+ "\n",
+ " # Draw grid lines by using minor ticks\n",
+ " ax.set_xticks([x - 0.5 for x in range(cols + 1)], minor=True)\n",
+ " ax.set_yticks([y - 0.5 for y in range(rows + 1)], minor=True)\n",
+ " ax.grid(which='minor', color='black', linestyle='-', linewidth=1)\n",
+ "\n",
+ " ax.tick_params(which='major', top=False, left=False)\n",
+ "\n",
+ " return ax\n",
+ "\n",
+ "def draw_rectangle(ax, top_left, bottom_right, color, width=6, fill_only=False):\n",
+ " start_row, start_col = top_left\n",
+ " end_row, end_col = bottom_right\n",
+ " rect = mpl.patches.Rectangle(\n",
+ " (start_col - 0.5, start_row - 0.5), # lower-left corner\n",
+ " (end_col - start_col + 1), # width\n",
+ " (end_row - start_row + 1), # height\n",
+ " edgecolor=color if not fill_only else 'none',\n",
+ " facecolor=color if fill_only else 'none',\n",
+ " linewidth=width, alpha=0.3 if fill_only else 1\n",
+ " )\n",
+ " ax.add_patch(rect)\n",
+ "\n",
+ "ax = draw_matrix(matrix)\n",
+ "\n",
+ "a, b = top_left\n",
+ "c, d = bottom_right\n",
+ "draw_rectangle(ax, top_left, bottom_right, 'k', fill_only=True)\n",
+ "draw_rectangle(ax, [0, 0], bottom_right, 'r', width=16)\n",
+ "draw_rectangle(ax, bottom_right, bottom_right, 'r', fill_only=True)\n",
+ "draw_rectangle(ax, [0, 0], [a-1, d], 'b', width=12)\n",
+ "draw_rectangle(ax, [a-1, d], [a-1, d], 'b', fill_only=True)\n",
+ "draw_rectangle(ax, [0, 0], [c, b-1], 'g', width=8)\n",
+ "draw_rectangle(ax, [c, b-1], [c, b-1], 'g', fill_only=True)\n",
+ "draw_rectangle(ax, [0, 0], [a-1, b-1], 'orange', width=6)\n",
+ "draw_rectangle(ax, [a-1, b-1], [a-1, b-1], 'orange', fill_only=True)\n",
+ "\n",
+ "\n",
+ "plt.tight_layout()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's first compute the sum \"manually\":"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#| code-fold: true\n",
+ "#| echo: false\n",
+ "sum_string_op = ''\n",
+ "sum_string_value = ''\n",
+ "sum_value = 0\n",
+ "for i in range(a, c + 1):\n",
+ " for j in range(b, d + 1):\n",
+ " sum_value += matrix[i][j]\n",
+ " sum_string_value += f' + {matrix[i][j]}'\n",
+ " sum_string_op += f' + [{i}, {j}]'\n",
+ "\n",
+ "sum_string_op = sum_string_op[3:] # remove the first \" + \"\n",
+ "sum_string_value = sum_string_value[3:]\n",
+ "\n",
+ "def print_str_nicely(string, max_line_length=80):\n",
+ " for i in range(0, len(string), max_line_length):\n",
+ " print(string[i:i+max_line_length])\n",
+ "\n",
+ "print(f'Sum (\"manual\" approach): {sum_value}')\n",
+ "print(\"Summed elements:\")\n",
+ "print_str_nicely(sum_string_op)\n",
+ "print(\"Summed values:\")\n",
+ "print_str_nicely(sum_string_value)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now, let's compute the sum from the integral image:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sum_value = ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]\n",
+ "print(f'Sum (integral image approach): {sum_value}')"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
},
- "nbformat": 4,
- "nbformat_minor": 4
-}
\ No newline at end of file
+ "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.10"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/src/pge_lectures/lecture_02/l2_matrices_processing.qmd b/src/pge_lectures/lecture_02/l2_matrices_processing.qmd
index 1aedc45..c1c1d9b 100644
--- a/src/pge_lectures/lecture_02/l2_matrices_processing.qmd
+++ b/src/pge_lectures/lecture_02/l2_matrices_processing.qmd
@@ -192,7 +192,8 @@ def shuffle_matrix(matrix):
for r in range(n_row):
for c in range(m_cols):
flat_index = r * m_cols + c
- shuffled_matrix[r][c] = matrix[int(flat_indices[flat_index] / m_cols)][flat_indices[flat_index] % m_cols]
+ shuffled_matrix[r][c] = matrix[int(flat_indices[flat_index]
+ / m_cols)][flat_indices[flat_index] % m_cols]
return shuffled_matrix
@@ -270,8 +271,9 @@ def multiply_two_matrices(mat_A, mat_B):
n_col_B = len(mat_B[0])
if n_col_A != n_row_B: # the matrices need to have the correct shapes
- raise ValueError(f"Matrices A and B cannot be multiplied. "
- f"Matrix B needs to have the same number of columns (has {n_col_B}) as matrix A has rows (has {n_row_A}).")
+ raise ValueError("Matrices A and B cannot be multiplied. "
+ "Matrix B needs to have the same number of columns "
+ f"(has {n_col_B}) as matrix A has rows (has {n_row_A}).")
# this function will not be able to deal with broadcasting!
mat_C = create_empty(n_row_A, n_col_B)
for r in range(n_row_A):
@@ -323,7 +325,8 @@ def find_equal_values(mat_A, mat_B, unique=False):
for r in range(n_row_B):
flat_values_B.extend(mat_B[r])
- if unique: # remove duplicates, otherwise the `remove` method below might not be enough
+ if unique: # remove duplicates,
+ # otherwise the `remove` method below might not be enough
remove_duplicates(flat_values_B)
equal_values = []
@@ -417,7 +420,7 @@ img = [
[ 0, 0, 0, 0, 100, 100, 0, 0, 0, 0],
]
-print(img)
+print(pretty_print_as_rows(img))
from matplotlib import pyplot as plt
@@ -432,7 +435,7 @@ plt.show()
#### Linear search
-Requires iterating over the array, asymptotic runtime is $O(n)$.
+Requires iterating over the array, until the item with the correct value is found.
```{python}
import numpy as np
@@ -452,7 +455,7 @@ print("Linear search run time:")
##### Computational complexity side-quest
-Even though the 'asymptotic' time of linear search is $O(n)$, the 'actual' time depends on the 'properties' of the array - where the item is located in the array.
+The asymptotic time of the linear search is $O(n)$, the 'actual' time depends on the 'properties' of the array - where the item is located in the array.
```{python}
N = 100
@@ -472,7 +475,7 @@ print("Searching for an item at the end of the array:")
#### Binary search
-Requires pre-sorted array but asymptotic runtime is $O(\log n)$ - we only look which "half" of the array might contain the item.
+Requires pre-sorted array but asymptotic runtime is $O(\log n)$ - we only look which "half" of the array might contain the item, then discard the other half.
```{python}
def binary_search_loop(a, x):
@@ -503,7 +506,8 @@ def binary_search_recursion(a, x):
# otherwise, we look at the left half of the array
return binary_search_recursion(a[:mid_idx], x)
-# recursion is slower in this case, as we are copying the whole array, instead of just using the index.
+# recursion is slower in this case, as we are copying the whole array,
+# instead of just using the index.
sa = np.sort(a)
print("Binary search run time:")
%timeit binary_search_loop(sa, 5)
@@ -518,22 +522,25 @@ print("'Fair' binary search run time (includes sorting):")
The advantage of binary search is relatively constant search time. As you can see from the following, the search time is shorter than "average" search time for linear search: $O(\log n) < O(n)$. This includes even sorting the array.
```{python}
-print("Searching for an item in the middle of the array, using binary search:")
+print("Searching for an item in the middle of the array, "
+"using binary search:")
searched_item = sa[N // 2]
%timeit binary_search_loop(np.sort(a), searched_item)
-print("Searching for an item at the beginning of the array, using binary search:")
+print("Searching for an item at the beginning of the array, "
+"using binary search:")
searched_item = sa[0]
%timeit binary_search_loop(np.sort(a), searched_item)
-print("Searching for an item at the end of the array, using binary search:")
+print("Searching for an item at the end of the array, "
+"using binary search:")
searched_item = sa[-1]
%timeit binary_search_loop(np.sort(a), searched_item)
```
#### Interpolation search
-Even more restrictive requirements than binary search - the array needs to be 'uniformly' (equal distribution of values) sorted. It uses similar approach as binary search but uses "an educated guess" of where the item might be (this is where the uniformity comes into play). The is $O(\log n)$
+Even more restrictive requirements than binary search - the array needs to be 'uniformly' (equal distribution of values) sorted. It uses similar approach as binary search but uses "an educated guess" of where the item might be (this is where the uniformity comes into play). The asymptotic run time is $O(\log(log n))$.
```{python}
def interpolation_search(a, x):
@@ -615,7 +622,8 @@ N = 100
a = np.random.choice(N**3, size=N, replace=True)
searched_item = a[N // 2]
occurrences = find_all(a, searched_item)
-print("The value", searched_item, "occurs", len(occurrences), "times in the array at indices", occurrences)
+print("The value", searched_item, "occurs", len(occurrences),
+ "times in the array at indices", occurrences)
```
### 'Statistical' computations on arrays
@@ -687,7 +695,8 @@ Summing items between two indices is useful for many tasks. For example, computi
```{python}
import numpy as np
temp_measurements = (np.random.rand(200) * 10).tolist()
-print("Input array: " + ', '.join([f"{x:.2f}" for x in temp_measurements[:10]]) + ", ...")
+print("Input array: " + ', '.join([f"{x:.2f}"
+ for x in temp_measurements[:10]]) + ", ...")
def range_sum(a, start, end):
s = 0
@@ -756,9 +765,11 @@ def range_sum_cs(a_cumsum, start, end):
print(my_list)
start, end = 3, 7
-print(f"Sum between indices {start} and {end}:", range_sum_cs(cumulative_sum(my_list), start, end))
+print(f"Sum between indices {start} and {end}:",
+ range_sum_cs(cumulative_sum(my_list), start, end))
# sanity check with the "simple" method:
-print("This is the same as with the simple `range_sum` method:", range_sum(my_list, start, end) == range_sum_cs(cumulative_sum(my_list), start, end))
+print("This is the same as with the simple `range_sum` method:",
+ range_sum(my_list, start, end) == range_sum_cs(cumulative_sum(my_list), start, end))
```
Now, let's say we have the cumulative sum precomputed and we want to compute the range sum for 100 different ranges. How long will it take?
@@ -793,7 +804,8 @@ def range_average(csa, start, end):
a = list(range(10))
csa = cumulative_sum(a)
print("Input array:", a)
-print("Average value of elements between indices 2 and 5 (inclusive):", range_average(csa, 2, 5))
+print("Average value of elements between indices 2 and 5 (inclusive):",
+ range_average(csa, 2, 5))
```
### Integral image (summed-area table, cumulative sum in 2D)
@@ -831,11 +843,11 @@ bottom_right = (7, 6) # c, d
Let's break it down: We want to compute the sum of the gray area - area of interest (AOI) in @fig-integral_image_full. To do that, we need to take the "total sum" (delimited by red rectangle in @fig-integral_image_full) of the area from $[0, 0]$ to $[c, d]$. That is, the sum from the origin of the image to the bottom_right corner of the AOI. This sum has the value at $ii[c][d]$. Then, we subtract the two areas from origin to the top-right corner of the AOI (delimited by blue rectangle in @fig-integral_image_full) and the area from the origin to the bottom-left corner of the AOI (delimited by green rectangle in @fig-integral_image_full). The sums of these areas are located at $ii[a-1][d]$ and $ii[c][b-1]$. This way, we subtracted twice the sum of the area from the origin to just above the top-left corner of the AOI. We need to add it once back-in. Therefore, the last step is to add the sum of this area (delimited by orange rectangle in @fig-integral_image_full) that is located at $ii[a-1][b-1]$. To recap, the full equation is:
-$$area_sum[[a, b], [c, d]] = ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]$$
+$$area\_sum[[a, b], [c, d]] = ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]$$
Specifically, in our case:
-$$area_sum[[5, 3], [7, 6]] = ii[7][6] - ii[4][6] - ii[7][3] + ii[4][3]$$
+$$area\_sum[[5, 3], [7, 6]] = ii[7][6] - ii[4][6] - ii[7][3] + ii[4][3]$$
| general term | current case term | color in @fig-integral_image_full |
|:-----------------|:------------------|:----------------------------------|
@@ -897,16 +909,16 @@ def draw_rectangle(ax, top_left, bottom_right, color, width=6, fill_only=False):
)
ax.add_patch(rect)
-ax = draw_matrix(matrix)
+ax = draw_matrix(ii)
a, b = top_left
c, d = bottom_right
draw_rectangle(ax, top_left, bottom_right, 'k', fill_only=True)
-draw_rectangle(ax, [0, 0], bottom_right, 'r', width=16)
+draw_rectangle(ax, [0, 0], bottom_right, 'r', width=22)
draw_rectangle(ax, bottom_right, bottom_right, 'r', fill_only=True)
-draw_rectangle(ax, [0, 0], [a-1, d], 'b', width=12)
+draw_rectangle(ax, [0, 0], [a-1, d], 'b', width=16)
draw_rectangle(ax, [a-1, d], [a-1, d], 'b', fill_only=True)
-draw_rectangle(ax, [0, 0], [c, b-1], 'g', width=8)
+draw_rectangle(ax, [0, 0], [c, b-1], 'g', width=10)
draw_rectangle(ax, [c, b-1], [c, b-1], 'g', fill_only=True)
draw_rectangle(ax, [0, 0], [a-1, b-1], 'orange', width=6)
draw_rectangle(ax, [a-1, b-1], [a-1, b-1], 'orange', fill_only=True)
@@ -916,7 +928,10 @@ plt.tight_layout()
plt.show()
```
+{{< pagebreak >}}
+
Let's first compute the sum "manually":
+
```{python}
#| code-fold: true
#| echo: false
@@ -936,15 +951,19 @@ def print_str_nicely(string, max_line_length=80):
for i in range(0, len(string), max_line_length):
print(string[i:i+max_line_length])
-print(f'Sum ("manual" approach): {sum_value}')
print("Summed elements:")
print_str_nicely(sum_string_op)
print("Summed values:")
print_str_nicely(sum_string_value)
+print(f'Sum ("manual" approach): {sum_value}')
```
Now, let's compute the sum from the integral image:
+
```{python}
sum_value = ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]
-print(f'Sum (integral image approach): {sum_value}')
+
+print("Sum (integral image approach):\n\t"
+ f"{ii[c][d]} - {ii[a-1][d]} - {ii[c][b-1]} + {ii[a-1][b-1]}"
+ f" = {sum_value}")
```
\ No newline at end of file
--
GitLab