import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from matplotlib.patches import Patch

import rasterio
from rasterio.plot import show
from rasterio.transform import from_bounds
from rasterio.warp import reproject, Resampling

from scipy.interpolate import griddata

import geopandas as gpd
import os

# Adapte les chemins selon l'emplacement de tes fichiers
BAND_FILES = [
    'rasters/Landsat4_1990_194028_b1.tif',
    'rasters/Landsat4_1990_194028_b2.tif',
    'rasters/Landsat4_1990_194028_b3.tif',
    'rasters/Landsat4_1990_194028_b4.tif',
    'rasters/Landsat4_1990_194028_b5.tif',
    'rasters/Landsat4_1990_194028_b6.tif',
    'rasters/Landsat4_1990_194028_b7.tif',
]
LANDCOVER_FILE = 'rasters/MapTicino1990_8classes.tif'
MNT_FILE       = 'rasters/MNT25_Ticino.tif'
DHM25_PATH     = 'dhm25_p/dhm25_p.shp'   # semis de points altimétriques (Partie 6)

# Lecture des métadonnées de la bande 4 (proche infrarouge)
with rasterio.open(BAND_FILES[3]) as src:   # index 3 = bande 4
    print("=== Métadonnées — Bande 4 (NIR) ===")
    print(f"CRS          : {src.crs}")
    print(f"Dimensions   : {src.height} lignes × {src.width} colonnes")
    print(f"Résolution   : {src.res[0]:.1f} m × {src.res[1]:.1f} m")
    print(f"Emprise      : {src.bounds}")
    print(f"Type données : {src.dtypes[0]}")
    print(f"NoData       : {src.nodata}")
    meta_ref = src.meta.copy()   # métadonnées de référence pour les sorties

# Lecture des 7 bandes et empilement
bands = []
for path in BAND_FILES:
    with rasterio.open(path) as src:
        bands.append(src.read(1).astype(np.float32))

landsat = np.stack(bands)   # shape : (7, height, width)
print(f"Image multibandes : {landsat.shape}")
print(f"  7 bandes × {landsat.shape[1]} lignes × {landsat.shape[2]} colonnes")

# Sauvegarde GeoTIFF multibandes (équivalent du Merge QGIS)
os.makedirs('rasters/output', exist_ok=True)
meta_multi = meta_ref.copy()
meta_multi.update({'count': 7, 'dtype': 'float32'})

with rasterio.open('rasters/output/Landsat4_1990_multibande.tif', 'w', **meta_multi) as dst:
    dst.write(landsat)

print("\nGeoTIFF multibandes sauvegardé : rasters/output/Landsat4_1990_multibande.tif")

def normalize_band(arr, pmin=2, pmax=98):
    """Étirement de contraste : coupe les percentiles extrêmes et normalise en [0, 1]."""
    lo, hi = np.nanpercentile(arr, [pmin, pmax])
    return np.clip((arr - lo) / (hi - lo + 1e-6), 0, 1)

def make_rgb(r_idx, g_idx, b_idx):
    """Construit une image RGB normalisée depuis les indices de bandes (1-based)."""
    r = normalize_band(landsat[r_idx - 1])
    g = normalize_band(landsat[g_idx - 1])
    b = normalize_band(landsat[b_idx - 1])
    return np.dstack([r, g, b])

compositions = [
    ('Vraies couleurs (3-2-1)',             3, 2, 1),
    ('Fausses couleurs végétation (4-3-2)', 4, 3, 2),
    ('Agriculture (5-4-1)',                 5, 4, 1),
]

fig, axes = plt.subplots(1, 3, figsize=(18, 6))
for ax, (title, r, g, b) in zip(axes, compositions):
    ax.imshow(make_rgb(r, g, b), origin='upper')
    ax.set_title(title, fontsize=10)
    ax.axis('off')

plt.suptitle('Compositions colorées — Landsat 4, Tessin 1990',
             fontsize=13, fontweight='bold')
plt.tight_layout()
plt.show()

rouge = landsat[2].astype(np.float32)   # B3 = Rouge
nir   = landsat[3].astype(np.float32)   # B4 = NIR

