Shrinking the Matrix and Clusters in Python

Code along video

Learning objectives

In this code along, we will write the remaining subroutines of upgma() in functions.py: finding the minimum element of the matrix, adding a row and column for a new cluster, and deleting the rows, columns, and clusters that are no longer needed.

Code along summary

Setup

We will continue working in the upgma directory from the previous code along.

Finding the minimum element

Our find_min_element() function returns the row index, column index, and value of the minimum element of the matrix. We first make sure that the matrix has at least two rows and columns. Then we initialize our three variables to the element at row 0 and column 1, which we know exists, and we scan over the values of the matrix; any time we find a smaller value, we update our variables. Note that we can update all three variables in one line.

Our first attempt uses a nested for loop to range over every row r and every column c in that row.

def find_min_element(mtx: DistanceMatrix) -> tuple[int, int, float]:
    """
    Scan the upper triangle of mtx for the minimum value.

    Args:
        mtx (DistanceMatrix): A square distance matrix (at least 2x2).

    Returns:
        tuple[int, int, float]: (row, col, min_val) where row < col.

    Raises:
        ValueError: If the matrix has only one row or column.
    """
    if len(mtx) <= 1 or len(mtx[0]) == 1:
        raise ValueError("Matrix isn't large enough.")

    row = 0
    col = 1
    min_val = mtx[row][col]

    # range over values of the matrix, and any time I find a smaller value, I update these values
    for r in range(len(mtx)):
        for c in range(len(mtx[r])):
            if mtx[r][c] < min_val:
                row, col, min_val = r, c, mtx[r][c]

    return row, col, min_val

Every value on the main diagonal of a distance matrix is 0, and so our function would always return 0 as the minimum. A distance matrix is symmetric, so we only need to look at the values above the main diagonal (or below it; it does not matter which, as long as we are consistent). We will look in the upper right triangle, where the column index is greater than the row index, which is why the docstring promises that row < col. To do so, we start c at r+1.

for r in range(len(mtx)):
    for c in range(r+1, len(mtx[r])):

There is also a more minor issue. The final row of the matrix contains no values in the upper right triangle, so r only needs to range up to len(mtx)-1. Here is our final function.

def find_min_element(mtx: DistanceMatrix) -> tuple[int, int, float]:
    """
    Scan the upper triangle of mtx for the minimum value.

    Args:
        mtx (DistanceMatrix): A square distance matrix (at least 2x2).

    Returns:
        tuple[int, int, float]: (row, col, min_val) where row < col.

    Raises:
        ValueError: If the matrix has only one row or column.
    """
    if len(mtx) <= 1 or len(mtx[0]) == 1:
        raise ValueError("Matrix isn't large enough.")

    row = 0
    col = 1
    min_val = mtx[row][col]

    # range over values of the matrix, and any time I find a smaller value, I update these values
    for r in range(len(mtx)-1):
        for c in range(r+1, len(mtx[r])):
            if mtx[r][c] < min_val:
                row, col, min_val = r, c, mtx[r][c]

    # nice exercise: write a function that asserts that mtx is a distance matrix (non-negative values, symmetric across main diagonal, have the triangle ineq.)

    return row, col, min_val

Adding a row and column

Next, we write add_row_col(), which adds a row and column to the matrix for the cluster that we just created by merging the clusters at indices row and col. This is where the weighted average from our walkthrough of UPGMA comes in. If the merged clusters have sizes s1 and s2, then the distance from the new cluster to cluster k is (s1*D(k, row) + s2*D(k, col))/(s1 + s2). Let’s add this formula to the docstring as a guide for ourselves.

def add_row_col(row: int, col: int, cluster_size1: int, cluster_size2: int, mtx: DistanceMatrix) -> DistanceMatrix:
    """
    Append a new row and column for the merged cluster.

    Uses the idea that if we merge clusters with sizes s1, s2 corresponding to indices row/col of the matrix, then the distance from new cluster to cluster k is

    (s1 * D(k, row) + s2*D(k, col))/(s1 + s2)

    Args:
        row (int): Index of the first cluster to merge (row < col).
        col (int): Index of the second cluster to merge.
        cluster_size1 (int): Leaf count of the first cluster.
        cluster_size2 (int): Leaf count of the second cluster.
        mtx (DistanceMatrix): The current distance matrix.

    Returns:
        DistanceMatrix: The matrix with one new row and column added.
    """

