Clustering

5 problems from Rosalind — Bioinformatics Textbook Track. Press Run on any block to execute it in your browser.

# Rosalind: BA8A — Implement FarthestFirstTraversal
# https://rosalind.info/problems/ba8a/
#
# Given: Integers k and m, then a set of points in m-dimensional space.
# Return: The k centers FarthestFirstTraversal chooses, starting from the first
# point of the data.

let k = 3
let points = [
    [0.0, 0.0], [5.0, 5.0], [0.0, 5.0], [1.0, 1.0],
    [2.0, 2.0], [3.0, 3.0], [1.0, 2.0],
]

fn distance(a, b) {
    sqrt(range(0, len(a)) |> map(|i| (a[i] - b[i]) * (a[i] - b[i])) |> sum())
}

# Each new center is the point furthest from every center chosen so far. That
# makes it deterministic — unlike k-means, which depends on where it starts —
# and it is why this is used to seed clustering rather than to do it: the points
# it picks are the extremes, not the middles.
let centers = [points[0]]
while len(centers) < k {
    let spreads = points |> map(|p| centers |> map(|c| distance(p, c)) |> min())
    centers = push(centers, points[argmax(spreads)])
}

let result = centers |> map(|c| c |> map(|v| str(v)) |> join(" ")) |> join(" / ")

println("Result:   " + result)
println("Expected: 0 0 / 5 5 / 0 5")

fn test_ba8a_farthest_first_traversal() {
    assert len(centers) == k, "BA8A: expected k centers"
    assert centers[0] == points[0], "BA8A: the first point seeds the traversal"
    assert contains(centers, [5.0, 5.0]), "BA8A: 5.0 5.0 should be chosen"
    assert contains(centers, [0.0, 5.0]), "BA8A: 0.0 5.0 should be chosen"
}

# Rosalind: BA8B — Compute the Squared Error Distortion
# https://rosalind.info/problems/ba8b/
#
# Given: Integers k and m, a set of centers, and a set of points.
# Return: The squared error distortion.

let centers = [[2.31, 4.55], [5.96, 9.08]]
let points = [
    [3.42, 6.03], [6.23, 8.25], [4.76, 1.64], [4.47, 4.33], [3.95, 7.61],
    [8.93, 2.97], [9.74, 4.03], [1.73, 1.28], [9.72, 5.01], [7.27, 3.77],
]

fn distance(a, b) {
    sqrt(range(0, len(a)) |> map(|i| (a[i] - b[i]) * (a[i] - b[i])) |> sum())
}

# Distortion is the mean squared distance from each point to its nearest center
# — squared, so a single badly-placed point counts for much more than several
# slightly-off ones. That is what k-means is minimising.
let distortion = (points
  |> map(|p| { let d = centers |> map(|c| distance(p, c)) |> min()
               d * d })
  |> sum()) / float(len(points))

println("Result:   " + str(round(distortion, 3)))
println("Expected: 18.246")

fn test_ba8b_squared_error_distortion() {
    assert round(distortion, 3) == 18.246, "BA8B: got " + str(round(distortion, 3))
    # A point sitting on a center contributes nothing.
    assert distance([1.0, 1.0], [1.0, 1.0]) == 0.0, "BA8B: distance to itself should be zero"
}

# Rosalind: BA8C — Implement the Lloyd Algorithm for k-Means Clustering
# https://rosalind.info/problems/ba8c/
#
# Given: Integers k and m, then a set of points.
# Return: The centers Lloyd's algorithm converges to, seeded with the first k
# points of the data.

let k = 2
let points = [
    [1.3, 1.1], [1.3, 0.2], [0.6, 2.8], [3.0, 3.2], [1.2, 0.7], [1.4, 1.6],
    [1.2, 1.0], [1.2, 1.1], [0.6, 1.5], [1.8, 2.6], [1.2, 1.3], [1.2, 1.0],
    [0.0, 1.9],
]

fn distance(a, b) {
    sqrt(range(0, len(a)) |> map(|i| (a[i] - b[i]) * (a[i] - b[i])) |> sum())
}

