"""Generate all Paper 2 figures at publication quality"""
import json, numpy as np, matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.patches import FancyArrowPatch
from scipy.optimize import curve_fit

# ── Global style ──
FS=12; TS=11; LS=11; TTS=13
plt.rcParams.update({
    'font.family':'serif','font.size':FS,
    'axes.labelsize':FS,'axes.titlesize':TTS,
    'xtick.labelsize':TS,'ytick.labelsize':TS,
    'legend.fontsize':LS,'lines.linewidth':1.8,
    'lines.markersize':7,'figure.dpi':300,
    'savefig.dpi':300,'savefig.bbox':'tight',
})

CK='#2166ac'; CF='#1a9641'; CM='#d73027'; CS='#ff7f00'

def kr_model(K,C0,alpha,R): return C0/(K**alpha)+R
def fit_kr(x,y):
    K=1.0/np.array(x,dtype=float)
    y=np.array(y,dtype=float)
    try:
        popt,_=curve_fit(kr_model,K,y,p0=[0.1,0.5,0.01],
            bounds=([0,0.001,0],[50,15,1]),method='trf',maxfev=20000)
        yp=kr_model(K,*popt)
        ss_res=np.sum((y-yp)**2); ss_tot=np.sum((y-np.mean(y))**2)
        R2=1.0-ss_res/ss_tot if ss_tot>0 else 0.0
        return popt, R2
    except: return None, None

def kr_curve(x_range, popt):
    return kr_model(1.0/np.array(x_range), *popt)

# Load data
with open('/mnt/user-data/uploads/ibm_hardware_v2_results.json') as f: v2=json.load(f)
with open('/mnt/user-data/uploads/kr_marrakesh_results.json') as f: mar=json.load(f)
with open('/mnt/user-data/uploads/kr_hardware_experiments.json') as f: hw=json.load(f)
with open('/mnt/user-data/uploads/ibm_8_regimes_results.json') as f: reg8=json.load(f)
with open('/mnt/user-data/uploads/kr_beta_extraction_results.json') as f: beta=json.load(f)
with open('/mnt/user-data/uploads/KR_Quantum_Comprehensive_Results.json') as f: sim=json.load(f)
with open('/mnt/user-data/uploads/ibm_hardproblems_results.json') as f: hp=json.load(f)

# ── Extract data ──
k_ghz_n=[int(n) for n in v2['backends']['ibm_kingston']['exp1_ghz_hardware']]
k_ghz_inf=[v2['backends']['ibm_kingston']['exp1_ghz_hardware'][str(n)]['infidelity'] for n in k_ghz_n]
f_ghz_n=[int(n) for n in v2['backends']['ibm_fez']['exp1_ghz_hardware']]
f_ghz_inf=[v2['backends']['ibm_fez']['exp1_ghz_hardware'][str(n)]['infidelity'] for n in f_ghz_n]
m_ghz_n=[int(n) for n in mar['exp1_ghz_scaling']['data']]
m_ghz_inf=[mar['exp1_ghz_scaling']['data'][str(n)]['infidelity'] for n in m_ghz_n]

k_dep_d=[int(d) for d in v2['backends']['ibm_kingston']['exp2_depth_hardware']]
k_dep_inf=[v2['backends']['ibm_kingston']['exp2_depth_hardware'][str(d)]['infidelity'] for d in k_dep_d]
f_dep_d=[int(d) for d in v2['backends']['ibm_fez']['exp2_depth_hardware']]
f_dep_inf=[v2['backends']['ibm_fez']['exp2_depth_hardware'][str(d)]['infidelity'] for d in f_dep_d]
m_dep_d=[int(d) for d in mar['exp2_depth_scaling']['data']]
m_dep_inf=[mar['exp2_depth_scaling']['data'][str(d)]['infidelity'] for d in m_dep_d]

