#!/usr/bin/env python3
"""Figures du chapitre 19 — modèles spatiaux et temporels (cas 1, Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).
Toutes les cartes sont tracées sur les coordonnées synthétiques (km est/nord autour du centroïde des ventes), sans fond de carte.
Sorties (assets/figures/ch19/) :
  ch19-fig-carte-residus.pdf      résidus MCO sans indicatrices de secteur vs modèle de référence
  ch19-fig-correlogramme.pdf      I de Moran par bande de distance
  ch19-fig-gwr-coefficients.pdf   coefficients locaux de la GWR (localisation, superficie, âge)
  ch19-fig-gwr-bande.pdf          score de validation croisée selon la bande passante
  ch19-fig-variogramme.pdf        variogramme empirique et modèle exponentiel
  ch19-fig-krigeage.pdf           surface krigée des résidus
  ch19-fig-sensibilite-W.pdf      I, ρ et λ selon la matrice de poids
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch19_figures.py"""
import sys, pathlib
ROOT = pathlib.Path(__file__).resolve().parents[3]
sys.path.insert(0, str(ROOT / "analysis/python")); sys.path.insert(0, str(ROOT / "analysis/python/cases"))
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
from traite.style import setup, palette, savefig, LARGEUR, DEMI
from traite import hedonique as h
import ch19_spatial as c19

SECTEURS = ["Hull", "Gatineau", "Aylmer", "Buckingham", "Masson-Angers"]

def _axes_km(ax, xy):
    ax.set_aspect("equal"); ax.set_xlim(xy[:, 0].min() - 1.5, xy[:, 0].max() + 1.5); ax.set_ylim(xy[:, 1].min() - 1.5, xy[:, 1].max() + 1.5)
    ax.set_xlabel("Est (km, coordonnées synthétiques)"); ax.grid(True, lw=0.4)

def _etiquettes_secteurs(ax, v, xy, dy=0.0):
    for s in SECTEURS:
        m = (v["secteur"] == s).values
        cx, cy = np.median(xy[m, 0]), np.median(xy[m, 1])
        ax.text(cx, cy + dy, s, fontsize=7, color=palette.ink2, ha="center", va="center",
                bbox=dict(boxstyle="round,pad=0.15", fc=palette.surface, ec="none", alpha=0.75))

def fig_carte_residus(res):
    b = res["base"]; v, xy, cs = b["ventes"], b["xy"], b["cs"]
    series = [("Sans indicatrices de secteur (I = %.2f)" % res["moran"].loc[2, "I_knn"], b["mco_r"]["residus"]),
              ("Modèle de référence (I = %.2f)" % res["moran"].loc[3, "I_knn"], b["modele"]["residus"])]
    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.3), sharey=True)
    vmax = 0.25; cmap = palette.cmap_div()
    for ax, (titre, e) in zip(axes, series):
        sc = ax.scatter(xy[:, 0], xy[:, 1], c=np.clip(e, -vmax, vmax), cmap=cmap, vmin=-vmax, vmax=vmax, s=9, lw=0)
        ax.scatter(cs[:, 0], cs[:, 1], marker="*", s=110, color=palette.ink, zorder=5)
        _axes_km(ax, xy); ax.set_title(titre, fontsize=8.5); _etiquettes_secteurs(ax, v, xy, dy=4.5)
    axes[0].set_ylabel("Nord (km)")
    cb = fig.colorbar(sc, ax=axes, orientation="horizontal", fraction=0.05, pad=0.16, aspect=40)
    cb.set_label("Résidu du modèle semi-logarithmique (ln P observé − ln P ajusté)", fontsize=8)
    cb.ax.tick_params(labelsize=7.5)
    savefig(fig, "ch19/ch19-fig-carte-residus.pdf")

def fig_correlogramme(res):
    cg = res["correlogramme"]
    fig, ax = plt.subplots(figsize=(LARGEUR, 3.0))
    for i, (cle, nom) in enumerate([("sans_secteur", "Résidus sans indicatrices de secteur"), ("reference", "Résidus du modèle de référence")]):
        pts = cg[cle]; x = np.array([(p["d_min"] + p["d_max"]) / 2 for p in pts]); I = np.array([p["I"] for p in pts])
        sd = np.sqrt(np.array([p["Var_I"] for p in pts]))
        ax.errorbar(x + (i - 0.5) * 0.12, I, yerr=1.96 * sd, fmt="o-", color=palette.cat[i], ms=4.5, lw=1.4, capsize=2, label=nom)
    ax.axhline(cg["reference"][0]["E_I"], color=palette.ink2, lw=0.9, ls="--")
    ax.text(14.6, cg["reference"][0]["E_I"] + 0.012, "E[I] = −1/(n−1)", fontsize=7, color=palette.ink2, ha="right")
    ax.set_xlabel("Distance entre paires de ventes (km, centre de la bande)"); ax.set_ylabel("I de Moran de la bande")
    ax.legend(loc="upper right", fontsize=7.5)
    savefig(fig, "ch19/ch19-fig-correlogramme.pdf")

