# ============================================================================
# SCRIPT DIAGNOSTIC SOCIAL DES COMMUNES DES PYRENEES-ORIENTALES (66)
# ============================================================================
# Etapes du script :
#   1. Importer et préparer les six indicateurs
#   2. Fusionner les données communales
#   3. Estimer les revenus médians manquants et évaluer le modèle
#   4. Calculer les scores de vulnérabilité et le score global
#   5. Choisir le nombre de clusters (coude + silhouette)
#   6. Appliquer le K-means et produire les graphiques utiles
#   7. Exporter les tableaux et les graphiques

# ============================================================================

# ----------------------------------------------------------------------------
# 0. PARAMETRES ET DEPENDANCES
# ----------------------------------------------------------------------------

# Packages utilisés
packages_requis <- c(
  "readxl",
  "dplyr",
  "stringr",
  "tidyr",
  "ggplot2",
  "factoextra",
  "cluster",
  "plotly",
  "htmlwidgets",
  "jsonlite"
)

# Pour retrouver les mêmes résultats à chaque exécution
graine <- 66

departement <- "66"

dossier_projet <- normalizePath(".", mustWork = TRUE)
dossier_resultats <- file.path(dossier_projet, "resultats")

dir.create(
  dossier_resultats,
  showWarnings = FALSE,
  recursive = TRUE
)

# Fichiers sources
chemins <- list(
  chomage = file.path(
    dossier_projet,
    "base_cc_emploi_pop_active_2023_xlsx.xlsx"
  ),
  revenus = file.path(
    dossier_projet,
    "FILOSOFI_CC_FR.xlsx"
  ),
  familles = file.path(
    dossier_projet,
    "base_cc_coupl_fam-men_2023.xlsx"
  ),
  population = file.path(
    dossier_projet,
    "TD_POP1_2023.xlsx"
  ),
  rsa = file.path(
    dossier_projet,
    "rsa.xlsx"
  )
)

# Contours des communes pour les cartes
dossier_donnees <- file.path(dossier_projet, "data")
chemin_geojson <- file.path(dossier_donnees, "communes_66.geojson")
url_geojson <- paste0(
  "https://geo.api.gouv.fr/departements/",
  departement,
  "/communes?format=geojson&geometry=contour"
)

if (!file.exists(chemin_geojson)) {
  dir.create(
    dossier_donnees,
    showWarnings = FALSE,
    recursive = TRUE
  )

  download.file(
    url_geojson,
    chemin_geojson,
    mode = "wb",
    quiet = TRUE
  )
}

# ----------------------------------------------------------------------------
# 1. TAUX DE CHOMAGE ET DIPLOMES DES ACTIFS
# ----------------------------------------------------------------------------

# Code et nom des communes
geo_chomage <- read_excel(
  chemins$chomage,
  sheet = "COM_2023",
  range = cell_cols("A:B")
)

# Nombre d'actifs 15-64 ans
actifs <- read_excel(
  chemins$chomage,
  sheet = "COM_2023",
  range = cell_cols("O")
)

# Nombre de chômeurs
chomeurs <- read_excel(
  chemins$chomage,
  sheet = "COM_2023",
  range = cell_cols("AM")
)

# Diplômes des actifs
diplomes_actifs <- read_excel(
  chemins$chomage,
  sheet = "COM_2023",
  range = cell_cols("BA:BD")
)

chomage_66 <- bind_cols(
  geo_chomage,
  actifs,
  chomeurs,
  diplomes_actifs
) |>
  transmute(
    code_commune = as.character(`Code géographique`),
    commune = `Libellé géographique`,
    actifs_15_64 = `Actifs 15-64 ans (princ)`,
    chomeurs_15_64 = `Chômeurs 15-64 ans (princ)`,
    actifs_bac = `Actifs Bac, brevet pro. ou équiv.  (princ)`,
    actifs_bac5_plus =
      `Actifs Enseignement sup de niveau bac + 5 ou plus  (princ)`
  ) |>
  filter(str_starts(code_commune, departement)) |>
  mutate(
    # Taux calculés parmi les actifs de 15 à 64 ans
    taux_chomage = if_else(
      actifs_15_64 > 0,
      100 * chomeurs_15_64 / actifs_15_64,
      NA_real_
    ),
    part_actifs_bac = if_else(
      actifs_15_64 > 0,
      100 * actifs_bac / actifs_15_64,
      NA_real_
    ),
    part_actifs_bac5_plus = if_else(
      actifs_15_64 > 0,
      100 * actifs_bac5_plus / actifs_15_64,
      NA_real_
    )
  )