k_amp_r=[int(r) for r in v2['backends']['ibm_kingston']['exp3_noise_amplification']]
k_amp_inf=[v2['backends']['ibm_kingston']['exp3_noise_amplification'][str(r)]['infidelity'] for r in k_amp_r]
m_amp_r=[int(k) for k in hw['h3_noise_amplification']['data']]
m_amp_inf=[hw['h3_noise_amplification']['data'][str(r)]['infidelity'] for r in m_amp_r]

h1_runs=hw['h1_temporal_drift']['runs']
h1_infs=[r['infidelity'] for r in h1_runs]
h1_times=[r['timestamp'] for r in h1_runs]

k_shot_s=[int(s) for s in v2['backends']['ibm_kingston']['exp4_shot_noise']]
k_shot_inf=[v2['backends']['ibm_kingston']['exp4_shot_noise'][str(s)]['infidelity'] for s in k_shot_s]
m_shot_s=[int(k) for k in hw['h4_shot_noise']['data']]
m_shot_inf=[hw['h4_shot_noise']['data'][str(s)]['infidelity'] for s in m_shot_s]

# ── FIG 1: THE KEY DIVERGENCE FIGURE ──────────────────────
print("Generating Figure 1: Key Divergence...")
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(11,5))

# Left: IBM T2 vs K-R alpha timeline
dates=['Jul 3\n2026','Jul 5\n2026 (Start)','Jul 5\n2026 (Mid)','Jul 5\n2026 (End)']
t2_vals=[303.1,339.6,339.6,339.6]
alpha_vals=[1.406,3.811,3.1787,3.7792]

ax1b=ax1.twinx()
l1,=ax1.plot([0,1,2,3],t2_vals,'o-',color='steelblue',linewidth=2,
             markersize=9,label='IBM T₂ (μs)')
l2,=ax1b.plot([0,1,2,3],alpha_vals,'s--',color=CM,linewidth=2,
              markersize=9,label='K-R α (circuit-level)')
ax1.set_xlabel('Session',fontsize=FS)
ax1.set_ylabel('IBM T₂ (μs)',fontsize=FS,color='steelblue')
ax1b.set_ylabel('K-R Scaling Exponent α',fontsize=FS,color=CM)
ax1.tick_params(axis='y',labelcolor='steelblue')
ax1b.tick_params(axis='y',labelcolor=CM)
ax1.set_xticks([0,1,2,3]); ax1.set_xticklabels(dates,fontsize=10)
ax1.set_ylim(200,420); ax1b.set_ylim(0,5)
ax1.axhspan(200,420,alpha=0.05,color='steelblue')
ax1.annotate('T₂ improves\n(+12%)',xy=(1,339.6),xytext=(0.3,370),
    fontsize=10,color='steelblue',
    arrowprops=dict(arrowstyle='->',color='steelblue'))
ax1b.annotate('α increases\n(+171%)',xy=(1,3.811),xytext=(1.5,4.3),
    fontsize=10,color=CM,
    arrowprops=dict(arrowstyle='->',color=CM))
lines=[l1,l2]
labs=[l.get_label() for l in lines]
ax1.legend(lines,labs,loc='lower left',fontsize=10)
ax1.set_title('(a) Divergence: IBM T₂ vs K-R α',fontsize=TTS,fontweight='bold')
ax1.grid(True,alpha=0.3)

# Right: GHZ infidelity Jul3 vs Jul5
n_all=[2,5,10,15,20,25,30]
inf_jul3=[0.020,0.056,0.164,0.266,0.602,0.772,0.844]
inf_jul5_start=[0.204,0.295,0.412,0.531,0.720,0.820,0.900]

popt3,R2_3=fit_kr(n_all,inf_jul3)
popt5,R2_5=fit_kr(n_all,inf_jul5_start)

n_sm=np.linspace(2,30,200)
ax2.plot(n_all,inf_jul3,'o',color=CK,markersize=8,label=f'Jul 3 (α=1.41, R²=0.97)')
if popt3 is not None:
    ax2.plot(n_sm,kr_curve(n_sm,popt3),'--',color=CK,alpha=0.7)
