import geopandas as gpd
import pandas as pd
import sqlite3
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np

# Chemin vers les données du TP2 : Remplace par le chemin vers ton fichier tp2.gpkg
DATA_PATH = 'tp2.gpkg'

# --- Chargement des couches géographiques ---
communes  = gpd.read_file(DATA_PATH, layer='Communes')
geologie  = gpd.read_file(DATA_PATH, layer='geology_VD')
accidents = gpd.read_file(DATA_PATH, layer='LI_Accident_tecto')
batiments = gpd.read_file(DATA_PATH, layer='BatimentsUNIL')
parcelles = gpd.read_file(DATA_PATH, layer='Parcelles')


for nom, gdf in [('Communes', communes), ('geology_VD', geologie),
                 ('LI_Accident_tecto', accidents), ('BatimentsUNIL', batiments),
                 ('Parcelles', parcelles)]:
    print(f"{nom:<25}  {str(gdf.shape):<12}  "
          f"CRS: EPSG:{gdf.crs.to_epsg()}  geom: {gdf.geom_type.unique()[0]}")

# --- Vue d'ensemble : superposer les couches du canton de Vaud ---
fig, ax = plt.subplots(figsize=(10, 8))

communes.plot(ax=ax, color='#f5f0e8', edgecolor='#888', linewidth=0.5)
geologie.plot(ax=ax, column='TECTO_F', cmap='tab20', alpha=0.4, legend=False)
accidents.plot(ax=ax, color='darkred', linewidth=0.4, alpha=0.6)

patches = [
    mpatches.Patch(facecolor='#f5f0e8', edgecolor='#888', label='Communes'),
    mpatches.Patch(color='steelblue', alpha=0.6, label='Géologie (TECTO_F)'),
    mpatches.Patch(color='darkred', label='Accidents tectoniques'),
]
ax.legend(handles=patches, loc='lower right', fontsize=9)
ax.set_title('Canton de Vaud — données tp2.gpkg', fontsize=13)
ax.set_axis_off()
plt.tight_layout()
plt.show()

# --- Tables attributaires (sans géométrie) ---
# Ces tables sont stockées dans le GeoPackage mais sans colonne géométrique.
# On les lit via sqlite3 (ou pd.read_sql).
con = sqlite3.connect(DATA_PATH)
voitures      = pd.read_sql('SELECT * FROM VoituresPlus', con)
proprietaires = pd.read_sql('SELECT * FROM Proprietaires', con)
prop_parc     = pd.read_sql('SELECT * FROM Proprietaire_Parcelle', con)
con.close()

print(f"VoituresPlus          : {voitures.shape}  cols={voitures.columns.tolist()}")
print(f"Proprietaires         : {proprietaires.shape}  cols={proprietaires.columns.tolist()}")
print(f"Proprietaire_Parcelle : {prop_parc.shape}  cols={prop_parc.columns.tolist()}")
print()
voitures.head()

# --- 2.1 Sélection simple ---
nappe_morcles = geologie[
    geologie['LEG_TEC_3'] == 'Nappe de Morcles (Chaine des Aravis incl.)'
]

print(f"Entités sélectionnées : {len(nappe_morcles)}")
print(f"Surface totale        : {nappe_morcles['AREA'].sum() / 1e6:.1f} km²")

fig, ax = plt.subplots(figsize=(10, 7))
communes.plot(ax=ax, color='#f5f0e8', edgecolor='#aaa', linewidth=0.4)
nappe_morcles.plot(ax=ax, color='#e63946', edgecolor='#900', linewidth=0.5)
ax.set_title("Nappe de Morcles (Chaine des Aravis incl.)", fontsize=12)
ax.set_axis_off()
# Légende manuelle
legend_patches = [
    mpatches.Patch(color='#e63946', label='Nappe de Morcles'),
    mpatches.Patch(facecolor='#f5f0e8', edgecolor='#aaa', label='Communes'),
]
ax.legend(handles=legend_patches, fontsize=9, loc='lower right')
plt.tight_layout()
plt.show()

# --- 2.2 Sélection par intervalle ---
petites_unites = geologie[
    (geologie['AREA'] >= 100_000) &
    (geologie['AREA'] <= 500_000)
]

