Phylogeny

19 problems from Rosalind — Bioinformatics Stronghold. Press Run on any block to execute it in your browser.

# Rosalind: TREE — Completing a Tree
# https://rosalind.info/problems/tree/
#
# Given: A positive integer n and an adjacency list of a graph on n nodes that
# forms a forest.
# Return: The minimum number of edges needed to produce a tree.

let n = 10
let edge_list = [[1,2],[2,8],[4,10],[5,9],[6,10],[7,9]]

# A tree on n nodes has exactly n-1 edges, and a forest with e edges has
# n-e components — so joining them needs (n-1)-e more.
let result = (n - 1) - len(edge_list)

println("Result:   " + str(result))
println("Expected: 3")

fn test_tree_edges_to_add() {
    assert result == 3, "TREE: got " + str(result)
}

# Rosalind: INOD — Counting Phylogenetic Ancestors
# https://rosalind.info/problems/inod/
#
# Given: A positive integer n (3 <= n <= 10000).
# Return: The number of internal nodes of any unrooted binary tree with n leaves.

let n = 4

# Every internal node of an unrooted binary tree has degree 3. Counting edge
# endpoints two ways gives internal = n - 2, independent of the shape.
let result = n - 2

println("Result:   " + str(result))
println("Expected: 2")

fn test_inod_internal_nodes() {
    assert result == 2, "INOD: got " + str(result)
}

# Rosalind: PDST — Creating a Distance Matrix
# https://rosalind.info/problems/pdst/
#
# Given: A collection of DNA strings of equal length in FASTA format.
# Return: The matrix D of p-distances, where D[i][j] is the proportion of
# positions at which strings i and j differ.

let strings = [
    "TTTCCATTTA",
    "GATTCATTTC",
    "TTTCCATTTT",
    "GTTCCATTTA"
]

fn p_distance(a, b) {
    let differing = range(0, len(a)) |> count_if(|i| substr(a, i, 1) != substr(b, i, 1))
    float(differing) / float(len(a))
}

let distance_grid = strings |> map(|a| strings |> map(|b| p_distance(a, b)))

println("Result:")
distance_grid |> each(|row| println("  " + (row |> map(|d| str(round(d, 5))) |> join(" "))))
println("Expected first row: 0.00000 0.40000 0.10000 0.10000")

fn test_pdst_distance_matrix() {
    assert round(distance_grid[0][0], 5) == 0.0, "PDST: diagonal not zero"
    assert round(distance_grid[0][1], 5) == 0.4, "PDST: [0][1] was " + str(distance_grid[0][1])
    assert round(distance_grid[0][2], 5) == 0.1, "PDST: [0][2] was " + str(distance_grid[0][2])
    assert round(distance_grid[0][3], 5) == 0.1, "PDST: [0][3] was " + str(distance_grid[0][3])
    assert round(distance_grid[1][0], 5) == round(distance_grid[0][1], 5), "PDST: not symmetric"
}

# Rosalind: ROOT — Counting Rooted Binary Trees
# https://rosalind.info/problems/root/
#
# Given: A positive integer n (n <= 1000).
# Return: The number of rooted binary trees on n labeled leaves, modulo
# 1,000,000.

let n = 5
let modulus = 1000000

# Adding the k-th leaf splits any of the (2k-3) existing edges, plus the root
# edge — so the count is the double factorial (2n-3)!!, built up one leaf at a
# time and reduced as we go.
let result = range(2, n + 1) |> reduce(|acc, k| (acc * (2 * k - 3)) % modulus, 1)

println("Result:   " + str(result))
println("Expected: 105   (7!! = 7 x 5 x 3 x 1)")

fn test_root_rooted_binary_trees() {
    assert result == 105, "ROOT: got " + str(result)
}

# Rosalind: CUNR — Counting Unrooted Binary Trees
# https://rosalind.info/problems/cunr/
#
# Given: A positive integer n (n <= 1000).
# Return: The number of unrooted binary trees on n labeled leaves, modulo
# 1,000,000.

let n = 5
let modulus = 1000000

# An unrooted tree on n leaves is a rooted tree on n-1 leaves with the root
# edge removed, so the count drops one double-factorial step to (2n-5)!!.
let result = range(3, n + 1) |> reduce(|acc, k| (acc * (2 * k - 5)) % modulus, 1)

println("Result:   " + str(result))
println("Expected: 15   (5!! = 5 x 3 x 1)")

fn test_cunr_unrooted_binary_trees() {
    assert result == 15, "CUNR: got " + str(result)
}

# Rosalind: NWCK — Distances in Trees
# https://rosalind.info/problems/nwck/
#
# Given: A collection of Newick trees, each followed by a pair of node names.
# Return: For each pair, the number of edges on the path between them.

let queries = [
    { tree: "(cat)dog;", from_node: "dog", to_node: "cat" },
    { tree: "((cat)dog,robot);", from_node: "dog", to_node: "robot" }
]

# ── A Newick parser ──────────────────────────────────────────
#
# Nodes are held in flat lists indexed by id: `names[i]` and `parents[i]`.
# A '(' opens an internal node; the label after the matching ')' names it. A
# bare label is a leaf. That is all the structure these problems need — branch
# lengths are ignored here, and read separately by the weighted variants.
fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)

        if c == "(" {
            # Flush any label sitting before the bracket, then open a child.
            let parent = stack[len(stack) - 1]
            names = push(names, "")
            parents = push(parents, parent)
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                # A pending label either names the node just closed, or is a
                # new leaf under the current parent.
                if pending != "" {
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k])
                    } else {
                        names = push(names, pending)
                        parents = push(parents, stack[len(stack) - 1])
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }

    { names: names, parents: parents }
}

# Edges are undirected for path-length purposes.
fn neighbours_of(tree, node) {
    let up = if tree.parents[node] >= 0 then [tree.parents[node]] else []
    let down = range(0, len(tree.parents)) |> filter(|k| tree.parents[k] == node)
    concat(up, down)
}

fn index_of_name(tree, name) {
    (range(0, len(tree.names)) |> filter(|k| tree.names[k] == name))[0]
}