# Two steps, alternating until nothing moves: assign every point to its nearest
# center, then move each center to the mean of what it was assigned. Seeded with
# the first k points, because the problem says so — Lloyd's converges to
# different answers from different seeds, which is why BA8A exists.
let centers = slice(points, 0, k)
let settled = false
let rounds = 0
while settled == false and rounds < 100 {
    let assignment = points |> map(|p| argmin(centers |> map(|c| distance(p, c))))
    let moved = []
    for index in range(0, k) {
        let members = range(0, len(points)) |> filter(|i| assignment[i] == index) |> map(|i| points[i])
        if len(members) == 0 {
            moved = push(moved, centers[index])
        } else {
            moved = push(moved, range(0, len(members[0]))
              |> map(|d| (members |> map(|p| p[d]) |> sum()) / float(len(members))))
        }
    }
    if moved == centers then settled = true
    centers = moved
    rounds = rounds + 1
}

let result = centers |> map(|c| c |> map(|v| str(round(v, 3))) |> join(" ")) |> join(" / ")

println("Result:   " + result)
println("Expected: 1.8 2.867 / 1.06 1.14   (converged in " + str(rounds) + " rounds)")

fn test_ba8c_lloyd_k_means() {
    assert settled, "BA8C: Lloyd's did not converge"
    let flat = centers |> map(|c| c |> map(|v| round(v, 3)))
    assert contains(flat, [1.8, 2.867]), "BA8C: missing the upper center, got " + str(flat)
    assert contains(flat, [1.06, 1.14]), "BA8C: missing the lower center, got " + str(flat)
}

# Rosalind: BA8D — Implement the Soft k-Means Clustering Algorithm
# https://rosalind.info/problems/ba8d/
#
# Given: Integers k and m, a stiffness parameter beta, and n points in m
# dimensions.
# Return: The k centers after 100 rounds of soft k-means.

let k = 2
let beta = 2.7
let steps = 100
let points = [
    [1.3, 1.1], [1.3, 0.2], [0.6, 2.8], [3.0, 3.2], [1.2, 0.7],
    [1.4, 1.6], [1.2, 1.0], [1.2, 1.1], [0.6, 1.5], [1.8, 2.6],
    [1.2, 1.3], [1.2, 1.0], [0.0, 1.9],
]

# BA8C's Lloyd algorithm makes every point choose one cluster. Here each point is
# shared out among all of them, weighted by e^(-beta * distance) — so a point
# midway between two centers pulls on both instead of being forced to pick, and a
# center is the weighted mean of everything rather than of its own members.
#
# beta is how decisive that sharing is. Large beta approaches hard assignment and
# reproduces Lloyd; small beta pulls every center towards the overall mean.
fn distance(a, b) {
    sqrt(range(0, len(a)) |> map(|i| (a[i] - b[i]) * (a[i] - b[i])) |> sum())
}

let centers = range(0, k) |> map(|i| points[i])

for _ in range(0, steps) {
    # E-step: how much each center is responsible for each point.
    let responsibility = centers |> map(|center|
        points |> map(|point| exp(0 - beta * distance(point, center))))

    let totals = range(0, len(points))
        |> map(|j| range(0, k) |> map(|i| responsibility[i][j]) |> sum())

    # M-step: each center becomes the weighted mean of every point.
    centers = range(0, k) |> map(|i| {
        let weights = range(0, len(points)) |> map(|j| responsibility[i][j] / totals[j])
        let weight_sum = sum(weights)
        range(0, len(points[0])) |> map(|d|
            (range(0, len(points)) |> map(|j| weights[j] * points[j][d]) |> sum()) / weight_sum)
    })
}

println("Result:")
for center in centers { println("  " + (center |> map(|c| str(round(c, 3))) |> join(" "))) }
println("Expected:")
println("  1.662 2.623")
println("  1.075 1.148")

