Hmm

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

# Rosalind: BA10A — Compute the Probability of a Hidden Path
# https://rosalind.info/problems/ba10a/
#
# Given: A hidden path pi, the states of an HMM, and its transition matrix.
# Return: Pr(pi).

let hidden_path = "AABBBAABABAAAABBBBAABBABABBBAABBAAAABABAABBABABBAB"

# No emissions are involved, so the model needs only its states and how they
# follow one another.
let model = {
    states: ["A", "B"],
    transition: {
        A: { A: 0.194, B: 0.806 },
        B: { A: 0.273, B: 0.727 },
    },
}

# Every state is equally likely to start, then each step multiplies in one
# transition. Fifty of them take the answer down to 1e-19, which is what makes
# the log-space treatment in BA10C and BA10D necessary rather than fussy.
let probability = hmm_path_probability(hidden_path, model)

# Printed as a mantissa and an exponent; the plain decimal expansion of 1e-19
# is unreadable next to the published answer.
println("Result:   " + str(round(probability * 1e19, 6)) + "e-19")
println("Expected: 5.017329e-19")

fn test_ba10a_hidden_path_probability() {
    assert abs(probability - 5.01732865318e-19) < 1e-30,
        "BA10A: got " + str(probability)
    # The same thing computed by hand, to check the builtin agrees with the
    # definition rather than only with itself.
    let steps = chars(hidden_path)
    let by_hand = range(1, len(steps))
        |> reduce(|running, i| running * model.transition[steps[i - 1]][steps[i]], 0.5)
    assert abs(probability - by_hand) < 1e-30, "BA10A: builtin and hand computation disagree"
}

# Rosalind: BA10B — Compute the Probability of an Outcome Given a Hidden Path
# https://rosalind.info/problems/ba10b/
#
# Given: A string x, its alphabet, a hidden path pi, the states, and the
# emission matrix.
# Return: Pr(x | pi).

let observed = "xxyzyxzzxzxyxyyzxxzzxxyyxxyxyzzxxyzyzxzxxyxyyzxxzx"
let hidden_path = "BBBAAABABABBBBBBAAAAAABAAAABABABBBBBABAABABABABBBB"

# The path is given, so transitions never come into it — only what each state
# emitted.
let model = {
    states: ["A", "B"],
    symbols: ["x", "y", "z"],
    emission: {
        A: { x: 0.612, y: 0.314, z: 0.074 },
        B: { x: 0.346, y: 0.317, z: 0.336 },
    },
}

# Conditioning on the path makes the positions independent, so this is one
# product of fifty emissions and nothing more.
let probability = hmm_emission_probability(observed, hidden_path, model)

println("Result:   " + str(round(probability * 1e28, 6)) + "e-28")
println("Expected: 1.931571e-28")

fn test_ba10b_outcome_given_path() {
    assert abs(probability - 1.93157070893e-28) < 1e-38,
        "BA10B: got " + str(probability)
    assert len(observed) == len(hidden_path), "BA10B: the sample's string and path are both 50"
    # The same product taken by hand.
    let symbols = chars(observed)
    let states = chars(hidden_path)
    let by_hand = range(0, len(symbols))
        |> reduce(|running, i| running * model.emission[states[i]][symbols[i]], 1.0)
    assert abs(probability - by_hand) < 1e-38, "BA10B: builtin and hand computation disagree"
}

# Rosalind: BA10C — Implement the Viterbi Algorithm
# https://rosalind.info/problems/ba10c/
#
# Given: A string x, the alphabet it was emitted from, the states of an HMM, and
# its transition and emission matrices.
# Return: A path that maximises the probability of x over all hidden paths.

let observed = "xyxzzxyxyy"

# A model is a plain record, so it reads the way the problem states it. The
# matrices are keyed by name in both directions — a transposed transition matrix
# is otherwise a wrong answer rather than an error, and it is the mistake
# everyone makes first.
let model = {
    states: ["A", "B"],
    symbols: ["x", "y", "z"],
    transition: {
        A: { A: 0.641, B: 0.359 },
        B: { A: 0.729, B: 0.271 },
    },
    emission: {
        A: { x: 0.117, y: 0.691, z: 0.192 },
        B: { x: 0.097, y: 0.42,  z: 0.483 },
    },
}

# Viterbi keeps, for each state at each position, the best path reaching it —
# and every path that is not best is dropped immediately. That is what makes it
# linear in the string instead of exponential in it.
let path = viterbi(observed, model) |> join("")