# Breadth-first search: every edge counts as one step.
fn distance_between(tree, a, b) {
    let start = index_of_name(tree, a)
    let goal = index_of_name(tree, b)
    let frontier = [start]
    let seen = [start]
    let steps = 0
    let answer = -1
    while len(frontier) > 0 and answer < 0 {
        if frontier |> contains(goal) then answer = steps
        if answer < 0 {
            let next = frontier
                |> flat_map(|node| neighbours_of(tree, node))
                |> filter(|node| !(seen |> contains(node)))
                |> unique()
            seen = concat(seen, next)
            frontier = next
            steps = steps + 1
        }
    }
    answer
}

let results = queries |> map(|q| distance_between(parse_newick(q.tree), q.from_node, q.to_node))

println("Result:   " + (results |> map(|d| str(d)) |> join(" ")))
println("Expected: 1 2")

fn test_nwck_tree_distances() {
    assert results[0] == 1, "NWCK[0]: got " + str(results[0])
    assert results[1] == 2, "NWCK[1]: got " + str(results[1])
}

# Rosalind: NKEW — Newick Format with Edge Weights
# https://rosalind.info/problems/nkew/
#
# Given: Newick trees carrying branch lengths, each followed by a pair of node
# names.
# Return: For each pair, the total weight of the path between them.

let queries = [
    { tree: "(dog:42,cat:33);", from_node: "cat", to_node: "dog" },
    { tree: "((dog:4,cat:3):74,robot:98,elephant:58);", from_node: "dog", to_node: "cat" }
]

# The same parser as NWCK, extended to read the ":length" suffix. A label is
# "name:weight"; either half may be empty, so an unnamed internal node can
# still carry a length.
fn split_label(label) {
    let colon = index_of(label, ":")
    if colon < 0 {
        { name: label, weight: 0.0 }
    } else {
        { name: substr(label, 0, colon), weight: float(substr(label, colon + 1, len(label) - colon - 1)) }
    }
}

fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let weights = [0.0]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)
        if c == "(" {
            names = push(names, "")
            parents = push(parents, stack[len(stack) - 1])
            weights = push(weights, 0.0)
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                if pending != "" {
                    let parts = split_label(pending)
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then parts.name else names[k])
                        weights = range(0, len(weights)) |> map(|k| if k == just_closed then parts.weight else weights[k])
                    } else {
                        names = push(names, parts.name)
                        parents = push(parents, stack[len(stack) - 1])
                        weights = push(weights, parts.weight)
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }

    { names: names, parents: parents, weights: weights }
}

fn index_of_name(tree, name) {
    (range(0, len(tree.names)) |> filter(|k| tree.names[k] == name))[0]
}

# Path to the root, as a list of node ids.
fn ancestry(tree, node) {
    let path = [node]
    let at = node
    while tree.parents[at] >= 0 {
        at = tree.parents[at]
        path = push(path, at)
    }
    path
}

# The path between two nodes runs up to their lowest common ancestor and back
# down; each node contributes the weight of the edge to its own parent.
fn weighted_distance(tree, a, b) {
    let up_a = ancestry(tree, index_of_name(tree, a))
    let up_b = ancestry(tree, index_of_name(tree, b))
    let shared = up_a |> filter(|node| up_b |> contains(node))
    let meeting = shared[0]
    let side_a = up_a |> filter(|node| node != meeting and !(up_b |> contains(node)))
    let side_b = up_b |> filter(|node| node != meeting and !(up_a |> contains(node)))
    let total_a = side_a |> map(|node| tree.weights[node]) |> sum()
    let total_b = side_b |> map(|node| tree.weights[node]) |> sum()
    total_a + total_b
}

let results = queries |> map(|q| weighted_distance(parse_newick(q.tree), q.from_node, q.to_node))

println("Result:   " + (results |> map(|d| str(int(d))) |> join(" ")))
println("Expected: 75 7")

fn test_nkew_weighted_distances() {
    assert int(results[0]) == 75, "NKEW[0]: got " + str(results[0])
    assert int(results[1]) == 7, "NKEW[1]: got " + str(results[1])
}

# Rosalind: CTBL — Creating a Character Table
# https://rosalind.info/problems/ctbl/
#
# Given: An unrooted binary tree in Newick format with n leaves.
# Return: The nontrivial characters the tree induces, each as a bit string over
# the leaves in alphabetical order.

let newick = "(dog,((elephant,mouse),robot),cat);"

# Newick parser: nodes in flat lists, '(' opens an internal node, the label
# after the matching ')' names it, a bare label is a leaf.
fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)
        if c == "(" {
            names = push(names, "")
            parents = push(parents, stack[len(stack) - 1])
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                if pending != "" {
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k])
                    } else {
                        names = push(names, pending)
                        parents = push(parents, stack[len(stack) - 1])
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }
    { names: names, parents: parents }
}

let tree = parse_newick(newick)

# A leaf is any node that is nobody's parent.
let leaf_ids = range(0, len(tree.names))
    |> filter(|k| (range(0, len(tree.parents)) |> count_if(|c| tree.parents[c] == k)) == 0)
let taxa = leaf_ids |> map(|k| tree.names[k]) |> sort()

fn descends_from(tree, node, ancestor) {
    let at = node
    let found = false
    while at >= 0 and !found {
        if at == ancestor then found = true
        at = tree.parents[at]
    }
    found
}

# Each internal node defines an edge to its parent, and cutting that edge
# splits the leaves in two. Splits of size 1 or n-1 are trivial — every tree
# has them — so only the rest are characters.
let internal_ids = range(0, len(tree.names)) |> filter(|k| !(leaf_ids |> contains(k)))

let characters = internal_ids |> flat_map(|node| {
    let inside = taxa |> map(|name| {
        let leaf = (leaf_ids |> filter(|k| tree.names[k] == name))[0]
        if descends_from(tree, leaf, node) then "1" else "0"
    })
    let size = inside |> count_if(|b| b == "1")
    if size > 1 and size < len(taxa) - 1 then [inside |> join("")] else []
}) |> unique() |> sort()

println("Taxa:     " + (taxa |> join(" ")))
println("Result:")
characters |> each(|c| println("  " + c))
println("Expected: 00110 and 00111")

fn test_ctbl_character_table() {
    assert len(characters) == 2, "CTBL: got " + str(len(characters)) + " characters"
    assert characters |> contains("00110"), "CTBL: missing 00110"
    assert characters |> contains("00111"), "CTBL: missing 00111"
}

