analysis.py

"""Statistical and TabICLv2-based analysis of player pay versus on-field performance.

All numbers produced here are computed from the loaded dataset. Nothing is
looked up externally. The module is independent of the user interface so it
can be unit-tested directly.
"""

from __future__ import annotations

import hashlib
import json
import os
import threading
from typing import Any

import numpy as np
import pandas as pd
from scipy import stats

STAT_COLUMNS = [
    "PA", "HR", "R", "RBI", "SB", "BB%", "K%", "AVG", "OBP", "SLG",
    "wOBA", "wRC+", "BsR", "Off", "Def", "WAR",
]
LOWER_IS_BETTER = {"K%"}
PERCENT_COLUMNS = {"BB%", "K%"}
THREE_DECIMAL_COLUMNS = {"AVG", "OBP", "SLG", "wOBA"}
INTEGER_COLUMNS = {"PA", "HR", "R", "RBI", "SB", "wRC+"}
SALARY_COLUMN = "AAV"
NAME_COLUMN = "Name"
TEAM_COLUMN = "Team"
REQUIRED_COLUMNS = [NAME_COLUMN, TEAM_COLUMN, SALARY_COLUMN] + STAT_COLUMNS

# Players paid at or below this amount are on a league-minimum type salary in
# this file (the lowest value in the data is 780,000). Their pay is fixed by
# service-time rules rather than by performance, so the salary model is fitted
# on the players above the threshold ("market contracts") only. The split is
# also used for descriptive comparisons.
MARKET_CONTRACT_THRESHOLD = 1_000_000

DEFAULT_DATA_PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "players.csv")
VALIDATION_CACHE_PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "loo_validation.json")

QUANTILE_GRID = [round(a, 2) for a in np.arange(0.01, 1.00, 0.01)]
REPORT_QUANTILES = [0.10, 0.25, 0.50, 0.75, 0.90]
N_COMPARABLES = 5


# ---------------------------------------------------------------------------
# Data loading
# ---------------------------------------------------------------------------

def load_players(path: str = DEFAULT_DATA_PATH) -> pd.DataFrame:
    """Load and validate the player table."""
    df = pd.read_csv(path)
    missing = [c for c in REQUIRED_COLUMNS if c not in df.columns]
    if missing:
        raise ValueError(f"Dataset is missing required columns: {missing}")
    numeric = STAT_COLUMNS + [SALARY_COLUMN]
    for col in numeric:
        df[col] = pd.to_numeric(df[col], errors="coerce")
    bad = df[numeric].isna().any(axis=1)
    if bad.any():
        raise ValueError(f"{int(bad.sum())} rows have non-numeric or missing values in {numeric}")
    if (df[SALARY_COLUMN] <= 0).any():
        raise ValueError("All salaries must be positive")
    if df[NAME_COLUMN].duplicated().any():
        dupes = df.loc[df[NAME_COLUMN].duplicated(), NAME_COLUMN].tolist()
        raise ValueError(f"Duplicate player names: {dupes}")
    if len(df) < 10:
        raise ValueError("At least 10 players are required")
    df = df.reset_index(drop=True)
    return df


def dataset_fingerprint(df: pd.DataFrame) -> str:
    """Stable hash of the analysed columns, used to validate cached results."""
    payload = df[REQUIRED_COLUMNS].to_csv(index=False).encode("utf-8")
    return hashlib.sha256(payload).hexdigest()[:16]


def dataset_overview(df: pd.DataFrame) -> dict[str, Any]:
    """Whole-file facts used in the interface text."""
    market = df[df[SALARY_COLUMN] > MARKET_CONTRACT_THRESHOLD]
    rho_all, p_all = stats.spearmanr(df["WAR"], df[SALARY_COLUMN])
    rho_market, p_market = stats.spearmanr(market["WAR"], market[SALARY_COLUMN])
    return {
        "n_players": int(len(df)),
        "n_market_contracts": int(len(market)),
        "n_minimum_contracts": int(len(df) - len(market)),
        "min_pa": int(df["PA"].min()),
        "max_pa": int(df["PA"].max()),
        "spearman_war_salary_all": round(float(rho_all), 2),
        "spearman_war_salary_all_p": round(float(p_all), 3),
        "spearman_war_salary_market": round(float(rho_market), 2),
        "spearman_war_salary_market_p": round(float(p_market), 3),
    }