with np.errstate(divide='ignore', invalid='ignore'):
    ndvi = np.where(
        (nir + rouge) == 0,
        np.nan,
        (nir - rouge) / (nir + rouge)
    )

print(f"NDVI — min={np.nanmin(ndvi):.3f}, max={np.nanmax(ndvi):.3f}, "
      f"moy={np.nanmean(ndvi):.3f}")

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

im = axes[0].imshow(ndvi, cmap='RdYlGn', vmin=-0.3, vmax=0.8, origin='upper')
plt.colorbar(im, ax=axes[0], label='NDVI', shrink=0.8)
axes[0].set_title('NDVI — Tessin 1990', fontsize=11)
axes[0].axis('off')

axes[1].hist(ndvi[~np.isnan(ndvi)].ravel(), bins=80,
             color='#388e3c', alpha=0.8, edgecolor='white')
axes[1].axvline(0,   color='navy',      linestyle='--', label='NDVI = 0')
axes[1].axvline(0.3, color='darkgreen', linestyle='--', label='NDVI = 0.3')
axes[1].set_xlabel('NDVI')
axes[1].set_ylabel('Pixels')
axes[1].set_title('Distribution des valeurs NDVI')
axes[1].legend()

plt.tight_layout()
plt.show()

# Sauvegarde
meta_ndvi = meta_ref.copy()
meta_ndvi.update({'dtype': 'float32', 'count': 1})
with rasterio.open('rasters/output/NDVI_Ticino1990.tif', 'w', **meta_ndvi) as dst:
    dst.write(ndvi.astype(np.float32), 1)
print("NDVI sauvegardé : rasters/output/NDVI_Ticino1990.tif")

with rasterio.open(LANDCOVER_FILE) as src:
    landcover = src.read(1)
    meta_lc   = src.meta.copy()

# Le MNT est en EPSG:21781 (LV03) ; on le reprojette sur la grille de la carte d'occupation du sol
with rasterio.open(MNT_FILE) as src_mnt:
    mnt = np.empty((meta_lc['height'], meta_lc['width']), dtype=np.float32)
    reproject(
        source=src_mnt.read(1).astype(np.float32),
        destination=mnt,
        src_transform=src_mnt.transform,
        src_crs=src_mnt.crs,
        dst_transform=meta_lc['transform'],
        dst_crs=meta_lc['crs'],
        resampling=Resampling.bilinear,
    )

print(f"Occupation du sol : {landcover.shape}")
print(f"MNT reprojeté     : {mnt.shape}  (même grille que l'occupation du sol)")

CLASSES = {1: 'Forêt', 2: 'Prés', 3: 'Eau', 4: 'Neige',
           5: 'Sols nus', 6: 'Aires urbaines', 7: 'Nuages', 8: 'Ombres'}
COLORS  = ['#2e7d32', '#aed581', '#1565c0', '#e3f2fd',
           '#bcaaa4', '#e53935', '#b0bec5', '#616161']

print("\nDistribution des classes (avant reclassification) :")
n_total = landcover.size
for val, nom in CLASSES.items():
    n = (landcover == val).sum()
    print(f"  Classe {val} ({nom:15s}) : {n:8d} pixels  ({100*n/n_total:.1f} %)")

cmap_lc = mcolors.ListedColormap(COLORS)
norm_lc = mcolors.BoundaryNorm(range(1, 10), cmap_lc.N)
patches = [Patch(color=COLORS[i], label=f"{i+1} — {CLASSES[i+1]}") for i in range(8)]

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

axes[0].imshow(landcover, cmap=cmap_lc, norm=norm_lc, origin='upper')
axes[0].legend(handles=patches, loc='lower right', fontsize=8, title='Classes')
axes[0].set_title('Occupation du sol — Tessin 1990', fontsize=11)
axes[0].axis('off')

im = axes[1].imshow(mnt, cmap='terrain', origin='upper')
plt.colorbar(im, ax=axes[1], label='Altitude (m)', shrink=0.8)
axes[1].set_title('MNT — Tessin', fontsize=11)
axes[1].axis('off')

