"""
QDI Manuscript — All Figures
Generated from real IBM hardware data.
Author: Ramakrishna Pasupuleti, July 2026
"""
import json, numpy as np, matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.patches import FancyArrowPatch, Rectangle
from matplotlib.lines import Line2D
from scipy.optimize import curve_fit
from scipy.stats import linregress, t as t_dist

with open('/home/claude/fig_data.json') as f: D=json.load(f)

plt.rcParams.update({
    'font.family':'serif','font.serif':['Times New Roman','Times','DejaVu Serif'],
    'font.size':11,'axes.labelsize':11,'axes.titlesize':12,
    'xtick.labelsize':10,'ytick.labelsize':10,'legend.fontsize':9,
    'lines.linewidth':1.8,'lines.markersize':7,
    'figure.dpi':300,'savefig.dpi':300,'savefig.bbox':'tight',
    'axes.spines.top':False,'axes.spines.right':False,
})

CK='#2166ac'; CF='#1a9641'; CM='#d73027'; CG='#762a83'
LGRAY='#999999'

def kr_model(K,C0,a,R): return C0/K**a+R
def kr_curve(x_range, popt):
    return kr_model(1.0/np.array(x_range), *popt)

# ── FIGURE 1: QDI ARCHITECTURE (conceptual) ──────────────
print("Fig 1: QDI Architecture...")
fig,ax=plt.subplots(figsize=(10,5))
ax.set_xlim(0,10); ax.set_ylim(0,5); ax.axis('off')

boxes=[
    (0.3,2.0,1.4,1.2,'INPUT\n13 circuits\n3 minutes','#F0F8FF'),
    (2.1,2.0,1.8,1.2,'LAYER 1\nBiomarker\nextraction\n\u03c7=(α,CV,R²)','#E8F5E9'),
    (4.2,2.0,1.8,1.2,'LAYER 2\nDiagnostic\nstate D(\u03c7)','#FFF3E0'),
    (6.3,2.0,1.8,1.2,'LAYER 3\nHealth\nscore H(\u03c7)','#FCE4EC'),
    (8.3,3.0,1.4,1.0,'LAYER 4\nRecommend.\nR(\u03c7)','#EDE7F6'),
    (8.3,1.5,1.4,1.0,'LAYER 5\nForecast\nF(\u03c7,t)','#E1F5FE'),
]
for x,y,w,h,txt,color in boxes:
    rect=Rectangle((x,y),w,h,linewidth=1.5,edgecolor='#333333',
                   facecolor=color,zorder=3,alpha=0.95)
    ax.add_patch(rect)
    ax.text(x+w/2,y+h/2,txt,ha='center',va='center',fontsize=8.5,
            fontweight='bold' if 'INPUT' in txt else 'normal',zorder=4,
            multialignment='center')

# Arrows
arrow_kw=dict(arrowstyle='-|>',color='#333333',lw=1.5,mutation_scale=14)
ax.annotate('',xy=(2.1,2.6),xytext=(1.7,2.6),arrowprops=arrow_kw)
ax.annotate('',xy=(4.2,2.6),xytext=(3.9,2.6),arrowprops=arrow_kw)
ax.annotate('',xy=(6.3,2.6),xytext=(6.0,2.6),arrowprops=arrow_kw)
ax.annotate('',xy=(8.3,3.3),xytext=(8.1,3.3),arrowprops=dict(arrowstyle='-|>',color='#333333',lw=1.2,mutation_scale=12))
ax.annotate('',xy=(8.3,2.0),xytext=(8.1,2.3),arrowprops=dict(arrowstyle='-|>',color='#333333',lw=1.2,mutation_scale=12))
ax.plot([8.1,8.1],[2.3,3.3],color='#333333',lw=1.2)
ax.annotate('',xy=(8.1,2.8),xytext=(8.0,2.8),)
ax.plot([8.0,8.1],[2.6,2.6],color='#333333',lw=1.2)

