#!/usr/bin/env python3
"""
Construye las tablas derivadas del post "Los que no llegaron".

Fuentes (públicas, del repositorio jmtoral/admisiones_unam, que a su vez
raspa las páginas públicas de la DGAE):

  - data/consolidated/resultados_control_2026.csv  (una fila por aspirante
    citado al examen de control 2026)
  - data/consolidated/metadata_control_2026.csv    (encabezados de la DGAE
    por oferta, examen de control)
  - data/consolidated/metadata_carreras.csv        (encabezados por oferta y
    año, 2021-2026, examen en línea)

Salida: CSV agregados a nivel oferta (carrera-campus). Ningún archivo de
salida contiene información a nivel de aspirante.

Uso:  python3 build_data.py
Sólo biblioteca estándar. Idempotente.
"""

import csv
import io
import json
import math
import os
import statistics as st
import urllib.request

BASE = ("https://raw.githubusercontent.com/jmtoral/admisiones_unam/"
        "main/data/consolidated/")
FUENTES = {
    "resultados_control_2026.csv": BASE + "resultados_control_2026.csv",
    "metadata_control_2026.csv": BASE + "metadata_control_2026.csv",
    "metadata_carreras.csv": BASE + "metadata_carreras.csv",
    # Salida del análisis de asistencia de jmtoral, para contrastar su
    # denominador estimado contra el que publica la propia DGAE.
    "presentaron_control.csv": ("https://raw.githubusercontent.com/jmtoral/"
                                "admisiones_unam/main/analysis/output/"
                                "presentaron_control.csv"),
}
AQUI = os.path.dirname(os.path.abspath(__file__))
CACHE = os.path.join(AQUI, ".cache")
SALIDA = os.path.join(AQUI, "data")

# Umbrales del conjunto analizable. Declarados aquí, no escondidos en el texto.
MIN_CITADOS = 50   # ofertas chicas dan tasas inestables
MIN_ANIOS = 3      # años de historia 2021-2025 para fijar una línea base


def descarga(nombre, url):
    """Baja una vez y cachea. La fuente no cambia retroactivamente."""
    destino = os.path.join(CACHE, nombre)
    if not os.path.exists(destino):
        os.makedirs(CACHE, exist_ok=True)
        with urllib.request.urlopen(url, timeout=120) as r:
            datos = r.read()
        with open(destino, "wb") as f:
            f.write(datos)
    with open(destino, encoding="utf-8") as f:
        return list(csv.DictReader(f))


def num(x):
    x = (x or "").replace(",", "").strip()
    return int(x) if x.isdigit() else None


def wls(y, X, w):
    """Mínimos cuadrados ponderados. Pesos analíticos (se normalizan a media 1),
    no de frecuencia: de lo contrario el error estándar se infla con el tamaño
    de las ofertas y todo sale significativo."""
    k, N = len(X[0]), len(y)
    s = sum(w)
    w = [x * N / s for x in w]
    A = [[sum(w[i] * X[i][a] * X[i][b] for i in range(N)) for b in range(k)]
         for a in range(k)]
    B = [sum(w[i] * X[i][a] * y[i] for i in range(N)) for a in range(k)]
    M = [A[i][:] + [1.0 if j == i else 0.0 for j in range(k)] + [B[i]]
         for i in range(k)]
    for c in range(k):                      # Gauss-Jordan: resuelve e invierte
        p = max(range(c, k), key=lambda r: abs(M[r][c]))
        M[c], M[p] = M[p], M[c]
        d = M[c][c]
        for j in range(2 * k + 1):
            M[c][j] /= d
        for r in range(k):
            if r != c and M[r][c]:
                f = M[r][c]
                for j in range(2 * k + 1):
                    M[r][j] -= f * M[c][j]
    beta = [M[i][2 * k] for i in range(k)]
    inv = [[M[i][k + j] for j in range(k)] for i in range(k)]
    res = [y[i] - sum(b * x for b, x in zip(beta, X[i])) for i in range(N)]
    s2 = sum(w[i] * res[i] ** 2 for i in range(N)) / (N - k)
    se = [math.sqrt(max(s2 * inv[i][i], 0)) for i in range(k)]
    ybar = sum(w[i] * y[i] for i in range(N)) / sum(w)
    sst = sum(w[i] * (y[i] - ybar) ** 2 for i in range(N))
    r2 = 1 - sum(w[i] * res[i] ** 2 for i in range(N)) / sst
    return beta, se, r2


