"""Outils économétriques : hédonique semi-log, indices (dummies temporelles, ventes répétées BMN / Case–Shiller),
statistiques de ratio IAAO (COD, PRD, PRB), autocorrélation spatiale (I de Moran), matrices de poids kNN/distance."""
from __future__ import annotations
import numpy as np

# ----------------------------------------------------------------------------- MCO

def mco(X: np.ndarray, y: np.ndarray, robuste: bool = True) -> dict:
    """Moindres carrés ordinaires avec erreurs-types classiques et HC1 (White). X doit contenir la constante."""
    X = np.asarray(X, float); y = np.asarray(y, float)
    n, k = X.shape
    XtX_inv = np.linalg.pinv(X.T @ X)
    beta = XtX_inv @ X.T @ y
    res = y - X @ beta
    s2 = res @ res / (n - k)
    V_cl = s2 * XtX_inv
    meat = (X * (res[:, None] ** 2)).T @ X
    V_hc1 = XtX_inv @ meat @ XtX_inv * n / (n - k)
    tss = ((y - y.mean()) ** 2).sum()
    r2 = 1 - (res @ res) / tss
    return {"beta": beta, "se": np.sqrt(np.diag(V_cl)), "se_hc1": np.sqrt(np.diag(V_hc1)), "residus": res,
            "sigma2": s2, "r2": r2, "r2_ajuste": 1 - (1 - r2) * (n - 1) / (n - k), "n": n, "k": k, "V_hc1": V_hc1}

def indice_dummies_temps(beta_temps: np.ndarray, base: float = 100.0) -> np.ndarray:
    """Indice hédonique à partir des coefficients des indicatrices de période (période 0 = référence, coefficient 0)."""
    b = np.concatenate([[0.0], np.asarray(beta_temps, float)])
    return base * np.exp(b)

# ----------------------------------------------------------------------------- ventes répétées

def ventes_repetees_bmn(prix_1, prix_2, t_1, t_2, n_periodes: int, base: float = 100.0) -> dict:
    """Estimateur de Bailey–Muth–Nourse (1963) : ln(P2/P1) = Σ_t D_t β_t + ε, D = −1 à t1, +1 à t2, période 0 omise."""
    p1 = np.asarray(prix_1, float); p2 = np.asarray(prix_2, float)
    t1 = np.asarray(t_1, int); t2 = np.asarray(t_2, int)
    y = np.log(p2 / p1)
    D = np.zeros((len(y), n_periodes))
    for i in range(len(y)):
        D[i, t1[i]] -= 1; D[i, t2[i]] += 1
    D = D[:, 1:]
    beta, *_ = np.linalg.lstsq(D, y, rcond=None)
    res = y - D @ beta
    idx = base * np.exp(np.concatenate([[0.0], beta]))
    return {"beta": beta, "indice": idx, "residus": res, "D": D, "y": y}

def ventes_repetees_case_shiller(prix_1, prix_2, t_1, t_2, n_periodes: int, base: float = 100.0) -> dict:
    """Case–Shiller (1987, 1989) pondéré en trois étapes : (1) BMN, (2) régression de ε² sur l'intervalle,
    (3) moindres carrés pondérés par 1/√(variance ajustée)."""
    e1 = ventes_repetees_bmn(prix_1, prix_2, t_1, t_2, n_periodes, base)
    gap = (np.asarray(t_2) - np.asarray(t_1)).astype(float)
    Z = np.column_stack([np.ones_like(gap), gap])
    a, *_ = np.linalg.lstsq(Z, e1["residus"] ** 2, rcond=None)
    var_hat = np.maximum(Z @ a, 1e-8)
    w = 1 / np.sqrt(var_hat)
    beta, *_ = np.linalg.lstsq(e1["D"] * w[:, None], e1["y"] * w, rcond=None)
    idx = base * np.exp(np.concatenate([[0.0], beta]))
    return {"beta": beta, "indice": idx, "coef_variance": a, "poids": w}

# ----------------------------------------------------------------------------- ratios IAAO