# Labels below boxes
labels=['','GHZ, Depth,\nEcho, Amp\ncircuits','5 states\n(physics rules)','0\u2013100\nexplainable','Do/Avoid\nMitigate','α(t)±CI\nLOW conf.']
for i,(x,y,w,h,txt,c) in enumerate(boxes):
    if labels[i]:
        ax.text(x+w/2,y-0.25,labels[i],ha='center',va='top',fontsize=7.5,
                color='#555555',multialignment='center')

ax.set_title('Fig. 1  Quantum Diagnostic Intelligence (QDI) — Five-Layer Architecture\n'
             'χ=(α_GHZ, α_depth, α_echo, α_amp, CV_drift, CV_shot, R²_GHZ, R²_depth)',
             fontsize=11,fontweight='bold',pad=10)
plt.tight_layout()
plt.savefig('/home/claude/Fig1_QDI_Architecture.pdf',format='pdf')
plt.savefig('/home/claude/Fig1_QDI_Architecture.png')
plt.close()
print("  Saved.")

# ── FIGURE 2: CENTRAL FINDING — DIVERGENCE ────────────────
print("Fig 2: Central Finding (Divergence)...")
fig=plt.figure(figsize=(12,5))
gs=gridspec.GridSpec(1,2,figure=fig,wspace=0.32)
ax1=fig.add_subplot(gs[0]); ax2=fig.add_subplot(gs[1])

# (a) Timeline: IBM T2 vs Health Score
sessions=['Jul 3\n(C4)','Jul 5\nRun 1','Jul 5\nRun 2','Jul 5\nRun 3']
T2_vals=[303.1,339.6,339.6,339.6]
health_vals=[82.6,37.8,32.1,15.4]  # computed from actual alphas
idx=[0,1,2,3]

ax1b=ax1.twinx()
l1,=ax1.plot(idx,T2_vals,'o-',color='steelblue',lw=2,ms=9,label='IBM T₂ Q2 (μs)',zorder=5)
ax1.fill_between(idx,[t-5 for t in T2_vals],[t+5 for t in T2_vals],alpha=0.1,color='steelblue')
l2,=ax1b.plot(idx,health_vals,'s--',color=CM,lw=2,ms=9,label='QDI Health Score',zorder=5)
ax1b.fill_between(idx,[0,7.6,0,0],[100,67.9,100,100],alpha=0.06,color=CM)

ax1.axhline(y=250,color='steelblue',ls=':',alpha=0.4)
ax1b.axhline(y=55,color=CM,ls=':',alpha=0.4,label='Warning threshold (55)')

ax1.annotate('T₂: +12%\n(reported improved)',xy=(1,339.6),xytext=(1.5,295),
    fontsize=9,color='steelblue',ha='center',
    arrowprops=dict(arrowstyle='->',color='steelblue',lw=1))
ax1b.annotate('Score: −45 pts\n(CRITICAL)',xy=(1,37.8),xytext=(1.5,62),
    fontsize=9,color=CM,ha='center',
    arrowprops=dict(arrowstyle='->',color=CM,lw=1))

ax1.set_xticks(idx); ax1.set_xticklabels(sessions,fontsize=9)
ax1.set_ylabel('IBM T₂ (μs)',color='steelblue',fontsize=10)
ax1b.set_ylabel('QDI Health Score',color=CM,fontsize=10)
ax1.tick_params(axis='y',labelcolor='steelblue')
ax1b.tick_params(axis='y',labelcolor=CM)
ax1.set_ylim(200,420); ax1b.set_ylim(-5,110)
ax1b.axhspan(0,30,alpha=0.06,color='red')
ax1b.text(3.9,15,'CRITICAL',fontsize=7,color='red',ha='right',style='italic')
ax1b.axhspan(30,55,alpha=0.04,color='orange')
lines=[l1,l2]; labs=[l.get_label() for l in lines]
ax1.legend(lines,labs,loc='lower left',fontsize=8)
ax1.set_title('(a) IBM T₂ vs QDI Health Score\n(ibm_marrakesh, Jul 2026)',
              fontsize=10,fontweight='bold')
ax1.grid(True,alpha=0.25)

# (b) GHZ infidelity Jul3 vs Jul5
m_ghz_n=D['m_ghz_n']; m_ghz_inf=D['m_ghz_inf']
m5_ghz_n=D['m5_ghz_n']; m5_ghz_inf=D['m5_ghz_inf']
mpopt=D['m_ghz_popt']
m5popt=D['m5_ghz_popt']
n_sm=np.linspace(2,30,300)

