"""
Tutorial 99 - Principal component analysis from scratch: where a box score hides its shooting.

Ten per-36-minute box-score rates describe every NBA player who logged 1,000 minutes in
2025-26. Principal component analysis (PCA) asks how many independent numbers those ten
really hold, by finding the directions in which the players differ most. Two directions
carry most of it: how much of the offence a player carries, and whether he plays near the
basket or away from it. A shuffled-noise test keeps a third.

The lesson is in the component PCA ranks last. It holds 0.18% of the variance, and it is
where most of the shooting-efficiency signal lives, because points are almost a fixed
function of attempts and the only freedom left is how well the attempts convert. Keep the
top three components, the usual compression, and Stephen Curry's shooting disappears.

Written by hand: standardization, the correlation matrix, power iteration with deflation
for every eigenvalue and eigenvector, sign orientation, scores, parallel analysis with
shuffled columns, the decomposition of a shooting measure across the components, rank-k
reconstruction, next-season persistence and a bootstrap over players. numpy's own
eigensolver is used only to check the hand-written one and, once it agrees, to run the
22,000 shuffles quickly. No scikit-learn. Every number the tutorial page quotes is
asserted below.

Run:

    python downloads/99_principal_component_analysis_from_scratch.py

Data: NBA Stats API (stats.nba.com, LeagueDashPlayerStats Base and Advanced, regular-season
totals) via nba_api, retrieved September 2026, bundled as nba_player_seasons.csv.
"""
import math
import os
import sys

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.patches import Patch

import sdt_common as sdt

try:                                   # player names such as Jokić print on any console
    sys.stdout.reconfigure(encoding="utf-8")
except (AttributeError, ValueError):
    pass

sdt.init("principal-component-analysis-from-scratch")
HERE = os.path.dirname(os.path.abspath(__file__))
CSV = os.path.join(HERE, "nba_player_seasons.csv")
SEED = 20261007
SEASON = "2025-26"
MIN_MINUTES = 1000
SHUFFLES = 1000                   # parallel-analysis shuffles per season
BOOT = 2000                       # bootstrap resamples of players
STATS = ["pts", "fg2a", "fg3a", "fta", "oreb", "dreb", "ast", "tov", "stl", "blk"]
NAMES = ["points", "2-pt attempts", "3-pt attempts", "free-throw attempts", "off. rebounds",
         "def. rebounds", "assists", "turnovers", "steals", "blocks"]
COLS = [s + "36" for s in STATS]


def near(a, b, tol=5e-4):
    return abs(a - b) < tol


# --- the method, written out in full -------------------------------------------------------
def standardize(X):
    """z-scores per column: subtract the mean, divide by the standard deviation (ddof=0)."""
    mu, sd = X.mean(axis=0), X.std(axis=0)
    return (X - mu) / sd, mu, sd


def power_iteration(A, tol=1e-12, max_iter=100_000):
    """The dominant eigenvalue and eigenvector of a symmetric matrix A.

    Multiply a vector by A and rescale it, again and again: the part of the vector that
    lies along the top eigenvector grows fastest, so the vector turns towards it. Each step
    shrinks the leftover by the ratio of the second eigenvalue to the first.
    """
    v = np.random.default_rng(0).normal(size=len(A))      # any start not exactly perpendicular
    v /= np.linalg.norm(v)
    for it in range(1, max_iter + 1):
        w = A @ v
        w /= np.linalg.norm(w)
        if np.abs(w - v).max() < tol:
            return float(w @ A @ w), w, it
        v = w
    raise RuntimeError("power iteration did not converge")


def pca_by_hand(A):
    """Every eigenpair of a symmetric positive semi-definite matrix, largest first.

    Find the top pair by power iteration, subtract it (deflation: A - lambda v v^T removes
    that direction and leaves the rest of A untouched), and repeat on what is left.
    """
    A = A.copy()
    vals, vecs, iters = [], [], []
    for _ in range(len(A)):
        lam, v, it = power_iteration(A)
        vals.append(lam)
        vecs.append(v)
        iters.append(it)
        A = A - lam * np.outer(v, v)
    return np.array(vals), np.column_stack(vecs), iters


def orient(V):
    """An eigenvector's sign is arbitrary: flip each one so its largest loading is positive."""
    return V * np.sign(V[np.abs(V).argmax(axis=0), np.arange(V.shape[1])])


def shuffled_eigenvalues(Z, reps, seed):
    """Eigenvalues of the correlation matrix after shuffling every column on its own.

    Shuffling keeps each column's values and destroys every relationship between columns,
    so these are the eigenvalues that pure chance produces for data of this exact shape.
    """
    rng = np.random.default_rng(seed)
    out = np.empty((reps, Z.shape[1]))
    for i in range(reps):
        Zs = rng.permuted(Z, axis=0)          # each column permuted independently
        out[i] = np.linalg.eigvalsh(Zs.T @ Zs / len(Zs))[::-1]
    return out