fn test_ba8d_soft_k_means() {
    let expected = [[1.662, 2.623], [1.075, 1.148]]
    for i in range(0, k) {
        for d in range(0, 2) {
            assert abs(centers[i][d] - expected[i][d]) < 5e-4,
                "BA8D: center " + str(i) + " coordinate " + str(d)
                    + " is " + str(round(centers[i][d], 3))
                    + ", expected " + str(expected[i][d])
        }
    }
    # Soft assignment means every point contributes to every center, so no center
    # can sit outside the bounding box of the data.
    for center in centers {
        for d in range(0, 2) {
            let column = points |> map(|p| p[d])
            assert center[d] >= min(column) and center[d] <= max(column),
                "BA8D: a weighted mean cannot fall outside the data"
        }
    }
}

# Rosalind: BA8E — Implement Hierarchical Clustering
# https://rosalind.info/problems/ba8e/
#
# Given: An integer n and an n x n distance matrix.
# Return: The clusters formed, in the order they are merged.

let n = 7
let distances = [
    [0.00, 0.74, 0.85, 0.54, 0.83, 0.92, 0.89],
    [0.74, 0.00, 1.59, 1.35, 1.20, 1.48, 1.55],
    [0.85, 1.59, 0.00, 0.63, 1.13, 0.69, 0.73],
    [0.54, 1.35, 0.63, 0.00, 0.66, 0.43, 0.88],
    [0.83, 1.20, 1.13, 0.66, 0.00, 0.72, 0.55],
    [0.92, 1.48, 0.69, 0.43, 0.72, 0.00, 0.80],
    [0.89, 1.55, 0.73, 0.88, 0.55, 0.80, 0.00],
]

# Unlike k-means, nothing has to be told how many clusters there are. Every point
# starts alone, the two closest merge, and the process is repeated — producing not
# one partition but a whole nested family of them, which is what a phylogenetic
# tree is.
#
# "Closest" here is average linkage: the mean distance between all pairs across
# the two clusters. That choice matters — single linkage would chain elongated
# clusters together, complete linkage would insist on compactness — and it is the
# one this problem specifies.
let clusters = range(0, n) |> map(|i| [i])
let merged_order = []

fn average_distance(left, right, matrix) {
    let total = left |> flat_map(|a| right |> map(|b| matrix[a][b])) |> sum()
    total / (len(left) * len(right))
}

while len(clusters) > 1 {
    let best_i = 0
    let best_j = 1
    let best = average_distance(clusters[0], clusters[1], distances)
    for i in range(0, len(clusters)) {
        for j in range(i + 1, len(clusters)) {
            let candidate = average_distance(clusters[i], clusters[j], distances)
            if candidate < best {
                best = candidate
                best_i = i
                best_j = j
            }
        }
    }

    let joined = clusters[best_i] + clusters[best_j]
    merged_order = push(merged_order, joined)
    clusters = range(0, len(clusters))
        |> filter(|i| i != best_i and i != best_j)
        |> map(|i| clusters[i])
    clusters = push(clusters, joined)
}

# Reported one-based, which is how the sample is written.
let reported = merged_order |> map(|cluster| cluster |> map(|i| str(i + 1)) |> join(" "))

println("Result:")
for line in reported { println("  " + line) }
println("Expected:")
println("  4 6 / 5 7 / 3 4 6 / 1 2 / 5 7 3 4 6 / 1 2 5 7 3 4 6")

fn test_ba8e_hierarchical_clustering() {
    assert len(merged_order) == n - 1, "BA8E: n points take n-1 merges"
    # The first merge must be the closest pair in the matrix, which is 4 and 6 at
    # 0.43 — checked against the matrix rather than against the expected output.
    assert reported[0] == "4 6", "BA8E: first merge was " + reported[0]
    assert distances[3][5] == 0.43, "BA8E: 4 and 6 are 0.43 apart"
    let closest = range(0, n) |> flat_map(|i| range(i + 1, n) |> map(|j| distances[i][j])) |> min()
    assert closest == 0.43, "BA8E: and nothing is closer"
    # The last merge gathers everything.
    assert len(merged_order[n - 2]) == n, "BA8E: the final cluster holds every point"
    assert sort(merged_order[n - 2]) == range(0, n), "BA8E: and holds each point once"
    assert join(reported, " / ") == "4 6 / 5 7 / 3 4 6 / 1 2 / 5 7 3 4 6 / 1 2 5 7 3 4 6",
        "BA8E: got " + join(reported, " / ")
}