# 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"
}