Solution to Doomsday Fuel foobar challenge

doomsday_fuel.py

from fractions import Fraction

# Replace trials by probabilties of occurrences
def replace_probability(m):        
    for row in range(len(m)):
        total = 0
        for item in range(len(m[row])):
            total += m[row][item]
        if total != 0:
            for item in range(len(m[row])):
                m[row][item] /= float(total)
    return m

# R - non-terminal -> terminal
# Q - non-terminal -> non-terminal
def RQ(m, terminal_state, non_terminal_state):
    R = []
    Q = []
    for i in non_terminal_state:
        temp_t = []
        temp_n = []
        for j in terminal_state:
            temp_t.append(m[i][j])
        for j in non_terminal_state:
            temp_n.append(m[i][j])
        R.append(temp_t)
        Q.append(temp_n)
    return R, Q

# Get Identity Matrix - Q
def subtract_Q_from_identity(Q):
    """
    If Q = [
        [1,2,3],
        [4,5,6],
        [7,8,9],
    ]
    I - Q:
    [[1,0,0]            [[0,-2,-3]
     [0,1,0]   - Q =     [-4,-4,-6]
     [0,0,1]]            [-7,-8,-8]]
    """

    n = len(Q)
    for row in range(len(Q)):
        for item in range(len(Q[row])):
            if row == item:
                Q[row][item] = 1 - Q[row][item]
            else:
                Q[row][item] = -Q[row][item]
    return Q

# Get minor matrix
def get_minor_matrix(Q,i,j):
    """
    Q = [
        [1,2,3],
        [4,5,6],
        [7,8,9],
    ]
    Minor matrix corresponding to 0,0 is
    [
        [5,6],
        [8,9],
    ]
    """

    minor_matrix = []
    for row in Q[:i] + Q[i+1:]:
        temp = []
        for item in row[:j] + row[j+1:]:
            temp.append(item)
        minor_matrix.append(temp)
    return minor_matrix

# Get determinant of a square matrix
def get_determinant(Q):
    if len(Q) == 1:
        return Q[0][0]
    if len(Q) == 2:
        return Q[0][0]*Q[1][1] - Q[0][1]*Q[1][0]
    
    determinant = 0
    for first_row_item in range(len(Q[0])):
        minor_matrix = get_minor_matrix(Q, 0, first_row_item)
        determinant += (((-1)**first_row_item)*Q[0][first_row_item] * get_determinant(minor_matrix))

    return determinant

# Get transpose of a square matrix
def get_transpose_square_matrix(Q):
    for i in range(len(Q)):
        for j in range(i, len(Q), 1):
            Q[i][j], Q[j][i] = Q[j][i], Q[i][j]
    return Q


def get_inverse(Q):
    Q1 = []
    for row in range(len(Q)):
        temp = []
        for column in range(len(Q[row])):
            minor_matrix = get_minor_matrix(Q, row, column)
            determinant = get_determinant(minor_matrix)
            temp.append(((-1)**(row+column))*determinant)
        Q1.append(temp)
    main_determinant = get_determinant(Q)
    Q1 = get_transpose_square_matrix(Q1)
    for i in range(len(Q)):
        for j in range(len(Q[i])):
            Q1[i][j] /= float(main_determinant)
    return Q1

def multiply_matrix(A, B):
    result = []
    dimension = len(A)
    for row in range(len(A)):
        temp = []
        for column in range(len(B[0])):
            product = 0
            for selector in range(dimension):
                product += (A[row][selector]*B[selector][column])
            temp.append(product)
        result.append(temp)
    return result

def gcd(a ,b):
    if b==0:
        return a
    else:
        return gcd(b,a%b)   

def sanitize(M):
    needed = M[0]
    to_fraction = [Fraction(i).limit_denominator() for i in needed]
    lcm = 1
    for i in to_fraction:
        if i.denominator != 1:
            lcm = i.denominator
    for i in to_fraction:
        if i.denominator != 1:
            lcm = lcm*i.denominator/gcd(lcm, i.denominator)
    to_fraction = [(i*lcm).numerator for i in to_fraction]
    to_fraction.append(lcm)
    return to_fraction

