"""F5-conf-03: counterexample to 'sum of normals is normal' (proposition A),
and baseline check for the independent case (proposition B).

Run from workspace root:  ~/wp29/venv/bin/python out/analysis.py
Reads only data/; writes out/results.json and out/counterexample.csv deterministically.
"""
import json
import os

import numpy as np
import pandas as pd
from scipy import stats

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
DATA = os.path.join(ROOT, "data")
OUT = os.path.join(ROOT, "out")

props = pd.read_csv(os.path.join(DATA, "propositions.csv"))
design = pd.read_csv(os.path.join(DATA, "design.csv"))

data_issues = []
if sorted(props["id"].tolist()) != ["A", "B"]:
    data_issues.append({"issue": "unexpected proposition ids", "handling": "reported"})
dmap = dict(zip(design["parameter"], design["value"]))
seed = int(dmap["seed"])
n = int(dmap["n"])

# Single generator, fixed draw order: X, then S, then baseline independent Y.
rng = np.random.default_rng(seed)
x = rng.standard_normal(n)
s = rng.choice(np.array([-1.0, 1.0]), size=n)  # Rademacher sign, independent of X
y = s * x                                          # Y = S X ~ N(0,1) by symmetry
z = x + y                                          # = (1+S) X : 0 w.p. 1/2, else 2X

def ks_vs_fitted_normal(v):
    m, sd = float(np.mean(v)), float(np.std(v, ddof=1))
    return float(stats.kstest(v, "norm", args=(m, sd)).statistic)

ks_sum = ks_vs_fitted_normal(z)
ks_x = ks_vs_fitted_normal(x)
ks_y = ks_vs_fitted_normal(y)
frac_sum_zero = float(np.mean(z == 0.0))

# Baseline (proposition B): independent Y' ~ N(0,1), same evaluation rule.
y_ind = rng.standard_normal(n)
ks_sum_indep = ks_vs_fitted_normal(x + y_ind)

pd.DataFrame({"x": x, "y": y}).to_csv(os.path.join(OUT, "counterexample.csv"),
                                      index=False, float_format="%.17g")
# Re-read so that summary statistics describe exactly the written file.
cx = pd.read_csv(os.path.join(OUT, "counterexample.csv"))
xv, yv = cx["x"].to_numpy(), cx["y"].to_numpy()
zv = xv + yv

cv = {
    "construction": ("X ~ N(0,1); S independent of X with P(S=1)=P(S=-1)=1/2; Y = S*X. "
                     "Then Y ~ N(0,1) by symmetry of N(0,1), and X+Y = (1+S)X equals 0 with "
                     "probability 1/2 and 2X ~ N(0,4) otherwise (non-degenerate, not normal)."),
    "n": int(len(cx)),
    "seed": seed,
    "mean_x": float(np.mean(xv)), "sd_x": float(np.std(xv, ddof=1)),
    "mean_y": float(np.mean(yv)), "sd_y": float(np.std(yv, ddof=1)),
    "mean_sum": float(np.mean(zv)), "sd_sum": float(np.std(zv, ddof=1)),
    "ks_sum_vs_normal": ks_vs_fitted_normal(zv),
}

results = {
    "task_id": "F5-conf-03",
    "statement_true": False,
    "counterexample_provided": True,
    "counterexample_values": cv,
    "statement_b_true": True,
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": data_issues,
    "n_input_rows": {"propositions.csv": int(len(props)), "design.csv": int(len(design))},
    "n_analysed": int(len(cx)),
    "random_seed": seed,
    # supplementary diagnostics
    "ks_x_vs_normal": ks_vs_fitted_normal(xv),
    "ks_y_vs_normal": ks_vs_fitted_normal(yv),
    "frac_sum_exactly_zero": float(np.mean(zv == 0.0)),
    "ks_critical_value_alpha_0_01": float(1.628 / np.sqrt(n)),
    "baseline_ks_sum_independent": ks_sum_indep,
}
with open(os.path.join(OUT, "results.json"), "w") as f:
    json.dump(results, f, indent=2)
print(json.dumps({k: results[k] for k in ["statement_true", "statement_b_true"]}))
print("ks_sum_vs_normal =", cv["ks_sum_vs_normal"], "| baseline independent =", ks_sum_indep)