ax2.plot(n_all,inf_jul5_start,'s',color=CM,markersize=8,label=f'Jul 5 (α=3.81, R²=0.98)')
if popt5 is not None:
    ax2.plot(n_sm,kr_curve(n_sm,popt5),'--',color=CM,alpha=0.7)

ax2.annotate('+920%\nat n=2',xy=(2,0.204),xytext=(5,0.35),fontsize=10,
    color=CM,arrowprops=dict(arrowstyle='->',color=CM))
ax2.set_xlabel('Number of Qubits n',fontsize=FS)
ax2.set_ylabel('GHZ Infidelity',fontsize=FS)
ax2.set_title('(b) GHZ Infidelity: Jul 3 vs Jul 5',fontsize=TTS,fontweight='bold')
ax2.legend(fontsize=10); ax2.grid(True,alpha=0.3); ax2.set_xlim(0,32)

plt.suptitle('The IBM T₂ / K-R α Divergence on ibm_marrakesh',
             fontsize=TTS+1,fontweight='bold',y=1.02)
plt.tight_layout()
plt.savefig('/home/claude/Fig1_Divergence.pdf',format='pdf')
plt.savefig('/home/claude/Fig1_Divergence.png',dpi=300)
plt.close()
print("  Fig 1 saved.")

# ── FIG 2: THREE-DEVICE FINGERPRINT ───────────────────────
print("Generating Figure 2: Three-Device Fingerprint...")
fig,axes=plt.subplots(1,3,figsize=(14,5))

# (a) GHZ scaling all 3 devices
ax=axes[0]
poptk,_=fit_kr(k_ghz_n,k_ghz_inf)
poptf,_=fit_kr(f_ghz_n,f_ghz_inf)
poptm,_=fit_kr(m_ghz_n,m_ghz_inf)
n_sm=np.linspace(2,30,200)
ax.plot(k_ghz_n,k_ghz_inf,'o',color=CK,markersize=7,label=f'ibm\\_kingston (α=1.54)')
ax.plot(np.linspace(2,10,100),kr_curve(np.linspace(2,10,100),poptk),'--',color=CK,alpha=0.7)
ax.plot(f_ghz_n,f_ghz_inf,'s',color=CF,markersize=7,label=f'ibm\\_fez (α=0.91)')
ax.plot(np.linspace(2,10,100),kr_curve(np.linspace(2,10,100),poptf),'--',color=CF,alpha=0.7)
ax.plot(m_ghz_n,m_ghz_inf,'^',color=CM,markersize=8,label=f'ibm\\_marrakesh (α=1.41)')
ax.plot(n_sm,kr_curve(n_sm,poptm),'--',color=CM,alpha=0.7)
ax.set_xlabel('Number of Qubits n',fontsize=FS)
ax.set_ylabel('GHZ Infidelity',fontsize=FS)
ax.set_title('(a) GHZ Scaling: Three Devices\n(n=2 to 30)',fontsize=TTS,fontweight='bold')
ax.legend(fontsize=9); ax.grid(True,alpha=0.3)
ax.annotate('n=30',xy=(30,0.844),xytext=(24,0.65),fontsize=9,color=CM,
    arrowprops=dict(arrowstyle='->',color=CM))