println("Result:   " + path)
println("Expected: AAABBAAAAA")

fn test_ba10c_viterbi() {
    assert path == "AAABBAAAAA", "BA10C: got " + path
    assert len(path) == len(observed), "BA10C: the path must be as long as the string"
    # No other path can score higher — check against the joint probability of the
    # path itself, which is what Viterbi claims to maximise.
    let best = hmm_path_probability(path, model)
             * hmm_emission_probability(observed, path, model)
    let rival = "BBBBBBBBBB"
    let worse = hmm_path_probability(rival, model)
              * hmm_emission_probability(observed, rival, model)
    assert best > worse, "BA10C: an all-B path scores at least as well"
}

# Rosalind: BA10D — Compute the Probability of a String Emitted by an HMM
# https://rosalind.info/problems/ba10d/
#
# Given: A string x, its alphabet, the states of an HMM, and its transition and
# emission matrices.
# Return: Pr(x), summed over every hidden path.

let observed = "xzyyzzyzyy"

let model = {
    states: ["A", "B"],
    symbols: ["x", "y", "z"],
    transition: {
        A: { A: 0.303, B: 0.697 },
        B: { A: 0.831, B: 0.169 },
    },
    emission: {
        A: { x: 0.533, y: 0.065, z: 0.402 },
        B: { x: 0.342, y: 0.334, z: 0.324 },
    },
}

# There are 2^10 paths here and 2^n in general, so summing them one at a time is
# not an option for a real string. The forward algorithm collapses them: once two
# paths reach the same state at the same position, nothing that follows can tell
# them apart, so they can be added together and carried as one number.
let probability = hmm_likelihood(observed, model)

println("Result:   " + str(round(probability * 1e6, 6)) + "e-06")
println("Expected: 1.100551e-06")

fn test_ba10d_string_probability() {
    assert abs(probability - 1.1005510319694847e-06) < 1e-16,
        "BA10D: got " + str(probability)
    # Small enough to enumerate, so check the sum really is over all 1024 paths.
    let states = model.states
    let total = range(0, 1024) |> reduce(|running, mask| {
        let path = range(0, len(observed))
            |> map(|i| states[(mask / pow(2, i)) % 2 |> floor()])
            |> join("")
        running + hmm_path_probability(path, model)
                * hmm_emission_probability(observed, path, model)
    }, 0.0)
    assert abs(probability - total) < 1e-15,
        "BA10D: forward gives " + str(probability) + " but enumeration gives " + str(total)
    # And it must be at least as large as the single best path.
    let best = viterbi(observed, model) |> join("")
    let best_probability = hmm_path_probability(best, model)
                         * hmm_emission_probability(observed, best, model)
    assert probability > best_probability, "BA10D: the sum must exceed its largest term"
}

# Rosalind: BA10J — Solve the Soft Decoding Problem
# https://rosalind.info/problems/ba10j/
#
# Given: A string x, its alphabet, the states of an HMM, and its transition and
# emission matrices.
# Return: For each position, the probability of each state given all of x.

let observed = "zyxxxxyxzz"

let model = {
    states: ["A", "B"],
    symbols: ["x", "y", "z"],
    transition: {
        A: { A: 0.911, B: 0.089 },
        B: { A: 0.228, B: 0.772 },
    },
    emission: {
        A: { x: 0.356, y: 0.191, z: 0.453 },
        B: { x: 0.04,  y: 0.467, z: 0.493 },
    },
}

# Forward-backward, not Viterbi. Viterbi answers "which single path", this
# answers "which state here" — conditioned on the whole string, so a symbol near
# the end can revise a call made near the start. The two disagree in general, and
# the position-wise winners need not even form a path the model can produce.
let posterior = hmm_posterior(observed, model)

println("Result:")
println("  A       B")
for distribution in posterior {
    println("  " + str(round(distribution.A, 4)) + "  " + str(round(distribution.B, 4)))
}
println("Expected first row: 0.5438  0.4562")
println("Expected last row:  0.8167  0.1833")

fn test_ba10j_soft_decoding() {
    let expected_a = [0.5438, 0.6492, 0.9647, 0.9936, 0.9957,
                      0.9891, 0.9154, 0.964,  0.8737, 0.8167]
    assert len(posterior) == len(observed), "BA10J: one row per position"
    for i in range(0, len(expected_a)) {
        assert abs(posterior[i].A - expected_a[i]) < 5e-5,
            "BA10J: position " + str(i) + " gives " + str(posterior[i].A)
                + ", expected " + str(expected_a[i])
        # Each position is a distribution over the states.
        assert abs(posterior[i].A + posterior[i].B - 1.0) < 1e-12,
            "BA10J: row " + str(i) + " does not sum to one"
    }
}

