3. Análisis exploratorio de datos

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as st
import statsmodels.api as sm
import json, sys

sys.path.insert(0, ".")
from src.config import (SEED, DATOS_CRUDOS, DATOS_PROC, FIGURAS, ESTACION_OBJETIVO,
                        FECHA_CORTE_TEST, GAP_DIAS, HORIZONTE, TOLERANCIA,
                        EXCLUIR_MODELO, fijar_semillas, estilo_graficos)

fijar_semillas(SEED)
estilo_graficos()

df = pd.read_csv(DATOS_PROC / "caudales.csv", parse_dates=["fecha"])
ESTACIONES = [c for c in df.columns if c != "fecha"]
PREDICTORAS = [e for e in ESTACIONES if e != ESTACION_OBJETIVO and e not in EXCLUIR_MODELO]

print(f"{len(df):,} filas, {len(ESTACIONES)} estaciones")

corte = pd.Timestamp(FECHA_CORTE_TEST)
inicio_test = corte + pd.Timedelta(days=GAP_DIAS)

train = df[df.fecha < corte].copy()
test = df[df.fecha >= inicio_test].copy()

print(f"Train: {train.fecha.min():%Y-%m-%d} a {train.fecha.max():%Y-%m-%d}   n = {len(train):,}")
print(f"Gap:   {GAP_DIAS} días")
print(f"Test:  {test.fecha.min():%Y-%m-%d} a {test.fecha.max():%Y-%m-%d}   n = {len(test):,}")

def curva_gasto(q, qs):
    """Ajusta Qs = a * Q^b por regresion lineal en log-log."""
    m = (q > 0) & (qs > 0) & q.notna() & qs.notna()
    b, loga, r, _, _ = st.linregress(np.log10(q[m]), np.log10(qs[m]))
    return 10**loga, b, r**2, int(m.sum())

# Pares (Q, Qs) de la estacion objetivo, solo en el periodo de entrenamiento
q_cal = pd.read_csv(sorted(DATOS_CRUDOS.glob("*CALAMAR*Q_MEDIA*.csv"))[0],
                    encoding="utf-8-sig", parse_dates=["Fecha"])[["Fecha", "Q_MEDIA_D_m3_s"]]
s_cal = pd.read_csv(sorted(DATOS_CRUDOS.glob("*CALAMAR*TR_KT*.csv"))[0],
                    encoding="utf-8-sig", parse_dates=["Fecha"])[["Fecha", "TR_KT_D_QS_D_kt_dia"]]
par = q_cal.merge(s_cal, on="Fecha").dropna()
par.columns = ["fecha", "Q", "Qs"]
par = par[par.fecha < corte]

a, b, r2, n = curva_gasto(par.Q, par.Qs)
print(f"Ajuste global:  Qs = {a:.4g} * Q^{b:.4f}   R² = {r2:.4f}   n = {n:,}")