def player_index(df: pd.DataFrame, name: str) -> int:
    matches = df.index[df[NAME_COLUMN] == name].tolist()
    if not matches:
        raise KeyError(f"Player not found: {name}")
    return int(matches[0])


# ---------------------------------------------------------------------------
# Descriptive statistics
# ---------------------------------------------------------------------------

def percentile_rank(others: np.ndarray, value: float, higher_is_better: bool = True) -> float:
    """Share of the other players this value beats, in percent (ties count half)."""
    others = np.asarray(others, dtype=float)
    if others.size == 0:
        return float("nan")
    if higher_is_better:
        below = np.sum(others < value)
    else:
        below = np.sum(others > value)
    ties = np.sum(others == value)
    return float(100.0 * (below + 0.5 * ties) / others.size)


def format_stat(col: str, value: float) -> str:
    if col in INTEGER_COLUMNS:
        return f"{int(round(value))}"
    if col in PERCENT_COLUMNS:
        return f"{100 * value:.1f}%"
    if col in THREE_DECIMAL_COLUMNS:
        return f"{value:.3f}"
    return f"{value:.1f}"


def format_money(value: float) -> str:
    return f"${value:,.0f}"


def format_money_short(value: float) -> str:
    if abs(value) >= 1_000_000:
        return f"${value / 1_000_000:.2f}M"
    return f"${value / 1_000:.0f}K"


def performance_table(df: pd.DataFrame, idx: int) -> pd.DataFrame:
    """Per-statistic comparison of one player against every other player."""
    rows = []
    others_mask = df.index != idx
    for col in STAT_COLUMNS:
        value = float(df.at[idx, col])
        others = df.loc[others_mask, col].to_numpy(dtype=float)
        higher = col not in LOWER_IS_BETTER
        pct = percentile_rank(others, value, higher_is_better=higher)
        sd = others.std(ddof=1)
        z = (value - others.mean()) / sd if sd > 0 else 0.0
        if not higher:
            z = -z
        rank_series = df[col].rank(ascending=not higher, method="min")
        rows.append({
            "Statistic": col,
            "Player": format_stat(col, value),
            "Sample median": format_stat(col, float(np.median(others))),
            "Sample min": format_stat(col, float(others.min())),
            "Sample max": format_stat(col, float(others.max())),
            "Rank": f"{int(rank_series.iloc[idx])} of {len(df)}",
            "Percentile": round(pct, 1),
            "z-score": round(float(z), 2),
        })
    return pd.DataFrame(rows)


def salary_summary(df: pd.DataFrame, idx: int) -> dict[str, Any]:
    salary = float(df.at[idx, SALARY_COLUMN])
    others = df.loc[df.index != idx, SALARY_COLUMN].to_numpy(dtype=float)
    market = df[df[SALARY_COLUMN] > MARKET_CONTRACT_THRESHOLD]
    minimum = df[df[SALARY_COLUMN] <= MARKET_CONTRACT_THRESHOLD]
    rank = int(df[SALARY_COLUMN].rank(ascending=False, method="min").iloc[idx])
    market_others = market.loc[market.index != idx, SALARY_COLUMN].to_numpy(dtype=float)
    return {
        "salary": salary,
        "salary_rank": rank,
        "salary_percentile": round(percentile_rank(others, salary), 1),
        "salary_percentile_among_market": (
            round(percentile_rank(market_others, salary), 1) if salary > MARKET_CONTRACT_THRESHOLD else None
        ),
        "n_market_others": int(len(market_others)),
        "sample_median_salary": float(np.median(others)),
        "sample_mean_salary": float(others.mean()),
        "n_players": int(len(df)),
        "n_market_contracts": int(len(market)),
        "n_minimum_contracts": int(len(minimum)),
        "market_median_salary": float(market[SALARY_COLUMN].median()) if len(market) else float("nan"),
        "is_market_contract": bool(salary > MARKET_CONTRACT_THRESHOLD),
    }