# ----------------------------------------------------------------------------
# 2. REVENU MEDIAN
# ----------------------------------------------------------------------------

# Les valeurs s et vm deviennent des NA
revenus_66 <- read_excel(
  chemins$revenus,
  sheet = "COM",
  skip = 5,
  col_names = c(
    "code_commune",
    "commune",
    "niveau_vie_median",
    "taux_pauvrete"
  ),
  col_types = c("text", "text", "numeric", "numeric"),
  na = c("s", "vm")
) |>
  filter(str_starts(code_commune, departement)) |>
  select(
    code_commune,
    niveau_vie_median
  )

# ----------------------------------------------------------------------------
# 3. PART DES FAMILLES MONOPARENTALES
# ----------------------------------------------------------------------------

# Données sur les familles
geo_familles <- read_excel(
  chemins$familles,
  sheet = "COM_2023",
  range = cell_cols("A:B")
)

mesures_familles <- read_excel(
  chemins$familles,
  sheet = "COM_2023",
  range = cell_cols("BP:BR")
)

familles_66 <- bind_cols(
  geo_familles,
  mesures_familles
) |>
  transmute(
    code_commune = as.character(`Code géographique`),
    familles_total = `Familles (compl)`,
    familles_monoparentales = `Fam Monoparentales (compl)`
  ) |>
  filter(str_starts(code_commune, departement)) |>
  mutate(
    # Part parmi l'ensemble des familles
    part_monoparentales = if_else(
      familles_total > 0,
      100 * familles_monoparentales / familles_total,
      NA_real_
    )
  )

# ----------------------------------------------------------------------------
# 4. POPULATION TOTALE POUR L'ESTIMATION DU REVENU
# ----------------------------------------------------------------------------

# Lecture des noms de colonnes
noms_population <- names(
  read_excel(
    chemins$population,
    sheet = "COM",
    range = "A1:GV1"
  )
)

# Lignes des communes du 66 uniquement
population_66_brute <- read_excel(
  chemins$population,
  sheet = "COM",
  range = "A25821:GV26046",
  col_names = noms_population,
  col_types = c("text", "text", rep("numeric", 202))
)

colonnes_population <- names(population_66_brute)[3:204]

# Somme de toutes les classes d'âge
population_totale_calculee <- rowSums(
  population_66_brute[, colonnes_population],
  na.rm = TRUE
)

population_66 <- population_66_brute |>
  transmute(
    code_commune = as.character(`Code géographique`),
    population_totale = population_totale_calculee
  )

# ----------------------------------------------------------------------------
# 5. ALLOCATAIRES DU RSA
# ----------------------------------------------------------------------------

# Le RSA est déjà donné pour 1 000 habitants de 15 à 64 ans.
# Il manque une commune, remplacée par la médiane du département.
rsa_66 <- suppressWarnings(
  read_excel(
    chemins$rsa,
    sheet = "Data",
    skip = 4,
    col_types = c("text", "text", "numeric", "numeric"),
    na = c("NA", "nd", "s")
  )
) |>
  transmute(
    code_commune = as.character(codgeo),
    annee = as.integer(an),
    rsa_pour_1000_observe = as.numeric(prsa)
  ) |>
  filter(
    str_starts(code_commune, departement),
    annee == 2023
  ) |>
  select(
    code_commune,
    rsa_pour_1000_observe
  ) |>
  mutate(
    statut_rsa = if_else(
      is.na(rsa_pour_1000_observe),
      "Estimé (médiane départementale)",
      "Observé"
    ),
    rsa_pour_1000 = coalesce(
      rsa_pour_1000_observe,
      median(rsa_pour_1000_observe, na.rm = TRUE)
    )
  )

# ----------------------------------------------------------------------------
# 6. FUSION DES CINQ SOURCES
# ----------------------------------------------------------------------------