plt.tight_layout()
plt.show()

# Reclassification : Aires urbaines (6) au-dessus de 1400 m → Sols nus (5)
# Équivalent QGIS : if( (MNT > 1400) AND (classes = 6), 5, classes )
condition = (mnt > 1400) & (landcover == 6)

landcover_reclass = np.where(condition, 5, landcover)

n_corriges = condition.sum()
print(f"Pixels reclassifiés (Urbain→Sol nu, altitude >1400 m) : {n_corriges}")
print()
print("Distribution après reclassification :")
for val, nom in CLASSES.items():
    n_avant = (landcover == val).sum()
    n_apres = (landcover_reclass == val).sum()
    diff = n_apres - n_avant
    marker = f"  ({diff:+d})" if diff != 0 else ""
    print(f"  Classe {val} ({nom:15s}) : {n_apres:8d} pixels{marker}")

fig, axes = plt.subplots(1, 2, figsize=(16, 7))
for ax, data, title in [
    (axes[0], landcover,         'Original'),
    (axes[1], landcover_reclass, 'Après reclassification (Urbain→Sol nu >1400 m)'),
]:
    ax.imshow(data, cmap=cmap_lc, norm=norm_lc, origin='upper')
    ax.legend(handles=patches, loc='lower right', fontsize=7, title='Classes')
    ax.set_title(title, fontsize=11)
    ax.axis('off')

plt.suptitle("Reclassification de l'occupation du sol", fontsize=13, fontweight='bold')
plt.tight_layout()
plt.show()

# Sauvegarde
with rasterio.open('rasters/output/ReClass_MapTicino1990_8classes.tif', 'w', **meta_lc) as dst:
    dst.write(landcover_reclass.astype(meta_lc['dtype']), 1)
print("Reclassifié sauvegardé : rasters/output/ReClass_MapTicino1990_8classes.tif")

pts = gpd.read_file(DHM25_PATH)
print(f"Points chargés : {len(pts)}")
print(f"Colonnes       : {pts.columns.tolist()}")
print(f"CRS            : {pts.crs}")
print(pts.head(3))

# L'altitude est encodée dans la coordonnée Z de la géométrie POINT Z
xy_pts = np.column_stack([pts.geometry.x, pts.geometry.y])
z_pts  = pts.geometry.z.values

# Grille cible (résolution 900 m)
RES = 900
bounds = pts.total_bounds   # [xmin, ymin, xmax, ymax]
xi = np.arange(bounds[0], bounds[2], RES)
yi = np.arange(bounds[1], bounds[3], RES)
XX, YY = np.meshgrid(xi, yi)

# Interpolation linéaire (équivalent IDW local)
mnt_idw = griddata(xy_pts, z_pts, (XX, YY), method='linear')
mnt_idw = np.flipud(mnt_idw)   # convention raster : ligne 0 en haut

print(f"MNT interpolé : {mnt_idw.shape}  —  résolution {RES} m")
print(f"Altitude : min={np.nanmin(mnt_idw):.0f} m, max={np.nanmax(mnt_idw):.0f} m")

fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(mnt_idw, cmap='terrain', origin='upper')
plt.colorbar(im, ax=ax, label='Altitude (m)', shrink=0.8)
ax.set_title(f'MNT interpolé IDW — résolution {RES} m\n(depuis points DHM25 Swisstopo)', fontsize=11)
ax.axis('off')
plt.tight_layout()
plt.show()

# Sauvegarde
idw_transform = from_bounds(bounds[0], bounds[1], bounds[2], bounds[3],
                             mnt_idw.shape[1], mnt_idw.shape[0])
with rasterio.open(
    'rasters/output/MNT_IDW_900m.tif', 'w',
    driver='GTiff', height=mnt_idw.shape[0], width=mnt_idw.shape[1],
    count=1, dtype='float32', crs=pts.crs, transform=idw_transform,
) as dst:
    dst.write(mnt_idw.astype(np.float32), 1)
print("MNT IDW sauvegardé : rasters/output/MNT_IDW_900m.tif")

