"""
MULTIPLE-COMPARISONS CORRECTION (Section 4.5, Section 5.2, Section 5.4).

This script takes the eight raw McNemar's exact-test p-values already
computed and reported by real_data_pipeline.py (four aggregation rules,
pooled-model peer groups) and same_model_robustness_check.py combined with
real_data_pipeline.py's run_full_statistical_analysis() (four aggregation
rules, same-model-only peer groups), and applies the Holm-Bonferroni
step-down procedure to control the family-wise error rate across all
eight comparisons jointly, since interpreting each p-value against an
uncorrected alpha=0.05 threshold in isolation inflates the true false-
positive rate for the family of tests as a whole. A Benjamini-Hochberg
false-discovery-rate correction is also reported for comparison; for
this particular set of eight p-values, both procedures yield identical
significance decisions, though this need not be true in general.

USAGE: after running real_data_pipeline.py and same_model_robustness_check.py
(which will have populated real_data_summary.json and
same_model_statistical_analysis.json with the raw McNemar p-values used
here), run: python3 multiple_comparisons_correction.py
"""
import json
from statsmodels.stats.multitest import multipletests


def load_raw_pvalues():
    with open("real_data_summary.json") as f:
        pooled = json.load(f)
    with open("same_model_statistical_analysis.json") as f:
        same_model = json.load(f)

    pooled_mcnemar = pooled["statistical_analysis"]["mcnemar_solo_vs_aggregation"]
    same_model_mcnemar = same_model["mcnemar_solo_vs_aggregation"]

    labels, pvals = [], []
    for key, val in pooled_mcnemar.items():
        labels.append(f"pooled: {key}")
        pvals.append(val["p_value"])
    for key, val in same_model_mcnemar.items():
        labels.append(f"same-model: {key}")
        pvals.append(val["p_value"])
    return labels, pvals


def main():
    labels, pvals = load_raw_pvalues()
    print(f"Loaded {len(pvals)} raw McNemar's p-values from real_data_summary.json "
          f"and same_model_statistical_analysis.json.\n")

    reject_holm, p_holm, _, _ = multipletests(pvals, alpha=0.05, method="holm")
    reject_bh, p_bh, _, _ = multipletests(pvals, alpha=0.05, method="fdr_bh")

    print(f"{'Comparison':32s} {'raw_p':>10s} {'holm_p':>10s} {'sig(holm)':>10s} "
          f"{'bh_p':>10s} {'sig(bh)':>10s}")
    results = []
    for lab, p, ph, rh, pb, rb in zip(labels, pvals, p_holm, reject_holm, p_bh, reject_bh):
        print(f"{lab:32s} {p:10.6f} {ph:10.6f} {str(rh):>10s} {pb:10.6f} {str(rb):>10s}")
        results.append({
            "comparison": lab, "raw_p": p,
            "holm_adjusted_p": float(ph), "significant_holm": bool(rh),
            "bh_adjusted_p": float(pb), "significant_bh": bool(rb),
        })

    agree = all(r["significant_holm"] == r["significant_bh"] for r in results)
    print(f"\nHolm and Benjamini-Hochberg agree on all {len(results)} significance "
          f"decisions: {agree}")

    with open("multiple_comparisons_correction_summary.json", "w") as f:
        json.dump({"n_tests": len(results), "results": results,
                   "holm_bh_agree": agree}, f, indent=2)
    print("\nSaved: multiple_comparisons_correction_summary.json")


if __name__ == "__main__":
    main()