# El mismo ajuste, pero por bloques de cinco años
bloques = []
for ini, sub in par.groupby(par.fecha.dt.year // 5 * 5):
    if len(sub) < 150:
        continue
    a_i, b_i, r2_i, n_i = curva_gasto(sub.Q, sub.Qs)
    bloques.append({"bloque": f"{ini}-{ini+4}", "n": n_i, "a": a_i, "b": b_i, "R2": r2_i})

auditoria = pd.DataFrame(bloques)
print(auditoria.round(6).to_string(index=False))
print(f"\nBloques con R² > 0.999: {(auditoria.R2 > 0.999).sum()} de {len(auditoria)}")

# Dos verificaciones adicionales
from statsmodels.tsa.stattools import acf

residuo = np.log10(par.Qs) - np.log10(a * par.Q**b)
ac = acf(residuo.values, nlags=30, fft=True)
print(f"ACF del residuo:  lag1 = {ac[1]:.4f}   lag7 = {ac[7]:.4f}   lag30 = {ac[30]:.4f}")
print()

# Histeresis: un transporte medido da exponentes distintos al subir y al bajar
p = par.sort_values("fecha").copy()
p["dQ"] = p.Q.diff()
for rama, cond in [("ascenso", p.dQ > 0), ("descenso", p.dQ < 0)]:
    a_r, b_r, _, n_r = curva_gasto(p.loc[cond, "Q"], p.loc[cond, "Qs"])
    print(f"rama de {rama:9} n = {n_r:>6,}   b = {b_r:.4f}   a = {a_r:.4g}")

def descarga_efectiva(caudal, a, b, n_clases=25):
    """Analisis magnitud-frecuencia. Devuelve Q_ef y la tabla por clases."""
    x = caudal.dropna()
    x = x[x > 0]
    bordes = np.logspace(np.log10(x.min()), np.log10(x.max()), n_clases + 1)
    frec, _ = np.histogram(x, bins=bordes)
    centros = np.sqrt(bordes[:-1] * bordes[1:])
    carga = frec * (a * centros**b)
    t = pd.DataFrame({"q_inf": bordes[:-1], "q_sup": bordes[1:], "q_centro": centros,
                      "dias": frec, "frec_rel": frec / frec.sum(),
                      "carga_rel": carga / carga.sum()})
    i = t.carga_rel.idxmax()
    return t.loc[i, "q_centro"], (t.loc[i, "q_inf"], t.loc[i, "q_sup"]), t

# Curva vigente: la del ultimo bloque de la auditoria
A_OF, B_OF = auditoria.iloc[-1][["a", "b"]]
Q_EF, CLASE_EF, tabla_ef = descarga_efectiva(train[ESTACION_OBJETIVO], A_OF, B_OF)

qm = train[ESTACION_OBJETIVO].mean()
print(f"Curva utilizada:    Qs = {A_OF:.4g} * Q^{B_OF:.4f}")
print(f"Descarga efectiva:  {Q_EF:,.0f} m³/s")
print(f"Clase modal:        [{CLASE_EF[0]:,.0f} - {CLASE_EF[1]:,.0f}] m³/s")
print(f"Caudal medio:       {qm:,.0f} m³/s    Q_ef/Q_medio = {Q_EF/qm:.2f}")
print(f"Excedencia:         {100*(train[ESTACION_OBJETIVO] >= Q_EF).mean():.1f}% del tiempo")
print()
print("Robustez frente a la curva elegida:")
for _, r in auditoria.iterrows():
    qe, _, _ = descarga_efectiva(train[ESTACION_OBJETIVO], r.a, r.b)
    print(f"  {r.bloque} (b = {r.b:.3f})  ->  Q_ef = {qe:,.0f} m³/s")

fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
i_max = tabla_ef.carga_rel.idxmax()

ax[0].bar(range(len(tabla_ef)), 100 * tabla_ef.frec_rel, color="steelblue")
ax[0].axvline(i_max, color="red", ls="--", lw=1.2)
ax[0].set(title="(a) Frecuencia de ocurrencia", xlabel="Clase de caudal", ylabel="% de los días")

ax[1].bar(range(len(tabla_ef)), 100 * tabla_ef.carga_rel, color="sienna")
ax[1].axvline(i_max, color="red", ls="--", lw=1.2, label=f"$Q_{{ef}}$ = {Q_EF:,.0f} m$^3$/s")
ax[1].set(title="(b) Carga de sedimentos transportada", xlabel="Clase de caudal",
          ylabel="% de la carga total")
ax[1].legend()

fig.suptitle("Análisis magnitud-frecuencia, Calamar 1940-2009", y=1.02)
plt.tight_layout()
plt.savefig(FIGURAS / "magnitud_frecuencia.png", dpi=150)
plt.show()

LIM = (Q_EF * (1 - TOLERANCIA), Q_EF * (1 + TOLERANCIA))
df["y"] = df[ESTACION_OBJETIVO].between(*LIM).astype(float).shift(-HORIZONTE)

print(f"Banda: [{LIM[0]:,.0f} - {LIM[1]:,.0f}] m³/s")
print()
print("Prevalencia según la definición del objetivo (en entrenamiento):")
qtr = train[ESTACION_OBJETIVO]
print(f"  clase modal exacta      [{CLASE_EF[0]:>7,.0f} - {CLASE_EF[1]:>8,.0f}]   {qtr.between(*CLASE_EF).mean():>7.2%}")
for tol in (0.10, 0.15, 0.25, 0.35):
    lo, hi = Q_EF * (1 - tol), Q_EF * (1 + tol)
    print(f"  banda ±{tol:.0%}{'':14}[{lo:>7,.0f} - {hi:>8,.0f}]   {qtr.between(lo, hi).mean():>7.2%}")
print(f"  excedencia Q >= Q_ef    [{Q_EF:>7,.0f} -      inf]   {(qtr >= Q_EF).mean():>7.2%}")
print()
ytr = df.loc[df.fecha < corte, "y"]
vc = ytr.value_counts().sort_index()
print("Distribución del objetivo en entrenamiento:")
for k, v in vc.items():
    print(f"  clase {int(k)}: {v:>7,}  ({v/vc.sum():6.2%})")
print(f"  desbalance: {vc.max()/vc.min():.2f} : 1")

d = df.dropna(subset=["y"]).copy()
d["mes"] = d.fecha.dt.month
d["decada"] = d.fecha.dt.year // 10 * 10

fig, ax = plt.subplots(1, 2, figsize=(13, 3.8))
d.groupby("mes").y.mean().mul(100).plot(kind="bar", ax=ax[0], color="teal")
ax[0].set(title="(a) Prevalencia por mes", xlabel="Mes", ylabel="% de días positivos")
d.groupby("decada").y.mean().mul(100).plot(kind="bar", ax=ax[1], color="darkslategray")
ax[1].set(title="(b) Prevalencia por década", xlabel="Década", ylabel="% de días positivos")
plt.tight_layout()
plt.savefig(FIGURAS / "objetivo_tiempo.png", dpi=150)
plt.show()

filas = []
for c in ESTACIONES:
    s = train[c].dropna()
    q1, q3 = s.quantile([.25, .75])
    iqr = q3 - q1
    filas.append({"estacion": c, "n": len(s), "faltante_%": 100 * train[c].isna().mean(),
                  "media": s.mean(), "mediana": s.median(), "desv": s.std(),
                  "p01": s.quantile(.01), "p99": s.quantile(.99),
                  "min": s.min(), "max": s.max(),
                  "asimetria": s.skew(), "curtosis": s.kurtosis(),
                  "outliers_IQR_%": 100 * ((s < q1 - 1.5*iqr) | (s > q3 + 1.5*iqr)).mean()})

pd.DataFrame(filas).set_index("estacion").round(2)

fig, axes = plt.subplots(2, 4, figsize=(15, 6.5))
for ax, e in zip(axes.ravel(), ESTACIONES):
    s = train[e].dropna()
    ax.hist(s, bins=60, color="steelblue")
    ax.axvline(s.median(), color="red", ls="--", lw=1, label=f"mediana {s.median():,.0f}")
    ax.set_title(f"{e}  (asimetría {s.skew():.2f})", fontsize=10)
    ax.set_xlabel("m$^3$/s")
    ax.legend(fontsize=7)
axes.ravel()[-1].axis("off")
fig.suptitle("Caudal medio diario por estación", y=1.01)
plt.tight_layout()
plt.savefig(FIGURAS / "histogramas.png", dpi=150)
plt.show()

# Pruebas de normalidad. Con n de este orden rechazan casi siempre,
# asi que la decision se apoya en la asimetria y en los histogramas.
filas = []
for e in ESTACIONES:
    s = train[e].dropna()
    m = s.sample(min(5000, len(s)), random_state=SEED)
    filas.append({"estacion": e, "shapiro_p": st.shapiro(m)[1],
                  "jarque_bera_p": st.jarque_bera(s)[1], "asimetria": s.skew()})

pd.DataFrame(filas).round(4)

fig, ax = plt.subplots(figsize=(8.5, 4.5))
for e in ESTACIONES:
    s = np.sort(train[e].dropna().values)[::-1]
    p = np.arange(1, len(s) + 1) / (len(s) + 1)
    ax.plot(100 * p, s, lw=1.3, label=e)

ax.axhline(Q_EF, color="red", ls="--", lw=1.2, label=f"$Q_{{ef}}$ = {Q_EF:,.0f}")
ax.axhspan(*LIM, color="red", alpha=0.10)
ax.set(yscale="log", xlabel="Probabilidad de excedencia (%)", ylabel="Caudal (m$^3$/s)",
       title="Curvas de duración de caudales")
ax.legend(fontsize=8, ncol=2)
ax.grid(True, which="both", alpha=.3)
plt.tight_layout()
plt.savefig(FIGURAS / "curva_duracion.png", dpi=150)
plt.show()

corr_p = train[ESTACIONES].corr()
corr_s = train[ESTACIONES].corr(method="spearman")

fig, ax = plt.subplots(1, 2, figsize=(13.5, 4.8))
for a_, (m, t) in zip(ax, [(corr_p, "Pearson"), (corr_s, "Spearman")]):
    sns.heatmap(m, annot=True, fmt=".2f", cmap="RdYlBu_r", center=0, vmin=-1, vmax=1,
                square=True, ax=a_, cbar_kws={"shrink": .8}, annot_kws={"size": 8})
    a_.set_title(f"Correlación de {t}")
plt.tight_layout()
plt.savefig(FIGURAS / "correlaciones.png", dpi=150)
plt.show()

print("Correlación con la estación objetivo:")
print(corr_p[ESTACION_OBJETIVO].drop(ESTACION_OBJETIVO).sort_values(ascending=False).round(3).to_string())

from statsmodels.stats.outliers_influence import variance_inflation_factor

sub = train.dropna(subset=ESTACIONES)
X = sm.add_constant(sub[ESTACIONES])
vif = pd.DataFrame({"estacion": X.columns,
                    "VIF": [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]})
vif = vif[vif.estacion != "const"].sort_values("VIF", ascending=False)

print(f"Periodo con todas las estaciones: {len(sub):,} días")
print(vif.round(2).to_string(index=False))

from statsmodels.stats.multitest import multipletests
from sklearn.feature_selection import mutual_info_classif

dtr = df[df.fecha < corte].dropna(subset=["y"])

filas = []
for e in ESTACIONES:
    s = dtr[[e, "y"]].dropna()
    a_, b_ = s.loc[s.y == 1, e], s.loc[s.y == 0, e]
    u, pv = st.mannwhitneyu(a_, b_)
    filas.append({"estacion": e, "mediana_y1": a_.median(), "mediana_y0": b_.median(),
                  "p_valor": pv, "efecto_A": u / (len(a_) * len(b_)),
                  "cohen_d": (a_.mean() - b_.mean()) / np.sqrt((a_.var() + b_.var()) / 2)})

aso = pd.DataFrame(filas)
aso["p_holm"] = multipletests(aso.p_valor, method="holm")[1]

mi = dtr[ESTACIONES + ["y"]].dropna()
aso["info_mutua"] = mutual_info_classif(mi[ESTACIONES], mi.y, random_state=SEED)

print("Con n de este orden casi toda diferencia resulta significativa, por lo que")
print("se reporta el tamaño de efecto junto al p-valor, con corrección de Holm.")
aso.round(4)

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.covariance import MinCovDet

Xs = StandardScaler().fit_transform(sub[ESTACIONES])
pca = PCA().fit(Xs)

varianza = pd.DataFrame({
    "componente": [f"PC{i+1}" for i in range(len(ESTACIONES))],
    "varianza_explicada": pca.explained_variance_ratio_,
    "acumulada": pca.explained_variance_ratio_.cumsum()})

n90 = int((varianza.acumulada < 0.90).sum() + 1)
print(varianza.round(4).to_string(index=False))
print(f"\nComponentes para el 90% de la varianza: {n90} de {len(ESTACIONES)}")

pesos = pd.DataFrame(pca.components_[:3].T, index=ESTACIONES, columns=["PC1", "PC2", "PC3"])
print("\nPesos de las tres primeras componentes:")
print(pesos.round(3).to_string())

fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
ax[0].bar(range(1, len(ESTACIONES) + 1), 100 * pca.explained_variance_ratio_, color="steelblue")
ax[0].plot(range(1, len(ESTACIONES) + 1), 100 * varianza.acumulada, "o-", color="darkred")
ax[0].axhline(90, ls="--", color="gray", lw=1)
ax[0].set(xlabel="Componente principal", ylabel="% de varianza", title="(a) Varianza explicada")

sns.heatmap(pesos, annot=True, fmt=".2f", cmap="RdBu_r", center=0, ax=ax[1], cbar_kws={"shrink": .8})
ax[1].set_title("(b) Pesos de las componentes")
plt.tight_layout()
plt.savefig(FIGURAS / "pca.png", dpi=150)
plt.show()

# Dias en que las estaciones se comportan de forma incoherente entre si
mcd = MinCovDet(random_state=SEED).fit(sub[ESTACIONES])
d2 = mcd.mahalanobis(sub[ESTACIONES])
umbral = st.chi2.ppf(0.999, df=len(ESTACIONES))

atipicos = sub.loc[d2 > umbral, ["fecha"] + ESTACIONES].copy()
atipicos["mahalanobis"] = d2[d2 > umbral]

print(f"Umbral chi² (p = 0.001, gl = {len(ESTACIONES)}): {umbral:.1f}")
print(f"Días atípicos: {len(atipicos):,} de {len(sub):,} ({len(atipicos)/len(sub):.2%})")
atipicos.nlargest(8, "mahalanobis").round(0)

from sklearn.metrics import roc_auc_score

filas = []
for e in ESTACIONES:
    s = dtr[[e, "y"]].dropna()
    auc = roc_auc_score(s.y, s[e])
    filas.append({"estacion": e, "AUC_univariado": max(auc, 1 - auc)})

fuga = pd.DataFrame(filas).sort_values("AUC_univariado", ascending=False)
print(fuga.round(4).to_string(index=False))
print(f"\nVariables con AUC > 0.95: {(fuga.AUC_univariado > 0.95).sum()}")
print(f"Fechas comunes entre train y test: {len(set(train.fecha) & set(test.fecha))}")

f = train.fecha
print(f"Cobertura: {f.min():%Y-%m-%d} a {f.max():%Y-%m-%d}  ({(f.max()-f.min()).days/365.25:.1f} años)")
print(f"Frecuencia: diaria regular    Fechas duplicadas: {f.duplicated().sum()}")
print()

s = train[ESTACION_OBJETIVO]
sin_dato = train.loc[s.isna(), "fecha"]
print(f"Días sin dato en la estación objetivo: {len(sin_dato):,} ({s.isna().mean():.2%})")

if len(sin_dato):
    g = sin_dato.diff().dt.days.ne(1).cumsum()
    huecos = sin_dato.groupby(g).agg(["size", "min", "max"])
    huecos.columns = ["dias", "desde", "hasta"]
    huecos = huecos.sort_values("dias", ascending=False)
    print(f"Número de huecos: {len(huecos)}")
    print("\nLos ocho más largos:")
    for _, r in huecos.head(8).iterrows():
        print(f"  {r.dias:>4} días   {r.desde:%Y-%m-%d} a {r.hasta:%Y-%m-%d}")

fig, ax = plt.subplots(2, 1, figsize=(13, 6.5))

ax[0].plot(train.fecha, train[ESTACION_OBJETIVO], lw=.35, color="steelblue")
ax[0].axhspan(*LIM, color="red", alpha=.13, label="banda de descarga efectiva")
ax[0].set(ylabel="Caudal (m$^3$/s)", title="(a) Serie completa")
ax[0].legend(fontsize=8)

mensual = train.set_index("fecha")[ESTACION_OBJETIVO].resample("MS").mean()
sns.boxplot(x=mensual.index.month, y=mensual.values, ax=ax[1], color="lightsteelblue")
ax[1].axhline(Q_EF, color="red", ls="--", lw=1.2, label="$Q_{ef}$")
ax[1].set(xlabel="Mes", ylabel="Caudal medio mensual (m$^3$/s)",
          title="(b) Estacionalidad por ciclo anual")
ax[1].legend(fontsize=8)

plt.tight_layout()
plt.savefig(FIGURAS / "serie_estacionalidad.png", dpi=150)
plt.show()

from statsmodels.tsa.seasonal import STL

serie = mensual.interpolate(limit=3).dropna()
stl = STL(serie, period=12, robust=True).fit()

fig = stl.plot()
fig.set_size_inches(11, 7)
fig.suptitle("Descomposición STL del caudal mensual", y=1.01)
plt.tight_layout()
plt.savefig(FIGURAS / "stl.png", dpi=150)
plt.show()

total = np.var(stl.trend) + np.var(stl.seasonal) + np.var(stl.resid)
print(f"Tendencia {np.var(stl.trend)/total:.1%} | "
      f"Estacionalidad {np.var(stl.seasonal)/total:.1%} | "
      f"Residuo {np.var(stl.resid)/total:.1%}")

from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

sd = train[ESTACION_OBJETIVO].dropna()
p_adf = adfuller(sd, autolag="AIC")[1]
p_kpss = kpss(sd, regression="c", nlags="auto")[1]

# Las dos pruebas tienen hipotesis nulas opuestas
print(f"ADF   H0: no estacionaria   p = {p_adf:.4f}")
print(f"KPSS  H0: sí estacionaria   p = {p_kpss:.4f}")
print("Conclusión:", "la serie es estacionaria en nivel" if p_adf < .05 and p_kpss >= .05
      else "las pruebas discrepan, interpretar con cautela")

fig, ax = plt.subplots(1, 2, figsize=(13, 3.6))
plot_acf(sd, lags=90, ax=ax[0])
ax[0].set_title("ACF del caudal diario")
plot_pacf(sd, lags=40, ax=ax[1], method="ywm")
ax[1].set_title("PACF del caudal diario")
plt.tight_layout()
plt.savefig(FIGURAS / "acf_pacf.png", dpi=150)
plt.show()

ac = acf(sd.values, nlags=60, fft=True)
print(f"\nACF: lag1 = {ac[1]:.4f} | lag7 = {ac[7]:.4f} | lag15 = {ac[15]:.4f} | lag30 = {ac[30]:.4f}")

# Tiempo de transito: rezago que maximiza la correlacion con la estacion objetivo
TRANSITOS, curvas = {}, {}
filas = []
for e in ESTACIONES:
    if e == ESTACION_OBJETIVO:
        continue
    r = [(k, train[e].shift(k).corr(train[ESTACION_OBJETIVO])) for k in range(31)]
    t = pd.DataFrame(r, columns=["rezago", "correlacion"])
    k_opt = int(t.loc[t.correlacion.idxmax(), "rezago"])
    TRANSITOS[e] = k_opt
    curvas[e] = t
    filas.append({"estacion": e, "transito_dias": k_opt,
                  "r_con_desfase": t.correlacion.max(),
                  "r_sin_desfase": t.loc[0, "correlacion"],
                  "ganancia": t.correlacion.max() - t.loc[0, "correlacion"]})

pd.DataFrame(filas).sort_values("transito_dias", ascending=False).round(4)

fig, ax = plt.subplots(figsize=(9, 4.5))
for e, t in curvas.items():
    ax.plot(t.rezago, t.correlacion, lw=1.5, label=f"{e} ({TRANSITOS[e]} d)")
    ax.plot(TRANSITOS[e], t.correlacion.max(), "o", ms=6)

ax.set(xlabel="Rezago aplicado a la estación de aguas arriba (días)",
       ylabel=f"Correlación con {ESTACION_OBJETIVO}",
       title="Correlación cruzada y tiempo de tránsito")
ax.legend(fontsize=8)
ax.grid(alpha=.3)
plt.tight_layout()
plt.savefig(FIGURAS / "transito.png", dpi=150)
plt.show()

# Los rezagos se estimaron sin usar informacion geografica.
# Contrastarlos con la ubicacion real es una verificacion externa.
ubic = pd.read_csv(DATOS_CRUDOS.parent / "UBICACION_ESTACIONES_DESCARGADAS_IDEAM.csv",
                   sep=";", encoding="utf-8-sig")
ubic.columns = ["estacion", "codigo", "departamento", "municipio", "lat", "lon", "altitud"]

MAPA = {"ARRANCAPLUMAS - AUT": "arrancaplumas", "PUERTO SALGAR - AUT": "pto_salgar",
        "TRES CRUCES": "tres_cruces", "COYONGAL": "coyongal", "EL BANCO - AUT": "el_banco",
        "SAN ROQUE": "san_roque", "CALAMAR": "calamar"}
ubic["clave"] = ubic.estacion.map(MAPA)
ubic["transito"] = ubic.clave.map({**TRANSITOS, ESTACION_OBJETIVO: 0})

g = ubic.dropna(subset=["transito"]).sort_values("transito", ascending=False)
print(g[["estacion", "municipio", "altitud", "lat", "transito"]].to_string(index=False))
print()
for col, nom in [("altitud", "altitud"), ("lat", "latitud")]:
    r, p = st.pearsonr(g[col], g.transito)
    rs, _ = st.spearmanr(g[col], g.transito)
    print(f"tránsito vs {nom:8}  Pearson r = {r:+.4f} (p = {p:.4f})   Spearman = {rs:+.4f}")

fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))

