"""
Tutorial 97 - Bradley-Terry ratings from scratch.

A win-loss record treats every opponent as the same. Bradley-Terry does not: it gives
every team a strength s, says team i beats team j with probability s_i / (s_i + s_j),
and finds the strengths that make the season's actual results most likely. Beat good
teams and your strength rises more than your record does.

The fit is Hunter's MM algorithm (2004): a one-line fixed-point update that is
guaranteed to climb the likelihood every iteration. A home-field factor theta goes in
the same way. Two things the textbook version leaves out are handled here: an unbeaten
or winless team has no finite maximum-likelihood strength, so a one-game prior keeps
every estimate finite; and the model is tested out of sample, where it has to beat the
plain win-loss record to be worth the trouble.

Written by hand: the MM update with and without home field, the prior, log5 from a
smoothed win percentage for comparison, log loss and Brier score, and a paired
bootstrap of the difference. No scipy, no statsmodels, no choix.

Run:

    python downloads/97_bradley_terry_ratings_from_scratch.py

Data: nflverse games table (github.com/nflverse/nfldata), June 2026 snapshot, bundled as
nfl_games_lines.csv.
"""
import math
import os

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

import sdt_common as sdt

sdt.init("bradley-terry-ratings-from-scratch")
HERE = os.path.dirname(os.path.abspath(__file__))
CSV = os.path.join(HERE, "nfl_games_lines.csv")
FRANCHISE = {"STL": "LA", "SD": "LAC", "OAK": "LV"}   # a move is not a new team
SEED = 20260923


# --- the model, written out in full ---------------------------------------------------
def fit_bt(games, home=True, prior=1.0, iters=5000, tol=1e-12):
    """Bradley-Terry by Hunter's MM algorithm.

    P(home team h beats away team a) = theta * s_h / (theta * s_h + s_a).

    The update for each team is  s_i <- W_i / sum over i's games of (weight / denominator),
    where W_i is i's wins, the weight is theta when i was at home and 1 when away, and the
    denominator is theta * s_home + s_away for that game. Theta updates the same way from
    the home wins. A tie counts half a win to each side.

    prior: every team also gets `prior` virtual wins and `prior` virtual losses against a
    reference team of strength 1, so an unbeaten or winless team still has a finite
    maximum. prior=0 is the textbook model.
    """
    teams = sorted(set(games.home_team) | set(games.away_team))
    ix = {t: i for i, t in enumerate(teams)}
    h = games.home_team.map(ix).to_numpy()
    a = games.away_team.map(ix).to_numpy()
    diff = games.home_score.to_numpy() - games.away_score.to_numpy()
    hw = np.where(diff > 0, 1.0, np.where(diff < 0, 0.0, 0.5))
    wins = np.zeros(len(teams))
    np.add.at(wins, h, hw)
    np.add.at(wins, a, 1 - hw)
    wins += prior
    s, theta, home_wins = np.ones(len(teams)), 1.0, hw.sum()
    for it in range(1, iters + 1):
        t = theta if home else 1.0
        den = t * s[h] + s[a]
        d = np.zeros(len(teams))
        np.add.at(d, h, t / den)
        np.add.at(d, a, 1 / den)
        d += 2 * prior / (s + 1.0)
        new = wins / d
        new /= np.exp(np.mean(np.log(new)))           # strengths are only defined up to scale
        new_theta = home_wins / np.sum(new[h] / (theta * new[h] + new[a])) if home else 1.0
        done = np.max(np.abs(np.log(new / s))) < tol and abs(math.log(new_theta / theta)) < tol
        s, theta = new, new_theta
        if done:
            break
    return pd.Series(s, index=teams), theta, it


def p_home(s, theta, home_team, away_team):
    return theta * s[home_team] / (theta * s[home_team] + s[away_team])


def log_loss(p, y):
    return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))


def brier(p, y):
    return float(np.mean((p - y) ** 2))


def team_records(games):
    diff = games.home_score - games.away_score
    home = pd.DataFrame({"team": games.home_team, "w": (diff > 0) + 0.5 * (diff == 0)})
    away = pd.DataFrame({"team": games.away_team, "w": (diff < 0) + 0.5 * (diff == 0)})
    both = pd.concat([home, away])
    return both.groupby("team").w.agg(wins="sum", games="count")


