# Standard library
import warnings

# Third-party libraries
import holidays
import joblib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import xgboost as xgb

from sklearn.compose import ColumnTransformer
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# Local imports
warnings.filterwarnings('ignore')
import viz_utils as viz # funciones de visualización personalizadas para el proyecto
viz.set_naturgy_style()

# Comprobración optativa: para gráficas como RTX Laptop 2000 que no son compatibles de forma nativa con xgboost.
import os
print(os.environ.get('LD_LIBRARY_PATH', 'NO DEFINIDA'))

viz.show_palette()

df_energy = pd.read_parquet('../data/energy_silver.parquet/')
viz.display_dataframe_summary(df_energy, name='Dataset de Energí­a')
# Se han transpuesto las filas por columnas para una visualización más legible.

df_weather = pd.read_parquet('../data/weather_features_silver.parquet')
viz.display_dataframe_summary(df_weather, name='Dataset de Clima')

viz.show_unique_values(df_weather, name='Dataset de Clima')

cols_to_drop_energy = [
    'generation hydro pumped storage aggregated',
    'forecast wind offshore eday ahead',
    'generation fossil coal-derived gas',
    'generation fossil oil shale',
    'generation fossil peat',
    'generation geothermal',
    'generation marine',
    'generation wind offshore',
    'price day ahead',
    
    # Producen data leakage
    'generation biomass',
    'generation fossil brown coal/lignite',
    'generation fossil gas',
    'generation fossil hard coal',
    'generation fossil oil',
    'generation hydro run-of-river and poundage',
    'generation hydro water reservoir',
    'generation nuclear',
    'generation other',
    'generation other renewable',
    'generation solar',
    'generation waste',
    'generation wind onshore'
]

df_energy_clean = df_energy.drop(columns=cols_to_drop_energy, errors='ignore')

print(f"Dataset original: {df_energy.shape[1]} columnas")
print(f"Dataset limpio:   {df_energy_clean.shape[1]} columnas")

print(f"Columnas en clima antes: {df_weather.shape[1]}")
df_weather = df_weather.drop(columns=["weather_description"])
print(f"Columnas en clima después: {df_weather.shape[1]}")
print("Columna 'weather_description' eliminada correctamente.")

numeric_cols_weather = df_weather.select_dtypes(include=np.number).columns

# Conversión a componente vectorial del viento (componente este-oeste)
df_weather['wind_x'] = df_weather['wind_speed'] * np.cos(np.radians((90 - df_weather['wind_deg']) % 360))

df_weather['city_name'] = df_weather['city_name'].str.strip()

weather_mapping = {
    'clear': 'Clear',
    'clouds': 'Cloudy',
    'rain': 'Rain',
    'drizzle': 'Rain',
    'thunderstorm': 'Storm',
    'snow': 'Snow',
    'squall': 'Wind',  # Ráfagas de viento como categoría separada
    'mist': 'Visibility_Issues',
    'fog': 'Visibility_Issues',
    'haze': 'Visibility_Issues',
    'smoke': 'Visibility_Issues',
    'dust': 'Visibility_Issues'
}

df_weather["weather_grouped"] = df_weather["weather_main"].str.lower().replace(weather_mapping)

#  One-Hot Encoding
clima_dummies = pd.get_dummies(df_weather["weather_grouped"], prefix="weather")

# Unión al dataset original
df_weather = pd.concat([df_weather, clima_dummies], axis=1)

pesos_ciudades = {
    'Madrid': 0.28,
    'Barcelona': 0.28,
    'Seville': 0.18,
    'Valencia': 0.16,
    'Bilbao': 0.10
}
df_weather["peso"] = df_weather["city_name"].map(pesos_ciudades)

columnas_numericas = [
    'temp', 'humidity', 'pressure', 'wind_speed', 'wind_x',
    'clouds_all', 'rain_1h',
    'weather_Clear', 'weather_Cloudy', 'weather_Rain', 'weather_Storm', 
    'weather_Visibility_Issues', 'weather_Snow', 'weather_Wind'
]


