"""
Extended Data Figure 2 — Cardiac differentiation lengthens the supra-compartment
                          dominant period 2.7-fold (H7 hESC → cardiac muscle).

REAL DATA ONLY. No rng. Reviewer-runnable.

Source CSV (BIOMNI 2026-04-04 H7 cardiac differentiation deliverable):
  - differentiation_downsampling_results.csv (5 days × 17 columns: Day 0-2-5-9-14)
       chrX dominant period + bootstrap-downsampled null at each timepoint
       (100 random subsamples of the Day 0 chrX peak set matched to target n_peaks).
  - cwt_h7_differentiation_summary.csv (Paul m=4 cross-check; secondary)

ENCODE accession ENCSR167QXA (H7 differentiation series):
  - Day 0: H7 hESC (ENCFF939OEJ; 21,422 chrX peaks)
  - Day 2: Mesodermal cell (ENCFF636AUG; 1,460 peaks)
  - Day 5: Cardiovascular progenitor (ENCFF496AXJ; 4,025 peaks)
  - Day 9: Cardiac myoblast (ENCFF874ORY; 4,815 peaks)
  - Day 14: Cardiac muscle cell (ENCFF701QBZ; 2,695 peaks)

Headline numbers reproduced (must match manuscript caption + Methods L236):
  - Day 0 → Day 14 dominant period: 15.63 → 42.24 Mbp (2.70-fold) ✓
  - Bootstrap downsampling z-scores: 3.67σ (Day 5) – 17.4σ (Day 2) ✓
  - All bootstrap p-values (one-tailed empirical): 0.00 (actual period > entire bootstrap distribution) ✓
  - sig95 fraction: Day 0 = 1.00, Day 14 = 0.984 (all timepoints exceed 5% AR(1) null)

Replaces 2026-04-21 v73 synthetic placeholder
(`figure_regen_v73/scripts/extended_data_figure_3_cardiac_tcell.py`, baked constant
24-26 Mbp across stages — does not match manuscript's 2.7-fold-increase claim).
"""
from __future__ import annotations
import sys, os
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

sys.path.insert(0, str(Path(__file__).resolve().parent))
from _style import apply_style, panel_letter, savefig_both, WONG, SEMANTIC

BUNDLE_ROOT = Path(__file__).resolve().parents[2]
DATA_DIR = BUNDLE_ROOT / '04_supplementary_data' / 'figure_data'
DOWNSAMPLING_CSV = DATA_DIR / 'differentiation_downsampling_results.csv'
SUMMARY_CSV = DATA_DIR / 'cwt_h7_differentiation_summary.csv'

OUT_BASENAMES = [
    # canonical name + alias matching the manuscript .tex \includegraphics path
    BUNDLE_ROOT / '02_figures' / 'extended_data' / 'extended_data_figure_2_differentiation',
    BUNDLE_ROOT / '02_figures' / 'extended_data' / 'extended_data_figure_3_differentiation',
    BUNDLE_ROOT / '03_figures_to_edit' / 'extended_data' / 'extended_data_figure_2_differentiation',
    BUNDLE_ROOT / '03_figures_to_edit' / 'extended_data' / 'extended_data_figure_3_differentiation',
]
DRAFT_BASENAME = BUNDLE_ROOT / '05_code' / 'figure_regen' / '_drafts' / 'ed2_differentiation_v1'

# ---------- Load ----------
df = pd.read_csv(DOWNSAMPLING_CSV)
# Coerce numeric columns
for col in ['actual_dom_period_Mbp', 'sub_mean_dom_period_Mbp', 'sub_sd_dom_period_Mbp',
            'sub_p2_5_Mbp', 'sub_p97_5_Mbp', 'actual_sig95_frac', 'z_score',
            'sub_mean_sig95_frac', 'sub_sd_sig95_frac']:
    df[col] = pd.to_numeric(df[col], errors='coerce')

cell_types = {0: 'H7 hESC', 2: 'Mesoderm', 5: 'Cardiac\nprogenitor', 9: 'Cardiac\nmyoblast', 14: 'Cardiac\nmuscle'}

# Verify
day0_period = df[df['day']==0]['actual_dom_period_Mbp'].iloc[0]
day14_period = df[df['day']==14]['actual_dom_period_Mbp'].iloc[0]
fold = day14_period / day0_period
z_min = df[df['day']>0]['z_score'].min()
z_max = df[df['day']>0]['z_score'].max()
print(f"[verify] Day 0 dominant period = {day0_period:.2f} Mbp (manuscript: 15.6)")
print(f"[verify] Day 14 dominant period = {day14_period:.2f} Mbp (manuscript: 42.2)")
print(f"[verify] fold change = {fold:.2f}× (manuscript: 2.7-fold)")
print(f"[verify] z-score range = {z_min:.1f}σ – {z_max:.1f}σ (manuscript: 3.7σ – 17.4σ)")
print(f"[verify] all empirical p < 0.01: {(df[df['day']>0]['p_value_one_tailed'].astype(float) < 0.01).all()}")