# Rosalind: SPTD — Phylogeny Comparison with Split Distance
# https://rosalind.info/problems/sptd/
#
# Given: A list of n taxa and two unrooted binary trees over them.
# Return: The split distance between the trees.

let taxa = ["dog", "rat", "elephant", "mouse", "cat", "rabbit"] |> sort()
let first_tree = "(rat,(dog,cat),(rabbit,(elephant,mouse)));"
let second_tree = "(rat,(cat,(dog,mouse)),(elephant,rabbit));"

fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)
        if c == "(" {
            names = push(names, "")
            parents = push(parents, stack[len(stack) - 1])
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                if pending != "" {
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k])
                    } else {
                        names = push(names, pending)
                        parents = push(parents, stack[len(stack) - 1])
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }
    { names: names, parents: parents }
}

fn descends_from(tree, node, ancestor) {
    let at = node
    let found = false
    while at >= 0 and !found {
        if at == ancestor then found = true
        at = tree.parents[at]
    }
    found
}

# The nontrivial splits of a tree, each as a bit string over `taxa`.
#
# A split and its complement describe the same edge, so each is stored in a
# canonical orientation — the one where the first taxon is a 0 — otherwise the
# same split could be counted as two different ones.
fn splits_of(newick, all_taxa) {
    let tree = parse_newick(newick)
    let leaf_ids = range(0, len(tree.names))
        |> filter(|k| (range(0, len(tree.parents)) |> count_if(|c| tree.parents[c] == k)) == 0)
    let internal_ids = range(0, len(tree.names)) |> filter(|k| !(leaf_ids |> contains(k)))

    internal_ids |> flat_map(|node| {
        let bits = all_taxa |> map(|name| {
            let leaf = (leaf_ids |> filter(|k| tree.names[k] == name))[0]
            if descends_from(tree, leaf, node) then "1" else "0"
        })
        let size = bits |> count_if(|b| b == "1")
        if size > 1 and size < len(all_taxa) - 1 {
            let canonical = if bits[0] == "1" {
                bits |> map(|b| if b == "1" then "0" else "1")
            } else {
                bits
            }
            [canonical |> join("")]
        } else {
            []
        }
    }) |> unique()
}

let first_splits = splits_of(first_tree, taxa)
let second_splits = splits_of(second_tree, taxa)
let shared = first_splits |> count_if(|s| second_splits |> contains(s))

# An unrooted binary tree on n taxa has n-3 nontrivial splits. The distance
# counts the splits unique to each tree, so it is 2(n-3) minus twice the
# shared ones.
let n = len(taxa)
let result = 2 * (n - 3) - 2 * shared

println("Splits in tree 1: " + (first_splits |> join(" ")))
println("Splits in tree 2: " + (second_splits |> join(" ")))
println("Shared:   " + str(shared))
println("Result:   " + str(result))
println("Expected: 6 — these two trees share no nontrivial split:")
println("  tree 1: {cat,dog} {elephant,mouse} {elephant,mouse,rabbit}")
println("  tree 2: {cat,dog,mouse} {dog,mouse} {elephant,rabbit}")

# Comparing a tree with itself must give 0, which checks the canonical
# orientation as much as the arithmetic: without it the same split could be
# recorded two ways and fail to match itself.
let self_shared = first_splits |> count_if(|s| first_splits |> contains(s))
let self_distance = 2 * (n - 3) - 2 * self_shared

fn test_sptd_split_distance() {
    assert len(first_splits) == n - 3, "SPTD: tree 1 has " + str(len(first_splits)) + " splits"
    assert len(second_splits) == n - 3, "SPTD: tree 2 has " + str(len(second_splits)) + " splits"
    assert self_distance == 0, "SPTD: a tree compared with itself gave " + str(self_distance)
    assert result == 6, "SPTD: got " + str(result)
    # The distance is even and cannot exceed twice the split count.
    assert result % 2 == 0, "SPTD: distance is odd"
    assert result >= 0 and result <= 2 * (n - 3), "SPTD: distance out of range"
}

# Rosalind: MEND — Inferring Genotype from a Pedigree
# https://rosalind.info/problems/mend/
#
# Given: A rooted binary tree in Newick format whose leaves are the genotypes
# of an individual's ancestors, for a factor with alleles A and a.
# Return: The probabilities that the root individual is AA, Aa and aa.

let newick = "((((Aa,aa),(Aa,Aa)),((aa,aa),(AA,Aa))),(((aa,AA),(aa,Aa)),((aa,aa),(aa,AA))));"

# Newick parser: nodes in flat lists, '(' opens an internal node, a bare label
# is a leaf.
fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)
        if c == "(" {
            names = push(names, "")
            parents = push(parents, stack[len(stack) - 1])
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                if pending != "" {
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k])
                    } else {
                        names = push(names, pending)
                        parents = push(parents, stack[len(stack) - 1])
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }
    { names: names, parents: parents }
}

let tree = parse_newick(newick)

fn children_of(node) {
    range(0, len(tree.parents)) |> filter(|k| tree.parents[k] == node)
}

# A genotype distribution is [P(AA), P(Aa), P(aa)].
fn genotype_of(label) {
    if label == "AA" then [1.0, 0.0, 0.0] else
    if label == "Aa" then [0.0, 1.0, 0.0] else [0.0, 0.0, 1.0]
}

# A parent passes allele A with probability P(AA) + P(Aa)/2, and the two
# parents contribute independently — so the child's distribution follows from
# just those two numbers.
fn cross(left, right) {
    let a = left[0] + left[1] / 2.0
    let b = right[0] + right[1] / 2.0
    [a * b, a * (1.0 - b) + (1.0 - a) * b, (1.0 - a) * (1.0 - b)]
}

fn distribution_at(node) {
    let kids = children_of(node)
    if len(kids) == 0 {
        genotype_of(tree.names[node])
    } else {
        cross(distribution_at(kids[0]), distribution_at(kids[1]))
    }
}

# Node 0 is the implicit root the parser starts from; the tree's own root is
# its single child.
let root = children_of(0)[0]
let result = distribution_at(root)

