Lecture 2 - matrices and array processing

Programming for Engineers

1 2D arrays (matrices), array processing (search, cumulative sum)

1.1 Prelude: Visualizing (printing) arrays

For debugging purposes, it is useful to see what is in an array. Smaller arrays are easy to print “whole”:

import numpy as np
a = np.random.permutation(10).tolist()
print(a)
[0, 3, 7, 1, 4, 9, 5, 6, 8, 2]

You can use the join method of strings to make the output nicer. More on that later, when we get to strings.

def pretty_print_array(a):
    print('[' + ', '.join([str(x) for x in a]) + ']')

pretty_print_array(a)
[0, 3, 7, 1, 4, 9, 5, 6, 8, 2]

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”):

def print_as_rows(a, row_length=10):
    for i in range(0, len(a), row_length):
        print(a[i:i+row_length])

a = np.random.permutation(100).tolist()
print_as_rows(a)
[86, 68, 53, 83, 92, 29, 98, 58, 11, 39]
[17, 95, 13, 84, 48, 71, 54, 81, 72, 99]
[67, 47, 12, 88, 90, 31, 34, 9, 19, 97]
[85, 3, 18, 57, 35, 6, 45, 40, 38, 26]
[91, 43, 20, 94, 51, 27, 23, 77, 21, 63]
[14, 15, 55, 66, 78, 64, 49, 7, 42, 96]
[30, 28, 32, 10, 74, 50, 56, 79, 65, 75]
[46, 52, 62, 16, 70, 24, 2, 5, 69, 8]
[80, 37, 76, 82, 1, 22, 87, 33, 4, 25]
[44, 93, 0, 41, 36, 73, 89, 59, 61, 60]

You can even improve it to have nicer output. Again, more on that next time, when we get to strings.

def pretty_print_as_rows(a, row_length=10):
    print('[')
    for i in range(0, len(a), row_length):
        print(', '.join([str(x) for x in a[i:i+row_length]]), end=',\n')
    print(']')

a = np.random.permutation(100).tolist()
pretty_print_as_rows(a)
[
31, 11, 56, 33, 52, 86, 96, 9, 59, 49,
14, 82, 68, 87, 63, 53, 19, 1, 20, 92,
42, 97, 67, 35, 4, 66, 39, 32, 80, 40,
41, 2, 95, 43, 75, 15, 29, 51, 3, 36,
47, 74, 13, 99, 73, 72, 37, 27, 55, 54,
17, 38, 88, 24, 30, 83, 16, 71, 48, 28,
58, 44, 94, 69, 70, 5, 46, 76, 85, 89,
18, 50, 79, 45, 0, 25, 78, 6, 65, 57,
64, 7, 22, 62, 61, 98, 23, 34, 93, 10,
91, 26, 60, 8, 77, 12, 81, 84, 90, 21,
]

1.1.0.1 Printing matrices:

(remember to define/run cell with the pretty_print function above)

Visualizing a matrix using the build-in print function might not be the best approach:

mat = np.random.randint(0, 20, size=(17, 11)).tolist()
print(mat)
[[9, 8, 16, 1, 6, 8, 6, 19, 1, 11, 2], [19, 12, 18, 7, 0, 8, 19, 5, 12, 16, 5], [19, 10, 7, 14, 16, 8, 8, 4, 1, 11, 15], [8, 4, 16, 3, 10, 3, 5, 6, 11, 8, 8], [18, 6, 0, 3, 13, 0, 15, 5, 13, 15, 7], [8, 4, 19, 9, 12, 17, 8, 6, 17, 7, 4], [9, 17, 6, 0, 6, 11, 1, 3, 12, 12, 16], [4, 3, 4, 0, 1, 10, 10, 8, 7, 3, 15], [17, 8, 2, 1, 13, 8, 18, 13, 19, 6, 10], [12, 16, 2, 4, 4, 19, 13, 19, 5, 18, 9], [18, 17, 12, 8, 15, 18, 2, 9, 10, 3, 2], [17, 15, 6, 3, 14, 6, 11, 12, 14, 4, 12], [6, 0, 19, 4, 17, 10, 7, 1, 12, 1, 13], [16, 11, 19, 11, 13, 15, 11, 1, 12, 17, 15], [16, 15, 7, 6, 13, 12, 6, 4, 4, 19, 6], [8, 2, 13, 3, 9, 4, 19, 6, 1, 11, 14], [19, 5, 18, 15, 6, 19, 16, 13, 17, 11, 9]]
def print_matrix(m):
    print('[', end='')
    n_rows = len(m)
    for ri, row in enumerate(m):
        print(row, end=',\n' if ri < n_rows - 1 else '')
    print(']')

