# Reviewer diagnostics on results/per_map.csv: who the non-positive maps are, and slopes that
# don't depend on the 1e-4 floor (median regression; rank correlation within collections).
import sys, numpy as np, pandas as pd, statsmodels.formula.api as smf
from scipy.stats import spearmanr
D = pd.read_csv(sys.argv[1])
A = D[D.status == "ok"].copy()
A["nonpos"] = A.delta2 <= 0
A["log_N"] = np.log(A.N)
A["fwhm_missing"] = A.fwhm_mm.isna().astype(int); A["fwhm_c"] = A.fwhm_mm.fillna(0)
print("analysed", len(A), "non-positive", int(A.nonpos.sum()))
A["Nbin"] = pd.cut(A.N, [9, 20, 40, 100, 10000], labels=["10-20", "21-40", "41-100", ">100"])
print(A.groupby("Nbin", observed=True).agg(maps=("N", "size"), nonpos_share=("nonpos", "mean"), median_log_d=("log_d", "median")).round(3).to_string())
print("non-positive by map type:"); print(A.groupby("map_type").nonpos.agg(["sum", "size", "mean"]).round(3).to_string())
np_ = A[A.nonpos]
print("non-positive maps: mean_sq quantiles", np_.mean_sq.quantile([0, .1, .5, .9, 1]).round(3).tolist())
print("non-positive maps: collections", np_.collection_id.nunique(), "largest share from one collection", np_.collection_id.value_counts().iloc[0])
print("top collections among non-positive:", np_.collection_id.value_counts().head(5).to_dict())
print("non-positive: share with baseline contrast", round(np_.baseline_contrast.mean(), 3), "vs all", round(A.baseline_contrast.mean(), 3))
# Expected sampling SD of delta2 if the map were a pure t field: sd(mean t^2) ~ sqrt(2/n_eff)/N; even with n_eff = 100, at N=10:
print("delta2 among positive maps: median", round(A.delta2[A.delta2 > 0].median(), 4))
rhs = "log_N + C(map_type) + fwhm_c + fwhm_missing + C(upload_year)"
ols = smf.ols("log_d ~ " + rhs, A).fit(cov_type="cluster", cov_kwds={"groups": A.collection_id})
q = smf.quantreg("log_d ~ " + rhs, A).fit(q=0.5)
print("OLS (no paradigm FE, clustered):", round(ols.params.log_N, 4), ols.conf_int().loc["log_N"].round(4).tolist())
print("median regression (same RHS, floor-invariant):", round(q.params.log_N, 4), q.conf_int().loc["log_N"].round(4).tolist())
q2 = smf.quantreg("log_d ~ " + rhs + " + C(task)", A).fit(q=0.5, max_iter=5000)
print("median regression with paradigm FE:", round(q2.params.log_N, 4), q2.conf_int().loc["log_N"].round(4).tolist())
# cluster bootstrap for the median regression with paradigm FE
rng = np.random.default_rng(1); cols = A.collection_id.unique(); bs = []
for _ in range(300):
    pick = rng.choice(cols, len(cols)); B = pd.concat([A[A.collection_id == c] for c in pick])
    try: bs.append(smf.quantreg("log_d ~ " + rhs + " + C(task)", B).fit(q=0.5, max_iter=5000).params.log_N)
    except Exception: pass
print("  collection-bootstrap 95% interval (", len(bs), "draws):", np.round(np.quantile(bs, [.025, .975]), 4).tolist())
# within-collection rank correlation: collections whose maps vary in N
w = [spearmanr(g.log_N, g.log_d).statistic for _, g in A.groupby("collection_id") if g.N.nunique() > 2]
print("collections with >2 distinct N:", len(w), "median within-collection rho(logN, log_d):", round(float(np.nanmedian(w)), 3))