# Fusion par code commune
tableau_base <- chomage_66 |>
  left_join(
    revenus_66,
    by = "code_commune"
  ) |>
  left_join(
    familles_66 |>
      select(
        code_commune,
        familles_total,
        familles_monoparentales,
        part_monoparentales
      ),
    by = "code_commune"
  ) |>
  left_join(
    population_66,
    by = "code_commune"
  ) |>
  left_join(
    rsa_66,
    by = "code_commune"
  ) |>
  arrange(code_commune)

# ----------------------------------------------------------------------------
# 7. ESTIMATION DES REVENUS MEDIANS MANQUANTS
# ----------------------------------------------------------------------------

# Modèle utilisé pour compléter les revenus manquants
# log1p réduit l'écart entre les petites communes et les grandes villes
formule_revenu <- niveau_vie_median ~
  taux_chomage +
  part_monoparentales +
  part_actifs_bac +
  part_actifs_bac5_plus +
  log1p(population_totale)

# Communes complètes pour entraîner le modèle
donnees_modele <- tableau_base |>
  filter(
    !is.na(niveau_vie_median),
    !is.na(taux_chomage),
    !is.na(part_monoparentales),
    !is.na(part_actifs_bac),
    !is.na(part_actifs_bac5_plus),
    !is.na(population_totale)
  )

# Validation croisée en 5 groupes
set.seed(graine)

donnees_modele <- donnees_modele |>
  mutate(
    groupe_validation = sample(
      rep(1:5, length.out = n())
    )
  )

predictions_validation <- lapply(
  1:5,
  function(groupe_test) {
    apprentissage <- donnees_modele |>
      filter(groupe_validation != groupe_test)

    test <- donnees_modele |>
      filter(groupe_validation == groupe_test)

    modele_temporaire <- lm(
      formule_revenu,
      data = apprentissage
    )

    test |>
      mutate(
        revenu_predit_validation = predict(
          modele_temporaire,
          newdata = test
        )
      )
  }
) |>
  bind_rows()

# Mesures de la qualité du modèle
qualite_modele <- predictions_validation |>
  summarise(
    MAE = mean(
      abs(niveau_vie_median - revenu_predit_validation),
      na.rm = TRUE
    ),
    RMSE = sqrt(
      mean(
        (niveau_vie_median - revenu_predit_validation)^2,
        na.rm = TRUE
      )
    ),
    R2 = 1 -
      sum(
        (niveau_vie_median - revenu_predit_validation)^2,
        na.rm = TRUE
      ) /
      sum(
        (niveau_vie_median - mean(niveau_vie_median))^2,
        na.rm = TRUE
      )
  )

graphique_validation <- ggplot(
  predictions_validation,
  aes(
    x = niveau_vie_median,
    y = revenu_predit_validation
  )
) +
  geom_abline(
    slope = 1,
    intercept = 0,
    linetype = "dashed",
    color = "#6B7280"
  ) +
  geom_point(
    color = "#2563EB",
    alpha = 0.75
  ) +
  coord_equal() +
  labs(
    title = "Validation du modèle de revenu médian",
    subtitle = paste0(
      "MAE = ", round(qualite_modele$MAE), " EUR ; ",
      "RMSE = ", round(qualite_modele$RMSE), " EUR ; ",
      "R2 = ", round(qualite_modele$R2, 2)
    ),
    x = "Revenu médian observé (EUR)",
    y = "Revenu médian prédit (EUR)"
  ) +
  theme_minimal()

# Modèle final avec toutes les communes connues
modele_revenu_final <- lm(
  formule_revenu,
  data = donnees_modele
)

tableau_final <- tableau_base |>
  mutate(
    statut_revenu = if_else(
      is.na(niveau_vie_median),
      "Estimé",
      "Observé"
    ),
    revenu_predit = NA_real_,
    revenu_min_estime = NA_real_,
    revenu_max_estime = NA_real_
  )

# Lignes pour lesquelles le revenu doit être estimé
lignes_a_predire <- which(
  is.na(tableau_final$niveau_vie_median) &
    complete.cases(
      tableau_final[, c(
        "taux_chomage",
        "part_monoparentales",
        "part_actifs_bac",
        "part_actifs_bac5_plus",
        "population_totale"
      )]
    )
)

