#!/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]