# (b) Depth scaling
ax=axes[1]
poptk,_=fit_kr(k_dep_d,k_dep_inf)
poptf,_=fit_kr(f_dep_d,f_dep_inf)
poptm,_=fit_kr(m_dep_d,m_dep_inf)
d_sm=np.linspace(1,64,200)
ax.plot(k_dep_d,k_dep_inf,'o',color=CK,markersize=7,label=f'ibm\\_kingston (α=0.51)')
ax.plot(d_sm,kr_curve(d_sm,poptk),'--',color=CK,alpha=0.7)
ax.plot(f_dep_d,f_dep_inf,'s',color=CF,markersize=7,label=f'ibm\\_fez (α=0.32)')
ax.plot(d_sm,kr_curve(d_sm,poptf),'--',color=CF,alpha=0.7)
ax.plot(m_dep_d,m_dep_inf,'^',color=CM,markersize=8,label=f'ibm\\_marrakesh (α=1.49)')
ax.plot(d_sm,kr_curve(d_sm,poptm),'--',color=CM,alpha=0.7)
ax.set_xlabel('Circuit Depth d',fontsize=FS)
ax.set_ylabel('Clifford Infidelity',fontsize=FS)
ax.set_title('(b) Clifford Depth Scaling',fontsize=TTS,fontweight='bold')
ax.legend(fontsize=9); ax.grid(True,alpha=0.3)

# (c) Fingerprint comparison bar chart
ax=axes[2]
devices=['kingston','fez','marrakesh']
ghz_alphas=[1.541,0.907,1.406]
dep_alphas=[0.507,0.321,1.493]
x=np.arange(3); w=0.35
colors=[CK,CF,CM]
b1=ax.bar(x-w/2,dep_alphas,w,label='Depth α',color=colors,alpha=0.85,edgecolor='black',linewidth=0.8)
b2=ax.bar(x+w/2,ghz_alphas,w,label='GHZ α',color=colors,alpha=0.55,edgecolor='black',linewidth=0.8,hatch='//')
for bar,val in zip(b1,dep_alphas):
    ax.text(bar.get_x()+bar.get_width()/2,bar.get_height()+0.03,f'{val:.3f}',
            ha='center',va='bottom',fontsize=9,fontweight='bold')
for bar,val in zip(b2,ghz_alphas):
    ax.text(bar.get_x()+bar.get_width()/2,bar.get_height()+0.03,f'{val:.3f}',
            ha='center',va='bottom',fontsize=9,fontweight='bold')
ax.axhline(y=1.0,color='red',linestyle='--',linewidth=1.5,alpha=0.7,label='α=1 boundary')
ax.set_xticks(x)
ax.set_xticklabels([f'ibm\\_\n{d}' for d in devices],fontsize=9)
ax.set_ylabel('K-R Scaling Exponent α',fontsize=FS)
ax.set_title('(c) Device Fingerprint\n(Three Heron r2 Processors)',fontsize=TTS,fontweight='bold')
ax.legend(fontsize=9); ax.grid(True,alpha=0.3,axis='y')

plt.suptitle('Three-Device K-R Fingerprint: ibm_kingston, ibm_fez, ibm_marrakesh',
             fontsize=TTS+1,fontweight='bold',y=1.02)
plt.tight_layout()
plt.savefig('/home/claude/Fig2_ThreeDevice.pdf',format='pdf')
plt.savefig('/home/claude/Fig2_ThreeDevice.png',dpi=300)
plt.close()
print("  Fig 2 saved.")

# ── FIG 3: INTRA-SESSION DRIFT + PROOF 4 ──────────────────
print("Generating Figure 3: Drift and Proof 4...")
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(11,5))

# (a) Intra-session drift
runs=[1,2,3]
drift_alphas=[2.9868,3.1787,3.7792]
drift_infs=h1_infs
ax1b=ax1.twinx()
l1,=ax1.plot(runs,drift_alphas,'o-',color=CM,linewidth=2.5,markersize=10,
    label='GHZ α (5-point suite)')
ax1.fill_between(runs,[a-0.1 for a in drift_alphas],
                 [a+0.1 for a in drift_alphas],alpha=0.2,color=CM)
mean_inf=np.mean(drift_infs)
l2,=ax1b.plot([1,2,3,4,5],drift_infs,'s--',color='gray',linewidth=1.5,markersize=8,
    label='Single GHZ-5 inf.')