print(f"Entités sélectionnées : {len(petites_unites)}")
print(f"Surface min / max     : {petites_unites['AREA'].min()/1e6:.3f} – "
      f"{petites_unites['AREA'].max()/1e6:.3f} km²")
print()
print("Top 8 unités les plus grandes dans la sélection :")
print(petites_unites[['OBJECTID', 'LEG_TEC_3', 'AREA']]
      .sort_values('AREA', ascending=False).head(8).to_string(index=False))

# --- 2.3 Sélection multiple avec .isin() ---
productivites = [
    'Peu productifs, dans les moraines',
    'Productif, a productivite variable ou faible'
]
peu_productifs = geologie[geologie['PRODUCTIV'].isin(productivites)]

print(f"Entités sélectionnées : {len(peu_productifs)}")
print()
print("Répartition par classe :")
print(peu_productifs['PRODUCTIV'].value_counts().to_string())

fig, ax = plt.subplots(figsize=(10, 7))
communes.plot(ax=ax, color='#f5f0e8', edgecolor='#aaa', linewidth=0.4)
peu_productifs.plot(ax=ax, column='PRODUCTIV', cmap='Set2', legend=True,
                    legend_kwds={'loc': 'lower right', 'fontsize': 8})
ax.set_title("Aquifères peu productifs — classes sélectionnées", fontsize=12)
ax.set_axis_off()
plt.tight_layout()
plt.show()

# --- 2.4 Requêtes avec .query() ---

# Sélection simple
nappe_q = geologie.query("LEG_TEC_3 == 'Nappe de Morcles (Chaine des Aravis incl.)'")

# Sélection par intervalle
petites_q = geologie.query("AREA >= 100_000 and AREA <= 500_000")

# Variable externe (préfixe @)
seuil_bas, seuil_haut = 1_000_000, 10_000_000
grandes_unites = geologie.query("AREA >= @seuil_bas and AREA <= @seuil_haut")

# Condition combinée avec méthode de chaîne (nécessite engine='python')
molasse = geologie.query(
    "LEG_TEC_3.str.contains('Molasse', na=False) and AREA > 1_000_000",
    engine='python'
)

print(f"Nappe de Morcles         : {len(nappe_q)} entités")
print(f"Petites unités (0.1–0.5 km²) : {len(petites_q)} entités")
print(f"Grandes unités (1–10 km²): {len(grandes_unites)} entités")
print(f"Molasse (> 1 km²)        : {len(molasse)} entités")

# --- 3. Requête spatiale : chevauchements alpins dans Vaud ---

# Étape 1 : sélection attributaire
types_chev = [
    'Chevauchement principal alpin (certain)',
    'Chevauchement principal alpin (probable)'
]
chev = accidents[accidents['Type'].isin(types_chev)].copy()
print(f"Chevauchements alpins (total dataset)  : {len(chev)}")
print(chev['Type'].value_counts().to_string())

# Étape 2 : sélection spatiale
# gpd.sjoin(predicate='intersects') → garde les lignes qui croisent un polygone communal
chev_vaud = gpd.sjoin(
    chev,
    communes[['geometry']],
    how='inner',
    predicate='intersects'
)
chev_vaud = chev_vaud[~chev_vaud.index.duplicated()]   # supprimer les doublons
print(f"\nChevauchements dans Vaud               : {len(chev_vaud)}")

# Visualisation
fig, ax = plt.subplots(figsize=(10, 7))
communes.plot(ax=ax, color='#f5f0e8', edgecolor='#bbb', linewidth=0.4)
chev_vaud.plot(ax=ax, column='Type', cmap='Set1', linewidth=1.2,
               legend=True, legend_kwds={'loc': 'lower right', 'fontsize': 8})
ax.set_title('Chevauchements principaux alpins — canton de Vaud', fontsize=12)
ax.set_axis_off()
plt.tight_layout()
plt.show()

voitures.head()

communes.head()

# --- 4.1 Jointure attributaire ---
communes_voit = communes.merge(
    voitures[['Communes', 'Voitures_de_tourisme', 'Voitures_de_tourisme__1000_hab']],
    left_on='NAME',
    right_on='Communes',
    how='left'
)

print(f"Résultat : {communes_voit.shape}")
print(f"Communes sans données  : "
      f"{communes_voit['Voitures_de_tourisme__1000_hab'].isna().sum()} / {len(communes_voit)}")
