"""Compute Benjamini-Hochberg FDR for all pre-specified tests."""
import sys, io
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')
import numpy as np

# All pre-specified hypothesis tests from the paper
tests = [
    # Llama baseline (Section 4.4)
    ("loop ↔ autocorr (baseline)", 0.002),
    ("surge ↔ max_norm (baseline)", 0.002),
    # Llama steered (Section 4.5)  
    ("shimmer ↔ norm_std (steered)", 0.005),
    ("shimmer paired Δ", 0.002),
    ("surge ↔ max_norm (steered)", 0.0005),
    ("surge paired Δ", 0.001),
    # Qwen baseline (Section 4.6)
    ("mirror ↔ spectral (Qwen)", 0.0001),  # p < 0.0001
    ("expand ↔ spectral (Qwen)", 0.0001),
    ("resonance ↔ max_norm (Qwen)", 0.0001),
]

names = [t[0] for t in tests]
pvals = np.array([t[1] for t in tests])
m = len(pvals)

# BH procedure
sorted_idx = np.argsort(pvals)
sorted_p = pvals[sorted_idx]
sorted_names = [names[i] for i in sorted_idx]

# q-values
qvals = np.zeros(m)
for i in range(m):
    rank = i + 1
    qvals[i] = sorted_p[i] * m / rank

# Enforce monotonicity (from bottom up)
for i in range(m-2, -1, -1):
    qvals[i] = min(qvals[i], qvals[i+1])

print(f"Benjamini-Hochberg FDR correction across {m} pre-specified tests")
print(f"{'Test':<40} {'p-value':>10} {'rank':>5} {'q-value':>10} {'Sig (q<0.05)':>12}")
print("-" * 80)
for i in range(m):
    sig = "YES" if qvals[i] < 0.05 else "NO"
    print(f"{sorted_names[i]:<40} {sorted_p[i]:>10.4f} {i+1:>5} {qvals[i]:>10.4f} {sig:>12}")