def hillshade(dem, res=25, azimuth=315, altitude=45):
    """Calcule l'ombrage depuis un MNT NumPy.
    azimuth  : direction de la lumière (degrés depuis le Nord, sens horaire)
    altitude : angle d'élévation de la source lumineuse (degrés)
    """
    az_rad  = np.radians(360 - azimuth + 90)
    alt_rad = np.radians(altitude)

    dy, dx = np.gradient(dem, res, res)
    slope  = np.arctan(np.sqrt(dx**2 + dy**2))
    aspect = np.arctan2(-dy, dx)

    hs = (np.sin(alt_rad) * np.cos(slope) +
          np.cos(alt_rad) * np.sin(slope) * np.cos(az_rad - aspect))
    return np.clip(hs, 0, 1)

# Le MNT reprojeté a la résolution de la carte d'occupation du sol (~31 m)
mnt_res = meta_lc['transform'].a   # taille du pixel en x
hs = hillshade(mnt, res=mnt_res)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

axes[0].imshow(mnt, cmap='terrain', origin='upper')
axes[0].set_title('MNT brut', fontsize=11)
axes[0].axis('off')

axes[1].imshow(hs, cmap='gray', vmin=0, vmax=1, origin='upper')
axes[1].set_title('Ombrage (Hillshade)\nazimuth=315°, élévation=45°', fontsize=11)
axes[1].axis('off')