ax2.plot(m_ghz_n,m_ghz_inf,'o',color=CK,ms=8,zorder=5,
         label=f'Jul 3 (α={D["m_ghz_alpha"]:.3f}, R²={D["m_ghz_R2"]:.3f})')
ax2.plot(n_sm,kr_curve(n_sm,mpopt),'--',color=CK,alpha=0.7,lw=1.5)

n5=np.linspace(2,20,200)
ax2.plot(m5_ghz_n,m5_ghz_inf,'s',color=CM,ms=8,zorder=5,
         label=f'Jul 5 (α={D["m5_ghz_alpha"]:.3f}, R²={D["m5_ghz_R2"]:.4f})')
ax2.plot(n5,kr_curve(n5,m5popt),'--',color=CM,alpha=0.7,lw=1.5)

# Annotate n=2 jump
ax2.annotate('',xy=(2,0.204),xytext=(2,0.020),
    arrowprops=dict(arrowstyle='<->',color='gray',lw=1.5))
ax2.text(2.4,0.11,'+920%\nat n=2',fontsize=9,color='gray',
         ha='left',va='center',fontweight='bold')

ax2.set_xlabel('Number of Qubits n',fontsize=10)
ax2.set_ylabel('GHZ Infidelity',fontsize=10)
ax2.set_title('(b) GHZ Scaling: Jul 3 vs Jul 5\n(ibm_marrakesh)',
              fontsize=10,fontweight='bold')
ax2.legend(fontsize=8,loc='upper left'); ax2.grid(True,alpha=0.25)
ax2.set_xlim(0,32); ax2.set_ylim(-0.02,1.0)
ax2.text(0.97,0.03,'IBM T₂ improved +12%\nbetween these sessions',
    transform=ax2.transAxes,ha='right',va='bottom',fontsize=8,
    color='steelblue',style='italic',
    bbox=dict(boxstyle='round',facecolor='lightblue',alpha=0.3))

plt.suptitle('Fig. 2  Central Finding: Pulse-Level vs Circuit-Level Divergence\n'
             'ibm_marrakesh (156-qubit Heron r2)',
             fontsize=11,fontweight='bold',y=1.02)
plt.savefig('/home/claude/Fig2_Divergence.pdf',format='pdf')
plt.savefig('/home/claude/Fig2_Divergence.png')
plt.close()
print("  Saved.")

# ── FIGURE 3: THREE-DEVICE FINGERPRINT ───────────────────
print("Fig 3: Three-Device Fingerprint...")
fig,axes=plt.subplots(1,3,figsize=(13,5))

# (a) GHZ scaling
ax=axes[0]
kpopt=D['k_ghz_popt']; fpopt=D['f_ghz_popt']
k_n=D['k_ghz_n']; k_inf=D['k_ghz_inf']
f_n=D['f_ghz_n']; f_inf=D['f_ghz_inf']
n10=np.linspace(2,10,100)
n30=np.linspace(2,30,300)
ax.plot(k_n,k_inf,'o',color=CK,ms=8,label=f'ibm_kingston (α={D["k_ghz_alpha"]:.3f})')
ax.plot(n10,kr_curve(n10,kpopt),'--',color=CK,alpha=0.7)
ax.plot(f_n,f_inf,'s',color=CF,ms=8,label=f'ibm_fez (α={D["f_ghz_alpha"]:.3f})')
ax.plot(n10,kr_curve(n10,fpopt),'--',color=CF,alpha=0.7)
ax.plot(m_ghz_n,m_ghz_inf,'^',color=CM,ms=8,label=f'ibm_marrakesh Jul3 (α={D["m_ghz_alpha"]:.3f})')
ax.plot(n30,kr_curve(n30,mpopt),'--',color=CM,alpha=0.7)
ax.set_xlabel('Number of Qubits n',fontsize=10)
ax.set_ylabel('GHZ Infidelity',fontsize=10)
ax.set_title('(a) GHZ Qubit Scaling',fontsize=10,fontweight='bold')
ax.legend(fontsize=7.5); ax.grid(True,alpha=0.25)
ax.annotate('n=30',xy=(30,0.844),xytext=(26,0.65),fontsize=8,color=CM,
    arrowprops=dict(arrowstyle='->',color=CM,lw=1))

