"""Outils géospatiaux génériques du traité (chapitre 38 ; réutilisables aux chapitres 43 et 48).
Fonctions pures NumPy, sans dépendance à geopandas : distances (haversine, équirectangulaire locale), facteur d'échelle
d'une projection transverse de Mercator, point dans polygone (lancer de rayon), interpolation par l'inverse de la
distance, grille hexagonale, classification des cartes choroplèthes (intervalles égaux, quantiles, Fisher–Jenks, GVF),
accessibilité gravitaire, conversion d'une probabilité sur 25 ans en probabilité annuelle, effet minimal détectable.
Tests : analysis/python/tests/test_ch38.py. Ajouté le 2026-10-31 ; n'altère aucune fonction existante."""
from __future__ import annotations
import numpy as np

R_TERRE_KM = 6371.0088          # rayon moyen de la Terre (IUGG, rayon volumétrique), km
MTM_K0 = 0.9999                 # facteur d'échelle au méridien central des zones MTM du Québec (EPSG 32187-32189)
UTM_K0 = 0.9996


def haversine_km(lat1, lon1, lat2, lon2, rayon: float = R_TERRE_KM) -> np.ndarray:
    """Distance orthodromique (formule du haversine) sur une sphère de rayon `rayon`, en km. Entrées en degrés."""
    p1, p2 = np.radians(np.asarray(lat1, float)), np.radians(np.asarray(lat2, float))
    dphi = p2 - p1; dlam = np.radians(np.asarray(lon2, float) - np.asarray(lon1, float))
    a = np.sin(dphi / 2) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(dlam / 2) ** 2
    return 2 * rayon * np.arcsin(np.sqrt(np.clip(a, 0.0, 1.0)))


def equirectangulaire_km(lat, lon, lat0: float, lon0: float) -> np.ndarray:
    """Projection équirectangulaire locale autour de (lat0, lon0) : x = R cos(lat0) Δλ, y = R Δφ (km). Même rayon que
    haversine_km, de sorte que les deux mesures coïncident au premier ordre autour du point de référence."""
    lat = np.radians(np.asarray(lat, float)); lon = np.radians(np.asarray(lon, float))
    x = R_TERRE_KM * np.cos(np.radians(lat0)) * (lon - np.radians(lon0)); y = R_TERRE_KM * (lat - np.radians(lat0))
    return np.column_stack([x, y])


def facteur_echelle_tm(lat, lon, lon0: float, k0: float = MTM_K0) -> np.ndarray:
    """Facteur d'échelle ponctuel d'une projection transverse de Mercator (approximation sphérique au second ordre) :
    k ≈ k0 [1 + ½ (Δλ cos φ)²], Δλ en radians. Une distance au sol de 1 km est représentée par k km sur la carte."""
    dl = np.radians(np.asarray(lon, float) - lon0); phi = np.radians(np.asarray(lat, float))
    return k0 * (1.0 + 0.5 * (dl * np.cos(phi)) ** 2)


def point_dans_polygone(points: np.ndarray, polygone: np.ndarray) -> np.ndarray:
    """Test point-dans-polygone par lancer de rayon (parité du nombre de croisements d'une demi-droite horizontale).
    points : (m, 2) ; polygone : (n, 2) sommets ordonnés (fermeture implicite). Retourne un tableau booléen (m,)."""
    P = np.asarray(points, float); V = np.asarray(polygone, float)
    x, y = P[:, 0], P[:, 1]; n = len(V); dedans = np.zeros(len(P), bool)
    for i in range(n):
        x1, y1 = V[i]; x2, y2 = V[(i + 1) % n]
        croise = (y1 > y) != (y2 > y)
        with np.errstate(divide="ignore", invalid="ignore"):
            xi = x1 + (y - y1) * (x2 - x1) / (y2 - y1)
        dedans ^= croise & (x < xi)
    return dedans


def idw(coords_obs: np.ndarray, z: np.ndarray, coords_cibles: np.ndarray, puissance: float = 2.0, voisins: int | None = None) -> np.ndarray:
    """Interpolation par l'inverse de la distance : ẑ(s0) = Σ w_i z_i / Σ w_i, w_i = d_i^(−p) (valeur exacte si d_i = 0)."""
    co = np.asarray(coords_obs, float); z = np.asarray(z, float); ct = np.asarray(coords_cibles, float)
    d = np.sqrt(((ct[:, None, :] - co[None, :, :]) ** 2).sum(-1)); out = np.zeros(len(ct))
    for i in range(len(ct)):
        idx = np.argsort(d[i])[:voisins] if voisins else np.arange(len(co))
        di = d[i, idx]
        if di.min() == 0:
            out[i] = z[idx][di.argmin()]; continue
        w = di ** (-puissance); out[i] = (w @ z[idx]) / w.sum()
    return out