def statistiques_ratio(valeurs_estimees, prix) -> dict:
    """Statistiques d'étude de ratios (IAAO, Standard on Ratio Studies) :
    ratio médian, COD = 100 × moyenne|r − médiane| / médiane, PRD = moyenne(r) / (Σ V / Σ P),
    PRB = pente de (r − médiane)/médiane sur log2 de la valeur proxy ((V/médiane + P)/2)."""
    V = np.asarray(valeurs_estimees, float); P = np.asarray(prix, float)
    r = V / P
    med = np.median(r)
    cod = 100.0 * np.mean(np.abs(r - med)) / med
    prd = np.mean(r) / (V.sum() / P.sum())
    proxy = 0.5 * (V / med + P)
    x = np.log2(proxy)
    yv = (r - med) / med
    X = np.column_stack([np.ones_like(x), x])
    b, *_ = np.linalg.lstsq(X, yv, rcond=None)
    res = yv - X @ b
    se = np.sqrt((res @ res) / (len(x) - 2) / ((x - x.mean()) ** 2).sum())
    return {"n": len(r), "ratio_median": med, "ratio_moyen": r.mean(), "COD": cod, "PRD": prd, "PRB": b[1], "PRB_se": se,
            "ratio_pondere": V.sum() / P.sum()}

# ----------------------------------------------------------------------------- spatial

def matrice_poids_knn(coords: np.ndarray, k: int = 8, standardiser: bool = True) -> np.ndarray:
    """Matrice W (n×n) des k plus proches voisins (distance euclidienne), standardisée en ligne."""
    c = np.asarray(coords, float)
    d = np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))
    np.fill_diagonal(d, np.inf)
    W = np.zeros_like(d)
    idx = np.argsort(d, axis=1)[:, :k]
    for i in range(len(c)):
        W[i, idx[i]] = 1.0
    if standardiser:
        W = W / W.sum(1, keepdims=True)
    return W

def matrice_poids_distance_inverse(coords: np.ndarray, portee: float, puissance: float = 1.0, standardiser: bool = True) -> np.ndarray:
    c = np.asarray(coords, float)
    d = np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))
    with np.errstate(divide="ignore"):
        W = np.where((d > 0) & (d <= portee), 1.0 / d ** puissance, 0.0)
    if standardiser:
        s = W.sum(1, keepdims=True); s[s == 0] = 1
        W = W / s
    return W

def moran_i(x: np.ndarray, W: np.ndarray) -> dict:
    """I de Moran et test sous hypothèse de normalité : I = (n/S0) · (z'Wz)/(z'z)."""
    z = np.asarray(x, float); z = z - z.mean()
    n = len(z); S0 = W.sum()
    I = (n / S0) * (z @ W @ z) / (z @ z)
    EI = -1.0 / (n - 1)
    S1 = 0.5 * ((W + W.T) ** 2).sum()
    S2 = ((W.sum(0) + W.sum(1)) ** 2).sum()
    VI = (n ** 2 * S1 - n * S2 + 3 * S0 ** 2) / ((n ** 2 - 1) * S0 ** 2) - EI ** 2
    zI = (I - EI) / np.sqrt(VI)
    from math import erf, sqrt
    p = 2 * (1 - 0.5 * (1 + erf(abs(zI) / sqrt(2))))
    return {"I": I, "E_I": EI, "Var_I": VI, "z": zI, "p": p}

def decalage_spatial(W: np.ndarray, x: np.ndarray) -> np.ndarray:
    return W @ np.asarray(x, float)

def sar_lag_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray, grille_rho=None) -> dict:
    """Modèle autorégressif spatial y = ρWy + Xβ + ε estimé par maximum de vraisemblance concentrée (recherche sur ρ)."""
    y = np.asarray(y, float); X = np.asarray(X, float); n = len(y)
    Wy = W @ y
    b0, *_ = np.linalg.lstsq(X, y, rcond=None); e0 = y - X @ b0
    bL, *_ = np.linalg.lstsq(X, Wy, rcond=None); eL = Wy - X @ bL
    eig = np.linalg.eigvals(W).real
    rhos = grille_rho if grille_rho is not None else np.linspace(-0.9, 0.99, 380)
    best = None
    for rho in rhos:
        e = e0 - rho * eL
        s2 = e @ e / n
        ll = np.log(np.abs(1 - rho * eig)).sum() - n / 2 * np.log(2 * np.pi * s2) - n / 2
        if best is None or ll > best[0]:
            best = (ll, rho, s2)
    ll, rho, s2 = best
    beta = b0 - rho * bL
    return {"rho": rho, "beta": beta, "sigma2": s2, "loglik": ll}