def parallel_keep(vals, null):
    """Leading components whose eigenvalue beats the 95th percentile of the shuffled ones."""
    cut = np.percentile(null, 95, axis=0)
    keep = 0
    while keep < len(vals) and vals[keep] > cut[keep]:
        keep += 1
    return keep, cut


def shooting_points(frame):
    """Points per 36 minus what the same attempts score at the group's own conversion rates."""
    per2 = 2 * (frame.fgm - frame.fg3m).sum() / frame.fg2a.sum()
    per3 = 3 * frame.fg3m.sum() / frame.fg3a.sum()
    perft = frame.ftm.sum() / frame.fta.sum()
    expected = per2 * frame.fg2a36 + per3 * frame.fg3a36 + perft * frame.fta36
    return (frame.pts36 - expected).to_numpy(), (per2, per3, perft)


def fit_season(df, season, min_minutes=MIN_MINUTES):
    """Standardize one season's qualifying players and decompose their correlation matrix."""
    d = df[(df.season == season) & (df["min"] >= min_minutes)].reset_index(drop=True)
    Z, mu, sd = standardize(d[COLS].to_numpy(float))
    w, W = np.linalg.eigh(Z.T @ Z / len(Z))
    order = np.argsort(w)[::-1]
    return d, Z, w[order], orient(W[:, order])


# --- 1. ten rates per 36 minutes, standardized ------------------------------------------------
df = pd.read_csv(CSV)
df["fg2a"] = df.fga - df.fg3a
for s in STATS:
    df[s + "36"] = 36 * df[s] / df["min"]

d = df[(df.season == SEASON) & (df["min"] >= MIN_MINUTES)].reset_index(drop=True)
X = d[COLS].to_numpy(float)
Z, mu, sd = standardize(X)
R = Z.T @ Z / len(Z)                       # the correlation matrix of the ten rates
raw_share = sd ** 2 / (sd ** 2).sum()

with sdt.snippet("data"):
    print(f"bundled file: {len(df):,} player-seasons, {df.season.nunique()} seasons "
          f"({df.season.min()} to {df.season.max()})")
    print(f"{SEASON} regular season: {int((df.season == SEASON).sum())} players, "
          f"{len(d)} with {MIN_MINUTES:,}+ minutes")
    print(f"{'per 36 minutes':22s}{'mean':>8s}{'sd':>8s}   share of the raw variance")
    for j, name in enumerate(NAMES):
        print(f"{name:22s}{mu[j]:8.3f}{sd[j]:8.3f}   {raw_share[j]:6.1%}")
    print("standardized: every column has mean 0 and sd 1, so each holds 10% of the variance")

assert len(df) == 11059 and df.season.nunique() == 22
assert df.season.min() == "2004-05" and df.season.max() == "2025-26"
assert not df.duplicated(["season", "player_id"]).any()
assert int((df.season == SEASON).sum()) == 582 and len(d) == 279
assert (df.pts == 2 * df.fgm + df.fg3m + df.ftm).all()            # the scoring identity, every row
assert np.allclose(Z.mean(axis=0), 0, atol=1e-12) and np.allclose(Z.std(axis=0), 1)
assert np.allclose(np.diag(R), 1) and np.allclose(R, R.T)
assert np.allclose(R, np.corrcoef(X, rowvar=False))                # same as numpy's correlation
assert near(mu[0], 17.510) and near(sd[0], 5.186)
assert near(raw_share[0], 0.4572, 5e-5) and int(np.argmax(raw_share)) == 0
assert f"{raw_share[0]:.1%}" == "45.7%" and f"{raw_share[8]:.1%}" == "0.3%"


# --- 2. two columns first: PCA on paper -------------------------------------------------------
r = R[4, 9]                                # offensive rebounds and blocks per 36
C2 = np.array([[1.0, r], [r, 1.0]])
q = (1 - r) / (1 + r)                      # how fast power iteration closes in
steps = []
v = np.array([1.0, 0.0])
for k in range(1, 7):
    v = C2 @ v
    v /= np.linalg.norm(v)
    angle = math.degrees(math.atan2(v[1], v[0]))
    steps.append((k, v.copy(), angle, 45 - angle, math.tan(math.radians(45 - angle))))
w2, W2 = np.linalg.eigh(C2)

with sdt.snippet("two"):
    print(f"offensive rebounds vs blocks per 36, {len(d)} players: r = {r:.4f}")
    print(f"eigenvalues 1 + r = {1 + r:.4f} and 1 - r = {1 - r:.4f}; "
          f"PC1 holds (1 + r) / 2 = {(1 + r) / 2:.2%}")
    print("eigenvectors (0.7071, 0.7071) and (0.7071, -0.7071): the two 45-degree lines")
    print(f"power iteration from (1, 0); q = (1 - r) / (1 + r) = {q:.5f}")
    print("step  vector               angle   gap to 45   tan(gap)   q^step")
    for k, vec, angle, gap, tg in steps:
        print(f"{k:4d}  ({vec[0]:.5f}, {vec[1]:.5f})  {angle:6.3f}  {gap:9.3f}  {tg:9.5f}  {q ** k:8.5f}")