ax[0].scatter(g.altitud, g.transito, s=70, color="sienna", zorder=3)
for _, r in g.iterrows():
    ax[0].annotate(r.clave, (r.altitud, r.transito), fontsize=7.5,
                   xytext=(5, 4), textcoords="offset points")
ax[0].set(xlabel="Altitud (m s.n.m.)", ylabel="Tiempo de tránsito (días)",
          title="(a) Tránsito frente a altitud")
ax[0].grid(alpha=.3)

ax[1].scatter(g.lon, g.lat, s=90, c=g.transito, cmap="viridis_r", zorder=3)
for _, r in g.iterrows():
    ax[1].annotate(f"{r.clave}\n{int(r.transito)} d", (r.lon, r.lat), fontsize=7,
                   xytext=(6, -3), textcoords="offset points")
ax[1].set(xlabel="Longitud (°W)", ylabel="Latitud (°N)",
          title="(b) Disposición geográfica")
ax[1].grid(alpha=.3)

plt.tight_layout()
plt.savefig(FIGURAS / "validacion_geografica.png", dpi=150)
plt.show()

# Deriva entre decadas
d_ = train.copy()
d_["decada"] = d_.fecha.dt.year // 10 * 10
decadas = sorted(d_.decada.unique())
base = decadas[0]