if (length(lignes_a_predire) > 0) {
  # Prédiction avec un intervalle à 80 %
  predictions_revenus <- predict(
    modele_revenu_final,
    newdata = tableau_final[lignes_a_predire, ],
    interval = "prediction",
    level = 0.80
  )

  tableau_final$revenu_predit[lignes_a_predire] <-
    predictions_revenus[, "fit"]

  tableau_final$revenu_min_estime[lignes_a_predire] <-
    predictions_revenus[, "lwr"]

  tableau_final$revenu_max_estime[lignes_a_predire] <-
    predictions_revenus[, "upr"]
}

tableau_final <- tableau_final |>
  mutate(
    # On garde le revenu observé quand il existe
    revenu_median_complet = coalesce(
      niveau_vie_median,
      revenu_predit
    )
  )

# ----------------------------------------------------------------------------
# 8. SCORES DE VULNERABILITE
# ----------------------------------------------------------------------------

# Les scores vont de 0 à 100 et comparent les communes entre elles
tableau_final <- tableau_final |>
  mutate(
    score_chomage = 100 * percent_rank(taux_chomage),

    # Le signe moins donne un score élevé aux faibles revenus
    score_revenu = 100 * percent_rank(-revenu_median_complet),

    score_monoparentalite = 100 * percent_rank(part_monoparentales),
    score_rsa = 100 * percent_rank(rsa_pour_1000),
    score_bac = 100 * percent_rank(part_actifs_bac),

    # Même principe pour la part de bac+5
    score_faible_bac5 = 100 * percent_rank(-part_actifs_bac5_plus),

    # Moyenne simple des 6 scores
    score_global = rowMeans(
      pick(
        score_chomage,
        score_revenu,
        score_monoparentalite,
        score_rsa,
        score_bac,
        score_faible_bac5
      ),
      na.rm = TRUE
    )
  )

# ----------------------------------------------------------------------------
# 9. PREPARATION DU CLUSTERING
# ----------------------------------------------------------------------------

# Données utilisées par le K-means
donnees_cluster <- tableau_final |>
  select(
    code_commune,
    commune,
    taux_chomage,
    revenu_median_complet,
    part_monoparentales,
    rsa_pour_1000,
    part_actifs_bac,
    part_actifs_bac5_plus,
    score_global
  ) |>
  drop_na(
    taux_chomage,
    revenu_median_complet,
    part_monoparentales,
    rsa_pour_1000,
    part_actifs_bac,
    part_actifs_bac5_plus
  )

variables_cluster <- donnees_cluster |>
  transmute(
    # Les signes sont inversés pour garder le même sens de lecture
    taux_chomage,
    faible_revenu_median = -revenu_median_complet,
    part_monoparentales,
    rsa_pour_1000,
    part_actifs_bac,
    faible_part_bac5_plus = -part_actifs_bac5_plus
  )

# Standardisation pour comparer les variables malgré leurs unités différentes
variables_standardisees <- scale(variables_cluster)

# ----------------------------------------------------------------------------
# 10. CHOIX DU NOMBRE DE CLUSTERS
# ----------------------------------------------------------------------------

k_max_coude <- min(15, nrow(donnees_cluster) - 1)
k_max_silhouette <- min(10, nrow(donnees_cluster) - 1)

set.seed(graine)

graphique_coude <- fviz_nbclust(
  variables_standardisees,
  FUNcluster = kmeans,
  method = "wss",
  k.max = k_max_coude,
  nstart = 100,
  iter.max = 100
) +
  labs(
    title = "Choix du nombre de clusters",
    subtitle = "Méthode du coude",
    x = "Nombre de clusters (k)",
    y = "WSS : inertie intra-classe"
  ) +
  theme_minimal()

set.seed(graine)

graphique_silhouette <- fviz_nbclust(
  variables_standardisees,
  FUNcluster = kmeans,
  method = "silhouette",
  k.max = k_max_silhouette,
  nstart = 100,
  iter.max = 100
) +
  labs(
    title = "Choix du nombre de clusters",
    subtitle = "Critère de silhouette",
    x = "Nombre de clusters (k)",
    y = "Silhouette moyenne"
  ) +
  theme_minimal()

# Sélection automatique : meilleure silhouette moyenne
valeurs_k <- 2:k_max_silhouette
distance_communes <- dist(variables_standardisees)

silhouettes_moyennes <- vapply(
  valeurs_k,
  function(k) {
    set.seed(graine)

    modele_temporaire <- kmeans(
      variables_standardisees,
      centers = k,
      nstart = 100,
      iter.max = 100
    )

    mean(
      silhouette(
        modele_temporaire$cluster,
        distance_communes
      )[, 3]
    )
  },
  numeric(1)
)

