Modelo FIA Optimización

Modelo FIA de optimización

Este documento presenta la fase final del proyecto FIA Speed Camera Enforcement, centrada en la optimización de la localización de nuevas cámaras de fotodetección en Bogotá. Partiendo de un conjunto de modelos predictivos de Machine Learning previamente entrenados, se emplea una simulación de Montecarlo a gran escala para evaluar miles de escenarios de intervención y, finalmente, identificar un conjunto de soluciones óptimas que maximicen la seguridad vial en la red.

El análisis se divide en dos grandes fases metodológicas:

• Modelo de Montecarlo: Un enfoque de simulación que nos permite explorar el vasto espacio de posibles configuraciones de la red y entender la distribución de los resultados esperados.

• Modelo U-NSGA-III: Una propuesta a futuro que utiliza algoritmos genéticos multiobjetivo para encontrar soluciones óptimas de manera más directa, basándose en las funciones de predicción calibradas en la fase anterior.

El punto de partida de este análisis es un conjunto de cuatro modelos de Random Forest (RF), entrenados para predecir indicadores clave de la operación y seguridad vial en los tramos de la ciudad. Estos modelos, cuya construcción se detalló en el informe Modelo RF FIA, fueron diseñados para estimar:

1. Velocidad promedio: La velocidad ponderada por la longitud del tramo en toda la red, un indicador de eficiencia.
2. Índice de exceso: La proporción de tiempo que, en promedio, los vehículos de la red circulan por encima del límite de velocidad. Es un indicador clave de incumplimiento.
3. Índice de riesgo: Un indicador ponderado de la siniestralidad, que da más peso a las fatalidades que a los heridos. Mide la peligrosidad de la red.
4. Probabilidad de siniestro: La probabilidad promedio de que ocurra un siniestro en un tramo de la red.

Estos modelos fueron entrenados utilizando una rica base de datos que incluye variables geométricas, operacionales, de entorno y de seguridad. Mediante un proceso de Recursive Feature Elimination (RFE), se seleccionó un subconjunto óptimo de predictores para cada modelo, garantizando su parsimonia y poder predictivo. Estos modelos entrenados y guardados en el archivo modelos_entrenados.RData.

Modelo Montecarlo

La simulación de Montecarlo es una técnica computacional que nos permite modelar la probabilidad de diferentes resultados en un proceso que no puede predecirse fácilmente debido a la intervención de variables aleatorias. En nuestro caso, la “variable aleatoria” es la decisión de dónde ubicar 100 nuevas cámaras de fotodetección.

Generación de escenarios en Python

El primer paso es crear los “universos posibles” que queremos evaluar. Para ello, se desarrolló un script en Python (Montecarlo.py) que realiza las siguientes tareas de manera altamente eficiente:

• Pre-carga de Datos Geoespaciales: Se carga el shapefile de la red vial de Bogotá una sola vez en memoria.

• Identificación de Candidatos: Se identifican todos los tramos de vía que actualmente no cuentan con una cámara, convirtiéndolos en candidatos para una nueva instalación.

• Generación de Escenarios: Se ejecuta un bucle para cada una de las 6,000 simulaciones. En cada iteración:

1. Se asigna un número único de identificación (seed) para garantizar la reproducibilidad.

2. Se seleccionan al azar 100 tramos viales de la lista de candidatos.

3. Se simula la “instalación” de estas cámaras y, utilizando algoritmos geoespaciales optimizados (cKDTree), se recalcula la variable dist_min_camara (distancia a la cámara más cercana) para todos los tramos de la red.

• Paralelización: Este proceso, que sería computacionalmente prohibitivo si se hiciera de forma secuencial, se paraleliza utilizando la librería multiprocessing de Python. Esto permite generar miles de escenarios en cuestión de minutos, aprovechando todos los núcleos de la CPU.

• Salida: El resultado de esta fase es un único archivo batch_scenarios_py.feather que contiene los datos de los 11,561 tramos viales replicados 6,000 veces, una por cada escenario, con los valores de seed y dist_min correspondientes a cada simulación.

 # -*- coding: utf-8 -*-
"""
@author: Universidad de los Andes - Dayling Páez
Adaptado para Simulación de Montecarlo (Versión Optimizada y Paralelizada)

==============================================================================
SCRIPT: simulacion_montecarlo_optimizado.py

DESCRIPCIÓN:
Versión final que implementa procesamiento por lotes para manejar grandes
simulaciones sin agotar la memoria, e incluye un diagnóstico de outliers
para un análisis de resultados más robusto.
==============================================================================
"""

import os
import subprocess
import pandas as pd
import numpy as np
import geopandas as gpd
from pyproj import CRS
from scipy.spatial import cKDTree
import matplotlib.pyplot as plt
import seaborn as sns
from tqdm import tqdm
import multiprocessing
from functools import partial
import time

# ============================================================
# --- CONFIGURACIÓN DE LA SIMULACIÓN Y RUTAS ---
# ============================================================

N_SIMULACIONES = 6000 # Objetivo final
N_NUEVAS_CAMARAS = 100 

R_EXECUTABLE_PATH = r"C:\Program Files\R\R-4.4.2\bin\Rscript.exe"  
BASE_WORK_DIR = r"C:\Users\tomas\Documents\UniAndes\FIA\Modelo\Modelo FIA Enforcement"
SHAPEFILE_DIR = r"C:\users\tomas\Documents\UniAndes\FIA\Modelo\Modelo FIA Enforcement\shp"

R_SCRIPT_NAME = "optim_optparse_batch.R"
VIAS_FILE_NAME = "Tabla_corregida_poligonos.xlsx"
MODEL_FILE_NAME = "modelos_entrenados.RData"
SHAPEFILE_NAME = "capa_modelo_opti.shp"

R_SCRIPT_PATH = os.path.join(BASE_WORK_DIR, R_SCRIPT_NAME)
VIAS_FILE_PATH = os.path.join(BASE_WORK_DIR, VIAS_FILE_NAME)
MODEL_FILE_PATH = os.path.join(BASE_WORK_DIR, MODEL_FILE_NAME)
POLY_PATH = os.path.join(SHAPEFILE_DIR, SHAPEFILE_NAME)

FIELD_ID = "PK_ID_CALZ"
FIELD_VIA = "TIPOS_VIA"
FIELD_SENT = "MVISVia"
FIELD_TIENE_CAM = "tiene_cama"
FIELD_LONG = "Longitudtr"

# ... (Todas las funciones auxiliares se mantienen igual, las he omitido por brevedad) ...
# (Asegúrate de mantener tus funciones init_worker, worker_function_for_pool, etc.)
def init_worker(gdf_from_main, n_cams_from_main):
    global worker_gdf_base, worker_n_cams
    worker_gdf_base = gdf_from_main
    worker_n_cams = n_cams_from_main

def worker_function_for_pool(seed):
    return generar_escenario_df(gdf_base=worker_gdf_base, 
                                n_new_cams=worker_n_cams, 
                                random_seed=seed)

def opposite(code: str):
    return {"FT": "TF", "TF": "FT"}.get(code, code)

def ensure_meters(gdf, epsg=3116):
    if gdf.crs is None: raise ValueError("⚠️ La capa no tiene CRS definido.")
    if not gdf.crs.is_projected: return gdf.to_crs(epsg=epsg)
    return gdf

def _to_num(df, col, default, dtype):
    if col not in df.columns: df[col] = default
    df[col] = pd.to_numeric(df[col], errors="coerce").fillna(default).astype(dtype)

def calcular_distancias_optimizadas(tramos_gdf, cams_gdf):
    if cams_gdf.empty: return np.full(len(tramos_gdf), np.inf)
    tramos_coords = np.array(list(zip(tramos_gdf.geometry.x, tramos_gdf.geometry.y)))
    cams_coords = np.array(list(zip(cams_gdf.geometry.x, cams_gdf.geometry.y)))
    cam_tree = cKDTree(cams_coords)
    k_vecinos = min(5, len(cams_gdf))
    distancias, indices = cam_tree.query(tramos_coords, k=k_vecinos, workers=-1)
    if distancias.ndim == 1:
        distancias, indices = distancias.reshape(-1, 1), indices.reshape(-1, 1)
    tramo_vias = tramos_gdf[FIELD_VIA].values[:, np.newaxis]
    tramo_sents_opuestos = np.vectorize(opposite)(tramos_gdf["SENT"].values)[:, np.newaxis]
    vecino_vias = cams_gdf[FIELD_VIA].values[indices]
    vecino_sents = cams_gdf["SENT"].values[indices]
    es_mismo_corredor_mat = (vecino_vias == tramo_vias)
    es_sentido_opuesto_mat = (vecino_sents == tramo_sents_opuestos)
    es_invalida_mat = es_mismo_corredor_mat & es_sentido_opuesto_mat
    distancias_validas = np.where(es_invalida_mat, np.inf, distancias)
    return np.min(distancias_validas, axis=1)