ax1b.axhline(y=mean_inf,color='gray',linestyle=':',alpha=0.7)
ax1.set_xlabel('Run Number',fontsize=FS)
ax1.set_ylabel('GHZ Scaling α (5-point)',fontsize=FS,color=CM)
ax1b.set_ylabel('GHZ-5 Infidelity (single-point)',fontsize=FS,color='gray')
ax1.tick_params(axis='y',labelcolor=CM)
ax1b.tick_params(axis='y',labelcolor='gray')
ax1.set_title('(a) Intra-Session Drift Detection\n(ibm_marrakesh, Jul 5)',
              fontsize=TTS,fontweight='bold')
ax1.annotate('Monotonic drift\ndetected by K-R',xy=(3,3.7792),xytext=(1.5,3.85),
    fontsize=10,color=CM,arrowprops=dict(arrowstyle='->',color=CM))
ax1.annotate('Single-point:\nCV=6.2% (stable?)',xy=(3,0.063),xytext=(3.2,0.073),
    fontsize=9,color='gray',transform=ax1.transData,
    xycoords='data',textcoords='data')
lines=[l1,l2]
ax1.legend(lines,[l.get_label() for l in lines],loc='upper left',fontsize=10)
ax1.grid(True,alpha=0.3)

# (b) Proof 4: alpha_cliff + alpha_echo = 2
ax2.barh(['Clifford depth\n(α_cliff)','Echo circuit\n(α_echo)'],
         [1.7853,0.3544],color=[CK,CF],edgecolor='black',height=0.5)
ax2.axvline(x=1.0,color='red',linestyle='--',linewidth=1.5,label='α=1 boundary')
ax2.text(1.85,0,'α_cliff=1.785',ha='center',va='center',fontsize=11,
         fontweight='bold',color=CK)
ax2.text(0.18,1,'α_echo=0.354',ha='center',va='center',fontsize=11,
         fontweight='bold',color=CF)
ax2.text(2.5,0.5,f'α_cliff + α_echo\n= {1.7853+0.3544:.4f}\n≈ 2.0 (theory)',
    ha='center',va='center',fontsize=12,fontweight='bold',color='darkgreen',
    bbox=dict(boxstyle='round',facecolor='lightgreen',alpha=0.4))
ax2.set_xlabel('K-R Scaling Exponent α',fontsize=FS)
ax2.set_title('(b) Proof 4 Confirmation\nα_cliff + α_echo = 2.1397 (7% error)',
              fontsize=TTS,fontweight='bold')
ax2.legend(fontsize=10); ax2.grid(True,alpha=0.3,axis='x')
ax2.set_xlim(0,3.2)

plt.tight_layout()
plt.savefig('/home/claude/Fig3_DriftProof4.pdf',format='pdf')
plt.savefig('/home/claude/Fig3_DriftProof4.png',dpi=300)
plt.close()
print("  Fig 3 saved.")

# ── FIG 4: CIRCUIT vs PULSE T2 + BETA EXTRACTION ──────────
print("Generating Figure 4: T2 Discrepancy...")
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(11,5))

# (a) Per-qubit T2 comparison
m1=beta['method1_backend_properties']
qubits=list(range(10))
ibm_t2s=[m1['qubit_data'][str(q)]['T2_us'] for q in qubits]
ibm_t1s=[m1['qubit_data'][str(q)]['T1_us'] for q in qubits]
circuit_t2_q0=9.0
circuit_t2_q2=30.0

x=np.arange(10); w=0.35
ax1.bar(x-w/2,ibm_t2s,w,label='IBM T₂ (pulse-level)',
        color='steelblue',alpha=0.8,edgecolor='navy')
circuit_vals=[circuit_t2_q0 if q==0 else circuit_t2_q2 if q==2 else np.nan
              for q in qubits]
ax1.bar([0,2],[circuit_t2_q0,circuit_t2_q2],w+0.1,
        label='Circuit T₂ (Ramsey)',color=CM,alpha=0.8,
        edgecolor='darkred',linewidth=2)