k_retenu <- valeurs_k[which.max(silhouettes_moyennes)]

# ----------------------------------------------------------------------------
# 11. COMPARAISON DES SOLUTIONS K = 2 ET K = 3
# ----------------------------------------------------------------------------

# Fonction utilisée pour tester k = 2 puis k = 3
analyser_kmeans <- function(k) {
  set.seed(graine)

  modele <- kmeans(
    variables_standardisees,
    centers = k,
    nstart = 100,
    iter.max = 100
  )

  donnees_k <- donnees_cluster |>
    mutate(
      cluster = factor(modele$cluster)
    )

  # Moyennes des indicateurs par cluster
  profil <- donnees_k |>
    group_by(cluster) |>
    summarise(
      nombre_communes = n(),
      taux_chomage_moyen = mean(taux_chomage),
      revenu_median_moyen = mean(revenu_median_complet),
      part_monoparentales_moyenne = mean(part_monoparentales),
      rsa_pour_1000_moyen = mean(rsa_pour_1000),
      part_actifs_bac_moyenne = mean(part_actifs_bac),
      part_actifs_bac5_plus_moyenne = mean(part_actifs_bac5_plus),
      score_global_moyen = mean(score_global),
      .groups = "drop"
    ) |>
    arrange(desc(score_global_moyen))

  graphique_clusters <- fviz_cluster(
    modele,
    data = variables_standardisees,
    geom = "point",
    ellipse.type = "convex",
    palette = "jco",
    ggtheme = theme_minimal()
  ) +
    labs(
      title = "Typologie sociale des communes",
      subtitle = paste0("K-means avec k = ", k),
      color = "Cluster"
    )

  # Profil des clusters par rapport à la moyenne départementale
  profils_standardises <- as_tibble(
    variables_standardisees
  ) |>
    mutate(
      cluster = donnees_k$cluster
    ) |>
    group_by(cluster) |>
    summarise(
      across(everything(), mean),
      .groups = "drop"
    ) |>
    pivot_longer(
      cols = -cluster,
      names_to = "indicateur",
      values_to = "valeur_standardisee"
    ) |>
    mutate(
      indicateur = recode(
        indicateur,
        taux_chomage = "Taux de chômage",
        faible_revenu_median = "Faible revenu médian",
        part_monoparentales = "Familles monoparentales",
        rsa_pour_1000 = "Allocataires du RSA",
        part_actifs_bac = "Actifs ayant uniquement le bac",
        faible_part_bac5_plus = "Faible part d'actifs bac+5"
      )
    )

  graphique_profils <- ggplot(
    profils_standardises,
    aes(
      x = indicateur,
      y = valeur_standardisee,
      fill = cluster
    )
  ) +
    geom_hline(
      yintercept = 0,
      color = "#6B7280",
      linetype = "dashed"
    ) +
    geom_col(
      position = position_dodge(width = 0.8),
      width = 0.7
    ) +
    labs(
      title = paste0("Profil moyen des clusters (k = ", k, ")"),
      subtitle = paste(
        "Écart à la moyenne départementale",
        "en unités standardisées"
      ),
      x = NULL,
      y = "Valeur standardisée",
      fill = "Cluster"
    ) +
    theme_minimal() +
    theme(
      axis.text.x = element_text(
        angle = 25,
        hjust = 1
      )
    )

  list(
    modele = modele,
    donnees = donnees_k,
    profil = profil,
    graphique_clusters = graphique_clusters,
    graphique_profils = graphique_profils
  )
}

resultat_k2 <- analyser_kmeans(2)
resultat_k3 <- analyser_kmeans(3)

modele_kmeans_k2 <- resultat_k2$modele
modele_kmeans_k3 <- resultat_k3$modele
profil_clusters_k2 <- resultat_k2$profil
profil_clusters_k3 <- resultat_k3$profil
graphique_clusters_k2 <- resultat_k2$graphique_clusters
graphique_clusters_k3 <- resultat_k3$graphique_clusters
graphique_profils_k2 <- resultat_k2$graphique_profils
graphique_profils_k3 <- resultat_k3$graphique_profils