def fig_gwr_coefficients(res):
    b = res["base"]; v, xy, cs = b["ventes"], b["xy"], b["cs"]; g = res["gwr"]["gwr"]; cols = list(b["Xr"].columns)
    j = cols.index("dist_centre"); dist = b["Xr"]["dist_centre"].values
    sect = v["secteur"].values
    loc = g["beta"][:, 0] + g["beta"][:, j] * dist - np.median(g["beta"][sect == "Hull", 0])
    panneaux = [("Effet de localisation local\n(écart à Hull, en log)", loc, palette.cmap_div(), -0.35, 0.35),
                ("Coefficient local de ln(superficie)\n(élasticité)", g["beta"][:, cols.index("ln_S")], palette.cmap_seq(), None, None),
                ("Coefficient local de l'âge\n(× 100, % par année)", 100 * g["beta"][:, cols.index("age")], palette.cmap_seq(), None, None)]
    fig, axes = plt.subplots(1, 3, figsize=(LARGEUR, 2.75))
    for ax, (titre, z, cmap, vmin, vmax) in zip(axes, panneaux):
        sc = ax.scatter(xy[:, 0], xy[:, 1], c=z, cmap=cmap, vmin=vmin, vmax=vmax, s=6, lw=0)
        ax.scatter(cs[:, 0], cs[:, 1], marker="*", s=70, color=palette.ink, zorder=5)
        _axes_km(ax, xy); ax.set_title(titre, fontsize=7.5); ax.tick_params(labelsize=6.5); ax.set_xlabel("Est (km)", fontsize=7)
        cb = fig.colorbar(sc, ax=ax, orientation="horizontal", fraction=0.06, pad=0.22, aspect=28); cb.ax.tick_params(labelsize=6.5)
    axes[0].set_ylabel("Nord (km)", fontsize=7)
    savefig(fig, "ch19/ch19-fig-gwr-coefficients.pdf")

def fig_gwr_bande(res):
    sel = res["gwr"]["selection"]
    fig, ax = plt.subplots(figsize=(DEMI, 2.5))
    ax.plot(sel["grille"], sel["cv_grille"], marker="o", ms=3.5, color=palette.cat[0])
    ax.axvline(sel["bande"], color=palette.ink2, ls="--", lw=0.9)
    ax.text(sel["bande"] * 1.08, ax.get_ylim()[1] * 0.985, "h* = %.2f km" % sel["bande"], fontsize=7, color=palette.ink2, va="top")
    ax.set_xscale("log"); ax.set_xlabel("Bande passante h (km, échelle log)"); ax.set_ylabel("Score de validation croisée\n(somme des carrés, log)")
    ax.set_xticks([0.3, 0.5, 1, 2, 5, 10]); ax.get_xaxis().set_major_formatter(FuncFormatter(lambda x, p: ("%g" % x)))
    savefig(fig, "ch19/ch19-fig-gwr-bande.pdf")