Rather than thinking about how to grow the matrix all at once, we first make a new row of zeros, which needs one more element than the current rows. We do not need to set the final value of the new row, which is the distance from the new cluster to itself and is already 0, so we range over all but the last index. We also do not need to set new_row[row] or new_row[col], since those are the clusters that are going away.

For every other index r, we apply our formula. What are the two distances that we are averaging? In our four-species example, the first value of the new row (3.5) was the average of the distances from row 0 to rows 2 and 3, where 2 and 3 were the row and column of the minimum element. So the two distances are mtx[r][row] and mtx[r][col].

num_rows = len(mtx)

# idea: make a new row and then I'll fill in the other values
new_row = [0.0]*(num_rows+1)  #we'll worry about the extra column values later
for r in range(len(new_row)-1):
    # I don't need to set new_row[row] or new_row[col]
    if r != row and r != col:
        new_row[r] = (cluster_size1 * mtx[r][row] + cluster_size2 * mtx[r][col])/(cluster_size1 + cluster_size2)

Now we do two things. We append the new row to the matrix, which leaves the matrix oddly shaped, since every other row is missing its last value. Then we use symmetry to append new_row[r] to the end of each original row r. (We could also access this value as mtx[num_rows][r], but new_row[r] is clearer.) Note that num_rows is still the original number of rows, which is exactly the number of rows that need a new value.

# do two things: append this row to matrix and add the extra values with symmetry
mtx.append(new_row)

# now set the missing values due to symmetry (don't need to set bottom right)
for r in range(num_rows):  # num_rows = original n
    mtx[r].append(new_row[r])

return mtx

Deleting clusters

Because of how we defined clusters, the row and column indices of the minimum element are exactly the indices of the two clusters that we want to remove. We can remove an element of a list with Python’s built-in del keyword, so our function seems easy to write.

def delete_clusters(clusters: list[Node], row: int, col: int) -> list[Node]:
    """
    Remove two merged clusters from the cluster list.

    Args:
        clusters (list[Node]): The active cluster list.
        row (int): Index of the first cluster to remove (row < col).
        col (int): Index of the second cluster to remove.

    Returns:
        list[Node]: The updated cluster list.
    """
    del clusters[row]
    del clusters[col]

    return clusters

To see the problem, let’s run a quick demonstration at the end of main(). Say that we have a list containing 1 through 6 and that we don’t like even numbers, so we delete the elements at indices 1 and 3.

c = [1, 2, 3, 4, 5, 6]
del c[1]  # c is now [1, 3, 4, 5, 6]
del c[3]  # we removed the 5, not the 4!
print("c is", c)

c is [1, 3, 4, 6]

Think of ejecting the people sitting in seats 1 and 3 of a row in a movie theater. After the first person leaves, everyone else shifts down a seat, so when we delete c[3], we are deleting the element that is now in seat 3, which is 5. The fix is to always delete the higher index first. This is why we insisted that col is bigger than row when we found the minimum element in the upper right triangle of the matrix.

del clusters[col]
del clusters[row]
# always delete the element with higher index first!

return clusters

Deleting rows and columns

Deleting rows and columns from the matrix is the two-dimensional version of the same idea. First, we delete the two rows, col first because it is bigger. Then n – 2 rows remain; we set num_rows to the current length of the matrix so that we do not need to think about this value, and for each remaining row, we delete the elements at indices col and row, again in that order. A statement like del mtx[r][col] might look tricky, but it is the same as deleting an element of a list: it says to look at the list mtx[r] and get rid of whatever is at index col.

def delete_row_col(mtx: DistanceMatrix, row: int, col: int) -> DistanceMatrix:
    """
    Remove two rows and columns from the distance matrix.

    Args:
        mtx (DistanceMatrix): The current distance matrix.
        row (int): Row/column index of the first cluster (row < col).
        col (int): Row/column index of the second cluster.

    Returns:
        DistanceMatrix: The matrix with two rows and columns removed.
    """
    # remember: col > row

    # first, delete 2 rows
    del mtx[col]
    del mtx[row]

    # then, n-2 rows left; range over all of these rows and remove the elements at indices row and col (col first)
    num_rows = len(mtx)   # now n -2
    for r in range(num_rows):
        del mtx[r][col]
        del mtx[r][row]

    return mtx

Looking ahead

Because we planned our code carefully and modularized it, the functions in functions.py have now all snapped into place. But we cannot run UPGMA yet, because upgma() calls the Node method count_leaves(), which counts the number of leaves in the tree beneath a given node. The most natural way of writing this method depends on a concept that we have not yet seen: recursion. We will introduce recursion in the next code along.