# 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
"""Transformations compositionnelles (parts de vote ↔ espace log-ratio).

Les intentions de vote sont compositionnelles (somme = 100 %). On modélise
l'état latent en espace « additive log-ratio » (alr) : y_k = log(p_k / p_ref),
ce qui garantit des parts positives sommant à 1 après transformation inverse
et induit la corrélation négative naturelle entre partis.
"""
from __future__ import annotations

import numpy as np

FLOOR = 0.0025  # plancher de part (0,25 %) pour stabilité numérique


def close(p: np.ndarray) -> np.ndarray:
    """Renormalise en composition stricte (somme 1), avec plancher."""
    p = np.asarray(p, dtype=float)
    p = np.clip(p, FLOOR, None)
    return p / p.sum(axis=-1, keepdims=True)


def alr(p: np.ndarray, ref: int = 0) -> np.ndarray:
    """shares (…, K) → log-ratios (…, K-1), référence = composante `ref`."""
    p = close(p)
    others = [i for i in range(p.shape[-1]) if i != ref]
    return np.log(p[..., others]) - np.log(p[..., ref:ref + 1])


def inv_alr(y: np.ndarray, ref: int = 0) -> np.ndarray:
    """log-ratios (…, K-1) → shares (…, K) via softmax."""
    y = np.asarray(y, dtype=float)
    K = y.shape[-1] + 1
    full = np.zeros(y.shape[:-1] + (K,), dtype=float)
    others = [i for i in range(K) if i != ref]
    full[..., others] = y
    full -= full.max(axis=-1, keepdims=True)
    e = np.exp(full)
    return e / e.sum(axis=-1, keepdims=True)


def alr_obs_cov(p: np.ndarray, n_eff: float, excess_pp: float = 0.0,
                ref: int = 0) -> np.ndarray:
    """Covariance d'observation en espace alr (méthode delta sur la multinomiale).

    Cov(y_i, y_j) ≈ (1/n)·(δ_ij / p_i + 1 / p_ref)  pour i,j ≠ ref,
    plus une variance excédentaire (erreur non-échantillonnale) diagonale.
    """
    p = close(p)
    others = [i for i in range(p.shape[-1]) if i != ref]
    po = p[others]
    pr = float(p[ref])
    n_eff = max(float(n_eff), 50.0)
    cov = (np.diag(1.0 / po) + 1.0 / pr) / n_eff
    if excess_pp > 0:
        e = excess_pp / 100.0
        cov = cov + np.diag((e / po) ** 2 + (e / pr) ** 2)
    return cov