# (b) Depth scaling
ax=axes[1]
k_d=D['k_dep_d']; k_di=D['k_dep_inf']
f_d=D['f_dep_d']; f_di=D['f_dep_inf']
m_d=D['m_dep_d']; m_di=D['m_dep_inf']
d_sm=np.linspace(1,64,200)
from scipy.optimize import curve_fit as cf2
def fit_for_plot(x,y):
    K=1/np.array(x,float)
    try:
        popt,_=cf2(kr_model,K,np.array(y,float),
            p0=[0.1,0.5,0.01],bounds=([0,0.001,0],[100,15,1]),method='trf',maxfev=50000)
        return popt
    except: return None
kpp=fit_for_plot(k_d,k_di)
fpp=fit_for_plot(f_d,f_di)
mpp=fit_for_plot(m_d,m_di)
ax.plot(k_d,k_di,'o',color=CK,ms=8,label=f'ibm_kingston (α={D["k_dep_alpha"]:.3f})')
if kpp is not None: ax.plot(d_sm,kr_curve(d_sm,kpp),'--',color=CK,alpha=0.7)
ax.plot(f_d,f_di,'s',color=CF,ms=8,label=f'ibm_fez (α={D["f_dep_alpha"]:.3f})')
if fpp is not None: ax.plot(d_sm,kr_curve(d_sm,fpp),'--',color=CF,alpha=0.7)
ax.plot(m_d,m_di,'^',color=CM,ms=8,label=f'ibm_marrakesh Jul3 (α={D["m_dep_alpha"]:.3f})')
if mpp is not None: ax.plot(d_sm,kr_curve(d_sm,mpp),'--',color=CM,alpha=0.7)
ax.set_xlabel('Circuit Depth d',fontsize=10)
ax.set_ylabel('Clifford Infidelity',fontsize=10)
ax.set_title('(b) Clifford Depth Scaling',fontsize=10,fontweight='bold')
ax.legend(fontsize=7.5); ax.grid(True,alpha=0.25)

# (c) Fingerprint radar/bar
ax=axes[2]
devs=['kingston','fez','mar_Jul3','mar_Jul5']
labels_bar=['α_GHZ','α_depth','α_echo','α_amp','CV_drift/5']
vals={
    'kingston':[D['k_ghz_alpha'],D['k_dep_alpha'],D['k_echo_alpha'],D['k_amp_alpha'],D['k_stab_CV']/5],
    'fez':     [D['f_ghz_alpha'],D['f_dep_alpha'],D['f_echo_alpha'],0,D['f_stab_CV']/5],
    'mar_Jul3':[D['m_ghz_alpha'],D['m_dep_alpha'],D['m5_Echo_alpha'],D['m_amp_alpha'],D['m_drift_CV']/5],
    'mar_Jul5':[D['m5_ghz_alpha'],D['m5_Cliff_alpha'],D['m5_Echo_alpha'],D['m5_Amp_alpha'],D['m5_drift_CV']/5],
}
x=np.arange(5); w=0.18
colors_bar=[CK,CF,CM,'#8B0000']
for i,(dev,color) in enumerate(zip(devs,colors_bar)):
    v=vals[dev]
    ax.bar(x+i*w,v,w,color=color,alpha=0.85,edgecolor='black',lw=0.5,
           label=dev.replace('_',' '))

ax.axhline(y=1.0,color='red',ls='--',lw=1.2,alpha=0.7,label='α=1 boundary')
ax.axhline(y=2.0,color='darkred',ls=':',lw=1.0,alpha=0.5,label='α=2 boundary')
ax.set_xticks(x+1.5*w)
ax.set_xticklabels(labels_bar,fontsize=8.5)
ax.set_ylabel('Biomarker Value',fontsize=10)
ax.set_title('(c) QDI Fingerprint\n(4 biomarkers + CV/5)',fontsize=10,fontweight='bold')
ax.legend(fontsize=7,ncol=2); ax.grid(True,alpha=0.25,axis='y')
ax.set_ylim(0,4.2)