println("Result:   " + (result |> map(|p| str(round(p, 3))) |> join(" ")))
println("Expected: 0.117 0.453 0.430")
println("")
println("Worked up from the leaves: the left half of the pedigree reaches")
println("[0.1406, 0.4688, 0.3906] and the right half [0.0938, 0.4375, 0.4688],")
println("so the root parents pass A with probability 0.375 and 0.3125 —")
println("giving 0.375 x 0.3125 = 0.1172 for AA.")

fn test_mend_pedigree_genotypes() {
    assert round(result[0], 3) == 0.117, "MEND: P(AA) = " + str(result[0])
    assert round(result[1], 3) == 0.453, "MEND: P(Aa) = " + str(result[1])
    assert round(result[2], 3) == 0.43, "MEND: P(aa) = " + str(result[2])
    # A probability distribution must sum to one.
    let total = result |> sum()
    assert abs(total - 1.0) < 0.000001, "MEND: distribution sums to " + str(total)
}

# Rosalind: CSET — Fixing an Inconsistent Character Set
# https://rosalind.info/problems/cset/
#
# Given: A collection of characters over the same taxa, at most one of which is
# inconsistent with the others.
# Return: The remaining characters, which form a consistent set.

let characters = ["100001", "000110", "111000", "100111"]

# Two splits are compatible when one of the four ways of intersecting their
# halves is empty — the four-gamete condition. If all four appear, no tree can
# carry both characters at once.
fn compatible(a, b) {
    let seen = range(0, len(a)) |> map(|i| substr(a, i, 1) ++ substr(b, i, 1)) |> unique()
    len(seen) < 4
}

fn all_compatible(rows) {
    let bad = range(0, len(rows)) |> flat_map(|i| {
        range(i + 1, len(rows)) |> filter(|j| !compatible(rows[i], rows[j]))
    })
    len(bad) == 0
}

# Drop each character in turn and keep the first set that becomes consistent.
let kept = (range(0, len(characters))
    |> map(|i| concat(take(characters, i), drop(characters, i + 1)))
    |> filter(|rows| all_compatible(rows)))[0]

let removed = characters |> filter(|c| !(kept |> contains(c)))

println("Result:")
kept |> each(|c| println("  " + c))
println("Removed:  " + (removed |> join(" ")))
println("Expected: one character removed, leaving a consistent set")
println("          (100001 and 111000 are the incompatible pair here)")

fn test_cset_consistent_character_set() {
    assert len(kept) == len(characters) - 1, "CSET: kept " + str(len(kept)) + " characters"
    # The point of the answer is consistency, so check it rather than a string.
    assert all_compatible(kept), "CSET: the kept set is still inconsistent"
    # The original set really was inconsistent, otherwise the problem is empty.
    assert !all_compatible(characters), "CSET: the input was already consistent"
}

# Rosalind: CSTR — Creating a Character Table from Genetic Strings
# https://rosalind.info/problems/cstr/
#
# Given: A collection of characterizable DNA strings of equal length.
# Return: A character table whose rows are the nontrivial characters, one per
# position where the strings disagree.

let strings = [
    "ATGCTACC",
    "CGTTTACC",
    "ATTCGACC",
    "AGTCTCCC",
    "CGTCTATC",
]

let width = len(strings[0])

# A position splits the collection into two groups. The split is nontrivial only
# when both groups have at least two members: a position where one string
# differs from all the others separates nothing, and neither does a position
# where every string agrees.
let characters = []
for column in range(0, width) {
    let symbols = strings |> map(|s| substr(s, column, 1))
    let distinct = unique(symbols)
    if len(distinct) == 2 {
        let first_count = symbols |> filter(|c| c == distinct[0]) |> len()
        let second_count = len(symbols) - first_count
        if first_count >= 2 and second_count >= 2 {
            # Which state gets 1 is arbitrary, so take the first symbol seen.
            characters = push(characters, symbols |> map(|c| if c == distinct[0] then "1" else "0") |> join(""))
        }
    }
}

println("Result:")
for character in characters {
    println("  " + character)
}
println("Expected:")
println("  10110")
println("  10100")

fn test_cstr_character_table() {
    assert len(characters) == 2, "CSTR: expected 2 nontrivial characters, got " + str(len(characters))
    # The 0/1 assignment is arbitrary, so a row and its complement are the same
    # character. Compare each row against the expected one either way round.
    fn flipped(row) {
        range(0, len(row)) |> map(|i| if substr(row, i, 1) == "1" then "0" else "1") |> join("")
    }
    assert characters[0] == "10110" or flipped(characters[0]) == "10110", "CSTR: first row " + characters[0]
    assert characters[1] == "10100" or flipped(characters[1]) == "10100", "CSTR: second row " + characters[1]
}

# Rosalind: ALPH — Alignment-Based Phylogeny
# https://rosalind.info/problems/alph/
#
# Given: A rooted binary tree in Newick format and a multiple alignment of its
# leaves.
# Return: The minimum total Hamming distance over the tree's edges, and internal
# node labels achieving it.

# (((ostrich,cat)rat,(duck,fly)mouse)dog,(elephant,pikachu)hamster)robot;
let children = {
    "rat": ["ostrich", "cat"],
    "mouse": ["duck", "fly"],
    "dog": ["rat", "mouse"],
    "hamster": ["elephant", "pikachu"],
    "robot": ["dog", "hamster"],
}
let root = "robot"
let leaves = {
    "ostrich": "AC", "cat": "CA", "duck": "T-",
    "fly": "GC", "elephant": "-T", "pikachu": "AA",
}

# The gap is a fifth symbol here, not a missing value. That is the whole reason
# this problem is separate from the ordinary parsimony ones: in an alignment a
# gap is evidence — an insertion or deletion actually happened — and treating it
# as unknown would make every column containing one free.
let symbols = ["A", "C", "G", "T", "-"]
let width = len(leaves["ostrich"])

fn is_leaf(node, leaf_map) { contains(keys(leaf_map), node) }

