"""Statistics for 'World-line divergence in Steins;Gate' (exact tests, n=5)."""
import itertools, math
fwd={'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)}
rev={'F5':(0.337187,0.409431),'F4':(0.409431,0.456914),'F3':(0.456914,0.523307),
     'F2':(0.523307,0.571046),'F0':(0.571046,1.130205)}
dF={k:b-a for k,(a,b) in fwd.items()}
for k,(a,b) in rev.items():
    e=(b-a)+dF[k]; print(f'{k}: fwd {dF[k]:+.6f} rev {b-a:+.6f} eps {e:+.6f} rel {abs(e)/abs(dF[k]):.2e}')
def rho(x,y):
    n=len(x); return 1-6*sum((a-b)**2 for a,b in zip(x,y))/(n*(n*n-1))
dv=[1,3,4,2,5]                      # ranks of |dD| for F1..F5
cases={'tau (F1 after F2)':[3,2,5,4,1],'tau (F1 before F2)':[2,3,5,4,1],'R':[1,2,4,3,5]}
P=list(itertools.permutations(range(1,6)))
for name,r in cases.items():
    obs=rho(r,dv); p=sum(rho(q,dv)>=obs-1e-12 for q in P)/len(P)
    print(f'{name}: rho={obs:+.2f}, exact P(rho>=obs)={p:.4f}')
print('sign test, 5/5 negative: one-sided',0.5**5,'two-sided',2*0.5**5)
print('gap log10(0.047483/9e-6)=',round(math.log10(0.047483/9e-6),2))
# Sensitivity: exclude F1 (n=4; F2..F5). |dD| ranks: F2=2, F3=3, F4=1, F5=4
dv4=[2,3,1,4]; P4=list(itertools.permutations(range(1,5)))
for name,r in {'R, no F1':[1,3,2,4],'tau, no F1':[2,4,3,1]}.items():
    obs=rho(r,dv4); p=sum(rho(q,dv4)>=obs-1e-12 for q in P4)/len(P4)
    print(f'{name}: rho={obs:+.2f}, exact P(rho>=obs)={p:.4f}')
print('sign test excluding F1 (4/4): one-sided',0.5**4,'two-sided',2*0.5**4)
print('R0 alternative 1.130238 (Mayuri route): eps =',round((1.130238-0.571046)+(0.571024-1.130426),6))
import matplotlib as mpl, matplotlib.pyplot as plt, numpy as np
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.16,y=1.04): ax.text(x,y,s,transform=ax.transAxes,fontsize=8,fontweight='bold',va='bottom')
# ---------- Fig 1 trajectory
seq=[('Start',1.130426),('F0',0.571024),('F1',0.571015),('F2',0.523299),('F3',0.456903),('F4',0.409420),('F5',0.337187),
     ('R5',0.409431),('R4',0.456914),('R3',0.523307),('R2',0.571046),('R0',1.130205),('SG',1.048596)]
