diff options
Diffstat (limited to '3.2.py')
| -rw-r--r-- | 3.2.py | 302 |
1 files changed, 302 insertions, 0 deletions
@@ -0,0 +1,302 @@ +#!/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] |