Doomsday Fuel

Created: May 01, 2024

Last modified: Oct 11, 2026

TL, DR.

Show challenge description

Welcome to level 3. The LAMBCHOP's fuel comes from a really weird ore that, while it's being processed, keeps switching forms at random until it lands on a stable one, and there it stays. We get a matrix where m[i][j] counts how many times state i was seen turning into state j. Everything starts at s0, and we have to say how likely it is to end up in each stable state, as exact fractions over a common denominator.

If you've ever heard of Markov chains, this is the textbook problem: states, transition probabilities that never change, and a few states you can never leave. If you haven't, no worries, we'll take it one step at a time.

# Reading the matrix

The rows hold counts, not probabilities, but converting them is easy: divide each number by its row's total. In the statement's example, s1 is [4, 0, 0, 3, 2, 0], which adds up to 9, so it moves to s0 with probability 4/9, to s3 with 3/9, and to s4 with 2/9.

What about the stable states? They're the rows full of zeros: nobody ever saw them turn into anything else. In Markov lingo they're called absorbing states, and the rest are transient.

# The trouble is the loops

If the ore only ever moved forward, we could just multiply probabilities along each path. But in the example, s0 can go to s1 and s1 can go back to s0, as many times as it likes. There are infinitely many paths.

Let's look at s5. From s0 you get there directly with probability 1/2. But you could also take one full lap s0 → s1 → s0 (probability $\tfrac{1}{2} \cdot \tfrac{4}{9} = \tfrac{2}{9}$) and then head to s5, or take two laps, or three...

A geometric series! Infinitely many paths, yet they add up to an exact number. With a single loop you can work it out by hand; with several states tangled together, you need the matrix version.

# Same idea, with matrices

Sort the states so the transient ones are separated from the absorbing ones, and keep two pieces of the probability matrix:

$Q^k$ tells you how likely you are to still be wandering around after $k$ steps. Add up every possible number of steps and you get the same geometric series as before, just with matrices:

$N$ is called the fundamental matrix: $N_{ij}$ is how many times, on average, you pass through state $j$ when starting from $i$. And since every visit to a transient state is a chance to jump into an absorbing one, multiplying by $R$ gives exactly what we're after:

$B_{ij}$ is the probability of ending up in absorbing state $j$ when starting from $i$. We always start at s0, so the answer is the first row of $B$.

# Drop some ore samples

Here are both examples from the statement as graphs: transient states on top, stable ones at the bottom (with a thick border). Each sample starts at s0 and hops around at random following the arrows until it gets stuck. Drop a single one to see its path, or a thousand at once: the bars show what fraction ended up in each state, and the black tick marks the exact probability. How many samples does it take for the bars to hit their marks?

No samples dropped yet.

# The code

My solution is the recipe above, written by hand (Foobar didn't allow numpy). I based it on this article on absorption probabilities.

# https://complex-systems-ai.com/en/markov-process/probability-of-state-absorption/
import copy
from functools import reduce
from fractions import Fraction

def inverse(mat):
    if len(mat) != len(mat[0]):
        return 0
    aux_mat = copy.deepcopy(mat)
    mat = [[Fraction(int(i == j)) for i in range(len(aux_mat))]
            for j in range(len(aux_mat))]
    for pivot, val in enumerate(aux_mat):
        pivotVal = aux_mat[pivot][pivot]
        for j, (a, b) in enumerate(zip(val, mat[pivot])):
            aux_mat[pivot][j] = Fraction(a, pivotVal)
            mat[pivot][j] = Fraction(b, pivotVal)
        for i in range(len(aux_mat)):
            if i != pivot:
                aux_piv = aux_mat[i][pivot]
                for j in range(len(aux_mat)):
                    aux_mat[i][j] = aux_mat[i][j] - aux_piv * aux_mat[pivot][j]
                    mat[i][j] = mat[i][j] - aux_piv * mat[pivot][j]
    return mat

def matrix_multiplication(a, b):
    zip_b = list(zip(*b))
    return [[sum(ele_a * ele_b for ele_a, ele_b in zip(row_a, col_b))
             for col_b in zip_b] for row_a in a]

def get_r_qt(mat):
    """
    This function convert the input matrix to Markov cannonical form in order to get the R matrix and (I-Q) matrix
             I | 0
        S = ---+---
             R | Q
        Qt = (I - Q)
    :param mat:
    :return:
    """
    absorb_states = {
        index for index, state in enumerate(mat) if sum(state) == 0
    }
    q = list()
    r = list()
    for i, states in enumerate(mat):
        if i in absorb_states:
            continue
        denominator = sum(states)
        _q = list()
        _r = list()
        for j, state in enumerate(states):
            current_fraction = Fraction(state, denominator)
            if j in absorb_states:
                _r.append(current_fraction)
                continue
            if i == j:
                _q.append(1 - current_fraction)
            else:
                _q.append(-current_fraction)
        q.append(_q)
        r.append(_r)
    return r, q

def lcmm(l):
    def gcd(a, b):
        while b:
            a, b = b, a % b
        return a
    return reduce(lambda x, y: (x * y) // gcd(x, y), l)

def solution(m):
    if len(m) < 2:
        return [1, 1]
    # s0 is already stable: it ends there with probability 1
    if sum(m[0]) == 0:
        terminals = sum(1 for row in m if sum(row) == 0)
        return [1] + [0] * (terminals - 1) + [1]
    r, q = get_r_qt(m)
    sol = matrix_multiplication(inverse(q), r)[0]
    lcm = lcmm([i.denominator for i in sol])
    output = [lcm // i.denominator * i.numerator for i in sol] + [lcm]
    return output

# The case that almost slips through

What if s0 is already stable? Say, [[0, 0, 0], [1, 0, 1], [0, 0, 0]]. The obvious answer is that the ore just stays put: [1, 0, 1]. But the matrix recipe has no idea. s0 is absorbing, so it has no row in $Q$ or $R$, and the "first row of $B$" ends up being s1's: it returns [1, 1, 2], the odds for a sample that started somewhere else. That's why solution starts with a check: if s0's row adds up to zero, it answers right away.

# Wrapping up

Infinitely many paths sounds like an impossible problem, until you notice they repeat in a pattern and add up to a geometric series. Absorbing Markov chains are exactly that, in matrix form: $(I - Q)^{-1}$ sums every lap in one go. Next time you see "this could repeat forever", check whether there's a series hiding in there.