# ACP utilisée uniquement pour afficher les points en 2 dimensions
projection_acp <- prcomp(
  variables_standardisees,
  center = FALSE,
  scale. = FALSE
)

coordonnees_acp <- as_tibble(
  projection_acp$x[, 1:2, drop = FALSE]
)

variance_acp <- 100 * summary(projection_acp)$importance[
  "Proportion of Variance",
  1:2
]

# Graphique interactif avec le nom de la commune au survol
creer_graphique_clusters_interactif <- function(resultat, k) {
  donnees_interactives <- bind_cols(
    coordonnees_acp,
    resultat$donnees |>
      select(
        code_commune,
        commune,
        taux_chomage,
        revenu_median_complet,
        part_monoparentales,
        rsa_pour_1000,
        part_actifs_bac,
        part_actifs_bac5_plus
      )
  )

  donnees_interactives <- donnees_interactives |>
    left_join(
      tableau_final |>
        select(code_commune, statut_revenu, statut_rsa),
      by = "code_commune"
    ) |>
    mutate(
      cluster = resultat$donnees$cluster,
      texte_survol = paste0(
        "<b>", commune, "</b>",
        "<br>Code : ", code_commune,
        "<br>Cluster : ", cluster,
        "<br>Chômage : ", round(taux_chomage, 1), " %",
        "<br>Revenu médian : ", round(revenu_median_complet), " EUR",
        " (", statut_revenu, ")",
        "<br>Familles monoparentales : ",
        round(part_monoparentales, 1), " %",
        "<br>Allocataires du RSA : ",
        round(rsa_pour_1000, 1), " pour 1 000",
        " (", statut_rsa, ")",
        "<br>Actifs ayant uniquement le bac : ",
        round(part_actifs_bac, 1), " %",
        "<br>Actifs ayant un bac+5 ou plus : ",
        round(part_actifs_bac5_plus, 1), " %"
      )
    )

  plot_ly(
    data = donnees_interactives,
    x = ~PC1,
    y = ~PC2,
    color = ~cluster,
    colors = c("#0072B2", "#E69F00", "#009E73"),
    type = "scatter",
    mode = "markers",
    text = ~texte_survol,
    hoverinfo = "text",
    marker = list(
      size = 9,
      opacity = 0.8,
      line = list(width = 0.5, color = "white")
    )
  ) |>
    layout(
      title = list(
        text = paste0(
          "Typologie sociale interactive des communes (k = ",
          k,
          ")"
        )
      ),
      xaxis = list(
        title = paste0(
          "Dimension 1 (",
          round(variance_acp[1], 1),
          " %)"
        )
      ),
      yaxis = list(
        title = paste0(
          "Dimension 2 (",
          round(variance_acp[2], 1),
          " %)"
        )
      ),
      legend = list(title = list(text = "Cluster")),
      hoverlabel = list(align = "left")
    )
}

graphique_clusters_interactif_k2 <-
  creer_graphique_clusters_interactif(resultat_k2, 2)

graphique_clusters_interactif_k3 <-
  creer_graphique_clusters_interactif(resultat_k3, 3)

# Versions HTML interactives des autres graphiques.
graphique_validation_interactif <- ggplotly(graphique_validation)
graphique_coude_interactif <- ggplotly(graphique_coude)
graphique_silhouette_interactif <- ggplotly(graphique_silhouette)
graphique_profils_interactif_k2 <- ggplotly(graphique_profils_k2)
graphique_profils_interactif_k3 <- ggplotly(graphique_profils_k3)

affectations_clusters <- donnees_cluster |>
  transmute(
    code_commune,
    cluster_k2 = resultat_k2$donnees$cluster,
    cluster_k3 = resultat_k3$donnees$cluster
  )

tableau_final <- tableau_final |>
  left_join(
    affectations_clusters,
    by = "code_commune"
  )

# ----------------------------------------------------------------------------
# 12. CARTES INTERACTIVES DES COMMUNES
# ----------------------------------------------------------------------------

geojson_communes <- fromJSON(
  chemin_geojson,
  simplifyVector = FALSE
)

# Echelle de couleurs sans dégradé entre les clusters
creer_echelle_discrete <- function(couleurs) {
  bornes <- seq(
    0,
    1,
    length.out = length(couleurs) + 1
  )

  unlist(
    lapply(
      seq_along(couleurs),
      function(i) {
        list(
          c(bornes[i], couleurs[i]),
          c(bornes[i + 1], couleurs[i])
        )
      }
    ),
    recursive = FALSE
  )
}

