"""Validation diagnostics (Phase 3): Monte Carlo simulation of the listing process against the analytic solution,
residual checks, bounds, reproducibility across seeds and sample sizes. Writes results/logs/validation.json and
results/processed/val_*.csv.

Run: python -m src.diagnostics
"""
from __future__ import annotations

import json
import sys
import time
from pathlib import Path

import numpy as np
import pandas as pd

from .model import load_params, Params
from .solver import solve_equilibrium, solve_sellers

ROOT = Path(__file__).resolve().parents[1]
PRO = ROOT / "results" / "processed"; LOG = ROOT / "results" / "logs"


def simulate_listings(eq, p: Params, n_listings: int, seed: int) -> pd.DataFrame:
    """Simulate n_listings sellers: type c ~ f, meetings ~ Poisson(Lam_s), withdrawal ~ Exp(delta), match w ~ G,
    acceptance w >= v*(c), price = (1-beta)V - beta U + beta w. Returns one row per listing."""
    rng = np.random.default_rng(seed)
    types = rng.choice(len(p.c), size=n_listings, p=p.f)
    Lam = eq.Lam_s
    out = np.empty((n_listings, 5))
    for i, ti in enumerate(types):
        vstar, V = eq.sellers.vstar[ti], eq.sellers.V[ti]
        t = 0.0; sold = 0; price = np.nan; n_meet = 0
        t_withdraw = rng.exponential(1.0 / p.delta) if p.delta > 0 else np.inf
        while True:
            t += rng.exponential(1.0 / Lam)
            if t > t_withdraw:
                t = t_withdraw; break
            n_meet += 1
            w = float(p.G.sample(rng, 1)[0])
            if w >= vstar:
                sold = 1; price = (1 - p.beta) * V - p.beta * eq.U + p.beta * w; break
        out[i] = (ti, t, sold, price, n_meet)
    return pd.DataFrame(out, columns=["type", "duration", "sold", "price", "n_meetings"]).astype({"type": int, "sold": int})


def main():
    p = load_params(ROOT / "config" / "baseline.yaml")
    eq = solve_equilibrium(p)
    report = {"config_hash": p.config_hash(), "timestamp": time.strftime("%Y-%m-%dT%H:%M:%S")}
    sb = eq.sellers
    report["residuals"] = {"seller_max_abs": float(np.abs(sb.resid).max()), "free_entry": float(eq.free_entry_resid),
                           "stock_eq_max_abs": float(np.abs(p.s * p.f - (sb.eta + p.delta) * eq.n).max()),
                           "omega_sum": float(eq.omega.sum()), "N_sum": float(eq.N.sum())}
    report["bounds"] = {"acc_in_[0,1]": bool(np.all((sb.acc >= 0) & (sb.acc <= 1))), "eta_nonneg": bool(np.all(sb.eta >= 0)),
                        "prices_in_support": bool(np.all((sb.pbar >= p.G.lo) & (sb.pbar <= p.G.hi))),
                        "V_monotone_dec": bool(np.all(np.diff(sb.V) < 0)), "eta_monotone_inc": bool(np.all(np.diff(sb.eta) > 0)),
                        "participation": bool(np.all(sb.V >= p.w0))}
    # Monte Carlo cross-check
    rows = []
    for n in (2000, 20000, 100000):
        for seed in (1, 2, 3):
            t0 = time.time(); df = simulate_listings(eq, p, n, seed)
            sold = df[df.sold == 1]
            an_sold_share = float((p.f * sb.eta / (sb.eta + p.delta)).sum())
            an_dur = float((p.f / (sb.eta + p.delta)).sum())
            an_price = float((eq.omega * sb.pbar).sum())
            an_meet = float((p.f * eq.Lam_s / (sb.eta + p.delta)).sum())
            rows.append({"n": n, "seed": seed, "sold_share_mc": df.sold.mean(), "sold_share_an": an_sold_share,
                         "mean_duration_mc": df.duration.mean(), "mean_duration_an": an_dur,
                         "mean_price_mc": sold.price.mean(), "mean_price_an": an_price,
                         "sd_price_mc": sold.price.std(), "sd_price_an": eq.moments(p)["price_sd_total"],
                         "mean_meetings_mc": df.n_meetings.mean(), "mean_meetings_an": an_meet,
                         "mean_tom_sold_mc": sold.duration.mean(), "mean_tom_sold_an": eq.moments(p)["mean_tom_sold"],
                         "seconds": time.time() - t0})
    mc = pd.DataFrame(rows); mc.to_csv(PRO / "val_montecarlo.csv", index=False)
    big = mc[mc.n == 100000]
    report["montecarlo"] = {k: {"max_abs_rel_err": float((abs(big[f"{k}_mc"] - big[f"{k}_an"]) / big[f"{k}_an"]).max())}
                            for k in ("sold_share", "mean_duration", "mean_price", "sd_price", "mean_meetings", "mean_tom_sold")}
    # per-type check with the largest sample
    df = simulate_listings(eq, p, 100000, 7)
    g = df.groupby("type").agg(sold_share=("sold", "mean"), mean_dur=("duration", "mean"), n=("sold", "size"),
                               mean_price=("price", "mean"))
    g["sold_share_an"] = sb.eta / (sb.eta + p.delta); g["mean_dur_an"] = 1 / (sb.eta + p.delta); g["mean_price_an"] = sb.pbar
    g["c"] = p.c; g.to_csv(PRO / "val_montecarlo_types.csv")
    report["montecarlo_types_max_rel_err_price"] = float((abs(g.mean_price - g.mean_price_an) / g.mean_price_an).max())
    # reproducibility: same seed -> identical results
    d1 = simulate_listings(eq, p, 5000, 11); d2 = simulate_listings(eq, p, 5000, 11)
    report["reproducible_same_seed"] = bool(d1.equals(d2))
    # duration-dependence check (Lemma 3): exposure-based piecewise-constant hazard vs analytic mixture hazard
    bins = np.linspace(0, 0.6, 13); emp = []; ana = []
    mu = sb.eta + p.delta
    for a, b in zip(bins[:-1], bins[1:]):
        exposure = (np.clip(df.duration, a, b) - a).sum(); exits = ((df.duration >= a) & (df.duration < b)).sum()
        emp.append(exits / exposure if exposure > 0 else np.nan)
        mid = 0.5 * (a + b); wts = p.f * np.exp(-mu * mid); ana.append(float((wts * mu).sum() / wts.sum()))
    report["hazard_by_duration"] = {"bin_upper_years": bins[1:].tolist(), "empirical": [float(h) for h in emp], "analytic_mixture": ana,
                                    "max_abs_rel_err": float(np.nanmax(np.abs(np.array(emp) - np.array(ana)) / np.array(ana)))}
    report["hazard_decreasing_analytic"] = bool(np.all(np.diff(ana) < 0))
    report["hazard_decreasing_empirical_first3_vs_last3"] = bool(np.nanmean(emp[:3]) > np.nanmean(emp[-3:]))
    (LOG / "validation.json").write_text(json.dumps(report, indent=2))
    print(json.dumps(report, indent=2))


if __name__ == "__main__":
    main()