print()
top10 = (communes_voit.dropna(subset=['Voitures_de_tourisme__1000_hab'])
         .sort_values('Voitures_de_tourisme__1000_hab', ascending=False)
         .head(10))
print("Top 10 communes — véhicules / 1 000 hab. :")
print(top10[['NAME', 'Voitures_de_tourisme', 'Voitures_de_tourisme__1000_hab']]
      .to_string(index=False))

# --- 4.2 Jointure spatiale : bâtiments UNIL → commune ---
# predicate='intersects' pour capturer aussi les bâtiments sur les limites
bat_communes = gpd.sjoin(
    batiments,
    communes[['NAME', 'BEZIRK', 'geometry']],
    how='left',
    predicate='intersects'
)
# Un bâtiment peut intersecter plusieurs communes : garder la première occurrence
bat_communes = bat_communes[~bat_communes.index.duplicated(keep='first')]

print(f"Résultat : {bat_communes.shape}")
print()
print("Commune de chaque bâtiment UNIL :")
print(bat_communes[['name', 'amenity', 'NAME', 'BEZIRK']]
      .sort_values('NAME').to_string(index=False))

# --- 4.3 Relations N-N : Parcelles ↔ Propriétaires ---
# Jointure en chaîne : Parcelles → Proprietaire_Parcelle → Proprietaires
parcelles_propri = (
    parcelles[['NO_IMM', 'NO_COMMUNE', 'SURFACE_RF', 'geometry']]
    .merge(prop_parc[['NO_IMM', 'NO_PROPRI']], on='NO_IMM', how='left')
    .merge(proprietaires[['NO_PROPRI', 'NOM_PROPRI']], on='NO_PROPRI', how='left')
)

print(f"Résultat : {parcelles_propri.shape}")
print(f"(une ligne par couple parcelle–propriétaire)")
print()

# Parcelles avec plusieurs propriétaires
multi_propri = (parcelles_propri.groupby('NO_IMM')
                .filter(lambda x: len(x) > 1)
                .sort_values('NO_IMM'))
print("Parcelles avec plusieurs propriétaires :")
print(multi_propri[['NO_IMM', 'SURFACE_RF', 'NO_PROPRI', 'NOM_PROPRI']]
      .to_string(index=False))

# --- 5.1 Choroplèthe : véhicules / 1 000 hab. par commune ---
fig, ax = plt.subplots(figsize=(12, 9))

communes_voit.plot(
    ax=ax,
    column='Voitures_de_tourisme__1000_hab',
    cmap='YlOrRd',
    scheme='quantiles',
    k=5,
    legend=True,
    legend_kwds={'title': 'Véhicules/1 000 hab.', 'loc': 'lower left'},
    edgecolor='white',
    linewidth=0.3,
    missing_kwds={'color': 'lightgrey', 'label': 'Données manquantes'}
)

ax.set_title('Voitures de tourisme pour 1 000 habitants\nCommunes vaudoises (2010)',
             fontsize=13, pad=10)
ax.set_axis_off()
plt.tight_layout()
plt.show()

# --- 5.2 Comparaison des méthodes de discrétisation ---
methodes = ['quantiles', 'equal_interval', 'natural_breaks']
titres   = ['Quantiles (5 classes)', 'Intervalles égaux (5 classes)', 'Jenks (5 classes)']

fig, axes = plt.subplots(1, 3, figsize=(18, 6))

for ax, methode, titre in zip(axes, methodes, titres):
    communes_voit.plot(
        ax=ax,
        column='Voitures_de_tourisme__1000_hab',
        cmap='YlOrRd', scheme=methode, k=5,
        legend=True, legend_kwds={'fontsize': 7},
        edgecolor='white', linewidth=0.3,
        missing_kwds={'color': 'lightgrey'}
    )
    ax.set_title(titre, fontsize=10)
    ax.set_axis_off()

plt.suptitle("Véhicules / 1 000 hab. — méthodes de discrétisation", fontsize=13, y=1.01)
plt.tight_layout()
plt.show()

# --- Carte choroplèthe complète avec éléments cartographiques ---
fig, ax = plt.subplots(figsize=(12, 10))