plt.suptitle('Relief ombré — Tessin', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.show()

# Sauvegarde
meta_hs = meta_ref.copy()
meta_hs.update({'dtype': 'float32', 'count': 1})
with rasterio.open('rasters/output/Hillshade_Ticino.tif', 'w', **meta_hs) as dst:
    dst.write(hs.astype(np.float32), 1)
print("Hillshade sauvegardé : rasters/output/Hillshade_Ticino.tif")

from sklearn.ensemble import RandomForestClassifier
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, ConfusionMatrixDisplay
import xgboost as xgb

LC_NAMES  = {1: 'Forêt', 2: 'Prés', 3: 'Eau', 4: 'Neige', 5: 'Sols nus', 6: 'Urbain'}
LC_COLORS = {1: '#2e7d32', 2: '#aed581', 3: '#1565c0', 4: '#e3f2fd', 5: '#bcaaa4', 6: '#e53935'}
BAND_NAMES = ['B1-Bleu', 'B2-Vert', 'B3-Rouge', 'B4-PIR', 'B5-SWIR1', 'B6-Therm', 'B7-SWIR2']

# 1. Rééchantillonner les 7 bandes Landsat sur la grille de l'occupation du sol
print("Rééchantillonnage des 7 bandes Landsat sur la grille de l'occupation du sol...")
bands_on_lc = []
for i, path in enumerate(BAND_FILES):
    band = np.empty((meta_lc['height'], meta_lc['width']), dtype=np.float32)
    with rasterio.open(path) as src:
        reproject(
            source=src.read(1).astype(np.float32),
            destination=band,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=meta_lc['transform'],
            dst_crs=meta_lc['crs'],
            resampling=Resampling.bilinear,
        )
    bands_on_lc.append(band)
    print(f"  {BAND_NAMES[i]:12s} rééchantillonné")

ls_on_lc = np.stack(bands_on_lc, axis=0)   # (7, H, W)
print(f"\nImage rééchantillonnée : {ls_on_lc.shape}  (bandes × lignes × colonnes)")

# 2. Masque : classes valides 1–6, toutes bandes finies et > 0, MNT positif
valid_clf = (
    np.isin(landcover, [1, 2, 3, 4, 5, 6]) &
    np.all(np.isfinite(ls_on_lc) & (ls_on_lc > 0), axis=0) &
    np.isfinite(mnt) & (mnt > 0)
)
print(f"Pixels valides : {valid_clf.sum():,} / {valid_clf.size:,} ({100*valid_clf.mean():.1f} %)")

# 3. Sous-échantillonnage aléatoire : 20 000 pixels
rng = np.random.default_rng(42)
idx_valid_clf = np.argwhere(valid_clf)
chosen = rng.choice(len(idx_valid_clf), 20_000, replace=False)
rows_clf, cols_clf = idx_valid_clf[chosen, 0], idx_valid_clf[chosen, 1]

# 4. Features : 7 bandes + MNT (8 features par pixel)
X_clf = np.column_stack([
    ls_on_lc[:, rows_clf, cols_clf].T,   # (20000, 7)
    mnt[rows_clf, cols_clf],              # (20000,)
])
y_clf = landcover[rows_clf, cols_clf]    # labels 1–6

feat_names_clf = BAND_NAMES + ['MNT']
print(f"\nFeatures ({X_clf.shape[1]}) : {feat_names_clf}")
print(f"Taille du jeu de données : {len(X_clf):,} pixels\n")
print("Distribution des classes :")
for v in sorted(LC_NAMES):
    n = (y_clf == v).sum()
    print(f"  {v} – {LC_NAMES[v]:15s} : {n:5d} pixels ({100*n/len(y_clf):.1f} %)")

from sklearn.neural_network import MLPClassifier
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

# --- Séparation 50 / 15 / 15 / 20 (entraînement / validation / calibration / test) ---
frac_train, frac_val, frac_cal, frac_test = 0.50, 0.15, 0.15, 0.20
idx_all = np.arange(len(X_clf))

idx_tr, idx_rest  = train_test_split(idx_all,  test_size=1 - frac_train, random_state=42)
idx_val, idx_rest2 = train_test_split(idx_rest, test_size=(frac_cal + frac_test) / (1 - frac_train), random_state=0)
idx_cal, idx_te    = train_test_split(idx_rest2, test_size=frac_test / (frac_cal + frac_test), random_state=1)

X_tr, X_val, X_cal, X_te = X_clf[idx_tr], X_clf[idx_val], X_clf[idx_cal], X_clf[idx_te]
y_tr, y_val, y_cal, y_te = y_clf[idx_tr], y_clf[idx_val], y_clf[idx_cal], y_clf[idx_te]
rows_te, cols_te = rows_clf[idx_te], cols_clf[idx_te]

print(f"Entraînement : {len(X_tr):,}   Validation : {len(X_val):,}   "
      f"Calibration : {len(X_cal):,}   Test : {len(X_te):,}")

# --- Random Forest ---
rf = RandomForestClassifier(n_estimators=200, max_depth=15, random_state=42, n_jobs=-1).fit(X_tr, y_tr)

# --- XGBoost (labels obligatoirement 0-indexés en version ≥ 3) ---
le = LabelEncoder()
y_tr_enc = le.fit_transform(y_tr)                          # 1-6 → 0-5
xgb_clf = xgb.XGBClassifier(
    n_estimators=200, max_depth=6, learning_rate=0.1,
    random_state=42, eval_metric='mlogloss', verbosity=0
).fit(X_tr, y_tr_enc)

# --- MLP : Pipeline avec standardisation (voir §7.2) ---
mlp = make_pipeline(
    StandardScaler(),
    MLPClassifier(hidden_layer_sizes=(32, 16), max_iter=500, random_state=42)
).fit(X_tr, y_tr)

# models[name] = (modèle, encodeur ou None) — l'encodeur ramène XGBoost à l'échelle 1-6
models = {'Random Forest': (rf, None), 'XGBoost': (xgb_clf, le), 'MLP': (mlp, None)}

def predict_labels(name, X):
    model, enc = models[name]
    pred = model.predict(X)
    return enc.inverse_transform(pred) if enc is not None else pred

# --- Sélection de modèle sur le jeu de VALIDATION (jamais utilisé pour l'entraînement) ---
print("\nSélection de modèle — précision sur le jeu de validation :")
val_acc = {}
for name in models:
    val_acc[name] = accuracy_score(y_val, predict_labels(name, X_val))
    print(f"  {name:15s} : {val_acc[name]:.3f}")

best_name = max(val_acc, key=val_acc.get)
print(f"\n→ Modèle retenu (validation la plus élevée) : {best_name}")

# --- Évaluation finale sur le jeu de TEST, pour les trois modèles (comparaison) ---
y_pred_rf  = predict_labels('Random Forest', X_te)
y_pred_xgb = predict_labels('XGBoost', X_te)
y_pred_mlp = predict_labels('MLP', X_te)
acc_rf  = accuracy_score(y_te, y_pred_rf)
acc_xgb = accuracy_score(y_te, y_pred_xgb)
acc_mlp = accuracy_score(y_te, y_pred_mlp)

class_labels = [LC_NAMES[v] for v in sorted(LC_NAMES)]

fig, axes = plt.subplots(1, 3, figsize=(19, 5))
for ax, preds, title in [
    (axes[0], y_pred_rf,  f'Random Forest  (test = {acc_rf:.3f})'),
    (axes[1], y_pred_xgb, f'XGBoost        (test = {acc_xgb:.3f})'),
    (axes[2], y_pred_mlp, f'MLP            (test = {acc_mlp:.3f})'),
]:
    ConfusionMatrixDisplay.from_predictions(
        y_te, preds, display_labels=class_labels,
        cmap='Blues', ax=ax, colorbar=False
    )
    ax.set_title(title, fontsize=11)
    ax.set_xlabel('Classe prédite')
    ax.set_ylabel('Classe réelle')

plt.suptitle("Matrices de confusion — classification de l'occupation du sol\n"
             "(lignes = classe réelle, colonnes = classe prédite)", fontsize=12, fontweight='bold')
plt.tight_layout()
plt.show()

# --- Importances des variables ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for ax, imp, title, color in [
    (axes[0], rf.feature_importances_,      f'Random Forest  (acc = {acc_rf:.3f})',  '#42a5f5'),
    (axes[1], xgb_clf.feature_importances_, f'XGBoost        (acc = {acc_xgb:.3f})', '#ef6c00'),
]:
    order = np.argsort(imp)
    ax.barh(
        [feat_names_clf[i] for i in order], imp[order],
        color=color, alpha=0.85, edgecolor='white'
    )
    ax.set_xlabel('Importance relative')
    ax.set_title(title, fontsize=11)
    ax.grid(True, alpha=0.3, axis='x')