assert near(r, 0.5772, 5e-5) and abs(r - np.corrcoef(d.oreb36, d.blk36)[0, 1]) < 1e-12
assert np.allclose(w2, [1 - r, 1 + r], atol=1e-12)                # eigenvalues 1 - r and 1 + r
assert np.allclose(np.abs(W2), 1 / math.sqrt(2), atol=1e-12)      # 45-degree axes, whatever r is
assert near((1 + r) / 2, 0.7886, 5e-5) and near(q, 0.26805, 5e-6)
assert all(abs(tg - q ** k) < 1e-12 for k, _v, _a, _g, tg in steps)   # tan(gap) = q^k exactly
assert [f"{a:.3f}" for _k, _v, a, _g, _t in steps[:4]] == ["29.995", "40.890", "43.897", "44.704"]
assert f"{steps[0][3]:.3f}" == "15.005" and f"{steps[1][3]:.2f}" == "4.11"


# --- 3. ten columns: power iteration with deflation -------------------------------------------
vals, V, iters = pca_by_hand(R)
share = vals / vals.sum()
np_vals = np.linalg.eigvalsh(R)[::-1]
sv = np.linalg.svd(Z, compute_uv=False)
raw_vals, raw_V, _ = pca_by_hand(np.cov(X, rowvar=False, bias=True))
raw_V = orient(raw_V)

with sdt.snippet("solver"):
    print("power iteration with deflation on the 10 x 10 correlation matrix")
    print("component  eigenvalue    share   cumulative   iterations")
    for k in range(10):
        print(f"{k + 1:9d}  {vals[k]:10.5f}  {share[k]:7.2%}  {share[:k + 1].sum():10.2%}  {iters[k]:11d}")
    print(f"eigenvalues sum to {vals.sum():.6f}, the number of standardized columns")
    print(f"largest gap to numpy.linalg.eigvalsh: {np.abs(vals - np_vals).max():.1e}; "
          f"to the SVD of Z: {np.abs(vals - sv ** 2 / len(Z)).max():.1e}")
    print(f"the same solver on the raw rates: PC1 takes {raw_vals[0] / raw_vals.sum():.2%} "
          "of the variance,")
    print(f"  and its loading on points is {raw_V[0, 0]:.3f}")

assert np.abs(vals - np_vals).max() < 1e-15 and np.abs(vals - sv ** 2 / len(Z)).max() < 1e-15
assert np.allclose(V.T @ V, np.eye(10), atol=1e-9)                # orthonormal axes
assert np.allclose(R @ V, V * vals, atol=1e-9)                    # R v = lambda v, all ten
assert abs(vals.sum() - 10) < 1e-12
assert all(vals[k] > vals[k + 1] for k in range(9))
assert [f"{x:.5f}" for x in vals[:3]] == ["3.50811", "2.83836", "1.20765"]
assert f"{vals[9]:.5f}" == "0.01801"
assert f"{share[0]:.2%}" == "35.08%" and f"{share[1]:.2%}" == "28.38%" and f"{share[9]:.2%}" == "0.18%"
assert f"{share[:2].sum():.2%}" == "63.46%" and f"{share[:3].sum():.2%}" == "75.54%"
assert max(iters) == iters[7] and 2000 < iters[7] < 2200          # PC8 and PC9 are nearly tied
assert all(it < 600 for k, it in enumerate(iters) if k != 7)
assert f"{vals[7]:.5f}" == "0.19592" and f"{vals[8]:.5f}" == "0.19382"
assert near(raw_vals[0] / raw_vals.sum(), 0.6539, 5e-5) and near(raw_V[0, 0], 0.824)


# --- 4. read the two big axes --------------------------------------------------------------------
V = orient(V)
S = Z @ V                                  # every player's score on every component
usage_r = [float(np.corrcoef(S[:, k], d.usg_pct)[0, 1]) for k in range(3)]
d["pc1"], d["pc2"] = S[:, 0], S[:, 1]
by1, by2 = d.sort_values("pc1"), d.sort_values("pc2")
i_curry = int(d.index[d.player == "Stephen Curry"][0])
z_show, l_show = np.round(Z[i_curry], 4), np.round(V[:, 0], 4)
products = np.round(z_show * l_show, 4)