# ----------------------------------------------------------------------------- spatial, suite (chapitre 19)
# Fonctions ajoutées le 2026-10-31 pour le chapitre 19 (tests : analysis/python/tests/test_ch19.py).
# Elles n'altèrent pas les fonctions précédentes. Déterminants jacobiens calculés avec les valeurs propres
# complexes de W (|det(I − ρW)| = Π |1 − ρ λ_i|), ce qui est exact pour une matrice standardisée en ligne.

def coordonnees_km(lat, lon, lat0: float | None = None, lon0: float | None = None) -> np.ndarray:
    """Projection équirectangulaire locale : (lat, lon) → (est, nord) en kilomètres autour de (lat0, lon0)."""
    lat = np.asarray(lat, float); lon = np.asarray(lon, float)
    lat0 = float(np.mean(lat)) if lat0 is None else lat0
    lon0 = float(np.mean(lon)) if lon0 is None else lon0
    x = (lon - lon0) * 111.32 * np.cos(np.radians(lat0)); y = (lat - lat0) * 110.57
    return np.column_stack([x, y])

def _distances(c: np.ndarray) -> np.ndarray:
    c = np.asarray(c, float)
    return np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))

def _logdet(eig: np.ndarray, rho: float) -> float:
    return float(np.log(np.abs(1.0 - rho * eig)).sum())

def _intervalle_rho(eig: np.ndarray) -> tuple[float, float]:
    """Intervalle admissible (1/λ_min, 1/λ_max) sur les valeurs propres réelles de W (λ_max = 1 si standardisée)."""
    re = eig.real[np.abs(eig.imag) < 1e-8]
    lo = 1.0 / re.min() if len(re) and re.min() < -1e-8 else -0.999
    hi = 1.0 / re.max() if len(re) and re.max() > 1e-8 else 0.999
    return max(lo, -0.999) + 1e-6, min(hi, 0.999) - 1e-6

def _optimiser_concentree(f, bornes, n_grille: int = 60):
    """Maximise une vraisemblance concentrée scalaire : grille puis raffinement (section dorée, scipy)."""
    from scipy.optimize import minimize_scalar
    g = np.linspace(bornes[0], bornes[1], n_grille)
    vals = np.array([f(r) for r in g]); i = int(np.argmax(vals))
    lo = g[max(i - 1, 0)]; hi = g[min(i + 1, n_grille - 1)]
    r = minimize_scalar(lambda t: -f(t), bounds=(lo, hi), method="bounded", options={"xatol": 1e-7})
    return float(r.x), float(-r.fun)

def sar_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray, durbin: bool = False) -> dict:
    """SAR (y = ρWy + Xβ + ε) ou SDM (durbin=True : y = ρWy + Xβ + WXθ + ε) par maximum de vraisemblance concentrée.
    X contient la constante en première colonne (non décalée dans le SDM). Erreurs-types asymptotiques (matrice d'information
    analytique, Anselin 1988, ch. 6). Retourne aussi la log-vraisemblance, l'AIC et les résidus."""
    y = np.asarray(y, float); X0 = np.asarray(X, float); n = len(y)
    Xr = np.column_stack([X0, W @ X0[:, 1:]]) if durbin else X0
    k = Xr.shape[1]
    eig = np.linalg.eigvals(W)
    Wy = W @ y
    XtX_inv = np.linalg.pinv(Xr.T @ Xr)
    b0 = XtX_inv @ Xr.T @ y; e0 = y - Xr @ b0
    bL = XtX_inv @ Xr.T @ Wy; eL = Wy - Xr @ bL
    def ll_c(rho):
        e = e0 - rho * eL; s2 = e @ e / n
        return _logdet(eig, rho) - n / 2 * np.log(2 * np.pi * s2) - n / 2
    rho, ll = _optimiser_concentree(ll_c, _intervalle_rho(eig))
    beta = b0 - rho * bL; e = e0 - rho * eL; s2 = float(e @ e / n)
    # matrice d'information (β, ρ, σ²)
    A_inv = np.linalg.inv(np.eye(n) - rho * W); WA = W @ A_inv; WAXb = WA @ (Xr @ beta)
    I = np.zeros((k + 2, k + 2))
    I[:k, :k] = Xr.T @ Xr / s2
    I[:k, k] = I[k, :k] = Xr.T @ WAXb / s2
    I[k, k] = np.trace(WA.T @ WA) + np.trace(WA @ WA) + WAXb @ WAXb / s2
    I[k, k + 1] = I[k + 1, k] = np.trace(WA) / s2
    I[k + 1, k + 1] = n / (2 * s2 ** 2)
    V = np.linalg.inv(I); se = np.sqrt(np.diag(V))
    return {"rho": rho, "beta": beta, "se": se[:k], "se_rho": float(se[k]), "sigma2": s2, "loglik": ll,
            "aic": 2 * (k + 2) - 2 * ll, "residus": e, "n": n, "k": k, "X": Xr, "durbin": durbin, "ajuste": y - e}