columnas_existentes = [col for col in columnas_numericas if col in df_weather.columns]

def promedio_ponderado(df_group):
    weights = df_group["peso"].values
    result = {}
    for col in columnas_existentes: 
        result[col] = np.average(df_group[col], weights=weights)
    return pd.Series(result)

# Agrupación final: 1 fila por timestamp
df_weather_nacional = df_weather.groupby("time").apply(promedio_ponderado).reset_index()

viz.display_describe(df_weather_nacional, name="Clima Nacional Ponderado")

print("Resumen")
print(f"Dataset de energía:  {df_energy_clean.shape}")
print(f"Dataset de clima:    {df_weather_nacional.shape}")

print(f"Rango temporal energía: {df_energy_clean['time'].min()} -> {df_energy_clean['time'].max()}")
print(f"Rango temporal clima:    {df_weather_nacional['time'].min()} -> {df_weather_nacional['time'].max()}")

print(f"Columnas energía: {list(df_energy_clean.columns)}")
print(f"Columnas clima: {list(df_weather_nacional.columns)}")

df_energy['time'] = pd.to_datetime(df_energy['time'])

# Filas donde el salto con la anterior sea mayor de 1 hora
saltos_temporales = df_energy[df_energy['time'].diff() > pd.Timedelta(hours=1)]

print(saltos_temporales['time'])

df_master = pd.merge(df_weather_nacional, df_energy_clean, on='time', how='outer')

df_master = df_master.sort_values('time').reset_index(drop=True)

print("Nulos antes de interpolar:")
print(df_master.isnull().sum()[df_master.isnull().sum() > 0])

# Interpolación lineal
columnas_a_interpolar = df_master.select_dtypes(include=[np.number]).columns
df_master[columnas_a_interpolar] = df_master[columnas_a_interpolar].interpolate(method='linear')

print("\nTabla final")
print(f"Filas totales: {df_master.shape[0]}")
print(f"Nulos restantes en el dataset: {df_master.isnull().sum().sum()}")

viz.display_dataframe_summary(df_master, name='Dataset Maestro')

viz.plot_target_correlation(df_master, target_col='price actual', figsize=(6, 8))

viz.plot_target_scatter(
    df_master,
    target_col='price actual',
    n_cols=3,
    figsize=(18, None),
    alpha=0.1
)

viz.plot_target_boxplot_by_time(
    df_master,
    time_col='time',
    target_col='price actual',
    figsize=(10, 5),
    title="Distribución del precio de la Luz por hora del día",
    show_stats=True
)

cols_legales = [
    'temp', 'humidity', 'pressure', 'wind_speed', 'clouds_all',
    'total load forecast', 'forecast solar day ahead', 'forecast wind onshore day ahead',
    'price actual'
]

viz.plot_correlation_heatmap(
    df_master,
    columns=cols_legales,
    figsize=(12, 10),
    title="Correlación entre predictores legales y precio"
)

cols_clima = ['temp', 'humidity', 'clouds_all', 'weather_Clear', 'weather_Cloudy', 'weather_Rain']

viz.plot_correlation_heatmap(
    df_master,
    columns=cols_clima,
    figsize=(8, 6),
    title="Multicolinealidad en Variables Climáticas"
)

columnas_a_borrar = [
    'weather_Wind', 
    # 'wind_x',
    'weather_Snow', 
    'weather_Storm', 
    'weather_Visibility_Issues', 
    'weather_Clear',
    'weather_Cloudy'
]

df_master = df_master.drop(columns=columnas_a_borrar)

# Protección contra el KeyError: si 'time' sigue siendo columna, la pasamos a índice.
if 'time' in df_master.columns:
    df_master = df_master.set_index('time')

# Índice como fecha sin zona horaria (para que cuadre con los festivos)
df_master.index = pd.to_datetime(df_master.index, utc=True).tz_localize(None)