for q,v in [(0,circuit_t2_q0),(2,circuit_t2_q2)]:
    ratio=ibm_t2s[q]/v
    ax1.text(q,max(ibm_t2s[q],v)+10,f'{ratio:.1f}×',
             ha='center',va='bottom',fontsize=10,
             fontweight='bold',color='red')
ax1.set_xlabel('Qubit Index',fontsize=FS)
ax1.set_ylabel('T₂ (μs)',fontsize=FS)
ax1.set_title('(a) IBM T₂ vs Circuit T₂\nDiscrepancy up to 11×',
              fontsize=TTS,fontweight='bold')
ax1.set_xticks(x); ax1.legend(fontsize=10); ax1.grid(True,alpha=0.3,axis='y')

# (b) Delay Ramsey coherence decay showing oscillation
m2=beta['method2_delay_ramsey']
taus_all=[1,5,10,20,50,100,200]
cohs=[m2['data'][str(t)]['coherence'] for t in taus_all]
t_fine=np.linspace(1,50,200)

ax2.plot(taus_all,cohs,'o',color=CM,markersize=9,zorder=5,label='Measured coherence')
ax2.axhline(y=0,color='black',linestyle='-',linewidth=0.8,alpha=0.5)
ax2.axhline(y=1/np.e,color='gray',linestyle='--',linewidth=1.5,
    alpha=0.7,label='1/e (T₂ threshold)')
# Fit Ramsey fringe to first 6 points
from scipy.optimize import curve_fit as cf2
def ramsey(t,A,T2,f,phi):
    return A*np.exp(-t/T2)*np.cos(2*np.pi*f*t+phi)
try:
    popt_r,_=cf2(ramsey,[1,5,10,20,30,50],cohs[:6],
        p0=[1.0,15.0,0.025,0.0],
        bounds=([0.5,1,0.001,-np.pi],[1.5,200,0.5,np.pi]),maxfev=10000)
    ax2.plot(t_fine,ramsey(t_fine,*popt_r),'-',color='navy',
        linewidth=2,label=f'Ramsey fit: T₂≈{popt_r[1]:.0f}μs, f≈{popt_r[2]*1000:.0f}kHz')
    ax2.annotate(f'Detuning:\n~{popt_r[2]*1000:.0f} kHz',
        xy=(20,-0.1),xytext=(40,0.2),fontsize=10,color='navy',
        arrowprops=dict(arrowstyle='->',color='navy'))
except: pass
ax2.annotate(f'IBM T₂=32μs\n(Q0 pulse-level)',xy=(32,1/np.e),
    xytext=(80,0.6),fontsize=10,color='steelblue',
    arrowprops=dict(arrowstyle='->',color='steelblue'))
ax2.set_xlabel('Delay τ (μs)',fontsize=FS)
ax2.set_ylabel('Ramsey Coherence C(τ)',fontsize=FS)
ax2.set_title('(b) Ramsey Circuit Coherence Decay\n(Oscillation = ZZ crosstalk detuning)',
              fontsize=TTS,fontweight='bold')
ax2.legend(fontsize=9); ax2.grid(True,alpha=0.3)
ax2.set_xlim(0,210)

plt.tight_layout()
plt.savefig('/home/claude/Fig4_T2Discrepancy.pdf',format='pdf')
plt.savefig('/home/claude/Fig4_T2Discrepancy.png',dpi=300)
plt.close()
print("  Fig 4 saved.")

# ── FIG 5: PHASE DIAGRAM (updated for paper 2) ────────────
print("Generating Figure 5: Phase Diagram...")
fig,ax=plt.subplots(figsize=(8,7))

# Zones
ax.axhspan(0.8,1.06,xmin=0,xmax=1/4.5,color='#d1e8ff',alpha=0.35,zorder=0)
ax.axhspan(0.8,1.06,xmin=1/4.5,xmax=1.0,color='#d4edda',alpha=0.35,zorder=0)
ax.axhspan(-0.05,0.5,color='#fde0dc',alpha=0.35,zorder=0)
ax.text(0.25,0.95,'Stochastic',style='italic',fontsize=11,
    color='#1565c0',ha='center',transform=ax.transAxes)