def escribe(nombre, filas, campos):
    os.makedirs(SALIDA, exist_ok=True)
    with open(os.path.join(SALIDA, nombre), "w", newline="",
              encoding="utf-8") as f:
        w = csv.DictWriter(f, fieldnames=campos)
        w.writeheader()
        for r in filas:
            w.writerow({k: r[k] for k in campos})
    print(f"  escrito data/{nombre} ({len(filas)} filas)")


def main():
    res = descarga("resultados_control_2026.csv",
                   FUENTES["resultados_control_2026.csv"])
    mctl = descarga("metadata_control_2026.csv",
                    FUENTES["metadata_control_2026.csv"])
    mcar = descarga("metadata_carreras.csv", FUENTES["metadata_carreras.csv"])

    # ---- 1. Validación: ¿'acreditado=N' sin aciertos es "no se presentó"? ----
    # La DGAE publica en el encabezado de cada tabla cuántos citados y cuántos
    # presentaron. Si la resta coincide con el conteo de filas N-sin-aciertos,
    # la interpretación queda anclada en la propia fuente y no en una suposición.
    # Ojo con la llave: `codigo` NO es único. Nueve códigos aparecen a la vez
    # en abierta y en SUAyED (Derecho 30500076, Contaduría 30400066, ...).
    # Unir sólo por `codigo` funde dos poblaciones distintas. La llave es
    # (modalidad, codigo).
    n_sin_por_llave = {}
    for r in res:
        if r["acreditado"] == "N" and not r["aciertos"]:
            k = (r["modalidad"], r["codigo"])
            n_sin_por_llave[k] = n_sin_por_llave.get(k, 0) + 1

    val = []
    for m in mctl:
        cod = m["codigo"]
        n_sin = n_sin_por_llave.get((m["modalidad"], cod), 0)
        cit, pres = num(m["aspirantes"]), num(m["presentaron_examen"])
        val.append(dict(codigo=cod, modalidad=m["modalidad"],
                        carrera=m["carrera"], campus=m["campus"],
                        citados=cit, presentaron=pres,
                        ausentes_encabezado=cit - pres,
                        filas_N_sin_aciertos=n_sin,
                        coincide=int(n_sin == cit - pres)))
    escribe("validacion_ausentes.csv", val,
            ["codigo", "modalidad", "carrera", "campus", "citados",
             "presentaron", "ausentes_encabezado", "filas_N_sin_aciertos",
             "coincide"])
    ok = sum(v["coincide"] for v in val)
    print(f"  validación: {ok}/{len(val)} ofertas cuadran al dígito")

    # ---- 2. Resumen por modalidad ----
    resumen = []
    for mod in ("escolarizado", "abierta", "suayed"):
        sub = [m for m in mctl if m["modalidad"] == mod]
        cit = sum(num(m["aspirantes"]) for m in sub)
        pres = sum(num(m["presentaron_examen"]) for m in sub)
        resumen.append(dict(modalidad=mod, ofertas=len(sub), citados=cit,
                            presentaron=pres, ausentes=cit - pres,
                            tasa_ausentismo=round((cit - pres) / cit, 4)))
    escribe("resumen_modalidad.csv", resumen,
            ["modalidad", "ofertas", "citados", "presentaron", "ausentes",
             "tasa_ausentismo"])

    # ---- 3. Tabla por oferta (escolarizado) ----
    # 'salto' = cuánto subió el mínimo del examen en línea 2026 respecto de la
    # mediana 2021-2025 de esa misma oferta. Es el proxy de qué tan anómalo fue
    # el examen en línea ahí. La mediana, no el promedio: 2021 fue atípico por
    # la pandemia en varias carreras.
    en_linea, hist = {}, {}
    for m in mcar:
        if m["modalidad"] != "escolarizado":
            continue
        if m["year"] == "2026":
            en_linea[m["codigo"]] = m
        else:
            v = num(m["aciertos_minimos"])
            if v is not None:
                hist.setdefault(m["codigo"], []).append(v)
    # (en escolarizado `codigo` sí es único: los nueve duplicados son de
    #  abierta/SUAyED, que este bloque ya excluyó)

    ofertas = []
    for m in mctl:
        if m["modalidad"] != "escolarizado":
            continue
        cod = m["codigo"]
        o, h = en_linea.get(cod), hist.get(cod, [])
        cit, pres = num(m["aspirantes"]), num(m["presentaron_examen"])
        min_ctl = num(m["aciertos_minimos"])
        min_lin = num(o["aciertos_minimos"]) if o else None
        base = st.median(h) if len(h) >= MIN_ANIOS else None
        asp_lin = num(o["aspirantes"]) if o else None
        lugares = num(o["oferta"]) if o else None
        fila = dict(
            codigo=cod, area=m["area"], carrera=m["carrera"],
            campus=m["campus"], citados=cit, presentaron=pres,
            ausentes=cit - pres, tasa_ausentismo=round((cit - pres) / cit, 4),
            min_linea_2026=min_lin, min_control_2026=min_ctl,
            min_historico=base, anios_historia=len(h),
            salto=(min_lin - base) if (min_lin is not None and base is not None)
                  else None,
            caida=(min_lin - min_ctl) if (min_lin is not None
                                          and min_ctl is not None) else None,
            lugares=lugares, aspirantes_linea=asp_lin,
            presion=round(asp_lin / lugares, 2)
                    if (asp_lin and lugares) else None,
        )
        fila["analizable"] = int(cit >= MIN_CITADOS and len(h) >= MIN_ANIOS
                                 and fila["salto"] is not None
                                 and fila["presion"] is not None)
        ofertas.append(fila)
    ofertas.sort(key=lambda r: -r["tasa_ausentismo"])
    escribe("ofertas_escolarizado.csv", ofertas,
            ["codigo", "area", "carrera", "campus", "citados", "presentaron",
             "ausentes", "tasa_ausentismo", "min_historico", "min_linea_2026",
             "min_control_2026", "salto", "caida", "lugares",
             "aspirantes_linea", "presion", "anios_historia", "analizable"])

    an = [r for r in ofertas if r["analizable"]]
    print(f"  conjunto analizable: {len(an)} ofertas, "
          f"{sum(r['citados'] for r in an)} citados")

    # ---- 4. Quintiles de salto ----
    s = sorted(an, key=lambda r: r["salto"])
    n, q = len(s), len(s) // 5
    quint = []
    for i in range(5):
        g = s[i * q:(i + 1) * q] if i < 4 else s[4 * q:]
        cit = sum(r["citados"] for r in g)
        aus = sum(r["ausentes"] for r in g)
        quint.append(dict(quintil=i + 1, salto_min=g[0]["salto"],
                          salto_max=g[-1]["salto"], ofertas=len(g),
                          citados=cit, ausentes=aus,
                          tasa_ausentismo=round(aus / cit, 4)))
    escribe("quintiles_salto.csv", quint,
            ["quintil", "salto_min", "salto_max", "ofertas", "citados",
             "ausentes", "tasa_ausentismo"])

    # ---- 5. Denominador: estimado vs. publicado ----
    # jmtoral estima a los convocados aplicando el umbral mínimo sobre el
    # examen en línea (lo documenta como estimación). El número exacto está en
    # el encabezado de la DGAE del propio examen de control.
    tor = descarga("presentaron_control.csv", FUENTES["presentaron_control.csv"])
    oficial = {}
    for m in mctl:
        k = (m["carrera"], m["campus"], m["modalidad"])
        a, p = num(m["aspirantes"]), num(m["presentaron_examen"])
        prev = oficial.get(k, (0, 0))
        oficial[k] = (prev[0] + a, prev[1] + p)

    den = []
    for r in tor:
        k = (r["carrera"], r["campus"], r["modalidad"])
        if k not in oficial:
            continue
        conv_est = round(float(r["convocados"]))
        pres_est = round(float(r["presc"]))
        conv_ofi, pres_ofi = oficial[k]
        den.append(dict(
            carrera=r["carrera"], campus=r["campus"], modalidad=r["modalidad"],
            convocados_estimado=conv_est, citados_publicados=conv_ofi,
            diferencia=conv_est - conv_ofi,
            presentaron_estimado=pres_est, presentaron_publicados=pres_ofi,
            asistencia_estimada=round(pres_est / conv_est, 4) if conv_est else "",
            asistencia_publicada=round(pres_ofi / conv_ofi, 4) if conv_ofi else "",
            coincide_convocados=int(conv_est == conv_ofi),
            coincide_presentaron=int(pres_est == pres_ofi)))
    den.sort(key=lambda r: r["diferencia"])
    escribe("denominador.csv", den,
            ["carrera", "campus", "modalidad", "convocados_estimado",
             "citados_publicados", "diferencia", "presentaron_estimado",
             "presentaron_publicados", "asistencia_estimada",
             "asistencia_publicada", "coincide_convocados",
             "coincide_presentaron"])
    print(f"  denominador: convocados coinciden en "
          f"{sum(r['coincide_convocados'] for r in den)}/{len(den)}; "
          f"presentaron en {sum(r['coincide_presentaron'] for r in den)}/{len(den)}")

    # ---- 6. Modelo ----
    # Controles: el mínimo histórico (qué tan competida es de suyo la oferta),
    # el tamaño, y la presión de demanda. Si el salto sobrevive a los tres,
    # no es un reflejo de que las carreras difíciles tengan más ausentismo.
    y = [r["tasa_ausentismo"] * 100 for r in an]
    w = [r["citados"] for r in an]
    nombres = ["(intercepto)", "salto del mínimo en línea",
               "mínimo histórico 2021-2025", "log(citados)",
               "presión (aspirantes/lugar)"]
    X = [[1.0, r["salto"], r["min_historico"], math.log(r["citados"]),
          r["presion"]] for r in an]
    b, se, r2 = wls(y, X, w)
    _, _, r2_sin = wls(y, [[1.0, r["min_historico"], math.log(r["citados"]),
                            r["presion"]] for r in an], w)
    modelo = [dict(variable=nm, coeficiente=round(bb, 3),
                   error_estandar=round(ss, 3), t=round(bb / ss, 2))
              for nm, bb, ss in zip(nombres, b, se)]
    escribe("modelo.csv", modelo,
            ["variable", "coeficiente", "error_estandar", "t"])

    cifras = dict(n_ofertas_modelo=len(an), r2=round(r2, 3),
                  r2_sin_salto=round(r2_sin, 3), min_citados=MIN_CITADOS,
                  min_anios=MIN_ANIOS, ofertas_validadas=ok,
                  ofertas_totales=len(val), den_ofertas=len(den),
                  den_coincide_conv=sum(r["coincide_convocados"] for r in den),
                  den_coincide_pres=sum(r["coincide_presentaron"] for r in den))
    # En CSV para que el .qmd no necesite jsonlite; en JSON por comodidad.
    escribe("cifras.csv", [cifras], list(cifras))
    with open(os.path.join(SALIDA, "cifras.json"), "w", encoding="utf-8") as f:
        json.dump(cifras, f, ensure_ascii=False, indent=2)
    print(f"  modelo: R2 = {r2:.3f} (sin el salto: {r2_sin:.3f})")


if __name__ == "__main__":
    main()