# Sankoff: the cheapest cost of a subtree given its root's symbol is the sum over
# children of their cheapest cost plus one if the symbol has to change. Columns
# are independent, so each is solved separately and the scores added.
fn score_column(node, position, leaf_map, child_map, alphabet) {
    if is_leaf(node, leaf_map) {
        let here = substr(leaf_map[node], position, 1)
        return {
            costs: alphabet |> map(|s| if s == here then 0 else 1000000),
            picks: {}, solved: [], kids: [],
        }
    }
    let kids = child_map[node]
    let solved = kids |> map(|kid| score_column(kid, position, leaf_map, child_map, alphabet))
    let picks = {}
    let costs = range(0, len(alphabet)) |> map(|mine| {
        range(0, len(kids)) |> map(|k| {
            let options = range(0, len(alphabet))
                |> map(|theirs| solved[k].costs[theirs] + (if theirs == mine then 0 else 1))
            let chosen = argmin(options)
            picks[str(mine) + "," + str(k)] = chosen
            options[chosen]
        }) |> sum()
    })
    { costs: costs, picks: picks, solved: solved, kids: kids }
}

fn assign(node, tree, chosen_index, labels, alphabet) {
    labels[node] = labels[node] + alphabet[chosen_index]
    for k in range(0, len(tree.kids)) {
        if len(tree.solved[k].kids) > 0 {
            labels = assign(tree.kids[k], tree.solved[k],
                            tree.picks[str(chosen_index) + "," + str(k)], labels, alphabet)
        }
    }
    labels
}

let labels = {}
for node in keys(children) { labels[node] = "" }
let total = 0
for position in range(0, width) {
    let solved = score_column(root, position, leaves, children, symbols)
    let best = argmin(solved.costs)
    total = total + solved.costs[best]
    labels = assign(root, solved, best, labels, symbols)
}

let named = {}
for node in keys(leaves) { named[node] = leaves[node] }
for node in keys(children) { named[node] = labels[node] }

fn hamming(a, b) { range(0, len(a)) |> count_if(|i| substr(a, i, 1) != substr(b, i, 1)) }

println("Result:   " + str(total))
for node in sort(keys(children)) { println("  " + node + ": " + named[node]) }
println("Expected: 8  (rat AC, mouse TC, dog AC, hamster AT, robot AC)")
println("This labelling differs and costs the same — several are optimal, and the")
println("assertion recounts the changes rather than trusting the printed one.")

fn test_alph_alignment_based_phylogeny() {
    assert total == 8, "ALPH: scored " + str(total)
    # The score has to be the changes the labelling actually shows — a number
    # without a matching labelling is the usual bug here.
    let shown = sort(keys(children)) |> flat_map(|parent|
        children[parent] |> map(|kid| hamming(named[parent], named[kid]))) |> sum()
    assert shown == total,
        "ALPH: the labelling shows " + str(shown) + " changes but the score is " + str(total)
    # Every internal label is the right length and drawn from the alphabet.
    for node in keys(children) {
        assert len(named[node]) == width, "ALPH: " + node + " has the wrong length"
        assert (chars(named[node]) |> count_if(|c| contains(symbols, c) == false)) == 0,
            "ALPH: " + named[node] + " uses a symbol outside the alphabet"
    }
    # Treating the gap as free would score lower, which is exactly what must not
    # happen — it is a real evolutionary event.
    assert (chars(leaves["duck"]) |> count_if(|c| c == "-")) == 1, "ALPH: duck carries a gap"
}

# Rosalind: CNTQ — Counting Quartets
# https://rosalind.info/problems/cntq/
#
# Given: n and an unrooted binary tree on n taxa in Newick format.
# Return: The number of quartets consistent with the tree, modulo 1,000,000.

let taxa = ["lobster", "cat", "dog", "caterpillar", "elephant", "mouse"]
# (lobster,(cat,dog),(caterpillar,(elephant,mouse)));
# Each internal edge splits the taxa in two; these are the non-trivial sides.
let clades = [
    ["cat", "dog"],
    ["caterpillar", "elephant", "mouse"],
    ["elephant", "mouse"],
]

# The answer is simply C(n,4), and the reason is worth more than the number: in a
# *fully resolved* binary tree, any four taxa are separated by some internal edge
# into two pairs, so every four-taxon set contributes exactly one consistent
# quartet. Nothing about the tree's shape enters into it.
#
# That stops being true the moment the tree has an unresolved node, which is why
# QRT and CHBP — where splits are partial — are real problems and this one is
# not.
let n = len(taxa)
let answer = choose(n, 4) % 1000000

println("Result:   " + str(answer))
println("Expected: 15")

fn test_cntq_counting_quartets() {
    assert answer == 15, "CNTQ: got " + str(answer)

    # Demonstrated rather than asserted: enumerate all four-taxon subsets and
    # confirm each really is separated into two pairs by some edge.
    fn subsets_of_four(items) {
        range(0, len(items)) |> flat_map(|a|
            range(a + 1, len(items)) |> flat_map(|b|
                range(b + 1, len(items)) |> flat_map(|c|
                    range(c + 1, len(items)) |> map(|d|
                        [items[a], items[b], items[c], items[d]]))))
    }

    let quartets = subsets_of_four(taxa)
    assert len(quartets) == choose(n, 4), "CNTQ: C(6,4) is 15 subsets"

    let consistent = quartets |> count_if(|four|
        (clades |> count_if(|clade|
            (four |> count_if(|t| contains(clade, t))) == 2)) > 0)
    assert consistent == len(quartets),
        "CNTQ: " + str(consistent) + " of " + str(len(quartets)) + " subsets are separated"
    assert consistent == answer, "CNTQ: the count and the formula must agree"
}

# Rosalind: QRT — Quartets
# https://rosalind.info/problems/qrt/
#
# Given: A partial character table.
# Return: Every quartet inferable from the splits its characters describe.

let taxa = ["cat", "dog", "elephant", "ostrich", "mouse", "rabbit", "robot"]
let characters = [
    "01xxx00",
    "x11xx00",
    "111x00x",
]

# A partial character is one that could not be scored for every taxon — an 'x'
# means "not known", not "third state". So a character constrains only the taxa
# it actually saw, and the quartets it supports are pairs drawn from its two
# scored groups.
#
# The whole point is that missing data does not make a character useless. It
# still separates the taxa it did score, and those separations are what a tree
# gets built from.
fn group_of(character, symbol, names) {
    range(0, len(names)) |> filter(|i| substr(character, i, 1) == symbol) |> map(|i| names[i])
}

fn pairs_of(items) {
    range(0, len(items)) |> flat_map(|a|
        range(a + 1, len(items)) |> map(|b| [items[a], items[b]]))
}