# --- 1. the data ------------------------------------------------------------------------
games = pd.read_csv(CSV)
games["home_team"] = games.home_team.replace(FRANCHISE)
games["away_team"] = games.away_team.replace(FRANCHISE)
reg = games[(games.game_type == "REG") & games.home_score.notna()].copy()
g25 = reg[reg.season == 2025]

with sdt.snippet("data"):
    print(f"regular-season games {reg.season.min()}-{reg.season.max()}: {len(reg):,}")
    print(f"2025: {len(g25)} games, {g25.home_team.nunique()} teams, "
          f"{int((g25.home_score == g25.away_score).sum())} ties")
    print(f"home teams won {(g25.home_score > g25.away_score).mean():.4f} of 2025 games")

# --- 2. the fit -------------------------------------------------------------------------
s25, th25, it25 = fit_bt(g25)
with sdt.snippet("fit"):
    print(f"converged in {it25} MM iterations")
    print(f"home-field factor theta = {th25:.4f}")
    print(f"two equal teams: the home side wins {th25 / (1 + th25):.4f}")
    rec = team_records(g25)
    tbl = pd.DataFrame({"strength": s25, "wins": rec.wins, "games": rec.games})
    tbl["win_pct"] = tbl.wins / tbl.games
    tbl = tbl.sort_values("strength", ascending=False)
    tbl["bt_rank"] = range(1, len(tbl) + 1)
    tbl["pct_rank"] = tbl.win_pct.rank(ascending=False, method="min").astype(int)
    print(tbl.head(8).round(4).to_string())
    print("...")
    print(tbl.tail(3).round(4).to_string())

# --- 3. same record, different strength --------------------------------------------------
with sdt.snippet("schedule"):
    for team in ("SEA", "DEN", "NE"):
        mine = g25[(g25.home_team == team) | (g25.away_team == team)]
        opps = np.where(mine.home_team == team, mine.away_team, mine.home_team)
        print(f"{team}: {tbl.loc[team, 'wins']:.0f}-{tbl.loc[team, 'games'] - tbl.loc[team, 'wins']:.0f}, "
              f"strength {s25[team]:.3f} (rank {tbl.loc[team, 'bt_rank']}), "
              f"opponents' mean strength {np.exp(np.mean(np.log(s25[opps]))):.3f}")
    se, ne = s25["SEA"], s25["NE"]
    print(f"SEA over NE on a neutral field: {se / (se + ne):.4f}")

# --- 4. why the prior exists --------------------------------------------------------------
early = g25[g25.week <= 4]
rec4 = team_records(early)
unbeaten = sorted(rec4[rec4.wins == rec4.games].index)
winless = sorted(rec4[rec4.wins == 0].index)
def log_lik(s, theta, games):
    """Log-likelihood of the actual results under strengths s (ties count half each way)."""
    p = np.array([p_home(s, theta, h, a) for h, a in zip(games.home_team, games.away_team)])
    diff = (games.home_score - games.away_score).to_numpy()
    w = np.where(diff > 0, 1.0, np.where(diff < 0, 0.0, 0.5))
    return float(np.sum(w * np.log(p) + (1 - w) * np.log(1 - p)))


def one_mm_step_no_prior(games):
    """The first textbook update from equal strengths, theta held at 1."""
    teams = sorted(set(games.home_team) | set(games.away_team))
    rec = team_records(games)
    games_played = rec.games.reindex(teams).to_numpy(float)
    return pd.Series(rec.wins.reindex(teams).to_numpy(float) / (games_played / 2.0), index=teams)


with sdt.snippet("prior"):
    print(f"after 4 weeks: {len(unbeaten)} unbeaten {unbeaten}, {len(winless)} winless {winless}")
    step = one_mm_step_no_prior(early)
    print("prior 0, first update from equal strengths:",
          ", ".join(f"{t} {step[t]:.3f}" for t in winless + unbeaten))
    s1, th1, it1 = fit_bt(early, prior=1.0)
    print(f"prior 1: converged in {it1} iterations, strengths {s1.min():.3f} to {s1.max():.3f}")
    for k in (1, 10, 100, 1000):
        s_k = s1.copy()
        s_k["BUF"] = s1["BUF"] * k
        print(f"BUF's strength x{k:>4}: log-likelihood of the 64 real results {log_lik(s_k, th1, early):.4f}")