def solution(m):
    n = len(m)
    if n==1:
        if len(m[0]) == 1 and m[0][0] == 0:
            return [1, 1]
    terminal_state = []
    non_terminal_state = []

    # Get terminal and non-terminal states
    for row in range(len(m)):
        count = 0
        for item in range(len(m[row])):
            if m[row][item] == 0:
                count += 1
        if count == n:
            terminal_state.append(row)
        else:
            non_terminal_state.append(row)
    # Replace trials by probabilties
    probabilities = replace_probability(m)
    # Get R and Q matrix
    R, Q = RQ(probabilities, terminal_state, non_terminal_state)
    IQ = subtract_Q_from_identity(Q)
    # Get Fundamental Matrix (F)
    IQ1 = get_inverse(IQ)
    product_IQ1_R = multiply_matrix(IQ1, R)
    return sanitize(product_IQ1_R)

# Case where state 0 itself is a terminal state
assert(solution(
    [
        [0],
    ]
)) == [1, 1]

assert(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]

assert(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]

@KaustubhRakhade do you mind clarifying why you divide by the determinant of N in the second to last line of your code?

It's been a while so I don't really remember it clearly.
But the elements of the Probability Matrix (B) are converted from their fractional values to integers by multiplying it by comm_denom / np.linalg.det(N)

This is driving me crazy. My solution gives the correct output for the two given cases, but even those fail. Would someone be able to dissect the reason why the verification is failing, even though the output is exactly as expected when running it locally?

`import numpy as np
import math

from fractions import Fraction

#reference for absorbing Markov Chains: https://math.libretexts.org/Bookshelves/Applied_Mathematics/Applied_Finite_Mathematics_(Sekhon_and_Bloom)/10%3A_Markov_Chains/10.04%3A_Absorbing_Markov_Chains

#code for inverting matrix from https://stackoverflow.com/questions/32114054/matrix-inversion-without-numpy

def transposeMatrix(m):
return map(list,zip(*m))

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

def getMatrixDeternminant(m):
#base case for 2x2 matrix
if len(m) == 2:
return m[0][0]*m[1][1]-m[0][1]*m[1][0]

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

def getMatrixInverse(m):
determinant = getMatrixDeternminant(m)
#special case for 2x2 matrix:
if len(m) == 2:
return [[m[1][1]/determinant, -1m[0][1]/determinant],
[-1m[1][0]/determinant, m[0][0]/determinant]]

#find matrix of cofactors
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]/determinant
return cofactors

def makeCommonDenominator(lcm, fractions):
result = []
for i in range(len(fractions)):
multiplier = lcm/fractions[i].denominator
result.append(int(multiplier*fractions[i].numerator))
result.append(int(lcm))
print(result)
return result

def getLCM(fractions):
denominators = []
for i in range(len(fractions)):
denominators.append(fractions[i].denominator)
lcm = math.lcm(*denominators)
return lcm

def convertToFraction(probabilities):
result = []
for i in range(len(probabilities)):
result.append(Fraction(probabilities[i]))
return result

def extractProbabilities(fa):
num_rows, num_cols = fa.shape
result = []
for i in range(num_cols):
result.append(fa[0,i])
return result

def fa(f, a):
return np.dot(f,a)

def findA(m, numberOfTerminalStates):
num_rows = len(m) - numberOfTerminalStates
a = np.zeros((num_rows, numberOfTerminalStates))
a = a + Fraction()
for i in range(num_rows):
for j in range(numberOfTerminalStates):
a[i,j] = m[i][j+num_rows]
return a

def calculateFundamentalMatrix(b):
num_rows, num_cols = b.shape
i = np.identity(num_rows, int)
i = i + Fraction()
#print(i)
difference = np.subtract(i, b)
#print(difference)
#inverse = np.linalg.inv(difference)
#inverse = difference**(-1)
inverse = getMatrixInverse(difference)
print(inverse)
return inverse

def findB(m, numberOfTerminalStates):
dimensionForB = len(m) - numberOfTerminalStates
b = np.zeros((dimensionForB, dimensionForB))
b = b + Fraction()
for i in range(dimensionForB):
for j in range(dimensionForB):
b[i, j] = m[i][j]
return b

def convertMatrixToProbabilities(m):
result = np.zeros((len(m), len(m)))
result = result + Fraction()

print(result)

for i in range(len(m)):
    denominator = 0
    for j in range(len(m)):
        denominator += m[i][j]
    for j in range(len(m)):
        if(denominator == 0):
            result[i,j] = Fraction(0, 1)
        else:
            result[i,j] = Fraction(m[i][j], denominator) 
return result

def countTerminalStates(m):
numberOfTerminalStates = 0
isTerminal = True
firstTime = True
for row in m:
if(isTerminal and (not firstTime)):
numberOfTerminalStates+=1
isTerminal = True
firstTime = False
for col in row:
if(col != 0):
isTerminal = False
if(isTerminal):
numberOfTerminalStates+=1
return numberOfTerminalStates

def helper(m):
numberOfTerminalStates = countTerminalStates(m)
convertedMatrix = convertMatrixToProbabilities(m)

print(convertedMatrix)

b = findB(convertedMatrix, numberOfTerminalStates)
print(b)
f = calculateFundamentalMatrix(b)
a = findA(convertedMatrix, numberOfTerminalStates)
fA = fa(f,a)
probabilities = extractProbabilities(fA)
fractions = convertToFraction(probabilities)
lcm = getLCM(fractions)
return makeCommonDenominator(lcm, fractions)            

def solution(m):
return helper(m)`