def generar_escenario_df(gdf_base, n_new_cams=100, random_seed=None):
    cent = gdf_base.copy()
    eligible = cent.loc[cent["HAS_BASE"] == 0].copy()
    n_take  = min(n_new_cams, len(eligible))
    rng = np.random.default_rng(random_seed)
    picked_idx = rng.choice(eligible.index.values, size=n_take, replace=False) if n_take > 0 else []
    cent["HAS_SCEN"] = cent["HAS_BASE"]
    cent.loc[picked_idx, "HAS_SCEN"] = 1
    cams_scen = cent.loc[cent["HAS_SCEN"] == 1].copy()
    cent["DIST_MIN"] = calcular_distancias_optimizadas(cent, cams_scen)
    out = cent[[FIELD_ID, "DIST_MIN", "HAS_SCEN"]].copy()
    out.rename(columns={FIELD_ID: "pk_id_calz", "HAS_SCEN": "hay_camara"}, inplace=True)
    out["DIST_MIN"] = out["DIST_MIN"].round(2)
    out['seed'] = random_seed
    return out

# ==========================================================================
# --- BLOQUE PRINCIPAL (VERSIÓN FINAL CON PROCESAMIENTO POR LOTES Y DIAGNÓSTICO) ---
# ==========================================================================
if __name__ == "__main__":
    
    start_time_total = time.time()
    print("--- Iniciando Simulación de Montecarlo (VERSIÓN BATCH POR LOTES) ---")

    CHUNK_SIZE = 500 # <-- Tu elección de tamaño de lote
    print(f"Configuración: Se procesarán lotes de {CHUNK_SIZE} simulaciones.")

    # --- 1. PRE-CARGA DE DATOS (Sin cambios) ---
    print("\nPaso 1: Cargando y preparando datos base...")
    gdf = gpd.read_file(POLY_PATH)
    gdf = ensure_meters(gdf)
    gdf_base = gdf.copy()
    gdf_base["geometry"] = gdf_base.geometry.centroid
    _to_num(gdf_base, FIELD_TIENE_CAM, 0, int)
    _to_num(gdf_base, FIELD_LONG, 0, float)
    gdf_base["SENT"] = gdf_base[FIELD_SENT].astype(str).str.upper().fillna("B")
    gdf_base.loc[~gdf_base["SENT"].isin(["FT","TF","B","N","SD"]), "SENT"] = "B"
    gdf_base["HAS_BASE"] = (gdf_base[FIELD_TIENE_CAM] > 0).astype(int)
    print("Datos base listos.")

    # --- 2. GENERACIÓN DE TODOS LOS ESCENARIOS (Sin cambios) ---
    print(f"\nPaso 2: Generando {N_SIMULACIONES} escenarios en memoria (en paralelo)...")
    num_procesos = max(1, multiprocessing.cpu_count() - 1)
    iterable_seeds = range(N_SIMULACIONES)
    lista_escenarios_df = []
    with multiprocessing.Pool(processes=num_procesos, initializer=init_worker, initargs=(gdf_base, N_NUEVAS_CAMARAS)) as pool:
        for escenario_df in tqdm(pool.imap_unordered(worker_function_for_pool, iterable_seeds), total=N_SIMULACIONES):
            if escenario_df is not None:
                lista_escenarios_df.append(escenario_df)
    
    batch_df = pd.concat(lista_escenarios_df, ignore_index=True)
    print("Generación de escenarios completada.")

    # --- 3. PROCESAMIENTO POR LOTES CON R ---
    print("\nPaso 3: Procesando escenarios en R por lotes para optimizar memoria...")
    
    ruta_salida_r_batch = os.path.join(BASE_WORK_DIR, "batch_resultados_r.csv")
    
    if os.path.exists(ruta_salida_r_batch):
        os.remove(ruta_salida_r_batch)

    num_chunks = int(np.ceil(N_SIMULACIONES / CHUNK_SIZE))
    all_chunks_successful = True
    
    for i in range(num_chunks):
        print(f"\n--- Procesando Lote {i+1} de {num_chunks} ---")
        
        start_seed = i * CHUNK_SIZE
        end_seed = min((i + 1) * CHUNK_SIZE, N_SIMULACIONES)
        
        chunk_df = batch_df[batch_df['seed'].between(start_seed, end_seed - 1)]
        
        ruta_chunk_py = os.path.join(BASE_WORK_DIR, "temp_chunk_py.feather")
        chunk_df.to_feather(ruta_chunk_py)
        
        command = [
            R_EXECUTABLE_PATH, R_SCRIPT_PATH,
            "--vector_file", ruta_chunk_py,
            "--vias_file", VIAS_FILE_PATH,
            "--model_file", MODEL_FILE_PATH,
            "--output_file", ruta_salida_r_batch
        ]
        
        try:
            result = subprocess.run(command, check=True, capture_output=True, text=True, encoding='utf-8')
            print(f"Lote {i+1} procesado exitosamente por R.")
            print("Salida de R (stdout):")
            print(result.stdout)
        
        except subprocess.CalledProcessError as e:
            print(f"\n❌ ERROR: El script de R falló durante el procesamiento del lote {i+1}.")
            print(e.stderr)
            all_chunks_successful = False
            break
            
        finally:
            if os.path.exists(ruta_chunk_py):
                os.remove(ruta_chunk_py)

    print("\n--- Todos los lotes han sido procesados ---")

    # ==============================================================================
    # --- PASO 4: ANÁLISIS, DIAGNÓSTICO, PARETO Y VISUALIZACIÓN DE RESULTADOS ---
    # ==============================================================================
    if all_chunks_successful and os.path.exists(ruta_salida_r_batch):
        print(f"\nPaso 4: Leyendo, diagnosticando y visualizando resultados desde: {ruta_salida_r_batch}")
        resultados_finales = pd.read_csv(ruta_salida_r_batch)
        
        # --- 4.1 Diagnóstico de Outliers ---
        print("\n--- Diagnóstico de Outliers ---")
        umbral_outlier = resultados_finales['velocidad_promedio_red'].quantile(0.99)
        outliers = resultados_finales[resultados_finales['velocidad_promedio_red'] > umbral_outlier]

        if not outliers.empty:
            print(f"Se han identificado {len(outliers)} outlier(s) por encima del percentil 99 ({umbral_outlier:.4f}).")
            print("Detalles del outlier principal (por velocidad):")
            print(outliers.sort_values('velocidad_promedio_red', ascending=False).head(1))
            
            print("\nCreando un set de datos filtrado para análisis y gráficos...")
            resultados_filtrados = resultados_finales[resultados_finales['velocidad_promedio_red'] <= umbral_outlier]
        else:
            print("No se encontraron outliers significativos.")
            resultados_filtrados = resultados_finales
        print("---------------------------------")

        # --- 4.2 Resumen Estadístico Completo (Filtrado) ---
        print("\n--- Resumen Estadístico (Filtrado) de las 4 Funciones Objetivo ---")
        with pd.option_context('display.max_columns', None, 'display.width', 1000):
            print(resultados_filtrados.drop(columns=['seed']).describe().round(4))

        # --- 4.3 Generación de Gráficas (con datos filtrados) ---
        print("\n--- Generando Gráficas de Resultados (con datos filtrados) ---")
        
        COL_VEL = 'velocidad_promedio_red'
        COL_EXC = 'indice_exceso_red'
        COL_RSK = 'indice_riesgo_red'
        COL_SIN = 'probabilidad_siniestro_red'
        num_escenarios_filtrados = len(resultados_filtrados)
        sns.set_theme(style="whitegrid")

        # Gráfico 1: Distribución de las 4 variables
        fig, axes = plt.subplots(2, 2, figsize=(18, 14))
        fig.suptitle(f'Distribución de Resultados en {num_escenarios_filtrados} Escenarios (Filtrados)', fontsize=20)
        
        sns.violinplot(ax=axes[0, 0], y=resultados_filtrados[COL_VEL], color="lightblue").set_title('Velocidad Promedio')
        sns.violinplot(ax=axes[0, 1], y=resultados_filtrados[COL_EXC], color="lightcoral").set_title('Índice de Exceso')
        sns.violinplot(ax=axes[1, 0], y=resultados_filtrados[COL_RSK], color="lightsalmon").set_title('Índice de Riesgo')
        sns.violinplot(ax=axes[1, 1], y=resultados_filtrados[COL_SIN], color="lightgreen").set_title('Prob. de Siniestro')
        plt.tight_layout(rect=[0, 0.03, 1, 0.95])
        plt.savefig(os.path.join(BASE_WORK_DIR, "grafica_distribucion_todas_filtrada.png"), dpi=300)
        plt.show()

        # Gráfico 2: Relaciones cruzadas
        g = sns.pairplot(
            data=resultados_filtrados,
            vars=[COL_VEL, COL_EXC, COL_RSK, COL_SIN],
            kind='reg',
            plot_kws={'scatter_kws': {'alpha': 0.1, 's': 20}},
            diag_kind='kde'
        )
        g.fig.suptitle("Relaciones Cruzadas entre Variables Objetivo (Filtrado)", y=1.02, fontsize=18)
        plt.savefig(os.path.join(BASE_WORK_DIR, "grafica_relaciones_cruzadas_filtrada.png"), dpi=300)
        plt.show()

        # --- 4.4 Análisis de Frente de Pareto ---
        print("\n--- Análisis de Frente de Pareto (3 Objetivos de Minimización) ---")

        pareto_data = resultados_filtrados.copy()

        # Los 3 objetivos a MINIMIZAR son: Exceso, Riesgo y Siniestralidad
        objetivos = ['indice_exceso_red', 'indice_riesgo_red', 'probabilidad_siniestro_red']
        
        is_pareto = np.ones(pareto_data.shape[0], dtype=bool)
        for i, row in pareto_data.iterrows():
            # Compara la fila 'i' con todas las demás
            # Si cualquier otra fila es mejor o igual en todos los objetivos Y estrictamente mejor en al menos uno...
            if np.any(np.all(pareto_data[objetivos] <= row[objetivos], axis=1) & 
                      np.any(pareto_data[objetivos] < row[objetivos], axis=1)):
                # ... entonces la fila 'i' está dominada.
                is_pareto[i] = False

        frente_pareto = resultados_filtrados.loc[is_pareto].copy()

        print(f"Se encontraron {len(frente_pareto)} escenarios en el Frente de Pareto (soluciones óptimas de seguridad).")
        
        if not frente_pareto.empty:
            print("\nMostrando los escenarios del Frente de Pareto con mejor rendimiento individual:")
            
            # Escenario con el MENOR índice de exceso
            mejor_exceso = frente_pareto.loc[frente_pareto[COL_EXC].idxmin()]
            print(f"\n--- Mejor para ÍNDICE DE EXCESO (Seed: {mejor_exceso['seed']:.0f}) ---")
            print(mejor_exceso)

            # Escenario con el MENOR índice de riesgo
            mejor_riesgo = frente_pareto.loc[frente_pareto[COL_RSK].idxmin()]
            print(f"\n--- Mejor para ÍNDICE DE RIESGO (Seed: {mejor_riesgo['seed']:.0f}) ---")
            print(mejor_riesgo)
            
            # Escenario con la MENOR probabilidad de siniestro
            mejor_siniestro = frente_pareto.loc[frente_pareto[COL_SIN].idxmin()]
            print(f"\n--- Mejor para PROB. DE SINIESTRO (Seed: {mejor_siniestro['seed']:.0f}) ---")
            print(mejor_siniestro)
            
            # Escenario "Balanceado": el que está más cerca del "punto utópico" (0,0,0)
            # Normalizamos los datos para que la distancia sea justa
            frente_pareto_norm = frente_pareto.copy()
            for col in objetivos:
                frente_pareto_norm[col] = (frente_pareto[col] - frente_pareto[col].min()) / (frente_pareto[col].max() - frente_pareto[col].min())
            
            frente_pareto_norm['distancia'] = np.sqrt(frente_pareto_norm[objetivos[0]]**2 + frente_pareto_norm[objetivos[1]]**2 + frente_pareto_norm[objetivos[2]]**2)
            mejor_balanceado_idx = frente_pareto_norm['distancia'].idxmin()
            mejor_balanceado = frente_pareto.loc[mejor_balanceado_idx]
            
            print(f"\n--- Escenario más BALANCEADO (Seed: {mejor_balanceado['seed']:.0f}) ---")
            print(mejor_balanceado)
            
    else:
        print("❌ No se generó el archivo de resultados. La simulación pudo haber fallado.")
        
    # --- Limpieza Final ---
    #if os.path.exists(ruta_salida_r_batch):
     #   os.remove(ruta_salida_r_batch)
        
    end_time_total = time.time()
    total_seconds = end_time_total - start_time_total
    mins, secs = divmod(total_seconds, 60)
    
    print("\n--- Proceso Finalizado ---")
    print(f"⏱️ Tiempo de ejecución total: {int(mins)} minutos y {secs:.2f} segundos.")