# Rosalind: BA10H — Estimate the Parameters of an HMM
# https://rosalind.info/problems/ba10h/
#
# Given: A string x, its alphabet, a path pi, and the states of an HMM whose
# transition and emission probabilities are unknown.
# Return: The matrices that maximise Pr(x, pi).

let observed = "yzzzyxzxxx"
let hidden_path = "BBABABABAB"

# Only the shape of the model is known, so that is all the skeleton carries. The
# matrices are what comes back.
let skeleton = {
    states: ["A", "B", "C"],
    symbols: ["x", "y", "z"],
}

# With the path given there is nothing to infer: the estimate that maximises
# Pr(x, pi) is just the fraction of the time each transition was taken and each
# symbol emitted. C never appears in the path, so its rows have nothing to count
# and come back uniform — the only choice that leaves them a distribution.
let learned = hmm_estimate(observed, hidden_path, skeleton)

fn show(matrix, rows, columns) {
    println("      " + join(columns, "       "))
    for name in rows {
        println("  " + name + "   " + (columns |> map(|c| str(round(matrix[name][c], 3))) |> join("   ")))
    }
}

println("Transition:")
show(learned.transition, learned.states, learned.states)
println("Emission:")
show(learned.emission, learned.states, learned.symbols)
println("Expected transition row B: 0.8 0.2 0.0")
println("Expected emission row A:   0.25 0.25 0.5")

fn test_ba10h_parameter_estimation() {
    # B is followed by A four times out of five.
    assert abs(learned.transition.B.A - 0.8) < 5e-4, "BA10H: B->A is " + str(learned.transition.B.A)
    assert abs(learned.transition.B.B - 0.2) < 5e-4, "BA10H: B->B is " + str(learned.transition.B.B)
    assert abs(learned.transition.A.B - 1.0) < 5e-4, "BA10H: A->B is " + str(learned.transition.A.B)
    assert abs(learned.emission.A.z - 0.5) < 5e-4, "BA10H: A emits z at " + str(learned.emission.A.z)
    assert abs(learned.emission.B.y - 0.167) < 5e-4, "BA10H: B emits y at " + str(learned.emission.B.y)
    # An unvisited state keeps a usable row rather than a row of NaNs.
    assert abs(learned.transition.C.A - 0.333) < 5e-4, "BA10H: C's row should be uniform"
    # Every row is still a distribution.
    for name in learned.states {
        let total = learned.states |> map(|to| learned.transition[name][to]) |> sum()
        assert abs(total - 1.0) < 1e-9, "BA10H: transition row " + name + " sums to " + str(total)
    }
}

# Rosalind: BA10I — Implement Viterbi Learning
# https://rosalind.info/problems/ba10i/
#
# Given: A number of iterations i, a string x, its alphabet, the states of an
# HMM, and initial transition and emission matrices.
# Return: Matrices that maximise Pr(x, pi) over all matrices and all paths pi.

let iterations = 100
let observed = "xxxzyzzxxzxyzxzxyxxzyzyzyyyyzzxxxzzxzyzzzxyxzzzxyzzxxxxzzzxyyxzzzzzyzzzxxzzxxxyxyzzyxzxxxyxzyxxyzyxz"

# Where this starts is part of the problem, not an implementation detail: the
# procedure climbs to a local optimum, and a different starting point reaches a
# different one.
let start = {
    states: ["A", "B"],
    symbols: ["x", "y", "z"],
    transition: {
        A: { A: 0.582, B: 0.418 },
        B: { A: 0.272, B: 0.728 },
    },
    emission: {
        A: { x: 0.129, y: 0.35,  z: 0.52  },
        B: { x: 0.422, y: 0.151, z: 0.426 },
    },
}

# Decode the most likely path under the current model, then re-estimate the model
# as if that path were the truth — which is BA10H — and repeat. Each round can
# only raise Pr(x, pi), so it settles.
let learned = hmm_viterbi_learning(observed, start, iterations)

