Recomendador de productos financieros

"""#9 Sistema de recomendación de productos financieros (next best offer).

Datos: Santander Product Recommendation (Kaggle; requiere aceptar las reglas y KAGGLE_API_TOKEN).
~950 mil clientes con la tenencia mensual de 24 productos (ene-2015 a may-2016). La pregunta de negocio:
¿qué producto NUEVO contratará cada cliente el próximo mes? Se evalúa en mayo 2016 (mes nunca visto),
entrenando con las altas de marzo y abril 2016.
Modelos: popularidad, filtrado colaborativo ítem-ítem, LightGBM multiclase y LightGBM LambdaRank; MAP@7.
"""
import os
import subprocess
import sys
import time
import zipfile
from pathlib import Path

import duckdb
import numpy as np
import pandas as pd
from lightgbm import LGBMClassifier, LGBMRanker

from common import RAW, SEED, bars, export, heatmap, line, table
from validacion import model_card

DIR = RAW / "kaggle_santander"
PRODUCTOS = {
    "ind_ahor_fin_ult1": "Cuenta de ahorro", "ind_aval_fin_ult1": "Garantías", "ind_cco_fin_ult1": "Cuenta corriente",
    "ind_cder_fin_ult1": "Cuenta derivados", "ind_cno_fin_ult1": "Cuenta nómina", "ind_ctju_fin_ult1": "Cuenta junior",
    "ind_ctma_fin_ult1": "Cuenta Más", "ind_ctop_fin_ult1": "Cuenta particular", "ind_ctpp_fin_ult1": "Cuenta particular plus",
    "ind_deco_fin_ult1": "Depósito corto plazo", "ind_deme_fin_ult1": "Depósito mediano plazo",
    "ind_dela_fin_ult1": "Depósito largo plazo", "ind_ecue_fin_ult1": "e-cuenta", "ind_fond_fin_ult1": "Fondos",
    "ind_hip_fin_ult1": "Hipoteca", "ind_plan_fin_ult1": "Plan de pensiones", "ind_pres_fin_ult1": "Préstamo",
    "ind_reca_fin_ult1": "Pago de impuestos", "ind_tjcr_fin_ult1": "Tarjeta de crédito", "ind_valo_fin_ult1": "Valores",
    "ind_viv_fin_ult1": "Cuenta hogar", "ind_nomina_ult1": "Domiciliación de nómina",
    "ind_nom_pens_ult1": "Domiciliación de pensión", "ind_recibo_ult1": "Domiciliación de recibos",
}
P = list(PRODUCTOS)
NOMBRES = list(PRODUCTOS.values())
MESES = ["2015-11-28", "2015-12-28", "2016-01-28", "2016-02-28", "2016-03-28", "2016-04-28", "2016-05-28"]
ENTRENO, PRUEBA = ["2016-03-28", "2016-04-28"], "2016-05-28"
K = 7


def datos() -> pd.DataFrame:
    cache = DIR / "santander_7m.parquet"
    if cache.exists():
        return pd.read_parquet(cache)
    csv = DIR / "train_ver2.csv"
    if not csv.exists():
        if not (DIR / "train_ver2.csv.zip").exists():
            subprocess.run([str(Path(sys.executable).parent / "kaggle"), "competitions", "download",
                            "santander-product-recommendation", "-f", "train_ver2.csv.zip", "-p", str(DIR)], check=True)
        with zipfile.ZipFile(DIR / "train_ver2.csv.zip") as z:
            z.extractall(DIR)
    prods = ", ".join(f"COALESCE(TRY_CAST({p} AS INTEGER), 0)::TINYINT AS {p}" for p in P)
    meses = ", ".join(f"DATE '{m}'" for m in MESES)
    df = duckdb.connect().execute(f"""
        SELECT fecha_dato::DATE AS fecha, ncodpers::INTEGER AS cliente,
               TRY_CAST(TRIM(age) AS INTEGER) AS edad, TRY_CAST(TRIM(antiguedad) AS INTEGER) AS antiguedad,
               TRY_CAST(renta AS DOUBLE) AS renta, segmento, canal_entrada AS canal,
               TRY_CAST(ind_actividad_cliente AS INTEGER) AS activo, TRY_CAST(cod_prov AS INTEGER) AS provincia,
               ind_empleado AS empleado, TRY_CAST(ind_nuevo AS INTEGER) AS nuevo, tiprel_1mes AS relacion,
               {prods}
        FROM read_csv_auto('{csv}', all_varchar=true)
        WHERE fecha_dato::DATE IN ({meses})
    """).df()
    df.to_parquet(cache)
    return df