def sem_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray) -> dict:
    """Modèle à erreur spatiale (SEM) y = Xβ + u, u = λWu + ε, par maximum de vraisemblance concentrée sur λ.
    Retourne λ, β (MCG spatiaux), σ², log-vraisemblance, AIC, résidus u (= y − Xβ) et innovations ε (= u − λWu)."""
    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape
    eig = np.linalg.eigvals(W); Wy = W @ y; WX = W @ X
    def ajuster(lam):
        ys = y - lam * Wy; Xs = X - lam * WX
        b = np.linalg.pinv(Xs.T @ Xs) @ Xs.T @ ys; eps = ys - Xs @ b
        return b, eps, float(eps @ eps / n), Xs
    def ll_c(lam):
        _, eps, s2, _ = ajuster(lam)
        return _logdet(eig, lam) - n / 2 * np.log(2 * np.pi * s2) - n / 2
    lam, ll = _optimiser_concentree(ll_c, _intervalle_rho(eig))
    beta, eps, s2, Xs = ajuster(lam)
    B_inv = np.linalg.inv(np.eye(n) - lam * W); WB = W @ B_inv
    I = np.zeros((k + 2, k + 2))
    I[:k, :k] = Xs.T @ Xs / s2
    I[k, k] = np.trace(WB.T @ WB) + np.trace(WB @ WB)
    I[k, k + 1] = I[k + 1, k] = np.trace(WB) / s2
    I[k + 1, k + 1] = n / (2 * s2 ** 2)
    V = np.linalg.inv(I); se = np.sqrt(np.diag(V))
    u = y - X @ beta
    return {"lambda": lam, "beta": beta, "se": se[:k], "se_lambda": float(se[k]), "sigma2": s2, "loglik": ll,
            "aic": 2 * (k + 2) - 2 * ll, "residus": u, "innovations": eps, "n": n, "k": k, "ajuste": X @ beta}

def tests_lm_spatial(y: np.ndarray, X: np.ndarray, W: np.ndarray) -> dict:
    """Tests du multiplicateur de Lagrange sur les résidus MCO (Anselin 1988 ; versions robustes d'Anselin, Bera, Florax et Yoon 1996) :
    LM-erreur, LM-décalage, RLM-erreur, RLM-décalage (χ² à 1 degré de liberté) et SARMA (2 d.l.)."""
    from scipy.stats import chi2
    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape
    XtX_inv = np.linalg.pinv(X.T @ X); b = XtX_inv @ X.T @ y; e = y - X @ b; s2 = e @ e / n
    M = np.eye(n) - X @ XtX_inv @ X.T
    T = np.trace((W.T + W) @ W)
    WXb = W @ (X @ b)
    RJ = (WXb @ M @ WXb) / s2 + T
    d_rho = (e @ W @ y) / s2; d_lam = (e @ W @ e) / s2
    lm_err = d_lam ** 2 / T; lm_lag = d_rho ** 2 / RJ
    rlm_lag = (d_rho - d_lam) ** 2 / (RJ - T)
    rlm_err = (d_lam - T / RJ * d_rho) ** 2 / (T * (1 - T / RJ))
    sarma = lm_err + rlm_lag
    p = lambda s, df=1: float(chi2.sf(s, df))
    return {"LM_erreur": (float(lm_err), p(lm_err)), "LM_decalage": (float(lm_lag), p(lm_lag)),
            "RLM_erreur": (float(rlm_err), p(rlm_err)), "RLM_decalage": (float(rlm_lag), p(rlm_lag)), "SARMA": (float(sarma), p(sarma, 2))}