fig,(a,b)=plt.subplots(1,2,figsize=(183*MM,62*MM),gridspec_kw={'width_ratios':[1.25,1],'wspace':0.28})
x=np.arange(len(seq)); y=[v for _,v in seq]
a.axhspan(1,1.2,color=RED,alpha=0.07,lw=0); a.axhspan(0,1,color=BLUE,alpha=0.07,lw=0)
a.axhspan(-0.4,0,color=GOLD,alpha=0.08,lw=0)
a.step(x,y,where='post',color='k',lw=0.8); a.plot([12,12],[1.130205,1.048596],color='k',lw=0.8); a.plot(x,y,'o',ms=3,color='k',zorder=3)
a.plot([7,8],[0.409431,-0.275349],ls='--',color=GOLD,lw=0.8); a.plot(8,-0.275349,'s',ms=3,color=GOLD)
a.text(8.2,-0.275349,'Faris route\n−0.275349 (Ω)',va='center',fontsize=6,color='#8a6500')
a.plot(12,1.048596,'*',ms=7,color='k',mfc='white',zorder=4)
a.set_xticks(x); a.set_xticklabels([s for s,_ in seq],fontsize=6)
a.set_ylabel('Divergence $D$ (%)'); a.set_ylim(-0.4,1.2); a.set_xlim(-0.5,12.7)
for t,xx,yy,c in [('β field',1.6,1.19,RED),('α field',-0.35,0.97,BLUE),('Ω field',-0.35,-0.03,'#8a6500')]: a.text(xx,yy,t,color=c,fontsize=6.5,va='top')
a.text(11.75,0.93,'Steins Gate 1.048596\n(listed separately,\nnot β; ref. 4)',ha='center',va='top',fontsize=5.5,color='k')
a.set_xlabel('Transition (chronological)'); lab(a,'a',-0.1)
# b zoom alpha with forward/back
fw=[0.571024,0.571015,0.523299,0.456903,0.409420,0.337187]; fl=['F0','F1','F2','F3','F4','F5']
bw=[0.337187,0.409431,0.456914,0.523307,0.571046]; bl=['F5','R5','R4','R3','R2']
b.plot(range(6),fw,'-o',color=BLUE,ms=3,lw=0.8,label='Sending D-mails (F1–F5)')
b.plot([5,4,3,2,1.0],bw,'-s',color=RED,ms=3,lw=0.8,mfc='white',label='Cancelling D-mails (R5–R2)')
b.set_xticks(range(6)); b.set_xticklabels(['after F0','F1','F2','F3','F4','F5'],fontsize=6)
b.set_ylabel('Divergence $D$ (%)'); b.set_xlabel('Set of active D-mails, $\\mathcal{A}$')
b.legend(loc='lower left',fontsize=6); lab(b,'b',-0.2)
ax2=b.inset_axes([0.62,0.58,0.35,0.36]); ax2.plot([0],[0.571024],'o',color=BLUE,ms=3); ax2.plot([0],[0.571015],'o',color=BLUE,ms=3,mfc='none')
ax2.plot([0],[0.571046],'s',color=RED,ms=3,mfc='white'); ax2.set_xticks([]); ax2.set_ylim(0.57100,0.57106)
ax2.ticklabel_format(useOffset=False); ax2.tick_params(labelsize=5); ax2.set_title('revisit offset ≤3.1×10$^{-5}$',fontsize=5.5)
fig.savefig('fig1.png',bbox_inches='tight'); plt.close()
# ---------- Fig 2 tau vs |dD|, R vs |dD|
names=['F1','F2','F3','F4','F5']; dD=np.array([9e-6,0.047716,0.066396,0.047483,0.072233])
tau=np.array([4,4,17.5*365,10*365,1]); tlo=np.array([2,4,17*365,10*365,1]); thi=np.array([7,4,18*365,10*365,1])
R=np.array([0.5/365,4/365,17,10,35]); Rlo=R.copy(); Rhi=R.copy()
fig,(a,b)=plt.subplots(1,2,figsize=(183*MM,58*MM),gridspec_kw={'wspace':0.32})
for ax,X,lo,hi,xl in [(a,tau,tlo,thi,'Look-back time $\\tau$ (days)'),(b,R,Rlo,Rhi,'Causal span of rewritten macro-history $R$ (years)')]:
  ax.errorbar(X,dD,xerr=[X-lo,hi-X],fmt='o',ms=4,color=BLUE,ecolor=GREY,elinewidth=0.6,capsize=1.5)
  ax.plot(X[0],dD[0],'o',ms=4,color=RED)
  for n,xx,yy,o in zip(names,X,dD,[(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.axhspan(3e-6,4e-5,color=GREY,alpha=0.12,lw=0)
a.text(0.25,0.55,'Spearman $\\rho=-0.30$ (or $-0.10$)\nexact $P(\\rho\\geq\\rho_{\\rm obs})=0.74$ (0.61)',transform=a.transAxes,va='top',fontsize=6)
b.text(0.25,0.55,'Spearman $\\rho=0.90$\nexact one-sided $P=0.042$',transform=b.transAxes,va='top',fontsize=6)
a.text(0.97,0.06,'$|\\Delta D|\\leq\\varepsilon_0=4\\times10^{-5}$',transform=a.transAxes,ha='right',fontsize=5.5,color='#555')
b.annotate('F1: no macro change\n(plotted at nominal $R$)',(R[0],dD[0]),xytext=(25,18),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()
# ---------- Fig 3 reversibility
pairs=['F0/R0','F2/R2','F3/R3','F4/R4','F5/R5']; f=np.array([0.559402,0.047716,0.066396,0.047483,0.072233]); r=np.array([0.559159,0.047739,0.066393,0.047483,0.072244])
eps=r-f; rel=np.abs(eps)/f*100
fig,(a,b)=plt.subplots(1,2,figsize=(183*MM,58*MM),gridspec_kw={'wspace':0.35})
g=np.logspace(-2,0,10); a.plot(g,g,color=GREY,lw=0.6,ls='--'); a.plot(f,r,'o',ms=4,color=BLUE)
for n,xx,yy,o in zip(['F0/R0','F2/R2 & F4/R4','F3/R3','',' F5/R5'],f,r,[(-10,-11),(6,-10),(-5,3),(6,-9),(4,-2)]): a.annotate(n,(xx,yy),xytext=o,textcoords='offset points',fontsize=6,ha='right' if n=='F3/R3' else 'left')
a.set_xscale('log'); a.set_yscale('log'); a.set_xlabel('$|\\Delta D_{\\rm fwd}|$ (%)'); a.set_ylabel('$|\\Delta D_{\\rm rev}|$ (%)'); a.text(0.04,0.9,'$y=x$',transform=a.transAxes,color=GREY,fontsize=6)
b.bar(range(5),rel,color=[RED if e<0 else BLUE for e in eps],width=0.6)
for i,(e,v) in enumerate(zip(eps,rel)): b.text(i,v+0.0015,('0' if abs(e)<1e-9 else (('+' if e>0 else '−')+f'{abs(e)*1e6:.0f}×10$^{{-6}}$')),ha='center',fontsize=5.5)
b.axhline(0.0482,color=GREY,lw=0.6,ls=':'); b.text(4.4,0.0495,'$\\hat\\eta=4.8\\times10^{-4}$',ha='right',fontsize=6,color='#555')
b.set_xticks(range(5)); b.set_xticklabels(pairs,fontsize=6); b.set_ylabel('Relative residual $|\\varepsilon|/|\\Delta D_{\\rm fwd}|$ (%)'); b.set_ylim(0,0.06)
lab(a,'a',-0.2); lab(b,'b',-0.2)
fig.savefig('fig3.png',bbox_inches='tight'); plt.close()
# ---------- Fig 4 schematic
fig,(a,b)=plt.subplots(1,2,figsize=(183*MM,55*MM),gridspec_kw={'wspace':0.25})
xs=np.linspace(-2.2,2.2,400); V=(xs**2-1)**2*1.0+0.15*np.cos(9*xs)*0.25
a.plot(xs,V,color='k',lw=0.9); a.axvline(0,color=GREY,lw=0.5,ls=':')
a.text(-1,-0.32,'α basin\nf = α',ha='center',color=BLUE,fontsize=6.5); a.text(1,-0.32,'β basin\nf = β',ha='center',color=RED,fontsize=6.5)
a.annotate('',xy=(-1.18,0.12),xytext=(-0.88,0.12),arrowprops=dict(arrowstyle='->',color=BLUE,lw=0.8)); a.text(-2.05,0.35,'within-field\nshift $\\delta$',fontsize=6,color=BLUE)
a.annotate('',xy=(0.95,0.1),xytext=(-0.8,0.1),arrowprops=dict(arrowstyle='->',color=RED,lw=0.8,connectionstyle='arc3,rad=-0.45'))
a.text(0,1.2,'flip of convergence-\ndetermining variable $c$',ha='center',fontsize=6,color=RED)
a.set_ylim(-0.45,1.6); a.set_xticks([]); a.set_yticks([]); a.set_xlabel('Macroscopic history state'); a.set_ylabel('Effective potential (schematic)'); lab(a,'a',-0.08)
t=np.linspace(0,1,200); 
for k,(amp,col,ls) in enumerate([(1,BLUE,'-'),(1,GREY,'--')]): pass
b.plot(t,np.exp(-8*t),color=GREY,lw=0.9,label='micro component $P_S\\,\\delta x$ (absorbed)')
b.plot(t,0.25+0.75*(1-np.exp(-6*t)),color=BLUE,lw=0.9,label='macro component $P_M\\,\\delta x$ (amplified)')
b.text(0.5,1.02,'conceptual schematic — not fitted to data',transform=b.transAxes,ha='center',fontsize=6,color='#555'); a.text(0.5,1.02,'conceptual schematic — not fitted to data',transform=a.transAxes,ha='center',fontsize=6,color='#555')
b.set_xlabel('Time since injection, $(t-t_0)/(T-t_0)$'); b.set_ylabel('Perturbation norm (a.u.)'); b.set_yticks([]); b.legend(fontsize=6,loc='center right'); lab(b,'b',-0.08)
fig.savefig('fig4.png',bbox_inches='tight'); plt.close()