println("Transition:")
println("      A       B")
for name in learned.states {
    println("  " + name + "   " + str(round(learned.transition[name].A, 3))
                     + "   " + str(round(learned.transition[name].B, 3)))
}
println("Emission:")
println("      x       y       z")
for name in learned.states {
    println("  " + name + "   " + (learned.symbols |> map(|s| str(round(learned.emission[name][s], 3))) |> join("   ")))
}
println("Expected transition: 0.875 0.125 / 0.011 0.989")
println("Expected emission:   0.0 0.75 0.25 / 0.402 0.174 0.424")

fn test_ba10i_viterbi_learning() {
    assert abs(learned.transition.A.A - 0.875) < 5e-4, "BA10I: A->A is " + str(learned.transition.A.A)
    assert abs(learned.transition.B.B - 0.989) < 5e-4, "BA10I: B->B is " + str(learned.transition.B.B)
    assert abs(learned.emission.A.x - 0.0)   < 5e-4, "BA10I: A emits x at " + str(learned.emission.A.x)
    assert abs(learned.emission.A.y - 0.75)  < 5e-4, "BA10I: A emits y at " + str(learned.emission.A.y)
    assert abs(learned.emission.B.z - 0.424) < 5e-4, "BA10I: B emits z at " + str(learned.emission.B.z)
    # Learning has to explain the data at least as well as the model it started
    # from — that is the property the whole procedure rests on.
    assert hmm_likelihood(observed, learned) >= hmm_likelihood(observed, start),
        "BA10I: learning made the observation less likely"
}

# Rosalind: BA10K — Implement Baum-Welch Learning
# https://rosalind.info/problems/ba10k/
#
# Given: A number of iterations i, a string x, its alphabet, the states of an
# HMM, and initial transition and emission matrices.
# Return: Matrices estimated after i rounds of Baum-Welch learning.

let iterations = 10
let observed = "xzyyzyzyxy"

let start = {
    states: ["A", "B"],
    symbols: ["x", "y", "z"],
    transition: {
        A: { A: 0.019, B: 0.981 },
        B: { A: 0.668, B: 0.332 },
    },
    emission: {
        A: { x: 0.175, y: 0.003, z: 0.821 },
        B: { x: 0.196, y: 0.512, z: 0.293 },
    },
}

# Baum-Welch is expectation-maximisation for an HMM. Where BA10I commits to the
# single best path and counts along it, this counts the *expected* number of
# times each transition was taken across every path at once — which
# forward-backward supplies without enumerating any of them. Keeping the paths it
# would otherwise discard is what makes it the better estimator.
let learned = hmm_baum_welch(observed, start, iterations)

println("Transition:")
println("      A       B")
for name in learned.states {
    println("  " + name + "   " + str(round(learned.transition[name].A, 3))
                     + "   " + str(round(learned.transition[name].B, 3)))
}
println("Emission:")
println("      x       y       z")
for name in learned.states {
    println("  " + name + "   " + (learned.symbols |> map(|s| str(round(learned.emission[name][s], 3))) |> join("   ")))
}
println("Expected transition: 0.0 1.0 / 0.786 0.214")
println("Expected emission:   0.242 0.0 0.758 / 0.172 0.828 0.0")

fn test_ba10k_baum_welch() {
    assert abs(learned.transition.A.B - 1.0)   < 5e-4, "BA10K: A->B is " + str(learned.transition.A.B)
    assert abs(learned.transition.B.A - 0.786) < 5e-4, "BA10K: B->A is " + str(learned.transition.B.A)
    assert abs(learned.emission.A.x - 0.242) < 5e-4, "BA10K: A emits x at " + str(learned.emission.A.x)
    assert abs(learned.emission.A.z - 0.758) < 5e-4, "BA10K: A emits z at " + str(learned.emission.A.z)
    assert abs(learned.emission.B.y - 0.828) < 5e-4, "BA10K: B emits y at " + str(learned.emission.B.y)
    # Each round can only raise Pr(x). Checked round by round rather than only
    # end to end, because a single bad update can be hidden by later good ones.
    let running = range(1, 6) |> map(|n| hmm_likelihood(observed, hmm_baum_welch(observed, start, n)))
    for i in range(1, len(running)) {
        assert running[i] >= running[i - 1] - 1e-12,
            "BA10K: round " + str(i + 1) + " lowered the likelihood"
    }
    assert hmm_likelihood(observed, learned) > hmm_likelihood(observed, start),
        "BA10K: learning should explain the observation better than the start did"
}

# Rosalind: BA10E — Construct a Profile HMM
# https://rosalind.info/problems/ba10e/
#
# Given: A threshold theta, an alphabet, and a multiple alignment.
# Return: The transition and emission probabilities of HMM(Alignment, theta).

