#!/usr/bin/env python3
"""
cross_domain_significance.py — Statistical significance test for
RRC cross-domain 1/n residual signatures.

Tests whether observed deviations from predicted 1/n scaling or
1/7 meta-solid threshold are statistically significant at ≥6σ.
"""

import json, math
from pathlib import Path
from datetime import datetime, timezone

ROOT = Path(__file__).resolve().parents[1]

def gaussian_p_value(z: float) -> float:
    """Two-tailed p-value from z-score via complementary error function."""
    return math.erfc(abs(z) / math.sqrt(2))

def sigma_from_p(p: float) -> float:
    """Convert p-value to sigma level (two-tailed) via inverse erfc approximation."""
    if p <= 0:
        return float('inf')
    if p >= 1:
        return 0.0
    t = math.sqrt(-2 * math.log(p if p <= 0.5 else 1 - p))
    c = (2.515517, 0.802853, 0.010328)
    d = (1.432788, 0.189269, 0.001308)
    return t - (c[0] + c[1]*t + c[2]*t**2) / (1 + d[0]*t + d[1]*t**2 + d[2]*t**3)

def main():
    sig_path = ROOT / "signatures" / "cross_domain_signatures.json"
    sigs = json.loads(sig_path.read_text())
    entries = sigs["signatures"]

    phases = {}
    for e in entries:
        phases.setdefault(e.get("phase", 0), []).append(e)

    results = []
    for phase_num in sorted(phases):
        group = phases[phase_num]
        domain = group[0].get("domain", "Unknown")
        deviations = [abs(e["deviation"]) for e in group if "deviation" in e]
        residuals = [e["residual_mhz"] for e in group if "residual_mhz" in e]

        values = deviations or residuals
        n = len(values)

        if n < 2:
            results.append({
                "phase": phase_num, "domain": domain,
                "n": n, "status": "insufficient_data"
            })
            continue

        mean_val = sum(values) / n
        var = sum((v - mean_val)**2 for v in values) / (n - 1) if n > 1 else 0
        se = math.sqrt(var / n)
        z = mean_val / se if se > 0 else 0.0
        p = gaussian_p_value(z)
        sigma = sigma_from_p(p)

        results.append({
            "phase": phase_num,
            "domain": domain,
            "n": n,
            "mean_deviation": round(mean_val, 6),
            "std_err": round(se, 6),
            "z_score": round(z, 4),
            "sigma_level": round(sigma, 2) if sigma < 1e6 else "inf",
            "passes_6sigma": bool(sigma >= 6.0),
        })

    summary = {
        "schema": "cross_domain_significance_v1",
        "generated_at": datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"),
        "source": "cross_domain_signatures.json",
        "total_entries": len(entries),
        "phases": results,
    }

    out_path = ROOT / "signatures" / "cross_domain_significance.json"
    out_path.write_text(json.dumps(summary, indent=2) + "\n")
    print(f"Wrote {out_path}")
    for r in results:
        s = r.get("sigma_level")
        if s is not None and s != "inf":
            tag = "✅" if r["passes_6sigma"] else "❌"
            print(f"  Phase {r['phase']} ({r['domain']}): σ={s} {tag}")
        else:
            tag = "✅" if r.get("status") == "insufficient_data" else "❌"
            print(f"  Phase {r['phase']} ({r['domain']}): σ={s} {tag}")

if __name__ == "__main__":
    main()