import numpy as np 
import math

from fractions import Fraction

#reference for absorbing Markov Chains: https://math.libretexts.org/Bookshelves/Applied_Mathematics/Applied_Finite_Mathematics_(Sekhon_and_Bloom)/10%3A_Markov_Chains/10.04%3A_Absorbing_Markov_Chains

#code for inverting matrix from https://stackoverflow.com/questions/32114054/matrix-inversion-without-numpy

def transposeMatrix(m):
    return map(list,zip(*m))

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

def getMatrixDeternminant(m):
    #base case for 2x2 matrix
    if len(m) == 2:
        return m[0][0]*m[1][1]-m[0][1]*m[1][0]

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

def getMatrixInverse(m):
    determinant = getMatrixDeternminant(m)
    #special case for 2x2 matrix:
    if len(m) == 2:
        return [[m[1][1]/determinant, -1*m[0][1]/determinant],
                [-1*m[1][0]/determinant, m[0][0]/determinant]]

    #find matrix of cofactors
    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]/determinant
    return cofactors

def makeCommonDenominator(lcm, fractions):
    result = []
    for i in range(len(fractions)):
        multiplier = lcm/fractions[i].denominator
        result.append(int(multiplier*fractions[i].numerator))
    result.append(int(lcm))
    print(result)
    return result

def getLCM(fractions):
    denominators = []
    for i in range(len(fractions)):
        denominators.append(fractions[i].denominator)
    lcm = math.lcm(*denominators)
    return lcm
        

def convertToFraction(probabilities):
    result = []
    for i in range(len(probabilities)):
        result.append(Fraction(probabilities[i]))
    return result

def extractProbabilities(fa):
    num_rows, num_cols = fa.shape
    result = []
    for i in range(num_cols):
        result.append(fa[0,i])
    return result
        

def fa(f, a):
    return np.dot(f,a)

def findA(m, numberOfTerminalStates):
    num_rows = len(m) - numberOfTerminalStates
    a = np.zeros((num_rows, numberOfTerminalStates))
    a = a + Fraction()
    for i in range(num_rows):
        for j in range(numberOfTerminalStates):
            a[i,j] = m[i][j+num_rows]
    return a

def calculateFundamentalMatrix(b):
    num_rows, num_cols = b.shape
    i = np.identity(num_rows, int)
    i = i + Fraction()
    #print(i)
    difference = np.subtract(i, b)
    #print(difference)
    #inverse = np.linalg.inv(difference)
    #inverse = difference**(-1)
    inverse = getMatrixInverse(difference)
    print(inverse)
    return inverse

def findB(m, numberOfTerminalStates):
    dimensionForB = len(m) - numberOfTerminalStates
    b = np.zeros((dimensionForB, dimensionForB))
    b = b + Fraction()
    for i in range(dimensionForB):
        for j in range(dimensionForB):
            b[i, j] = m[i][j]
    return b