# A quartet has no orientation — {a,b}|{c,d} is the same as {c,d}|{a,b} — so each
# is recorded in a canonical form and duplicates dropped. Two characters often
# support the same quartet, and counting it twice would overstate the evidence.
fn canonical(left, right) {
    let one = join(sort(left), ",")
    let two = join(sort(right), ",")
    if one < two then one + "|" + two else two + "|" + one
}

let seen = {}
let quartets = []
for character in characters {
    let zeros = group_of(character, "0", taxa)
    let ones = group_of(character, "1", taxa)
    for left in pairs_of(zeros) {
        for right in pairs_of(ones) {
            let key = canonical(left, right)
            if contains(keys(seen), key) == false {
                seen[key] = true
                quartets = push(quartets, { left: left, right: right })
            }
        }
    }
}

let written = quartets
    |> map(|q| "{" + join(q.left, ", ") + "} {" + join(q.right, ", ") + "}")
    |> sort()

println("Result:")
for line in written { println("  " + line) }
println("Expected (any order): {elephant, dog} {rabbit, robot} / {cat, dog} {mouse, rabbit}")
println("                      {mouse, rabbit} {cat, elephant} / {dog, elephant} {mouse, rabbit}")

fn test_qrt_quartets() {
    assert len(quartets) == 4, "QRT: expected 4 quartets, got " + str(len(quartets))
    # Compared as unordered pairs of unordered pairs, since neither side nor the
    # order within a side carries meaning.
    let expected = [
        canonical(["elephant", "dog"], ["rabbit", "robot"]),
        canonical(["cat", "dog"], ["mouse", "rabbit"]),
        canonical(["mouse", "rabbit"], ["cat", "elephant"]),
        canonical(["dog", "elephant"], ["mouse", "rabbit"]),
    ]
    let got = quartets |> map(|q| canonical(q.left, q.right))
    assert sort(got) == sort(expected), "QRT: got " + join(sort(got), " ")
    # Every quartet's four taxa are distinct and really were scored by some
    # character — an 'x' can never appear in one.
    for q in quartets {
        let four = q.left + q.right
        assert len(unique(four)) == 4, "QRT: a quartet must name four different taxa"
        let supporting = characters |> count_if(|c|
            (q.left |> count_if(|t| substr(c, (range(0, len(taxa))
                |> filter(|i| taxa[i] == t))[0], 1) == "0")) == 2
            and (q.right |> count_if(|t| substr(c, (range(0, len(taxa))
                |> filter(|i| taxa[i] == t))[0], 1) == "1")) == 2)
        assert supporting > 0, "QRT: no character supports " + canonical(q.left, q.right)
    }
}

# Rosalind: CHBP — Character-Based Phylogeny
# https://rosalind.info/problems/chbp/
#
# Given: Species names and a consistent character table.
# Return: An unrooted binary tree in Newick format modelling the table.

let taxa = ["cat", "dog", "elephant", "mouse", "rabbit", "rat"]
let characters = ["011101", "001101", "001100"]

# The inverse of CTBL: there, a tree produced splits; here the splits have to
# produce the tree. It works because a consistent table's splits are *laminar* —
# any two are nested or disjoint, never crossing — and a laminar family is
# exactly a tree.
#
# An unrooted tree has no first node, so a reference taxon is picked and every
# split is taken from the side away from it. That makes the clades nest, and the
# reference becomes the outermost branch.
let reference = taxa[0]

fn side_away_from(character, names, anchor) {
    let anchor_index = (range(0, len(names)) |> filter(|i| names[i] == anchor))[0]
    let anchor_state = substr(character, anchor_index, 1)
    range(0, len(names))
        |> filter(|i| substr(character, i, 1) != anchor_state)
        |> map(|i| names[i])
}

let clades = characters |> map(|c| side_away_from(c, taxa, reference))

fn is_subset(small, big) { (small |> count_if(|x| contains(big, x) == false)) == 0 }

fn build(members, groups) {
    # The groups strictly inside this one, and of those the maximal — anything
    # contained in another is handled a level deeper.
    let inside = groups |> filter(|g| is_subset(g, members) and len(g) < len(members))
    let maximal = inside |> filter(|g|
        (inside |> count_if(|other| len(other) > len(g) and is_subset(g, other))) == 0)
    let covered = maximal |> flat_map(|g| g)
    let loose = members |> filter(|m| contains(covered, m) == false)
    let parts = (maximal |> map(|g| build(g, groups))) + loose
    if len(parts) == 1 then parts[0] else "(" + join(parts, ",") + ")"
}

let rest = taxa |> filter(|t| t != reference)
let inner = build(rest, clades)
# The reference sits alongside the rest at the unrooted centre, so its branch is
# spliced in rather than wrapped around.
let newick = "(" + reference + "," + substr(inner, 1, len(inner) - 2) + ");"

println("Result:   " + newick)
println("Expected: (dog,(cat,rabbit),(rat,(elephant,mouse)));")
println("Both describe one unrooted tree — they differ only in which branch is")
println("written first, which an unrooted tree does not fix.")

fn test_chbp_character_based_phylogeny() {
    # The real requirement is that the tree induces exactly the table's splits.
    # Newick is not canonical for an unrooted tree, so comparing strings would be
    # comparing a formatting choice.
    fn splits_of(tree_text, names) {
        # Every parenthesised group is a clade; take each as a split.
        let found = []
        let stack = []
        for i in range(0, len(tree_text)) {
            let symbol = substr(tree_text, i, 1)
            if symbol == "(" { stack = push(stack, i) }
            if symbol == ")" {
                let opened = stack[len(stack) - 1]
                stack = slice(stack, 0, len(stack) - 1)
                let inner_text = substr(tree_text, opened + 1, i - opened - 1)
                let members = names |> filter(|t| contains(inner_text, t))
                if len(members) > 1 and len(members) < len(names) {
                    found = push(found, join(sort(members), ","))
                }
            }
        }
        unique(found)
    }

    let mine = splits_of(newick, taxa)
    let published = splits_of("(dog,(cat,rabbit),(rat,(elephant,mouse)));", taxa)

    # A split and its complement are the same split, so compare canonically.
    fn canonical_splits(raw, names) {
        raw |> map(|s| {
            let members = split(s, ",")
            let other = names |> filter(|t| contains(members, t) == false)
            let one = join(sort(members), ",")
            let two = join(sort(other), ",")
            if one < two then one else two
        }) |> unique() |> sort()
    }

    assert canonical_splits(mine, taxa) == canonical_splits(published, taxa),
        "CHBP: splits " + join(canonical_splits(mine, taxa), " | ")
            + " against " + join(canonical_splits(published, taxa), " | ")
    # And those splits are exactly the ones the characters describe.
    let from_characters = canonical_splits(
        clades |> map(|c| join(sort(c), ",")), taxa)
    assert canonical_splits(mine, taxa) == from_characters,
        "CHBP: the tree must induce the table's splits and no others"
}