let threshold = 0.289
let alphabet = ["A", "B", "C", "D", "E"]
let alignment = [
    "EBA",
    "EBD",
    "EB-",
    "EED",
    "EBD",
    "EBE",
    "E-D",
    "EBD",
]

# A profile HMM turns a family of related sequences into something a new sequence
# can be scored against — it is what Pfam and HMMER search with. Each conserved
# column becomes a match state; the gappy ones become insertions, which is what
# keeps the model's length at the family's length rather than the alignment's.
#
# Here no column is gappy enough to cross the threshold, so all three are match
# columns and the model is three layers long.
let profile = hmm_profile(alignment, alphabet, threshold)

println("States: " + join(profile.states, " "))
println("")
println("Transitions that carry any probability:")
for source in profile.states {
    for target in profile.states {
        if profile.transition[source][target] > 0 {
            println("  " + source + " -> " + target + "   "
                    + str(round(profile.transition[source][target], 3)))
        }
    }
}
println("Expected: S->M1 1.0, M1->M2 0.875, M1->D2 0.125, M2->M3 0.857, M2->D3 0.143")

fn test_ba10e_profile_hmm() {
    assert len(profile.states) == 12, "BA10E: expected 12 states, got " + str(len(profile.states))
    assert profile.states[0] == "S" and profile.states[11] == "E", "BA10E: S and E bracket the model"

    assert abs(profile.transition.S.M1 - 1.0) < 5e-4, "BA10E: S->M1"
    assert abs(profile.transition.M1.M2 - 0.875) < 5e-4, "BA10E: M1->M2"
    assert abs(profile.transition.M1.D2 - 0.125) < 5e-4, "BA10E: M1->D2"
    assert abs(profile.transition.M2.M3 - 0.857) < 5e-4, "BA10E: M2->M3"
    assert abs(profile.transition.M2.D3 - 0.143) < 5e-4, "BA10E: M2->D3"
    assert abs(profile.transition.D2.M3 - 1.0) < 5e-4, "BA10E: D2->M3"
    assert abs(profile.transition.M3.E - 1.0) < 5e-4, "BA10E: M3->E"
    assert abs(profile.transition.D3.E - 1.0) < 5e-4, "BA10E: D3->E"

    assert abs(profile.emission.M1.E - 1.0) < 5e-4, "BA10E: M1 emits E"
    assert abs(profile.emission.M2.B - 0.857) < 5e-4, "BA10E: M2 emits B"
    assert abs(profile.emission.M3.D - 0.714) < 5e-4, "BA10E: M3 emits D"

    # Without pseudocounts, a state the alignment never reaches keeps an empty
    # row instead of a uniform one.
    let out_of_i0 = profile.states |> map(|to| profile.transition.I0[to]) |> sum()
    assert out_of_i0 == 0, "BA10E: I0 is never used, so it should have no transitions"
    # Deletion states are silent.
    let emitted_by_d1 = alphabet |> map(|s| profile.emission.D1[s]) |> sum()
    assert emitted_by_d1 == 0, "BA10E: D1 is a deletion state and cannot emit"
}

# Rosalind: BA10F — Construct a Profile HMM with Pseudocounts
# https://rosalind.info/problems/ba10f/
#
# Given: A threshold theta, a pseudocount sigma, an alphabet, and a multiple
# alignment.
# Return: The transition and emission probabilities of HMM(Alignment, theta, sigma).

let threshold = 0.358
let pseudocount = 0.01
let alphabet = ["A", "B", "C", "D", "E"]
let alignment = [
    "ADA",
    "ADA",
    "AAA",
    "ADC",
    "-DA",
    "D-A",
]

# Without a pseudocount, anything the alignment never happened to do is scored as
# impossible — so a new sequence differing in one position gets probability zero
# rather than a low score. Six sequences cannot rule out the rest of the family,
# and that is what the pseudocount corrects.
#
# It is added after the counts are turned into probabilities, not before: added
# to raw counts its influence would depend on how many sequences the alignment
# happens to contain, which is not something anyone means to tune.
let profile = hmm_profile(alignment, alphabet, threshold, pseudocount)