def tabla_mes(df: pd.DataFrame, mes: str) -> pd.DataFrame:
    """Clientes del mes con: tenencia del mes anterior, rezagos (2–4 meses), perfil y las ALTAS del mes (objetivo)."""
    i = MESES.index(mes)
    actual = df[df.fecha == pd.Timestamp(mes)].set_index("cliente")
    prev = df[df.fecha == pd.Timestamp(MESES[i - 1])].set_index("cliente")
    # Perfil y tenencia del mes ANTERIOR: es lo que se conoce al decidir la oferta (el perfil del mismo mes,
    # p. ej. "cliente activo", puede cambiar justamente por la contratación y filtraría el resultado)
    base = prev.drop(columns=P + ["fecha"]).join(prev[P].add_suffix("_t1"))
    base = base[base.index.isin(actual.index)]
    for r in (2, 3, 4):
        lag = df[df.fecha == pd.Timestamp(MESES[i - r])].set_index("cliente")[P]
        base = base.join(lag.add_suffix(f"_t{r}"), how="left")
    for p in P:   # ¿hace cuánto lo dejó? (1 si lo tuvo en algún rezago pero no en t-1)
        base[f"{p}_tuvo_antes"] = ((base[[f"{p}_t{r}" for r in (2, 3, 4)]].max(axis=1) == 1) &
                                   (base[f"{p}_t1"] == 0)).astype(int)
    base["n_productos_t1"] = base[[f"{p}_t1" for p in P]].sum(axis=1)
    altas = (actual[P].reindex(base.index).fillna(0).values == 1) & (base[[f"{p}_t1" for p in P]].values == 0)
    return base, pd.DataFrame(altas.astype(int), index=base.index, columns=P)


def ap_at_k(reales: np.ndarray, ranking: np.ndarray, k: int = K) -> float:
    """Average precision @k de un cliente: reales = vector binario de altas, ranking = índices ordenados."""
    if reales.sum() == 0:
        return 0.0
    aciertos, suma = 0, 0.0
    for pos, j in enumerate(ranking[:k]):
        if reales[j]:
            aciertos += 1
            suma += aciertos / (pos + 1)
    return suma / min(reales.sum(), k)


def evaluar(scores: np.ndarray, poseidos: np.ndarray, altas: np.ndarray) -> np.ndarray:
    s = np.where(poseidos == 1, -np.inf, scores)            # nunca recomendar lo que ya tiene
    rank = np.argsort(-s, axis=1)[:, :K]
    return np.array([ap_at_k(a, r) for a, r in zip(altas, rank)])