communes_voit.plot(
    ax=ax,
    column='Voitures_de_tourisme__1000_hab',
    cmap='YlOrRd',
    scheme='natural_breaks', k=5,
    legend=True,
    legend_kwds={
        'title': 'Véhicules\n/1 000 hab.',
        'loc': 'lower right',
        'fontsize': 9,
        'title_fontsize': 10,
    },
    edgecolor='white', linewidth=0.3,
    missing_kwds={'color': '#cccccc', 'label': 'Données manquantes'}
)

ax.set_title("Voitures de tourisme pour 1 000 habitants\nCommunes vaudoises — 2010",
             fontsize=13, pad=12)
ax.set_axis_off()
ax.annotate('Source : OFS / tp2.gpkg', xy=(0.01, 0.01), xycoords='axes fraction',
            fontsize=8, color='grey')
plt.tight_layout()
plt.show()

import pandas as pd
import numpy as np
from scipy import stats
import statsmodels.formula.api as smf
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_absolute_error

# --- Chargement des résultats de votation (tables du GeoPackage tp2.gpkg) ---
con = sqlite3.connect(DATA_PATH)
v23_raw = pd.read_sql('SELECT * FROM votation_2023', con)
v26_raw = pd.read_sql('SELECT * FROM votation_2026', con)
con.close()

# Dans un GeoPackage, les colonnes numériques peuvent être stockées en texte :
# on les convertit explicitement avant de filtrer et de calculer des pourcentages.
for d in (v23_raw, v26_raw):
    for col in ['proposal_id', 'yes_votes_count', 'cast_votes_count']:
        d[col] = pd.to_numeric(d[col], errors='coerce')

# Niveau communal uniquement, filtre par identifiant de scrutin
lci = v23_raw[(v23_raw['region_type'] == 'municipality') &
              (v23_raw['proposal_id'] == 6630)].copy()   # Loi Climat
m10 = v26_raw[(v26_raw['region_type'] == 'municipality') &
              (v26_raw['proposal_id'] == 6860)].copy()   # Initiative 10M

lci['oui_pct_lci'] = lci['yes_votes_count'] / lci['cast_votes_count'] * 100
m10['oui_pct_m10'] = m10['yes_votes_count'] / m10['cast_votes_count'] * 100

# Jointure interne sur l'identifiant communal : seules les communes présentes
# dans les deux scrutins sont conservées
df = (lci[['region_id', 'region_name', 'parent_canton_name', 'oui_pct_lci']]
      .merge(m10[['region_id', 'oui_pct_m10']], on='region_id')
      .dropna())

print(f"Communes retenues : {len(df)}")
print(f"  LCI — taux OUI : moy={df['oui_pct_lci'].mean():.1f}%  "
      f"[{df['oui_pct_lci'].min():.1f} – {df['oui_pct_lci'].max():.1f}%]")
print(f"  10M — taux OUI : moy={df['oui_pct_m10'].mean():.1f}%  "
      f"[{df['oui_pct_m10'].min():.1f} – {df['oui_pct_m10'].max():.1f}%]")
df.head()

r, p_val = stats.pearsonr(df['oui_pct_lci'], df['oui_pct_m10'])
print(f"r = {r:.3f}   p-value = {p_val:.2e}")
print(f"→ Corrélation {'forte' if abs(r)>0.7 else 'modérée'} et "
      f"{'négative' if r < 0 else 'positive'} (R² = {r**2:.3f})")

fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(df['oui_pct_lci'], df['oui_pct_m10'],
           alpha=0.25, s=6, color='steelblue', linewidths=0)
ax.set_xlabel('% OUI — Loi Climat 2023 (LCI)', fontsize=11)
ax.set_ylabel('% OUI — Initiative 10 millions 2026', fontsize=11)
ax.set_title(f'Corrélation entre deux votations fédérales suisses\n'
             f'n = {len(df)} communes   r = {r:.3f}   R² = {r**2:.3f}', fontsize=11)
plt.tight_layout()
plt.show()

# Ajustement OLS avec statsmodels (formule R-like)
model_ols = smf.ols('oui_pct_m10 ~ oui_pct_lci', data=df).fit()
print(model_ols.summary().tables[1])   # tableau des coefficients seulement
print(f"\nR² = {model_ols.rsquared:.3f}")