def formatted_print_matrix(m):  # works only for integers
    n_rows = len(m)
    maximum_value = abs(max([max(row) for row in m])) + 1e-6
    if maximum_value > 0:
        max_digits =  int(np.ceil(np.log10(maximum_value)))
    else:  # all zeros
        max_digits = 1
    print('[', end='')
    for ri, row in enumerate(m):
        row_prefix = ('' if ri == 0 else ' ') + '['
        row_numbers = ', '.join([f"{x:{max_digits}d}" for x in row])
        row_end='],\n' if ri < n_rows - 1 else ']'
        print(row_prefix + row_numbers, end=row_end)
    print(']')

mat = np.random.randint(0, 111, size=(17, 11)).tolist()
print("Row-by-row print:")
print_matrix(mat)
print()
print("Formatted print:")
formatted_print_matrix(mat)
Row-by-row print:
[[69, 71, 41, 26, 106, 101, 24, 103, 63, 63, 34],
[49, 23, 95, 106, 39, 25, 73, 89, 0, 72, 41],
[55, 64, 81, 10, 44, 82, 103, 105, 34, 27, 21],
[47, 85, 31, 38, 106, 92, 76, 18, 59, 44, 76],
[79, 50, 86, 25, 71, 66, 96, 41, 19, 66, 97],
[26, 74, 86, 14, 25, 32, 98, 3, 60, 71, 44],
[1, 26, 37, 105, 103, 52, 93, 39, 110, 19, 21],
[52, 15, 102, 37, 17, 91, 29, 36, 34, 73, 44],
[62, 74, 12, 28, 24, 74, 13, 88, 69, 32, 58],
[37, 72, 68, 0, 48, 1, 36, 107, 50, 79, 66],
[11, 103, 59, 7, 84, 42, 32, 53, 66, 85, 75],
[77, 85, 81, 53, 84, 4, 39, 84, 22, 33, 19],
[42, 80, 24, 59, 81, 25, 34, 70, 85, 71, 89],
[91, 6, 52, 59, 43, 72, 73, 77, 6, 88, 58],
[109, 98, 102, 93, 52, 40, 28, 25, 3, 3, 100],
[58, 57, 75, 79, 93, 70, 32, 12, 110, 10, 33],
[78, 17, 52, 62, 42, 25, 93, 12, 10, 74, 15]]

Formatted print:
[[ 69,  71,  41,  26, 106, 101,  24, 103,  63,  63,  34],
 [ 49,  23,  95, 106,  39,  25,  73,  89,   0,  72,  41],
 [ 55,  64,  81,  10,  44,  82, 103, 105,  34,  27,  21],
 [ 47,  85,  31,  38, 106,  92,  76,  18,  59,  44,  76],
 [ 79,  50,  86,  25,  71,  66,  96,  41,  19,  66,  97],
 [ 26,  74,  86,  14,  25,  32,  98,   3,  60,  71,  44],
 [  1,  26,  37, 105, 103,  52,  93,  39, 110,  19,  21],
 [ 52,  15, 102,  37,  17,  91,  29,  36,  34,  73,  44],
 [ 62,  74,  12,  28,  24,  74,  13,  88,  69,  32,  58],
 [ 37,  72,  68,   0,  48,   1,  36, 107,  50,  79,  66],
 [ 11, 103,  59,   7,  84,  42,  32,  53,  66,  85,  75],
 [ 77,  85,  81,  53,  84,   4,  39,  84,  22,  33,  19],
 [ 42,  80,  24,  59,  81,  25,  34,  70,  85,  71,  89],
 [ 91,   6,  52,  59,  43,  72,  73,  77,   6,  88,  58],
 [109,  98, 102,  93,  52,  40,  28,  25,   3,   3, 100],
 [ 58,  57,  75,  79,  93,  70,  32,  12, 110,  10,  33],
 [ 78,  17,  52,  62,  42,  25,  93,  12,  10,  74,  15]]

1.2 Matrices

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).

m = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
formatted_print_matrix(m)
also_m = [
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
]
formatted_print_matrix(also_m)
m_as_well = [[j + i * 3 for j in range(1, 4)] for i in range(3)]
formatted_print_matrix(m_as_well)
[[1, 2, 3],
 [4, 5, 6],
 [7, 8, 9]]
[[1, 2, 3],
 [4, 5, 6],
 [7, 8, 9]]
[[1, 2, 3],
 [4, 5, 6],
 [7, 8, 9]]

Creating a matrix in a loop (the following two functions are equivalent):

def create_counting_matrix(height, width):
    return [[j + i * width for j in range(1, width + 1)] for i in range(height)]

def create_counting_matrix_unpacked(height, width):
    mat = []
    for i in range(height):
        row = []
        for j in range(1, width + 1):
            row.append(j + i * width)
        mat.append(row)
    return mat

formatted_print_matrix(create_counting_matrix(3, 4))
[[ 1,  2,  3,  4],
 [ 5,  6,  7,  8],
 [ 9, 10, 11, 12]]

Matrix pre-allocation:

def create_empty(n_rows, n_cols):
    return [[0] * n_cols for _ in range(n_rows)]
M, N = 3, 4
mat = create_empty(M, N)
print(f"Empty {M} by {N} matrix:")
formatted_print_matrix(mat)
Empty 3 by 4 matrix:
[[     0,      0,      0,      0],
 [     0,      0,      0,      0],
 [     0,      0,      0,      0]]

1.2.1 Filling matrices with values

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.

def fill_matrix_with_value(matrix, value):
    n_row = len(matrix)
    if n_row == 0:  # let's check if matrix is empty to avoid some errors
        print("Matrix is empty!")
        return
    n_col = len(matrix[0])
    for r in range(n_row):
        for c in range(n_col):
            matrix[r][c] = value


m = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
fill_matrix_with_value(m, 11)
formatted_print_matrix(m)
[[11, 11, 11],
 [11, 11, 11],
 [11, 11, 11]]

Random fill:

import random

def fill_matrix_with_randint(matrix, min_value=0, max_value=100):
    n_row = len(matrix)
    if n_row == 0:  # let's check if matrix is empty to avoid some errors
        print("Matrix is empty!")
        return
    n_col = len(matrix[0])

    for r in range(n_row):
        for c in range(n_col):
            matrix[r][c] = random.randint(min_value, max_value)

mat = [[1] * 5 for _ in range(4)]
fill_matrix_with_randint(mat)
formatted_print_matrix(mat)
[[30, 54, 32, 79, 20],
 [75, 94, 35, 72,  0],
 [39, 94, 12, 65, 85],
 [67, 89,  6, 47, 37]]

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.

def shuffle_matrix(matrix):
    n_row = len(matrix)
    m_cols = len(matrix[0])
    flat_indices = list(range(n_row * m_cols))  # compute flat indices
    random.shuffle(flat_indices)  # shuffle the flat indices
    shuffled_matrix = create_empty(n_row, m_cols)  # preallocate empty matrix
    for r in range(n_row):
        for c in range(m_cols):
            flat_index = r * m_cols + c
            new_ind = flat_indices[flat_index]  # get shuffled index
            # '//' is integer division - indices must be integers
            shuffled_matrix[r][c] = matrix[new_ind // m_cols][new_ind % m_cols]

    return shuffled_matrix

m = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
sm = shuffle_matrix(m)
print("Original matrix:")
formatted_print_matrix(m)
print("Shuffled matrix:")
formatted_print_matrix(sm)
Original matrix:
[[1, 2, 3],
 [4, 5, 6],
 [7, 8, 9]]
Shuffled matrix:
[[7, 2, 6],
 [8, 1, 3],
 [4, 9, 5]]

1.2.2 Matrix operations

First, let’s create some matrices:

mat_A = create_empty(3, 4)
fill_matrix_with_randint(mat_A, 0, 20)
mat_B = create_empty(3, 4)
fill_matrix_with_randint(mat_B, 0, 20)
print("Matrix A:")
formatted_print_matrix(mat_A)
print("Matrix B:")
formatted_print_matrix(mat_B)
Matrix A:
[[10,  5,  7,  6],
 [18, 18,  5, 18],
 [ 2, 17, 12, 12]]
Matrix B:
[[ 1,  8,  1,  4],
 [16, 11,  4,  0],
 [ 0, 12, 15,  0]]

1.2.2.1 Addition

Compute \(C = A + B\). This one is easy, we just need to loop through elements of the two matrices and add them.

def add_two_matrices(mat_A, mat_B):
    n_row = len(mat_A)
    n_col = len(mat_A[0])
    mat_C = create_empty(n_row, n_col)
    for r in range(n_row):
        for c in range(n_col):
            mat_C[r][c] = mat_A[r][c] + mat_B[r][c]
    return mat_C

mat_C = add_two_matrices(mat_A, mat_B)
print("Result of addition of two matrices:")
formatted_print_matrix(mat_C)
Result of addition of two matrices:
[[11, 13,  8, 10],
 [34, 29,  9, 18],
 [ 2, 29, 27, 12]]

1.2.2.2 Transposition

We want to compute \(A^T\). We just swap the indices for the output matrix.

def transpose(mat):
    n_row = len(mat)  # get the matrix dimensions
    n_col = len(mat[0])
    mat_T = create_empty(n_col, n_row)  # create a new matrix
    for r in range(n_row):
        for c in range(n_col):
            mat_T[c][r] = mat[r][c]  # swap the row and column indices
    return mat_T

print("Original matrix:")
formatted_print_matrix(mat_A)
mat_T = transpose(mat_A)
print("Transposed matrix:")
formatted_print_matrix(mat_T)
Original matrix:
[[10,  5,  7,  6],
 [18, 18,  5, 18],
 [ 2, 17, 12, 12]]
Transposed matrix:
[[10, 18,  2],
 [ 5, 18, 17],
 [ 7,  5, 12],
 [ 6, 18, 12]]

1.2.2.3 Multiplication

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\).