Simulación con modelos entrenados en R

Una vez generados los escenarios, la siguiente fase consiste en ejecutar nuestros modelos de Random Forest sobre ellos para predecir los resultados. Esta tarea se orquesta mediante un script de R (optim_optparse_batch.R) que está diseñado para ser llamado desde Python y procesar los datos por lotes, evitando así problemas de memoria.

El script de R realiza los siguientes pasos para cada lote de simulaciones:

• Carga de datos: Carga los datos maestros de las vías, los modelos de RF entrenados desde modelos_entrenados.RData, y el lote de escenarios generado por Python.

• Replicación del pre-procesamiento: Es el paso más crítico. El script replica exactamente la misma secuencia de ingeniería de variables que se utilizó para entrenar los modelos originales (documentada en el informe Modelo RF FIA). Esto incluye la limpieza de nombres con janitor, la conversión de tipos de datos, la creación de variables como indice_exceso y, fundamentalmente, la conversión de variables categóricas a factores con los mismos niveles de referencia (forcats::fct_relevel). Esto garantiza que los datos que se pasan a la función predict() tienen la estructura y formato exactos que el modelo espera.

• Optimización por deduplicación: Antes de predecir, el script implementa una optimización crucial. En lugar de ejecutar los modelos sobre el DataFrame masivo (ej. 500 simulaciones * 11,561 vías), primero identifica las filas “únicas” basadas en el conjunto de todas las variables predictoras requeridas por los cuatro modelos. Esto reduce drásticamente el número de predicciones necesarias (a menudo en más de un 90%), convirtiendo un proceso de horas en uno de minutos.

• Predicción y agregación: Se ejecutan los cuatro modelos de RF sobre el DataFrame deduplicado. Luego, los resultados de las predicciones se unen de vuelta al DataFrame completo. Finalmente, se agrupan los resultados por seed y se calcula el promedio ponderado por longitud de tramo (longitudtr) para cada una de las cuatro variables objetivo.

• Salida: Los resultados agregados de cada lote se añaden a un único archivo CSV (batch_resultados_r.csv).

#!/usr/bin/env Rscript

# ==============================================================================
# SCRIPT: optim_optparse_batch.R (VERSIÓN FINAL, CORRECCIÓN DE LOGGING)
# FECHA: 2025-10-10
#
# DESCRIPCIÓN:
# Versión final y funcional que corrige un error menor en una línea de
# impresión (logging). La lógica principal es correcta.
# ==============================================================================

# --- 1. Carga de Librerías ---
suppressPackageStartupMessages(library(optparse))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(readxl))
suppressPackageStartupMessages(library(janitor))
suppressPackageStartupMessages(library(caret))
suppressPackageStartupMessages(library(ranger))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(tictoc))

# --- 2. Definición de Opciones de Línea de Comandos ---
option_list <- list(
  make_option(c("-v", "--vector_file"), action = 'store', default = NULL, type = 'character', help = "Ruta al archivo .feather del lote de escenarios."),
  make_option(c("-d", "--vias_file"), action = 'store', default = NULL, type = 'character', help = "Ruta al archivo .xlsx con la base de datos maestra de tramos viales."),
  make_option(c("-m", "--model_file"), action = 'store', default = NULL, type = 'character', help = "Ruta al archivo .RData con los modelos entrenados."),
  make_option(c("-o", "--output_file"), action = 'store', default = NULL, type = 'character', help = "Ruta del archivo CSV de salida para los resultados.")
)
opt <- parse_args(OptionParser(option_list = option_list))
if (is.null(opt$vector_file) || is.null(opt$vias_file) || is.null(opt$model_file) || is.null(opt$output_file)) {
  stop("Todos los argumentos (...) son obligatorios.", call. = FALSE)
}