ax.text(0.70,0.95,'Structured\n(coherent-like)',style='italic',fontsize=11,
    color='#1b5e20',ha='center',transform=ax.transAxes)
ax.text(0.15,0.20,'Model\nBreakdown',style='italic',fontsize=11,
    color='#b71c1c',ha='center',transform=ax.transAxes)
ax.axhline(y=0.8,color='gray',linestyle='--',linewidth=1.0,alpha=0.6)
ax.axhline(y=0.5,color='gray',linestyle='--',linewidth=1.0,alpha=0.6)
ax.axvline(x=1.0,color='gray',linestyle='--',linewidth=1.0,alpha=0.6)

# Hardware points
k_pts=[(1.541,0.975,'GHZ C1'),(0.507,0.813,'Dep C1'),
       (1.636,0.977,'GHZ C2'),(0.533,0.998,'Clif C2'),
       (1.190,1.000,'Amp C2'),(0.609,0.894,'Echo C2')]
f_pts=[(0.907,0.986,'GHZ C1'),(0.321,0.812,'Dep C1'),
       (2.582,0.912,'GHZ C2'),(0.690,0.975,'Echo')]
m_pts_jul3=[(1.406,0.967,'GHZ Jul3'),(1.493,0.900,'Dep Jul3'),
            (0.599,0.996,'Amp Jul3')]
m_pts_jul5=[(3.811,0.980,'GHZ Jul5'),(1.785,0.968,'Clif Jul5'),
            (0.354,0.946,'Echo Jul5')]

for alpha,r2,lbl in k_pts:
    ax.scatter(alpha,r2,marker='*',s=130,color=CK,zorder=6,
               edgecolors='navy',linewidths=0.5)
for alpha,r2,lbl in f_pts:
    ax.scatter(alpha,r2,marker='X',s=90,color=CF,zorder=6,
               edgecolors='darkgreen',linewidths=0.5)
for alpha,r2,lbl in m_pts_jul3:
    ax.scatter(alpha,r2,marker='^',s=100,color=CM,zorder=7,
               edgecolors='darkred',linewidths=0.8,alpha=0.6)
for alpha,r2,lbl in m_pts_jul5:
    ax.scatter(alpha,r2,marker='p',s=130,color='darkred',zorder=8,
               edgecolors='black',linewidths=1.0)

# Arrow showing Jul3→Jul5 drift
ax.annotate('',xy=(3.811,0.980),xytext=(1.406,0.967),
    arrowprops=dict(arrowstyle='->>',color='darkred',lw=2.0))
ax.text(2.5,0.90,'Jul 3→Jul 5\nα: 1.41→3.81',fontsize=10,
    color='darkred',ha='center',fontweight='bold')

# Simulator
sim_pts=[(1.825,0.9999),(0.522,0.961),(2.998,0.710),(0.010,0.0),(0.049,0.305)]
for alpha,r2 in sim_pts:
    ax.scatter(alpha,r2,marker='D',s=70,color=CS,zorder=5,
               edgecolors='darkorange',linewidths=0.5)

from matplotlib.lines import Line2D
legend_elements=[
    Line2D([0],[0],marker='*',color='w',markerfacecolor=CK,markersize=12,
           markeredgecolor='navy',label='ibm\\_kingston'),
    Line2D([0],[0],marker='X',color='w',markerfacecolor=CF,markersize=11,
           markeredgecolor='darkgreen',label='ibm\\_fez'),
    Line2D([0],[0],marker='^',color='w',markerfacecolor=CM,markersize=11,
           markeredgecolor='darkred',alpha=0.6,label='ibm\\_marrakesh (Jul 3)'),
    Line2D([0],[0],marker='p',color='w',markerfacecolor='darkred',markersize=12,
           markeredgecolor='black',label='ibm\\_marrakesh (Jul 5)'),
    Line2D([0],[0],marker='D',color='w',markerfacecolor=CS,markersize=10,
           markeredgecolor='darkorange',label='Simulator'),
]
ax.legend(handles=legend_elements,loc='lower right',fontsize=10)
ax.set_xlabel('K-R Scaling Exponent α',fontsize=FS)
ax.set_ylabel('Goodness of Fit R²',fontsize=FS)
ax.set_title('(α, R²) Phase Diagram: Hardware + Simulator\n'
             'Arrow shows ibm_marrakesh temporal drift Jul 3→Jul 5',
             fontsize=TTS,fontweight='bold')