# Rosalind: EUBT — Enumerating Unrooted Binary Trees
# https://rosalind.info/problems/eubt/
#
# Given: n taxa.
# Return: Every unrooted binary tree on those taxa, in Newick format.

let taxa = ["dog", "cat", "mouse", "elephant"]

# Built by insertion: start with the only tree on three taxa — a star — then add
# each remaining taxon by splitting one existing edge in two and hanging it off
# the new node. Every unrooted binary tree arises exactly once this way, which is
# why the count is the double factorial (2n-5)!! and why it explodes: fifteen
# trees at n=5, over two million at n=10.
#
# That growth is the reason nobody enumerates trees to find the best one. It is
# also why the parsimony and distance methods in this pack exist.
let edges_start = [
    { a: "center", b: taxa[0] },
    { a: "center", b: taxa[1] },
    { a: "center", b: taxa[2] },
]

let trees = [edges_start]
for step in range(3, len(taxa)) {
    let leaf = taxa[step]
    let grown = []
    for tree in trees {
        for i in range(0, len(tree)) {
            let fresh = "node" + str(step) + "_" + str(i)
            let split_edge = tree[i]
            let others = range(0, len(tree)) |> filter(|j| j != i) |> map(|j| tree[j])
            grown = push(grown, others + [
                { a: split_edge.a, b: fresh },
                { a: fresh, b: split_edge.b },
                { a: fresh, b: leaf },
            ])
        }
    }
    trees = grown
}

fn neighbours_of(tree, node) {
    (tree |> filter(|e| e.a == node) |> map(|e| e.b))
        + (tree |> filter(|e| e.b == node) |> map(|e| e.a))
}

fn is_taxon(node, names) { contains(names, node) }

fn walk(tree, node, parent, names) {
    if is_taxon(node, names) { return node }
    let children = neighbours_of(tree, node) |> filter(|c| c != parent)
    "(" + (children |> map(|c| walk(tree, c, node, names)) |> join(",")) + ")"
}

# Written from the first taxon, which an unrooted tree does not privilege — it
# is just somewhere to start reading.
fn to_newick(tree, names) {
    let anchor = names[0]
    let hub = neighbours_of(tree, anchor)[0]
    let branches = neighbours_of(tree, hub) |> filter(|c| c != anchor)
    "(" + anchor + "," + (branches |> map(|c| walk(tree, c, hub, names)) |> join(",")) + ");"
}

let written = trees |> map(|t| to_newick(t, taxa))

println("Result:")
for line in written { println("  " + line) }
println("Expected: three trees — (mouse,cat)|(elephant,dog), (elephant,mouse)|(cat,dog),")
println("          (elephant,cat)|(mouse,dog), written from a different anchor.")

fn test_eubt_enumerating_unrooted_binary_trees() {
    # (2n-5)!! trees: 1 x 3 for n = 4.
    assert len(trees) == 3, "EUBT: expected 3 trees, got " + str(len(trees))
    assert len(unique(written)) == 3, "EUBT: the trees must be distinct"

    # A tree on n taxa has 2n-3 edges and n-2 internal nodes.
    for tree in trees {
        assert len(tree) == 2 * len(taxa) - 3,
            "EUBT: expected " + str(2 * len(taxa) - 3) + " edges, got " + str(len(tree))
        for taxon in taxa {
            assert len(neighbours_of(tree, taxon)) == 1, "EUBT: " + taxon + " must be a leaf"
        }
    }

    # The three topologies are exactly the three ways of pairing four taxa. That
    # is the whole content of the answer, and it does not depend on how the
    # Newick is rooted.
    fn pairing_of(tree, names) {
        # The one internal edge separates the taxa into two pairs. Which side the
        # edge happens to be stored from is arbitrary, so the smaller of the two
        # names the topology.
        let internal = tree |> filter(|e|
            is_taxon(e.a, names) == false and is_taxon(e.b, names) == false)
        let side = neighbours_of(tree, internal[0].a) |> filter(|x| is_taxon(x, names))
        let other = names |> filter(|t| contains(side, t) == false)
        let one = join(sort(side), ",")
        let two = join(sort(other), ",")
        if one < two then one else two
    }
    let pairings = trees |> map(|t| pairing_of(t, taxa)) |> sort()
    assert len(unique(pairings)) == 3, "EUBT: each tree must pair the taxa differently"
    # With four taxa the topology is fixed by which one cat is paired with, and
    # all three possibilities appear exactly once.
    assert pairings == sort(["cat,dog", "cat,mouse", "cat,elephant"]),
        "EUBT: got pairings " + join(pairings, " | ")
}

# Rosalind: QRTD — Quartet Distance
# https://rosalind.info/problems/qrtd/
#
# Given: n taxa and two unrooted binary trees.
# Return: The quartet distance between them.

let taxa = ["A", "B", "C", "D", "E"]
# (A,C,((B,D),E));   and   (C,(B,D),(A,E));
# Each internal edge is a split; these are the non-trivial sides.
let first_clades = [["B", "D"], ["B", "D", "E"]]
let second_clades = [["B", "D"], ["A", "E"]]

# Two trees can share every taxon and still disagree about almost everything, so
# comparing them needs a measure. Quartets give one: each set of four taxa is
# separated into two pairs by exactly one edge, and two trees either agree about
# that pairing or they do not.
#
# It beats counting shared splits because it degrades gracefully — moving one
# taxon changes a handful of quartets, where it can destroy every split at once.
fn pairing_in(four, clades) {
    # The clade cutting these four 2-and-2 names the quartet.
    let cutting = clades |> filter(|c| (four |> count_if(|t| contains(c, t))) == 2)
    if len(cutting) == 0 { return "" }
    let side = four |> filter(|t| contains(cutting[0], t))
    let other = four |> filter(|t| contains(cutting[0], t) == false)
    let one = join(sort(side), ",")
    let two = join(sort(other), ",")
    if one < two then one + "|" + two else two + "|" + one
}