plt.suptitle('Importances des variables prédictives\n'
             '(quelle bande aide le plus à distinguer les classes ?)',
             fontsize=12, fontweight='bold')
plt.tight_layout()
plt.show()

print("Classement RF — de la plus à la moins importante :")
for i in np.argsort(rf.feature_importances_)[::-1]:
    print(f"  {feat_names_clf[i]:12s} : {rf.feature_importances_[i]:.4f}")

from math import ceil

alpha = 0.10   # couverture cible : 90 %
best_model, best_enc = models[best_name]

# 1. Scores de non-conformité sur le jeu de calibration
#    Méthode LAC (Least Ambiguous Classifier) : score = 1 – P(vraie classe | x)
#    Un score proche de 0 signifie que le modèle est confiant et correct.
#    Un score proche de 1 signifie que le modèle a attribué très peu de probabilité
#    à la vraie classe (erreur ou incertitude élevée).
proba_cal = best_model.predict_proba(X_cal)                                  # (n_cal, 6)
# best_enc.inverse_transform(best_model.classes_) remet XGBoost à l'échelle 1-6 ;
# pour RF/MLP, .classes_ est déjà dans cette échelle (best_enc = None).
classes_best = best_enc.inverse_transform(best_model.classes_) if best_enc is not None else best_model.classes_
true_idx_cal = np.array([np.where(classes_best == yc)[0][0] for yc in y_cal])
scores_cal = 1.0 - proba_cal[np.arange(len(y_cal)), true_idx_cal]

# 2. Quantile q̂ avec correction finie — garantie de couverture exacte
n_cal = len(scores_cal)
q_hat = np.quantile(scores_cal, ceil((n_cal + 1) * (1 - alpha)) / n_cal)
print(f"Modèle utilisé : {best_name}")
print(f"n_cal = {n_cal:,}")
print(f"q̂    = {q_hat:.4f}  →  on retient toutes les classes avec P(classe|x) ≥ {1-q_hat:.4f}")