def multiply_two_matrices(mat_A, mat_B):
    n_row_A = len(mat_A)
    n_col_A = len(mat_A[0])
    n_row_B = len(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("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):
        for c in range(n_col_B):
            mat_C[r][c] = 0
            for k in range(n_col_A):
                mat_C[r][c] += mat_A[r][k] * mat_B[k][c]
    return mat_C

mat_C = multiply_two_matrices(mat_A, transpose(mat_B))
print("Result of multiplication of two matrices:")
formatted_print_matrix(mat_C)
Result of multiplication of two matrices:
[[ 81, 243, 165],
 [239, 506, 291],
 [198, 267, 384]]

1.2.3 Miscellaneous operations on matrices

1.2.3.1 Comparison

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 find_equal_values function requires us to remove duplicates from the flattened matrix first (you can remove them from the 2D 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.

def compare_two_matrices(mat_A, mat_B):
    n_row = len(mat_A)
    n_col = len(mat_A[0])
    if n_row != len(mat_B) or n_col != len(mat_B[0]):
        return False  # matrices don't have the same shape => cannot be equal
    for r in range(n_row):
        for c in range(n_col):
            if mat_A[r][c] != mat_B[r][c]:
                return False  # "early" termination for efficiency
    return True

def remove_duplicates(arr):
    i = 0
    while i < len(arr) - 1:
        j = i + 1
        while j < len(arr):
            if arr[i] == arr[j]:
                del arr[j]
            j += 1
        i += 1

def find_equal_values(mat_A, mat_B, unique=False):
    n_row_A = len(mat_A)
    n_col_A = len(mat_A[0])
    n_row_B = len(mat_B)
    n_col_B = len(mat_B[0])
    flat_values_B = []
    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
        remove_duplicates(flat_values_B)

    equal_values = []
    for r in range(n_row_A):
        for c in range(n_col_A):
            item_A = mat_A[r][c]
            if item_A in flat_values_B:
                equal_values.append(item_A)
                if unique:
                    flat_values_B.remove(item_A)
    return equal_values

mat_A = create_empty(3, 4)
fill_matrix_with_randint(mat_A, 0, 10)
mat_B = create_empty(3, 4)
fill_matrix_with_randint(mat_B, 0, 10)
print("Matrix A:")
formatted_print_matrix(mat_A)
print("Matrix B:")
formatted_print_matrix(mat_B)

print("Matrices are equal:", compare_two_matrices(mat_A, mat_B))
print("Equal values:", find_equal_values(mat_A, mat_B))
print("Uniquely equal values:", find_equal_values(mat_A, mat_B, unique=True))
Matrix A:
[[1, 6, 8, 0],
 [7, 2, 9, 3],
 [1, 7, 9, 8]]
Matrix B:
[[ 0, 10,  5, 10],
 [ 8,  3,  1,  5],
 [ 9,  7,  3,  2]]
Matrices are equal: False
Equal values: [1, 8, 0, 7, 2, 9, 3, 1, 7, 9, 8]
Uniquely equal values: [1, 8, 0, 7, 2, 9, 3]

Of course, for duplicate removal, set also works:

def remove_duplicates_set(mat):
    # of course, this also works:
    arr = []
    for r in mat:
        arr.extend(r)
    return list(set(arr))

flat_A = []
for r in mat_A:
    flat_A.extend(r)
remove_duplicates(flat_A)
print("Duplicates removed from matrix A:", flat_A)
print("Duplicates removed from matrix A:", remove_duplicates_set(mat_A))

%timeit remove_duplicates(flat_A)
%timeit remove_duplicates_set(mat_A)
Duplicates removed from matrix A: [1, 6, 8, 0, 7, 2, 9, 3]
Duplicates removed from matrix A: [0, 1, 2, 3, 6, 7, 8, 9]
1.84 μs ± 110 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
537 ns ± 10.3 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

1.2.3.2 All different

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.

# print("Matrix contains only unique values:", len() == 0)
mat_U = create_empty(3, 4)
fill_matrix_with_randint(mat_U, 0, 50)
formatted_print_matrix(mat_U)
num_equal_values = len(find_equal_values(mat_U, mat_U, unique=True))
if num_equal_values == len(mat_U) * len(mat_U[0]):
    print("Matrix contains only unique values.")
else:
    print("Matrix does not contain only unique values.")
[[47, 12, 38, 49],
 [40, 19, 33,  5],
 [ 5, 30,  0, 46]]
Matrix does not contain only unique values.

Alternatively, a custom method with similar approach can be made. The easiest approach is to flatten the matrix and then basically compare “two” arrays.

def all_different(mat):
    # flatten the matrix
    flat_values = []
    for r in range(len(mat)):
        flat_values.extend(mat[r])
    n = len(flat_values)
    for i in range(n - 1):
        for j in range(i + 1, n):
            if flat_values[i] == flat_values[j]:
                return False
    return True

mat_U = create_empty(3, 4)
fill_matrix_with_randint(mat_U, 0, 50)
formatted_print_matrix(mat_U)
if all_different(mat_U):
    print("Matrix contains only unique values.")
else:
    print("Matrix does not contain only unique values.")
[[33, 11, 20,  8],
 [43,  6, 42, 23],
 [ 8, 17, 46, 45]]
Matrix does not contain only unique values.

1.2.4 “Graphical” matrices (images)

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.

However, we can have a little fun with gray-scale images.

img = [
    [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],
]

print(pretty_print_as_rows(img))

from matplotlib import pyplot as plt

plt.imshow(img, cmap='gray')
plt.axis('off')
plt.show()
[
[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],
]
None

1.3 Array & matrix processing

1.3.2 ‘Statistical’ computations on arrays

1.3.2.1 Simple operations

1.3.2.1.1 Minimum, maximum

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.

import numpy as np
my_list = np.random.permutation(10).tolist()
print("Input array:", my_list)

def minimum(a):
    min_value = a[0]
    for i in range(1, len(a)):
        if a[i] < min_value:
            min_value = a[i]
        # Alternatively:
        # min_value = min(min_value, a[i])
    return min_value

def maximum(a):
    max_value = a[0]
    for i in range(1, len(a)):
        if a[i] > max_value:
            max_value = a[i]
        # Alternatively:
        # max_value = max(max_value, a[i])
    return max_value

print("Minimum: ", minimum(my_list))
print(minimum(my_list) == min(my_list))
print("Maximum: ", maximum(my_list))
print(maximum(my_list) == max(my_list))
Input array: [6, 0, 8, 7, 2, 1, 3, 5, 9, 4]
Minimum:  0
True
Maximum:  9
True
1.3.2.1.2 Mean

First, we need to compute the sum (there is also a built-in sum() function):

import numpy as np
my_list = np.random.permutation(10).tolist()
print("Input array:", my_list)

def arr_sum(a):
    s = 0
    for i in range(len(a)):
        s += a[i]
    return s

print("Sum: ", arr_sum(my_list))
Input array: [5, 6, 0, 9, 4, 2, 1, 3, 8, 7]
Sum:  45

Mean is then simply the sum divided by the number of items:

def mean(a):
    return arr_sum(a) / len(a)

print("Mean: ", mean(my_list))
Mean:  4.5

1.3.2.2 Cumulative sum (prefix sum)

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).

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]]) + ", ...")