es_holidays = holidays.Spain(years=range(2015, 2020))
festivos_dt = pd.to_datetime(sorted(es_holidays.keys()))

# festivo?
df_master['is_holiday'] = df_master.index.normalize().isin(festivos_dt).astype(int)

# Es víspera de festivo?
visperas = festivos_dt - pd.Timedelta(days=1)
df_master['is_holiday_eve'] = df_master.index.normalize().isin(visperas).astype(int)

# Feature numérica: días hasta/desde el festivo
fechas_actuales = df_master.index.normalize()

def calc_days_to(fecha):
    diffs = (festivos_dt - fecha).days
    futuros = diffs[diffs >= 0]
    return futuros.min() if len(futuros) > 0 else 0

def calc_days_from(fecha):
    diffs = (fecha - festivos_dt).days
    pasados = diffs[diffs >= 0]
    return pasados.min() if len(pasados) > 0 else 0

df_master['days_to_holiday'] = fechas_actuales.map(calc_days_to)
df_master['days_from_holiday'] = fechas_actuales.map(calc_days_from)

print(f"Festivos cargados: {len(festivos_dt)}")
print(f"Registros en festivo: {df_master['is_holiday'].sum()}")
print(f"Registros en víspera: {df_master['is_holiday_eve'].sum()}")

# Extracción de componentes temporales
df_master['hour'] = df_master.index.hour
df_master['dayofweek'] = df_master.index.dayofweek
df_master['month'] = df_master.index.month
df_master['dayofyear'] = df_master.index.dayofyear

# Transformaciones cíclicas
ciclos = {
    'hour': 24,
    'dayofweek': 7,
    'month': 12
}

for col, periodo in ciclos.items():
    df_master[f'{col}_sin'] = np.sin(2 * np.pi * df_master[col] / periodo)
    df_master[f'{col}_cos'] = np.cos(2 * np.pi * df_master[col] / periodo)

# Borrado de las columnas temporales
df_master = df_master.drop(columns=['hour', 'dayofweek', 'month', 'dayofyear'])

# TENDENCIA: Media móvil de 24 horas y 168 horas (desplazadas 1h)
df_master['price_roll_mean_24h'] = df_master['price actual'].shift(1).rolling(window=24).mean()
df_master['price_roll_mean_168h'] = df_master['price actual'].shift(1).rolling(window=168).mean()

# VOLATILIDAD: Desviación estándar de las últimas 24 horas (desplazada 1h)
df_master['price_roll_std_24h'] = df_master['price actual'].shift(1).rolling(window=24).std()

# EXTREMOS: Máximos y mínimos recientes (desplazados 1h)
df_master['price_roll_max_24h'] = df_master['price actual'].shift(1).rolling(window=24).max()
df_master['price_roll_min_24h'] = df_master['price actual'].shift(1).rolling(window=24).min()

df_master = df_master.dropna()

lags_clave = [1, 24, 168]

def add_lags(df, col, lags):
    for lag in lags:
        df[f'{col}_lag_{lag}h'] = df[col].shift(lag)
    return df

df_master = add_lags(df_master, 'price actual', lags_clave)
df_master = add_lags(df_master, 'total load actual', lags_clave)

df_master = df_master.dropna()

print(f"Dimensiones: {df_master.shape}")

esta_ordenado = df_master.index.is_monotonic_increasing

diferencia_tiempo = df_master.index.to_series().diff().dropna().value_counts()

print(f"¿Está el índice ordenado cronológicamente?: {esta_ordenado}")
print("Frecuencia de los saltos de tiempo:")
print(diferencia_tiempo)

viz.display_dataframe_summary(df_master, name='Dataset Maestro (tras Feature Engineering)')

models_path = './bestmodels/'

# Variable objetivo: el precio de dentro de 24 horas.
df_master['target_price_24h'] = df_master['price actual'].shift(-24)