println("Transitions out of S, I0, M1, D1:")
for source in ["S", "I0", "M1", "D1"] {
    let row = ["I0", "M1", "D1", "I1", "M2", "D2"]
        |> map(|target| target + "=" + str(round(profile.transition[source][target], 3)))
        |> join("  ")
    println("  " + source + ":  " + row)
}
println("Expected S:   I0=0.01  M1=0.819  D1=0.172")
println("Expected I0:  I0=0.333  M1=0.333  D1=0.333")
println("Expected M1:  I1=0.01  M2=0.786")
println("Expected D1:  I1=0.01  M2=0.981")

fn test_ba10f_profile_hmm_with_pseudocounts() {
    assert abs(profile.transition.S.I0 - 0.01)  < 5e-4, "BA10F: S->I0"
    assert abs(profile.transition.S.M1 - 0.819) < 5e-4, "BA10F: S->M1"
    assert abs(profile.transition.S.D1 - 0.172) < 5e-4, "BA10F: S->D1"
    # A row with no counts smooths to uniform over what the topology allows —
    # three states, not all twelve.
    assert abs(profile.transition.I0.I0 - 0.333) < 5e-4, "BA10F: I0->I0"
    assert abs(profile.transition.I0.M1 - 0.333) < 5e-4, "BA10F: I0->M1"
    assert abs(profile.transition.I0.D1 - 0.333) < 5e-4, "BA10F: I0->D1"
    assert abs(profile.transition.M1.I1 - 0.01)  < 5e-4, "BA10F: M1->I1"
    assert abs(profile.transition.M1.M2 - 0.786) < 5e-4, "BA10F: M1->M2"
    assert abs(profile.transition.D1.M2 - 0.981) < 5e-4, "BA10F: D1->M2"

    assert abs(profile.emission.I0.A - 0.2)   < 5e-4, "BA10F: I0 emits A"
    assert abs(profile.emission.M1.A - 0.771) < 5e-4, "BA10F: M1 emits A"
    assert abs(profile.emission.M1.B - 0.01)  < 5e-4, "BA10F: M1 emits B"
    assert abs(profile.emission.M2.D - 0.771) < 5e-4, "BA10F: M2 emits D"

    # A pseudocount does not make a silent state emit, and does not open
    # transitions the topology forbids — smoothing those would invent paths the
    # model does not have.
    let emitted_by_d1 = alphabet |> map(|s| profile.emission.D1[s]) |> sum()
    assert emitted_by_d1 == 0, "BA10F: D1 is a deletion state and cannot emit"
    assert profile.transition.S.M2 == 0, "BA10F: S cannot skip a layer to M2"
    assert profile.transition.M2.M1 == 0, "BA10F: a profile HMM cannot go backwards"
}

# Rosalind: BA10G — Perform a Multiple Sequence Alignment with a Profile HMM
# https://rosalind.info/problems/ba10g/
#
# Given: A string Text, a multiple alignment, a threshold theta, and a
# pseudocount sigma.
# Return: An optimal hidden path emitting Text in HMM(Alignment, theta, sigma).

let text = "AEFDFDC"
let threshold = 0.4
let pseudocount = 0.01
let alphabet = ["A", "B", "C", "D", "E", "F"]
let alignment = [
    "ACDEFACADF",
    "AFDA---CCF",
    "A--EFD-FDC",
    "ACAEF--A-C",
    "ADDEFAAADF",
]

# This is what a profile HMM is for: having learned a family from an alignment,
# align a new sequence to the family rather than to any one of its members.
let profile = hmm_profile(alignment, alphabet, threshold, pseudocount)

# Not `viterbi`, and not because of a naming preference. Deletion states emit
# nothing, so the path is longer than the string it explains — nine states here
# for seven symbols. Ordinary Viterbi advances one state per symbol and cannot
# express that, so the silent states have to be settled in layer order at each
# position before the emitting ones look back at them.
let path = hmm_profile_align(text, profile)

println("Result:   " + join(path, " "))
println("Expected: M1 D2 D3 M4 M5 I5 M6 M7 M8")

fn test_ba10g_align_to_profile() {
    assert join(path, " ") == "M1 D2 D3 M4 M5 I5 M6 M7 M8", "BA10G: got " + join(path, " ")
    # Exactly the emitting states account for the string; the two deletions are
    # the difference between the path's length and the text's.
    let emitting = path |> filter(|state| starts_with(state, "M") or starts_with(state, "I"))
    assert len(emitting) == len(text),
        "BA10G: " + str(len(emitting)) + " emitting states for " + str(len(text)) + " symbols"
    assert len(path) > len(text), "BA10G: the silent states should make the path longer"
}