#!/usr/bin/python2

# Oppsieee https://github.com/mkutny/absorbing-markov-chains

# Doomsday Fuel
# =============

# Making fuel for the LAMBCHOP's reactor core is a tricky process because of
# the exotic matter involved. It starts as raw ore, then during processing,
# begins randomly changing between forms, eventually reaching a stable form.
# There may be multiple stable forms that a sample could ultimately reach, not
# all of which are useful as fuel.

# Commander Lambda has tasked you to help the scientists increase fuel creation
# efficiency by predicting the end state of a given ore sample. You have
# carefully studied the different structures that the ore can take and which
# transitions it undergoes. It appears that, while random, the probability of
# each structure transforming is fixed. That is, each time the ore is in 1
# state, it has the same probabilities of entering the next state (which might
# be the same state). You have recorded the observed transitions in a matrix.
# The others in the lab have hypothesized more exotic forms that the ore can
# become, but you haven't seen all of them.

# Write a function solution(m) that takes an array of array of nonnegative ints
# representing how many times that state has gone to the next state and return
# an array of ints for each terminal state giving the exact probabilities of
# each terminal state, represented as the numerator for each state, then the
# denominator for all of them at the end and in simplest form. The matrix is at
# most 10 by 10. It is guaranteed that no matter which state the ore is in,
# there is a path from that state to a terminal state. That is, the processing
# will always eventually end in a stable state. The ore starts in state 0. The
# denominator will fit within a signed 32-bit integer during the calculation,
# as long as the fraction is simplified regularly.

# For example, consider the matrix m:
# [
#   [0,1,0,0,0,1],  # s0, the initial state, goes to s1 and s5 with equal probability
#   [4,0,0,3,2,0],  # s1 can become s0, s3, or s4, but with different probabilities
#   [0,0,0,0,0,0],  # s2 is terminal, and unreachable (never observed in practice)
#   [0,0,0,0,0,0],  # s3 is terminal
#   [0,0,0,0,0,0],  # s4 is terminal
#   [0,0,0,0,0,0],  # s5 is terminal
# ]
# So, we can consider different paths to terminal states, such as:
# s0 -> s1 -> s3
# s0 -> s1 -> s0 -> s1 -> s0 -> s1 -> s4
# s0 -> s1 -> s0 -> s5
# Tracing the probabilities of each, we find that
# s2 has probability 0
# s3 has probability 3/14
# s4 has probability 1/7
# s5 has probability 9/14
# So, putting that together, and making a common denominator, gives an answer in the form of
# [s2.numerator, s3.numerator, s4.numerator, s5.numerator, denominator] which is
# [0, 3, 2, 9, 14].


from fractions import Fraction
from fractions import gcd

def num_of_transients(m):
    if len(m) == 0:
        raise Exception("Can't get transient states of empty matrix")

    for r in range(len(m)):
        for c in range(len(m[r])):
            if m[r][c] != 0:
                break
        else:
            return r
    raise Exception("Not a valid AMC matrix: no absorbing (terminal) states")

def decompose(m):
    t = num_of_transients(m)
    if t == 0:
        raise Exception("No transient states. At least initial state is needed.")

    Q = []
    for r in range(t):
        qRow = []
        for c in range(t):
            qRow.append(m[r][c])
        Q.append(qRow)
    if Q == []:
        raise Exception("Not a valid AMC matrix: no transient states")

    R = []
    for r in range(t):
        rRow = []
        for c in range(t, len(m[r])):
            rRow.append(m[r][c])
        R.append(rRow)
    if R == []:
        raise Exception("Not a valid AMC matrix: missing absorbing states")
    return Q, R

def identity(t):
    m = []
    for i in range(t):
        r = []
        for j in range(t):
            r.append(int(i == j))
        m.append(r)
    return m

def isZero(m):
    for r in range(len(m)):
        for c in range(len(m[r])):
            if m[r][c] != 0:
                return False
    return True

def swap(m, i, j):
    n = []
    s = len(m)

    if s != len(m[0]):
        raise Exception("Cannot swap non-square matrix")

    if i == j:
        return m

    for r in range(s):
        nRow = []
        tmpRow = m[r]
        if r == i:
            tmpRow = m[j]
        if r == j:
            tmpRow = m[i]
        for c in range(s):
            tmpEl = tmpRow[c]
            if c == i:
                tmpEl = tmpRow[j]
            if c == j:
                tmpEl = tmpRow[i]
            nRow.append(tmpEl)
        n.append(nRow)
    return n