# Prédictions sur une grille régulière pour tracer IC et IP
x_grid = pd.DataFrame({'oui_pct_lci': np.linspace(df['oui_pct_lci'].min(),
                                                    df['oui_pct_lci'].max(), 300)})
pred = model_ols.get_prediction(x_grid).summary_frame(alpha=0.05)
# Colonnes : mean, mean_ci_lower, mean_ci_upper, obs_ci_lower, obs_ci_upper

fig, ax = plt.subplots(figsize=(9, 6))
ax.scatter(df['oui_pct_lci'], df['oui_pct_m10'],
           alpha=0.2, s=6, color='steelblue', linewidths=0, zorder=1)
ax.plot(x_grid['oui_pct_lci'], pred['mean'],
        color='crimson', linewidth=2, label='Droite OLS', zorder=3)
ax.fill_between(x_grid['oui_pct_lci'], pred['mean_ci_lower'], pred['mean_ci_upper'],
                color='crimson', alpha=0.25, label='IC 95 % (moyenne)', zorder=2)
ax.fill_between(x_grid['oui_pct_lci'], pred['obs_ci_lower'], pred['obs_ci_upper'],
                color='orange', alpha=0.15, label='IP 95 % (observation)', zorder=2)
ax.set_xlabel('% OUI — Loi Climat (LCI)', fontsize=11)
ax.set_ylabel('% OUI — Initiative 10 millions', fontsize=11)
ax.set_title('Régression OLS avec intervalles de confiance et de prédiction', fontsize=11)
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

rng = np.random.default_rng(0)
sample_idx = rng.choice(len(df), size=15, replace=False)
rest_idx   = np.setdiff1d(np.arange(len(df)), sample_idx)

x_s, y_s = df['oui_pct_lci'].values[sample_idx], df['oui_pct_m10'].values[sample_idx]
x_r, y_r = df['oui_pct_lci'].values[rest_idx],   df['oui_pct_m10'].values[rest_idx]

degrees = [1, 6]
x_line = np.linspace(df['oui_pct_lci'].min(), df['oui_pct_lci'].max(), 300)

fig, ax = plt.subplots(figsize=(9, 6))
ax.scatter(x_r, y_r, alpha=0.08, s=6, color='grey', label='Reste des communes (jamais vues)')
ax.scatter(x_s, y_s, s=40, color='black', zorder=5, label="Échantillon d'entraînement (n=15)")

results = []
for deg, color in zip(degrees, ['crimson', 'orange']):
    coeffs = np.polyfit(x_s, y_s, deg)
    ax.plot(x_line, np.polyval(coeffs, x_line), color=color, linewidth=2, label=f'Degré {deg}')

    r2_train   = r2_score(y_s, np.polyval(coeffs, x_s))
    r2_heldout = r2_score(y_r, np.polyval(coeffs, x_r))
    results.append((deg, r2_train, r2_heldout))

ax.set_ylim(y_r.min() - 10, y_r.max() + 10)   # les oscillations du degré 6 sortent vite du cadre
ax.set_xlabel('% OUI — LCI'); ax.set_ylabel('% OUI — 10M')
ax.legend(fontsize=9)
ax.set_title('Surapprentissage : droite vs polynôme de degré 6 (ajustés sur n=15)')
plt.tight_layout()
plt.show()

print(f"{'Degré':<8}{'R² (entraînement, n=15)':<26}{'R² (reste, jamais vu)':<24}")
for deg, r2_tr, r2_ho in results:
    print(f"{deg:<8}{r2_tr:<26.3f}{r2_ho:<24.3f}")
print("\nLe polynôme de degré 6 explique mieux les 15 points d'entraînement que la droite,")
print("mais généralise beaucoup plus mal (R² très négatif, bien pire que de prédire la")
print("moyenne) : il a appris le bruit de cet échantillon précis, pas la tendance réelle.")

X = df[['oui_pct_lci']].values
y = df['oui_pct_m10'].values

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42)

model_val = LinearRegression().fit(X_train, y_train)
y_pred_test = model_val.predict(X_test)

r2 = r2_score(y_test, y_pred_test)
print(f"Entraînement : {len(X_train)} communes   Test : {len(X_test)} communes")
print(f"R² (test) = {r2:.3f}")