with sdt.snippet("axes"):
    print("loadings, each axis flipped so its largest loading is positive")
    print(f"{'':22s}{'PC1':>8s}{'PC2':>8s}{'PC3':>8s}")
    for j, name in enumerate(NAMES):
        print(f"{name:22s}{V[j, 0]:+8.3f}{V[j, 1]:+8.3f}{V[j, 2]:+8.3f}")
    print("correlation with the league's usage rate, which is not one of the ten columns:")
    print(f"  PC1 {usage_r[0]:+.3f}, PC2 {usage_r[1]:+.3f}, PC3 {usage_r[2]:+.3f}")
    for label, frame, col in (("PC1 highest", by1[::-1], "pc1"), ("PC1 lowest ", by1, "pc1"),
                              ("PC2 highest", by2[::-1], "pc2"), ("PC2 lowest ", by2, "pc2")):
        print(f"{label}: " + ", ".join(f"{p} {s:+.2f}" for p, s in zip(frame.player[:3], frame[col][:3])))
    print("worked example, Stephen Curry's PC1 score = his ten z-scores times the PC1 loadings")
    print(f"{'':22s}{'z':>9s}{'loading':>10s}{'product':>10s}")
    for j, name in enumerate(NAMES):
        print(f"{name:22s}{z_show[j]:+9.4f}{l_show[j]:+10.4f}{products[j]:+10.4f}")
    print(f"{'sum':22s}{'':19s}{products.sum():+10.4f}   (unrounded: {S[i_curry, 0]:+.4f})")

assert near(usage_r[0], 0.910) and abs(usage_r[1]) < 0.25 and abs(usage_r[2]) < 0.2
assert near(V[0, 0], 0.469) and near(V[3, 0], 0.470) and near(V[7, 0], 0.448)   # volume
assert near(V[4, 1], 0.536) and near(V[9, 1], 0.457) and near(V[2, 1], -0.452)  # inside vs outside
assert near(V[8, 2], 0.760) and near(V[6, 2], 0.402)                             # steals, assists
assert list(by1.player[::-1][:3]) == ["Giannis Antetokounmpo", "Luka Dončić", "Nikola Jokić"]
assert list(by1.player[:3]) == ["Nicolas Batum", "Spencer Jones", "Dean Wade"]
assert (d.usg_pct.rank()[by1.index[:3]] <= 5).all()                   # three of the five lowest usage rates
assert list(by2.player[::-1][:3]) == ["Mitchell Robinson", "Robert Williams III", "Yves Missi"]
assert list(by2.player[:3]) == ["Stephen Curry", "LaMelo Ball", "Darius Garland"]
assert near(S[i_curry, 0], 2.8068) and near(S[i_curry, 1], -2.7409)
assert near(float(products.sum()), S[i_curry, 0], 5e-4)          # the rounded hand sum holds
assert f"{products.sum():.4f}" == "2.8067" and f"{S[i_curry, 0]:.4f}" == "2.8068"
assert near(by1.pc1.iloc[-1], 7.07, 5e-3) and near(by2.pc2.iloc[-1], 5.84, 5e-3)


# --- 5. how many components are real? -----------------------------------------------------------
kaiser = int((vals > 1).sum())
null = shuffled_eigenvalues(Z, SHUFFLES, SEED)
keep, cut = parallel_keep(vals, null)

seasons = sorted(df.season.unique())
fits, sweep = {}, []
for s in seasons:
    ds, Zs, ws, Ws = fit_season(df, s)
    k_s, cut_s = parallel_keep(ws, shuffled_eigenvalues(Zs, SHUFFLES, SEED))
    fits[s] = (ds, Zs, ws, Ws)
    sweep.append((s, len(ds), ws[2], cut_s[2], int((ws > 1).sum()), k_s))
sweep = pd.DataFrame(sweep, columns=["season", "n", "third", "cut3", "kaiser", "parallel"])

# a real trait should come back next season: score both seasons on the fixed 2025-26 axes
W26 = fits[SEASON][3]
pers, raw_pers, ts_pers, pair_n = [], [], [], []
for a, b in zip(seasons[:-1], seasons[1:]):
    da, Za, _, _ = fits[a]
    db, Zb, _, _ = fits[b]
    both = np.intersect1d(da.player_id, db.player_id)
    ia = pd.Series(da.index, index=da.player_id)[both].to_numpy()
    ib = pd.Series(db.index, index=db.player_id)[both].to_numpy()
    Sa, Sb = Za[ia] @ W26, Zb[ib] @ W26
    pers.append([np.corrcoef(Sa[:, k], Sb[:, k])[0, 1] for k in range(10)])
    raw_pers.append([np.corrcoef(da[c].to_numpy()[ia], db[c].to_numpy()[ib])[0, 1] for c in COLS])
    ts_a = (da.pts / (2 * (da.fga + 0.44 * da.fta))).to_numpy()[ia]
    ts_b = (db.pts / (2 * (db.fga + 0.44 * db.fta))).to_numpy()[ib]
    ts_pers.append(np.corrcoef(ts_a, ts_b)[0, 1])
    pair_n.append(len(both))