# --- 3. Proceso Principal ---
tryCatch({
  
  tictoc::tic("Total R Script")
  
  # --- Carga de Datos ---
  tictoc::tic("1. Carga de datos")
  cat("INFO: Cargando archivos de entrada...\n")
  vector_batch <- arrow::read_feather(opt$vector_file) %>% clean_names()
  vias_maestra <- read_xlsx(opt$vias_file) %>% 
    clean_names() %>%
    mutate(across(where(is.character), ~na_if(.x, "SIN DATO")))
  load(opt$model_file)
  tictoc::toc()
  
  # --- Preparación del Lote Completo (Replicando Lógica de Entrenamiento) ---
  tictoc::tic("2. Join y Mutate (Creación Batch)")
  cat("INFO: [BATCH] Uniendo y preparando el lote completo...\n")
  datos_completos <- vias_maestra %>%
    left_join(vector_batch, by = "pk_id_calz")
  
  datos_procesados <- datos_completos %>%
    mutate(
      dist_min_camara = dist_min,
      indice_exceso = case_when(
        "exceso" %in% names(.) ~ pmin(pmax(as.numeric(exceso), 0), 1),
        "hora_mas_5" %in% names(.) ~ pmin(pmax(as.numeric(hora_mas_5)/24, 0), 1),
        TRUE ~ NA_real_
      ),
      n_muer = as.integer(n_muer), n_her = as.integer(n_her), estrato_fi = as.integer(estrato_fi),
      estaashto = if_else(is.na(estaashto), "BUENO", estaashto),
      tipomallae = forcats::fct_relevel(as.factor(tipomallae), "Arterial"),
      estaashto = forcats::fct_relevel(as.factor(estaashto), "FALLADO"),
      aashto_bue = forcats::fct_relevel(as.factor(aashto_bue), "0"),
      tipo = forcats::fct_relevel(as.factor(tipo), "Sin Construir"),
      flujo_1 = forcats::fct_relevel(as.factor(flujo_1), "BAJO"),
      pres_anden = forcats::fct_relevel(as.factor(pres_anden), "0"),
      ciclorruta = forcats::fct_relevel(as.factor(ciclorruta), "0"),
      sitp = forcats::fct_relevel(as.factor(sitp), "0"),
      ancho_carril = ifelse(carriles > 0, ancho / carriles, 0)
    )
  rm(datos_completos, vias_maestra, vector_batch); gc()
  tictoc::toc()
  
  # --- Optimización: Deduplicación para Predicciones ---
  tictoc::tic("3. Deduplicación de Datos")
  cat("INFO: Optimizando: extrayendo predictores de perfiles RFE y deduplicando...\n")
  
  preds_vel  <- predictors(rfe_profile_vel)
  preds_exc  <- predictors(rfe_profile_exc)
  preds_risk <- predictors(rfe_profile_risk)
  preds_log  <- predictors(rfe_profile_log)
  
  predictor_vars <- unique(c(preds_vel, preds_exc, preds_risk, preds_log))
  cat(paste("INFO: Se identificaron", length(predictor_vars), "predictores únicos en total en los 4 modelos.\n"))
  
  datos_unicos <- datos_procesados %>%
    distinct(across(all_of(predictor_vars)))
  
  # --- CORRECCIÓN DE LA LÍNEA DE LOGGING ---
  cat(paste0("INFO: Reducción de datos para predicción: de ", nrow(datos_procesados), " a ", nrow(datos_unicos), " filas únicas.\n"))
  tictoc::toc()
  
  # --- Predicciones sobre Dataset Único ---
  cat("INFO: [BATCH] Ejecutando predicciones optimizadas...\n")
  
  tictoc::tic("4.1 Predicción Velocidad"); datos_unicos$vel_promed_pred <- predict(rf_vel_final, newdata = datos_unicos, na.action = na.pass); tictoc::toc()
  tictoc::tic("4.2 Predicción Exceso"); datos_unicos$indice_exceso_pred <- predict(rf_exc_final, newdata = datos_unicos, na.action = na.pass); datos_unicos$indice_exceso_pred <- pmin(pmax(datos_unicos$indice_exceso_pred, 0), 1); tictoc::toc()
  tictoc::tic("4.3 Predicción Riesgo"); datos_unicos$indice_riesgo_pred <- predict(rf_risk_final, newdata = datos_unicos, na.action = na.pass); tictoc::toc()
  tictoc::tic("4.4 Predicción Siniestro"); pred_siniestros_prob <- predict(rf_log_final, newdata = datos_unicos, type = "prob", na.action = na.pass); datos_unicos$prob_siniestro_pred <- pred_siniestros_prob[, "Sí"]; tictoc::toc()
  
  # --- Unir Predicciones al Dataset Completo ---
  tictoc::tic("5. Join de Predicciones")
  cat("INFO: Uniendo los resultados de las predicciones al lote completo...\n")
  
  datos_procesados <- datos_procesados %>%
    left_join(datos_unicos, by = predictor_vars)
  rm(datos_unicos); gc()
  tictoc::toc()
  
  # --- Cálculo de Métricas Agregadas ---
  tictoc::tic("6. Agregación Final")
  cat("INFO: [BATCH] Calculando métricas agregadas por escenario...\n")
  
  resultados_finales <- datos_procesados %>%
    mutate(longitudtr = as.numeric(longitudtr)) %>%
    filter(!is.na(longitudtr) & longitudtr > 0) %>%
    group_by(seed) %>%
    summarise(
      longitud_total_red = sum(longitudtr, na.rm = TRUE),
      velocidad_promedio_red = sum(vel_promed_pred * longitudtr, na.rm = TRUE) / longitud_total_red,
      indice_exceso_red = sum(indice_exceso_pred * longitudtr, na.rm = TRUE) / longitud_total_red,
      indice_riesgo_red = sum(indice_riesgo_pred * longitudtr, na.rm = TRUE) / longitud_total_red,
      probabilidad_siniestro_red = sum(prob_siniestro_pred * longitudtr, na.rm = TRUE) / longitud_total_red,
      .groups = 'drop'
    ) %>%
    select(-longitud_total_red)
  tictoc::toc()
  
  # --- Guardado de Resultados ---
  tictoc::tic("7. Guardado de Resultados")
  cat(paste0("INFO: Guardando resultados consolidados en: ", opt$output_file, "\n"))
  write_csv(resultados_finales, opt$output_file, append = file.exists(opt$output_file))
  tictoc::toc()
  
  cat("\n✅ Proceso en R finalizado exitosamente.\n")
  tictoc::toc(log = TRUE)
  
}, error = function(e) {
  cat(paste0("\n❌ ERROR: Ha ocurrido un error durante la ejecución del script de R.\n",
             "Mensaje de error: ", e$message, "\n"))
  quit(status = 1)
})

Análisis de resultados y frente de Pareto

Con los 6,000 resultados agregados en un solo archivo CSV, pasamos a la fase final de análisis, que se realiza en este mismo documento de R Markdown.

Carga y limpieza de datos

# 1. Cargar librerías

suppressPackageStartupMessages({
  library(tidyverse)
  library(sf)
  library(skimr)
  library(GGally)
  library(plotly)
  library(rgl)
  library(kableExtra)
  library(here) 
  library(plot3D)
  library(rayshader)
})

# 2. Definir rutas a los archivos
ruta_resultados_csv <- "C:/Users/tomas/Documents/UniAndes/FIA/Modelo/Modelo FIA Enforcement/batch_resultados_r.csv"
ruta_shapefile_vias <- "C:/Users/tomas/Documents/UniAndes/FIA/Modelo/Modelo FIA Enforcement/shp/capa_modelo_opti.shp"

# 3. Cargar y limpiar datos

resultados_limpios <- read_csv(ruta_resultados_csv, show_col_types = FALSE) %>%
  filter(!is.na(seed))

# Cargamos los datos espaciales una sola vez aquí para usarlos después.
vias_sf <- st_read(ruta_shapefile_vias, quiet = TRUE) %>%
  janitor::clean_names()

Análisis estadístico y gráficos

Una vez limpios los datos, se realiza un análisis exploratorio para entender la naturaleza de los resultados. Se genera un resumen estadístico completo y dos visualizaciones clave:

1. Gráficos de Violín: Muestran la distribución de la densidad de cada una de las cuatro variables objetivo. Nos permiten ver dónde se concentra la mayoría de los resultados y cuál es el rango de los valores más extremos.

2. Gráfico de Pares (Pair Plot): Es una matriz de gráficos de dispersión que visualiza la relación entre cada par de variables objetivo. Es fundamental para entender los “trade-offs”: por ejemplo, ¿los escenarios que mejoran la velocidad tienden a empeorar el índice de riesgo?

# 1. Resumen estadístico
resultados_limpios %>%
  select(-seed) %>%
  skim() %>%
  print()

## ── Data Summary ────────────────────────
##                            Values    
## Name                       Piped data
## Number of rows             6000      
## Number of columns          4         
## _______________________              
## Column type frequency:               
##   numeric                  4         
## ________________________             
## Group variables            None      
## 
## ── Variable type: numeric ──────────────────────────────────────────────────────
##   skim_variable              n_missing complete_rate   mean       sd     p0
## 1 velocidad_promedio_red             0             1 40.4   0.0298   40.3  
## 2 indice_exceso_red                  0             1  0.179 0.000709  0.177
## 3 indice_riesgo_red                  0             1  4.09  0.00771   4.06 
## 4 probabilidad_siniestro_red         0             1  0.568 0.000415  0.567
##      p25    p50    p75   p100 hist 
## 1 40.4   40.4   40.4   40.6   ▁▃▇▂▁
## 2  0.178  0.179  0.179  0.184 ▁▇▂▁▁
## 3  4.08   4.09   4.09   4.11  ▁▂▇▅▁
## 4  0.568  0.568  0.569  0.570 ▁▃▇▅▁