def run():
    if not ((DIR / "santander_7m.parquet").exists() or (DIR / "train_ver2.csv.zip").exists()
            or os.environ.get("KAGGLE_API_TOKEN") or (Path.home() / ".kaggle").exists()):
        print("  · p09: sin datos de Kaggle ni KAGGLE_API_TOKEN; se omite")
        return None
    df = datos()
    for c in ["segmento", "canal", "empleado", "relacion"]:
        df[c] = df[c].fillna("NA").astype("category")

    tr_partes = [tabla_mes(df, m) for m in ENTRENO]
    X_te, Y_te = tabla_mes(df, PRUEBA)
    X_tr = pd.concat([x for x, _ in tr_partes])
    Y_tr = pd.concat([y for _, y in tr_partes])
    pos_te = X_te[[f"{p}_t1" for p in P]].values
    con_alta = Y_te.values.sum(axis=1) > 0
    n_clientes, n_alta = len(X_te), int(con_alta.sum())
    Ya, posa = Y_te.values[con_alta], pos_te[con_alta]
    resultados, ap_por_modelo, tiempos = [], {}, {}

    # 1) Popularidad: frecuencia de altas en entrenamiento
    pop = Y_tr.values.mean(axis=0)
    ap_por_modelo["Popularidad"] = evaluar(np.tile(pop, (n_alta, 1)), posa, Ya)

    # 2) Filtrado colaborativo ítem-ítem: P(alta j | posee i), suavizado y combinado sobre lo que posee
    pos_tr = X_tr[[f"{p}_t1" for p in P]].values
    co = (pos_tr.T @ Y_tr.values + pop * 20) / (pos_tr.sum(axis=0)[:, None] + 20)   # 24 × 24
    sin_nada = (posa.sum(axis=1, keepdims=True) == 0)
    cf = np.where(sin_nada, pop, posa @ co / np.maximum(posa.sum(axis=1, keepdims=True), 1))
    ap_por_modelo["Filtrado colaborativo ítem-ítem"] = evaluar(cf, posa, Ya)

    # 3) LightGBM multiclase: una fila por (cliente, producto dado de alta), objetivo = producto
    feats = list(X_tr.columns)
    filas = np.repeat(np.arange(len(X_tr)), Y_tr.values.sum(axis=1))
    clase = np.concatenate([np.flatnonzero(r) for r in Y_tr.values if r.any()])
    t = time.time()
    mc = LGBMClassifier(objective="multiclass", n_estimators=300, learning_rate=0.05, num_leaves=63,
                        min_child_samples=50, colsample_bytree=0.7, random_state=SEED, verbose=-1)
    mc.fit(X_tr.iloc[filas][feats], clase)
    tiempos["LightGBM multiclase"] = time.time() - t
    proba = np.zeros((n_alta, len(P)))
    proba[:, mc.classes_] = mc.predict_proba(X_te[con_alta][feats])
    ap_por_modelo["LightGBM multiclase"] = evaluar(proba, posa, Ya)

    # 4) LightGBM LambdaRank: pares (cliente, producto candidato no poseído), optimiza directamente el orden
    rng = np.random.default_rng(SEED)
    idx_tr = np.flatnonzero(Y_tr.values.sum(axis=1) > 0)
    idx_tr = rng.choice(idx_tr, min(len(idx_tr), 60_000), replace=False)

    def pares(X, Y, poseidos, idx):
        reps = []
        for j, p in enumerate(P):
            ok = poseidos[idx, j] == 0
            d = X.iloc[idx[ok]][feats].copy()
            d["producto"] = j
            d["pop_producto"] = pop[j]
            d["cf_producto"] = (poseidos[idx[ok]] @ co[:, j]) / np.maximum(poseidos[idx[ok]].sum(axis=1), 1)
            d["y"] = Y[idx[ok], j]
            d["grupo"] = idx[ok]
            reps.append(d)
        out = pd.concat(reps).sort_values("grupo", kind="stable")
        out["producto"] = out["producto"].astype("category")
        return out
    D_tr = pares(X_tr, Y_tr.values, pos_tr, idx_tr)
    t = time.time()
    rk = LGBMRanker(objective="lambdarank", n_estimators=300, learning_rate=0.05, num_leaves=63,
                    min_child_samples=50, colsample_bytree=0.7, random_state=SEED, verbose=-1,
                    lambdarank_truncation_level=K)
    rk.fit(D_tr.drop(columns=["y", "grupo"]), D_tr.y, group=D_tr.groupby("grupo", sort=False).size().values)
    tiempos["LightGBM LambdaRank"] = time.time() - t
    idx_te = np.flatnonzero(con_alta)
    D_te = pares(X_te, Y_te.values, pos_te, idx_te)
    D_te["s"] = rk.predict(D_te.drop(columns=["y", "grupo"]))
    S = np.full((n_clientes, len(P)), -np.inf)
    S[D_te.grupo.values, D_te.producto.astype(int).values] = D_te.s.values
    ap_por_modelo["LightGBM LambdaRank"] = evaluar(S[con_alta], posa, Ya)

    # Métricas con IC bootstrap (por cliente) y comparación pareada contra el mejor
    boots = [rng.integers(0, n_alta, n_alta) for _ in range(500)]
    for nombre, ap in ap_por_modelo.items():
        bs = [ap[i].mean() for i in boots]
        resultados.append({"Modelo": nombre, "MAP@7 (clientes con alta)": round(ap.mean(), 4),
                           "IC 95 %": f"{np.percentile(bs, 2.5):.4f}–{np.percentile(bs, 97.5):.4f}",
                           "MAP@7 estilo Kaggle (todos)": round(ap.sum() / n_clientes, 5),
                           "Entrenamiento (s)": round(tiempos.get(nombre, 0), 1)})
    res = pd.DataFrame(resultados)
    mejor = res.loc[res["MAP@7 (clientes con alta)"].idxmax(), "Modelo"]
    pareadas = []
    for otro in ap_por_modelo:
        if otro == mejor:
            continue
        d = ap_por_modelo[mejor] - ap_por_modelo[otro]
        bs = [d[i].mean() for i in boots]
        pareadas.append({"Comparación": f"{mejor} − {otro}", "Δ MAP@7": round(d.mean(), 4),
                         "IC 95 %": f"{np.percentile(bs, 2.5):+.4f} a {np.percentile(bs, 97.5):+.4f}"})
    pareadas = pd.DataFrame(pareadas)
    seg_row = pareadas.loc[pareadas["Δ MAP@7"].idxmin()]
    segundo = seg_row["Comparación"].split(" − ")[1]
    lo_seg = float(seg_row["IC 95 %"].split(" a ")[0])
    empate = lo_seg <= 0

    # Diagnóstico del mejor: precisión de la 1.ª recomendación por producto y curva de aciertos por posición
    S_best = {"LightGBM multiclase": proba, "LightGBM LambdaRank": S[con_alta],
              "Filtrado colaborativo ítem-ítem": cf, "Popularidad": np.tile(pop, (n_alta, 1))}[mejor]
    s = np.where(posa == 1, -np.inf, S_best)
    top = np.argsort(-s, axis=1)[:, :K]
    hit_k = {nombre: [] for nombre in ap_por_modelo}
    for nombre, M in {"Popularidad": np.tile(pop, (n_alta, 1)), "Filtrado colaborativo ítem-ítem": cf,
                      "LightGBM multiclase": proba, "LightGBM LambdaRank": S[con_alta]}.items():
        tk = np.argsort(-np.where(posa == 1, -np.inf, M), axis=1)[:, :K]
        aciertos = np.take_along_axis(Ya, tk, axis=1)
        hit_k[nombre] = [[k, round(float((aciertos[:, :k].sum(axis=1) > 0).mean() * 100), 2)] for k in range(1, K + 1)]
    primera = top[:, 0]
    acierto1 = Ya[np.arange(n_alta), primera]
    por_prod = pd.DataFrame({"producto": [NOMBRES[j] for j in primera], "acierto": acierto1}).groupby("producto") \
        .agg(recomendado=("acierto", "size"), precision=("acierto", "mean")).sort_values("recomendado", ascending=False)
    por_prod = por_prod[por_prod.recomendado >= 50].head(10)
    altas_prod = pd.Series(Y_te.values.sum(axis=0), index=NOMBRES).sort_values(ascending=False)
    # Matriz de co-ocurrencia (colaborativo): probabilidad de alta del producto j dado que posee i (top productos)
    top_p = [P.index(p) for p in P if PRODUCTOS[p] in altas_prod.index[:8]]
    r_ = res.set_index("Modelo")
    m1 = r_.loc[mejor, "MAP@7 (clientes con alta)"]
    pop_v = r_.loc["Popularidad", "MAP@7 (clientes con alta)"]
    hit1 = dict(hit_k[mejor])[1]

    result = {
        "title": "Recomendador de productos financieros (next best offer)",
        "summary": [
            f"{df.cliente.nunique():,} clientes reales de un banco con 24 productos. Objetivo: el producto NUEVO que "
            f"contratará el próximo mes. Prueba en mayo 2016 (mes no visto): {n_alta:,} de {n_clientes:,} clientes "
            f"({n_alta / n_clientes:.1%}) dieron de alta algo; el resto no se puede acertar.",
            f"Mejor modelo: {mejor}, MAP@7 {m1:.3f} (IC {r_.loc[mejor, 'IC 95 %']}) entre clientes con alta, vs. "
            f"popularidad {pop_v:.3f} ({m1 / pop_v - 1:+.0%}). La primera recomendación acierta en {hit1:.1f} % de los casos. "
            + (f"Empata con {segundo} (Δ {seg_row['Δ MAP@7']:+.4f}, IC {seg_row['IC 95 %']}): ante un empate conviene el "
               "más simple de operar." if empate else ""),
            f"El filtrado colaborativo puro ({r_.loc['Filtrado colaborativo ítem-ítem', 'MAP@7 (clientes con alta)']:.3f}) "
            "mejora poco la popularidad: con 24 productos el catálogo es pequeño y lo que más predice es el perfil "
            "(edad, renta, antigüedad, canal) y lo que el cliente tuvo y dejó, que el boosting sí usa.",
            f"En el estilo de la competencia (promediando también a quienes no compraron) el MAP@7 es "
            f"{r_.loc[mejor, 'MAP@7 estilo Kaggle (todos)']:.4f}; el ganador de Kaggle obtuvo ~0.031 con muchos más "
            f"meses de historia y features (y en junio, otro mes). La métrica refleja que solo {n_alta / n_clientes:.1%} de los clientes compra algo en el mes.",
            "Uso en una financiera automotriz: ordenar a quién ofrecer seguro, garantía extendida o refinanciamiento "
            "en cada contacto, con la misma lógica (perfil + tenencia + historial).",
        ],
        "metrics": [
            {"label": "MAP@7 (mejor)", "value": f"{m1:.3f}"},
            {"label": "Mejora vs. popularidad", "value": f"{m1 / pop_v - 1:+.0%}"},
            {"label": "Acierto 1.ª recomendación", "value": f"{hit1:.1f} %"},
            {"label": "Clientes evaluados", "value": f"{n_alta:,}"},
        ],
        "charts": [
            bars("map", "MAP@7 por modelo (mayo 2016, clientes con alta)", res.Modelo,
                 res["MAP@7 (clientes con alta)"].values, fmt="dec"),
            line("hit", "Clientes con al menos un acierto entre las primeras k recomendaciones", hit_k,
                 x_label="k (número de recomendaciones)", y_label="% de clientes", fmt="pct"),
            bars("altas", "Altas por producto en mayo 2016", altas_prod.index[:12], altas_prod.values[:12], fmt="int"),
            bars("precision1", f"Precisión de la 1.ª recomendación por producto ({mejor})", por_prod.index,
                 (por_prod.precision * 100).round(1).values, fmt="pct"),
            heatmap("co", "Colaborativo: % de alta del producto (columna) si ya posee el de la fila",
                    co[np.ix_(top_p, top_p)] * 100, [NOMBRES[j] for j in top_p], [NOMBRES[j] for j in top_p], fmt="pct"),
        ],
        "tables": [table("Comparación de modelos (MAP@7 con IC bootstrap por cliente)", res),
                   table(f"Diferencias pareadas contra {mejor}", pareadas),
                   table("Primera recomendación por producto", por_prod.reset_index().rename(
                       columns={"producto": "Producto", "recomendado": "Veces recomendado", "precision": "Precisión"})
                       .round(3))],
        "card": model_card(
            "Popularidad, filtrado colaborativo ítem-ítem (probabilidad condicional suavizada), LightGBM multiclase "
            "sobre altas y LightGBM LambdaRank sobre pares cliente-producto; nunca se recomienda algo que ya tiene.",
            ["Validación fuera de tiempo: entrenamiento con altas de mar–abr 2016, prueba en may 2016.",
             "MAP@7 con IC bootstrap por cliente y diferencias pareadas; aciertos por posición (hit@k)."],
            ["La alta del próximo mes depende de la tenencia actual y del perfil, no de la campaña (sin datos de oferta).",
             "Sexo no se usa como variable."],
            ["Datos de un banco español de 2015–2016; el ranking óptimo cambia con campañas y estacionalidad.",
             "Sin datos de respuesta a ofertas: optimiza 'qué contratará' y no 'qué cambia con la oferta' (uplift, #3)."]),
    }
    muestra = X_te[con_alta][["edad", "antiguedad", "renta"] + [f"{p}_t1" for p in P[:12]]].copy()
    muestra.columns = ["edad", "antiguedad", "renta"] + NOMBRES[:12]
    muestra["alta"] = [NOMBRES[np.flatnonzero(r)[0]] for r in Ya]
    export("p09-recomendador", result, sample=muestra.sample(4000, random_state=SEED).round(1))
    return result


if __name__ == "__main__":
    run()