def effets_sar(rho: float, beta: np.ndarray, W: np.ndarray, theta: np.ndarray | None = None) -> dict:
    """Effets direct, indirect et total (LeSage et Pace 2009) : S_r = (I − ρW)^{-1}(β_r I + θ_r W) ;
    direct = tr(S_r)/n, total = 1'S_r 1/n, indirect = total − direct. theta = coefficients des WX (SDM), None pour le SAR."""
    n = W.shape[0]; A_inv = np.linalg.inv(np.eye(n) - rho * W)
    beta = np.asarray(beta, float); theta = np.zeros_like(beta) if theta is None else np.asarray(theta, float)
    direct, total = [], []
    for br, tr in zip(beta, theta):
        S = A_inv @ (br * np.eye(n) + tr * W)
        direct.append(np.trace(S) / n); total.append(S.sum() / n)
    direct = np.array(direct); total = np.array(total)
    return {"direct": direct, "indirect": total - direct, "total": total}

def matrice_poids_spatio_temporelle(coords: np.ndarray, temps: np.ndarray, portee: float, fenetre: float,
                                     puissance: float = 1.0, unidirectionnel: bool = True, standardiser: bool = True) -> np.ndarray:
    """Matrice spatio-temporelle à la Dubé et Legros (2013) : w_ij = s_ij · t_ij, avec s_ij = 1/d_ij^p pour 0 < d_ij ≤ portée
    et t_ij = 1 si la vente j précède la vente i de moins de `fenetre` unités de temps (unidirectionnel : seules les ventes
    antérieures informent). temps : dates en unités numériques (mois). Standardisée en ligne (lignes vides laissées nulles)."""
    d = _distances(coords); t = np.asarray(temps, float)
    dt = t[:, None] - t[None, :]           # > 0 si j est antérieure à i
    with np.errstate(divide="ignore"):
        S = np.where((d > 0) & (d <= portee), 1.0 / d ** puissance, 0.0)
    if unidirectionnel:
        Tm = ((dt > 0) & (dt <= fenetre)).astype(float)
    else:
        Tm = ((dt != 0) & (np.abs(dt) <= fenetre)).astype(float)
    W = S * Tm
    if standardiser:
        s = W.sum(1, keepdims=True); s[s == 0] = 1.0; W = W / s
    return W

def moran_par_distance(x: np.ndarray, coords: np.ndarray, bornes) -> list[dict]:
    """Corrélogramme : I de Moran calculé avec une matrice binaire des paires dont la distance tombe dans chaque bande [b_k, b_{k+1})."""
    d = _distances(coords); out = []
    for lo, hi in zip(bornes[:-1], bornes[1:]):
        W = ((d >= lo) & (d < hi)).astype(float); np.fill_diagonal(W, 0.0)
        npaires = int(W.sum() / 2)
        if npaires < 10:
            continue
        r = moran_i(x, W); r.update({"d_min": lo, "d_max": hi, "paires": npaires}); out.append(r)
    return out

