"""
Toetst de bewering:

    D = (azen + tienen) - (vrouwen + boeren)
    D >  3  ->  hand is 1 punt meer waard
    D < -3  ->  hand is 1 punt minder waard

Methode: regressie van het dubbeldummy-aantal slagen op HCP + vormcontroles,
en kijken of D (of de aanpassing zelf) nog iets toevoegt. De coefficient van
D wordt teruggerekend naar HCP-equivalenten via de marginale waarde van 1 HCP.

Gebruik: python3 analyse.py <csv> [nt|suit]
"""

import sys

import numpy as np
import pandas as pd
import statsmodels.api as sm

pd.set_option("display.width", 200)


# ---------------------------------------------------------------- hulpfuncties
def ols(df, y, xs):
    X = sm.add_constant(df[xs].astype(float))
    return sm.OLS(df[y].astype(float), X).fit()


def hcp_marginal(res, df):
    """Marginale waarde van 1 HCP in slagen, op het gemiddelde van de steekproef."""
    b = res.params
    m = b["hcp"]
    if "hcp2" in b:
        m += 2 * b["hcp2"] * df["hcp"].mean()
    return m


def build(df, mode):
    d = pd.DataFrame(index=df.index)

    for seat in ("N", "S", "E", "W"):
        df[f"{seat}_D"] = (df[f"{seat}_a"] + df[f"{seat}_t"]) - (
            df[f"{seat}_q"] + df[f"{seat}_j"]
        )
        df[f"{seat}_adj"] = np.where(
            df[f"{seat}_D"] > 3, 1, np.where(df[f"{seat}_D"] < -3, -1, 0)
        )

    d["hcp"] = df["N_hcp"] + df["S_hcp"]
    d["hcp2"] = d["hcp"] ** 2
    for k in ("a", "k", "q", "j", "t"):
        d[k] = df[f"N_{k}"] + df[f"S_{k}"]
    d["D"] = df["N_D"] + df["S_D"]
    d["adj"] = df["N_adj"] + df["S_adj"]
    d["D_N"], d["D_S"] = df["N_D"], df["S_D"]
    d["AQJ"] = d["a"] - d["q"] - d["j"]  # helft 1 van de regel, zonder tienen
    d["tens"] = d["t"]  # helft 2 van de regel

    # vormcontroles
    suits = ("sp", "he", "ru", "kl")
    fits = np.column_stack([df[f"N_{s}"] + df[f"S_{s}"] for s in suits])
    d["fit"] = fits.max(axis=1)
    lengths = np.column_stack(
        [df[f"{seat}_{s}"] for seat in ("N", "S") for s in suits]
    )
    d["lang"] = np.clip(lengths - 4, 0, None).sum(axis=1)
    d["kort"] = np.clip(3 - lengths, 0, None).sum(axis=1)
    d["fit2"] = d["fit"] ** 2

    if mode == "nt":
        d["y"] = df[["tr_NN", "tr_SN"]].max(axis=1)
    else:
        cols = [f"tr_{s}{st}" for s in ("N", "S") for st in ("C", "D", "H", "S")]
        d["y"] = df[cols].max(axis=1)
    return d


CONTROLS = ["hcp", "hcp2", "fit", "fit2", "lang", "kort"]