def cost_per_war(df: pd.DataFrame, idx: int) -> dict[str, Any]:
    """Dollars of AAV per unit of WAR, compared with the market-contract players."""
    war = float(df.at[idx, "WAR"])
    salary = float(df.at[idx, SALARY_COLUMN])
    market = df[(df[SALARY_COLUMN] > MARKET_CONTRACT_THRESHOLD) & (df["WAR"] > 0) & (df.index != idx)]
    market_ratio = (market[SALARY_COLUMN] / market["WAR"]).to_numpy(dtype=float)
    result: dict[str, Any] = {
        "war": war,
        "market_median_cost_per_war": float(np.median(market_ratio)) if market_ratio.size else float("nan"),
        "n_market_positive_war": int(market_ratio.size),
    }
    if war > 0:
        ratio = salary / war
        result["cost_per_war"] = ratio
        result["cost_per_war_percentile_among_market"] = round(
            percentile_rank(market_ratio, ratio, higher_is_better=False), 1
        ) if market_ratio.size else float("nan")
    else:
        result["cost_per_war"] = None
        result["cost_per_war_percentile_among_market"] = None
    return result


def standardized_features(df: pd.DataFrame) -> np.ndarray:
    X = df[STAT_COLUMNS].to_numpy(dtype=float)
    mean = X.mean(axis=0)
    sd = X.std(axis=0, ddof=1)
    sd[sd == 0] = 1.0
    Z = (X - mean) / sd
    Z[:, [STAT_COLUMNS.index(c) for c in LOWER_IS_BETTER]] *= -1
    return Z


def comparable_players(df: pd.DataFrame, idx: int, k: int = N_COMPARABLES) -> tuple[pd.DataFrame, dict[str, Any]]:
    """The k market-contract players closest to the selected player in standardized statistics.

    Returns the display table and a summary of the comparables' salaries.
    """
    Z = standardized_features(df)
    dist = np.sqrt(((Z - Z[idx]) ** 2).sum(axis=1))
    dist[~model_training_mask(df, idx)] = np.inf
    order = np.argsort(dist)[:k]
    order = order[np.isfinite(dist[order])]
    rows = []
    for j in order:
        rows.append({
            "Player": df.at[j, NAME_COLUMN],
            "Team": df.at[j, TEAM_COLUMN],
            "AAV": format_money(float(df.at[j, SALARY_COLUMN])),
            "WAR": format_stat("WAR", float(df.at[j, "WAR"])),
            "wRC+": format_stat("wRC+", float(df.at[j, "wRC+"])),
            "OPS": f"{float(df.at[j, 'OBP'] + df.at[j, 'SLG']):.3f}",
            "Distance": f"{float(dist[j]):.2f}",
        })
    salaries = df.loc[order, SALARY_COLUMN].to_numpy(dtype=float)
    summary = {
        "comparable_rule": f"nearest players with AAV above {MARKET_CONTRACT_THRESHOLD:,} by standardized statistics",
        "comparable_names": df.loc[order, NAME_COLUMN].tolist(),
        "comparable_median_salary": float(np.median(salaries)) if salaries.size else float("nan"),
        "comparable_min_salary": float(salaries.min()) if salaries.size else float("nan"),
        "comparable_max_salary": float(salaries.max()) if salaries.size else float("nan"),
    }
    return pd.DataFrame(rows), summary


# ---------------------------------------------------------------------------
# TabICLv2 salary model
# ---------------------------------------------------------------------------

def _make_regressor(n_estimators: int, random_state: int = 42):
    from tabicl import TabICLRegressor

    return TabICLRegressor(
        n_estimators=n_estimators,
        device=os.environ.get("TABICL_DEVICE", "cpu"),
        random_state=random_state,
        verbose=False,
    )


def _model_n_estimators() -> int:
    return int(os.environ.get("TABICL_N_ESTIMATORS", "8"))


def fit_predict_log_salary(
    X_train: np.ndarray,
    y_train_log: np.ndarray,
    X_test: np.ndarray,
    n_estimators: int | None = None,
) -> np.ndarray:
    """Fit TabICLv2 on the training rows and return predicted log-salary quantiles.

    The result has shape (n_test, len(QUANTILE_GRID)); quantiles are made
    monotone so that interpolation across them is well defined.
    """
    reg = _make_regressor(n_estimators or _model_n_estimators())
    reg.fit(X_train, y_train_log)
    quantiles = np.asarray(reg.predict(X_test, output_type="quantiles", alphas=QUANTILE_GRID), dtype=float)
    return np.maximum.accumulate(quantiles, axis=1)