# Création d'une carte pour k = 2 ou k = 3
creer_carte_clusters <- function(colonne_cluster, k) {
  couleurs <- c(
    "#0072B2",
    "#E69F00",
    "#009E73"
  )[1:k]

  donnees_carte <- tableau_final |>
    filter(!is.na(.data[[colonne_cluster]])) |>
    mutate(
      cluster_carte = as.integer(.data[[colonne_cluster]]),
      texte_survol = paste0(
        "<b>", commune, "</b>",
        "<br>Code : ", code_commune,
        "<br>Cluster : ", cluster_carte,
        "<br>Score global : ", round(score_global, 1),
        "<br>Chômage : ", round(taux_chomage, 1), " %",
        "<br>Revenu médian : ", round(revenu_median_complet), " EUR",
        " (", statut_revenu, ")",
        "<br>Familles monoparentales : ",
        round(part_monoparentales, 1), " %",
        "<br>Allocataires du RSA : ",
        round(rsa_pour_1000, 1), " pour 1 000",
        " (", statut_rsa, ")",
        "<br>Actifs ayant uniquement le bac : ",
        round(part_actifs_bac, 1), " %",
        "<br>Actifs ayant un bac+5 ou plus : ",
        round(part_actifs_bac5_plus, 1), " %"
      )
    )

  plot_ly(
    data = donnees_carte,
    type = "choroplethmapbox",
    geojson = geojson_communes,
    locations = ~code_commune,
    z = ~cluster_carte,
    featureidkey = "properties.code",
    zmin = 0.5,
    zmax = k + 0.5,
    colorscale = creer_echelle_discrete(couleurs),
    text = ~texte_survol,
    hoverinfo = "text",
    marker = list(
      line = list(
        color = "white",
        width = 0.7
      )
    ),
    colorbar = list(
      title = "Cluster",
      tickmode = "array",
      tickvals = 1:k,
      ticktext = paste("Cluster", 1:k),
      len = 0.45
    )
  ) |>
    layout(
      title = list(
        text = paste0(
          "Typologie sociale des communes des Pyrénées-Orientales (k = ",
          k,
          ")"
        )
      ),
      mapbox = list(
        style = "carto-positron",
        center = list(
          lon = 2.45,
          lat = 42.62
        ),
        zoom = 7.6
      ),
      margin = list(
        l = 0,
        r = 0,
        t = 60,
        b = 0
      )
    )
}

carte_clusters_k2 <- creer_carte_clusters("cluster_k2", 2)
carte_clusters_k3 <- creer_carte_clusters("cluster_k3", 3)

# ----------------------------------------------------------------------------
# 13. TABLEAU LISIBLE ET EXPORTS
# ----------------------------------------------------------------------------

# Tableau final classé par score
tableau_lisible <- tableau_final |>
  arrange(desc(score_global)) |>
  transmute(
    Code = code_commune,
    Commune = commune,
    `Chômage (%)` = round(taux_chomage, 1),
    `Revenu médian observé (€)` = round(niveau_vie_median),
    `Revenu médian utilisé (€)` = round(revenu_median_complet),
    `Statut du revenu` = statut_revenu,
    `Familles monoparentales (%)` = round(part_monoparentales, 1),
    `Allocataires RSA (pour 1 000 habitants de 15-64 ans)` =
      round(rsa_pour_1000, 1),
    `Statut du RSA` = statut_rsa,
    `Actifs ayant uniquement le bac (%)` = round(part_actifs_bac, 1),
    `Actifs ayant un bac+5 ou plus (%)` =
      round(part_actifs_bac5_plus, 1),
    `Score global` = round(score_global, 1),
    `Cluster (k = 2)` = cluster_k2,
    `Cluster (k = 3)` = cluster_k3
  )

# Suppression des anciens résultats avant de recréer les nouveaux
anciennes_sorties <- list.files(
  dossier_resultats,
  pattern = paste0(
    "^(0[1-9]_.*|",
    "diagnostic_communes_66\\.csv$|",
    "profil_clusters.*\\.csv$|",
    "qualite_modele_revenu\\.csv$)"
  ),
  full.names = TRUE,
  all.files = FALSE
)