pers, raw_pers = np.array(pers), np.array(raw_pers)
med, raw_med = np.median(pers, axis=0), np.median(raw_pers, axis=0)

with sdt.snippet("howmany"):
    print(f"Kaiser's rule (eigenvalue above 1) keeps {kaiser}")
    print(f"parallel analysis, {SHUFFLES:,} shuffles of {SEASON}:")
    print("component  eigenvalue  95th pct of shuffled  kept")
    for k in range(10):
        print(f"{k + 1:9d}  {vals[k]:10.4f}  {cut[k]:20.4f}  {'yes' if k < keep else 'no'}")
    print(f"22 seasons, {MIN_MINUTES:,}+ minutes: Kaiser keeps 3 in {int((sweep.kaiser == 3).sum())}; "
          f"parallel analysis keeps 3 in {int((sweep.parallel == 3).sum())} and 2 in "
          f"{int((sweep.parallel == 2).sum())}")
    print("season   3rd eig  cut   keeps | season   3rd eig  cut   keeps")
    half = len(sweep) // 2
    for i in range(half):
        a, b = sweep.iloc[i], sweep.iloc[i + half]
        print(f"{a.season}  {a.third:6.3f} {a.cut3:6.3f}  {a.parallel:3d}   | "
              f"{b.season}  {b.third:6.3f} {b.cut3:6.3f}  {b.parallel:3d}")
    print(f"next-season persistence on the {SEASON} axes, median r over {len(pers)} season pairs")
    print(f"  ({min(pair_n)}-{max(pair_n)} players each):")
    print("  " + "  ".join(f"PC{k + 1} {med[k]:.3f}" for k in range(5)))
    print("  " + "  ".join(f"PC{k + 1} {med[k]:.3f}" for k in range(5, 10)))
    print(f"  steadiest single stat: {NAMES[int(np.argmax(raw_med))]} {raw_med.max():.3f}; "
          f"points {raw_med[0]:.3f}; true shooting % {np.median(ts_pers):.3f}")
    print(f"  PC2 beats all ten single stats in {int((pers[:, 1] > raw_pers.max(axis=1)).sum())} "
          f"of {len(pers)} pairs; PC1 beats points in {int((pers[:, 0] > raw_pers[:, 0]).sum())}")

assert kaiser == 3 and keep == 3
assert near(cut[0], 1.3946, 5e-4) and near(cut[2], 1.1884, 5e-4)
assert vals[2] > cut[2] and vals[2] - cut[2] < 0.03                # PC3 clears it by a hair
assert vals[3] < cut[3]
assert list(sweep.season) == seasons and (sweep.kaiser == 3).all()
assert int((sweep.parallel == 3).sum()) == 12 and int((sweep.parallel == 2).sum()) == 10
assert (sweep[sweep.season >= "2016-17"].parallel == 3).all()      # every season since 2016-17
assert int((sweep[sweep.season < "2016-17"].parallel == 2).sum()) == 10
assert abs(sweep.iloc[-1].cut3 - cut[2]) < 1e-12                   # the sweep reruns 2025-26 exactly
assert len(pers) == 21 and min(pair_n) == 177 and max(pair_n) == 215
assert near(med[0], 0.897) and near(med[1], 0.968) and near(med[9], 0.584)
assert int(np.argmax(med)) == 1 and int(np.argmin(med)) == 9
assert NAMES[int(np.argmax(raw_med))] == "off. rebounds" and near(raw_med.max(), 0.942)
assert near(raw_med[0], 0.872) and near(float(np.median(ts_pers)), 0.629)
assert (pers[:, 1] > raw_pers.max(axis=1)).all()                   # PC2, 21 of 21
assert (pers[:, 0] > raw_pers[:, 0]).all()                         # PC1 beats points, 21 of 21
assert med[1] > raw_med.max() and med[9] < float(np.median(ts_pers))
assert all(m > 0.67 for m in med[3:9])                             # small, and still repeatable


# --- 6. the smallest component is the shooting ----------------------------------------------------
shot, (per2, per3, perft) = shooting_points(d)
A = np.column_stack([np.ones(len(d)), d.fg2a36, d.fg3a36, d.fta36])
coef, *_ = np.linalg.lstsq(A, d.pts36.to_numpy(), rcond=None)
resid = d.pts36.to_numpy() - A @ coef
r2_attempts = 1 - resid.var() / d.pts36.var(ddof=0)
corr_shot = np.array([np.corrcoef(S[:, k], shot)[0, 1] for k in range(10)])
shot_share = corr_shot ** 2
ts = (d.pts / (2 * (d.fga + 0.44 * d.fta))).to_numpy()
K = 3
X3 = mu + sd * (S[:, :K] @ V[:, :K].T)    # every player rebuilt from three numbers
shot3 = X3[:, 0] - (per2 * X3[:, 1] + per3 * X3[:, 2] + perft * X3[:, 3])
SHOW = ["Stephen Curry", "Nikola Jokić", "Bam Adebayo"]
idx_show = [int(d.index[d.player == p][0]) for p in SHOW]