filas = []
for dd in decadas[1:]:
    a_ = d_.loc[d_.decada == base, ESTACION_OBJETIVO].dropna()
    b_ = d_.loc[d_.decada == dd, ESTACION_OBJETIVO].dropna()
    if len(b_) < 300:
        continue
    ks, pv = st.ks_2samp(a_, b_)
    filas.append({"decada": dd, "n": len(b_), "media": b_.mean(), "KS": ks, "p_valor": pv})

print("Prueba de Kolmogorov-Smirnov frente a la primera década:")
print(pd.DataFrame(filas).round(4).to_string(index=False))
print()
dy = df[df.fecha < corte].dropna(subset=["y"]).copy()
dy["decada"] = dy.fecha.dt.year // 10 * 10
print("Prevalencia del objetivo por década:")
print(dy.groupby("decada").y.agg(["mean", "size"]).round(4).to_string())

def construir_features(d):
    """Predictores. Todos usan informacion hasta t; el objetivo esta en t+h."""
    x = d.sort_values("fecha").copy()
    o = ESTACION_OBJETIVO

    for k in (1, 2, 3, 5, 7, 10, 15, 30, 60):
        x[f"{o}_lag{k}"] = x[o].shift(k)

    for w in (7, 30, 90, 365):
        base = x[o].shift(1)          # shift(1) excluye el dia actual
        x[f"{o}_media{w}"] = base.rolling(w, min_periods=max(2, w//3)).mean()
        x[f"{o}_desv{w}"] = base.rolling(w, min_periods=max(2, w//3)).std()

    x[f"{o}_d1"] = x[o].diff(1)
    x[f"{o}_d7"] = x[o].diff(7)

    # Cada estacion entra con su tiempo de transito
    for e in PREDICTORAS:
        k = TRANSITOS.get(e, 0)
        x[f"{e}_t{k}"] = x[e].shift(k)
        x[f"{e}_d3"] = x[e].shift(k).diff(3)
        x[f"{e}_rel"] = x[e].shift(k) / x[e].shift(1).rolling(365, min_periods=120).mean()

    # Distancia a Q_ef: convierte la banda en un problema monotono
    x["dist_qef"] = (x[o] - Q_EF).abs()
    x["dist_qef_rel"] = x["dist_qef"] / Q_EF
    x[f"{o}_cuad"] = x[o] ** 2
    for k in (7, 15):
        x[f"dist_qef_lag{k}"] = (x[f"{o}_lag{k}"] - Q_EF).abs()

    # Regimen bimodal: primer y segundo armonico
    doy = x.fecha.dt.dayofyear
    x["sen_anual"] = np.sin(2 * np.pi * doy / 365.25)
    x["cos_anual"] = np.cos(2 * np.pi * doy / 365.25)
    x["sen_semi"] = np.sin(4 * np.pi * doy / 365.25)
    x["cos_semi"] = np.cos(4 * np.pi * doy / 365.25)
    return x

dfx = construir_features(df)

BASE_COLS = ([c for c in dfx.columns if c.startswith(f"{ESTACION_OBJETIVO}_")
              and not c.endswith("_cuad")]
             + [c for e in PREDICTORAS for c in dfx.columns if c.startswith(e + "_")]
             + ["sen_anual", "cos_anual", "sen_semi", "cos_semi"])
DOM_COLS = ["dist_qef", "dist_qef_rel", f"{ESTACION_OBJETIVO}_cuad",
            "dist_qef_lag7", "dist_qef_lag15"]
TODAS = BASE_COLS + DOM_COLS

print(f"Predictores: {len(TODAS)}")
print(f"  derivados del caudal objetivo y calendario: {len([c for c in BASE_COLS if not any(c.startswith(e) for e in PREDICTORAS)])}")
print(f"  derivados de estaciones aguas arriba:       {len([c for c in BASE_COLS if any(c.startswith(e) for e in PREDICTORAS)])}")
print(f"  de dominio:                                 {len(DOM_COLS)}")

# Verificacion anti-fuga sobre las variables construidas
dd = dfx[dfx.fecha < corte].dropna(subset=["y"])
filas = []
for c in TODAS:
    s = dd[[c, "y"]].dropna()
    if len(s) < 500:
        continue
    auc = roc_auc_score(s.y, s[c])
    filas.append({"variable": c, "AUC": max(auc, 1 - auc)})

vf = pd.DataFrame(filas).sort_values("AUC", ascending=False)
print("Las diez variables con mayor AUC univariado:")
print(vf.head(10).round(4).to_string(index=False))
print()
print("dist_qef tiene AUC alto por construcción: es la distancia al centro de la banda")
print("evaluada en t. Usa solo información de t y el objetivo está en t+15.")

# El capitulo 4 parte de estos archivos
dfx[["fecha", "y"] + TODAS + [ESTACION_OBJETIVO]].to_csv(DATOS_PROC / "features.csv", index=False)

parametros = {
    "Q_EF": float(Q_EF), "LIM": [float(LIM[0]), float(LIM[1])],
    "A_OF": float(A_OF), "B_OF": float(B_OF),
    "TRANSITOS": TRANSITOS, "BASE_COLS": BASE_COLS, "DOM_COLS": DOM_COLS,
}
with open(DATOS_PROC / "parametros.json", "w", encoding="utf-8") as f:
    json.dump(parametros, f, indent=2, ensure_ascii=False)

print(f"Guardados: features.csv ({len(dfx):,} filas) y parametros.json")