df_master = df_master.dropna()

# Fecha de corte
cutoff_date = pd.to_datetime('2018-01-01')

target_col = 'target_price_24h'
cols_no_features = [target_col, 'price actual', 'total load actual']

# Sets de enrenamiento y test
df_train_val = df_master[df_master.index < cutoff_date]
df_final_test = df_master[df_master.index >= cutoff_date]

# Conjunto de evaluacion final
X_final_test = df_final_test.drop(columns=cols_no_features, errors='ignore')
y_final_test = df_final_test[target_col]

# Conjunto de entrenamiento
X_dev_features = df_train_val.drop(columns=cols_no_features, errors='ignore')
y_dev_targets = df_train_val[target_col]

X_train, X_val, y_train, y_val = train_test_split(
    X_dev_features, y_dev_targets, test_size=0.2, random_state=42, shuffle=False
)

numeric_features = X_train.select_dtypes(include=np.number).columns.tolist()

preprocessor = ColumnTransformer(
    transformers=[
        ('num', StandardScaler(), numeric_features)
    ]
)

preprocessor.fit(X_train)

xgb_pipeline = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('regressor', xgb.XGBRegressor(
        objective='reg:squarederror',
        random_state=42,
        n_estimators=6000,
        early_stopping_rounds=100,
        max_depth=4,
        min_child_weight=8,
        learning_rate=0.003,
        reg_alpha=2.0,
        reg_lambda=2.0,
        subsample=0.8,
        colsample_bytree=0.8,
        tree_method='hist',
        device='cuda'
    ))
])

# Transformamos train y val para el seguimiento del error
X_train_t = preprocessor.transform(X_train)
X_val_t = preprocessor.transform(X_val)

xgb_pipeline.fit(
    X_train, y_train,
    regressor__eval_set=[(X_train_t, y_train), (X_val_t, y_val)],
    regressor__verbose=500
)

print(f"Entrenamiento finalizado. Árbol óptimo: {xgb_pipeline.named_steps['regressor'].best_iteration}")

model = xgb_pipeline.named_steps['regressor']
results = model.evals_result()

plt.figure(figsize=(10, 5))
plt.plot(results['validation_0']['rmse'], label='Entrenamiento (Train)', color=viz.NOBLUE[0])
plt.plot(results['validation_1']['rmse'], label='Validación (Val)', color=viz.NORANGE[0])
plt.title('Evolución del error (RMSE) durante el entrenamiento')
plt.xlabel('Iteraciones (n_estimators)')
plt.ylabel('RMSE')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

model_filename_xgb = models_path + '/xgb_model-energy_weather.joblib'
joblib.dump(xgb_pipeline, model_filename_xgb)
print(f"Modelo XGBoost guardado en: {model_filename_xgb}")

xgb_pipeline = joblib.load(model_filename_xgb)

y_pred = xgb_pipeline.predict(X_final_test)

# Ventana de una semana
plt.figure(figsize=(15, 6))
plt.plot(y_final_test.index[:168], y_final_test[:168], label='Precio Real', color=viz.NOBLUE[0], alpha=0.8)
plt.plot(y_final_test.index[:168], y_pred[:168], label='Predicción XGBoost', color=viz.NORANGE[0], linestyle='--')