last = []
for s in seasons:
    ds, Zs, ws, Ws = fits[s]
    sh, _ = shooting_points(ds)
    c2 = np.array([np.corrcoef(Zs @ Ws[:, k], sh)[0, 1] for k in range(10)]) ** 2
    last.append((s, c2[9], int(np.argmax(c2)), ws[9] / 10, c2.sum()))
last = pd.DataFrame(last, columns=["season", "share", "best", "var", "total"])

floors = {}
for floor in (500, 1500):                  # does the minutes floor matter?
    ds, Zs, ws, Ws = fit_season(df, SEASON, floor)
    sh, _ = shooting_points(ds)
    c2 = np.array([np.corrcoef(Zs @ Ws[:, k], sh)[0, 1] for k in range(10)]) ** 2
    floors[floor] = (len(ds), c2[9], int(np.argmax(c2)))

rng = np.random.default_rng(SEED)
boot_share, boot_gap = np.empty(BOOT), np.empty(BOOT)
for i in range(BOOT):
    rows = rng.integers(0, len(d), len(d))
    Zb, _, _ = standardize(X[rows])
    wb, Wb = np.linalg.eigh(Zb.T @ Zb / len(Zb))                    # ascending: column 0 is smallest
    sb, _ = shooting_points(d.iloc[rows])
    boot_share[i] = np.corrcoef(Zb @ Wb[:, 0], sb)[0, 1] ** 2
    boot_gap[i] = wb[-1] - wb[-2]
share_ci = np.percentile(boot_share, [2.5, 97.5])
gap_ci = np.percentile(boot_gap, [2.5, 97.5])

with sdt.snippet("shooting"):
    print(f"PC10: eigenvalue {vals[9]:.5f}, {share[9]:.2%} of the variance")
    print("  loadings: " + ", ".join(f"{NAMES[j]} {V[j, 9]:+.3f}" for j in range(2)) + ",")
    print("            " + ", ".join(f"{NAMES[j]} {V[j, 9]:+.3f}" for j in range(2, 4)))
    print(f"  the other six loadings are all within {np.abs(V[4:, 9]).max():.3f} of zero")
    print(f"points per 36 regressed on the three attempt rates: R^2 = {r2_attempts:.4f}")
    print(f"group conversion: {per2:.4f} points per 2-pt attempt, {per3:.4f} per 3-pt attempt,")
    print(f"  {perft:.4f} per free throw")
    print(f"shooting points per 36 = points - what those attempts score at those rates "
          f"(sd {shot.std():.3f})")
    print("component  share of variance  r with shooting points  share of shooting signal")
    for k in range(10):
        print(f"{k + 1:9d}  {share[k]:17.2%}  {corr_shot[k]:+22.3f}  {shot_share[k]:24.2%}")
    print(f"{'total':9s}  {share.sum():17.2%}  {'':22s}  {shot_share.sum():24.2%}")
    print(f"keep {K} components ({share[:K].sum():.2%} of the variance): "
          f"{shot_share[:K].sum():.2%} of the shooting signal survives")
    print(f"{'':22s}{'pts/36':>8s}{'from 3':>8s}{'shooting':>10s}{'from 3':>8s}")
    for p, i in zip(SHOW, idx_show):
        print(f"{p:22s}{d.pts36[i]:8.2f}{X3[i, 0]:8.2f}{shot[i]:+10.2f}{shot3[i]:+8.2f}")
    print(f"true shooting % vs PC10: r = {np.corrcoef(S[:, 9], ts)[0, 1]:+.3f}")
    print(f"22 seasons: the last component holds the largest share in {int((last.best == 9).sum())};")
    print(f"  its share runs {last.share.min():.1%} to {last.share.max():.1%}, its variance never "
          f"above {last['var'].max():.2%}")
    print(f"{SEASON} with other minute floors:")
    for f, (n, s2, _b) in floors.items():
        print(f"  {f:,}+ minutes, {n} players: PC10 holds {s2:.1%}")
    print(f"bootstrap, {BOOT:,} resamples of the {len(d)} players: PC10's share "
          f"{share_ci[0]:.1%} to {share_ci[1]:.1%}")