def grille_hexagonale(xmin: float, ymin: float, xmax: float, ymax: float, rayon: float) -> list[np.ndarray]:
    """Pavage hexagonal (hexagones « pointe en haut » de rayon circonscrit `rayon`) couvrant la boîte englobante.
    Retourne la liste des sommets (6, 2) de chaque cellule ; aire d'une cellule = (3√3/2) rayon²."""
    w = np.sqrt(3) * rayon; hstep = 1.5 * rayon; cells = []
    ang = np.radians(np.arange(30, 390, 60))
    j = 0; y = ymin
    while y <= ymax + rayon:
        off = (w / 2) if (j % 2) else 0.0; x = xmin + off
        while x <= xmax + w:
            cells.append(np.column_stack([x + rayon * np.cos(ang), y + rayon * np.sin(ang)])); x += w
        y += hstep; j += 1
    return cells


def bornes_intervalles_egaux(z, k: int) -> np.ndarray:
    z = np.asarray(z, float); return np.linspace(z.min(), z.max(), k + 1)


def bornes_quantiles(z, k: int) -> np.ndarray:
    z = np.asarray(z, float); return np.quantile(z, np.linspace(0, 1, k + 1))


def bornes_jenks(z, k: int) -> np.ndarray:
    """Classification de Fisher–Jenks (programmation dynamique exacte) : minimise la somme intra-classe des carrés des écarts.
    Retourne k + 1 bornes (min, ruptures, max)."""
    z = np.sort(np.asarray(z, float)); n = len(z)
    if k >= n:
        return np.concatenate([[z[0]], z[1:], [z[-1]]])
    cs = np.concatenate([[0.0], np.cumsum(z)]); cs2 = np.concatenate([[0.0], np.cumsum(z ** 2)])
    def ssd(i, j):  # somme des carrés des écarts de z[i:j]
        m = j - i; s = cs[j] - cs[i]; return (cs2[j] - cs2[i]) - s * s / m
    INF = np.inf
    cout = np.full((k + 1, n + 1), INF); arg = np.zeros((k + 1, n + 1), int); cout[0, 0] = 0.0
    for c in range(1, k + 1):
        for j in range(c, n + 1):
            best, bi = INF, -1
            for i in range(c - 1, j):
                v = cout[c - 1, i] + ssd(i, j)
                if v < best:
                    best, bi = v, i
            cout[c, j] = best; arg[c, j] = bi
    bornes = [z[-1]]; j = n
    for c in range(k, 0, -1):
        i = arg[c, j]; bornes.append(z[i]); j = i
    return np.array(bornes[::-1])


def gvf(z, bornes) -> float:
    """Goodness of variance fit = 1 − SDCM/SDAM (somme des carrés intra-classe / somme totale des carrés)."""
    z = np.asarray(z, float); b = np.asarray(bornes, float)
    cl = np.clip(np.searchsorted(b, z, side="right") - 1, 0, len(b) - 2)
    sdcm = sum(((z[cl == c] - z[cl == c].mean()) ** 2).sum() for c in range(len(b) - 1) if (cl == c).any())
    return 1.0 - sdcm / ((z - z.mean()) ** 2).sum()


def classer(z, bornes) -> np.ndarray:
    """Indice de classe (0 … k−1) de chaque valeur selon des bornes croissantes ; la dernière borne est incluse."""
    z = np.asarray(z, float); b = np.asarray(bornes, float)
    return np.clip(np.searchsorted(b, z, side="right") - 1, 0, len(b) - 2)


def accessibilite_gravitaire(coords: np.ndarray, coords_poles: np.ndarray, masses, friction: float, forme: str = "exp") -> np.ndarray:
    """Indice de Hansen A_i = Σ_k E_k f(d_ik), f(d) = exp(−λ d) (forme « exp ») ou d^(−λ) (forme « puissance », d en km, d ≥ 0,1)."""
    c = np.asarray(coords, float); p = np.asarray(coords_poles, float); E = np.asarray(masses, float)
    d = np.sqrt(((c[:, None, :] - p[None, :, :]) ** 2).sum(-1))
    f = np.exp(-friction * d) if forme == "exp" else np.maximum(d, 0.1) ** (-friction)
    return f @ E


def probabilite_annuelle(p_horizon: float, horizon: int = 25) -> float:
    """Probabilité annuelle d'au moins une occurrence équivalente à une probabilité p sur `horizon` années (événements indépendants) :
    p_1 = 1 − (1 − p)^(1/horizon)."""
    return 1.0 - (1.0 - p_horizon) ** (1.0 / horizon)


def effet_minimal_detectable(se: float, alpha: float = 0.05, puissance: float = 0.80) -> float:
    """Effet minimal détectable (bilatéral) : (z_{1−α/2} + z_{puissance}) · se."""
    from scipy.stats import norm
    return float((norm.ppf(1 - alpha / 2) + norm.ppf(puissance)) * se)