# 3. Ensembles de prédiction sur le test
proba_te  = best_model.predict_proba(X_te)   # (n_te, 6)
pred_sets = proba_te >= (1.0 - q_hat)        # booléen (n_te, 6)
set_sizes = pred_sets.sum(axis=1)            # taille de l'ensemble pour chaque pixel

# 4. Couverture empirique : la vraie classe est-elle dans l'ensemble ?
true_idx_te = np.array([np.where(classes_best == yc)[0][0] for yc in y_te])
covered     = pred_sets[np.arange(len(y_te)), true_idx_te]

print(f"\nCouverture empirique  : {covered.mean():.3f}  (cible ≥ {1-alpha:.2f})")
print(f"Taille moy. ensemble  : {set_sizes.mean():.2f} classes / pixel")
print(f"\nDistribution des tailles d'ensemble :")
for sz in sorted(np.unique(set_sizes)):
    n = (set_sizes == sz).sum()
    print(f"  {sz} classe(s) : {n:5d} pixels ({100*n/len(set_sizes):.1f} %)")

# 5. Visualisation : carte des prédictions + carte d'incertitude conforme
y_pred_best = predict_labels(best_name, X_te)
acc_best = accuracy_score(y_te, y_pred_best)

fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# Panneau gauche : classe la plus probable du modèle retenu (fond MNT)
axes[0].imshow(mnt, cmap='gray', origin='upper', alpha=0.35)
axes[0].scatter(
    cols_te, rows_te,
    c=y_pred_best,
    cmap=mcolors.ListedColormap([LC_COLORS[v] for v in sorted(LC_NAMES)]),
    vmin=0.5, vmax=6.5, s=6, alpha=0.8, linewidths=0
)
patches_lc = [Patch(color=LC_COLORS[v], label=f'{v} – {LC_NAMES[v]}')
              for v in sorted(LC_NAMES)]
axes[0].legend(handles=patches_lc, loc='lower right', fontsize=8, title='Classe prédite')
axes[0].set_title(f'Classe prédite ({best_name})\nacc = {acc_best:.3f}', fontsize=11)
axes[0].set_xlim(0, mnt.shape[1])
axes[0].set_ylim(mnt.shape[0], 0)
axes[0].axis('off')

# Panneau droit : incertitude = taille de l'ensemble de prédiction conforme
size_cmap = mcolors.ListedColormap(['#1a9850', '#fdae61', '#d73027', '#7b2d8b'])
size_norm = mcolors.BoundaryNorm([0.5, 1.5, 2.5, 3.5, 6.5], size_cmap.N)
axes[1].imshow(mnt, cmap='gray', origin='upper', alpha=0.35)
axes[1].scatter(
    cols_te, rows_te,
    c=set_sizes, cmap=size_cmap, norm=size_norm,
    s=6, alpha=0.8, linewidths=0
)
size_patches = [
    Patch(color='#1a9850', label='1 classe  (très certain)'),
    Patch(color='#fdae61', label='2 classes (ambigu)'),
    Patch(color='#d73027', label='3 classes (incertain)'),
    Patch(color='#7b2d8b', label='4+ classes (très incertain)'),
]
axes[1].legend(handles=size_patches, loc='lower right', fontsize=8, title='Incertitude')
axes[1].set_title(
    f'Incertitude conforme (α = {alpha})\n'
    f'taille moy. = {set_sizes.mean():.2f}  |  couverture = {covered.mean():.3f}',
    fontsize=11
)
axes[1].set_xlim(0, mnt.shape[1])
axes[1].set_ylim(mnt.shape[0], 0)
axes[1].axis('off')

plt.suptitle(
    f'Prédiction conforme pour la classification ({best_name}, LAC, α = {alpha})\n'
    f'q̂ = {q_hat:.4f} → ensembles à ≥ {(1-alpha)*100:.0f} % de couverture garantie',
    fontsize=12, fontweight='bold'
)
plt.tight_layout()
plt.show()