def percentile_within_distribution(actual_log: float, quantiles_log: np.ndarray) -> float:
    """Where the actual value falls inside the predicted distribution, in percent."""
    alphas = np.asarray(QUANTILE_GRID, dtype=float)
    if actual_log <= quantiles_log[0]:
        return float(100 * alphas[0])
    if actual_log >= quantiles_log[-1]:
        return float(100 * alphas[-1])
    return float(100 * np.interp(actual_log, quantiles_log, alphas))


def model_training_mask(df: pd.DataFrame, idx: int | None = None) -> np.ndarray:
    """Rows used to fit the salary model: market contracts, excluding the selected player."""
    mask = np.array(df[SALARY_COLUMN] > MARKET_CONTRACT_THRESHOLD, dtype=bool)
    if idx is not None:
        mask[idx] = False
    return mask


def model_salary_for_player(df: pd.DataFrame, idx: int, n_estimators: int | None = None) -> dict[str, Any]:
    """Fit TabICLv2 on the other market-contract players and predict the selected player's salary."""
    X = df[STAT_COLUMNS].to_numpy(dtype=float)
    y_log = np.log(df[SALARY_COLUMN].to_numpy(dtype=float))
    train = model_training_mask(df, idx)
    if train.sum() < 10:
        raise ValueError("Fewer than 10 market-contract players are available to fit the model")
    q = fit_predict_log_salary(X[train], y_log[train], X[[idx]], n_estimators)[0]
    actual = float(df.at[idx, SALARY_COLUMN])
    actual_log = float(np.log(actual))
    report = {}
    for alpha in REPORT_QUANTILES:
        report[f"q{int(round(alpha * 100)):02d}"] = float(np.exp(q[QUANTILE_GRID.index(round(alpha, 2))]))
    predicted_median = float(report["q50"])
    return {
        "n_train": int(train.sum()),
        "training_rule": f"players with AAV above {MARKET_CONTRACT_THRESHOLD:,}, excluding the selected player",
        "selected_player_is_market_contract": bool(actual > MARKET_CONTRACT_THRESHOLD),
        "n_estimators": n_estimators or _model_n_estimators(),
        "predicted_median_salary": predicted_median,
        "predicted_quantiles": report,
        "actual_salary": actual,
        "actual_percentile_in_predicted": round(percentile_within_distribution(actual_log, q), 1),
        "actual_over_predicted_median": actual / predicted_median,
        "inside_10_90": bool(report["q10"] <= actual <= report["q90"]),
        "inside_25_75": bool(report["q25"] <= actual <= report["q75"]),
    }


# ---------------------------------------------------------------------------
# Leave-one-out validation of the model on the whole file
# ---------------------------------------------------------------------------

def leave_one_out_validation(df: pd.DataFrame, n_estimators: int | None = None,
                             progress: Any = None) -> dict[str, Any]:
    """Hold out each market-contract player in turn, fit on the others, and score the predictions."""
    X_all = df[STAT_COLUMNS].to_numpy(dtype=float)
    y_all = np.log(df[SALARY_COLUMN].to_numpy(dtype=float))
    market_idx = np.flatnonzero(model_training_mask(df))
    n = len(market_idx)
    y_log = y_all[market_idx]
    pred_median = np.zeros(n)
    q10 = np.zeros(n)
    q25 = np.zeros(n)
    q75 = np.zeros(n)
    q90 = np.zeros(n)
    for k, i in enumerate(market_idx):
        train = model_training_mask(df, int(i))
        q = fit_predict_log_salary(X_all[train], y_all[train], X_all[[i]], n_estimators)[0]
        pred_median[k] = q[QUANTILE_GRID.index(0.50)]
        q10[k] = q[QUANTILE_GRID.index(0.10)]
        q25[k] = q[QUANTILE_GRID.index(0.25)]
        q75[k] = q[QUANTILE_GRID.index(0.75)]
        q90[k] = q[QUANTILE_GRID.index(0.90)]
        if progress is not None:
            progress(k + 1, n)
    ss_res = float(((y_log - pred_median) ** 2).sum())
    ss_tot = float(((y_log - y_log.mean()) ** 2).sum())
    rho, p_value = stats.spearmanr(y_log, pred_median)
    return {
        "fingerprint": dataset_fingerprint(df),
        "n_players": int(n),
        "training_rule": f"players with AAV above {MARKET_CONTRACT_THRESHOLD:,}",
        "n_estimators": n_estimators or _model_n_estimators(),
        "r2_log_salary": round(1 - ss_res / ss_tot, 3),
        "spearman_rho": round(float(rho), 3),
        "spearman_p_value": round(float(p_value), 4),
        "mean_abs_error_log": round(float(np.abs(y_log - pred_median).mean()), 3),
        "median_abs_pct_error": round(float(np.median(np.abs(np.exp(pred_median - y_log) - 1)) * 100), 1),
        "coverage_10_90": round(float(np.mean((y_log >= q10) & (y_log <= q90))) * 100, 1),
        "coverage_25_75": round(float(np.mean((y_log >= q25) & (y_log <= q75))) * 100, 1),
    }