plt.suptitle('Fig. 3  Three-Device K–R Fingerprint: ibm_kingston, ibm_fez, ibm_marrakesh',
             fontsize=11,fontweight='bold',y=1.02)
plt.tight_layout()
plt.savefig('/home/claude/Fig3_ThreeDevice.pdf',format='pdf')
plt.savefig('/home/claude/Fig3_ThreeDevice.png')
plt.close()
print("  Saved.")

# ── FIGURE 4: HEALTH SCORE BREAKDOWN ─────────────────────
print("Fig 4: Health Score Breakdown...")
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(11,5))

# (a) Health scores comparison
devices=['ibm_kingston','ibm_fez','ibm_marrakesh\nJul 3','ibm_marrakesh\nJul 5']
scores=[94.8,95.1,82.6,37.8]
ci_lo=[80.8,91.1,32.3,7.6]
ci_hi=[100,99.0,100,67.9]
colors_h=[CK,CF,CM,'#8B0000']

bars=ax1.bar(range(4),scores,color=colors_h,alpha=0.85,edgecolor='black',lw=0.8,width=0.5)
for i,(s,lo,hi) in enumerate(zip(scores,ci_lo,ci_hi)):
    ax1.errorbar(i,s,yerr=[[s-lo],[hi-s]],fmt='none',color='black',capsize=6,lw=1.5,capthick=1.5)
    ax1.text(i,hi+2,f'{s:.1f}',ha='center',va='bottom',fontsize=9,fontweight='bold')

ax1.axhspan(80,100,alpha=0.05,color='green'); ax1.text(3.7,87,'HEALTHY',fontsize=7.5,color='green',style='italic')
ax1.axhspan(55,80,alpha=0.05,color='orange'); ax1.text(3.7,65,'WARNING',fontsize=7.5,color='orange',style='italic')
ax1.axhspan(0,55,alpha=0.05,color='red');    ax1.text(3.7,30,'CRITICAL',fontsize=7.5,color='red',style='italic')
ax1.axhline(y=55,color='orange',ls='--',lw=1,alpha=0.6)
ax1.axhline(y=80,color='green',ls='--',lw=1,alpha=0.6)

# T2 annotation
T2s=[147.1,303.1,303.1,339.6]
for i,t2 in enumerate(T2s):
    ax1.text(i,-8,f'T₂={t2:.0f}μs',ha='center',va='top',fontsize=8,color='steelblue')

ax1.set_xticks(range(4)); ax1.set_xticklabels(devices,fontsize=8.5)
ax1.set_ylabel('QDI Health Score H(χ)',fontsize=10)
ax1.set_title('(a) Health Scores with 95% CI\n(Error bars = uncertainty from biomarker σ)',
              fontsize=10,fontweight='bold')
ax1.set_ylim(-15,112); ax1.grid(True,alpha=0.25,axis='y')
ax1.text(0,-15,'IBM T₂ (μs):',ha='right',va='top',fontsize=8,color='steelblue')

# (b) Breakdown for marrakesh Jul5
components=['α_GHZ\n(3.811>1.5)','CV_drift\n(10.2%>3%)','R² quality\n(satisfactory)','Proof4\n(dev=0.14)']
penalties=[40.0,21.5,0.0,0.7]
pcts=[64,35,0,1]
colorsb=['#d73027','#fc8d59','#91bfdb','#e0f3f8']

wedges,texts,autotexts=ax2.pie(
    [40.0,21.5,0.7,0.001],
    labels=None,
    colors=['#d73027','#fc8d59','#fdae61','#91bfdb'],
    autopct=lambda p: f'{p:.0f}%' if p>2 else '',
    startangle=90, pctdistance=0.75,
    wedgeprops={'edgecolor':'white','linewidth':2}
)
for at in autotexts: at.set_fontsize(10); at.set_fontweight('bold')

legend_labels=[f'α_GHZ penalty: −40.0 pts (64%)',
               f'CV_drift penalty: −21.5 pts (35%)',
               f'Proof4 penalty: −0.7 pts (1%)',
               f'R² penalty: 0.0 pts (0%)']