def gwr(y: np.ndarray, X: np.ndarray, coords: np.ndarray, bande: float, cibles: np.ndarray | None = None) -> dict:
    """Régression géographiquement pondérée à noyau gaussien fixe : w_ij = exp(−½ (d_ij / h)²).
    Coefficients locaux à chaque cible (par défaut les observations), valeurs ajustées, trace de la matrice chapeau,
    score de validation croisée (somme des carrés des erreurs leave-one-out) et AICc (Fotheringham, Brunsdon et Charlton 2002)."""
    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape
    c = np.asarray(coords, float); cib = c if cibles is None else np.asarray(cibles, float)
    d = np.sqrt(((cib[:, None, :] - c[None, :, :]) ** 2).sum(-1))
    Wk = np.exp(-0.5 * (d / bande) ** 2)
    betas = np.zeros((len(cib), k))
    for i in range(len(cib)):
        Xw = X * Wk[i][:, None]
        betas[i] = np.linalg.solve(X.T @ Xw + 1e-10 * np.eye(k), Xw.T @ y)
    out = {"beta": betas, "bande": bande}
    if cibles is None:
        yhat = (X * betas).sum(1)
        # matrice chapeau et validation croisée leave-one-out
        S_diag = np.zeros(n); cv = 0.0
        for i in range(n):
            Xw = X * Wk[i][:, None]
            Ainv = np.linalg.inv(X.T @ Xw + 1e-10 * np.eye(k))
            S_diag[i] = X[i] @ Ainv @ Xw[i]
            w_loo = Wk[i].copy(); w_loo[i] = 0.0
            b_loo = np.linalg.solve(X.T @ (X * w_loo[:, None]) + 1e-10 * np.eye(k), (X * w_loo[:, None]).T @ y)
            cv += (y[i] - X[i] @ b_loo) ** 2
        res = y - yhat; rss = float(res @ res); tr_S = float(S_diag.sum())
        aicc = n * np.log(rss / n) + n * np.log(2 * np.pi) + n * (n + tr_S) / (n - 2 - tr_S)
        out.update({"ajuste": yhat, "residus": res, "trace_S": tr_S, "cv": float(cv), "aicc": float(aicc),
                    "r2": 1 - rss / float(((y - y.mean()) ** 2).sum()), "sigma": float(np.sqrt(rss / (n - tr_S)))})
    return out

def gwr_bande_cv(y: np.ndarray, X: np.ndarray, coords: np.ndarray, bornes=(0.3, 15.0), n_grille: int = 25) -> dict:
    """Bande passante (km) minimisant le score de validation croisée du noyau gaussien : grille log puis section dorée."""
    from scipy.optimize import minimize_scalar
    g = np.exp(np.linspace(np.log(bornes[0]), np.log(bornes[1]), n_grille))
    cv = np.array([gwr(y, X, coords, h)["cv"] for h in g]); i = int(np.argmin(cv))
    lo = g[max(i - 1, 0)]; hi = g[min(i + 1, n_grille - 1)]
    r = minimize_scalar(lambda h: gwr(y, X, coords, h)["cv"], bounds=(lo, hi), method="bounded", options={"xatol": 1e-3})
    return {"bande": float(r.x), "cv": float(r.fun), "grille": g, "cv_grille": cv}

def variogramme_empirique(coords: np.ndarray, z: np.ndarray, bornes) -> dict:
    """Semi-variogramme empirique (Matheron) : γ(h) = Σ_{paires dans la bande} (z_i − z_j)² / (2 N(h))."""
    d = _distances(coords); z = np.asarray(z, float)
    dz2 = (z[:, None] - z[None, :]) ** 2
    iu = np.triu_indices(len(z), 1); dd = d[iu]; zz = dz2[iu]
    h, g, nb = [], [], []
    for lo, hi in zip(bornes[:-1], bornes[1:]):
        m = (dd >= lo) & (dd < hi)
        if m.sum() < 30:
            continue
        h.append(dd[m].mean()); g.append(zz[m].sum() / (2 * m.sum())); nb.append(int(m.sum()))
    return {"h": np.array(h), "gamma": np.array(g), "n": np.array(nb)}

def variogramme_exponentiel(h, pepite: float, palier_partiel: float, portee: float) -> np.ndarray:
    """γ(h) = c₀ + c₁ (1 − exp(−h/a)) pour h > 0 ; portée pratique ≈ 3a."""
    h = np.asarray(h, float)
    return np.where(h > 0, pepite + palier_partiel * (1 - np.exp(-h / portee)), 0.0)