# --- 5. out of sample ---------------------------------------------------------------------
rows, per_game = [], []
for yr in range(2010, 2026):
    season = reg[reg.season == yr]
    train, test = season[season.week <= 12], season[season.week > 12]
    test = test[test.home_score != test.away_score]
    y = (test.home_score > test.away_score).to_numpy(float)
    s_h, th, _ = fit_bt(train)
    s_n, _, _ = fit_bt(train, home=False)
    p_bt = np.array([p_home(s_h, th, h, a) for h, a in zip(test.home_team, test.away_team)])
    p_bn = np.array([p_home(s_n, 1.0, h, a) for h, a in zip(test.home_team, test.away_team)])
    rec = team_records(train)
    wp = (rec.wins + 1) / (rec.games + 2)              # the same one-game prior, on the record
    pa, pb = wp[test.home_team].to_numpy(), wp[test.away_team].to_numpy()
    p_l5 = (pa - pa * pb) / (pa + pb - 2 * pa * pb)     # log5
    home_rate = (train.home_score > train.away_score).mean()
    rows.append({"season": yr, "games": len(test),
                 "coin": log_loss(np.full(len(y), 0.5), y),
                 "home_rate": log_loss(np.full(len(y), home_rate), y),
                 "log5": log_loss(p_l5, y), "bt_no_home": log_loss(p_bn, y),
                 "bt": log_loss(p_bt, y), "brier_log5": brier(p_l5, y), "brier_bt": brier(p_bt, y)})
    per_game.append(pd.DataFrame({"y": y, "bt": p_bt, "log5": p_l5}))
oos = pd.DataFrame(rows)
pg = pd.concat(per_game, ignore_index=True)
ll_bt = -(pg.y * np.log(pg.bt) + (1 - pg.y) * np.log(1 - pg.bt))
ll_l5 = -(pg.y * np.log(pg.log5) + (1 - pg.y) * np.log(1 - pg.log5))
d = (ll_l5 - ll_bt).to_numpy()
rng = np.random.default_rng(SEED)
boot = np.array([d[rng.integers(0, len(d), len(d))].mean() for _ in range(10000)])
lo, hi = np.percentile(boot, [2.5, 97.5])

with sdt.snippet("oos"):
    print(oos.round(4).to_string(index=False))
    print()
    m = oos.mean(numeric_only=True)
    print(f"mean log loss  coin {m.coin:.4f}  home rate {m.home_rate:.4f}  log5 {m.log5:.4f}  "
          f"BT no home {m.bt_no_home:.4f}  BT {m.bt:.4f}")
    print(f"mean Brier     log5 {m.brier_log5:.4f}  BT {m.brier_bt:.4f}")
    print(f"BT beats log5 in {(oos.bt < oos.log5).sum()} of {len(oos)} seasons")
    print(f"per game over {len(d):,} games: log5 minus BT = {d.mean():+.5f}, "
          f"95% bootstrap interval {lo:+.5f} to {hi:+.5f}")

# --- 6. the chart -------------------------------------------------------------------------
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.5, 5.2))
ax1.scatter(tbl.win_pct, tbl.strength, s=26, color=sdt.sport_color("football"))
for team in ("SEA", "DEN", "NE", "JAX", "LV", "NYJ"):
    ax1.annotate(team, (tbl.loc[team, "win_pct"], tbl.loc[team, "strength"]), textcoords="offset points",
                 xytext=(5, 3), fontsize=9)
ax1.set_yscale("log")
ax1.set_xlabel("2025 win percentage")
ax1.set_ylabel("Bradley-Terry strength (log scale)")
ax1.set_title("Same record, different strength: SEA, DEN and NE all went 14-3")
x = np.arange(len(oos))
ax2.plot(x, oos.bt, marker="o", label="Bradley-Terry with home field")
ax2.plot(x, oos.log5, marker="s", label="log5 from win %")
ax2.plot(x, oos.home_rate, ls="--", color="grey", label="home-win rate only")
ax2.axhline(math.log(2), color="black", lw=0.8, ls=":", label="coin flip")
ax2.set_xticks(x[::3])
ax2.set_xticklabels(oos.season[::3])
ax2.set_ylabel("log loss on weeks 13+ (lower is better)")
ax2.set_title("Out of sample, the two finish level")
ax2.legend(fontsize=8)
sdt.save_fig(fig, "bt_strength_and_oos",
             source="nflverse games table (June 2026 snapshot), regular seasons 2010-2025")