def report(d, label):
    n = len(d)
    print(f"\n{'=' * 78}\n{label}   (n = {n:,})\n{'=' * 78}")

    base = ols(d, "y", CONTROLS)
    mh = hcp_marginal(base, d)
    print(f"Basismodel  R2 = {base.rsquared:.4f}   RMSE = {np.sqrt(base.mse_resid):.4f} slagen")
    print(f"1 HCP is aan de marge {mh:.4f} slagen waard  (dus 1 slag ~ {1/mh:.2f} HCP)")

    # ---- D continu
    m = ols(d, "y", CONTROLS + ["D"])
    b, se = m.params["D"], m.bse["D"]
    print(f"\n[1] y ~ controles + D")
    print(f"    coef D = {b:+.5f} slagen  (se {se:.5f}, t = {b/se:+.1f})")
    print(f"    -> {b/mh:+.4f} HCP per eenheid D   ({4*b/mh:+.2f} HCP bij D = 4)")
    print(f"    R2 {base.rsquared:.4f} -> {m.rsquared:.4f}")

    # ---- de regel zelf
    m2 = ols(d, "y", CONTROLS + ["adj"])
    b2, se2 = m2.params["adj"], m2.bse["adj"]
    print(f"\n[2] y ~ controles + aanpassing volgens de regel (-1/0/+1 per hand)")
    print(f"    coef = {b2:+.5f} slagen (se {se2:.5f}, t = {b2/se2:+.1f})")
    print(f"    -> de aanpassing is in werkelijkheid {b2/mh:+.3f} HCP waard")
    print(f"    R2 {base.rsquared:.4f} -> {m2.rsquared:.4f}")

    # ---- de twee helften apart
    m3 = ols(d, "y", CONTROLS + ["AQJ", "tens"])
    print(f"\n[3] de regel opgesplitst")
    for name, lab in (("AQJ", "azen min (vrouwen+boeren)"), ("tens", "tienen")):
        bb, ss = m3.params[name], m3.bse[name]
        print(f"    {lab:28s} {bb:+.5f} slagen -> {bb/mh:+.3f} HCP per kaart (t {bb/ss:+.1f})")

    # ---- empirische kaartwaarden
    m4 = ols(d, "y", ["a", "k", "q", "j", "t", "fit", "fit2", "lang", "kort"])
    cA, cK, cQ, cJ, cT = (m4.params[k] for k in ("a", "k", "q", "j", "t"))
    scale = 40.0 / (4 * (cA + cK + cQ + cJ))
    print(f"\n[4] empirische kaartwaarden, geschaald zodat A+H+V+B samen 40 is")
    print("    kaart :   aas    heer   vrouw   boer    tien")
    print(f"    waarde: {cA*scale:6.2f} {cK*scale:6.2f} {cQ*scale:6.2f} {cJ*scale:6.2f} {cT*scale:6.2f}")
    print("    HCP   :   4.00   3.00   2.00   1.00   0.00")

    # ---- welke drempel hoort bij precies 1 punt?
    print(f"\n[4b] welke drempel levert precies 1 punt op?")
    print("    drempel   waarde van de aanpassing (HCP)   extra verklaarde variantie")
    for thr in (1, 2, 3, 4, 5):
        x = np.where(d["D_N"] > thr, 1, np.where(d["D_N"] < -thr, -1, 0)) + np.where(
            d["D_S"] > thr, 1, np.where(d["D_S"] < -thr, -1, 0)
        )
        dd = d.assign(x=x)
        mt = ols(dd, "y", CONTROLS + ["x"])
        print(f"      > {thr}          {mt.params['x']/mh:+.3f}"
              f"                        {mt.rsquared - base.rsquared:.4f}")
    print(f"      lineair in D                                    "
          f"{m.rsquared - base.rsquared:.4f}")

    # ---- restanten per D-waarde
    d = d.copy()
    d["resid_hcp"] = base.resid / mh
    print(f"\n[5] gemiddelde onderwaardering (in HCP) per waarde van D (paar)")
    g = d.groupby("D")["resid_hcp"].agg(["count", "mean"])
    g = g[g["count"] >= max(30, n // 500)]
    for dv, row in g.iterrows():
        print(f"    D = {dv:+3.0f}  n = {int(row['count']):6d}   {row['mean']:+.3f} HCP")


def main():
    path = sys.argv[1]
    mode = sys.argv[2] if len(sys.argv) > 2 else "nt"
    raw = pd.read_csv(path)
    d = build(raw, mode)

    # frequentie van de regel per hand
    allD = pd.concat([(raw[f"{s}_a"] + raw[f"{s}_t"]) - (raw[f"{s}_q"] + raw[f"{s}_j"])
                      for s in ("N", "E", "S", "W")])
    print(f"Verdeling van D per losse hand (n = {len(allD):,})")
    print(f"  gemiddelde {allD.mean():+.3f}, sd {allD.std():.3f}")
    print(f"  D >  3 (+1 punt): {(allD > 3).mean()*100:.2f} % van de handen")
    print(f"  D < -3 (-1 punt): {(allD < -3).mean()*100:.2f} % van de handen")
    print(f"  regel doet iets : {((allD > 3) | (allD < -3)).mean()*100:.2f} % van de handen")

    report(d, f"ALLE SPELLEN  ({'SA' if mode == 'nt' else 'beste kleurcontract'})")
    sub = d[(d["hcp"] >= 20) & (d["hcp"] <= 28)]
    report(sub, f"BIEDBARE ZONE 20-28 HCP  ({'SA' if mode == 'nt' else 'beste kleurcontract'})")


if __name__ == "__main__":
    main()