def load_cached_validation(df: pd.DataFrame, path: str = VALIDATION_CACHE_PATH) -> dict[str, Any] | None:
    """Return the cached validation only if it was computed on exactly this dataset."""
    if not os.path.exists(path):
        return None
    try:
        with open(path, "r", encoding="utf-8") as fh:
            cached = json.load(fh)
    except (OSError, ValueError):
        return None
    if cached.get("fingerprint") != dataset_fingerprint(df):
        return None
    if cached.get("n_estimators") != _model_n_estimators():
        return None
    return cached


def save_validation(result: dict[str, Any], path: str = VALIDATION_CACHE_PATH) -> None:
    os.makedirs(os.path.dirname(path), exist_ok=True)
    with open(path, "w", encoding="utf-8") as fh:
        json.dump(result, fh, indent=2)


class ValidationRunner:
    """Runs the leave-one-out validation in a background thread and caches the result."""

    def __init__(self, df: pd.DataFrame, cache_path: str = VALIDATION_CACHE_PATH):
        self.df = df
        self.cache_path = cache_path
        self.result = load_cached_validation(df, cache_path)
        self.error: str | None = None
        self.completed = 0
        self.total = int(model_training_mask(df).sum())
        self._lock = threading.Lock()
        self._thread: threading.Thread | None = None

    def start(self) -> None:
        if self.result is not None or self._thread is not None:
            return
        self._thread = threading.Thread(target=self._run, name="loo-validation", daemon=True)
        self._thread.start()

    def _progress(self, done: int, total: int) -> None:
        with self._lock:
            self.completed = done
            self.total = total

    def _run(self) -> None:
        try:
            result = leave_one_out_validation(self.df, progress=self._progress)
            try:
                save_validation(result, self.cache_path)
            except OSError:
                pass
            with self._lock:
                self.result = result
        except Exception as exc:  # noqa: BLE001 - surfaced to the UI
            with self._lock:
                self.error = f"{type(exc).__name__}: {exc}"

    def status(self) -> dict[str, Any]:
        with self._lock:
            return {
                "result": self.result,
                "error": self.error,
                "completed": self.completed,
                "total": self.total,
                "running": self._thread is not None and self._thread.is_alive(),
            }


# ---------------------------------------------------------------------------
# Verdict
# ---------------------------------------------------------------------------