plt.title('Comparativa: Precio Real vs Predicción (Primera semana de 2018)')
plt.xlabel('Fecha')
plt.ylabel('Precio (€/MWh)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Prodicciones sobre el conjunto de validación
y_val_pred = xgb_pipeline.predict(X_val)

rmse_val = np.sqrt(mean_squared_error(y_val, y_val_pred))
mae_val = mean_absolute_error(y_val, y_val_pred)
r2_val = r2_score(y_val, y_val_pred)

print("--- Métricas de Rendimiento en Validación ---")
print(f"RMSE (Error cuadrático medio): {rmse_val:.2f} €/MWh")
print(f"MAE (Error absoluto medio):    {mae_val:.2f} €/MWh")
print(f"R² (Coeficiente de determinación): {r2_val:.4f}")

residuos = y_val - y_val_pred

plt.figure(figsize=(10, 5))
plt.scatter(y_val, residuos, alpha=0.3, color=viz.NOBLUE[1])
plt.axhline(0, color=viz.NORANGE[0], linestyle='--')
plt.title('Distribución de residuos (¿El modelo falla más cuando el precio es alto?)')
plt.xlabel('Precio Real')
plt.ylabel('Error (Real - Pred)')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

#Evaluación final sobre el conjunto de test aislado (2018)
# (xgb_pipeline ya se cargó en la sección de evaluación; el pipeline escala internamente)
y_final_pred = xgb_pipeline.predict(X_final_test)

rmse_final = np.sqrt(mean_squared_error(y_final_test, y_final_pred))
mae_final = mean_absolute_error(y_final_test, y_final_pred)
r2_final = r2_score(y_final_test, y_final_pred)

print(f"RMSE (Error cuadrático medio): {rmse_final:.2f} €/MWh")
print(f"MAE  (Error absoluto medio):   {mae_final:.2f} €/MWh")
print(f"R²   (Coeficiente de determinación): {r2_final:.4f}")
print(f"\nNº de registros en test: {len(y_final_test)}")

plt.figure(figsize=(16, 6))
plt.plot(y_final_test.index[:720], y_final_test[:720],
         label='Precio Real (€/MWh)', color=viz.NOBLUE[0], linewidth=1.5)
plt.plot(y_final_test.index[:720], y_final_pred[:720],
         label='Predicción XGBoost', color=viz.NORANGE[0], linestyle='--', linewidth=1.5)
plt.title('Evaluación Final: Real vs Predicción — Enero 2018 (Test Aislado)', fontsize=14)
plt.xlabel('Fecha')
plt.ylabel('Precio (€/MWh)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

residuos_final = y_final_test - y_final_pred

plt.figure(figsize=(14, 5))

plt.subplot(1, 2, 1)
plt.scatter(y_final_test, residuos_final, alpha=0.4, color=viz.NOBLUE[1])
plt.axhline(0, color=viz.NORANGE[0], linestyle='--')
plt.title('Residuos vs Precio Real (Test 2018)')
plt.xlabel('Precio Real (€/MWh)')
plt.ylabel('Residuo (Real - Pred)')

plt.subplot(1, 2, 2)
plt.hist(residuos_final, bins=60, color=viz.NOBLUE[1], edgecolor='white', alpha=0.8)
plt.title('Distribución de Residuos (Test 2018)')
plt.xlabel('Error (€/MWh)')
plt.ylabel('Frecuencia')

plt.tight_layout()
plt.show()

print(f"Error medio (bias): {residuos_final.mean():.3f} €/MWh")
print(f"Error máximo absoluto: {residuos_final.abs().max():.2f} €/MWh")
print(f"Percentil 95 del |error|: {residuos_final.abs().quantile(0.95):.2f} €/MWh")

# Importancia de variables del modelo XGBoost
model = xgb_pipeline.named_steps['regressor']
feat_names = [n.replace('num__', '') for n in xgb_pipeline.named_steps['preprocessor'].get_feature_names_out()]
importancias = pd.Series(model.feature_importances_, index=feat_names).sort_values(ascending=False)

top_n = 20
top = importancias.head(top_n).iloc[::-1]

fig, ax = plt.subplots(figsize=(10, 8))
ax.barh(top.index, top.values, color=viz.NOBLUE[0])
ax.set_title(f'Top {top_n} variables más importantes (XGBoost)')
ax.set_xlabel('Importancia relativa')
plt.tight_layout()
plt.show()

print("Top 10 variables:")
print(importancias.head(10).round(4))

energy_v2 = pd.read_parquet('../data/energy_silver.parquet/')
cols_energia = ['time', 'generation hydro pumped storage consumption',
                'forecast solar day ahead', 'forecast wind onshore day ahead',
                'total load forecast', 'total load actual', 'price actual', 'price day ahead']
energy_v2 = energy_v2[[c for c in cols_energia if c in energy_v2.columns]].copy()

weather_v2 = pd.read_parquet('../data/weather_features_silver.parquet').drop(columns=['weather_description'])
weather_v2['city_name'] = weather_v2['city_name'].str.strip()
weather_v2['peso'] = weather_v2['city_name'].map(pesos_ciudades)

ang = np.radians((90 - weather_v2['wind_deg']) % 360)
weather_v2['wind_x'] = weather_v2['wind_speed'] * np.cos(ang)

# Agregación climática ponderada por ciudad
cols_clima_v2 = ['temp', 'humidity', 'pressure', 'wind_speed', 'wind_x', 'clouds_all', 'rain_1h']

def _pond_v2(g):
    w = g['peso'].values
    return pd.Series({c: np.average(g[c], weights=w) for c in cols_clima_v2})

# include_groups=False evita el FutureWarning de pandas y no mete 'time' en el grupo
weather_nac_v2 = weather_v2.groupby('time').apply(_pond_v2, include_groups=False).reset_index()

# Merge + índice temporal
dfv = pd.merge(weather_nac_v2, energy_v2, on='time', how='outer').sort_values('time').reset_index(drop=True)
dfv['time'] = pd.to_datetime(dfv['time'], utc=True).dt.tz_localize(None)
dfv = dfv.set_index('time')

# Serie horaria continua y ordenada
dfv = dfv.asfreq('h')

# Interpolación SOLO de exógenas; el precio NO se interpola
precio_raw = dfv['price actual'].copy()
exog = [c for c in dfv.select_dtypes(include=[np.number]).columns if c != 'price actual']
dfv[exog] = dfv[exog].interpolate('linear', limit_direction='both')

# festivos + víspera
dfv['is_holiday'] = dfv.index.normalize().isin(festivos_dt).astype(int)
dfv['is_holiday_eve'] = dfv.index.normalize().isin(festivos_dt - pd.Timedelta(days=1)).astype(int)
dfv['weekend'] = (dfv.index.dayofweek >= 5).astype(int)

# Variables cíclicas
for col, periodo in {'hour': 24, 'dayofweek': 7, 'month': 12}.items():
    v = getattr(dfv.index, col)
    dfv[f'{col}_sin'] = np.sin(2 * np.pi * v / periodo)
    dfv[f'{col}_cos'] = np.cos(2 * np.pi * v / periodo)

# Demanda residual prevista
dfv['residual_load'] = dfv['total load forecast'] - dfv['forecast solar day ahead'] - dfv['forecast wind onshore day ahead']
dfv['renew_ratio']   = (dfv['forecast solar day ahead'] + dfv['forecast wind onshore day ahead']) / dfv['total load forecast']

# Es el mínimo limpio (deja 1 h de margen sobre el corte).
GAP = 12
pg = precio_raw.shift(GAP)  # precio más reciente garantizado al cierre de las 12:00

dfv['price_roll_mean_24h']  = pg.rolling(24).mean()
dfv['price_roll_mean_168h'] = pg.rolling(168).mean()
dfv['price_roll_std_24h']   = pg.rolling(24).std()
dfv['price_roll_min_24h']   = pg.rolling(24).min()
dfv['price_roll_max_24h']   = pg.rolling(24).max()

# Lags de precio. shift(k) sobre t = (k + 24) h antes del target.
for lag in [12, 36, 168]:
    dfv[f'price actual_lag_{lag}h'] = precio_raw.shift(lag)

# Predictores alineados t+24
dfv['load_fc_target']  = dfv['total load forecast'].shift(-24)
dfv['solar_fc_target'] = dfv['forecast solar day ahead'].shift(-24)
dfv['wind_fc_target']  = dfv['forecast wind onshore day ahead'].shift(-24)
dfv['resid_fc_target'] = dfv['load_fc_target'] - dfv['solar_fc_target'] - dfv['wind_fc_target']

# Clima (shift -24).
for c in ['temp', 'clouds_all', 'humidity', 'wind_speed', 'wind_x', 'pressure']:
    dfv[f'{c}_target'] = dfv[c].shift(-24)

# Target
dfv['target_price_24h'] = precio_raw.shift(-24)

cols_fuera = ['target_price_24h', 'price actual', 'total load actual',
              'price day ahead', 'pda_target',
              'generation hydro pumped storage consumption',
              'temp', 'humidity', 'pressure', 'wind_speed', 'wind_x',
              'clouds_all', 'rain_1h']

# dropna SOLO sobre lo que el modelo usa (X + target), como dice el comentario
X_cols = [c for c in dfv.columns if c not in cols_fuera]
dfv = dfv.dropna(subset=[*X_cols, 'target_price_24h'])

X_v2 = dfv[X_cols]
y_v2 = dfv['target_price_24h']

print(f"Dataset V2 listo: {dfv.shape} | nº features en X: {X_v2.shape[1]}")

tr_v2 = dfv.index < cutoff_date
te_v2 = dfv.index >= cutoff_date
X_dev_v2, y_dev_v2 = X_v2[tr_v2], y_v2[tr_v2]
X_test_v2, y_test_v2 = X_v2[te_v2], y_v2[te_v2]

X_tr2, X_val2, y_tr2, y_val2 = train_test_split(
    X_dev_v2, y_dev_v2, test_size=0.2, random_state=42, shuffle=False)

num_feat_v2 = X_tr2.select_dtypes(include=np.number).columns.tolist()
pre_v2 = ColumnTransformer([('num', StandardScaler(), num_feat_v2)])
pre_v2.fit(X_tr2)

xgb_v2 = xgb.XGBRegressor(
    objective='reg:squarederror',
    random_state=32,
    n_estimators=15000,
    early_stopping_rounds=150,
    max_depth=4,
    min_child_weight=5,
    learning_rate=0.003,
    reg_alpha=1.0,
    reg_lambda=1.0,
    subsample=0.8,
    colsample_bytree=0.8,
    tree_method='hist', device='cuda')

print("Entrenando Modelo V2 ...")
xgb_v2.fit(pre_v2.transform(X_tr2), y_tr2,
           eval_set=[(pre_v2.transform(X_val2), y_val2)], verbose=500)
print(f"Árbol óptimo: {xgb_v2.best_iteration}")

y_pred_v2 = xgb_v2.predict(pre_v2.transform(X_test_v2))
rmse_v2 = np.sqrt(mean_squared_error(y_test_v2, y_pred_v2))
mae_v2 = mean_absolute_error(y_test_v2, y_pred_v2)
r2_v2 = r2_score(y_test_v2, y_pred_v2)

print(f"RMSE: {rmse_v2:.2f} €/MWh  |  MAE: {mae_v2:.2f} €/MWh  |  R²: {r2_v2:.4f}")
print(f"(Modelo base V1 era R² ≈ 0.72)")

# Comparativa real vs predicción (primer mes de 2018)
plt.figure(figsize=(16, 6))
plt.plot(y_test_v2.index[:720], y_test_v2[:720], label='Precio Real',
         color=viz.NOBLUE[0], linewidth=1.5)
plt.plot(y_test_v2.index[:720], y_pred_v2[:720], label='Predicción V2 (purista)',
         color=viz.NORANGE[0], linestyle='--', linewidth=1.5)
plt.title(f'Modelo V2 purista: Real vs Predicción — Enero 2018  (R²={r2_v2:.3f})', fontsize=14)
plt.xlabel('Fecha')
plt.ylabel('Precio (€/MWh)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()