# 2. Gráficos
# Gráfico 1: Distribuciones (Violin Plots)
resultados_limpios %>%
  pivot_longer(cols = -seed, names_to = "variable", values_to = "valor") %>%
  ggplot(aes(x = variable, y = valor, fill = variable)) +
  geom_violin(trim = FALSE, alpha = 0.8) +
  geom_boxplot(width = 0.1, fill = "white") +
  facet_wrap(~variable, scales = "free") +
  theme_minimal() +
  theme(legend.position = "none", axis.text.x = element_blank()) +
  labs(title = paste("Distribución de resultados en", nrow(resultados_limpios), "Escenarios"),
       x = "Variable objetivo", y = "Valor")

# Gráfico 2: Relaciones Cruzadas (Pair Plot)
GGally::ggpairs(
  resultados_limpios,
  columns = c("velocidad_promedio_red", "indice_exceso_red", "indice_riesgo_red", "probabilidad_siniestro_red"),
  upper = list(continuous = GGally::wrap("cor", size = 4)),
  lower = list(continuous = GGally::wrap("points", alpha = 0.1, size = 0.5)),
  title = "Relaciones cruzadas entre variables objetivo"
)

Identificación de frentes de Pareto

Dado que tenemos múltiples objetivos de seguridad que queremos minimizar (exceso, riesgo y siniestralidad), no existe un único “mejor” escenario, sino un conjunto de soluciones óptimas conocido como el Frente de Pareto. Un escenario pertenece a este frente si es “no dominado”, lo que significa que no existe ningún otro escenario que sea mejor en al menos uno de los objetivos sin ser peor en los demás.

Estos escenarios representan los mejores compromisos posibles que se pueden lograr. Para identificarlos, se implementa una función que compara cada escenario con todos los demás y filtra aquellos que no son dominados.

# Función para determinar si un punto es dominado
is_dominated <- function(punto, todos_los_puntos) {
  # Un punto es dominado si existe CUALQUIER otro punto que es
  # mejor o igual en todos los objetivos Y estrictamente mejor en al menos uno.
   any(
    colSums(t(todos_los_puntos) <= punto) == ncol(todos_los_puntos) &
    colSums(t(todos_los_puntos) < punto) > 0
  )
}

# Objetivos a MINIMIZAR
objetivos_min <- c("indice_exceso_red", "probabilidad_siniestro_red")
valores_objetivos <- resultados_limpios %>% select(all_of(objetivos_min))

# Encontrar los puntos no dominados (el Frente de Pareto)
# Usamos un bucle, que es más fácil de leer. Para datasets muy grandes, hay métodos más rápidos.
es_pareto <- !sapply(1:nrow(valores_objetivos), function(i) {
  is_dominated(as.numeric(valores_objetivos[i, ]), valores_objetivos[-i, ])
})

frente_pareto <- resultados_limpios[es_pareto, ]

cat(paste("Se encontraron", nrow(frente_pareto), "escenarios en el Frente de Pareto.\n\n"))

## Se encontraron 13 escenarios en el Frente de Pareto.

Visualización del frente de Pareto

En esta sección, se identifica el frente de pareto y la solución de compromiso. Se ha decidido enfocar el análisis en minimizar dos objetivos clave: el índice de exceso y la probabilidad de siniestro, ya que representan las metas de seguridad más directas.

Para pasar de los resultados de la simulación a una recomendación estratégica, realizamos un análisis de “compromiso”. El objetivo es identificar un único escenario que represente el mejor balance posible entre los tres objetivos de seguridad que buscamos minimizar: índice de exceso, índice de riesgo y probabilidad de siniestro.

Además, para cuantificar la mejora, incluimos la “condición actual” o punto base, que representa el rendimiento de la red sin ninguna de las 100 cámaras nuevas.

El método seleccionado es el de programación por compromiso (compromise programming), que consiste en:

1. Normalizar los valores de los tres objetivos para que todos estén en una escala comparable (de 0 a 1).
2. Calcular la distancia euclidiana de cada solución en el Frente de Pareto al “punto ideal” utópico (0, 0, 0).
3. Seleccionar el escenario que tenga la mínima distancia a este punto ideal. Esta solución es, por definición, la más balanceada.

# --- 1. DEFINIR EL PUNTO BASE (CONDICIÓN ACTUAL) ---
punto_base <- tibble(
  seed = -1,
  velocidad_promedio_red = 40.5916,
  indice_exceso_red = 0.1832,
  indice_riesgo_red = 4.0238, # Se mantiene para consistencia de datos
  probabilidad_siniestro_red = 0.5691
)


# Normalizar los datos del FRENTE DE PARETO
punto_compromiso <- frente_pareto %>%
  mutate(
    exceso_norm = scales::rescale(indice_exceso_red),
    siniestro_norm = scales::rescale(probabilidad_siniestro_red)
  ) %>%
  mutate(distancia_euclidiana = sqrt(exceso_norm^2 + siniestro_norm^2)) %>%
  slice_min(order_by = distancia_euclidiana, n = 1)

cat(paste("La mejor solución balanceada es la del Seed:", punto_compromiso$seed, "\n\n"))

## La mejor solución balanceada es la del Seed: 4417

# --- 3. CREACIÓN DE LAS 4 FIGURAS ESTILO ACADÉMICO ---
# --- 1. PREPARAR UN ÚNICO DATAFRAME PARA GRAFICAR ---

# Añadimos una columna 'tipo' a cada conjunto de datos para identificarlos
datos_grafico <- bind_rows(
  resultados_limpios %>% mutate(tipo = "6000 escenarios simulados"),
  frente_pareto %>% mutate(tipo = "Frente de Pareto"),
  punto_base %>% mutate(tipo = "Situación actual"),
  punto_compromiso %>% mutate(tipo = "Solución de compromiso")
)

# Convertimos 'tipo' a un factor para controlar el orden de dibujado (y de la leyenda)
# Los primeros niveles se dibujan al fondo, los últimos encima.
datos_grafico$tipo <- factor(datos_grafico$tipo, levels = c(
  "6000 escenarios simulados", 
  "Frente de Pareto", 
  "Situación actual", 
  "Solución de compromiso"
))

# --- 2. DEFINIR LAS ESCALAS MANUALES (COLORES, FORMAS, TAMAÑOS) ---
# Definimos esto una vez para reutilizarlo en los tres gráficos.

# Paleta de colores nombrada
colores_manuales <- c(
  "6000 escenarios simulados"   = "grey80",
  "Frente de Pareto"       = "#0072B2",
  "Situación actual"       = "black",
  "Solución de compromiso" = "firebrick"
)

# Formas nombradas
formas_manuales <- c(
  "6000 escenarios simulados"   = 16,  # Círculo relleno
  "Frente de Pareto"       = 16,  # Círculo relleno
  "Situación actual"       = 4,   # Una 'X'
  "Solución de compromiso" = 17   # Triángulo
)

# Tamaños nombrados
tamanos_manuales <- c(
  "6000 escenarios simulados"   = 1.5,
  "Frente de Pareto"       = 3,
  "Situación actual"       = 5,
  "Solución de compromiso" = 5
)

# --- 3. CREAR LOS GRÁFICOS CON LEYENDA AUTOMÁTICA ---

# Proyección 1: Velocidad vs Exceso
ggplot(datos_grafico, aes(x = velocidad_promedio_red, y = indice_exceso_red, 
                           color = tipo, shape = tipo, size = tipo)) +
  geom_point(alpha = 0.7, stroke = 1.2) + # stroke es para el grosor de formas como la 'X'
  scale_color_manual(values = colores_manuales, name = "Leyenda") +
  scale_shape_manual(values = formas_manuales, name = "Leyenda") +
  scale_size_manual(values = tamanos_manuales, name = "Leyenda") +
  guides(color = guide_legend(override.aes = list(size = 4))) + # Hace los puntos de la leyenda más grandes
  labs(
    title = "Velocidad vs. índice de exceso",
    x = "Velocidad promedio (km/h)",
    y = "Índice de exceso"
  ) +
  theme_bw(base_size = 14)

# Proyección 2: Velocidad vs Siniestralidad
ggplot(datos_grafico, aes(x = velocidad_promedio_red, y = probabilidad_siniestro_red, 
                           color = tipo, shape = tipo, size = tipo)) +
  geom_point(alpha = 0.7, stroke = 1.2) +
  scale_color_manual(values = colores_manuales, name = "Leyenda") +
  scale_shape_manual(values = formas_manuales, name = "Leyenda") +
  scale_size_manual(values = tamanos_manuales, name = "Leyenda") +
  guides(color = guide_legend(override.aes = list(size = 4))) +
  labs(
    title = "Velocidad vs. probabilidad de siniestro",
    x = "Velocidad promedio (km/h)",
    y = "Probabilidad de siniestro"
  ) +
  theme_bw(base_size = 14)