def range_sum(a, start, end):
    s = 0
    for i in range(start, end + 1):  # +1 to include the last element
        s += a[i]
    return s

print("Sum between indices 5 and 15:", range_sum(temp_measurements, 5, 15))
Input array: 8.83, 3.47, 0.39, 7.92, 8.12, 5.39, 8.62, 8.55, 5.43, 0.47, ...
Sum between indices 5 and 15: 58.349988517893195

This is fine, if we need to to this once, but what if we need to do it many times?

def generate_range_queries(a, n, query_len):
    numel_a = len(a)
    assert query_len <= numel_a, "Query length must be less than the array length"
    max_pos = numel_a - query_len
    queries = []
    for _ in range(n):
        start = np.random.randint(0, max_pos)
        end = start + query_len - 1
        queries.append((start, end))
    return queries

print("Summing once run time:")
%timeit range_sum(temp_measurements, 50, 150)
print("Summing 100 times run time:")
queries = generate_range_queries(temp_measurements, 100, 50)
%timeit [range_sum(temp_measurements, start, end) for start, end in queries]
Summing once run time:
2.46 μs ± 56.8 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
Summing 100 times run time:
149 μs ± 6.33 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Is there a better way? Why yes! We will use cumulative sum (also known as prefix sum). It is done by adding all the previous values to the current value in the array.