ax2.legend(legend_labels,loc='lower center',bbox_to_anchor=(0.5,-0.2),
           fontsize=8,ncol=2)
ax2.set_title('(b) Score Breakdown: ibm_marrakesh Jul 5\n'
              'H(χ) = 37.8/100  [95% CI: 7.6, 67.9]',
              fontsize=10,fontweight='bold')
ax2.text(0,0,'Score\n37.8/100\nCRITICAL',ha='center',va='center',fontsize=10,
         fontweight='bold',color='#d73027')

plt.suptitle('Fig. 4  QDI Health Score H(χ): Results and Explainable Decomposition',
             fontsize=11,fontweight='bold',y=1.02)
plt.tight_layout()
plt.savefig('/home/claude/Fig4_HealthScore.pdf',format='pdf')
plt.savefig('/home/claude/Fig4_HealthScore.png')
plt.close()
print("  Saved.")

# ── FIGURE 5: PHASE DIAGRAM + DRIFT FORECAST ─────────────
print("Fig 5: Phase Diagram and Drift Forecast...")
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(12,5.5))

# (a) (α, R²) Phase Diagram
# Zones
ax1.axhspan(0.8,1.06,xmin=0,xmax=1/4.5,color='#d1e8ff',alpha=0.4,zorder=0)
ax1.axhspan(0.8,1.06,xmin=1/4.5,xmax=1.0,color='#d4edda',alpha=0.4,zorder=0)
ax1.axhspan(-0.05,0.5,color='#fde0dc',alpha=0.3,zorder=0)
ax1.axhspan(0.5,0.8,color='#fff9c4',alpha=0.3,zorder=0)
ax1.text(0.25,0.95,'Stochastic',style='italic',fontsize=10,color='#1565c0',ha='center',transform=ax1.transAxes)
ax1.text(0.72,0.95,'Structured',style='italic',fontsize=10,color='#1b5e20',ha='center',transform=ax1.transAxes)
ax1.text(0.10,0.18,'Model\nBreakdown',style='italic',fontsize=9,color='#b71c1c',ha='center',transform=ax1.transAxes)
ax1.axhline(y=0.8,color='gray',ls='--',lw=0.8,alpha=0.6)
ax1.axhline(y=0.5,color='gray',ls='--',lw=0.8,alpha=0.6)
ax1.axvline(x=1.0,color='gray',ls='--',lw=0.8,alpha=0.6)
ax1.axvline(x=2.0,color='darkred',ls=':',lw=0.8,alpha=0.5)
ax1.text(2.02,0.52,'α=2.0',fontsize=8,color='darkred',alpha=0.7)

# All hardware points
hw_pts=[
    (D['k_ghz_alpha'],D['k_ghz_R2'],CK,'o','ibm_kingston'),
    (D['k_dep_alpha'],D['k_dep_R2'],CK,'*',''),
    (D['k_echo_alpha'],0.894,CK,'D',''),
    (D['f_ghz_alpha'],D['f_ghz_R2'],CF,'o','ibm_fez'),
    (D['f_dep_alpha'],D['f_dep_R2'],CF,'*',''),
    (D['m_ghz_alpha'],D['m_ghz_R2'],CM,'^','ibm_marrakesh Jul3'),
    (D['m_dep_alpha'],D['m_dep_R2'],CM,'v',''),
    (D['m5_ghz_alpha'],D['m5_ghz_R2'],'#8B0000','p','ibm_marrakesh Jul5'),
    (D['m5_Cliff_alpha'],D['m5_Cliff_R2'],'#8B0000','H',''),
    (D['m5_Echo_alpha'],D['m5_Echo_R2'],'#8B0000','X',''),
]
for a,r2,c,mk,lbl in hw_pts:
    ax1.scatter(a,r2,marker=mk,s=100,color=c,zorder=6,
                edgecolors='black' if mk=='p' else c,linewidths=1.0,
                label=lbl if lbl else None)

