"""Steins;Gate divergence records - analysis script, v10 (main corpus, robustness, supplement, cross-source check,
exploratory cross-work analysis). Readings: Steins;Gate Wiki "List of Known World Lines", cross-checked against the
Japanese strategy wiki, Moegirl, the official Maniacs guidebook and fan-translated game scripts (see paper, Table S9).
Main analysis uses tau_F1 = 170 h (game text); the alternative ordering puts F1 before F2 (tau_F1 = 3 d).
Run: python3 sg_analysis_v10.py"""
import itertools, json, numpy as np
from decimal import Decimal as Dm
import matplotlib as mpl, matplotlib.pyplot as plt
mpl.rcParams.update({'font.family':'Arial','font.size':7,'axes.linewidth':0.6,'xtick.major.width':0.6,'ytick.major.width':0.6,
 'xtick.major.size':2.5,'ytick.major.size':2.5,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':400,'legend.frameon':False})
MM=1/25.4; BLUE='#3B6FB6'; RED='#C8443A'; GREY='#8C8C8C'; GOLD='#D9A21B'
def lab(ax,s,x=-0.18,y=1.04): ax.text(x,y,s,transform=ax.transAxes,fontsize=8,fontweight='bold',va='bottom')
out={}
# ---------- readings (exact decimals)
ops=[('F0','1.130426','0.571024'),('F1','0.571024','0.571015'),('F2','0.571015','0.523299'),('F3','0.523299','0.456903'),
     ('F4','0.456903','0.409420'),('F5','0.409420','0.337187'),('R5','0.337187','0.409431'),('R4','0.409431','0.456914'),
     ('R3','0.456914','0.523307'),('R2','0.523307','0.571046'),('R0','0.571046','1.130205'),('R0b','0.571046','1.130238'),
     ('Sa','0.409431','-0.275349'),('Sc','1.130205','1.048596')]
dD={k:Dm(b)-Dm(a) for k,a,b in ops}; out['dD']={k:str(v) for k,v in dD.items()}
pairs=[('F5','R5'),('F4','R4'),('F3','R3'),('F2','R2'),('F0','R0'),('F0','R0b')]
res={f'{f}/{r}':(dD[f]+dD[r], abs(dD[f]+dD[r])/abs(dD[f])) for f,r in pairs}
out['residuals']={k:[str(e),f'{float(q):.3e}'] for k,(e,q) in res.items()}
rev=[('0.409420','0.409431'),('0.456903','0.456914'),('0.523299','0.523307'),('0.571015','0.571046'),('1.130205','1.130238')]
out['revisits']={f'{a}/{b}':str(Dm(b)-Dm(a)) for a,b in rev}
# ---------- rank statistics
def ranks(x):
    x=np.asarray(x,float); o=np.argsort(x,kind='mergesort'); r=np.empty(len(x)); r[o]=np.arange(1,len(x)+1)
    for v in np.unique(x):  # average ties
        m=x==v; r[m]=r[m].mean()
    return r
def rho(x,y):
    a,b=ranks(x),ranks(y); a-=a.mean(); b-=b.mean(); return float((a*b).sum()/np.sqrt((a*a).sum()*(b*b).sum()))
def perm_dist(n):
    base=np.arange(1,n+1); vals=[round(rho(base,p),4) for p in itertools.permutations(base)]
    u,c=np.unique(vals,return_counts=True); return dict(zip([float(v) for v in u],[int(k) for k in c]))
def pge(r,n): d=perm_dist(n); return sum(c for v,c in d.items() if v>=r-1e-9)/sum(d.values())
out['perm5']=perm_dist(5); out['perm4']=perm_dist(4)
N=['F1','F2','F3','F4','F5']
A=np.array([float(abs(dD[k])) for k in N])           # |dD|
TAU={'F1':None,'F2':4,'F3':17.555*365.25,'F4':10.344*365.25,'F5':1}   # days
RR={'F1':0,'F2':4,'F3':16.934*365.25,'F4':10.344*365.25,'F5':35*365.25}
def tauvec(t1): return np.array([t1 if k=='F1' else TAU[k] for k in N],float)
# (1) F1 tau sweep
sweep=np.round(np.arange(2,7.01,0.25),2); out['tau_sweep']={float(t):round(rho(tauvec(t),A),3) for t in sweep}
# (2) F5 R sweep (days)
Rs=np.logspace(0,np.log10(60*365.25),400)
def Rvec(r5): v=np.array([RR[k] for k in N],float); v[4]=r5; return v
rsw=[rho(Rvec(r),A) for r in Rs]; out['R5_steps']=sorted({round(x,3) for x in rsw})
out['R5_25_35']=sorted({round(rho(Rvec(y*365.25),A),3) for y in np.linspace(25,35,41)})
out['R5_threshold_years']=float(min(r for r,x in zip(Rs,rsw) if x>0.85)/365.25)
# (3) leave-one-out
loo={}
for lbl,t1 in [('tau(F1=170h)',170/24),('tau(F1=3d)',3)]:
    t=tauvec(t1); loo[lbl]={N[i]:round(rho(np.delete(t,i),np.delete(A,i)),3) for i in range(5)}
r=Rvec(RR['F5']); loo['R']={N[i]:round(rho(np.delete(r,i),np.delete(A,i)),3) for i in range(5)}
out['loo']=loo
out['loo_P']={k:{n:round(pge(v,4),3) for n,v in d.items()} for k,d in loo.items()}
out['main']={'tau_main(F1=170h)':[rho(tauvec(170/24),A),pge(-0.3,5)],'tau_alt(F1=3d)':[rho(tauvec(3),A),pge(-0.1,5)],'tau4tie':rho(tauvec(4),A),'R':[rho(Rvec(RR['F5']),A),pge(0.9,5)]}
# (4) epsilon0 / eta definitions
alpha_res=[abs(float(res[k][0])) for k in ['F5/R5','F4/R4','F3/R3','F2/R2']]
alpha_rev=[abs(float(Dm(b)-Dm(a))) for a,b in rev[:4]]
defs={'max α residual':max(alpha_res),'max α revisit':max(alpha_rev),'β context pair':3.3e-5,'rounded (main text)':4e-5,
      'incl. F0/R0 (not same 𝒜)':2.43e-4}
minbig=float(min(abs(dD[k]) for k in ['F2','F3','F4','F5']))
eps_tab={}
for k,e in defs.items():
    eta_need=max(0.0,(2.43e-4-e)/0.559402)
    eps_tab[k]={'eps0':e,'10eps0':10*e,'F1_below_10eps0':9e-6<10*e,'gap_F2F5_over_10eps0_log10':float(np.log10(minbig/(10*e))),
                'eta_needed_for_F0R0':eta_need,'alpha_pairs_need_eta':any(x>e for x in alpha_res)}
out['eps_sens']=eps_tab
# (5) quantisation robustness: worst-case perturbation +-5e-7 per reading
gaps=sorted(abs(dD[a])-abs(dD[b]) for a,b in itertools.combinations(['F2','F3','F4','F5'],2))
out['min_gap_F2F5']=str(min(abs(g) for g in gaps)); out['worst_case_dD_error']=1e-6
# (6) identifiability: event-category design matrix for F2-F5
cats=['IBN 5100 去向','生理性别','人物生死','城市文化格局','时间旅行者行动/结局']
X=np.array([[1,0,0,0,0],[1,1,0,0,0],[1,0,1,1,0],[0,0,0,0,1]],float)
out['design_rank']=int(np.linalg.matrix_rank(X)); out['design_shape']=list(X.shape)
json.dump(out,open('results_v10.json','w'),ensure_ascii=False,indent=1,default=str)
# ---------- Fig 2
fig,(a,b)=plt.subplots(1,2,figsize=(183*MM,58*MM),gridspec_kw={'wspace':0.32})
t=tauvec(170/24)
for ax,X_,xl in [(a,t,'D-mail look-back time $\\tau$ (days)'),(b,np.where(Rvec(RR['F5'])==0,0.5,Rvec(RR['F5'])),'Constructed span $R$ (days; descriptor only)')]:
    ax.plot(X_,A,'o',ms=4,color=BLUE); ax.plot(X_[0],A[0],'o',ms=4,color=RED)
    for n,xx,yy,o in zip(N,X_,A,[(5,4),(5,-9),(4,5),(-12,6),(5,4)]): ax.annotate(n,(xx,yy),xytext=o,textcoords='offset points',fontsize=6)
    ax.set_xscale('log'); ax.set_yscale('log'); ax.set_ylim(3e-6,0.3); ax.set_xlabel(xl); ax.set_ylabel('$|\\Delta D|$ (%)')
    ax.axhline(4e-5,color=GREY,lw=0.6,ls=':')
a.text(0.35,0.62,'Spearman $\\rho=-0.30$\n(F1: 170 h, game text)',transform=a.transAxes,va='top',fontsize=6)
b.text(0.30,0.62,'$R\\approx\\tau$ for F2–F4 by causality;\ndiffers only at F1 ($R$=0) and F5 (†)',transform=b.transAxes,va='top',fontsize=6)
a.text(0.97,0.30,'$\\varepsilon_0$ scale (reversal/revisit),\nnot read-out error',transform=a.transAxes,ha='right',fontsize=5.5,color='#555')
b.annotate('F1: $R$=0, plotted at 0.5 d',(0.5,A[0]),xytext=(30,-4),textcoords='offset points',fontsize=5.5,color='#555',arrowprops=dict(arrowstyle='-',color=GREY,lw=0.5))
lab(a,'a',-0.2); lab(b,'b',-0.2); fig.savefig('fig2.png',bbox_inches='tight'); plt.close()
# ---------- Extended Data Fig 1: robustness
fig,ax=plt.subplots(1,4,figsize=(183*MM,50*MM),gridspec_kw={'wspace':0.55})
ts=np.linspace(2,7,201); ax[0].plot(ts,[rho(tauvec(x),A) for x in ts],color=BLUE,lw=1)
ax[0].axvline(4,color=GREY,ls=':',lw=0.6); ax[0].set_ylim(-1,1); ax[0].axhline(0.9,color=RED,lw=0.6,ls='--')
ax[0].set_xlabel('assumed $\\tau_{\\rm F1}$ (days)'); ax[0].set_ylabel('Spearman $\\rho(\\tau,|\\Delta D|)$')
ax[0].text(4.2,0.75,'5% extreme: $\\rho\\geq0.9$',fontsize=5.5,color=RED); lab(ax[0],'a',-0.32)
ax[1].plot(Rs/365.25,rsw,color=BLUE,lw=1); ax[1].set_xscale('log'); ax[1].axvspan(25,35,color=GOLD,alpha=0.2,lw=0)
ax[1].set_ylim(-1,1); ax[1].axhline(0.9,color=RED,lw=0.6,ls='--'); ax[1].set_xlabel('assumed $R_{\\rm F5}$ (years)'); ax[1].set_ylabel('$\\rho(R,|\\Delta D|)$')
ax[1].text(0.004,-0.55,'$R_{\\rm F5}=\\tau_{\\rm F5}$\n(no time-machine\nattribution)',fontsize=5.2); lab(ax[1],'b',-0.32)
xs=np.arange(5); w=0.27
for j,(k,c) in enumerate([('tau(F1=170h)',BLUE),('tau(F1=3d)','#7FA6D9'),('R',GOLD)]):
    ax[2].bar(xs+(j-1)*w,[loo[k][n] for n in N],w,color=c,label={'tau(F1=170h)':'$\\tau$, F1 after F2','tau(F1=3d)':'$\\tau$, F1 before F2','R':'$R$'}[k])
ax[2].set_xticks(xs); ax[2].set_xticklabels(['−'+n for n in N],fontsize=6); ax[2].axhline(0,color='k',lw=0.5)
ax[2].set_ylim(-1,1.05); ax[2].set_ylabel('$\\rho$ (n=4)'); ax[2].legend(fontsize=5,loc='lower left'); lab(ax[2],'c',-0.32)
ks=list(defs); e=[defs[k] for k in ks]
ax[3].scatter(range(len(ks)),[10*x for x in e],color=GREY,s=10,label='$10\\varepsilon_0$ (C1 threshold)')
ax[3].axhline(9e-6,color=RED,lw=0.7,label='F1 $|\\Delta D|$'); ax[3].axhline(minbig,color=BLUE,lw=0.7,label='min $|\\Delta D|$, F2–F5')
ax[3].set_yscale('log'); ax[3].set_ylim(3e-6,0.2); ax[3].set_xticks(range(len(ks))); ax[3].set_xticklabels(['α res','α rev','β ctx','4e-5','+F0/R0'],fontsize=5.5,rotation=30)
ax[3].set_ylabel('% '); ax[3].legend(fontsize=5,loc='upper left',bbox_to_anchor=(0,0.88)); ax[3].set_xlabel('$\\varepsilon_0$ definition'); lab(ax[3],'d',-0.32)
fig.savefig('edfig1_robustness.png',bbox_inches='tight'); plt.close()
print(json.dumps({k:out[k] for k in ['main','tau_sweep','R5_steps','R5_25_35','R5_threshold_years','loo','loo_P','min_gap_F2F5','design_rank','design_shape']},ensure_ascii=False,default=str))
print(json.dumps(out['eps_sens'],ensure_ascii=False))
print(json.dumps(out['perm5'])); print(json.dumps(out['perm4']))
print(json.dumps(out['dD'])); print(json.dumps(out['residuals'])); print(json.dumps(out['revisits']))

# ---------- v6: Kendall tau-a (no ties among |dD|) as alternative statistic
def kendall(x,y):
    c=d=0
    for i,j in itertools.combinations(range(len(x)),2):
        s=(x[i]-x[j])*(y[i]-y[j]); c+=s>0; d+=s<0
    return (c-d)/(len(x)*(len(x)-1)/2)
def kperm(n):
    b=list(range(n)); return [kendall(b,list(p)) for p in itertools.permutations(b)]
kout={}
for lbl,t1 in [('tau(F1=170h)',170/24),('tau(F1=3d)',3)]:
    t=list(tauvec(t1)); a=list(A)
    kout[lbl]={'all':round(kendall(t,a),3),
               'loo':{N[i]:round(kendall(t[:i]+t[i+1:],a[:i]+a[i+1:]),3) for i in range(5)}}
kout['R']=round(kendall(list(Rvec(RR['F5'])),list(A)),3)
k5,k4=kperm(5),kperm(4)
kout['P5']={str(v):f"{sum(z>=v-1e-3 for z in k5)}/120" for v in (kout['tau(F1=170h)']['all'],kout['tau(F1=3d)']['all'],kout['R'])}
kout['P4']={str(v):f"{sum(z>=v-1e-3 for z in k4)}/24" for v in (kout['tau(F1=170h)']['loo']['F5'],kout['tau(F1=3d)']['loo']['F5'])}
print(json.dumps({'kendall':kout},ensure_ascii=False))


# ---------- cross-source check: alternative values listed by the Japanese strategy wiki (probably novel versions); substitute one at a time
alts={'0.523299':'0.523291','0.456903':'0.456904','0.456914':'0.456919','1.130426':'1.130427','1.130205':'1.130204'}
alt_out={}
for old,new in alts.items():
    o2=[(k,new if a==old else a,new if b==old else b) for k,a,b in ops]; d2={k:Dm(b)-Dm(a) for k,a,b in o2}
    r2={f'{f}/{r}':str(d2[f]+d2[r]) for f,r in pairs[:4]}
    alt_out[f'{old}->{new}']={'alpha_res':r2,'max_abs_alpha_res':str(max(abs(d2[f]+d2[r]) for f,r in pairs[:4])),
        'min_F2F5':str(min(abs(d2[k]) for k in ('F2','F3','F4','F5'))),'F1':str(d2['F1'])}
print(json.dumps({'alt_values':alt_out}))

# ---------- exploratory cross-work analysis (NOT blind; not pooled with the main corpus)
ext=[('X1','1.129848','1.382733','third-party TT'),('X2','1.382733','1.129954','third-party undo'),
     ('X3','1.129848','1.064750','present'),('X4','1.129954','1.143688','present'),
     ('X5','1.053649','1.055821','present'),('X6','1.053649','1.064756','present'),
     ('X7','1.064750','1.081163','present'),('X8','1.097302','1.123581','D-RINE (anime)')]
xd={k:Dm(b)-Dm(a) for k,a,b,_ in ext}
x2_res=xd['X1']+xd['X2']; x2_rel=abs(x2_res)/abs(xd['X1'])
near={'1.064750/1.064756':str(Dm('1.064756')-Dm('1.064750')),'1.129848/1.129954':str(Dm('1.129954')-Dm('1.129848'))}
SG=Dm('1.048596'); toward={k:(abs(Dm(b)-SG)<abs(Dm(a)-SG)) for k,a,b,_ in ext}
print(json.dumps({'ext_dD':{k:str(v) for k,v in xd.items()},'X1X2_res':str(x2_res),'X1X2_rel':f'{x2_rel:.3e}',
  'abs_eps0_pass':float(abs(x2_res))<=4e-5,'rel_max_4.82e-4_pass':float(x2_rel)<=4.82e-4,
  'max_abs_res_passing_rel':f'{4.82e-4*float(abs(xd["X1"])):.2e}',
  'eta_needed_eps0=4e-5':f'{(float(abs(x2_res))-4e-5)/float(abs(xd["X1"])):.2e}','near_revisits':near,
  'signs':{k:('+' if v>0 else '-') for k,v in xd.items()},'toward_1.048596':toward,
  'branch_1.053649':str(xd['X6']-xd['X5'])}))
# Extended Data Fig 2
fig,(a,b)=plt.subplots(1,2,figsize=(183*MM,62*MM),gridspec_kw={'wspace':0.35})
groups=[('D-mail\n(SG VN, α)',[float(abs(dD[k])) for k in N],BLUE,N,'o'),
        ('D-RINE\n(SG0 anime, β)',[float(abs(xd['X8']))],BLUE,['X8'],'o'),
        ('present action\n(SG0 VN, β)',[float(abs(xd[k])) for k in ('X3','X4','X5','X6','X7')],RED,['X3','X4','X5','X6','X7'],'o'),
        ('3rd-party\ntime travel\n(SG0 VN, β)',[float(abs(xd['X1']))],GOLD,['X1'],'o')]
off={'F1':(5,2),'F2':(-14,-6),'F3':(-10,5),'F4':(4,-7),'F5':(5,3),'X3':(5,2),'X4':(-14,2),'X5':(5,-3),'X6':(3,-9),'X7':(5,2)}
for i,(g,v,c,labs,mk) in enumerate(groups):
    xs=np.full(len(v),i)+np.linspace(-0.18,0.18,len(v)) if len(v)>1 else [i]
    a.scatter(xs,v,s=14,color=c,zorder=3,marker=mk)
    for x,y,l in zip(xs,v,labs): a.annotate(l,(x,y),xytext=off.get(l,(4,3)),textcoords='offset points',fontsize=5)
a.axhline(4e-4,color=GREY,lw=0.6,ls=':'); a.text(3.45,4.6e-4,'10$\\varepsilon_0$ = 4×10$^{-4}$',ha='right',fontsize=5.5,color='#555')
a.set_yscale('log'); a.set_ylim(3e-6,1); a.set_xlim(-0.6,3.5); a.set_xticks(range(4)); a.set_xticklabels([g[0] for g in groups],fontsize=5.5)
a.set_ylabel('|$\\Delta D$| (%)'); lab(a,'a',-0.16)
pp=[('F5/R5',0.072233,1.1e-5,BLUE,'o'),('F4/R4',0.047483,2e-6,BLUE,'v'),('F3/R3',0.066396,3e-6,BLUE,'o'),('F2/R2',0.047716,2.3e-5,BLUE,'o'),
    ('F0/R0',0.559402,2.43e-4,GREY,'o'),('X1/X2',float(abs(xd['X1'])),float(abs(x2_res)),GOLD,'s')]
for n_,f,e,c,mk in pp:
    b.scatter(f,e,s=14,color=c,zorder=3,marker=mk); b.text(f*1.08,e*1.1,n_,fontsize=5)
xx=np.logspace(-2,0,50); b.plot(xx,4.82e-4*xx,color=RED,lw=0.7,ls='--')
b.text(0.13,1.6e-5,'max in-sample ratio\n$|\\varepsilon|/|\\Delta D|=4.8\\times10^{-4}$ (F2/R2)',fontsize=5.2,color=RED)
b.axhline(4e-5,color=GREY,lw=0.6,ls=':'); b.text(0.011,4.6e-5,'$\\varepsilon_0=4\\times10^{-5}$',fontsize=5.2,color='#555')
b.set_xscale('log'); b.set_yscale('log'); b.set_xlim(0.01,1); b.set_ylim(1e-6,1e-3)
b.set_xlabel('|$\\Delta D^{\\rm fwd}$| (%)'); b.set_ylabel('|$\\varepsilon$| (%)'); lab(b,'b',-0.2)
fig.savefig('edfig2_crosswork.png',bbox_inches='tight'); plt.close()