def cumulative_sum(a):
    cs = [0] * len(a)  # pre-allocate
    for i in range(1, len(cs)):
        cs[i] = a[i] + cs[i - 1]  # add previous sum to the current element
    return cs
my_list = list(range(10))
print("Cumulative sum of ordered numbers from 0 to 9: ")
print("Input array:   [" + ', '.join([f"{x:2d}" for x in my_list]) + "]")
print("Cumulative sum: " + str(cumulative_sum(my_list)))
print()
my_list = np.random.permutation(10).tolist()
print("Cumulative sum of randomly shuffled numbers: ")
print("Input array:   [" + ', '.join([f"{x:2d}" for x in my_list]) + "]")
print("Cumulative sum: " + str(cumulative_sum(my_list)))
Cumulative sum of ordered numbers from 0 to 9: 
Input array:   [ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9]
Cumulative sum: [0, 1, 3, 6, 10, 15, 21, 28, 36, 45]

Cumulative sum of randomly shuffled numbers: 
Input array:   [ 6,  3,  5,  2,  8,  7,  9,  1,  0,  4]
Cumulative sum: [0, 3, 8, 10, 18, 25, 34, 35, 35, 39]
“Visualizing” cumulative sum
index 0 1 2 3 4 5 6 7 8 9
value 4 2 0 8 1 5 7 9 3 6
cumsum 4 6 6 14 15 20 27 36 39 45
4 4+2 4+2+0 4+2+0+8 …

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.

my_list = np.random.permutation(10).tolist()

def range_sum_cs(a_cumsum, start, end):
    return a_cumsum[end] - a_cumsum[start - 1]

print(my_list)
start, end = 3, 7
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))
[9, 4, 8, 6, 5, 7, 2, 3, 0, 1]
Sum between indices 3 and 7: 23
This is the same as with the simple `range_sum` method: True

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?

print("Run time of cumulative sum computation:")
%timeit cumulative_sum(temp_measurements)

cumsum_measurements = cumulative_sum(temp_measurements)
print("Summing 100 times run time:")
%timeit [range_sum_cs(cumsum_measurements, start, end) for start, end in queries]
Run time of cumulative sum computation:
9.52 μs ± 398 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
Summing 100 times run time:
8.05 μs ± 461 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

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.

short_queries = generate_range_queries(temp_measurements, 10, 5)
print("Run time of simple range sum:")
%timeit [range_sum(temp_measurements, start, end) for start, end in short_queries]
print("Run time of cumulative sum & range sum:")
%timeit cumsum_measurements = cumulative_sum(temp_measurements); [range_sum_cs(cumsum_measurements, start, end) for start, end in short_queries]
Run time of simple range sum:
3.04 μs ± 173 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
Run time of cumulative sum & range sum:
10.4 μs ± 366 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Remember to always use the right tool for the job!

Now, we can compute average values over varying ranges:

