# QC Élection Forecast — Plateforme de prévision électorale du Québec 2026
# Auteur : Simon-Pierre Boucher
# Contact : contact@spboucher.ai
# https://www.qc-election.com
"""Scoring de prévisions électorales (v3 « 127 » — validation §20/§38).

Toutes les métriques standard au même endroit, réutilisées par le backtest,
le replay historique, l'ablation et le tournoi de modèles :
MAE/RMSE du vote, couverture d'intervalles multi-niveaux, Brier multiclasse,
log score, exactitude des vainqueurs, sharpness et biais signé.
"""
from __future__ import annotations

import numpy as np


def vote_errors(forecast: dict[str, float], actual: dict[str, float],
                parties: list[str]) -> dict:
    errs = [forecast[p] - actual[p] for p in parties if p in forecast and p in actual]
    return {"mae_pp": round(float(np.mean(np.abs(errs))), 3),
            "rmse_pp": round(float(np.sqrt(np.mean(np.square(errs)))), 3),
            "bias_pp": {p: round(forecast[p] - actual[p], 2)
                        for p in parties if p in forecast and p in actual}}


def interval_coverage(summary: dict, actual: dict[str, float],
                      parties: list[str]) -> dict:
    """Couverture des intervalles 50/80/95 % (résumé {party:{p05..p95,lo95,hi95}})."""
    out = {}
    for level, lo_k, hi_k in [("50", "p25", "p75"), ("90", "p05", "p95"),
                              ("95", "lo95", "hi95")]:
        hits = [1.0 if summary[p].get(lo_k) is not None
                and summary[p][lo_k] <= actual[p] <= summary[p][hi_k] else 0.0
                for p in parties if p in summary and p in actual
                and summary[p].get(lo_k) is not None]
        if hits:
            out[f"coverage_{level}"] = round(float(np.mean(hits)), 3)
    return out


def winner_scores(district_probs: list[dict], actual_winners: dict[str, str],
                  parties: list[str]) -> dict:
    """district_probs: [{district, win_probs{party:p}, favorite}]."""
    briers, lls, correct, n = [], [], 0, 0
    for d in district_probs:
        actual = actual_winners.get(d["district"])
        if actual is None:
            continue
        n += 1
        wp = d["win_probs"]
        briers.append(sum((wp.get(p, 0.0) - (1.0 if p == actual else 0.0)) ** 2
                          for p in parties))
        lls.append(-np.log(max(wp.get(actual, 0.0), 1e-6)))
        if d.get("favorite") == actual:
            correct += 1
    if not n:
        return {}
    return {"brier": round(float(np.mean(briers)), 4),
            "log_loss": round(float(np.mean(lls)), 4),
            "winner_accuracy": round(correct / n, 4), "n_districts": n}


def sharpness(summary: dict, parties: list[str]) -> float | None:
    """Largeur moyenne des IC 95 % (pp) — une prévision calibrée ET étroite
    domine; à calibration égale, plus petit = mieux."""
    widths = [summary[p]["hi95"] - summary[p]["lo95"]
              for p in parties if p in summary and "hi95" in summary[p]]
    return round(float(np.mean(widths)), 2) if widths else None


def calibration_bins(probs: np.ndarray, outcomes: np.ndarray,
                     width: float = 0.1, min_n: int = 5) -> list[dict]:
    """Diagramme de fiabilité (bacs de largeur `width`)."""
    bins = []
    for lo in np.arange(0.0, 1.0, width):
        hi = lo + width
        m = (probs >= lo) & ((probs < hi) if hi < 1.0 else (probs <= 1.0))
        if m.sum() >= min_n:
            bins.append({"bin": f"{int(lo * 100)}–{int(hi * 100)} %",
                         "predicted": round(float(probs[m].mean()), 3),
                         "observed": round(float(outcomes[m].mean()), 3),
                         "n": int(m.sum())})
    return bins


def crps_from_draws(draws: np.ndarray, actual: np.ndarray) -> float:
    """CRPS moyen par parti à partir d'échantillons du posterior (en pp).
    CRPS = E|X − y| − ½·E|X − X'| (estimateur par paires décalées)."""
    n = len(draws)
    term1 = np.abs(draws - actual[None, :]).mean(axis=0)
    perm = np.roll(np.arange(n), n // 2)
    term2 = 0.5 * np.abs(draws - draws[perm]).mean(axis=0)
    return round(float((term1 - term2).mean()), 3)


def log_score_alr(x: np.ndarray, P: np.ndarray, actual_alr: np.ndarray) -> float:
    """Log-densité prédictive de la composition réelle sous le posterior alr
    (gaussienne multivariée) — plus haut = mieux."""
    Mx = len(x)
    diff = actual_alr - x
    Pi = P + 1e-10 * np.eye(Mx)
    sign, logdet = np.linalg.slogdet(Pi)
    return round(float(-0.5 * (Mx * np.log(2 * np.pi) + logdet
                               + diff @ np.linalg.inv(Pi) @ diff)), 3)


def score_snapshot(forecast_summary: dict, districts: list[dict],
                   actual_national: dict[str, float],
                   actual_winners: dict[str, str],
                   parties: list[str]) -> dict:
    """Tableau de bord complet d'un snapshot vs les résultats réels."""
    fc_mean = {p: forecast_summary[p]["mean"] for p in parties
               if p in forecast_summary}
    return {**vote_errors(fc_mean, actual_national, parties),
            **interval_coverage(forecast_summary, actual_national, parties),
            **winner_scores(districts, actual_winners, parties),
            "sharpness_pp": sharpness(forecast_summary, parties)}