if (length(anciennes_sorties) > 0) {
  unlink(
    anciennes_sorties,
    recursive = TRUE,
    force = TRUE
  )
}

# Export des tableaux
write.csv2(
  tableau_lisible,
  file.path(dossier_resultats, "diagnostic_communes_66.csv"),
  row.names = FALSE,
  fileEncoding = "UTF-8"
)

write.csv2(
  profil_clusters_k2,
  file.path(dossier_resultats, "profil_clusters_k2.csv"),
  row.names = FALSE,
  fileEncoding = "UTF-8"
)

write.csv2(
  profil_clusters_k3,
  file.path(dossier_resultats, "profil_clusters_k3.csv"),
  row.names = FALSE,
  fileEncoding = "UTF-8"
)

write.csv2(
  qualite_modele,
  file.path(dossier_resultats, "qualite_modele_revenu.csv"),
  row.names = FALSE,
  fileEncoding = "UTF-8"
)

# Export des graphiques PNG
ggsave(
  file.path(dossier_resultats, "01_validation_revenu.png"),
  graphique_validation,
  width = 8,
  height = 6,
  dpi = 300
)

ggsave(
  file.path(dossier_resultats, "02_methode_coude.png"),
  graphique_coude,
  width = 8,
  height = 6,
  dpi = 300
)

ggsave(
  file.path(dossier_resultats, "03_silhouette.png"),
  graphique_silhouette,
  width = 8,
  height = 6,
  dpi = 300
)

ggsave(
  file.path(dossier_resultats, "04_clusters_kmeans_k2.png"),
  graphique_clusters_k2,
  width = 8,
  height = 6,
  dpi = 300
)

ggsave(
  file.path(dossier_resultats, "05_profils_clusters_k2.png"),
  graphique_profils_k2,
  width = 10,
  height = 6,
  dpi = 300
)

ggsave(
  file.path(dossier_resultats, "06_clusters_kmeans_k3.png"),
  graphique_clusters_k3,
  width = 8,
  height = 6,
  dpi = 300
)

ggsave(
  file.path(dossier_resultats, "07_profils_clusters_k3.png"),
  graphique_profils_k3,
  width = 10,
  height = 6,
  dpi = 300
)

# Export des graphiques interactifs
saveWidget(
  graphique_validation_interactif,
  file.path(dossier_resultats, "01_validation_revenu_interactif.html"),
  selfcontained = FALSE
)

saveWidget(
  graphique_coude_interactif,
  file.path(dossier_resultats, "02_methode_coude_interactive.html"),
  selfcontained = FALSE
)

saveWidget(
  graphique_silhouette_interactif,
  file.path(dossier_resultats, "03_silhouette_interactive.html"),
  selfcontained = FALSE
)

saveWidget(
  graphique_clusters_interactif_k2,
  file.path(dossier_resultats, "04_clusters_kmeans_k2_interactif.html"),
  selfcontained = FALSE
)

saveWidget(
  graphique_profils_interactif_k2,
  file.path(dossier_resultats, "05_profils_clusters_k2_interactif.html"),
  selfcontained = FALSE
)

saveWidget(
  graphique_clusters_interactif_k3,
  file.path(dossier_resultats, "06_clusters_kmeans_k3_interactif.html"),
  selfcontained = FALSE
)

saveWidget(
  graphique_profils_interactif_k3,
  file.path(dossier_resultats, "07_profils_clusters_k3_interactif.html"),
  selfcontained = FALSE
)

saveWidget(
  carte_clusters_k2,
  file.path(dossier_resultats, "08_carte_clusters_k2_interactive.html"),
  selfcontained = FALSE
)

saveWidget(
  carte_clusters_k3,
  file.path(dossier_resultats, "09_carte_clusters_k3_interactive.html"),
  selfcontained = FALSE
)

# Affichage dans RStudio.
if (interactive()) {
  print(graphique_validation)
  print(graphique_coude)
  print(graphique_silhouette)
  print(graphique_clusters_k2)
  print(graphique_profils_k2)
  print(graphique_clusters_k3)
  print(graphique_profils_k3)
  print(graphique_clusters_interactif_k2)
  print(graphique_clusters_interactif_k3)
  print(carte_clusters_k2)
  print(carte_clusters_k3)
  View(tableau_lisible)
}