mae  = mean_absolute_error(y_test, y_pred_test)
rmse = np.sqrt(np.mean((y_test - y_pred_test) ** 2))

# --- Benchmark : modèle naïf qui prédit toujours la moyenne d'entraînement ---
y_pred_baseline = np.full_like(y_test, y_train.mean())
r2_base   = r2_score(y_test, y_pred_baseline)
mae_base  = mean_absolute_error(y_test, y_pred_baseline)
rmse_base = np.sqrt(np.mean((y_test - y_pred_baseline) ** 2))

print(f"{'Métrique':<10}{'Benchmark (moyenne)':<22}{'Modèle OLS':<12}")
print(f"{'R²':<10}{r2_base:<22.3f}{r2:<12.3f}")
print(f"{'MAE':<10}{mae_base:<22.2f}{mae:<12.2f}")
print(f"{'RMSE':<10}{rmse_base:<22.2f}{rmse:<12.2f}")
print()
print(f"Le benchmark a un R² ≈ 0 : sans surprise, « toujours prédire la moyenne » n'explique")
print(f"aucune variance (§6.2). Le modèle OLS réduit le MAE de "
      f"{100 * (1 - mae / mae_base):.0f} % par rapport à ce benchmark.")
print(f"Interprétation : en moyenne, l'erreur de prédiction du modèle est ±{mae:.1f} points de %,")
print(f"contre ±{mae_base:.1f} points de % pour le benchmark — le modèle apporte une information réelle,")
print(f"pas seulement un ajustement qui a l'air raisonnable sur le papier.")

# Géométries des communes suisses (couche Communes_CH de tp2.gpkg) pour cartographier les résidus.
# GEM_TEIL == '0' : polygone principal de chaque commune (pas les enclaves)
communes_ch = (gpd.read_file(DATA_PATH, layer='Communes_CH')
               .to_crs('EPSG:21781'))
communes_ch = communes_ch[
    communes_ch['KANTONSNUM'].str.startswith('CH', na=False) &
    (communes_ch['GEM_TEIL'] == '0')
]

# Régions linguistiques (simplifié — cantons officiellement francophones ou italophones)
romand = {'Vaud', 'Genève', 'Neuchâtel', 'Jura', 'Fribourg', 'Valais'}
tessin = {'Tessin'}
def region(canton):
    if canton in romand: return 'Romand'
    if canton in tessin: return 'Tessin'
    return 'Alémanique'

df['region_ling'] = df['parent_canton_name'].map(region)
print("Distribution linguistique :")
print(df['region_ling'].value_counts().to_string())

mask_al = (df['region_ling'] == 'Alémanique').values
n_tr, n_te = mask_al.sum(), (~mask_al).sum()

# --- Split géographique : entraînement = Alémanique, test = Romandie + Tessin ---
model_geo = LinearRegression().fit(X[mask_al], y[mask_al])
resid_geo_arr = y[~mask_al] - model_geo.predict(X[~mask_al])
r2_geo  = r2_score(y[~mask_al], model_geo.predict(X[~mask_al]))
mae_geo = mean_absolute_error(y[~mask_al], model_geo.predict(X[~mask_al]))

# --- Split aléatoire à effectif IDENTIQUE (n_tr entraînement / n_te test) ---
# Un modèle réajusté sur toutes les communes (comme model_ols en §6.3) aurait déjà "vu"
# les communes de test : ce ne serait pas une évaluation hors échantillon valable.
rng_split = np.random.default_rng(42)
perm = rng_split.permutation(len(df))
idx_tr_rand, idx_te_rand = perm[:n_tr], perm[n_tr:n_tr + n_te]

model_rand = LinearRegression().fit(X[idx_tr_rand], y[idx_tr_rand])
resid_rand_arr = y[idx_te_rand] - model_rand.predict(X[idx_te_rand])
r2_rand  = r2_score(y[idx_te_rand], model_rand.predict(X[idx_te_rand]))
mae_rand = mean_absolute_error(y[idx_te_rand], model_rand.predict(X[idx_te_rand]))

bias_rand, bias_geo = resid_rand_arr.mean(), resid_geo_arr.mean()