def assessment(perf: pd.DataFrame, salary: dict[str, Any], model: dict[str, Any],
               comps: dict[str, Any]) -> dict[str, Any]:
    """Rule-based reading of the evidence. Every threshold is stated in the method notes.

    The TabICLv2 result leads: where the actual salary falls in the salary distribution the model
    predicts for these statistics. The rank comparison can turn a clear reading into a mixed one.
    The comparables are reported as context.
    """
    war_pct = float(perf.loc[perf["Statistic"] == "WAR", "Percentile"].iloc[0])
    mean_pct = float(perf["Percentile"].mean())
    salary_pct = float(salary["salary_percentile"])
    p_model = float(model["actual_percentile_in_predicted"])
    ratio_comps = salary["salary"] / comps["comparable_median_salary"] if comps["comparable_median_salary"] > 0 else float("nan")

    if p_model < 25:
        band = "below"
    elif p_model <= 75:
        band = "inside"
    else:
        band = "above"
    model_reading = f"{band} the central range the model expects for these statistics"

    gap = war_pct - salary_pct
    if gap >= 20:
        rank_reading = "performance rank is well ahead of pay rank"
    elif gap <= -20:
        rank_reading = "pay rank is well ahead of performance rank"
    else:
        rank_reading = "performance rank and pay rank are close"

    if not np.isfinite(ratio_comps):
        comps_reading = "no comparable salary available"
    elif ratio_comps < 0.8:
        comps_reading = "paid less than that median"
    elif ratio_comps <= 1.25:
        comps_reading = "paid about the same as that median"
    else:
        comps_reading = "paid more than that median"

    if band == "below":
        headline = ("Performance supports the salary. The pay is below the central range the TabICLv2 model "
                    "expects for these statistics, so the player produces more than the salary implies.")
    elif band == "inside":
        headline = ("Performance supports the salary. The pay is inside the central range the TabICLv2 model "
                    "expects for these statistics.")
    else:
        headline = ("Performance does not support the salary. The pay is above the central range the TabICLv2 "
                    "model expects for these statistics.")
    if band == "above" and gap >= 20:
        headline = ("The evidence is mixed. The pay is above the central range the TabICLv2 model expects for "
                    "these statistics, but the performance rank is well ahead of the pay rank.")
    elif band != "above" and gap <= -20:
        headline = (f"The evidence is mixed. The pay is {band} the central range the TabICLv2 model expects for "
                    "these statistics, but the pay rank is well ahead of the performance rank.")

    return {
        "headline": headline,
        "model_band": band,
        "war_percentile": round(war_pct, 1),
        "mean_stat_percentile": round(mean_pct, 1),
        "salary_percentile": round(salary_pct, 1),
        "war_minus_salary_percentile": round(gap, 1),
        "model_reading": model_reading,
        "rank_reading": rank_reading,
        "comparables_reading": comps_reading,
        "salary_over_comparable_median": round(float(ratio_comps), 2) if np.isfinite(ratio_comps) else None,
    }


def analyze_player(df: pd.DataFrame, name: str, n_estimators: int | None = None) -> dict[str, Any]:
    """Full analysis bundle for one player. Returns plain Python types plus two DataFrames."""
    idx = player_index(df, name)
    perf = performance_table(df, idx)
    salary = salary_summary(df, idx)
    model = model_salary_for_player(df, idx, n_estimators)
    comps_table, comps = comparable_players(df, idx)
    cost = cost_per_war(df, idx)
    verdict = assessment(perf, salary, model, comps)
    return {
        "player": {
            "name": str(df.at[idx, NAME_COLUMN]),
            "team": str(df.at[idx, TEAM_COLUMN]),
            "index": idx,
        },
        "performance_table": perf,
        "comparables_table": comps_table,
        "salary": salary,
        "model": model,
        "comparables": comps,
        "cost_per_war": cost,
        "assessment": verdict,
        "stat_columns": list(STAT_COLUMNS),
        "lower_is_better": sorted(LOWER_IS_BETTER),
    }


def analysis_to_json(result: dict[str, Any]) -> dict[str, Any]:
    """Plain-Python view of an analysis: DataFrames become lists of records,
    numpy scalars become Python numbers, and non-finite floats become None."""
    return _sanitize(result)


def _sanitize(obj: Any) -> Any:
    if isinstance(obj, pd.DataFrame):
        return [_sanitize(row) for row in obj.to_dict(orient="records")]
    if isinstance(obj, dict):
        return {str(k): _sanitize(v) for k, v in obj.items()}
    if isinstance(obj, (list, tuple)):
        return [_sanitize(v) for v in obj]
    if isinstance(obj, np.ndarray):
        return [_sanitize(v) for v in obj.tolist()]
    if isinstance(obj, (bool, np.bool_)):
        return bool(obj)
    if isinstance(obj, (int, np.integer)):
        return int(obj)
    if isinstance(obj, (float, np.floating)):
        value = float(obj)
        return value if np.isfinite(value) else None
    return obj