# Arrow Jul3→Jul5
ax1.annotate('',xy=(D['m5_ghz_alpha'],D['m5_ghz_R2']),
    xytext=(D['m_ghz_alpha'],D['m_ghz_R2']),
    arrowprops=dict(arrowstyle='-|>',color='darkred',lw=2.0,mutation_scale=16))
ax1.text(2.5,0.88,'Jul 3→Jul 5\n+171% α',fontsize=8.5,color='darkred',
         ha='center',fontweight='bold')

legend_elements=[
    Line2D([0],[0],marker='o',color='w',mfc=CK,ms=10,mec=CK,label='ibm_kingston'),
    Line2D([0],[0],marker='o',color='w',mfc=CF,ms=10,mec=CF,label='ibm_fez'),
    Line2D([0],[0],marker='^',color='w',mfc=CM,ms=10,mec=CM,label='ibm_marrakesh Jul3'),
    Line2D([0],[0],marker='p',color='w',mfc='#8B0000',ms=12,mec='black',label='ibm_marrakesh Jul5'),
]
ax1.legend(handles=legend_elements,loc='lower right',fontsize=8)
ax1.set_xlabel('K–R Scaling Exponent α',fontsize=10)
ax1.set_ylabel('Goodness-of-Fit R²',fontsize=10)
ax1.set_title('(a) (α, R²) Phase Diagram\nAll real hardware measurements',
              fontsize=10,fontweight='bold')
ax1.set_xlim(-0.1,4.5); ax1.set_ylim(-0.08,1.07)
ax1.grid(True,alpha=0.25)

# (b) Drift Forecast with CI
t_obs=np.array(D['jul5_times']); a_obs=np.array(D['jul5_alphas'])
sl=D['jul5_slope']; ic=D['jul5_intercept']; se=D['jul5_slope_se']
n=3; t_mean=np.mean(t_obs); S_xx=np.sum((t_obs-t_mean)**2)
s2_res=(1-D['jul5_R2'])*np.var(a_obs,ddof=1)*(n-1)/(n-2)
t_crit=t_dist.ppf(0.975,1)  # df=1

t_all=np.linspace(0,5,300)
mu_all=np.minimum(ic+sl*t_all,5.0)
PI_all=np.sqrt(np.maximum(s2_res*(1+1/n+(t_all-t_mean)**2/max(S_xx,1e-10)),0))
lo_all=np.maximum(0,mu_all-t_crit*PI_all)
hi_all=np.minimum(5,mu_all+t_crit*PI_all)

ax2.fill_between(t_all,lo_all,hi_all,alpha=0.15,color=CM,label='95% prediction interval (df=1)')
ax2.plot(t_all,mu_all,'--',color=CM,lw=1.5,label='Linear forecast')
ax2.plot(t_obs,a_obs,'o',color=CM,ms=10,zorder=6,label='Observed (REAL data)',mec='black',mew=1)

# Threshold line
ax2.axhline(y=2.0,color='orange',ls='--',lw=1.2,alpha=0.8,label='α=2.0 (Coherent threshold)')
ax2.axhline(y=4.0,color='red',ls=':',lw=1.0,alpha=0.7,label='α=4.0 (Extreme coherent)')
ax2.axhspan(4.0,5.1,alpha=0.05,color='red')

ax2.text(2.5,4.7,'PREDICTION INTERVAL\nt₀.₉₇₅(df=1) = 12.71\n→ CI spans [0, 5]',
    fontsize=8,color=CM,ha='center',
    bbox=dict(boxstyle='round',facecolor='white',alpha=0.8,edgecolor=CM))

ax2.set_xlabel('Time within session (hours)',fontsize=10)
ax2.set_ylabel('α_GHZ',fontsize=10)
ax2.set_title('(b) Intra-Session Drift Forecast\n'
              f'slope = {sl:.3f}±{se:.3f} α/hr  R²={D["jul5_R2"]:.3f}  n=3',
              fontsize=10,fontweight='bold')
ax2.legend(fontsize=8,loc='upper left'); ax2.grid(True,alpha=0.25)
ax2.set_xlim(-0.3,5.0); ax2.set_ylim(0,5.2)
ax2.text(0.97,0.03,'Confidence: LOW\nLinear model, n=3 only',
    transform=ax2.transAxes,ha='right',va='bottom',fontsize=8,
    color='gray',style='italic')

