P1 = [[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]]
P2 = [[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]]
import numpy as np
def solution(m):
m = np.matrix(m, dtype=float)
# Diagonal cases are not relevant
for i in range(len(m)):
m[i, i] = 0
active, stable = separate(m)
# Check if first state is stable
if 0 in stable:
return [1] + [0]*(len(m)-2) + [1]
# Take only active states
m = m[active]
# Get common denominator
c_denom = np.lcm.reduce(np.sum(m, axis=1).A1.astype(int))
row_metric = np.sum(m, axis=1)
prob_matrx = m / row_metric
# Take matrix Q
Q = prob_matrx[:, active]
# Take matrix R
R = prob_matrx[:, stable]
# Compute fundamental matrix
F = (np.identity(len(active)) - Q) ** -1
FR = F.dot(R)
# Re-apply metric and correct the inverse distortion
# Found after headaching trial and error
FR = np.multiply(FR, c_denom) / np.linalg.det(F)
print((np.multiply(FR, c_denom) / np.linalg.det(F)).dtype)
# Integer black magic
FR = FR.round().astype(np.int32)
ret = list(FR[0].A1) + [np.sum(FR[0]).astype(np.int32)]
for i in range(len(ret)):
print(type(ret[i]))
return list(FR[0].A1) + [np.sum(FR[0])]
def separate(m):
act_idx, stb_idx = [], []
for i, el in enumerate(m):
if np.any(el):
act_idx.append(i)
else:
stb_idx.append(i)
return act_idx, stb_idx
print(solution(P1))
print(solution(P2))
UDEgPSBbWzAsIDEsIDAsIDAsIDAsIDFdLCBbNCwgMCwgMCwgMywgMiwgMF0sIFswLCAwLCAwLCAwLCAwLCAwXSwgWzAsIDAsIDAsIDAsIDAsIDBdLCBbMCwgMCwgMCwgMCwgMCwgMF0sIFswLCAwLCAwLCAwLCAwLCAwXV0KUDIgPSBbWzAsIDIsIDEsIDAsIDBdLCBbMCwgMCwgMCwgMywgNF0sIFswLCAwLCAwLCAwLCAwXSwgWzAsIDAsIDAsIDAsMF0sIFswLCAwLCAwLCAwLCAwXV0KaW1wb3J0IG51bXB5IGFzIG5wCgpkZWYgc29sdXRpb24obSk6CiAgICBtID0gbnAubWF0cml4KG0sIGR0eXBlPWZsb2F0KQogICAgIyBEaWFnb25hbCBjYXNlcyBhcmUgbm90IHJlbGV2YW50CiAgICBmb3IgaSBpbiByYW5nZShsZW4obSkpOgogICAgICAgIG1baSwgaV0gPSAwCiAgICBhY3RpdmUsIHN0YWJsZSA9IHNlcGFyYXRlKG0pCiAgICAjIENoZWNrIGlmIGZpcnN0IHN0YXRlIGlzIHN0YWJsZQogICAgaWYgMCBpbiBzdGFibGU6CiAgICAgICAgcmV0dXJuIFsxXSArIFswXSoobGVuKG0pLTIpICsgWzFdCiAgICAjIFRha2Ugb25seSBhY3RpdmUgc3RhdGVzCiAgICBtID0gbVthY3RpdmVdCiAgICAjIEdldCBjb21tb24gZGVub21pbmF0b3IKICAgIGNfZGVub20gPSBucC5sY20ucmVkdWNlKG5wLnN1bShtLCBheGlzPTEpLkExLmFzdHlwZShpbnQpKQogICAgcm93X21ldHJpYyA9IG5wLnN1bShtLCBheGlzPTEpCiAgICBwcm9iX21hdHJ4ID0gbSAvIHJvd19tZXRyaWMKICAgICMgVGFrZSBtYXRyaXggUQogICAgUSA9IHByb2JfbWF0cnhbOiwgYWN0aXZlXQogICAgIyBUYWtlIG1hdHJpeCBSCiAgICBSID0gcHJvYl9tYXRyeFs6LCBzdGFibGVdCiAgICAKICAgICMgQ29tcHV0ZSBmdW5kYW1lbnRhbCBtYXRyaXgKICAgIEYgPSAobnAuaWRlbnRpdHkobGVuKGFjdGl2ZSkpIC0gUSkgKiogLTEKICAgIEZSID0gRi5kb3QoUikKICAgICMgUmUtYXBwbHkgbWV0cmljIGFuZCBjb3JyZWN0IHRoZSBpbnZlcnNlIGRpc3RvcnRpb24KICAgICMgRm91bmQgYWZ0ZXIgaGVhZGFjaGluZyB0cmlhbCBhbmQgZXJyb3IKICAgIEZSID0gbnAubXVsdGlwbHkoRlIsIGNfZGVub20pIC8gbnAubGluYWxnLmRldChGKQogICAgcHJpbnQoKG5wLm11bHRpcGx5KEZSLCBjX2Rlbm9tKSAvIG5wLmxpbmFsZy5kZXQoRikpLmR0eXBlKQogICAgIyBJbnRlZ2VyIGJsYWNrIG1hZ2ljCiAgICBGUiA9IEZSLnJvdW5kKCkuYXN0eXBlKG5wLmludDMyKQogICAgcmV0ID0gbGlzdChGUlswXS5BMSkgKyBbbnAuc3VtKEZSWzBdKS5hc3R5cGUobnAuaW50MzIpXQogICAgZm9yIGkgaW4gcmFuZ2UobGVuKHJldCkpOgogICAgCXByaW50KHR5cGUocmV0W2ldKSkKICAgIHJldHVybiBsaXN0KEZSWzBdLkExKSArIFtucC5zdW0oRlJbMF0pXQpkZWYgc2VwYXJhdGUobSk6CiAgICBhY3RfaWR4LCBzdGJfaWR4ID0gW10sIFtdCiAgICBmb3IgaSwgZWwgaW4gZW51bWVyYXRlKG0pOgogICAgICAgIGlmIG5wLmFueShlbCk6CiAgICAgICAgICAgIGFjdF9pZHguYXBwZW5kKGkpCiAgICAgICAgZWxzZToKICAgICAgICAgICAgc3RiX2lkeC5hcHBlbmQoaSkKICAgIHJldHVybiBhY3RfaWR4LCBzdGJfaWR4CnByaW50KHNvbHV0aW9uKFAxKSkKcHJpbnQoc29sdXRpb24oUDIpKQ==