def ajuster_variogramme(h, gamma, n_paires, portee_init: float = 2.0, portee_max: float = 10.0) -> dict:
    """Ajustement du modèle exponentiel par moindres carrés pondérés (poids N(h)/γ_modèle², Cressie 1985) ; le paramètre
    de portée a est borné par portee_max (portée pratique 3a) pour éviter les solutions dégénérées quand le palier n'est pas atteint."""
    from scipy.optimize import least_squares
    h = np.asarray(h, float); g = np.asarray(gamma, float); nb = np.asarray(n_paires, float)
    def res(p):
        gm = variogramme_exponentiel(h, *p); return np.sqrt(nb) * (g - gm) / np.maximum(gm, 1e-9)
    p0 = [0.5 * g[0], max(g[-1] - 0.5 * g[0], 1e-6), portee_init]
    r = least_squares(res, p0, bounds=([0, 0, 1e-3], [np.inf, np.inf, portee_max]))
    c0, c1, a = r.x
    return {"pepite": float(c0), "palier_partiel": float(c1), "portee": float(a), "portee_pratique": float(3 * a), "palier": float(c0 + c1)}

def krigeage_ordinaire(coords_obs: np.ndarray, z: np.ndarray, coords_cibles: np.ndarray, params: dict, voisins: int | None = 40) -> dict:
    """Krigeage ordinaire avec covariance exponentielle C(h) = c₁ exp(−h/a) (+ c₀ en h = 0), voisinage des `voisins` points
    les plus proches (None = tous). Retourne la prédiction et la variance de krigeage à chaque cible."""
    co = np.asarray(coords_obs, float); z = np.asarray(z, float); ct = np.asarray(coords_cibles, float)
    c0, c1, a = params["pepite"], params["palier_partiel"], params["portee"]
    cov = lambda h: c1 * np.exp(-h / a)
    dct = np.sqrt(((ct[:, None, :] - co[None, :, :]) ** 2).sum(-1))
    pred = np.zeros(len(ct)); var = np.zeros(len(ct))
    for i in range(len(ct)):
        idx = np.argsort(dct[i])[:voisins] if voisins else np.arange(len(co))
        sub = co[idx]; dd = _distances(sub); m = len(idx)
        K = np.zeros((m + 1, m + 1)); K[:m, :m] = cov(dd); K[np.arange(m), np.arange(m)] = c0 + c1
        K[m, :m] = 1.0; K[:m, m] = 1.0
        rhs = np.concatenate([cov(dct[i, idx]), [1.0]])
        sol = np.linalg.solve(K, rhs); w = sol[:m]; mu = sol[m]
        pred[i] = w @ z[idx]; var[i] = (c0 + c1) - w @ rhs[:m] - mu
    return {"prediction": pred, "variance": var}

def krigeage_validation_croisee(coords: np.ndarray, z: np.ndarray, params: dict, voisins: int = 40) -> np.ndarray:
    """Prédiction leave-one-out de chaque observation par krigeage ordinaire de ses voisins (hors elle-même)."""
    c = np.asarray(coords, float); z = np.asarray(z, float); d = _distances(c); out = np.zeros(len(z))
    c0, c1, a = params["pepite"], params["palier_partiel"], params["portee"]; cov = lambda h: c1 * np.exp(-h / a)
    for i in range(len(z)):
        idx = np.argsort(d[i])[1:voisins + 1]; m = len(idx); dd = d[np.ix_(idx, idx)]
        K = np.zeros((m + 1, m + 1)); K[:m, :m] = cov(dd); K[np.arange(m), np.arange(m)] = c0 + c1; K[m, :m] = 1; K[:m, m] = 1
        rhs = np.concatenate([cov(d[i, idx]), [1.0]]); sol = np.linalg.solve(K, rhs); out[i] = sol[:m] @ z[idx]
    return out

def sem_esperance_conditionnelle(lam: float, W_aug: np.ndarray, u_obs: np.ndarray, idx_cibles) -> np.ndarray:
    """Espérance conditionnelle des perturbations u aux cibles sachant les u observés, dans un SEM de matrice augmentée W_aug
    (observations puis cibles) : avec Q = (I − λW)'(I − λW) (précision à σ² près), E[u_c | u_o] = −Q_cc^{-1} Q_co u_o."""
    n = W_aug.shape[0]; B = np.eye(n) - lam * W_aug; Q = B.T @ B
    idx_c = np.asarray(idx_cibles); idx_o = np.setdiff1d(np.arange(n), idx_c)
    Qcc = Q[np.ix_(idx_c, idx_c)]; Qco = Q[np.ix_(idx_c, idx_o)]
    return -np.linalg.solve(Qcc, Qco @ np.asarray(u_obs, float))