"""
Tutorial 91 - Is xG calibrated? Checking a published model against 1,494 real shots.

Every xG number is a promise: "shots like this one score 12% of the time." This
tutorial takes StatsBomb's own xG for every shot of the 2022 World Cup and audits
the promise three ways - calibration-in-the-large (does total xG match total
goals?), a reliability table (do 0.05 shots score 5% of the time?), and a Brier
score against the laziest possible baseline. Along the way the data teaches the
first rule of model auditing: know what a row IS before you average it - this
file quietly contains 41 penalty-shootout kicks that are not match goals at all.

Runs entirely offline from the bundled wc2022_shots.csv next to this script
(no API, no key). Data provided by StatsBomb (https://github.com/statsbomb/
open-data), CC BY 4.0 - attribution required, gladly given.

Run:  python downloads/91_is_xg_calibrated.py
"""
import os

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import sdt_common as sdt

sdt.init("is-xg-calibrated")
HERE = os.path.dirname(os.path.abspath(__file__))
CSV = os.path.join(HERE, "wc2022_shots.csv")

shots = pd.read_csv(CSV)

# --- know what a row is --------------------------------------------------------------
# 1,494 rows and 195 "goals" - but 41 rows are period 5: the penalty shoot-outs.
# A shoot-out kick is a shot event with an xG value, yet its goal does not count
# in any match score. Own goals, meanwhile, are nobody's shot, so they are not
# in a shot table at all. Averaging without knowing this bakes both mistakes in.
with sdt.snippet("inventory"):
    print(f"{len(shots)} shots, {shots.is_goal.sum()} rows flagged as goals")
    print("\nby period (5 = penalty shoot-outs):")
    tbl = shots.groupby("period").agg(rows=("xg", "size"), goals=("is_goal", "sum"),
                                      penalties=("is_penalty", "sum"))
    sdt.show_df(tbl.reset_index())
    so = shots[shots.period == 5]
    print(f"\nshoot-out rows: {len(so)}, 'goals' among them: {so.is_goal.sum()} - "
          f"none of which appear in a match score")
    print("own goals: not shots, so never in this table at all")

match = shots[shots.period < 5].copy()          # real match play only
open_sp = match[match.is_penalty == 0]           # non-penalty match shots
pens = match[match.is_penalty == 1]              # in-game penalties

# --- calibration in the large --------------------------------------------------------
# If the model is honest, the sum of xG over many shots should land near the
# actual goal count. The standard error of that promise is sqrt(sum p*(1-p)) -
# each shot is a Bernoulli trial with its own p.
with sdt.snippet("large"):
    for name, df in [("non-penalty match shots", open_sp),
                     ("in-game penalties", pens),
                     ("all match shots", match)]:
        exp, act = df.xg.sum(), df.is_goal.sum()
        se = np.sqrt((df.xg * (1 - df.xg)).sum())
        z = (act - exp) / se
        print(f"{name:24}  n={len(df):5}  xG {exp:7.2f}  goals {act:4}  "
              f"diff {act - exp:+6.2f}  SE {se:5.2f}  z {z:+.2f}")
    naive = shots.xg.sum(), shots.is_goal.sum()
    print(f"\nnaive (shoot-outs left in): xG {naive[0]:.2f} vs 'goals' {naive[1]} - "
          f"a fake surplus built from rows that are not match goals")

# --- the reliability table -----------------------------------------------------------
# Bin non-penalty shots by their promised probability; compare the promise
# (mean xG in the bin) to what happened (goal rate), with a Wilson 95% interval
# on the observed rate so small bins cannot pretend to be precise.
def wilson(k, n, z=1.96):
    if n == 0:
        return (np.nan, np.nan)
    p = k / n
    denom = 1 + z * z / n
    centre = (p + z * z / (2 * n)) / denom
    half = z * np.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / denom
    return (centre - half, centre + half)

EDGES = [0, 0.02, 0.05, 0.10, 0.20, 0.40, 1.00]
LABELS = ["0-.02", ".02-.05", ".05-.10", ".10-.20", ".20-.40", ".40-1"]
open_sp = open_sp.assign(bin=pd.cut(open_sp.xg, EDGES, labels=LABELS, right=False))

rows = []
for lab, grp in open_sp.groupby("bin", observed=True):
    n, k = len(grp), int(grp.is_goal.sum())
    lo, hi = wilson(k, n)
    rows.append({"bin": lab, "shots": n, "mean_xg": grp.xg.mean(),
                 "goal_rate": k / n, "ci_lo": lo, "ci_hi": hi,
                 "promise_inside": lo <= grp.xg.mean() <= hi})
rel = pd.DataFrame(rows)

with sdt.snippet("bins"):
    out = rel.copy()
    for c in ["mean_xg", "goal_rate", "ci_lo", "ci_hi"]:
        out[c] = out[c].map("{:.3f}".format)
    sdt.show_df(out)
    inside = int(rel.promise_inside.sum())
    print(f"\nbins where the promised rate sits inside the observed 95% CI: "
          f"{inside} of {len(rel)}")

# --- Brier score vs the laziest baseline ---------------------------------------------
# The Brier score is the mean squared error of a probability. Any candidate
# model must beat "predict the base rate for every shot" or it has learned
# nothing about individual shots.
with sdt.snippet("brier"):
    p = open_sp.xg.to_numpy()
    y = open_sp.is_goal.to_numpy()
    base = y.mean()
    bs_model = np.mean((p - y) ** 2)
    bs_base = np.mean((base - y) ** 2)
    skill = 1 - bs_model / bs_base
    print(f"base rate: {base:.4f}  ({y.sum()} of {len(y)} non-penalty shots scored)")
    print(f"Brier, xG model:      {bs_model:.4f}")
    print(f"Brier, base rate:     {bs_base:.4f}")
    print(f"skill vs base rate:   {skill:+.1%}")
    print(f"mean xG when scored:  {p[y == 1].mean():.3f}")
    print(f"mean xG when missed:  {p[y == 0].mean():.3f}")

# --- the picture ---------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(7.6, 5.4))
ax.plot([0, 0.55], [0, 0.55], ls="--", lw=1, color="#8a8a8a",
        label="perfect calibration (y = x)")
ax.errorbar(rel.mean_xg, rel.goal_rate,
            yerr=[rel.goal_rate - rel.ci_lo, rel.ci_hi - rel.goal_rate],
            fmt="o", ms=7, lw=1.6, capsize=4, color=sdt.sport_color("soccer"),
            label="observed goal rate (Wilson 95% CI)")
for _, r in rel.iterrows():
    ax.annotate(f"n={r.shots}", (r.mean_xg, r.goal_rate),
                textcoords="offset points", xytext=(8, -11), fontsize=8, color="#666")
ax.set_xlabel("promised probability (mean xG in bin)")
ax.set_ylabel("what happened (goal rate in bin)")
ax.set_title("StatsBomb xG vs reality - 1,430 non-penalty shots, World Cup 2022")
ax.legend(loc="upper left", frameon=False)
ax.set_xlim(0, 0.55)
ax.set_ylim(0, 0.62)
fig.tight_layout()
sdt.save_fig(fig, "reliability", source="StatsBomb open data (CC BY 4.0), WC2022, 1,430 non-penalty shots")