def fig_variogramme(res):
    vk = res["variogramme"]; vg, vr, p = vk["variogramme"], vk["variogramme_reference"], vk["parametres"]
    fig, ax = plt.subplots(figsize=(LARGEUR, 3.1))
    taille = 12 + 60 * vg["n"] / vg["n"].max()
    ax.scatter(vg["h"], vg["gamma"], s=taille, color=palette.cat[0], alpha=0.85, lw=0, label="Résidus sans indicatrices de secteur (taille ∝ nombre de paires)", zorder=3)
    hh = np.linspace(0.05, 16, 200)
    ax.plot(hh, h.variogramme_exponentiel(hh, p["pepite"], p["palier_partiel"], p["portee"]), color=palette.cat[0], lw=1.6,
            label="Modèle exponentiel ajusté (pépite %.4f, palier %.4f, portée pratique %.0f km)" % (p["pepite"], p["palier"], p["portee_pratique"]))
    ax.scatter(vr["h"], vr["gamma"], s=14, marker="s", color=palette.cat[1], lw=0, label="Résidus du modèle de référence", zorder=3)
    ax.axhline(vk["sigma2_reference"], color=palette.cat[1], lw=0.9, ls=":")
    ax.axhline(vk["variance_residus"], color=palette.cat[0], lw=0.9, ls=":")
    ax.text(15.8, vk["variance_residus"] + 0.0004, "variance des résidus sans secteur", fontsize=7, color=palette.ink2, ha="right")
    ax.text(15.8, vk["sigma2_reference"] - 0.0012, "σ² du modèle de référence", fontsize=7, color=palette.ink2, ha="right")
    ax.set_xlabel("Distance h entre paires de ventes (km)"); ax.set_ylabel("Semi-variance γ(h) (log²)")
    ax.set_ylim(0, 0.02); ax.legend(loc="upper left", fontsize=7)
    savefig(fig, "ch19/ch19-fig-variogramme.pdf")

def fig_krigeage(res):
    b = res["base"]; v, xy, cs = b["ventes"], b["xy"], b["cs"]; gr = res["variogramme"]["grille"]
    GX, GY = np.meshgrid(gr["x"], gr["y"]); pts = np.column_stack([GX.ravel(), GY.ravel()])
    from scipy.spatial import cKDTree
    dmin = cKDTree(xy).query(pts)[0].reshape(GX.shape)
    Z = np.ma.masked_where(dmin > 2.5, gr["z"])
    fig, ax = plt.subplots(figsize=(LARGEUR, 3.6))
    vmax = 0.25
    cf = ax.contourf(GX, GY, np.clip(Z, -vmax, vmax), levels=np.linspace(-vmax, vmax, 21), cmap=palette.cmap_div(), vmin=-vmax, vmax=vmax)
    ax.scatter(xy[:, 0], xy[:, 1], s=3, color=palette.ink, alpha=0.45, lw=0)
    ax.scatter(cs[:, 0], cs[:, 1], marker="*", s=120, color=palette.ink, zorder=5)
    ax.annotate("sujet : résidu krigé %.3f" % res["variogramme"]["sujet"]["prediction"][0], (cs[0, 0], cs[0, 1]), xytext=(-70, -22), textcoords="offset points",
                fontsize=7, color=palette.ink, arrowprops=dict(arrowstyle="-", color=palette.ink2, lw=0.7))
    _axes_km(ax, xy); ax.set_ylabel("Nord (km)"); _etiquettes_secteurs(ax, v, xy, dy=5.0)
    cb = fig.colorbar(cf, ax=ax, fraction=0.035, pad=0.02); cb.set_label("Résidu krigé (log)", fontsize=8); cb.ax.tick_params(labelsize=7.5)
    savefig(fig, "ch19/ch19-fig-krigeage.pdf")

def fig_sensibilite(res):
    s = res["sensibilite"]; k = s[s["famille"].str.startswith("k")]
    fig, ax = plt.subplots(figsize=(LARGEUR, 2.9))
    ax.errorbar(k["parametre"], k["rho"], yerr=1.96 * k["se_rho"], fmt="o-", color=palette.cat[0], ms=4.5, capsize=2, label="ρ du SAR (± 1,96 é.-t.)")
    ax.errorbar(k["parametre"], k["lambda"], yerr=1.96 * k["se_lambda"], fmt="s-", color=palette.cat[1], ms=4.5, capsize=2, label="λ du SEM (± 1,96 é.-t.)")
    ax.plot(k["parametre"], k["I"], "^-", color=palette.cat[2], ms=4.5, label="I de Moran des résidus MCO")
    ax.set_xlabel("Nombre k de plus proches voisins dans W"); ax.set_ylabel("Coefficient"); ax.set_xticks(k["parametre"].astype(int))
    ax.set_ylim(0, 0.95); ax.legend(loc="upper left", fontsize=7.5)
    savefig(fig, "ch19/ch19-fig-sensibilite-W.pdf")

def main():
    setup()
    res = c19.main(verbose=False)
    for f in [fig_carte_residus, fig_correlogramme, fig_gwr_coefficients, fig_gwr_bande, fig_variogramme, fig_krigeage, fig_sensibilite]:
        f(res); print("OK", f.__name__)

if __name__ == "__main__":
    main()