ax.set_xlim(-0.1,4.5); ax.set_ylim(-0.1,1.07)
ax.grid(True,alpha=0.3)
plt.tight_layout()
plt.savefig('/home/claude/Fig5_PhaseDiagram.pdf',format='pdf')
plt.savefig('/home/claude/Fig5_PhaseDiagram.png',dpi=300)
plt.close()
print("  Fig 5 saved.")

# ── FIG 6: NOISE AMPLIFICATION + SHOT NOISE ───────────────
print("Generating Figure 6: Amplification and Shot Noise...")
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(11,5))

# (a) Noise amplification 3 devices
x=np.arange(4); w=0.35
ax1.bar(x-w/2,k_amp_inf,w,label=f'ibm\\_kingston (α=1.190)',
        color=CK,alpha=0.85,edgecolor='navy')
ax1.bar(x+w/2,m_amp_inf,w,label=f'ibm\\_marrakesh (α=0.599)',
        color=CM,alpha=0.85,edgecolor='darkred')
for i,(v1,v2) in enumerate(zip(k_amp_inf,m_amp_inf)):
    ax1.text(i-w/2,v1+0.01,f'{v1:.3f}',ha='center',va='bottom',fontsize=9,color='navy')
    ax1.text(i+w/2,v2+0.01,f'{v2:.3f}',ha='center',va='bottom',fontsize=9,color='darkred')
ax1.set_xlabel('Repetitions (Gate Fold Factor)',fontsize=FS)
ax1.set_ylabel('Infidelity',fontsize=FS)
ax1.set_xticks(x); ax1.set_xticklabels(['1×','3×','5×','7×'],fontsize=FS)
ax1.set_title('(a) Noise Amplification\nkingston (Structured) vs marrakesh (Stochastic)',
              fontsize=TTS,fontweight='bold')
ax1.legend(fontsize=10); ax1.grid(True,alpha=0.3,axis='y')

# (b) Shot noise stability
ax2.plot(k_shot_s,k_shot_inf,'o-',color=CK,linewidth=2,markersize=9,
    label=f'ibm\\_kingston (CV={100*np.std(k_shot_inf)/np.mean(k_shot_inf):.1f}%)')
ax2.plot(m_shot_s,m_shot_inf,'s-',color=CM,linewidth=2,markersize=9,
    label=f'ibm\\_marrakesh (CV=5.7%)')
ax2.axhline(y=np.mean(k_shot_inf),color=CK,linestyle=':',alpha=0.6)
ax2.axhline(y=np.mean(m_shot_inf),color=CM,linestyle=':',alpha=0.6)
ax2.set_xlabel('Shot Count',fontsize=FS)
ax2.set_ylabel('GHZ-5 Infidelity',fontsize=FS)
ax2.set_xscale('log',base=2)
ax2.set_xticks([512,1024,4096]); ax2.set_xticklabels(['512','1024','4096'])
ax2.set_title('(b) Shot Noise Stability\n(Monotonic increase on marrakesh = drift)',
              fontsize=TTS,fontweight='bold')
ax2.legend(fontsize=10); ax2.grid(True,alpha=0.3)
plt.tight_layout()
plt.savefig('/home/claude/Fig6_AmpShot.pdf',format='pdf')
plt.savefig('/home/claude/Fig6_AmpShot.png',dpi=300)
plt.close()
print("  Fig 6 saved.")

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