print()
print(f"Effectif (identique dans les deux cas) : {n_tr} entraînement / {n_te} test")
print(f"Split aléatoire     : R² = {r2_rand:.3f}   MAE = {mae_rand:.2f} pts%   biais = {bias_rand:+.2f} pts%")
print(f"Split géographique  : R² = {r2_geo:.3f}   MAE = {mae_geo:.2f} pts%   biais = {bias_geo:+.2f} pts%")
print("  (entraîné sur la Suisse alémanique → prédit Romandie + Tessin)")

# --- Rattacher chaque commune à sa géométrie, avec son rôle (entraînement / test) ---
geo = communes_ch[['NAME', 'geometry']]

role_rand = np.empty(len(df), dtype=object)
role_rand[idx_tr_rand] = 'train'
role_rand[idx_te_rand] = 'test'
df_split = df[['region_name']].copy()
df_split['role_rand'] = role_rand
df_split['role_geo']  = np.where(mask_al, 'train', 'test')

gdf_rand_all = geo.merge(df_split[['region_name', 'role_rand']],
                         left_on='NAME', right_on='region_name', how='inner')
gdf_geo_all  = geo.merge(df_split[['region_name', 'role_geo']],
                         left_on='NAME', right_on='region_name', how='inner')

# Résidus rattachés à la géométrie (jeu de test uniquement)
df_rand_te = df.iloc[idx_te_rand][['region_name']].copy()
df_rand_te['resid'] = resid_rand_arr
gdf_rand = geo.merge(df_rand_te, left_on='NAME', right_on='region_name', how='inner')

df_geo_te = df.loc[~mask_al, ['region_name']].copy()
df_geo_te['resid'] = resid_geo_arr
gdf_geo = geo.merge(df_geo_te, left_on='NAME', right_on='region_name', how='inner')

vabs = max(gdf_rand['resid'].abs().quantile(0.97), gdf_geo['resid'].abs().quantile(0.97))

# --- Figure 2×2 : en haut la stratégie de séparation, en bas les résidus ---
fig, axes = plt.subplots(2, 2, figsize=(13, 11))
C_TRAIN, C_TEST = '#cfd8dc', '#ef6c00'

# Rangée du haut : QUI est en entraînement (gris) vs en test (orange)
for ax, gdf_all, role_col, titre in [
    (axes[0, 0], gdf_rand_all, 'role_rand', f'Split aléatoire — test dispersé (n={n_te})'),
    (axes[0, 1], gdf_geo_all,  'role_geo',  f'Split géographique — test = Romandie + Tessin (n={n_te})'),
]:
    gdf_all[gdf_all[role_col] == 'train'].plot(ax=ax, color=C_TRAIN, edgecolor='white', linewidth=0.15)
    gdf_all[gdf_all[role_col] == 'test'].plot(ax=ax, color=C_TEST, edgecolor='white', linewidth=0.15)
    ax.set_title(titre, fontsize=10)
    ax.set_axis_off()
axes[0, 0].legend(handles=[mpatches.Patch(facecolor=C_TRAIN, label='Entraînement'),
                           mpatches.Patch(facecolor=C_TEST, label='Test')],
                  loc='lower left', fontsize=9, frameon=False)

# Rangée du bas : les résidus qui en résultent, même échelle de couleur des deux côtés
for ax, gdf_resid, r2v, maev, biasv, titre in [
    (axes[1, 0], gdf_rand, r2_rand, mae_rand, bias_rand, 'Résidus — split aléatoire'),
    (axes[1, 1], gdf_geo,  r2_geo,  mae_geo,  bias_geo,  'Résidus — split géographique'),
]:
    communes_ch.plot(ax=ax, color='#f0f0f0', edgecolor='white', linewidth=0.2)
    gdf_resid.plot(ax=ax, column='resid', cmap='PiYG', vmin=-vabs, vmax=vabs, legend=True,
                   legend_kwds={'shrink': 0.6, 'label': 'résidu (points %)'}, edgecolor='none')
    ax.set_title(f'{titre}\nR² = {r2v:.2f}   MAE = {maev:.2f}   biais = {biasv:+.2f} pts%', fontsize=10)
    ax.set_axis_off()

plt.suptitle("Global (aléatoire) vs localisé (géographique) : même effectif d'entraînement, résultats très différents",
             fontsize=13, y=0.99)
plt.tight_layout()
plt.show()