def range_average(csa, start, end):
    return range_sum_cs(csa, start, end) / (end - start + 1)

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))
Input array: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
Average value of elements between indices 2 and 5 (inclusive): 3.5

1.3.3 Integral image (summed-area table, cumulative sum in 2D)

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).

matrix = create_empty(11, 10)
fill_matrix_with_randint(matrix, 0, 255)

def integral_image(mat):
    ii = create_empty(len(mat), len(mat[0]))
    for i in range(len(mat)):
        for j in range(len(mat[0])):
            ii[i][j] = mat[i][j]
            if i > 0:
                ii[i][j] += ii[i - 1][j]
            if j > 0:
                ii[i][j] += ii[i][j - 1]
            if i > 0 and j > 0:
                ii[i][j] -= ii[i - 1][j - 1]
    return ii

ii = integral_image(matrix)
formatted_print_matrix(matrix)
formatted_print_matrix(ii)
[[120,  45, 118, 170, 217, 142, 105, 218, 196, 113],
 [ 51, 197,  93,  95, 227,  39,  66,  37, 220, 199],
 [172,  24, 221, 175, 239, 226, 219, 211,  11,  93],
 [140,  63, 118, 176, 205,  43, 158,  69, 148, 192],
 [248, 235, 114, 187,  76, 220, 205,  69, 225, 174],
 [156,  30, 178, 117, 134,   5,  90, 140,  86, 135],
 [199, 114, 179,  69,  30, 223,  72, 138, 202, 110],
 [ 15, 108, 255, 181,  25, 231,  95,  75, 142,  57],
 [ 48, 212,  35, 111,   9, 187, 214,  53,  76, 187],
 [  7, 157, 215, 114, 246, 240, 157, 109, 158, 221],
 [ 84, 110,  32, 148, 178, 121, 159,  94,  95, 157]]
[[  120,   165,   283,   453,   670,   812,   917,  1135,  1331,  1444],
 [  171,   413,   624,   889,  1333,  1514,  1685,  1940,  2356,  2668],
 [  343,   609,  1041,  1481,  2164,  2571,  2961,  3427,  3854,  4259],
 [  483,   812,  1362,  1978,  2866,  3316,  3864,  4399,  4974,  5571],
 [  731,  1295,  1959,  2762,  3726,  4396,  5149,  5753,  6553,  7324],
 [  887,  1481,  2323,  3243,  4341,  5016,  5859,  6603,  7489,  8395],
 [ 1086,  1794,  2815,  3804,  4932,  5830,  6745,  7627,  8715,  9731],
 [ 1101,  1917,  3193,  4363,  5516,  6645,  7655,  8612,  9842, 10915],
 [ 1149,  2177,  3488,  4769,  5931,  7247,  8471,  9481, 10787, 12047],
 [ 1156,  2341,  3867,  5262,  6670,  8226,  9607, 10726, 12190, 13671],
 [ 1240,  2535,  4093,  5636,  7222,  8899, 10439, 11652, 13211, 14849]]

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:

top_left = (5, 3)  # a, b
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 Figure 1. To do that, we need to take the “total sum” (delimited by red rectangle in Figure 1) 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 Figure 1) and the area from the origin to the bottom-left corner of the AOI (delimited by green rectangle in Figure 1). 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 Figure 1) 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]\]

Specifically, in our case:

\[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 Figure 1
\(+ ii[c][d]\) \(+ii[7][6]\) red
\(- ii[a-1][d]\) \(-ii[4][6]\) blue
\(- ii[c][b-1]\) \(-ii[7][3]\) green
\(+ ii[a-1][b-1]\) \(+ii[4][3]\) orange
Figure 1: Integral image with colored regions.

Let’s first compute the sum “manually”:

Summed elements:
[5, 3] + [5, 4] + [5, 5] + [5, 6] + [6, 3] + [6, 4] + [6, 5] + [6, 6] + [7, 3] +
 [7, 4] + [7, 5] + [7, 6]
Summed values:
117 + 134 + 5 + 90 + 69 + 30 + 223 + 72 + 181 + 25 + 231 + 95
Sum ("manual" approach): 1272

Now, let’s compute the sum from the integral image:

sum_value = ii[c][d] - ii[a-1][d] - ii[c][b-1] + ii[a-1][b-1]

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}")
Sum (integral image approach):
    7655 - 5149 - 3193 + 1959 = 1272