def convertMatrixToProbabilities(m):
    result = np.zeros((len(m), len(m)))
    result = result + Fraction()
   # print(result)
    for i in range(len(m)):
        denominator = 0
        for j in range(len(m)):
            denominator += m[i][j]
        for j in range(len(m)):
            if(denominator == 0):
                result[i,j] = Fraction(0, 1)
            else:
                result[i,j] = Fraction(m[i][j], denominator) 
    return result

def countTerminalStates(m): 
    numberOfTerminalStates = 0 
    isTerminal = True
    firstTime = True
    for row in m:
        if(isTerminal and (not firstTime)):
            numberOfTerminalStates+=1
        isTerminal = True
        firstTime = False
        for col in row:
            if(col != 0):
                isTerminal = False
    if(isTerminal):
        numberOfTerminalStates+=1
    return numberOfTerminalStates
    
def helper(m):
    numberOfTerminalStates = countTerminalStates(m)
    convertedMatrix = convertMatrixToProbabilities(m)
   # print(convertedMatrix)
    b = findB(convertedMatrix, numberOfTerminalStates)
    print(b)
    f = calculateFundamentalMatrix(b)
    a = findA(convertedMatrix, numberOfTerminalStates)
    fA = fa(f,a)
    probabilities = extractProbabilities(fA)
    fractions = convertToFraction(probabilities)
    lcm = getLCM(fractions)
    return makeCommonDenominator(lcm, fractions)            

def solution(m):
    return helper(m)

The first comment pasted the code in a really horrible way, please see above for my code in a cleaner format.

Thanks for sharing, currently I am stuck in test case #9, I am not sure why.
and I just saw someone was able to use numpy!!! I swear I was not able to use it in the beginning but It was able when I tried....

really thanks to KaustubhRakhade . I noticed test case 9 is about the special case where the source is terminal!

@kavinyudhitia

I just saw someone was able to use numpy!!! I swear I was not able to use it in the beginning but It was able when I tried....

Yeah, even I was not able to use it. But it turns out you can. Good luck with the submission.

Currently I am only failing test 3, I was decided to not search on foruns directly related to the foobar chalange, but with only one day to go I'm beginnig to get desperate LOL
Anyone has any insight regarding test 3?

Nevermind I got it, it was something on the way I was getting R and Q Matrices

Hi there!
I share my solution, I didnt used matrix stuff like inverse or traspose because I don't understand how to apply this concept to this excercise at all, but think it is useful for someone that is making it in the same way than me:

def solution(m):
    
    #case when 0 is terminal state
    if(not any(m[0])):
        return [1] + ([0] * (len(m)-1)) + [1]
   
    #diagonal values arent relevant 
    cleanDiagonal(m)
    
    #useful structures
    probabilitiesMatrix = generateProbabilityMatrix(m)
    terminals, notTerminals = getTerminalsAndNotTerminals(m)
    
    #we will remove one by one all of the not-terminals nodes 
    for i in notTerminals:
        absorbNode(probabilitiesMatrix, i)
    
    #now we can take the solution
    terminalsProbabilities = list(map(lambda x: probabilitiesMatrix[0][x], terminals))
    commonDenominator = getCommonDenominator(list(map(lambda x: x[1], terminalsProbabilities)))
    unsimplifiedNumerators = list(map(lambda x: fracUnsimplify(x, commonDenominator)[0], terminalsProbabilities))
    
    return unsimplifiedNumerators + [commonDenominator]
    
   
def cleanDiagonal(m):
    for i in range(len(m)):
        m[i][i] = 0


def generateProbabilityMatrix(m):
    result = []
    for i in range(len(m)):
        result.append([None] * len(m))
        for j in range(len(m)):
            result[i][j] = fracDiv([m[i][j],1], [sum(m[i]),1])
    return result
            
            
def getTerminalsAndNotTerminals(m):
    terminals = []
    notTerminals = list(range(1, len(m)))
    for i in range(len(m)):
        if(not any(m[i])):
            terminals.append(i)
            notTerminals.remove(i)
    return terminals, notTerminals
    
    
def absorbNode(pm, node):
    for i in range(len(pm)):
        for j in range(len(pm)):
            if(i != node and j != node):
                pm[i][j] = fracAdd(pm[i][j], fracMult(pm[i][node], pm[node][j]))
                
    for k in range(len(pm)):
        pm[k][node] = [0, 1]
        pm[node][k] = [0, 1]
        
    for i in range(len(pm)):
        if(pm[i][i] != [0, 1]):
            multiplier = solveGeometricSeries(pm[i][i])
            for j in range(len(pm)):
                if(i == j):
                    pm[i][j] = [0, 1]
                else:
                    pm[i][j] = fracMult(pm[i][j] ,multiplier)
                    
                    