fn subsets_of_four(items) {
    range(0, len(items)) |> flat_map(|a|
        range(a + 1, len(items)) |> flat_map(|b|
            range(b + 1, len(items)) |> flat_map(|c|
                range(c + 1, len(items)) |> map(|d|
                    [items[a], items[b], items[c], items[d]]))))
}

let quartets = subsets_of_four(taxa)
let shared = quartets |> count_if(|four| {
    let one = pairing_in(four, first_clades)
    let two = pairing_in(four, second_clades)
    one != "" and one == two
})

# Both trees are fully resolved, so each induces exactly one quartet per subset.
let first_count = quartets |> count_if(|four| pairing_in(four, first_clades) != "")
let second_count = quartets |> count_if(|four| pairing_in(four, second_clades) != "")
let distance = first_count + second_count - 2 * shared

println("Result:   " + str(distance))
println("Expected: 4")
println("(" + str(len(quartets)) + " quartets each, " + str(shared) + " agreeing)")

fn test_qrtd_quartet_distance() {
    assert distance == 4, "QRTD: got " + str(distance)
    # A resolved tree resolves every four-taxon subset, which CNTQ established.
    assert first_count == choose(len(taxa), 4), "QRTD: the first tree resolves all 5"
    assert second_count == choose(len(taxa), 4), "QRTD: and so does the second"
    assert shared == 3, "QRTD: 3 of the 5 quartets agree"

    # A tree against itself is distance zero — the identity any distance must
    # satisfy, and worth checking rather than assuming.
    let self_shared = quartets |> count_if(|four| pairing_in(four, first_clades) != "")
    assert first_count + first_count - 2 * self_shared == 0,
        "QRTD: a tree must be distance 0 from itself"
    # The two disagreeing quartets are the ones involving both B,D and A,E.
    let disagreeing = quartets |> filter(|four|
        pairing_in(four, first_clades) != pairing_in(four, second_clades))
    assert len(disagreeing) == 2, "QRTD: exactly two quartets differ"
}

# Rosalind: RSUB — Identifying Reversing Substitutions
# https://rosalind.info/problems/rsub/
#
# Given: A rooted binary tree with every node labelled by a string.
# Return: Every reversing substitution — a change that later changes back, with
# no other change in between.

# (((ostrich,cat)rat,mouse)dog,elephant)robot;
let children = {
    "robot": ["dog", "elephant"],
    "dog": ["rat", "mouse"],
    "rat": ["ostrich", "cat"],
}
let sequences = {
    "robot": "AATTG", "dog": "GGGCA", "mouse": "AAGAC", "rat": "GTTGT",
    "cat": "GAGGC", "ostrich": "GTGTC", "elephant": "AATTC",
}
let root = "robot"

# A site that mutates and later mutates back looks, from the tips alone, as
# though nothing ever happened — the ancestor and the descendant agree. Only the
# internal labels reveal it, which is why this problem hands them over rather
# than asking for them.
#
# It matters because reversing substitutions are exactly what makes distant
# relationships hard to recover: the signal erases itself, and two lineages look
# more similar than their history warrants.
fn kids_of(node) { if has_key(children, node) then children[node] else [] }

# Walk down from a node where a change occurred, following only lineages that
# still carry the new character, and report wherever it changes back.
fn reversions_below(start, position, was, became) {
    let found = []
    let frontier = kids_of(start)
    while len(frontier) > 0 {
        let node = frontier[0]
        frontier = slice(frontier, 1, len(frontier))
        let here = substr(sequences[node], position, 1)
        if here == was {
            # Changed back, with nothing else in between.
            found = push(found, node)
        } else {
            # Still carrying the substitution; keep descending. Any third
            # character ends the lineage's relevance.
            if here == became { frontier = frontier + kids_of(node) }
        }
    }
    found
}

let width = len(sequences[root])
let reported = []
let stack = [root]
while len(stack) > 0 {
    let parent = stack[len(stack) - 1]
    stack = slice(stack, 0, len(stack) - 1)
    for child in kids_of(parent) {
        stack = push(stack, child)
        for i in range(0, width) {
            let was = substr(sequences[parent], i, 1)
            let became = substr(sequences[child], i, 1)
            if was != became {
                for reverted in reversions_below(child, i, was, became) {
                    reported = push(reported, child + " " + reverted + " " + str(i + 1)
                        + " " + was + "->" + became + "->" + was)
                }
            }
        }
    }
}

println("Result:")
for line in sort(reported) { println("  " + line) }
println("Expected (any order): dog mouse 1 A->G->A / dog mouse 2 A->G->A")
println("                      rat ostrich 3 G->T->G / rat cat 3 G->T->G / dog rat 3 T->G->T")

fn test_rsub_reversing_substitutions() {
    let expected = ["dog mouse 1 A->G->A", "dog mouse 2 A->G->A",
                    "rat ostrich 3 G->T->G", "rat cat 3 G->T->G", "dog rat 3 T->G->T"]
    assert sort(reported) == sort(expected), "RSUB: got " + join(sort(reported), " / ")

    # Every report must describe a real reversion: the character genuinely
    # differs from the parent and genuinely returns at the named descendant.
    for line in reported {
        let parts = split(line, " ")
        let changed = parts[0]
        let reverted = parts[1]
        let position = int(parts[2]) - 1
        assert substr(sequences[changed], position, 1) != substr(sequences[reverted], position, 1),
            "RSUB: " + line + " does not actually revert"
    }
    # The tips alone hide this: robot and mouse agree at position 1 despite two
    # substitutions on the path between them.
    assert substr(sequences["robot"], 0, 1) == substr(sequences["mouse"], 0, 1),
        "RSUB: the endpoints agree, which is what makes the change invisible"
    assert substr(sequences["dog"], 0, 1) != substr(sequences["robot"], 0, 1),
        "RSUB: yet the intermediate differs from both"
}