assert near(V[0, 9], 0.695) and near(V[1, 9], -0.487) and near(V[2, 9], -0.438) and near(V[3, 9], -0.270)
assert np.abs(V[4:, 9]).max() < 0.10
assert near(r2_attempts, 0.9543, 5e-5)
assert near(per2, 1.1053, 5e-5) and near(per3, 1.0908, 5e-5) and near(perft, 0.7893, 5e-5)
assert near(shot.std(), 1.134)
assert abs(shot_share.sum() - 1) < 1e-9                            # shooting points sit in the span
assert int(np.argmax(shot_share)) == 9 and near(shot_share[9], 0.6986, 5e-5)
assert near(corr_shot[9], 0.836) and f"{shot_share[9]:.2%}" == "69.86%"
assert f"{shot_share[:K].sum():.2%}" == "15.80%" and f"{shot_share[:9].sum():.2%}" == "30.14%"
assert near(np.var(shot3) / np.var(shot), shot_share[:K].sum(), 1e-9)   # rebuilt = projected
assert near(shot[idx_show[0]], 2.478) and near(shot3[idx_show[0]], -0.229)
assert near(d.pts36[idx_show[0]], 30.94, 5e-3) and near(X3[idx_show[0], 0], 27.19, 5e-3)
assert near(shot[idx_show[1]], 2.754) and near(shot3[idx_show[1]], -0.386)
assert near(shot[idx_show[2]], -1.912) and near(shot3[idx_show[2]], 0.377)
assert near(float(np.corrcoef(S[:, 9], ts)[0, 1]), 0.718)
assert near(shot_share[1], 0.0636, 5e-5) and near(shot_share[2], 0.0924, 5e-5)
assert near(shot_share[8], 0.0997, 5e-5)
assert max(shot_share[k] for k in (0, 3, 4, 5, 6, 7)) < 0.02         # the rest, under 2% each
assert all(b == 9 and s2 > 0.7 for _n, s2, b in floors.values())   # the floor does not matter
assert floors[500][0] == 378 and floors[1500][0] == 164
assert f"{floors[500][1]:.1%}" == "71.7%" and f"{floors[1500][1]:.1%}" == "73.8%"
assert (last.best == 9).all() and np.allclose(last.total, 1, atol=1e-9)
assert near(last.share.min(), 0.6617, 5e-5) and near(last.share.max(), 0.8528, 5e-5)
assert last.loc[last.share.idxmin(), "season"] == "2018-19" and last["var"].max() < 0.0027
assert f"{share_ci[0]:.1%}" == "57.5%" and f"{share_ci[1]:.1%}" == "76.7%"
assert share_ci[0] > 0.5                                           # over half, even at the low end


# --- 7. do the axes hold still? ----------------------------------------------------------------
def pc1_kind(Ws):
    """Name a season's PC1 by which loads harder on it: offensive rebounds or points."""
    reb, pts = abs(Ws[4, 0]), abs(Ws[0, 0])
    if abs(reb - pts) < 0.01:
        return "blend"
    return "inside-outside" if reb > pts else "volume"


kinds = {s: pc1_kind(fits[s][3]) for s in seasons}
P05, P26 = fits["2004-05"][3][:, :2], fits[SEASON][3][:, :2]
cosines = np.linalg.svd(P05.T @ P26, compute_uv=False)        # cosines of the principal angles
angles = np.degrees(np.arccos(np.clip(cosines, -1, 1)))
w05, W05 = fits["2004-05"][2], fits["2004-05"][3]

with sdt.snippet("stability"):
    print(f"2004-05: PC1 {w05[0]:.3f} = inside-outside: "
          + ", ".join(f"{NAMES[j]} {W05[j, 0]:+.3f}" for j in (4, 5, 9)) + ",")
    print("                                  "
          + ", ".join(f"{NAMES[j]} {W05[j, 0]:+.3f}" for j in (6, 2)))
    print(f"         PC2 {w05[1]:.3f} = volume: "
          + ", ".join(f"{NAMES[j]} {W05[j, 1]:+.3f}" for j in (0, 3)) + ",")
    print("                          "
          + ", ".join(f"{NAMES[j]} {W05[j, 1]:+.3f}" for j in (1, 7)))
    print(f"{SEASON}: PC1 {vals[0]:.3f} = volume; PC2 {vals[1]:.3f} = inside-outside")
    print("what PC1 is, season by season (offensive-rebound vs points loading):")
    for kind in ("inside-outside", "blend", "volume"):
        ss = [s for s in seasons if kinds[s] == kind]
        print(f"  {kind:15s}{len(ss):3d}  {ss[0]} to {ss[-1]}" if len(ss) > 2 else
              f"  {kind:15s}{len(ss):3d}  {', '.join(ss)}")
    print(f"the plane of PC1 and PC2, 2004-05 vs {SEASON}:")
    print(f"  principal angles {angles[0]:.1f} and {angles[1]:.1f} degrees "
          f"(cosines {cosines[0]:.3f}, {cosines[1]:.3f})")
    print(f"{SEASON} gap between the first two eigenvalues: {vals[0] - vals[1]:.3f}, "
          f"bootstrap 95% interval {gap_ci[0]:.3f} to {gap_ci[1]:.3f}")