def sort(m):
    size = len(m)

    zero_row = -1
    for r in range(size):
        sum = 0
        for c in range(size):
            sum += m[r][c]
        if sum == 0:
            zero_row = r
        if sum != 0 and zero_row > -1:
            n = swap(m, r, zero_row)
            return sort(n)
    return m

def normalize(m, use_fractions=False ):
    n = []
    for r in range(len(m)):
        sum = 0
        cols = len(m[r])
        for c in range(cols):
            sum += m[r][c]

        nRow = []

        if sum == 0:
            nRow = m[r]
        else:
            for c in range(cols):
                if use_fractions:
                    nRow.append(Fraction(m[r][c], sum))
                else:
                    nRow.append(float(m[r][c])/sum)
        n.append(nRow)
    return n

def subtract(i, q):
    if len(i) != len(i[0]) or len(q) != len(q[0]):
        raise Exception("non-square matrices")

    if len(i) != len(q) or len(i[0]) != len(q[0]):
        raise Exception("Cannot subtract matrices of different sizes")

    s = []
    for r in range(len(i)):
        sRow = []
        for c in range(len(i[r])):
            sRow.append(i[r][c] - q[r][c])
        s.append(sRow)
    return s

def multiply(a, b):
    if a == [] or b == []:
        raise Exception("Cannot multiply empty matrices")

    if len(a[0]) != len(b):
        raise Exception("Cannot multiply matrices of incompatible sizes")

    m = []
    rows = len(a)
    cols = len(b[0])
    iters = len(a[0])

    for r in range(rows):
        mRow = []
        for c in range(cols):
            sum = 0
            for i in range(iters):
                sum += a[r][i]*b[i][c]
            mRow.append(sum)
        m.append(mRow)
    return m

def transposeMatrix(m):
    t = []
    for r in range(len(m)):
        tRow = []
        for c in range(len(m[r])):
            if c == r:
                tRow.append(m[r][c])
            else:
                tRow.append(m[c][r])
        t.append(tRow)
    return t

def getMatrixMinor(m,i,j):
    return [row[:j] + row[j+1:] for row in (m[:i]+m[i+1:])]

def getMatrixDeternminant(m):
    if len(m) == 2:
        return m[0][0]*m[1][1]-m[0][1]*m[1][0]

    d = 0
    for c in range(len(m)):
        d += ((-1)**c)*m[0][c]*getMatrixDeternminant(getMatrixMinor(m,0,c))

    return d

def getMatrixInverse(m):
    d = getMatrixDeternminant(m)

    if d == 0:
        raise Exception("Cannot get inverse of matrix with zero determinant")

    if len(m) == 2:
        return [[m[1][1]/d, -1*m[0][1]/d],
                [-1*m[1][0]/d, m[0][0]/d]]

    cofactors = []
    for r in range(len(m)):
        cofactorRow = []
        for c in range(len(m)):
            minor = getMatrixMinor(m,r,c)
            cofactorRow.append(((-1)**(r+c)) * getMatrixDeternminant(minor))
        cofactors.append(cofactorRow)
    cofactors = transposeMatrix(cofactors)
    for r in range(len(cofactors)):
        for c in range(len(cofactors)):
            cofactors[r][c] = cofactors[r][c]/d
    return cofactors


def find_lcm(a):
    lcm = a[0]
    for i in a[1:]:
        lcm = lcm*i//gcd(lcm, i)
    return lcm
    
def simplify(b):
    a = [i.denominator for i in b]
    lcm = find_lcm(a)

    c = []
    for i in range(len(b)):
        c.append(lcm * b[i].numerator // b[i].denominator)
    c.append(lcm)
    return c

def solution(m):
    m = sort(m)
    n = normalize(m,use_fractions=True)
    (q, r) = decompose(n)
    i = identity(len(q))
    s = subtract(i, q)
    v = getMatrixInverse(s)
    b = multiply(v, r)

    return simplify(b[0])

print(solution([[0, 2, 1, 0, 0],
                [0, 0, 0, 3, 4],
                [0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0]]))
#   [7, 6, 8, 21]

print(solution([[0, 1, 0, 0, 0, 1],
                [4, 0, 0, 3, 2, 0],
                [0, 0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0, 0]]))
#    [0, 3, 2, 9, 14]