plt.suptitle('Fig. 5  Phase Diagram and Drift Forecast F(χ, t)',
             fontsize=11,fontweight='bold',y=1.02)
plt.tight_layout()
plt.savefig('/home/claude/Fig5_PhaseDiagram_Forecast.pdf',format='pdf')
plt.savefig('/home/claude/Fig5_PhaseDiagram_Forecast.png')
plt.close()
print("  Saved.")

# ── FIGURE 6: VALIDATION ROADMAP ─────────────────────────
print("Fig 6: Validation Roadmap...")
fig,ax=plt.subplots(figsize=(11,5))
ax.set_xlim(0,11); ax.set_ylim(0,5.5); ax.axis('off')

stages=[
    (0.3,3.2,2.0,1.6,'CURRENT\n(this paper)\nn=4 sessions\n3 devices\nFeasibility\ndemonstration','#FFF3E0','black'),
    (3.0,3.2,2.0,1.6,'STAGE 2\nn=20 sessions\nConference\npaper\nWeak statistical\nclaim','#E8F5E9','darkgreen'),
    (5.7,3.2,2.0,1.6,'STAGE 3\nn=50 sessions\n+ algorithm\noutcomes\nPhys. Rev. Appl.','#E3F2FD','#1565c0'),
    (8.4,3.2,2.0,1.6,'STAGE 4\nn=100 sessions\nProspective\nvalidation\nnpj QI / PRX Q','#FCE4EC','darkred'),
]
for x,y,w,h,txt,fc,ec in stages:
    rect=Rectangle((x,y),w,h,lw=2,edgecolor=ec,facecolor=fc,alpha=0.9,zorder=3)
    ax.add_patch(rect)
    ax.text(x+w/2,y+h/2,txt,ha='center',va='center',fontsize=8,
            fontweight='bold' if 'CURRENT' in txt else 'normal',
            zorder=4,multialignment='center')

# Arrows
for xi in [2.3,5.0,7.7]:
    ax.annotate('',xy=(xi+0.7,4.0),xytext=(xi,4.0),
        arrowprops=dict(arrowstyle='-|>',color='#555',lw=1.5,mutation_scale=14))

# Required experiments below
exps=[
    (0.3,0.5,2.0,2.2,'What exists:\n\u2022 191 circuits\n\u2022 3 devices\n\u2022 3/3 TEST correct\n\u2022 Physics rules','#FFF9C4'),
    (3.0,0.5,2.0,2.2,'Needed:\n\u2022 Hourly α, 3 weeks\n\u2022 IBM calib. history\n\u2022 Exact binomial CI\n\u2022 >29% accuracy','#F1F8E9'),
    (5.7,0.5,2.0,2.2,'Needed:\n\u2022 50+ sessions\n\u2022 VQE/QAOA outcomes\n\u2022 ROC curves\n\u2022 AUC > 0.75','#E8EAF6'),
    (8.4,0.5,2.0,2.2,'Needed:\n\u2022 100+ sessions\n\u2022 Prospective design\n\u2022 Blind prediction\n\u2022 2nd hardware\nplatform','#FFEBEE'),
]
for x,y,w,h,txt,fc in exps:
    rect=Rectangle((x,y),w,h,lw=1,edgecolor='#999',facecolor=fc,alpha=0.85,zorder=3)
    ax.add_patch(rect)
    ax.text(x+w/2,y+h/2,txt,ha='center',va='center',fontsize=7.5,zorder=4,multialignment='center')
    ax.plot([x+w/2,x+w/2],[y+h,3.2],color='#aaa',lw=1,ls=':',zorder=2)

ax.text(5.5,5.3,'Validation Roadmap — Evidence Requirements for Publication',
    ha='center',va='center',fontsize=11,fontweight='bold')

plt.savefig('/home/claude/Fig6_ValidationRoadmap.pdf',format='pdf')
plt.savefig('/home/claude/Fig6_ValidationRoadmap.png')
plt.close()
print("  Saved.")

print("\nAll 6 figures generated successfully.")
print("PDFs ready for manuscript submission.")