ins = [s for s in seasons if kinds[s] == "inside-outside"]
vol = [s for s in seasons if kinds[s] == "volume"]
assert ins == seasons[:12] and ins[-1] == "2015-16"
assert [s for s in seasons if kinds[s] == "blend"] == ["2016-17", "2017-18"]
assert vol == seasons[14:] and len(vol) == 8 and vol[0] == "2018-19"
assert near(w05[0], 3.608) and near(w05[1], 2.949)
assert near(W05[4, 0], 0.476) and near(W05[0, 1], 0.511)
assert near(cosines[0], 0.967) and near(cosines[1], 0.952)
assert f"{angles[0]:.1f}" == "14.7" and f"{angles[1]:.1f}" == "17.8"
assert near(vals[0] - vals[1], 0.670) and gap_ci[0] > 0
assert near(gap_ci[0], 0.382) and near(gap_ci[1], 1.106)


# --- 8. the exhibit ------------------------------------------------------------------------------
ORANGE, BLUE, INK, GREY = sdt.sport_color("basketball"), "#2C5E8A", "#20242B", "#8A8577"
SHORT = ["PTS", "2PA", "3PA", "FTA", "OREB", "DREB", "AST", "TOV", "STL", "BLK"]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.6, 5.4), gridspec_kw={"width_ratios": [1.05, 1]})

ax1.scatter(S[:, 0], S[:, 1], s=14, color=ORANGE, alpha=0.55, linewidths=0, zorder=2)
ax1.axhline(0, color=GREY, lw=0.8, zorder=1)
ax1.axvline(0, color=GREY, lw=0.8, zorder=1)
SCALE = 6.0
NUDGE = {"PTS": (0.32, 0.30), "TOV": (0.32, -0.32), "BLK": (0.42, -0.12)}   # PTS and TOV nearly coincide
for j, lab in enumerate(SHORT):
    x, y = SCALE * V[j, 0], SCALE * V[j, 1]
    if math.hypot(x, y) < 1.0:
        continue                                       # steals sit almost on the origin here
    ax1.annotate("", xy=(x, y), xytext=(0, 0),
                 arrowprops=dict(arrowstyle="-|>", color=BLUE, lw=1.1, alpha=0.85), zorder=3)
    tx, ty = (x + NUDGE[lab][0], y + NUDGE[lab][1]) if lab in NUDGE else (x * 1.12, y * 1.12)
    ax1.text(tx, ty, lab, color=BLUE, fontsize=8, ha="center", va="center", zorder=4)
LABELS = {"Giannis Antetokounmpo": (-10, 9), "Luka Dončić": (6, -10), "Stephen Curry": (6, -9),
          "Mitchell Robinson": (8, 2), "Nicolas Batum": (-6, -11), "Rudy Gobert": (8, -4)}
for p, (dx, dy) in LABELS.items():
    i = int(d.index[d.player == p][0])
    ax1.scatter([S[i, 0]], [S[i, 1]], s=26, facecolor="none", edgecolor=INK, linewidths=1.0, zorder=5)
    ax1.annotate(p, (S[i, 0], S[i, 1]), xytext=(dx, dy), textcoords="offset points", fontsize=8,
                 color=INK, ha="left" if dx > 0 else "right", zorder=6)
ax1.set_xlabel(f"PC1, volume ({share[0]:.1%} of the variance)")
ax1.set_ylabel(f"PC2, inside-outside ({share[1]:.1%})")
ax1.set_title(f"The two big axes: {len(d)} players, {SEASON}", fontsize=12)
ax1.set_xlim(-6.2, 8.6)
ax1.set_ylim(-4.6, 7.2)

ks = np.arange(1, 11)
ax2.bar(ks - 0.2, 100 * share, width=0.4, color=GREY, label="share of the variance", zorder=2)
ax2.bar(ks + 0.2, 100 * shot_share, width=0.4, color=ORANGE, label="share of the shooting signal",
        zorder=2)
ax2.set_xticks(ks)
ax2.set_xticklabels([f"PC{k}" for k in ks], fontsize=8.5)
ax2.set_ylabel("percent")
ax2.set_ylim(0, 80)
ax2.set_title("Ranked by variance, the shooting comes last", fontsize=12)
ax2.annotate(f"PC10: {share[9]:.2%} of the variance,\n{shot_share[9]:.1%} of the shooting signal",
             xy=(10.2, 100 * shot_share[9]), xytext=(5.6, 62), fontsize=9, color=INK, ha="center",
             arrowprops=dict(arrowstyle="->", color=INK, lw=0.9))
ax2.legend(handles=[Patch(facecolor=GREY, label="share of the variance"),
                    Patch(facecolor=ORANGE, label="share of the shooting signal")],
           loc="upper left", fontsize=8.6, frameon=False)
fig.tight_layout()
sdt.save_fig(fig, "pca_box_score",
             source="NBA Stats API (stats.nba.com) via nba_api, 2025-26 regular season",
             asof="September 2026")
print("\nall asserts passed")