# ---------- Figure ----------
apply_style()
fig, axes = plt.subplots(1, 3, figsize=(13.0, 4.4),
                          gridspec_kw={'left': 0.06, 'right': 0.985, 'top': 0.86,
                                       'bottom': 0.21, 'wspace': 0.36})
ax_a, ax_b, ax_c = axes

xs = np.arange(len(df))
day_labels = [f'Day {int(d)}\n{cell_types[int(d)]}' for d in df['day']]

# === Panel a: actual dominant period across days ===
periods = df['actual_dom_period_Mbp'].values
boot_means = df['sub_mean_dom_period_Mbp'].fillna(np.nan).values
boot_sds = df['sub_sd_dom_period_Mbp'].fillna(np.nan).values
# Plot actual period as gold bars (key finding)
ax_a.bar(xs, periods, color=WONG['gold'], edgecolor='white', linewidth=0.5,
         label='actual dominant period', zorder=3)
# Overlay bootstrap mean ± SD as black error bars (where defined)
mask = ~np.isnan(boot_means)
ax_a.errorbar(xs[mask], boot_means[mask], yerr=boot_sds[mask],
              fmt='o', color=WONG['black'], markersize=7, capsize=4,
              markerfacecolor='white', markeredgewidth=1.4,
              elinewidth=1.0, zorder=5,
              label='bootstrap (Day 0 peak set\ndownsampled to Day-N count)')
# Annotate values
for x, p in zip(xs, periods):
    ax_a.text(x, p + 1.5, f'{p:.1f}', ha='center', va='bottom',
              fontsize=8.5, fontweight='bold')
ax_a.set_xticks(xs); ax_a.set_xticklabels(day_labels, fontsize=7.5)
ax_a.set_ylabel('chrX dominant period (Mbp)')
ax_a.set_ylim(0, max(periods)*1.18)
ax_a.legend(loc='upper left', fontsize=8)
ax_a.set_title(f'Period lengthens {fold:.1f}× (15.6 → 42.2 Mbp)',
               loc='left', pad=8, fontsize=10)
panel_letter(fig, ax_a, 'a')

# === Panel b: bootstrap z-scores ===
mask_b = df['day'] > 0
xs_b = xs[mask_b]
zs = df.loc[mask_b, 'z_score'].astype(float).values
ax_b.bar(xs_b, zs, color=WONG['blue'], edgecolor='white', linewidth=0.5)
ax_b.axhline(2, color='0.5', lw=0.7, ls='--', label='2σ threshold')
for x, z in zip(xs_b, zs):
    ax_b.text(x, z + 0.4, f'{z:.1f}σ', ha='center', va='bottom',
              fontsize=8.5, fontweight='bold')
ax_b.set_xticks(xs_b); ax_b.set_xticklabels([day_labels[i] for i in xs_b], fontsize=7.5)
ax_b.set_ylabel('z-score (actual vs bootstrap)')
ax_b.set_ylim(0, max(zs)*1.20)
ax_b.legend(loc='upper right', fontsize=8)
ax_b.set_title('Not a peak-count artefact (z = 3.7 – 17.4σ)',
               loc='left', pad=8, fontsize=10)
panel_letter(fig, ax_b, 'b')

# === Panel c: spectral concentration (sig95) per day ===
sig95 = df['actual_sig95_frac'].values
ax_c.bar(xs, sig95, color=WONG['vermil'], edgecolor='white', linewidth=0.5)
ax_c.axhline(0.05, color='0.5', lw=0.7, ls='--', label='5% AR(1) null')
for x, s in zip(xs, sig95):
    ax_c.text(x, s + 0.02, f'{s:.2f}', ha='center', va='bottom',
              fontsize=8.5, fontweight='bold')
ax_c.set_xticks(xs); ax_c.set_xticklabels(day_labels, fontsize=7.5)
ax_c.set_ylabel('sig95 fraction')
ax_c.set_ylim(0, 1.15)
ax_c.legend(loc='lower right', fontsize=8)
ax_c.set_title('All timepoints above 5% AR(1) null',
               loc='left', pad=8, fontsize=10)
panel_letter(fig, ax_c, 'c')

# Save
DRAFT_BASENAME.parent.mkdir(parents=True, exist_ok=True)
savefig_both(fig, DRAFT_BASENAME)
print(f"[wrote draft] {DRAFT_BASENAME}.png + .svg")
if os.environ.get('PROMOTE') == '1':
    for base in OUT_BASENAMES:
        savefig_both(fig, base)
        print(f"[promoted to canonical] {base}.png + .svg")
plt.close(fig)