# Proyección 3: Exceso vs Siniestralidad
ggplot(datos_grafico, aes(x = indice_exceso_red, y = probabilidad_siniestro_red, 
                           color = tipo, shape = tipo, size = tipo)) +
  geom_point(alpha = 0.7, stroke = 1.2) +
  scale_color_manual(values = colores_manuales, name = "Leyenda") +
  scale_shape_manual(values = formas_manuales, name = "Leyenda") +
  scale_size_manual(values = tamanos_manuales, name = "Leyenda") +
  guides(color = guide_legend(override.aes = list(size = 4))) +
  labs(
    title = "Frente de Pareto: Índice de exceso vs. probabilidad de siniestro",
    x = "Índice de exceso",
    y = "Probabilidad de siniestro"
  ) +
  theme_bw(base_size = 14)

Gráfico 3D del frente de pareto

# ==============================================================================
# --- Visualizaciones 3D del Espacio de Soluciones ---
# A continuación se presentan cuatro métodos para visualizar los mismos datos en 3D.
# Todos están configurados para funcionar con knitr y reutilizan la estética
# definida en la sección de gráficos 2D para mantener la consistencia.
# ==============================================================================

# --- Variables de ejes para consistencia ---
eje_x <- "Vel. promedio (km/h)"
eje_y <- "Índice de exceso"
eje_z <- "Prob. de siniestro"

# Separar los datos por tipo para aplicar estilos diferentes
escenarios <- datos_grafico %>% filter(tipo == "6000 escenarios simulados")
pareto_pts <- datos_grafico %>% filter(tipo == "Frente de Pareto")
base_pt    <- datos_grafico %>% filter(tipo == "Situación actual")
compromiso_pt <- datos_grafico %>% filter(tipo == "Solución de compromiso")


# ==============================================================================
# --- OPCIÓN 1: rgl (Interactivo, Rápido y Ligero) ---
# Ideal para exploración interactiva dentro del HTML. Usa WebGL.
# RECUERDA: Debes tener `rgl::setupKnitr()` en tu chunk de setup.
# ==============================================================================

datos_grafico_rgl <- bind_rows(
  resultados_limpios %>% mutate(tipo = "Escenarios simulados"),
  frente_pareto %>% mutate(tipo = "Frente de Pareto"),
  punto_base %>% mutate(tipo = "Situacion actual"), # SIN TILDE
  punto_compromiso %>% mutate(tipo = "Solucion de compromiso") # SIN TILDE
)

# Creamos el factor con los niveles ya limpios
datos_grafico_rgl$tipo <- factor(datos_grafico_rgl$tipo, levels = c(
  "Escenarios simulados", 
  "Frente de Pareto", 
  "Situacion actual", 
  "Solucion de compromiso"
))

# Ahora filtramos desde este dataframe 100% limpio
escenarios_rgl <- datos_grafico_rgl %>% filter(tipo == "Escenarios simulados")
pareto_rgl     <- datos_grafico_rgl %>% filter(tipo == "Frente de Pareto")
base_rgl       <- datos_grafico_rgl %>% filter(tipo == "Situacion actual")
compromiso_rgl <- datos_grafico_rgl %>% filter(tipo == "Solucion de compromiso")

# --- 2. DIBUJAR LA ESCENA ---
# Ahora que tanto los datos como las etiquetas están limpios, el gráfico funcionará.

bg3d("white")

aspect3d("iso")

# Capa 1: Nube de puntos
plot3d(
  x = escenarios_rgl$velocidad_promedio_red,
  y = escenarios_rgl$indice_exceso_red,
  z = escenarios_rgl$probabilidad_siniestro_red,
  type = "p",
  col = "grey80",
  size = 4,
  alpha = 0.3,
  main = "Opcion 1: Grafico 3D con rgl (Interactivo)",
  xlab = "Vel. promedio (kmh)", 
  ylab = "Indice de exceso", 
  zlab = "Prob. de siniestro"
)

# Capa 2: Frente de Pareto
points3d(
  x = pareto_rgl$velocidad_promedio_red, y = pareto_rgl$indice_exceso_red, z = pareto_rgl$probabilidad_siniestro_red,
  col = "#0072B2", size = 6
)

# Capa 3: Situación actual
points3d(
  x = base_rgl$velocidad_promedio_red, y = base_rgl$indice_exceso_red, z = base_rgl$probabilidad_siniestro_red,
  col = "black", pch = 4, size = 8, lwd = 2
)

# Capa 4: Solución de compromiso
points3d(
  x = compromiso_rgl$velocidad_promedio_red, y = compromiso_rgl$indice_exceso_red, z = compromiso_rgl$probabilidad_siniestro_red,
  col = "firebrick", pch = 17, size = 8
)

# --- 3. AÑADIR LEYENDA (También con texto limpio) ---
legend3d(
  "topright",
  legend = levels(datos_grafico_rgl$tipo),
  pch = c(16, 16, 4, 17),
  col = c("grey80", "#0072B2", "black", "firebrick"),
  cex = 1,
  inset = c(0.02)
)

# --- 4. ESTABLECER VISTA INICIAL ---
view3d(theta = 45, phi = 25, zoom = 0.8)

# ==============================================================================
# --- OPCIÓN 2: plotly (Interactivo, Moderno y con Tooltips) ---
# La mejor opción para reportes HTML modernos. Permite ver información al
# pasar el cursor sobre los puntos.
# ==============================================================================

plot_ly(data = datos_grafico) %>%
  # Añadimos las trazas en el orden correcto de dibujado (de fondo a frente)
  add_trace(
    data = escenarios,
    x = ~velocidad_promedio_red, y = ~indice_exceso_red, z = ~probabilidad_siniestro_red,
    type = 'scatter3d', mode = 'markers',
    marker = list(color = colores_manuales["6000 escenarios simulados"], opacity = 0.2, size = 3),
    hoverinfo = 'none', # Desactivamos tooltips para la nube de puntos
    name = "6000 escenarios simulados"
  ) %>%
  add_trace(
    data = pareto_pts,
    x = ~velocidad_promedio_red, y = ~indice_exceso_red, z = ~probabilidad_siniestro_red,
    type = 'scatter3d', mode = 'markers',
    marker = list(color = colores_manuales["Frente de Pareto"], size = 5),
    text = ~paste("Seed:", seed), hoverinfo = 'text',
    name = "Frente de Pareto"
  ) %>%
  add_trace(
    data = base_pt,
    x = ~velocidad_promedio_red, y = ~indice_exceso_red, z = ~probabilidad_siniestro_red,
    type = 'scatter3d', mode = 'markers',
    marker = list(color = colores_manuales["Situación actual"], size = 6, symbol = 'cross'),
    name = "Situación actual"
  ) %>%
  add_trace(
    data = compromiso_pt,
    x = ~velocidad_promedio_red, y = ~indice_exceso_red, z = ~probabilidad_siniestro_red,
    type = 'scatter3d', mode = 'markers',
    marker = list(color = colores_manuales["Solución de compromiso"], size = 7, symbol = 'diamond'),
    text = ~paste("Seed Óptimo:", seed), hoverinfo = 'text',
    name = "Solución de compromiso"
  ) %>%
  # Layout y estilo académico
  layout(
    title = "Opción 2: Gráfico 3D con Plotly (Interactivo)",
    scene = list(
      xaxis = list(title = eje_x, gridcolor = 'rgb(220, 220, 220)'),
      yaxis = list(title = eje_y, gridcolor = 'rgb(220, 220, 220)'),
      zaxis = list(title = eje_z, gridcolor = 'rgb(220, 220, 220)'),
      bgcolor = "white",
      camera = list(eye = list(x = 1.8, y = -1.8, z = 0.8)) # Buen ángulo inicial
    ),
    legend = list(orientation = 'h', xanchor = "center", x = 0.5, y = -0.1)
  )

# ==============================================================================
# --- OPCIÓN 3: plot3D (Estático, Estilo Base-R) ---
# Genera un gráfico 2D estático que simula una perspectiva 3D.
# Es menos flexible que rgl pero no requiere WebGL.
# ==============================================================================

# Se necesita llamar a scatter3D una vez para crear la escena, y luego a points3D
# para añadir las demás capas.
scatter3D(
  x = escenarios$velocidad_promedio_red,
  y = escenarios$indice_exceso_red,
  z = escenarios$probabilidad_siniestro_red,
  col = colores_manuales["6000 escenarios simulados"],
  alpha = 0.2, pch = 16, cex = 0.5,
  bty = "g", ticktype = "detailed",
  main = "Opción 3: Gráfico 3D con plot3D (Estático)",
  xlab = paste("\n\n\n", eje_x), 
  ylab = paste("\n\n\n", eje_y), 
  zlab = paste("\n\n\n", eje_z),
  theta = 40, phi = 25
)