#we will work with fractions, so lets create some functions 

def fracSimplify(a):
    if(a[0] == 0):
        a[1] = 1
    i=2
    while (i <= max(a)):
        if(a[0]%i == 0 and a[1]%i == 0):
            a[0] //= i
            a[1] //= i
        else:
            i += 1
    return a
    
def fracAdd(a, b):
    return fracSimplify([a[0]*b[1] + b[0]*a[1] , a[1]*b[1]])
    
def fracSubs(a, b):
    return fracSimplify([a[0]*b[1] - b[0]*a[1] , a[1]*b[1]])
    
def fracMult(a, b):
    return fracSimplify([a[0]*b[0], a[1]*b[1]])

def fracDiv(a, b):
    if(a[1] == 0 or b[1] == 0):
        return [0, 1]
    return fracSimplify([a[0]*b[1], a[1]*b[0]])

def solveGeometricSeries(r):
    if(r == [1,1]):
        return [1,1]
    n = [1,1]
    d = fracSubs([1,1], r)
    return fracDiv(n, d)
    
def getCommonDenominator(l):
    greater = min(l)
    while(not all(list(map(lambda x: greater % x == 0, l)))):
        greater += 1
    return greater

def fracUnsimplify(a, d):
    return [int(a[0]*(d/a[1])), d]

Does anybody know why my code works fine for the two given test cases in my IDE, but it fails all test cases in foobar?
import numpy as np
from fractions import Fraction
#from fractions import gcd #python 2.7
from math import gcd
def solution(m):
#ai represents the chance of starting from state i and going to certain terminal state
'''
[
[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
]
Terminal state 5:
a0 = 0.5 + 0.5a1
a1 = 4/9a0
=> a0 = 9/14, a1 = 2/7
terminal 5 = 9/14
Terminal state 3:
a0 = 0.5a1
a1 = 4/9a0 + 3/9
a0 = 3/14, a1 = 3/7
terminal 3 = 3/14
Terminal state 4:
a0 = 0.5a1
a1 = 4/9a0 + 2/9
a0 = 1/7, a1 = 2/7
Terminal 4: 1/7
ans = 0:3/14:1/7:9/14 = 0:3:2:9
'''
#find all terminal state
def lcm(x, y):
return x * y // gcd(x, y)
Terminal = set()
ans = []
res = []
l = 1 # lcm of demoninator
for i in range(1,len(m)):
isTerminal = True
for j in range(0, len(m[i])):
if m[i][j] !=0:
isTerminal = False
break
if isTerminal == True:
Terminal.add(i)
#convert each line into probability
terminal =[]
for elem in Terminal:
terminal.append(elem)
terminal.sort()
for i in range(0,len(m)):
ttl = 0
if i not in terminal:
for j in range(0, len(m[i])):
ttl += m[i][j]
for j in range(0, len(m[i])):
m[i][j] =m[i][j]/ttl
#loop for each terminal state, list out the linear equation and solve it
lc = 1#lcm
for t in terminal:
#define coefficient matrix and constant vector
vector = []
coei = []
for i in range(0,len(m)):
if i in terminal:#terminal state
continue
con = 0
c = [1] *(len(m)-len(terminal))# current state's coefficient is 1
for j in range(0, len(m[i])):
if j ==i:
continue
if j in terminal:
if j==t:
con = m[i][j]
else:#not terminal state
c[j]=-m[i][j]
vector.append(con)
coei.append(list(c))
y = np.linalg.solve(coei, vector)
# Convert the solution to fractions
y_fraction = [Fraction(num).limit_denominator() for num in y]
numerator = y_fraction[0].numerator
denominator = y_fraction[0].denominator
lc = lcm(lc,denominator)
res.append([numerator,denominator])
for r in res:
a = int(r[0]*lc//r[1])
ans.append(a)
ans.append(int(lc))
return ans

天才，我大概看懂了题解，实在是不知道怎么用保持分数的形式进行运算，比如说求逆