# Añadir las otras capas
points3D(
  x = pareto_pts$velocidad_promedio_red, y = pareto_pts$indice_exceso_red, z = pareto_pts$probabilidad_siniestro_red,
  add = TRUE, col = colores_manuales["Frente de Pareto"], pch = 16, cex = 1
)
points3D(
  x = base_pt$velocidad_promedio_red, y = base_pt$indice_exceso_red, z = base_pt$probabilidad_siniestro_red,
  add = TRUE, col = colores_manuales["Situación actual"], pch = formas_manuales["Situación actual"], cex = 2, lwd = 2
)
points3D(
  x = compromiso_pt$velocidad_promedio_red, y = compromiso_pt$indice_exceso_red, z = compromiso_pt$probabilidad_siniestro_red,
  add = TRUE, col = colores_manuales["Solución de compromiso"], pch = formas_manuales["Solución de compromiso"], cex = 2
)

# Añadir leyenda manualmente
legend(
  "topright", 
  legend = names(colores_manuales),
  pch = unname(formas_manuales),
  col = unname(colores_manuales),
  cex = 0.8,
  ncol = 2,      # <-- La clave para que no se corte
  bty = "n"      # <-- Quita la caja de la leyenda
)

# ==============================================================================
# --- OPCIÓN 4: rayshader (Estático, Alta Calidad para Publicación) ---
# Crea una renderización de alta calidad a partir de un ggplot 2D.
# El resultado es una imagen estática, ideal para PDFs o artículos.
# NOTA: Puede ser lento de compilar en el knitr.
# ==============================================================================

datos_grafico_rayshader <- datos_grafico %>%
  mutate(
    # Asignamos un tamaño numérico basado en el tipo
    size_numeric = case_when(
      tipo == "6000 escenarios simulados" ~ 1.5,
      tipo == "Frente de Pareto"       ~ 3.0,
      tipo == "Situación actual"       ~ 5.0,
      tipo == "Solución de compromiso" ~ 5.0
    )
  )

p_base_rayshader <- ggplot(datos_grafico_rayshader, 
                           # El eje Z (velocidad) se mapea a un color degradado
                           # que rayshader usa para la altura.
                           aes(x = indice_exceso_red, 
                               y = probabilidad_siniestro_red,
                               color = velocidad_promedio_red)) +
  # Dibujamos las capas de fondo a frente
  geom_point(data = . %>% filter(tipo == "6000 escenarios simulados"), alpha = 0.3, size = 1.5) +
  geom_point(data = . %>% filter(tipo == "Frente de Pareto"), size = 3) +
  geom_point(data = . %>% filter(tipo == "Situación actual"), color = "black", shape = 4, size = 5, stroke = 1.5) +
  geom_point(data = . %>% filter(tipo == "Solución de compromiso"), color = "firebrick", shape = 17, size = 5) +
  scale_color_viridis_c(name = "Velocidad Promedio") +
  labs(
    x = eje_y,
    y = eje_z,
    title = "Opción 4: Gráfico 3D con rayshader (Estático)"
  ) +
  theme_classic(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "bottom")


# 2. Renderizar la escena 3D con plot_gg
# NOTA: `plot_gg` puede ser lento, pero `setupKnitr` lo capturará.
plot_gg(
  p_base_rayshader,
  multicore = TRUE,
  width = 7, height = 7,
  scale = 300,
  zoom = 0.60, theta = 45, phi = 30,
  solid = TRUE, solidcolor = "white",
  shadow = TRUE
)

render_camera(theta = 45, phi = 30, zoom = 0.65)

Identificación de ubicaciones óptimas de cámaras

El resultado más accionable de este análisis es saber dónde se ubicaron las cámaras en los escenarios más prometedores.

Para lograr esto, se implementa una función en R que replica la lógica de selección aleatoria de Python. Esta función toma un seed como entrada y, usando ese mismo seed, reproduce la selección de los 100 tramos pk_id_calz que se hizo en la simulación original. Al ejecutar esta función para los seeds de los escenarios que identificamos como óptimos en el Frente de Pareto, podemos obtener la lista exacta de ubicaciones para cada estrategia ganadora.

Mapa 1: Ubicaciones para la Solución de Compromiso

A continuación, se muestra la distribución geográfica de las 100 nuevas cámaras propuestas por la solución de compromiso, el escenario que ofrece el mejor balance entre la reducción del índice de exceso y la probabilidad de siniestro.

# --- IDENTIFICAR EL SEED ÓPTIMO ---
seed_optimo <- punto_compromiso$seed

# --- RECREAR LAS UBICACIONES PARA ESE SEED ---
ubicaciones_optimas_final <- obtener_ubicaciones_por_seed(
  seed_id = seed_optimo, 
  datos_base_sf = vias_sf, 
  n_camaras_nuevas = 100
)

# --- PREPARAR DATOS Y GENERAR EL MAPA ---
vias_mapa_compromiso <- vias_sf %>%
  mutate(
    tipo_camara = case_when(
      pk_id_calz %in% ubicaciones_optimas_final ~ "Nueva (solución óptima)",
      tiene_cama == 1 ~ "Existente",
      TRUE ~ "Sin cámara"
    ),
    tipo_camara = factor(tipo_camara, levels = c("Nueva (solución óptima)", "Existente", "Sin cámara"))
  )

ggplot() +
  geom_sf(data = vias_mapa_compromiso, aes(color = tipo_camara, linewidth = tipo_camara), alpha = 0.9) +
  scale_color_manual(
    name = "Tipo de cámara",
    values = c("Nueva (solución óptima)" = "#0072B2", "Existente" = "firebrick", "Sin cámara" = "grey80")
  ) +
  scale_linewidth_manual(
    values = c("Nueva (solución óptima)" = 0.8, "Existente" = 0.8, "Sin cámara" = 0.15),
    guide = "none" # Ocultar la leyenda del grosor
  ) +
  labs(
    title = "Ubicación de cámaras para la solución de compromiso",
    subtitle = paste("Escenario:", seed_optimo)
  ) +
  theme_void() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    legend.position = "bottom"
  )

Mapa 2: Mapa de Calor de Ubicaciones Prioritarias

Para obtener una visión más robusta, se analizaron los escenarios más destacados (los 10 mejores para cada métrica y todos los puntos del Frente de Pareto). El siguiente mapa de calor muestra la frecuencia con la que cada tramo de vía fue seleccionado en este conjunto de escenarios óptimos. Las áreas más intensas representan las ubicaciones más consistentemente prioritarias.

# Selección de los mejores 10 escenarios por cada variable
top_seeds <- list(
  arrange(resultados_limpios, desc(velocidad_promedio_red)),
  arrange(resultados_limpios, indice_exceso_red),
  arrange(resultados_limpios, indice_riesgo_red),
  arrange(resultados_limpios, probabilidad_siniestro_red)
) %>%
  map(~ head(., 10) %>% pull(seed)) %>%
  unlist()

seeds_optimos <- unique(c(top_seeds, frente_pareto$seed))

# 2. Recreación de ubicación por tramos y conteo de frecuencias
# Cargar los datos geoespaciales base
vias_sf <- st_read(ruta_shapefile_vias, quiet = TRUE) %>%
  janitor::clean_names()

# Replicar la función de generación de escenarios en R
obtener_ubicaciones_por_seed <- function(seed_id, datos_base_sf, n_camaras_nuevas) {
  
  # Asegurar reproducibilidad
  set.seed(seed_id)
  
  # Identificar tramos elegibles (sin cámaras base)
  elegibles <- datos_base_sf %>% 
    filter(tiene_cama == 0) %>%
    pull(pk_id_calz)
  
  # Seleccionar n_camaras_nuevas al azar
  n_a_tomar <- min(n_camaras_nuevas, length(elegibles))
  
  ubicaciones_elegidas <- sample(elegibles, size = n_a_tomar, replace = FALSE)
  
  return(ubicaciones_elegidas)
}


# Obtener una lista de TODOS los tramos seleccionados en los escenarios óptimos
lista_completa_tramos <- unlist(lapply(seeds_optimos, function(s) {
  obtener_ubicaciones_por_seed(seed_id = s, datos_base_sf = vias_sf, n_camaras_nuevas = 100)
}))

# Contar la frecuencia de aparición de cada tramo
frecuencia_tramos <- as.data.frame(table(lista_completa_tramos)) %>%
  rename(pk_id_calz = lista_completa_tramos, frecuencia = Freq) %>%
  mutate(pk_id_calz = as.numeric(as.character(pk_id_calz)))

# 3. Unión geoespacial y categorización
# Unir la tabla de frecuencias al shapefile
vias_mapeo <- vias_sf %>%
  left_join(frecuencia_tramos, by = "pk_id_calz") %>%
  # Crear una columna categórica para el mapeo
  mutate(
    tipo_tramo = case_when(
      tiene_cama == 1 ~ "Cámara Existente",
      frecuencia > 0 ~ "Ubicación Óptima Sugerida",
      TRUE ~ "Red Vial Base"
    ),
    # Asegurar que la frecuencia sea numérica y que los NA sean 0
    frecuencia = replace_na(frecuencia, 0)
  )

# 4. Creación del mapa
# Extraer los centroides de los tramos con cámaras existentes para mapearlos como puntos
centroides_existentes <- vias_mapeo %>%
  filter(tipo_tramo == "Cámara Existente") %>%
  st_centroid()

## Warning: st_centroid assumes attributes are constant over geometries

ggplot() +
  # Capa 1: Toda la red vial en gris claro como base
  geom_sf(data = vias_mapeo, color = "grey80", size = 0.2, alpha = 0.5) +
  
  # Capa 2: Ubicaciones óptimas sugeridas, con estética mejorada
  geom_sf(data = vias_mapeo %>% filter(tipo_tramo == "Ubicación Óptima Sugerida"), 
          aes(color = frecuencia, 
              alpha = frecuencia,  # La transparencia también depende de la frecuencia
              size = frecuencia)) + # El grosor también depende de la frecuencia
  
  # Capa 3: Cámaras existentes como puntos rojos semi-transparentes
  geom_sf(data = centroides_existentes, color = "firebrick", size = 1, alpha = 0.7, shape = 16) +
  
  # Escala de colores para la frecuencia
  scale_color_viridis_c(
    name = "Frecuencia de instalación en escenarios óptimos",
    guide = guide_colorbar(
      direction = "horizontal", barheight = unit(2, units = "mm"),
      barwidth = unit(50, units = "mm"), title.position = 'top',
      title.hjust = 0.5, label.hjust = 0.5
    )
  ) +
  
  # Control manual de la transparencia y el grosor
  scale_alpha_continuous(range = c(0.4, 1.0), guide = "none") + # Desde apenas visible a totalmente opaco
  scale_size_continuous(range = c(0.4, 2), guide = "none") +  # Desde un poco más grueso a notablemente grueso
  
  labs(title = "Mapa de calor de ubicaciones óptimas para 100 nuevas cámaras",
       subtitle = "Basado en la frecuencia de aparición en los 40 mejores escenarios de simulación",
       caption = "Puntos rojos representan las cámaras existentes.") +
  
  theme_void() +
  theme(
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    plot.caption = element_text(size = 10, hjust = 0.5),
    legend.position = "bottom"
  )

Tabla de tramos y corredores prioritarios

Para obtener una visión más robusta, se analizaron los escenarios más destacados (los 10 mejores para cada métrica y todos los puntos del Frente de Pareto). El siguiente mapa de calor muestra la frecuencia con la que cada tramo de vía fue seleccionado en este conjunto de escenarios óptimos. Las áreas más intensas representan las ubicaciones más consistentemente prioritarias.

# Tablas con tramos y corredores más frecuentes

# Tabla 1: TOP 20 TRAMOS más frecuentes
# Seleccionamos los tramos con frecuencia > 0 y los unimos con la información del corredor
top_tramos <- vias_mapeo %>%
  filter(frecuencia > 0) %>%
  st_drop_geometry() %>% # Quitamos la geometría para trabajar con la tabla
  select(pk_id_calz, tipos_via, frecuencia) %>%
  arrange(desc(frecuencia)) %>%
  head(20)

knitr::kable(top_tramos,
             col.names = c("ID del tramo (pk_id_calz)", "Corredor vial (tipos_via)", "Frecuencia"),
             caption = "Top 20 tramos prioritarios para la ubicación de 100 nuevas cámaras",
             align = "lcr") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Top 20 tramos prioritarios para la ubicación de 100 nuevas cámaras

ID del tramo (pk_id_calz)	Corredor vial (tipos_via)	Frecuencia
503605	AVENIDA JORGE URIBE BOTERO ( KR 35 )	4
603810	AVENIDA LUIS CARLOS GALAN SARMIENTO ( AC 39 )	4
385342	AVENIDA FERROCARRIL DE OCCIDENTE ( CL 28 )	4
298272	AVENIDA CARACAS ( )	4
91025363	TRANSVERSAL DE SUBA ( KR 61 )	4
184483	AVENIDA FUCHA ( CL 8 SUR Y CL 11 SUR )	3
24122993	AVENIDA SANTA BARBARA ( KR 19 )	3
24122773	AVENIDA BOYACA ( )	3
91024696	AVENIDA BOYACA ( )	3
24119890	AVENIDA CIUDAD DE QUITO ( AK 30 )	3
24122154	AVENIDA ALBERTO LLERAS CAMARGO ( AK 7 )	3
24120850	AVENIDA CIUDAD DE QUITO ( AK 30 )	3
91018066	AVENIDA CIUDAD DE CALI ( )	3
24122940	AVENIDA DE LAS VILLAS ( KR 52- KR 56 Y CL 125A Y CL 190 )	3
91019449	AVENIDA LAUREANO GOMEZ ( AK 9 )	3
504137	AVENIDA CARACAS ( )	3
213434	AVENIDA DE LOS CERROS ( AV CIRCUNVALAR )	3
514327	AVENIDA CHILE ( CL 72 )	3
24119835	AVENIDA CALLE 92 ( CL 92 )	3
200456	AVENIDA FUCHA ( CL 8 SUR Y CL 11 SUR )	3

# Tabla 2: TOP 15 CORREDORES con mayor frecuencia acumulada
# Agrupamos por corredor para calcular las métricas totales y de selección
top_corredores <- vias_mapeo %>%
  st_drop_geometry() %>%
  # Primero, obtenemos los totales de cada corredor
  group_by(tipos_via) %>%
  summarise(
    longitud_total_km = first(longitud_k),
    numero_tramos_totales = n()
  ) %>%
  # Ahora, unimos la información de las selecciones
  left_join(
    vias_mapeo %>%
      st_drop_geometry() %>%
      filter(frecuencia > 0) %>%
      group_by(tipos_via) %>%
      summarise(
        frecuencia_absoluta = sum(frecuencia, na.rm = TRUE),
        numero_tramos_sugeridos = n()
      ),
    by = "tipos_via"
  ) %>%
  # Rellenar con 0 para corredores que nunca fueron seleccionados
  mutate(
    frecuencia_absoluta = replace_na(frecuencia_absoluta, 0),
    numero_tramos_sugeridos = replace_na(numero_tramos_sugeridos, 0)
  ) %>%
  # Calcular las métricas numéricas finales
  mutate(
    freq_rel_tramos = frecuencia_absoluta / numero_tramos_totales,
    pct_tramos_sugeridos = numero_tramos_sugeridos / numero_tramos_totales,
    freq_por_km = frecuencia_absoluta / longitud_total_km
  ) %>%
  # Ordenar por la métrica de interés y tomar el Top 15
  arrange(desc(frecuencia_absoluta)) %>%
  head(15) %>%
    mutate(
    # La función scales::percent convierte la proporción (ej. 0.25) a texto ("25.0%")
    `% Tramos sugeridos` = scales::percent(pct_tramos_sugeridos, accuracy = 0.1)
  ) %>%
  select(
    Corredor = tipos_via,
    `Frecuencia absoluta` = frecuencia_absoluta,
    `Frec. / Tramos totales` = freq_rel_tramos,
    `% Tramos sugeridos`, # Ya está formateado y nombrado
    `Frec. / km` = freq_por_km
  )

# --- Impresión de la Tabla con Formato (ahora más simple) ---
knitr::kable(
  top_corredores,
  caption = "Top 15 corredores prioritarios",
  digits = 2,
  align = "lcrrr" # Alineación de columnas: Izquierda, Centro, Derecha, Derecha, Derecha
) %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"), 
    full_width = FALSE, 
    position = "center"
  )

Top 15 corredores prioritarios

Corredor	Frecuencia absoluta	Frec. / Tramos totales	% Tramos sugeridos	Frec. / km
AVENIDA BOYACA ( )	315	0.35	29.5%	3.17
AVENIDA CARACAS ( )	236	0.40	32.1%	3.01
AVENIDA PASEO DE LOS LIBERTADORES ( AUTONORTE )	217	0.36	30.7%	2.56
AVENIDA CIUDAD DE QUITO ( AK 30 )	188	0.42	33.6%	2.68
AVENIDA ALBERTO LLERAS CAMARGO ( AK 7 )	174	0.39	33.1%	3.56
AVENIDA JORGE ELIECER GAITAN ( CL 26 )	169	0.40	33.3%	2.24
AVENIDA CIUDAD DE CALI ( )	145	0.37	30.3%	2.75
AVENIDA CIUDAD DE VILLAVICENCIO ( )	137	0.42	34.2%	4.19
AVENIDA DE LOS CERROS ( AV CIRCUNVALAR )	136	0.42	33.5%	3.51
AVENIDA PRIMERO DE MAYO ( )	122	0.37	31.4%	4.11
AVENIDA MEDELLIN ( CL 80 )	113	0.46	39.8%	2.63
AVENIDA DEL CONGRESO EUCARISTICO ( AK 68 )	101	0.39	32.6%	2.04
AVENIDA CHILE ( CL 72 )	96	0.42	35.7%	4.98
AVENIDA DE LAS AMERICAS ( )	94	0.39	32.6%	2.93
AVENIDA BATALLON CALDAS ( KR 50 Y TV 42 )	89	0.44	34.7%	4.22

Modelo UNGSA-III