feat: add 10 new ML models for auction optimization (Phases 1-6)

Phase 1 - Quick Wins:
- QuantileEnsemble: P10/P50/P90 predictions for risk-aware bidding
- MinutesSurvivalModel: Weibull AFT for minutes distribution modeling

Phase 2 - Adaptive Auction:
- BanditAuctionSolver: Thompson Sampling for live auction bids
- OpponentBidModel: Predict competitor bids via LightGBM
- BudgetOptimizer: Bayesian optimization for role-level allocation

Phase 3 - Deep Learning:
- RLAuctionPolicy: Double DQN agent for auction strategy
- SetTransformer: Team composition valuation via set-based ML

Phase 4 - Probabilistic:
- BayesianPlayerModel: Hierarchical pooling for rookie uncertainty
- ConformalPredictor: Calibrated prediction intervals

Phase 5 - Chemistry & Form:
- PlayerChemistryGAT: Graph attention network for player synergies
- PlayerFormModel: Hawkes process for form momentum

Phase 6 - Causal:
- TransferCausalModel: Causal forest for transfer effects
- AuctionEffectAnalyzer: Bid adjustment from causal analysis

81 tests passing

requirements.txt

@@ -26,6 +26,14 @@ torch-geometric>=2.5
 
 # Optimization
 pulp>=2.8
+scikit-optimize>=0.9
+
+# Bayesian modeling (optional)
+pymc>=5.0
+lifelines>=0.28
+
+# Causal inference (optional)
+econml>=0.15
 
 # Browser automation (Playwright fallback for Cloudflare)
 playwright>=1.42

src/models/__init__.py

@@ -1 +1,30 @@
 """Machine learning models."""
+
+from .base_model import BaseModel
+from .gbm_model import GBMEnsemble
+from .card_model import CardClassifier, PenaltyModel, GoalProbabilityModel
+from .tgcn_model import TemporalGNN
+from .distribution_head import (
+    SinhArcsinhDistribution,
+    BernoulliCleanSheet,
+    sinh_arcsinh_params,
+)
+from .train import ModelTrainer
+
+# Phase 1: Quantile regression + survival analysis
+from .quantile_model import QuantileEnsemble
+from .survival_model import MinutesSurvivalModel
+
+# Phase 4: Bayesian pooling + conformal prediction
+from .bayesian_pooling import BayesianPlayerModel
+from .conformal_predictor import ConformalPredictor
+
+# Phase 5: Player chemistry + form momentum
+from .gat_model import PlayerChemistryGAT
+from .hawkes_form import PlayerFormModel
+
+# Phase 3: Set-based team valuation
+from .set_transformer import SetTransformer
+
+# Phase 6: Causal inference for transfers
+from .causal_forest import TransferCausalModel, AuctionEffectAnalyzer

src/models/bayesian_pooling.py

@@ -0,0 +1,378 @@
+"""Hierarchical Bayesian partial pooling for player skill estimation.
+
+Implements two modes:
+- PyMC: Full MCMC-based hierarchical model with role-level priors
+- Scipy fallback: James-Stein-style shrinkage with empirical Bayes estimates
+"""
+
+import logging
+from typing import Optional, Tuple
+
+import numpy as np
+import pandas as pd
+
+from .base_model import BaseModel
+
+logger = logging.getLogger(__name__)
+
+VALID_ROLES = {"P", "D", "C", "A"}
+
+
+def _has_pymc() -> bool:
+    try:
+        import pymc as pm  # noqa: F401
+        return True
+    except ImportError:
+        return False
+
+
+class BayesianPlayerModel(BaseModel):
+    """Hierarchical Bayesian model with partial pooling by player role.
+
+    Two implementation modes:
+    - PyMC (if installed): Full MCMC hierarchical model
+    - Scipy fallback: Empirical Bayes with James-Stein shrinkage
+    """
+
+    def __init__(
+        self,
+        model_dir: str = "models_trained",
+        use_pymc: Optional[bool] = None,
+        samples: int = 2000,
+        tune: int = 1000,
+        chains: int = 2,
+        random_seed: int = 42,
+    ):
+        super().__init__(model_dir)
+        self.use_pymc = use_pymc if use_pymc is not None else _has_pymc()
+        self.samples = samples
+        self.tune = tune
+        self.chains = chains
+        self.random_seed = random_seed
+
+        self.trace = None
+        self.player_indices = {}
+        self.role_encoder = {}
+        self.role_reverse = {}
+        self.player_means = {}
+        self.player_vars = {}
+        self.role_means = {}
+        self.role_vars = {}
+        self.fitted = False
+        self._pymc_mode = False
+
+    def _validate_roles(self, X: pd.DataFrame):
+        if "role" not in X.columns:
+            raise ValueError("X must contain a 'role' column with values: P, D, C, A")
+        unknown = set(X["role"].unique()) - VALID_ROLES
+        if unknown:
+            raise ValueError(f"Unknown role values: {unknown}. Allowed: {VALID_ROLES}")
+
+    def _prepare_data(self, X: pd.DataFrame, y: pd.Series):
+        self._validate_roles(X)
+        roles = X["role"].values
+        unique_roles = sorted(VALID_ROLES)
+        self.role_encoder = {r: i for i, r in enumerate(unique_roles)}
+        self.role_reverse = {i: r for r, i in self.role_encoder.items()}
+
+        role_idx = np.array([self.role_encoder[r] for r in roles])
+        n_players = len(y)
+        n_roles = len(unique_roles)
+
+        player_map = {}
+        player_id = np.zeros(n_players, dtype=int)
+        for i in range(n_players):
+            key = (roles[i], i)
+            if key not in player_map:
+                player_map[key] = len(player_map)
+            player_id[i] = player_map[key]
+        n_unique = len(player_map)
+
+        return role_idx, player_id, n_players, n_roles, n_unique
+
+    def _fit_pymc(self, X: pd.DataFrame, y: pd.Series):
+        import pymc as pm
+
+        role_idx, player_id, n_players, n_roles, n_unique = self._prepare_data(X, y)
+
+        with pm.Model() as model:
+            mu_role = pm.Normal("mu_role", mu=6.0, sigma=2.0, shape=n_roles)
+            sigma_role = pm.HalfNormal("sigma_role", sigma=1.0, shape=n_roles)
+            sigma_obs = pm.HalfNormal("sigma_obs", sigma=1.0)
+
+            player_skill = pm.Normal(
+                "player_skill",
+                mu=mu_role[role_idx],
+                sigma=sigma_role[role_idx],
+                shape=n_players,
+            )
+
+            pm.Normal(
+                "observed_fv",
+                mu=player_skill,
+                sigma=sigma_obs,
+                observed=y.values,
+            )
+
+            self.trace = pm.sample(
+                draws=self.samples,
+                tune=self.tune,
+                chains=self.chains,
+                random_seed=self.random_seed,
+                progressbar=False,
+            )
+
+        logger.info(
+            f"PyMC model fitted: {n_players} players, {n_roles} roles, "
+            f"{len(self.trace.posterior.draw) * len(self.trace.posterior.chain)} posterior samples"
+        )
+        self._pymc_mode = True
+
+    def _fit_scipy(self, X: pd.DataFrame, y: pd.Series):
+        role_idx, player_id, n_players, n_roles, n_unique = self._prepare_data(X, y)
+        roles = X["role"].values
+
+        y_vals = y.values.astype(np.float64)
+
+        self.role_means = {}
+        self.role_vars = {}
+        for r_idx, r_name in self.role_reverse.items():
+            mask = role_idx == r_idx
+            if mask.sum() > 0:
+                self.role_means[r_name] = float(np.mean(y_vals[mask]))
+                role_var = float(np.var(y_vals[mask], ddof=1)) if mask.sum() > 1 else 0.0
+                self.role_vars[r_name] = role_var
+            else:
+                self.role_means[r_name] = 6.0
+                self.role_vars[r_name] = 2.0
+
+        player_data = {}
+        for i in range(n_players):
+            role = roles[i]
+            val = y_vals[i]
+            if role not in player_data:
+                player_data[role] = {}
+            player_data[role][i] = val
+
+        self.player_means = {}
+        self.player_vars = {}
+        for role, players_by_idx in player_data.items():
+            vals = list(players_by_idx.values())
+            role_mean = self.role_means[role]
+            role_var = max(self.role_vars[role], 1e-8)
+            for idx in players_by_idx:
+                self.player_means[idx] = vals[0]
+                self.player_vars[idx] = role_var
+
+        logger.info(
+            f"Scipy fallback fitted: {n_players} players, {n_roles} roles"
+        )
+        self._pymc_mode = False
+
+    def fit(self, X: pd.DataFrame, y: pd.Series, **kwargs):
+        if len(X) == 0:
+            raise ValueError("X cannot be empty")
+        if len(X) != len(y):
+            raise ValueError(f"X and y lengths must match: {len(X)} vs {len(y)}")
+
+        if self.use_pymc and _has_pymc():
+            self._fit_pymc(X, y)
+        else:
+            if self.use_pymc and not _has_pymc():
+                logger.warning("PyMC requested but not installed. Falling back to scipy.")
+            self.use_pymc = False
+            self._fit_scipy(X, y)
+
+        self.fitted = True
+        return self
+
+    def predict(self, X: pd.DataFrame) -> np.ndarray:
+        if not self.fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+        mean, _ = self.predict_with_uncertainty(X)
+        return mean
+
+    def predict_with_uncertainty(self, X: pd.DataFrame) -> Tuple[np.ndarray, np.ndarray]:
+        if not self.fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+        self._validate_roles(X)
+
+        if self._pymc_mode:
+            return self._predict_with_uncertainty_pymc(X)
+
+        means = np.zeros(len(X))
+        stds = np.zeros(len(X))
+        roles = X["role"].values
+
+        for i, role in enumerate(roles):
+            role_mean = self.role_means.get(role, 6.0)
+            role_var = self.role_vars.get(role, 2.0)
+            player_mean = self.player_means.get(i, role_mean)
+            player_var = self.player_vars.get(i, role_var)
+
+            shrinkage = role_var / max(role_var + player_var, 1e-8)
+            means[i] = role_mean + (1.0 - shrinkage) * (player_mean - role_mean)
+            stds[i] = np.sqrt(role_var * (1.0 - shrinkage))
+
+        return means, stds
+
+    def _predict_with_uncertainty_pymc(self, X: pd.DataFrame) -> Tuple[np.ndarray, np.ndarray]:
+        import pymc as pm
+        import arviz as az
+
+        n_players = len(X)
+        roles = X["role"].values
+
+        with pm.Model() as pred_model:
+            n_roles = len(self.role_encoder)
+            mu_role = pm.Normal("mu_role", mu=6.0, sigma=2.0, shape=n_roles)
+            sigma_role = pm.HalfNormal("sigma_role", sigma=1.0, shape=n_roles)
+            sigma_obs = pm.HalfNormal("sigma_obs", sigma=1.0)
+            player_skill = pm.Normal(
+                "player_skill",
+                mu=mu_role[[self.role_encoder.get(r, 0) for r in roles]],
+                sigma=sigma_role[[self.role_encoder.get(r, 0) for r in roles]],
+                shape=n_players,
+            )
+            pm.Normal("observed_fv", mu=player_skill, sigma=sigma_obs, shape=n_players)
+
+            ppc = pm.sample_posterior_predictive(
+                self.trace,
+                var_names=["observed_fv"],
+                random_seed=self.random_seed,
+                progressbar=False,
+            )
+
+        observed_samples = ppc.posterior_predictive["observed_fv"].values
+        draws_per_chain = observed_samples.shape[0]
+        n_chains = observed_samples.shape[1]
+        observed_flat = observed_samples.reshape(draws_per_chain * n_chains, n_players)
+
+        means = observed_flat.mean(axis=0)
+        stds = observed_flat.std(axis=0)
+
+        return means, stds
+
+    def posterior_predictive(self, X: pd.DataFrame, n_samples: int = 2000) -> np.ndarray:
+        if not self.fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+        self._validate_roles(X)
+
+        if self._pymc_mode:
+            return self._posterior_predictive_pymc(X, n_samples)
+
+        n_players = len(X)
+        means, stds = self.predict_with_uncertainty(X)
+        rng = np.random.RandomState(self.random_seed)
+        draws = rng.normal(
+            loc=means[np.newaxis, :],
+            scale=stds[np.newaxis, :] + 1e-6,
+            size=(n_samples, n_players),
+        )
+        return np.clip(draws, -10, 20)
+
+    def _posterior_predictive_pymc(self, X: pd.DataFrame, n_samples: int) -> np.ndarray:
+        import pymc as pm
+
+        n_players = len(X)
+        roles = X["role"].values
+
+        with pm.Model() as pred_model:
+            n_roles = len(self.role_encoder)
+            mu_role = pm.Normal("mu_role", mu=6.0, sigma=2.0, shape=n_roles)
+            sigma_role = pm.HalfNormal("sigma_role", sigma=1.0, shape=n_roles)
+            sigma_obs = pm.HalfNormal("sigma_obs", sigma=1.0)
+            player_skill = pm.Normal(
+                "player_skill",
+                mu=mu_role[[self.role_encoder.get(r, 0) for r in roles]],
+                sigma=sigma_role[[self.role_encoder.get(r, 0) for r in roles]],
+                shape=n_players,
+            )
+            pm.Normal("observed_fv", mu=player_skill, sigma=sigma_obs, shape=n_players)
+
+            ppc = pm.sample_posterior_predictive(
+                self.trace,
+                var_names=["observed_fv"],
+                random_seed=self.random_seed,
+                progressbar=False,
+            )
+
+        observed_samples = ppc.posterior_predictive["observed_fv"].values
+        draws_per_chain = observed_samples.shape[0]
+        n_chains = observed_samples.shape[1]
+        observed_flat = observed_samples.reshape(draws_per_chain * n_chains, n_players)
+
+        total = observed_flat.shape[0]
+        if total > n_samples:
+            rng = np.random.RandomState(self.random_seed)
+            idx = rng.choice(total, size=n_samples, replace=False)
+            return observed_flat[idx]
+        return observed_flat
+
+    def get_player_reliability(self, X: pd.DataFrame) -> np.ndarray:
+        if not self.fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+        self._validate_roles(X)
+
+        roles = X["role"].values
+        scores = np.zeros(len(X))
+
+        if self._pymc_mode:
+            means, stds = self.predict_with_uncertainty(X)
+            for i, role in enumerate(roles):
+                role_var = stds[i] ** 2
+                total_var = role_var + 1.0
+                scores[i] = np.clip(1.0 - (role_var / max(total_var, 1e-8)), 0.0, 1.0)
+            return np.clip(scores, 0.0, 1.0)
+
+        for i, role in enumerate(roles):
+            role_var = self.role_vars.get(role, 2.0)
+            player_var = self.player_vars.get(i, role_var)
+            total_var = role_var + player_var
+            scores[i] = np.clip(role_var / max(total_var, 1e-8), 0.0, 1.0)
+
+        return np.clip(scores, 0.0, 1.0)
+
+    def get_rookie_estimates(self, X: pd.DataFrame, min_observations: int = 5) -> pd.DataFrame:
+        if not self.fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+        self._validate_roles(X)
+
+        reliability = self.get_player_reliability(X)
+        rookie_mask = reliability < (1.0 / max(min_observations, 1))
+        means, stds = self.predict_with_uncertainty(X)
+        roles = X["role"].values
+
+        results = []
+        for i in range(len(X)):
+            if not rookie_mask[i]:
+                continue
+            role = roles[i]
+            role_mean = self.role_means.get(role, 6.0)
+            results.append({
+                "index": i,
+                "role": role,
+                "player_estimate": float(means[i]),
+                "role_mean": float(role_mean),
+                "naive_player_mean": self.player_means.get(i, role_mean),
+                "shrunken_estimate": float(means[i]),
+                "shrunken_std": float(stds[i]),
+                "reliability": float(reliability[i]),
+                "shrinkage_factor": float(
+                    (self.player_means.get(i, role_mean) - means[i])
+                    / max(abs(self.player_means.get(i, role_mean) - role_mean), 1e-8)
+                ) if abs(self.player_means.get(i, role_mean) - role_mean) > 1e-8 else 1.0,
+            })
+
+        if not results:
+            logger.info("No rookie players found (all have sufficient observations)")
+            return pd.DataFrame(columns=[
+                "index", "role", "player_estimate", "role_mean",
+                "naive_player_mean", "shrunken_estimate", "shrunken_std",
+                "reliability", "shrinkage_factor",
+            ])
+
+        df = pd.DataFrame(results)
+        df = df.sort_values("reliability")
+        logger.info(f"Found {len(results)} rookie players (heavy shrinkage toward role mean)")
+        return df

src/models/conformal_predictor.py

@@ -0,0 +1,184 @@
+"""Conformal prediction for calibrated prediction intervals.
+
+Wraps any sklearn-compatible predictor with distribution-free, finite-sample
+valid prediction bands via split conformal prediction with absolute residuals
+as the nonconformity score.
+"""
+
+import logging
+from typing import Optional, Tuple
+
+import numpy as np
+import pandas as pd
+
+logger = logging.getLogger(__name__)
+
+
+class ConformalPredictor:
+    """Split conformal prediction wrapper for calibrated uncertainty.
+
+    Computes nonconformity scores (absolute residuals) on a calibration set
+    and uses the (1-alpha) quantile to construct prediction bands with
+    guaranteed marginal coverage under exchangeability.
+    """
+
+    def __init__(self, base_model, alpha: float = 0.10):
+        if not 0 < alpha < 1:
+            raise ValueError(f"alpha must be in (0, 1), got {alpha}")
+        self.base_model = base_model
+        self.alpha = alpha
+        self.q_hat = None
+        self.is_calibrated = False
+        self._mad = False
+        self._cal_residuals = None
+        self._base_has_predict_std = self._check_predict_std()
+
+    def __repr__(self):
+        status = "calibrated" if self.is_calibrated else "uncalibrated"
+        mad_status = "MAD" if self._mad else "absolute"
+        return (
+            f"ConformalPredictor(base={type(self.base_model).__name__}, "
+            f"alpha={self.alpha:.2f}, status={status}, score={mad_status})"
+        )
+
+    def _check_predict_std(self) -> bool:
+        return hasattr(self.base_model, "predict_distribution") or hasattr(
+            self.base_model, "predict_std"
+        )
+
+    def _get_base_predictions(self, X: pd.DataFrame) -> np.ndarray:
+        preds = self.base_model.predict(X)
+        preds = np.asarray(preds, dtype=np.float64)
+        if preds.ndim == 2 and preds.shape[1] == 1:
+            preds = preds.ravel()
+        return preds
+
+    def _get_base_std(self, X: pd.DataFrame) -> Optional[np.ndarray]:
+        if hasattr(self.base_model, "predict_distribution"):
+            _, std = self.base_model.predict_distribution(X)
+            std = np.asarray(std, dtype=np.float64).ravel()
+            return std
+        if hasattr(self.base_model, "predict_std"):
+            std = self.base_model.predict_std(X)
+            std = np.asarray(std, dtype=np.float64).ravel()
+            return std
+        return None
+
+    def _compute_nonconformity(self, residuals: np.ndarray, std: Optional[np.ndarray] = None) -> np.ndarray:
+        if self._mad and std is not None and np.all(std > 0):
+            return np.abs(residuals) / np.maximum(std, 1e-8)
+        return np.abs(residuals)
+
+    def calibrate(self, X_cal: pd.DataFrame, y_cal: pd.Series, mad: bool = False):
+        if len(X_cal) == 0 or len(y_cal) == 0:
+            raise ValueError("Calibration set cannot be empty")
+
+        self._mad = mad
+        y_true = np.asarray(y_cal, dtype=np.float64).ravel()
+        y_pred = self._get_base_predictions(X_cal)
+        residuals = y_true - y_pred
+
+        if self._mad:
+            std = self._get_base_std(X_cal)
+            if std is None:
+                logger.warning(
+                    "MAD mode requested but base_model has no predict_std/predict_distribution. "
+                    "Falling back to absolute residuals."
+                )
+                self._mad = False
+                std = None
+        else:
+            std = None
+
+        self._cal_residuals = self._compute_nonconformity(residuals, std)
+        self._cal_preds = y_pred
+
+        n = len(self._cal_residuals)
+        correction = (1.0 + 1.0 / n)
+        q_idx = min(int(np.ceil((1.0 - self.alpha) * (n + 1))) - 1, n - 1)
+        if q_idx < 0:
+            q_idx = 0
+        self.q_hat = np.sort(self._cal_residuals)[q_idx]
+
+        self.is_calibrated = True
+        logger.info(
+            f"Conformal calibration complete: {n} samples, "
+            f"alpha={self.alpha:.2f}, quantile={self.q_hat:.4f}"
+        )
+        return self
+
+    def predict_with_band(self, X: pd.DataFrame) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+        if not self.is_calibrated:
+            raise RuntimeError("Model not calibrated. Call calibrate() first.")
+
+        y_pred = self._get_base_predictions(X)
+
+        if self._mad:
+            std = self._get_base_std(X)
+            if std is None:
+                self._mad = False
+                std = None
+        else:
+            std = None
+
+        if self._mad and std is not None:
+            half_width = self.q_hat * std
+        else:
+            half_width = self.q_hat
+
+        lower = y_pred - half_width
+        upper = y_pred + half_width
+
+        return y_pred, lower, upper
+
+    def predict(self, X: pd.DataFrame) -> np.ndarray:
+        return self._get_base_predictions(X)
+
+    def predict_interval_width(self, X: pd.DataFrame) -> float:
+        _, lower, upper = self.predict_with_band(X)
+        widths = upper - lower
+        return float(np.mean(widths))
+
+    def average_interval_width(self, X: pd.DataFrame) -> float:
+        widths = self.predict_interval_width(X)
+        return float(np.mean(widths))
+
+    def is_inside_band(self, X: pd.DataFrame, y_true: pd.Series) -> np.ndarray:
+        _, lower, upper = self.predict_with_band(X)
+        y = np.asarray(y_true, dtype=np.float64).ravel()
+        return (y >= lower) & (y <= upper)
+
+    def coverage(self, X: pd.DataFrame, y_true: pd.Series) -> float:
+        inside = self.is_inside_band(X, y_true)
+        return float(np.mean(inside))
+
+    def update(self, X_new: pd.DataFrame, y_new: pd.Series):
+        if not self.is_calibrated:
+            return self.calibrate(X_new, y_new)
+
+        y_true = np.asarray(y_new, dtype=np.float64).ravel()
+        y_pred = self._get_base_predictions(X_new)
+        residuals = y_true - y_pred
+
+        if self._mad:
+            std = self._get_base_std(X_new)
+            if std is not None:
+                new_scores = self._compute_nonconformity(residuals, std)
+            else:
+                new_scores = np.abs(residuals)
+        else:
+            new_scores = np.abs(residuals)
+
+        self._cal_residuals = np.concatenate([self._cal_residuals, new_scores])
+
+        n = len(self._cal_residuals)
+        q_idx = min(int(np.ceil((1.0 - self.alpha) * (n + 1))) - 1, n - 1)
+        if q_idx < 0:
+            q_idx = 0
+        self.q_hat = np.sort(self._cal_residuals)[q_idx]
+
+        logger.info(
+            f"Conformal update: +{len(new_scores)} samples, "
+            f"total={n}, alpha={self.alpha:.2f}, quantile={self.q_hat:.4f}"
+        )
+        return self

src/models/gat_model.py

@@ -0,0 +1,647 @@
+"""Graph Attention Network for player chemistry/interaction modeling.
+
+Models player-to-player synergy and redundancy using a GAT architecture when
+PyTorch Geometric is available, falling back to network-science feature
+extraction (PageRank, betweenness, clustering, eigenvector centrality) with a
+scikit-learn MLPRegressor otherwise.
+
+Replaces the stub tgcn_model.py.
+"""
+
+import logging
+from itertools import combinations
+from typing import Dict, List, Optional, Tuple
+
+import numpy as np
+import pandas as pd
+from sklearn.neural_network import MLPRegressor
+from sklearn.preprocessing import StandardScaler
+
+from .base_model import BaseModel
+
+logger = logging.getLogger(__name__)
+
+
+def _check_torch_geometric() -> bool:
+    try:
+        import torch                         # noqa: F401
+        from torch_geometric.nn import GATConv  # noqa: F401
+        return True
+    except ImportError:
+        return False
+
+
+def _build_graph_from_edges(edges: Dict[Tuple[str, str], float]) -> Tuple[Dict[str, int], np.ndarray, np.ndarray]:
+    """Convert edge dict to node mapping, adjacency matrix, and edge_index."""
+    nodes = set()
+    for (a, b) in edges:
+        nodes.add(a)
+        nodes.add(b)
+    node_list = sorted(nodes)
+    node_idx = {name: i for i, name in enumerate(node_list)}
+    n = len(node_list)
+    adj = np.zeros((n, n), dtype=np.float64)
+    for (a, b), w in edges.items():
+        i, j = node_idx[a], node_idx[b]
+        adj[i, j] = w
+    edge_list = [(node_idx[a], node_idx[b]) for (a, b) in edges]
+    edge_index = np.array(edge_list, dtype=np.int64).T if edge_list else np.empty((2, 0), dtype=np.int64)
+    return node_idx, adj, edge_index
+
+
+def _pagerank(adj: np.ndarray, alpha: float = 0.85, tol: float = 1e-6, max_iter: int = 100) -> np.ndarray:
+    n = adj.shape[0]
+    out_deg = adj.sum(axis=1)
+    out_deg[out_deg == 0] = 1.0
+    M = adj.T / out_deg
+    pr = np.ones(n) / n
+    for _ in range(max_iter):
+        pr_new = alpha * M.dot(pr) + (1.0 - alpha) / n
+        if np.abs(pr_new - pr).sum() < tol:
+            return pr_new
+        pr = pr_new
+    return pr
+
+
+def _eigenvector_centrality(adj: np.ndarray, tol: float = 1e-6, max_iter: int = 200) -> np.ndarray:
+    n = adj.shape[0]
+    x = np.ones(n)
+    for _ in range(max_iter):
+        x_new = adj.dot(x)
+        norm = np.linalg.norm(x_new, 2)
+        if norm < 1e-12:
+            break
+        x_new /= norm
+        if np.abs(x_new - x).sum() < tol:
+            return x_new
+        x = x_new
+    return x / np.linalg.norm(x, 2) if np.linalg.norm(x, 2) > 1e-12 else x
+
+
+def _betweenness_centrality(adj: np.ndarray) -> np.ndarray:
+    """Brandes algorithm for weighted undirected graphs."""
+    n = adj.shape[0]
+    if n <= 2:
+        return np.zeros(n)
+    dist = np.where(adj > 0, 1.0 / np.maximum(adj, 1e-12), np.inf)
+    np.fill_diagonal(dist, 0.0)
+    bc = np.zeros(n)
+    for s in range(n):
+        S = []
+        P = [[] for _ in range(n)]
+        sigma = np.zeros(n)
+        sigma[s] = 1.0
+        d = np.full(n, np.inf)
+        d[s] = 0.0
+        Q = [s]
+        for v in Q:
+            S.append(v)
+            for w in range(n):
+                if adj[v, w] <= 0 or w == v:
+                    continue
+                new_d = d[v] + dist[v, w]
+                if new_d < d[w] - 1e-12:
+                    d[w] = new_d
+                    sigma[w] = 0.0
+                    P[w] = [v]
+                    Q.append(w)
+                elif abs(new_d - d[w]) < 1e-12:
+                    sigma[w] += sigma[v]
+                    P[w].append(v)
+        delta = np.zeros(n)
+        while S:
+            w = S.pop()
+            for v in P[w]:
+                delta[v] += (sigma[v] / max(sigma[w], 1e-12)) * (1.0 + delta[w])
+            if w != s:
+                bc[w] += delta[w]
+    bc /= max(n - 1, 1) * max(n - 2, 1)
+    return bc
+
+
+def _clustering_coefficient(adj: np.ndarray) -> np.ndarray:
+    n = adj.shape[0]
+    cc = np.zeros(n)
+    binary = (adj > 0).astype(np.float64)
+    for i in range(n):
+        neighbors = np.where(binary[i] > 0)[0]
+        deg = len(neighbors)
+        if deg < 2:
+            cc[i] = 0.0
+            continue
+        sub = binary[np.ix_(neighbors, neighbors)]
+        triangles = np.sum(sub) / 2.0
+        cc[i] = 2.0 * triangles / (deg * (deg - 1))
+    return cc
+
+
+def _assortativity_by_position(adj: np.ndarray, role_groups: Dict[int, str]) -> float:
+    """Compute assortativity coefficient with respect to positional role."""
+    n = adj.shape[0]
+    if n <= 1:
+        return 0.0
+    binary = (adj > 0).astype(np.float64)
+    degrees = binary.sum(axis=1)
+    total_edges = degrees.sum()
+    if total_edges == 0:
+        return 0.0
+    same_type = 0.0
+    for i in range(n):
+        for j in range(n):
+            if binary[i, j] > 0 and role_groups.get(i) == role_groups.get(j):
+                same_type += 1.0
+    p_same = same_type / max(total_edges, 1e-12)
+    type_counts: Dict[str, float] = {}
+    for idx in range(n):
+        role = role_groups.get(idx, "unknown")
+        type_counts[role] = type_counts.get(role, 0.0) + degrees[idx]
+    total_deg = sum(type_counts.values())
+    expected = sum((t / max(total_deg, 1e-12)) ** 2 for t in type_counts.values())
+    max_possible = 1.0 - expected
+    if abs(max_possible) < 1e-12:
+        return 0.0
+    return (p_same - expected) / max_possible
+
+
+class PlayerChemistryGAT(BaseModel):
+    """Graph Attention Network for player chemistry/interaction modeling.
+
+    Two modes:
+    - PyTorch Geometric: full GATConv layers predicting delta-fantavoto.
+    - Fallback: scikit-learn MLPRegressor trained on network-science features
+      (PageRank, betweenness, clustering, eigenvector, assortativity).
+
+    Parameters
+    ----------
+    model_dir : str
+        Directory for persisting trained models.
+    hidden_dim : int
+        Hidden dimension for GAT / MLP.
+    num_layers : int
+        Number of GATConv layers.
+    heads : int
+        Number of attention heads.
+    dropout : float
+        Dropout probability.
+    """
+
+    def __init__(
+        self,
+        model_dir: str = "models_trained",
+        hidden_dim: int = 64,
+        num_layers: int = 2,
+        heads: int = 4,
+        dropout: float = 0.2,
+    ):
+        super().__init__(model_dir)
+        self.hidden_dim = hidden_dim
+        self.num_layers = num_layers
+        self.heads = heads
+        self.dropout = dropout
+
+        self._edges: Dict[Tuple[str, str], float] = {}
+        self.interaction_edges: Dict[Tuple[str, str], float] = {}
+        self._node_idx: Dict[str, int] = {}
+        self._idx_node: Dict[int, str] = {}
+        self._adj: Optional[np.ndarray] = None
+        self._edge_index: Optional[np.ndarray] = None
+
+        self._graph_features: Optional[np.ndarray] = None
+        self._feature_names: Optional[List[str]] = None
+        self._scaler = StandardScaler()
+
+        self._use_torch: bool = False
+        self._gat_module = None
+        self._mlp: Optional[MLPRegressor] = None
+
+    # ------------------------------------------------------------------
+    # Graph construction
+    # ------------------------------------------------------------------
+    def build_graph(self, historical_data: pd.DataFrame):
+        """Build player interaction graph from pass/assist/cross data.
+
+        Args:
+            historical_data: DataFrame with columns
+                [player, teammate, passes_to, assists_to, crosses_to, matchday]
+        """
+        required = {"player", "teammate", "passes_to", "assists_to", "crosses_to"}
+        missing = required - set(historical_data.columns)
+        if missing:
+            raise ValueError(f"Missing required columns: {missing}")
+
+        edges: Dict[Tuple[str, str], float] = {}
+        for (player, teammate), group in historical_data.groupby(["player", "teammate"]):
+            weight = (
+                group["passes_to"].sum()
+                + group["assists_to"].sum() * 5
+                + group["crosses_to"].sum() * 3
+            )
+            if weight > 0:
+                edges[(str(player), str(teammate))] = float(weight)
+
+        self._edges = edges
+        self.interaction_edges = dict(edges)
+        self._node_idx, self._adj, self._edge_index = _build_graph_from_edges(edges)
+        self._idx_node = {v: k for k, v in self._node_idx.items()}
+        logger.info(
+            "Built interaction graph: %d nodes, %d edges",
+            len(self._node_idx),
+            len(edges),
+        )
+        return self
+
+    # ------------------------------------------------------------------
+    # Graph feature extraction
+    # ------------------------------------------------------------------
+    def _extract_graph_features(self, role_map: Optional[Dict[str, str]] = None) -> np.ndarray:
+        n = len(self._node_idx)
+        if n == 0:
+            return np.empty((0, 6))
+
+        adj = self._adj.copy() if self._adj is not None else np.zeros((n, n))
+
+        pr = _pagerank(adj)
+        bc = _betweenness_centrality(adj)
+        cc = _clustering_coefficient(adj)
+        deg = adj.sum(axis=1) + adj.sum(axis=0)
+        ec = _eigenvector_centrality(adj)
+
+        role_groups: Dict[int, str] = {}
+        if role_map:
+            for i in range(n):
+                name = self._idx_node.get(i, "")
+                role_groups[i] = role_map.get(name, "unknown")
+        assortative = _assortativity_by_position(adj, role_groups) if role_map else 0.0
+
+        features = np.column_stack([pr, bc, cc, deg, ec, np.full(n, assortative)])
+        self._feature_names = [
+            "pagerank",
+            "betweenness",
+            "clustering_coef",
+            "degree",
+            "eigenvector",
+            "assortativity",
+        ]
+        return features
+
+    # ------------------------------------------------------------------
+    # PyTorch GAT module
+    # ------------------------------------------------------------------
+    def _init_torch_gat(self, in_channels: int):
+        try:
+            import torch
+            import torch.nn as nn
+            from torch_geometric.nn import GATConv
+
+            class GATModule(nn.Module):
+                def __init__(
+                    self,
+                    in_channels: int,
+                    hidden_dim: int,
+                    out_dim: int,
+                    heads: int,
+                    dropout: float,
+                ):
+                    super().__init__()
+                    self.conv1 = GATConv(in_channels, hidden_dim, heads=heads, dropout=dropout)
+                    self.conv2 = GATConv(hidden_dim * heads, out_dim, heads=1, concat=False, dropout=dropout)
+                    self.head = nn.Linear(out_dim, 1)
+
+                def forward(self, x, edge_index, edge_attr=None):
+                    x = self.conv1(x, edge_index)
+                    x = torch.relu(x)
+                    x = self.conv2(x, edge_index)
+                    x = torch.relu(x)
+                    return self.head(x).squeeze(-1)
+
+            self._gat_module = GATModule(
+                in_channels=in_channels,
+                hidden_dim=self.hidden_dim,
+                out_dim=self.hidden_dim // 2,
+                heads=self.heads,
+                dropout=self.dropout,
+            )
+            self._use_torch = True
+            logger.info("PyTorch Geometric GAT initialized (in=%d)", in_channels)
+        except ImportError:
+            self._use_torch = False
+            logger.info("torch_geometric not installed; using MLP fallback")
+
+    # ------------------------------------------------------------------
+    # Fit
+    # ------------------------------------------------------------------
+    def fit(self, X: pd.DataFrame, y: pd.Series, **kwargs):
+        """Fit the player chemistry model.
+
+        Args:
+            X: DataFrame. Must contain a 'player' column for node identification.
+               Optionally 'team' for team-based subgraph grouping and 'role' for
+               positional assortativity.
+            y: Target fantavoto scores per row.
+        """
+        if X.empty:
+            raise ValueError("X cannot be empty")
+
+        role_map: Optional[Dict[str, str]] = None
+        if "role" in X.columns:
+            role_map = dict(zip(X["player"].astype(str), X["role"].astype(str)))
+
+        team_groups = []
+        if "team" in X.columns:
+            for _, group in X.groupby("team"):
+                players = group["player"].unique().tolist()
+                if len(players) >= 2:
+                    for a, b in combinations(players, 2):
+                        team_groups.append({"player": str(a), "teammate": str(b),
+                                            "passes_to": 1, "assists_to": 0, "crosses_to": 0})
+
+        if "player" not in X.columns:
+            raise ValueError("X must contain a 'player' column")
+
+        if not self._edges and not team_groups:
+            if "player" in X.columns and "team" not in X.columns:
+                logger.warning("No interaction edges and no 'team' column; graph will be empty")
+            self._edges = {}
+            self.interaction_edges = {}
+            self._node_idx = {}
+            self._adj = np.array([[0.0]])
+            self._edge_index = np.empty((2, 0), dtype=np.int64)
+
+        if team_groups and not self._edges:
+            gdf = pd.DataFrame(team_groups)
+            if "passes_to" not in gdf.columns:
+                gdf["passes_to"] = 1
+            if "assists_to" not in gdf.columns:
+                gdf["assists_to"] = 0
+            if "crosses_to" not in gdf.columns:
+                gdf["crosses_to"] = 0
+            self.build_graph(gdf)
+
+        if "player" in X.columns:
+            player_set = set(self._node_idx.keys())
+            new_nodes = []
+            for row_player in X["player"]:
+                p = str(row_player)
+                if p and p not in player_set:
+                    idx = len(self._node_idx)
+                    self._node_idx[p] = idx
+                    self._idx_node[idx] = p
+                    player_set.add(p)
+                    new_nodes.append(p)
+            if self._adj is None:
+                n = len(self._node_idx)
+                self._adj = np.zeros((n, n))
+                self._edge_index = np.empty((2, 0), dtype=np.int64)
+            elif new_nodes:
+                old_n = self._adj.shape[0]
+                new_n = len(self._node_idx)
+                adj_new = np.zeros((new_n, new_n))
+                adj_new[:old_n, :old_n] = self._adj
+                self._adj = adj_new
+
+        self._graph_features = self._extract_graph_features(role_map)
+        n_nodes = len(self._node_idx)
+
+        if n_nodes == 0:
+            logger.warning("Zero nodes in graph; fitting dummy model")
+            self._mlp = MLPRegressor(
+                hidden_layer_sizes=(self.hidden_dim,),
+                max_iter=300,
+                random_state=42,
+            )
+            self._mlp.fit(np.zeros((1, 6)), np.zeros(1))
+            return self
+
+        node_to_player = {v: k for k, v in self._node_idx.items()}
+        target_deltas = np.zeros(n_nodes)
+        count_deltas = np.zeros(n_nodes)
+        player_indices = {
+            str(row["player"]): i
+            for i, (_, row) in enumerate(X.iterrows())
+        }
+        y_mean = float(np.mean(y)) if len(y) > 0 else 6.0
+
+        for i, (_, row) in enumerate(X.iterrows()):
+            p = str(row["player"])
+            if p in self._node_idx:
+                idx = self._node_idx[p]
+                target_deltas[idx] += (float(y.iloc[i]) - y_mean)
+                count_deltas[idx] += 1.0
+
+        for j in range(n_nodes):
+            if count_deltas[j] > 0:
+                target_deltas[j] /= count_deltas[j]
+
+        if _check_torch_geometric():
+            self._init_torch_gat(in_channels=self._graph_features.shape[1])
+            self._fit_torch_gat(self._graph_features, self._edge_index, target_deltas)
+        else:
+            self._fit_sklearn_fallback(self._graph_features, target_deltas)
+
+        logger.info(
+            "PlayerChemistryGAT fitted: %d nodes, %d edges, mode=%s",
+            n_nodes, len(self._edges),
+            "torch" if self._use_torch else "sklearn",
+        )
+        return self
+
+    def _fit_torch_gat(
+        self,
+        features: np.ndarray,
+        edge_index: np.ndarray,
+        targets: np.ndarray,
+    ):
+        import torch
+        import torch.nn as nn
+        import torch.optim as optim
+
+        device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+        x_tensor = torch.tensor(features, dtype=torch.float32).to(device)
+        edge_tensor = torch.tensor(edge_index, dtype=torch.long).to(device)
+        y_tensor = torch.tensor(targets, dtype=torch.float32).to(device)
+
+        self._gat_module = self._gat_module.to(device)
+        optimizer = optim.Adam(self._gat_module.parameters(), lr=0.01, weight_decay=1e-5)
+        criterion = nn.MSELoss()
+
+        n = features.shape[0]
+        self._gat_module.train()
+        for epoch in range(200):
+            optimizer.zero_grad()
+            pred = self._gat_module(x_tensor, edge_tensor)
+            loss = criterion(pred, y_tensor)
+            loss.backward()
+            optimizer.step()
+            if (epoch + 1) % 50 == 0:
+                logger.debug(f"GAT epoch {epoch + 1}: loss={loss.item():.6f}")
+
+        self._gat_module.eval()
+
+    def _fit_sklearn_fallback(self, features: np.ndarray, targets: np.ndarray):
+        self._scaler.fit(features)
+        features_scaled = self._scaler.transform(features)
+        n_samples = features_scaled.shape[0]
+        val_frac = 0.1 if n_samples >= 20 else 0.0
+        self._mlp = MLPRegressor(
+            hidden_layer_sizes=(self.hidden_dim, self.hidden_dim // 2),
+            activation="relu",
+            solver="adam",
+            max_iter=500,
+            random_state=42,
+            early_stopping=(n_samples >= 20),
+            validation_fraction=val_frac if val_frac > 0 else 0.1,
+            n_iter_no_change=20,
+        )
+        self._mlp.fit(features_scaled, targets)
+        self._use_torch = False
+
+    # ------------------------------------------------------------------
+    # Predict
+    # ------------------------------------------------------------------
+    def predict(self, X: pd.DataFrame) -> np.ndarray:
+        """Return chemistry-adjusted point projection bonuses.
+
+        These bonuses can be positive (synergy) or negative (redundancy).
+
+        Returns an array the same length as X with delta-fantavoto values.
+        """
+        if self._adj is None or len(self._node_idx) == 0:
+            return np.zeros(X.shape[0])
+
+        role_map: Optional[Dict[str, str]] = None
+        if "role" in X.columns:
+            role_map = dict(zip(X["player"].astype(str), X["role"].astype(str)))
+
+        features = self._extract_graph_features(role_map)
+        n_nodes = len(self._node_idx)
+        if features.shape[0] != n_nodes or n_nodes == 0:
+            return np.zeros(X.shape[0])
+
+        if self._use_torch and self._gat_module is not None:
+            return self._predict_torch(features, X)
+        elif self._mlp is not None:
+            return self._predict_sklearn(features, X)
+        return np.zeros(X.shape[0])
+
+    def _predict_torch(self, features: np.ndarray, X: pd.DataFrame) -> np.ndarray:
+        import torch
+
+        device = next(self._gat_module.parameters()).device
+        x_tensor = torch.tensor(features, dtype=torch.float32).to(device)
+        edge_tensor = torch.tensor(self._edge_index, dtype=torch.long).to(device)
+
+        self._gat_module.eval()
+        with torch.no_grad():
+            node_deltas = self._gat_module(x_tensor, edge_tensor).cpu().numpy()
+
+        result = np.zeros(X.shape[0])
+        for i, (_, row) in enumerate(X.iterrows()):
+            p = str(row.get("player", row.get("name", "")))
+            if p in self._node_idx:
+                result[i] = float(node_deltas[self._node_idx[p]])
+        return result
+
+    def _predict_sklearn(self, features: np.ndarray, X: pd.DataFrame) -> np.ndarray:
+        features_scaled = self._scaler.transform(features)
+        node_deltas = self._mlp.predict(features_scaled)
+
+        result = np.zeros(X.shape[0])
+        for i, (_, row) in enumerate(X.iterrows()):
+            p = str(row.get("player", row.get("name", "")))
+            if p in self._node_idx:
+                result[i] = float(node_deltas[self._node_idx[p]])
+        return result
+
+    # ------------------------------------------------------------------
+    # Interaction feature extraction (public API)
+    # ------------------------------------------------------------------
+    def extract_interaction_features(self, player: str, teammates: List[str]) -> Dict[str, float]:
+        """Extract interaction features between a player and their teammates.
+
+        Args:
+            player: Player name.
+            teammates: List of teammate names.
+
+        Returns:
+            Dict with keys: interaction_outgoing_sum, interaction_incoming_sum,
+            interaction_synergy.
+        """
+        outgoing = 0.0
+        incoming = 0.0
+        synergy = 0.0
+        count = 0
+        for teammate in teammates:
+            w = self._edges.get((str(player), str(teammate)), 0.0)
+            outgoing += w
+            w_in = self._edges.get((str(teammate), str(player)), 0.0)
+            incoming += w_in
+            if w > 0 or w_in > 0:
+                synergy += (w + w_in) / 2.0
+                count += 1
+        avg_synergy = synergy / max(count, 1)
+        return {
+            "interaction_outgoing_sum": float(outgoing),
+            "interaction_incoming_sum": float(incoming),
+            "interaction_synergy": float(avg_synergy),
+        }
+
+    # ------------------------------------------------------------------
+    # Interaction bonus / chemistry matrix
+    # ------------------------------------------------------------------
+    def compute_interaction_bonus(self, player_a: str, player_b: str) -> float:
+        """Bonus factor for two players in the same lineup.
+
+        Positive = synergy, negative = redundancy.
+        """
+        if not self._edges:
+            return 0.0
+        a = str(player_a)
+        b = str(player_b)
+        weight = self._edges.get((a, b), 0.0) + self._edges.get((b, a), 0.0)
+        max_weight = max(self._edges.values()) if self._edges else 1.0
+        return float(weight) / max(max_weight, 1.0) if weight > 0 else 0.0
+
+    def get_chemistry_matrix(self, players: List[str]) -> np.ndarray:
+        """Return NxN matrix of pairwise chemistry bonuses for given players.
+
+        Positive values = synergy, negative = redundancy.
+        """
+        n = len(players)
+        matrix = np.zeros((n, n))
+        for i in range(n):
+            for j in range(n):
+                if i != j:
+                    matrix[i, j] = self.compute_interaction_bonus(players[i], players[j])
+        return matrix
+
+    def get_redundancy_penalty(self, players: List[str]) -> Dict[Tuple[str, str], float]:
+        """Identify negative synergies between players on the same team.
+
+        Two players may compete for the same actions (both take corners,
+        both demand the ball in similar zones). Returns dict mapping
+        player-pairs to a negative penalty value (0 = no redundancy).
+
+        The penalty is derived from co-occurrence overlap: if two players
+        both have high outgoing edges to the same teammates and both
+        receive from the same sources, they're likely redundant.
+        """
+        penalty: Dict[Tuple[str, str], float] = {}
+        if not self._edges:
+            return penalty
+
+        for i, a in enumerate(players):
+            for b in players[i + 1:]:
+                a_out = set(t for (p, t) in self._edges if p == str(a))
+                b_out = set(t for (p, t) in self._edges if p == str(b))
+                a_in = set(p for (p, t) in self._edges if t == str(a))
+                b_in = set(p for (p, t) in self._edges if t == str(b))
+
+                out_overlap = len(a_out & b_out)
+                in_overlap = len(a_in & b_in)
+                total = max(len(a_out | b_out) + len(a_in | b_in), 1)
+                overlap_ratio = (out_overlap + in_overlap) / total
+
+                if overlap_ratio > 0.3:
+                    penalty[(str(a), str(b))] = -overlap_ratio
+
+        return penalty

src/models/hawkes_form.py

@@ -0,0 +1,503 @@
+"""Hawkes-process player form model for momentum modeling.
+
+Models player form as a self-exciting Hawkes process: good performances
+increase the probability of more good performances (momentum / hot streak).
+
+Uses scipy.optimize.minimize to fit per-player Hawkes parameters
+(mu, alpha, beta, s) via maximum likelihood.
+"""
+
+import logging
+from dataclasses import dataclass
+from typing import Dict, List, Optional, Tuple
+
+import numpy as np
+import pandas as pd
+from scipy.optimize import minimize
+
+from .base_model import BaseModel
+
+logger = logging.getLogger(__name__)
+
+FORM_STATUS_HOT = "HOT"
+FORM_STATUS_COLD = "COLD"
+FORM_STATUS_NEUTRAL = "NEUTRAL"
+FORM_STATUSES = {FORM_STATUS_HOT, FORM_STATUS_COLD, FORM_STATUS_NEUTRAL}
+
+_EPS = 1e-12
+_MIN_BETA = 1e-4
+_MAX_ALPHA = 20.0
+_MAX_MU = 50.0
+_MIN_S = -3.0
+_MAX_S = 3.0
+
+
+@dataclass
+class PlayerHawkesParams:
+    mu: float
+    alpha: float
+    beta: float
+    s: float
+    baseline_mean: float
+    baseline_std: float
+
+
+class PlayerFormModel(BaseModel):
+    """Self-exciting Hawkes process for player form / momentum.
+
+    Each player's match performances are modeled as a point process where
+    above-average games ("excitatory events") temporarily raise the
+    probability of subsequent above-average games.
+
+    Parameters
+    ----------
+    model_dir : str
+        Directory for persisting trained models.
+    decay_window : int
+        Maximum match-gap over which excitation persists (default 10).
+    """
+
+    def __init__(
+        self,
+        model_dir: str = "models_trained",
+        decay_window: int = 10,
+    ):
+        super().__init__(model_dir)
+        self.decay_window = int(decay_window)
+        self._player_params: Dict[str, PlayerHawkesParams] = {}
+        self._global_baseline: float = 6.0
+        self._global_std: float = 1.5
+
+    # ------------------------------------------------------------------
+    # Hawkes log-likelihood and intensity
+    # ------------------------------------------------------------------
+    @staticmethod
+    def _hawkes_intensity(times: np.ndarray, event_mask: np.ndarray, mu: float,
+                          alpha: float, beta: float) -> np.ndarray:
+        lam = np.full_like(times, mu, dtype=np.float64)
+        for i in range(1, len(times)):
+            if event_mask[i - 1]:
+                dt = times[i:] - times[i - 1]
+                mask = dt > 0
+                lam[i:] += alpha * np.exp(-beta * dt) * mask
+        return np.maximum(lam, _EPS)
+
+    @staticmethod
+    def _hawkes_integral(mu: float, alpha: float, beta: float,
+                         event_times: np.ndarray, T: float) -> float:
+        result = mu * T
+        for te in event_times:
+            remaining = T - te
+            if remaining > 0:
+                result += (alpha / beta) * (1.0 - np.exp(-beta * remaining))
+        return result
+
+    @staticmethod
+    def _hawkes_nll(params: np.ndarray, times: np.ndarray, event_mask: np.ndarray) -> float:
+        mu, alpha, beta = max(params[0], _EPS), max(params[1], _EPS), max(params[2], _MIN_BETA)
+        T = times[-1] if len(times) > 0 else 1.0
+        lam = PlayerFormModel._hawkes_intensity(times, event_mask, mu, alpha, beta)
+        log_lik = np.sum(np.log(lam))
+        integral = PlayerFormModel._hawkes_integral(mu, alpha, beta,
+                                                     times[event_mask.astype(bool)], T)
+        return -(log_lik - integral)
+
+    # ------------------------------------------------------------------
+    # Fit helper: tune per-player
+    # ------------------------------------------------------------------
+    def _fit_player(self, times: np.ndarray, scores: np.ndarray) -> Optional[PlayerHawkesParams]:
+        n = len(scores)
+        if n < 5:
+            return None
+
+        times_float = times.astype(np.float64)
+        scores_float = scores.astype(np.float64)
+        baseline_mean = float(np.mean(scores_float))
+        baseline_std = float(np.std(scores_float, ddof=1)) if n > 1 else 1.0
+
+        best_nll = float("inf")
+        best_params = None
+
+        for s_candidate in [-1.0, 0.0, 0.5, 1.0, 1.5]:
+            threshold = baseline_mean + s_candidate * max(baseline_std, 0.5)
+            event_mask = (scores_float > threshold).astype(np.float64)
+            n_events = event_mask.sum()
+            if n_events < 2:
+                continue
+
+            init_mu = max(max(n_events / max(times_float[-1] - times_float[0], 1.0), 0.05), _EPS)
+            init_alpha = min(n_events / max(n, 1) * 2.0, _MAX_ALPHA)
+            init_beta = 0.5
+
+            for init_scale in [0.5, 1.0, 2.0]:
+                x0 = np.array([
+                    init_mu * init_scale,
+                    init_alpha * init_scale,
+                    init_beta * init_scale,
+                ])
+
+                try:
+                    result = minimize(
+                        self._hawkes_nll,
+                        x0,
+                        args=(times_float, event_mask),
+                        method="L-BFGS-B",
+                        bounds=[(_EPS, _MAX_MU), (_EPS, _MAX_ALPHA), (_MIN_BETA, 10.0)],
+                        options={"maxiter": 500, "ftol": 1e-10},
+                    )
+                    if result.success and result.fun < best_nll:
+                        best_nll = result.fun
+                        best_params = PlayerHawkesParams(
+                            mu=float(max(result.x[0], _EPS)),
+                            alpha=float(max(result.x[1], _EPS)),
+                            beta=float(max(result.x[2], _MIN_BETA)),
+                            s=float(s_candidate),
+                            baseline_mean=baseline_mean,
+                            baseline_std=baseline_std,
+                        )
+                except Exception:
+                    continue
+
+        if best_params is None:
+            best_params = PlayerHawkesParams(
+                mu=0.1,
+                alpha=1.0,
+                beta=0.3,
+                s=0.0,
+                baseline_mean=baseline_mean,
+                baseline_std=baseline_std,
+            )
+
+        return best_params
+
+    # ------------------------------------------------------------------
+    # Fit
+    # ------------------------------------------------------------------
+    def fit(self, X: pd.DataFrame, y: pd.Series, **kwargs):
+        """Fit per-player Hawkes parameters.
+
+        X must contain:
+          - 'player' or 'name': player identifier.
+          - 'match_date' or 'matchday': temporal ordering column.
+
+        y: fantavoto scores.
+
+        Additional kwargs:
+          - 'match_date' column name override.
+        """
+        if X.empty:
+            raise ValueError("X cannot be empty")
+
+        player_col = None
+        for candidate in ["player", "name"]:
+            if candidate in X.columns:
+                player_col = candidate
+                break
+        if player_col is None:
+            raise ValueError("X must contain a 'player' or 'name' column")
+
+        date_col = kwargs.get("date_col", None)
+        if date_col is None:
+            for candidate in ["match_date", "matchday", "date", "giornata"]:
+                if candidate in X.columns:
+                    date_col = candidate
+                    break
+        if date_col is None:
+            logger.warning("No date column found; using row index as temporal order")
+            times = np.arange(len(X), dtype=np.float64)
+        else:
+            col_vals = X[date_col]
+            if pd.api.types.is_datetime64_any_dtype(col_vals):
+                times = col_vals.astype(np.int64).values.astype(np.float64) / 1e9 / 86400.0
+            else:
+                times = col_vals.astype(np.float64).values
+
+        players = X[player_col].astype(str).values
+        scores = y.values.astype(np.float64)
+
+        self._global_baseline = float(np.mean(scores)) if len(scores) > 0 else 6.0
+        self._global_std = float(np.std(scores, ddof=1)) if len(scores) > 1 else 1.5
+
+        self._player_params = {}
+        unique_players = np.unique(players)
+        fitted = 0
+        for player in unique_players:
+            mask = players == player
+            p_times = times[mask]
+            p_scores = scores[mask]
+            sort_idx = np.argsort(p_times)
+            p_times = p_times[sort_idx]
+            p_scores = p_scores[sort_idx]
+            params = self._fit_player(p_times, p_scores)
+            if params is not None:
+                self._player_params[str(player)] = params
+                fitted += 1
+
+        logger.info(
+            "Hawkes form model fitted: %d/%d players with sufficient history",
+            fitted, len(unique_players),
+        )
+        return self
+
+    # ------------------------------------------------------------------
+    # Predict
+    # ------------------------------------------------------------------
+    def predict(self, X: pd.DataFrame) -> np.ndarray:
+        """Return form-adjusted projections as additive bonuses to base.
+
+        Positive = player is in form (HOT), negative = out of form (COLD).
+        """
+        if not self._player_params:
+            return np.zeros(X.shape[0])
+
+        player_col = "player" if "player" in X.columns else "name"
+        multipliers = self._compute_multipliers(X)
+        base = X.get("base_prediction", pd.Series(np.full(X.shape[0], 6.0)))
+        base_vals = base.values.astype(np.float64)
+        return (multipliers - 1.0) * base_vals
+
+    def _compute_multipliers(self, X: pd.DataFrame) -> np.ndarray:
+        player_col = "player" if "player" in X.columns else "name"
+        multipliers = np.ones(X.shape[0], dtype=np.float64)
+        players = X[player_col].astype(str).values
+
+        for i, player in enumerate(players):
+            params = self._player_params.get(player)
+            if params is None:
+                continue
+            recent = self._compute_recent_intensity(params)
+            base_rate = params.mu
+            if base_rate > _EPS:
+                ratio = recent / base_rate
+                clamped = np.clip(ratio, 0.85, 1.15)
+                multipliers[i] = float(clamped)
+
+        return multipliers
+
+    def _compute_recent_intensity(self, params: PlayerHawkesParams) -> float:
+        return max(params.mu, _EPS)
+
+    # ------------------------------------------------------------------
+    # Momentum projection
+    # ------------------------------------------------------------------
+    def predict_momentum(
+        self,
+        X: pd.DataFrame,
+        player_history: pd.DataFrame,
+        n_future: int = 5,
+    ) -> np.ndarray:
+        """Project form trajectory for next *n_future* matches.
+
+        Returns (n_future, n_players) array of momentum multipliers.
+        """
+        if not self._player_params:
+            return np.ones((n_future, X.shape[0]))
+
+        player_col = "player" if "player" in X.columns else "name"
+        players = X[player_col].astype(str).values
+        n_players = X.shape[0]
+        trajectory = np.ones((n_future, n_players), dtype=np.float64)
+
+        history_player_col = None
+        for c in ["player", "name"]:
+            if c in player_history.columns:
+                history_player_col = c
+                break
+
+        history_date_col = None
+        for c in ["match_date", "matchday", "date"]:
+            if c in player_history.columns:
+                history_date_col = c
+                break
+
+        for j, player in enumerate(players):
+            params = self._player_params.get(player)
+            if params is None:
+                continue
+
+            event_times = []
+            if history_player_col and history_date_col:
+                p_hist = player_history[player_history[history_player_col].astype(str) == player]
+                if len(p_hist) > 0:
+                    target_col = None
+                    for c in ["fantavoto", "score", "fv"]:
+                        if c in p_hist.columns:
+                            target_col = c
+                            break
+                    if target_col and params.baseline_std > 0:
+                        threshold = params.baseline_mean + params.s * params.baseline_std
+                        p_sorted = p_hist.sort_values(history_date_col)
+                        times = p_sorted[history_date_col].values
+                        scores = p_sorted[target_col].values
+                        if pd.api.types.is_datetime64_any_dtype(p_sorted[history_date_col]):
+                            event_times_float = times.astype(np.int64).astype(np.float64) / 1e9 / 86400.0
+                        else:
+                            event_times_float = times.astype(np.float64)
+                        for ti, si in zip(event_times_float, scores):
+                            if float(si) > threshold:
+                                event_times.append(ti)
+
+            if not event_times:
+                continue
+
+            last_t = max(event_times)
+            for k in range(1, n_future + 1):
+                future_t = last_t + k
+                lam = params.mu
+                for te in event_times:
+                    dt = future_t - te
+                    if dt > 0 and dt <= self.decay_window:
+                        lam += params.alpha * np.exp(-params.beta * dt)
+                trajectory[k - 1, j] = float(np.clip(lam / max(params.mu, _EPS), 0.85, 1.15))
+
+        return trajectory
+
+    # ------------------------------------------------------------------
+    # Form status
+    # ------------------------------------------------------------------
+    def get_form_status(self, X: pd.DataFrame) -> List[str]:
+        """Return status string per player: HOT, COLD, or NEUTRAL."""
+        player_col = "player" if "player" in X.columns else "name"
+        players = X[player_col].astype(str).values
+        statuses: List[str] = []
+
+        for player in players:
+            params = self._player_params.get(player)
+            if params is None:
+                statuses.append(FORM_STATUS_NEUTRAL)
+                continue
+            intensity = self._compute_recent_intensity(params)
+            baseline = max(params.mu, _EPS)
+            ratio = intensity / baseline
+            if ratio > 1.1 and params.alpha > 0.1:
+                statuses.append(FORM_STATUS_HOT)
+            elif ratio < 0.9:
+                statuses.append(FORM_STATUS_COLD)
+            else:
+                statuses.append(FORM_STATUS_NEUTRAL)
+
+        return statuses
+
+    # ------------------------------------------------------------------
+    # Intensity curve
+    # ------------------------------------------------------------------
+    def compute_intensity_curve(
+        self,
+        player_name: str,
+        history: pd.DataFrame,
+        match_dates: np.ndarray,
+        future_dates: np.ndarray,
+    ) -> np.ndarray:
+        """Compute λ(t) over match_dates and future_dates for one player.
+
+        Returns an array of intensity values at each date.
+        """
+        params = self._player_params.get(str(player_name))
+        if params is None:
+            mu = self._global_baseline
+            all_dates = np.concatenate([match_dates, future_dates])
+            return np.full_like(all_dates, max(mu, _EPS), dtype=np.float64)
+
+        target_col = None
+        for c in ["fantavoto", "score", "fv"]:
+            if c in history.columns:
+                target_col = c
+                break
+
+        threshold = params.baseline_mean + params.s * max(params.baseline_std, 0.5)
+        event_times = []
+        if target_col:
+            for _, row in history.iterrows():
+                if float(row.get(target_col, 0)) > threshold:
+                    event_times.append(float(row.name) if isinstance(row.name, (int, float)) else 0.0)
+
+        if match_dates is not None and len(match_dates) > 0:
+            match_vals = match_dates.astype(np.float64)
+            for ti in match_vals:
+                if ti not in event_times:
+                    score = None
+                    for _, row in history.iterrows():
+                        d_val = float(row.name) if isinstance(row.name, (int, float)) else 0.0
+                        if abs(d_val - ti) < _EPS:
+                            score = row.get(target_col, 0) if target_col else 0
+                            break
+                    if score is not None and float(score) > threshold:
+                        event_times.append(ti)
+
+        all_dates = np.concatenate([
+            match_dates.astype(np.float64) if match_dates is not None and len(match_dates) > 0
+            else np.array([], dtype=np.float64),
+            future_dates.astype(np.float64) if future_dates is not None and len(future_dates) > 0
+            else np.array([], dtype=np.float64),
+        ])
+
+        if len(all_dates) == 0:
+            return np.array([params.mu])
+
+        intensity = np.full(len(all_dates), params.mu, dtype=np.float64)
+        for i, t in enumerate(all_dates):
+            lam = params.mu
+            for te in event_times:
+                dt = t - te
+                if dt > 0 and dt <= self.decay_window:
+                    lam += params.alpha * np.exp(-params.beta * dt)
+                elif dt > self.decay_window:
+                    pass
+            intensity[i] = max(lam, _EPS)
+
+        return intensity
+
+    # ------------------------------------------------------------------
+    # Streak detection
+    # ------------------------------------------------------------------
+    def detect_streak(
+        self,
+        player_history: pd.DataFrame,
+    ) -> Tuple[bool, int, str]:
+        """Detect whether a player is on a hot or cold streak.
+
+        Returns (is_streak, streak_length, streak_direction).
+        streak_direction is "HOT_STREAK" or "COLD_STREAK".
+        """
+        if len(player_history) < 3:
+            return (False, 0, "NO_STREAK")
+
+        target_col = None
+        for c in ["fantavoto", "score", "fv"]:
+            if c in player_history.columns:
+                target_col = c
+                break
+        if target_col is None:
+            return (False, 0, "NO_STREAK")
+
+        date_col = None
+        for c in ["match_date", "matchday", "date"]:
+            if c in player_history.columns:
+                date_col = c
+                break
+        if date_col:
+            sorted_hist = player_history.sort_values(date_col)
+        else:
+            sorted_hist = player_history
+
+        scores = sorted_hist[target_col].values.astype(np.float64)
+        mean_score = np.mean(scores)
+        std_score = max(np.std(scores, ddof=1), 0.5)
+
+        above = scores[-1] > mean_score + 0.5 * std_score
+        below = scores[-1] < mean_score - 0.5 * std_score
+
+        if not above and not below:
+            return (False, 0, "NO_STREAK")
+
+        direction = "HOT_STREAK" if above else "COLD_STREAK"
+        streak_len = 1
+        for j in range(len(scores) - 2, -1, -1):
+            if direction == "HOT_STREAK" and scores[j] > mean_score + 0.5 * std_score:
+                streak_len += 1
+            elif direction == "COLD_STREAK" and scores[j] < mean_score - 0.5 * std_score:
+                streak_len += 1
+            else:
+                break
+
+        return (streak_len >= 3, streak_len, direction)

src/models/quantile_model.py

@@ -0,0 +1,207 @@
+"""Quantile regression ensemble for probabilistic score prediction.
+
+Trains one LightGBM quantile regressor per target quantile (default P10, P50, P90)
+to output a full predictive distribution of Fantavoto scores. Supports downside risk,
+upside potential, and Value-at-Risk-safe estimates.
+"""
+
+import logging
+from typing import Optional, Tuple
+
+import numpy as np
+import pandas as pd
+from sklearn.preprocessing import StandardScaler
+
+from .base_model import BaseModel
+
+logger = logging.getLogger(__name__)
+
+
+class QuantileEnsemble(BaseModel):
+    """Quantile regression ensemble for multi-quantile score prediction.
+
+    Stores one LGBMRegressor per quantile, each with ``objective="quantile"``
+    and the corresponding ``alpha`` value.  Provides convenience accessors for
+    point predictions (P50), downside/upside risk, and VaR-safe floors.
+    """
+
+    def __init__(
+        self,
+        quantiles: Tuple[float, ...] = (0.10, 0.50, 0.90),
+        model_dir: str = "models_trained",
+        n_estimators: int = 300,
+        learning_rate: float = 0.05,
+        max_depth: int = 5,
+        num_leaves: int = 31,
+        subsample: float = 0.8,
+        colsample_bytree: float = 0.8,
+    ):
+        super().__init__(model_dir)
+        self.quantiles = tuple(quantiles)
+        self.n_estimators = n_estimators
+        self.learning_rate = learning_rate
+        self.max_depth = max_depth
+        self.num_leaves = num_leaves
+        self.subsample = subsample
+        self.colsample_bytree = colsample_bytree
+
+        self._models: dict = {}               # alpha → LGBMRegressor
+        self.feature_names = None
+        self.scaler = StandardScaler()
+
+    # ------------------------------------------------------------------
+    # Training
+    # ------------------------------------------------------------------
+    def fit(self, X: pd.DataFrame, y: pd.Series, **kwargs):
+        """Fit one LightGBM quantile regressor per target quantile.
+
+        Args:
+            X: Feature matrix.
+            y: Target values (fantavoto scores).
+        """
+        try:
+            import lightgbm as lgb
+        except ImportError:
+            raise ImportError("LightGBM is required.  pip install lightgbm")
+
+        self.feature_names = list(X.columns)
+        X_clean = X.select_dtypes(include=[np.number]).fillna(0)
+        X_scaled = self.scaler.fit_transform(X_clean)
+
+        self._models = {}
+        for alpha in self.quantiles:
+            model = lgb.LGBMRegressor(
+                n_estimators=self.n_estimators,
+                learning_rate=self.learning_rate,
+                max_depth=self.max_depth,
+                num_leaves=self.num_leaves,
+                subsample=self.subsample,
+                colsample_bytree=self.colsample_bytree,
+                objective="quantile",
+                alpha=alpha,
+                random_state=42,
+                verbose=-1,
+            )
+            model.fit(X_scaled, y)
+            label = f"P{int(alpha * 100)}"
+            self._models[label] = model
+            logger.info(f"Quantile model {label} (α={alpha:.2f}) trained")
+
+        logger.info(f"QuantileEnsemble fitted: {list(self._models.keys())}")
+        return self
+
+    # ------------------------------------------------------------------
+    # Prediction helpers
+    # ------------------------------------------------------------------
+    def _preprocess(self, X: pd.DataFrame) -> np.ndarray:
+        """Clean, impute, and scale features."""
+        if self.feature_names is None:
+            raise RuntimeError("Model not trained. Call fit() first.")
+        X_c = X[self.feature_names].select_dtypes(include=[np.number]).fillna(0)
+        return self.scaler.transform(X_c)
+
+    # ------------------------------------------------------------------
+    # Core predict
+    # ------------------------------------------------------------------
+    def predict(self, X: pd.DataFrame) -> dict:
+        """Return a dict mapping quantile label → prediction array.
+
+        Example:
+            {"P10": array([4.2, 5.1, ...]), "P50": array([5.8, ...]), ...}
+        """
+        if not self._models:
+            raise RuntimeError("Model not trained. Call fit() first.")
+        X_scaled = self._preprocess(X)
+        result = {}
+        for label, model in self._models.items():
+            result[label] = model.predict(X_scaled)
+        return result
+
+    def predict_points(self, X: pd.DataFrame) -> np.ndarray:
+        """Convenience: return the P50 (median) prediction array."""
+        preds = self.predict(X)
+        p50_key = "P50"
+        if p50_key not in preds:
+            available = min(preds.keys(), key=lambda k: abs(float(k[1:]) / 100 - 0.50))
+            logger.warning(f"P50 not trained; falling back to {available}")
+            return preds[available]
+        return preds[p50_key]
+
+    # ------------------------------------------------------------------
+    # Risk / reward utilities
+    # ------------------------------------------------------------------
+    def predict_downside_risk(
+        self, X: pd.DataFrame, threshold: float = 5.5
+    ) -> np.ndarray:
+        """Probability that player score falls below *threshold*.
+
+        Uses CDF interpolation across the trained quantiles.  The returned
+        probability is the fraction of the predictive distribution that lies
+        below the threshold.
+        """
+        preds = self.predict(X)
+        n = len(next(iter(preds.values())))
+        probs = np.zeros(n)
+
+        for i in range(n):
+            q_vals = [preds[label][i] for label in sorted(preds.keys())]
+            q_levels = sorted([float(k[1:]) / 100 for k in preds.keys()])
+
+            if threshold <= q_vals[0]:
+                probs[i] = q_levels[0]
+            elif threshold >= q_vals[-1]:
+                probs[i] = q_levels[-1]
+            else:
+                idx = np.searchsorted(q_vals, threshold)
+                lo_q, hi_q = q_vals[idx - 1], q_vals[idx]
+                lo_level, hi_level = q_levels[idx - 1], q_levels[idx]
+                frac = (threshold - lo_q) / (hi_q - lo_q + 1e-10)
+                probs[i] = lo_level + frac * (hi_level - lo_level)
+
+        return np.clip(probs, 0.0, 1.0)
+
+    def predict_upside(self, X: pd.DataFrame, threshold: float = 7.0) -> np.ndarray:
+        """Probability that player score exceeds *threshold*."""
+        downside = self.predict_downside_risk(X, threshold)
+        return 1.0 - downside
+
+    def value_at_risk_safe(
+        self, X: pd.DataFrame, confidence: float = 0.90
+    ) -> np.ndarray:
+        """VaR-safe estimate: floor that the player exceeds with given confidence.
+
+        For confidence=0.90 returns the score floor that the player exceeds 90%
+        of the time — i.e. the P(100-confidence) quantile.  A higher VaR-safe
+        means more reliable upside.
+        """
+        alpha = 1.0 - confidence
+        preds = self.predict(X)
+        levels = np.array(sorted([float(k[1:]) / 100 for k in preds.keys()]))
+
+        # If the exact alpha was trained return it directly
+        atol = 0.005
+        for label, q_pred in preds.items():
+            if abs(float(label[1:]) / 100 - alpha) <= atol:
+                return np.asarray(q_pred)
+
+        # Otherwise linearly interpolate
+        idx = np.searchsorted(levels, alpha)
+        if idx == 0:
+            label = sorted(preds.keys())[0]
+            return np.asarray(preds[label])
+        if idx >= len(levels):
+            label = sorted(preds.keys())[-1]
+            return np.asarray(preds[label])
+
+        lo_level = levels[idx - 1]
+        hi_level = levels[idx]
+        lo_label = f"P{int(round(lo_level * 100))}"
+        hi_label = f"P{int(round(hi_level * 100))}"
+        # Fall back to closest trained quantile keys
+        lo_label = min(preds.keys(), key=lambda k: abs(float(k[1:]) / 100 - lo_level))
+        hi_label = min(preds.keys(), key=lambda k: abs(float(k[1:]) / 100 - hi_level))
+
+        lo_vals = preds[lo_label]
+        hi_vals = preds[hi_label]
+        frac = (alpha - lo_level) / (hi_level - lo_level + 1e-10)
+        return lo_vals + frac * (hi_vals - lo_vals)

src/models/set_transformer.py

@@ -0,0 +1,726 @@
+"""Set Transformer for team-level valuation.
+
+Models a Fantacalcio roster as an unordered set of players rather than
+a simple sum of individual projections. Captures non-linear team composition
+effects: synergies between players, role balance, and diminishing returns
+from overlapping skill sets.
+
+Two implementation modes:
+- PyTorch: Full Set Transformer with Induced Set Attention Blocks (ISAB)
+  and Pooling by Multihead Attention (PMA).
+- sklearn: "Bag-of-Players" approximation using per-role aggregates,
+  pairwise interactions, and non-linear regression.
+"""
+
+import logging
+import math
+from collections import Counter
+from pathlib import Path
+from typing import Dict, List, Optional, Tuple
+
+import numpy as np
+import pandas as pd
+
+from .base_model import BaseModel
+
+logger = logging.getLogger(__name__)
+
+ROLES = ["P", "D", "C", "A"]
+ROLE_INDEX = {"P": 0, "D": 1, "C": 2, "A": 3}
+
+_EMBEDDING_FEATURES = [
+    "projected_points",
+    "market_value",
+    "ceiling_price",
+    "vor",  # value over replacement
+    "goals_per_game",
+    "assists_per_game",
+    "cs_prob",  # clean sheet probability
+    "minutes_avg",
+    "form_recent",
+    "injury_risk",
+    "starter_prob",
+]
+
+
+def _has_torch() -> bool:
+    try:
+        import torch  # noqa: F401
+        return True
+    except ImportError:
+        return False
+
+
+# ---------------------------------------------------------------------------
+# Torch Set Transformer modules (lazy-built to handle missing torch)
+# ---------------------------------------------------------------------------
+
+
+def _build_torch_classes():
+    """Factory that returns torch nn.Module classes for the Set Transformer.
+
+    Only called when torch is available. All returned classes are proper
+    nn.Module subclasses so they work with nn.ModuleList, nn.Sequential, etc.
+    """
+    import torch
+    import torch.nn as nn
+    import torch.nn.functional as F
+
+    class MAB(nn.Module):
+        """Multihead Attention Block.
+
+            H = LayerNorm(X + Multihead(Q=X, K=Y, V=Y))
+            Out = LayerNorm(H + FFN(H))
+        """
+
+        def __init__(self, dim, n_heads, rff_dim, dropout=0.1):
+            super().__init__()
+            self.dim = dim
+            self.n_heads = n_heads
+            self.head_dim = dim // n_heads
+            self.scale = self.head_dim ** 0.5
+
+            self.ln1 = nn.LayerNorm(dim)
+            self.ln2 = nn.LayerNorm(dim)
+            self.W_q = nn.Linear(dim, dim, bias=False)
+            self.W_k = nn.Linear(dim, dim, bias=False)
+            self.W_v = nn.Linear(dim, dim, bias=False)
+            self.W_o = nn.Linear(dim, dim, bias=False)
+            self.ffn = nn.Sequential(
+                nn.Linear(dim, rff_dim),
+                nn.ReLU(),
+                nn.Dropout(dropout),
+                nn.Linear(rff_dim, dim),
+                nn.Dropout(dropout),
+            )
+
+        def forward(self, X, Y):
+            H = self.ln1(X + self._mh_attention(X, Y, Y))
+            return self.ln2(H + self.ffn(H))
+
+        def _mh_attention(self, Q, K, V):
+            B, N_Q, _ = Q.shape
+            N_K = K.shape[1]
+
+            q = self.W_q(Q).view(B, N_Q, self.n_heads, self.head_dim).transpose(1, 2)
+            k = self.W_k(K).view(B, N_K, self.n_heads, self.head_dim).transpose(1, 2)
+            v = self.W_v(V).view(B, N_K, self.n_heads, self.head_dim).transpose(1, 2)
+
+            attn = torch.matmul(q, k.transpose(-2, -1)) / self.scale
+            attn = torch.softmax(attn, dim=-1)
+            out = torch.matmul(attn, v)
+            out = out.transpose(1, 2).contiguous().view(B, N_Q, self.dim)
+            return self.W_o(out)
+
+    class ISAB(nn.Module):
+        """Induced Set Attention Block.
+
+            X -> MAB(X, MAB(I, X))
+
+        with learnable inducing points I.
+        """
+
+        def __init__(self, dim, n_heads, n_inducing, rff_dim, dropout=0.1):
+            super().__init__()
+            self.mab1 = MAB(dim, n_heads, rff_dim, dropout)
+            self.mab2 = MAB(dim, n_heads, rff_dim, dropout)
+            self.I = nn.Parameter(torch.randn(1, n_inducing, dim) * 0.1)
+            self.n_inducing = n_inducing
+
+        def forward(self, X):
+            B = X.shape[0]
+            I_expanded = self.I.expand(B, -1, -1)
+            H = self.mab1(I_expanded, X)
+            return self.mab2(X, H)
+
+    class PMA(nn.Module):
+        """Pooling by Multihead Attention.
+
+            SEMA = MAB(S, X)  where S is a learnable seed vector.
+        """
+
+        def __init__(self, dim, n_heads, n_seeds, rff_dim, dropout=0.1):
+            super().__init__()
+            self.mab = MAB(dim, n_heads, rff_dim, dropout)
+            self.S = nn.Parameter(torch.randn(1, n_seeds, dim) * 0.1)
+            self.n_seeds = n_seeds
+
+        def forward(self, X):
+            B = X.shape[0]
+            S_expanded = self.S.expand(B, -1, -1)
+            return self.mab(S_expanded, X)
+
+    class SetTransformerTorch(nn.Module):
+        """Full Set Transformer with ISAB layers and PMA pooling."""
+
+        def __init__(self, d_input, d_model=128, n_heads=4, n_layers=3,
+                     n_inducing=32, dropout=0.1):
+            super().__init__()
+            self.d_model = d_model
+            self.input_proj = nn.Linear(d_input, d_model)
+            self.isabs = nn.ModuleList([
+                ISAB(d_model, n_heads, n_inducing, d_model * 2, dropout)
+                for _ in range(n_layers)
+            ])
+            self.pma = PMA(d_model, n_heads, 1, d_model * 2, dropout)
+            self.output_head = nn.Sequential(
+                nn.Linear(d_model, d_model // 2),
+                nn.ReLU(),
+                nn.Dropout(dropout),
+                nn.Linear(d_model // 2, 1),
+            )
+
+        def forward(self, x):
+            h = self.input_proj(x)
+            for isab in self.isabs:
+                h = isab(h)
+            pooled = self.pma(h)
+            pooled = pooled.squeeze(1)
+            return self.output_head(pooled).view(-1)
+
+    return MAB, ISAB, PMA, SetTransformerTorch
+
+
+# ---------------------------------------------------------------------------
+# sklearn fallback: Bag-of-Players aggregator
+# ---------------------------------------------------------------------------
+
+
+class _BagOfPlayers:
+    """Sklearn-based team valuation using set-level aggregate features.
+
+    Instead of learning over the entire permutation-invariant set structure,
+    we compute fixed aggregate statistics per team:
+    - Per-role: mean, max, min, std of each numeric feature
+    - Pairwise cosine similarities between player embeddings
+    - Position entropy
+    - Budget allocation fractions
+    """
+
+    def __init__(self, random_state: int = 42):
+        self.random_state = random_state
+        self.model = None
+        self.fitted = False
+        self._feature_names: List[str] = []
+        self._scaler = None
+        self._n_agg_features = 0
+
+    def _extract_features(self, team_df: pd.DataFrame) -> np.ndarray:
+        df = team_df.copy()
+
+        present_cols = [c for c in _EMBEDDING_FEATURES if c in df.columns]
+        if not present_cols:
+            present_cols = list(df.select_dtypes(include=[np.number]).columns)
+        if "role" not in df.columns:
+            df["role"] = "C"
+
+        df = df.fillna(0.0)
+        features: List[float] = []
+
+        # Per-role aggregates
+        for role in ROLES:
+            subset = df[df["role"] == role][present_cols]
+            if len(subset) == 0:
+                for col in present_cols:
+                    features.extend([0.0, 0.0, 0.0, 0.0])
+                continue
+            values = subset.values.astype(np.float64)
+            for col_idx, _col in enumerate(present_cols):
+                col_vals = values[:, col_idx]
+                features.append(float(np.mean(col_vals)))
+                features.append(float(np.max(col_vals)))
+                features.append(float(np.min(col_vals)))
+                features.append(float(np.std(col_vals)) if len(col_vals) > 1 else 0.0)
+
+        # Role count features
+        n_players = len(df)
+        for role in ROLES:
+            features.append(float((df["role"] == role).sum()) / max(n_players, 1))
+
+        # Position entropy
+        role_counts = Counter(df["role"].tolist())
+        entropy = 0.0
+        for role in ROLES:
+            count = role_counts.get(role, 0)
+            if count > 0:
+                p = count / max(n_players, 1)
+                entropy -= p * math.log(p + 1e-10)
+        features.append(entropy)
+        features.append(float(n_players))
+
+        # Pairwise cosine similarities
+        if "projected_points" in df.columns and "market_value" in df.columns:
+            emb_cols = [c for c in ["projected_points", "market_value", "vor", "starter_prob"]
+                        if c in df.columns]
+            if emb_cols:
+                emb = df[emb_cols].fillna(0.0).values.astype(np.float64)
+                emb_norm = emb / (np.linalg.norm(emb, axis=1, keepdims=True) + 1e-8)
+                sim_matrix = emb_norm @ emb_norm.T
+                n = sim_matrix.shape[0]
+                if n > 1:
+                    triu_idx = np.triu_indices(n, k=1)
+                    triu_vals = sim_matrix[triu_idx]
+                    features.append(float(np.mean(triu_vals)))
+                    features.append(float(np.max(triu_vals)))
+                    features.append(float(np.min(triu_vals)))
+                    features.append(float(np.std(triu_vals)))
+                else:
+                    features.extend([0.0, 0.0, 0.0, 0.0])
+        else:
+            features.extend([0.0, 0.0, 0.0, 0.0])
+
+        # Cross-role interaction: dot products between role-group means
+        role_means = {}
+        for role in ROLES:
+            subset = df[df["role"] == role]
+            if len(subset) == 0 or "projected_points" not in subset.columns:
+                role_means[role] = np.zeros(len(present_cols))
+            else:
+                role_means[role] = subset[present_cols].fillna(0.0).mean().values
+
+        for i, r1 in enumerate(ROLES):
+            for r2 in ROLES[i + 1:]:
+                if np.linalg.norm(role_means[r1]) > 1e-8 and np.linalg.norm(role_means[r2]) > 1e-8:
+                    dot = np.dot(role_means[r1], role_means[r2]) / max(
+                        np.linalg.norm(role_means[r1]) * np.linalg.norm(role_means[r2]), 1e-8
+                    )
+                else:
+                    dot = 0.0
+                features.append(float(dot))
+
+        # Budget allocation features
+        if "market_value" in df.columns:
+            total_team_value = df["market_value"].sum()
+            for role in ROLES:
+                subset = df[df["role"] == role]
+                role_value = subset["market_value"].sum() if len(subset) > 0 else 0.0
+                features.append(role_value / max(total_team_value, 1))
+        else:
+            for role in ROLES:
+                features.append(0.0)
+
+        return np.array(features, dtype=np.float64)
+
+    def fit(self, teams_data: List[pd.DataFrame], team_values: List[float]):
+        from sklearn.ensemble import RandomForestRegressor
+        from sklearn.preprocessing import StandardScaler
+
+        if len(teams_data) == 0:
+            raise ValueError("teams_data cannot be empty")
+
+        X_list = [self._extract_features(df) for df in teams_data]
+        X = np.vstack(X_list)
+        y = np.array(team_values, dtype=np.float64).ravel()
+
+        if len(X) != len(y):
+            raise ValueError(f"Mismatch: {len(X)} teams vs {len(y)} values")
+
+        self._scaler = StandardScaler()
+        X_scaled = self._scaler.fit_transform(X)
+
+        self.model = RandomForestRegressor(
+            n_estimators=200,
+            max_depth=12,
+            min_samples_leaf=5,
+            random_state=self.random_state,
+            n_jobs=-1,
+        )
+        self.model.fit(X_scaled, y)
+        self.fitted = True
+        self._n_agg_features = X.shape[1]
+
+        # Log feature importance
+        if hasattr(self.model, "feature_importances_"):
+            top_n = min(10, len(self.model.feature_importances_))
+            top_idx = np.argsort(self.model.feature_importances_)[::-1][:top_n]
+            logger.info(
+                f"BagOfPlayers fitted: {len(teams_data)} teams, "
+                f"{X.shape[1]} features. Top features: {top_idx.tolist()}"
+            )
+
+        return self
+
+    def predict(self, team_df: pd.DataFrame) -> float:
+        if not self.fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+        feats = self._extract_features(team_df).reshape(1, -1)
+        feats_scaled = self._scaler.transform(feats)
+        pred = self.model.predict(feats_scaled)
+        return float(pred[0])
+
+
+# ---------------------------------------------------------------------------
+# Set Transformer main model
+# ---------------------------------------------------------------------------
+
+
+class SetTransformer(BaseModel):
+    """Set-based team valuation model.
+
+    Values an entire Fantacalcio roster as a permutation-invariant set,
+    capturing non-additive synergies between players.
+
+    Args:
+        model_dir: directory for model persistence.
+        d_model: hidden dimension.
+        n_heads: attention heads.
+        n_layers: ISAB layers.
+        use_torch: force PyTorch mode (auto-detect by default).
+    """
+
+    def __init__(
+        self,
+        model_dir: str = "models_trained",
+        d_model: int = 128,
+        n_heads: int = 4,
+        n_layers: int = 3,
+        use_torch: Optional[bool] = None,
+    ):
+        super().__init__(model_dir)
+        self.d_model = d_model
+        self.n_heads = n_heads
+        self.n_layers = n_layers
+        self.use_torch = use_torch if use_torch is not None else _has_torch()
+
+        self._torch_model = None
+        self._sklearn_model: Optional[_BagOfPlayers] = None
+        self._d_input: Optional[int] = None
+        self._optimal_tau: float = 0.02
+        self._device: Optional[str] = None
+
+        self._trained = False
+        self._using_torch = False
+
+    def _init_torch_model(self, d_input: int):
+        import torch
+        _, _, _, SetTransformerTorch = _build_torch_classes()
+        self._torch_model = SetTransformerTorch(
+            d_input=d_input,
+            d_model=self.d_model,
+            n_heads=self.n_heads,
+            n_layers=self.n_layers,
+            n_inducing=32,
+            dropout=0.1,
+        )
+        self._device = "cuda" if torch.cuda.is_available() else "cpu"
+        self._torch_model.to(self._device)
+        self._d_input = d_input
+        self._using_torch = True
+        logger.info(f"SetTransformer (torch) initialized on {self._device}")
+
+    def _init_sklearn_model(self):
+        self._sklearn_model = _BagOfPlayers(random_state=42)
+        self._using_torch = False
+        logger.info("SetTransformer (sklearn fallback) initialized")
+
+    # ------------------------------------------------------------------
+    # Data preparation
+    # ------------------------------------------------------------------
+
+    def _prepare_set(self, team_df: pd.DataFrame) -> np.ndarray:
+        """Extract player feature vectors from a team DataFrame.
+
+        Returns array of shape (n_players, d_input).
+        """
+        df = team_df.copy()
+        present_cols = [c for c in _EMBEDDING_FEATURES if c in df.columns]
+        if not present_cols:
+            present_cols = list(df.select_dtypes(include=[np.number]).columns)
+
+        df = df.fillna(0.0)
+        if "projected_points" not in df.columns and len(present_cols) > 0:
+            df["projected_points"] = df[present_cols[0]]
+
+        for col in _EMBEDDING_FEATURES:
+            if col not in df.columns:
+                df[col] = 0.0
+
+        features = df[_EMBEDDING_FEATURES].values.astype(np.float64)
+        return features
+
+    def _collate_sets(self, teams_data: List[pd.DataFrame]) -> List[np.ndarray]:
+        return [self._prepare_set(df) for df in teams_data]
+
+    # ------------------------------------------------------------------
+    # Fit
+    # ------------------------------------------------------------------
+
+    def fit(self, teams_data: List[pd.DataFrame], team_values: List[float], **kwargs):
+        """Fit the Set Transformer on team-level data.
+
+        Args:
+            teams_data: list of DataFrames, one per team, each containing
+                        player-level features.
+            team_values: list of total season points for each team.
+        """
+        if len(teams_data) == 0:
+            raise ValueError("teams_data cannot be empty")
+        if len(teams_data) != len(team_values):
+            raise ValueError(
+                f"Length mismatch: {len(teams_data)} teams, {len(team_values)} values"
+            )
+
+        if self.use_torch and _has_torch():
+            return self._fit_torch(teams_data, team_values, **kwargs)
+        else:
+            if self.use_torch and not _has_torch():
+                logger.warning("torch requested but not installed. Falling back to sklearn.")
+            return self._fit_sklearn(teams_data, team_values)
+
+    def _fit_torch(
+        self,
+        teams_data: List[pd.DataFrame],
+        team_values: List[float],
+        epochs: int = 200,
+        batch_size: int = 16,
+        lr: float = 1e-3,
+        weight_decay: float = 1e-5,
+    ) -> "SetTransformer":
+        import torch
+        import torch.nn as nn
+        import torch.optim as optim
+
+        sets = self._collate_sets(teams_data)
+        d_input = sets[0].shape[1]
+        self._init_torch_model(d_input)
+
+        y = np.array(team_values, dtype=np.float64)
+        y_mean = y.mean()
+        y_std = max(y.std(), 1e-8)
+        y_norm = (y - y_mean) / y_std
+        self._y_mean = y_mean
+        self._y_std = y_std
+
+        model = self._torch_model
+        optimizer = optim.Adam(model.parameters(), lr=lr, weight_decay=weight_decay)
+        loss_fn = nn.MSELoss()
+
+        n_samples = len(sets)
+        best_loss = float("inf")
+        best_state = model.state_dict()
+
+        for epoch in range(epochs):
+            perm = torch.randperm(n_samples)
+            epoch_loss = 0.0
+
+            for start in range(0, n_samples, batch_size):
+                idx = perm[start:start + batch_size]
+                batch_loss = torch.tensor(0.0, device=self._device)
+
+                for i in idx:
+                    x_np = sets[i.item()]
+                    x_t = torch.tensor(
+                        x_np, dtype=torch.float32, device=self._device
+                    ).unsqueeze(0)
+                    y_t = torch.tensor(
+                        [float(y_norm[i.item()])], dtype=torch.float32, device=self._device
+                    )
+                    pred = model(x_t)
+                    loss = loss_fn(pred, y_t)
+                    batch_loss = batch_loss + loss
+
+                batch_loss = batch_loss / len(idx)
+                optimizer.zero_grad()
+                batch_loss.backward()
+                torch.nn.utils.clip_grad_norm_(model.parameters(), 5.0)
+                optimizer.step()
+                epoch_loss += batch_loss.item() * len(idx)
+
+            epoch_loss /= n_samples
+
+            if epoch_loss < best_loss:
+                best_loss = epoch_loss
+                best_state = {k: v.cpu().clone() for k, v in model.state_dict().items()}
+
+            if (epoch + 1) % max(epochs // 10, 1) == 0:
+                logger.info(f"Epoch {epoch + 1}/{epochs} | loss={epoch_loss:.6f}")
+
+        if best_state is not None:
+            model.load_state_dict(best_state)
+
+        self._trained = True
+        logger.info(
+            f"SetTransformer (torch) trained: {n_samples} teams, "
+            f"best_loss={best_loss:.6f}, d_model={self.d_model}"
+        )
+        return self
+
+    def _fit_sklearn(
+        self, teams_data: List[pd.DataFrame], team_values: List[float]
+    ) -> "SetTransformer":
+        self._init_sklearn_model()
+        self._sklearn_model.fit(teams_data, team_values)
+        self._trained = True
+        self._using_torch = False
+        return self
+
+    # ------------------------------------------------------------------
+    # Predict
+    # ------------------------------------------------------------------
+
+    def predict(self, team_roster_df: pd.DataFrame):
+        """Predict total team value (season-long points)."""
+        if not self._trained:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+
+        val = self._predict_single(team_roster_df)
+        return val
+
+    def _predict_single(self, team_roster_df: pd.DataFrame) -> float:
+        if self._using_torch and self._torch_model is not None:
+            import torch
+            x_np = self._prepare_set(team_roster_df)
+            if len(x_np) == 0:
+                return 0.0
+            x_t = torch.tensor(x_np, dtype=torch.float32, device=self._device).unsqueeze(0)
+            with torch.no_grad():
+                raw = self._torch_model(x_t)
+            val = raw.item() * getattr(self, "_y_std", 1.0) + getattr(self, "_y_mean", 0.0)
+            return float(val)
+        elif self._sklearn_model is not None:
+            return self._sklearn_model.predict(team_roster_df)
+        return 0.0
+
+    def predict_batch(self, teams: List[pd.DataFrame]) -> np.ndarray:
+        if not self._trained:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+        return np.array([self._predict_single(df) for df in teams])
+
+    # ------------------------------------------------------------------
+    # Marginal value analysis
+    # ------------------------------------------------------------------
+
+    def value_added(self, team_roster_df: pd.DataFrame, new_player) -> float:
+        """Marginal value: delta when adding new_player to the team."""
+        baseline = self._predict_single(team_roster_df)
+        if isinstance(new_player, pd.DataFrame):
+            augmented = pd.concat([team_roster_df, new_player], ignore_index=True)
+        else:
+            augmented = pd.concat(
+                [team_roster_df, pd.DataFrame([new_player])], ignore_index=True
+            )
+        augmented_val = self._predict_single(augmented)
+        return augmented_val - baseline
+
+    def value_removed(self, team_roster_df: pd.DataFrame, removed_player) -> float:
+        """Marginal loss: delta when removing a player (by index or name)."""
+        baseline = self._predict_single(team_roster_df)
+        if isinstance(removed_player, str):
+            reduced = team_roster_df[team_roster_df["name"] != removed_player]
+        else:
+            reduced = team_roster_df.drop(team_roster_df.index[removed_player])
+        reduced_val = self._predict_single(reduced)
+        return baseline - reduced_val
+
+    def optimal_replacement(
+        self,
+        team_roster_df: pd.DataFrame,
+        candidate_pool: pd.DataFrame,
+        to_replace: List,
+    ) -> Dict:
+        """For each player to replace, rank candidates by predicted team value delta.
+
+        Args:
+            team_roster_df: current team roster.
+            candidate_pool: DataFrame of free-agent candidates.
+            to_replace: list of player names (str) or indices (int) in team_roster_df
+                       to consider replacing.
+
+        Returns:
+            dict mapping player_name -> DataFrame of candidates ranked by delta.
+        """
+        results = {}
+        for rp in to_replace:
+            if isinstance(rp, str):
+                base_team = team_roster_df[team_roster_df["name"] != rp]
+                key = rp
+            else:
+                base_team = team_roster_df.drop(team_roster_df.index[rp])
+                key = rp
+            deltas = []
+            for _, cand in candidate_pool.iterrows():
+                cand_dict = cand.to_dict()
+                new_team = pd.concat(
+                    [base_team, pd.DataFrame([cand_dict])], ignore_index=True
+                )
+                new_val = self._predict_single(new_team)
+                current_val = self._predict_single(team_roster_df)
+                deltas.append({
+                    "candidate": cand.get("name", str(cand.name)),
+                    "role": cand.get("role", ""),
+                    "projected_points": cand.get("projected_points", 0),
+                    "team_value_delta": new_val - current_val,
+                })
+
+            results[key] = pd.DataFrame(deltas).sort_values(
+                "team_value_delta", ascending=False
+            )
+        return results
+
+    # ------------------------------------------------------------------
+    # Roster analysis
+    # ------------------------------------------------------------------
+
+    def get_role_synergies(self, team_roster_df: pd.DataFrame) -> pd.DataFrame:
+        """Analyze which role combinations maximize team value.
+
+        Systematically tests removing one role at a time and measures
+        the value contribution per role. Returns a DataFrame with
+        per-role synergy scores.
+        """
+        baseline = self._predict_single(team_roster_df)
+        results = []
+
+        for role in ROLES:
+            role_players = team_roster_df[team_roster_df["role"] == role]
+            if len(role_players) == 0:
+                synergy = 0.0
+                per_player = 0.0
+            else:
+                without_role = team_roster_df[team_roster_df["role"] != role]
+                without_val = self._predict_single(without_role) if len(without_role) > 0 else 0.0
+                synergy = baseline - without_val
+                per_player = synergy / len(role_players)
+
+            results.append({
+                "role": role,
+                "n_players": len(role_players),
+                "absolute_synergy": synergy,
+                "synergy_per_player": per_player,
+                "synergy_pct": (synergy / max(baseline, 1)) * 100,
+            })
+
+        return pd.DataFrame(results).sort_values("absolute_synergy", ascending=False)
+
+    def get_redundancy_score(self, team_roster_df: pd.DataFrame) -> float:
+        """Compute 0-1 redundancy score for the roster.
+
+        High redundancy = lots of overlapping skill sets = diminishing
+        returns beyond sum of individual projections.
+
+        Uses the ratio of (sum of individual values) / (team value) as
+        a proxy: if team value << sum of parts, players are redundant.
+
+        Returns:
+            float between 0 (perfect synergy) and 1 (maximum redundancy).
+        """
+        team_val = self._predict_single(team_roster_df)
+        if team_val <= 0:
+            return 0.0
+
+        individual_sum = 0.0
+        for _, player in team_roster_df.iterrows():
+            solo = pd.DataFrame([player.to_dict()])
+            individual_sum += self._predict_single(solo)
+
+        if individual_sum <= 0:
+            return 0.0
+
+        ratio = team_val / individual_sum
+
+        # Map: ratio ~1 means additive (no synergy, no redundancy)
+        #      ratio >1 means synergy
+        #      ratio <<1 means redundancy
+        redundancy = np.clip(1.0 - ratio, 0.0, 1.0)
+
+        return float(redundancy)

src/models/survival_model.py

@@ -0,0 +1,354 @@
+"""Minutes-played survival model via Weibull AFT.
+
+Models the distribution of minutes played per matchweek using Weibull
+Accelerated Failure Time.  Supports both lifelines (preferred) and a pure-scipy
+MLE fallback so the module works in minimal environments.
+
+Provides:
+- Expected minutes / confidence intervals
+- Starter probability (≥60 min)
+- Full-match probability (90 min)
+"""
+
+import logging
+import math
+from typing import Optional, Tuple
+
+import numpy as np
+import pandas as pd
+from sklearn.linear_model import LinearRegression
+from sklearn.preprocessing import StandardScaler
+
+from .base_model import BaseModel
+
+logger = logging.getLogger(__name__)
+
+# ---------------------------------------------------------------------------
+# Feature-selection keyword list
+# ---------------------------------------------------------------------------
+_MINUTE_KEYWORDS = [
+    "minute", "game", "rest", "fatigue", "age", "injury",
+    "played", "starter", "bench", "appearance",
+    "recovery", "rotation", "squad", "season",
+    "match", "form", "fitness",
+]
+
+
+def _select_survival_features(X: pd.DataFrame) -> list:
+    """Pick columns whose name contains any survival-relevant keyword."""
+    lower_cols = {c: str(c).lower() for c in X.columns}
+    selected = [
+        c for c, cl in lower_cols.items()
+        if any(kw in cl for kw in _MINUTE_KEYWORDS)
+    ]
+    if not selected:
+        selected = list(X.select_dtypes(include=[np.number]).columns[:20])
+        logger.info("No keyword-matched survival features; using first 20 numeric columns")
+    else:
+        logger.info(f"Selected {len(selected)} survival features via keyword matching")
+    return selected
+
+
+# ===================================================================
+# Weibull helper functions for the scipy fallback
+# ===================================================================
+
+def _weibull_log_likelihood(params, X, t, event, eps=1e-10):
+    """Negative log-likelihood for Weibull AFT model.
+
+    Parameters
+    ----------
+    params : ndarray (p_features + 1,)
+        First p entries: beta (coefficients for X).
+        Last entry: log_k (log shape parameter ensures k > 0).
+    X : ndarray (n, p)
+        Scaled feature matrix.
+    t : ndarray (n,)
+        Observed durations (minutes played).
+    event : ndarray (n,)
+        0 → exact failure (subbed off),  1 → right-censored (completed 90).
+    eps : float
+        Small epsilon for numerical stability.
+
+    Returns
+    -------
+    neg_ll : float
+        Negative log-likelihood (to be minimized).
+    """
+    p = X.shape[1]
+    beta = params[:p]
+    log_k = params[p]
+    k = np.exp(log_k) + eps
+
+    log_lambda = X.dot(beta)              # log(λ_i) = X_i * beta
+    lambda_ = np.exp(log_lambda) + eps
+    log_t = np.log(np.maximum(t, eps))
+    z = t / lambda_
+
+    # Log-PDF for uncensored (event == 0)
+    log_pdf = np.log(k) - log_lambda + (k - 1.0) * (log_t - log_lambda) - z ** k
+
+    # Log-SF for censored (event == 1)
+    log_sf = -(z ** k)
+
+    # event==1 → censored → use SF; event==0 → observed → use PDF
+    ll = np.where(event == 1, log_sf, log_pdf)
+    return -ll.sum()
+
+
+def _fit_weibull_mle(X, t, event):
+    """Fit Weibull AFT via scipy MLE.
+
+    Returns
+    -------
+    beta : ndarray (p,)
+        Feature coefficients (scaled to original duration range).
+    k : float
+        Shape parameter.
+    t_scale : float
+        Scale factor to convert normalized predictions back to minutes.
+    """
+    from scipy.optimize import minimize
+
+    n, p = X.shape
+    t_scale = max(t.max(), 1.0)
+    t_norm = np.clip(t / t_scale, 1e-6, 1.0)
+    log_t_norm = np.log(np.maximum(t_norm, 1e-9))
+
+    lr = LinearRegression(fit_intercept=False)
+    lr.fit(X, log_t_norm)
+    beta0 = np.clip(lr.coef_.copy(), -5, 5)
+
+    bounds = [(-10, 10)] * p + [(-5, 3)]
+    init = np.concatenate([beta0, [0.0]])
+
+    result = minimize(
+        _weibull_log_likelihood,
+        init,
+        args=(X, t_norm, event),
+        method="L-BFGS-B",
+        bounds=bounds,
+        options={"maxiter": 2000, "ftol": 1e-10},
+    )
+    if not result.success:
+        logger.warning(f"Weibull MLE did not converge: {result.message}")
+
+    beta = result.x[:p]
+    k = max(np.exp(result.x[p]), 1e-4)
+    return beta, k, t_scale
+
+
+# ===================================================================
+# MinutesSurvivalModel
+# ===================================================================
+
+class MinutesSurvivalModel(BaseModel):
+    """Weibull AFT model for minutes-played distribution.
+
+    Parameters
+    ----------
+    model_dir : str
+        Directory for persisting trained models.
+    force_scipy : bool
+        If True, use the pure-scipy MLE fallback even when lifelines
+        is installed.
+    """
+
+    def __init__(
+        self,
+        model_dir: str = "models_trained",
+        force_scipy: bool = False,
+    ):
+        super().__init__(model_dir)
+        self.force_scipy = force_scipy
+
+        self.scaler = StandardScaler()
+        self.feature_names = None
+
+        # Weibull parameters
+        self._beta = None           # feature coefficients → log(λ)
+        self._k = None              # shape parameter
+        self._afitter = None        # lifelines WeibullAFTFitter instance (if used)
+        self._t_scale = 90.0
+        self._use_lifelines = False
+
+    # ------------------------------------------------------------------
+    # Fit
+    # ------------------------------------------------------------------
+    def fit(
+        self,
+        X: pd.DataFrame,
+        durations: np.ndarray,
+        events: np.ndarray,
+        **kwargs,
+    ):
+        """Fit the Weibull AFT model.
+
+        Args:
+            X: Feature matrix (one row per player-match).
+            durations: Minutes played (0–90); `y` alias for BaseModel compat.
+            events:
+                0 → exact duration observed (subbed off before 90).
+                1 → right-censored (player completed the full 90 minutes).
+        """
+        self.feature_names = _select_survival_features(X)
+        X_clean = X[self.feature_names].select_dtypes(include=[np.number]).fillna(0)
+        X_scaled = self.scaler.fit_transform(X_clean)
+
+        durations = np.asarray(durations, dtype=float)
+        events = np.asarray(events, dtype=int)
+
+        # Try lifelines first --------------------------------------------------
+        if not self.force_scipy:
+            try:
+                import lifelines  # noqa: F401
+                from lifelines import WeibullAFTFitter
+
+                df = pd.DataFrame(X_scaled, columns=self.feature_names)
+                df["duration"] = durations
+                df["event"] = events
+
+                aft = WeibullAFTFitter()
+                aft.fit(df, duration_col="duration", event_col="event")
+                self._afitter = aft
+                self._use_lifelines = True
+                self._t_scale = 1.0  # lifelines works in original duration scale
+                logger.info(
+                    "MinutesSurvivalModel fitted via lifelines "
+                    f"(n={len(durations)}, features={len(self.feature_names)})"
+                )
+                return self
+            except ImportError:
+                logger.info("lifelines not installed; falling back to scipy MLE")
+            except Exception as exc:
+                logger.warning(f"lifelines failed ({exc}); falling back to scipy MLE")
+
+        # Scipy fallback -------------------------------------------------------
+        self._use_lifelines = False
+        self._beta, self._k, self._t_scale = _fit_weibull_mle(X_scaled, durations, events)
+        logger.info(
+            f"MinutesSurvivalModel fitted via scipy MLE "
+            f"(n={len(durations)}, features={len(self.feature_names)}, "
+            f"k={self._k:.3f})"
+        )
+        return self
+
+    # ------------------------------------------------------------------
+    # Predict (BaseModel interface — returns expected minutes)
+    # ------------------------------------------------------------------
+    def predict(self, X: pd.DataFrame) -> np.ndarray:
+        """Return expected minutes (E[T]) — BaseModel interface."""
+        return self.predict_expected_minutes(X)
+
+    # ------------------------------------------------------------------
+    # Preprocessing
+    # ------------------------------------------------------------------
+    def _preprocess(self, X: pd.DataFrame) -> np.ndarray:
+        if self.feature_names is None:
+            raise RuntimeError("Model not trained. Call fit() first.")
+        X_c = X[self.feature_names].select_dtypes(include=[np.number]).fillna(0)
+        return self.scaler.transform(X_c)
+
+    # ------------------------------------------------------------------
+    # Core distribution
+    # ------------------------------------------------------------------
+    def predict_distribution(
+        self, X: pd.DataFrame
+    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+        """Return (expected_minutes, lower_bound, upper_bound).
+
+        ``lower_bound`` and ``upper_bound`` are approximate 95 % confidence
+        intervals derived from the Weibull variance.
+        """
+        X_scaled = self._preprocess(X)
+
+        if self._use_lifelines and self._afitter is not None:
+            df = pd.DataFrame(X_scaled, columns=self.feature_names)
+            # Lifelines returns median survival in its summary; we approximate
+            # expected minutes using the median and estimated shape.
+            median = self._afitter.predict_median(df).values.flatten()
+            # Heuristic: for Weibull, E[T] ≈ median / (ln 2)^(1/k).
+            # Derive k from the lifelines summary if possible, else guess ~1.
+            try:
+                summary = self._afitter.summary
+                log_k = summary.loc["lambda_", "coef"]
+                k = 1.0 / np.exp(log_k) if abs(log_k) > 1e-8 else 1.0
+            except Exception:
+                k = 1.0
+            expected = median * np.exp(np.log(np.log(2)) / k)
+            # Std via coefficient of variation
+            coef_var = np.sqrt(np.exp(
+                np.log(math.gamma(1 + 2 / k)) - 2 * np.log(math.gamma(1 + 1 / k))
+            ))
+            std = expected * coef_var
+            lower = np.maximum(0, expected - 1.96 * std)
+            upper = np.minimum(90, expected + 1.96 * std)
+            return expected, lower, upper
+
+        # Scipy / stored parameters (normalized scale, convert to minutes)
+        log_lambda = X_scaled.dot(self._beta)
+        lambda_ = np.exp(log_lambda) * self._t_scale
+        k = self._k
+
+        # Expected value: λ * Γ(1 + 1/k)
+        gamma_1 = math.gamma(1.0 + 1.0 / k)
+        expected = lambda_ * gamma_1
+
+        # Variance = λ² * (Γ(1+2/k) - Γ²(1+1/k))
+        gamma_2 = math.gamma(1.0 + 2.0 / k)
+        var = (lambda_ ** 2) * (gamma_2 - gamma_1 ** 2)
+        std = np.sqrt(np.maximum(var, 0.01))
+
+        lower = np.maximum(0, expected - 1.96 * std)
+        upper = np.minimum(90, expected + 1.96 * std)
+        expected = np.clip(expected, 0, 90)
+        return expected, lower, upper
+
+    def predict_expected_minutes(self, X: pd.DataFrame) -> np.ndarray:
+        """Return E[minutes] for each row."""
+        expected, _, _ = self.predict_distribution(X)
+        return expected
+
+    def predict_full_match_probability(self, X: pd.DataFrame) -> np.ndarray:
+        """Probability the player completes 90 minutes: P(T ≥ 90)."""
+        X_scaled = self._preprocess(X)
+
+        if self._use_lifelines and self._afitter is not None:
+            try:
+                surv = self._afitter.predict_survival_function(
+                    pd.DataFrame(X_scaled, columns=self.feature_names),
+                    times=[90.0],
+                )
+                return 1.0 - surv.values.flatten()
+            except Exception:
+                pass
+
+        log_lambda = X_scaled.dot(self._beta)
+        lambda_ = np.exp(log_lambda) * self._t_scale
+        surv = np.exp(-((90.0 / lambda_) ** self._k))
+        prob = 1.0 - surv
+        return np.clip(prob, 0.0, 1.0)
+
+    def predict_starter_probability(
+        self, X: pd.DataFrame, min_minutes: float = 60.0
+    ) -> np.ndarray:
+        """Probability player plays at least *min_minutes* (default 60).
+
+        Useful as a "likely starter" proxy.
+        """
+        X_scaled = self._preprocess(X)
+
+        if self._use_lifelines and self._afitter is not None:
+            try:
+                surv = self._afitter.predict_survival_function(
+                    pd.DataFrame(X_scaled, columns=self.feature_names),
+                    times=[min_minutes],
+                )
+                return surv.values.flatten()
+            except Exception:
+                pass
+
+        log_lambda = X_scaled.dot(self._beta)
+        lambda_ = np.exp(log_lambda) * self._t_scale
+        surv = np.exp(-((min_minutes / lambda_) ** self._k))
+        return np.clip(surv, 0.0, 1.0)

src/optimization/__init__.py

@@ -1 +1,14 @@
 """Optimization modules for auction and lineup selection."""
+
+from .auction_solver import AuctionSolver, AuctionConfig, PlayerValuation
+from .lineup_solver import LineupSolver, LineupConstraints, PlayerScore, MCTSNode
+from .opponent_model import OpponentModel
+from .transfer_analyzer import TransferAnalyzer
+
+# Phase 2: Bandit, opponent bidding, budget optimization
+from .bandit_auction import BanditAuctionSolver
+from .opponent_bidding_model import OpponentBidModel
+from .budget_optimizer import BudgetOptimizer
+
+# Phase 3: Reinforcement learning auction agent
+from .rl_auction_agent import AuctionEnv, RLAuctionPolicy, RLAuctionTrainer, QNetwork

src/optimization/bandit_auction.py

@@ -0,0 +1,359 @@
+"""Contextual Multi-Armed Bandit for live auction bidding decisions.
+
+Uses Thompson Sampling with Beta-distributed posteriors over discrete bid levels.
+Falls back to UCB when exploration depth is insufficient.
+Integrates with AuctionConfig from auction_solver.py.
+"""
+
+import logging
+from dataclasses import dataclass, field
+from typing import Dict, List, Optional, Tuple
+
+import numpy as np
+import pandas as pd
+from scipy.stats import beta as beta_dist
+
+logger = logging.getLogger(__name__)
+
+# Default bid arms as fractions of total budget
+DEFAULT_BID_ARMS = np.array(
+    [0.0, 0.005, 0.01, 0.02, 0.03, 0.05, 0.08, 0.12, 0.18, 0.25],
+    dtype=np.float64,
+)
+
+ROLES = ["P", "D", "C", "A"]
+ROLE_SCARCITY = {"P": 3, "D": 8, "C": 8, "A": 6}
+ROLE_POOL_SIZE = {"P": 4, "D": 22, "C": 24, "A": 12}
+
+
+@dataclass
+class BanditArmState:
+    alpha: float = 1.0
+    beta: float = 1.0
+    trials: int = 0
+    wins: float = 0.0
+
+
+@dataclass
+class AuctionState:
+    budget_remaining: float = 500.0
+    total_budget: float = 500.0
+    slots_filled: Dict[str, int] = field(default_factory=lambda: {"P": 0, "D": 0, "C": 0, "A": 0})
+    slots_total: Dict[str, int] = field(default_factory=lambda: {"P": 3, "D": 8, "C": 8, "A": 6})
+    round_number: int = 1
+    opponent_budgets: List[float] = field(default_factory=list)
+
+
+@dataclass
+class PlayerContext:
+    name: str
+    role: str
+    projected_points: float
+    role_scarcity: float
+    value_over_replacement: float
+    budget_remaining_fraction: float
+    slots_remaining_in_role: int
+    round_number: int
+    opponent_budget_avg: float
+
+
+def _get_attr(obj, key, default=None):
+    """Get attribute or dict key from an object."""
+    if isinstance(obj, dict):
+        return obj.get(key, default)
+    return getattr(obj, key, default)
+
+
+class BanditAuctionSolver:
+    """Thompson Sampling bandit for live auction bid selection.
+
+    Arms are discrete bid fractions. Each (role, scarcity_level) maintains
+    independent Beta posteriors. Falls back to UCB when total observations
+    for a context group are < 50.
+    """
+
+    def __init__(
+        self,
+        bid_arms: Optional[np.ndarray] = None,
+        total_budget: float = 500.0,
+        min_obs_for_ts: int = 50,
+        ucb_exploration: float = 1.414,
+        config=None,
+    ):
+        if config is not None:
+            from .auction_solver import AuctionConfig
+            total_budget = config.total_budget
+        self.bid_arms = bid_arms if bid_arms is not None else DEFAULT_BID_ARMS
+        self.n_arms = len(self.bid_arms)
+        self.total_budget = total_budget
+        self.min_obs_for_ts = min_obs_for_ts
+        self.ucb_exploration = ucb_exploration
+        self.bid_fractions = self.bid_arms
+
+        self.posteriors: Dict[Tuple[str, int], List[BanditArmState]] = {}
+        for role in ROLES:
+            self._ensure_posteriors(role, 1)
+
+    def _context_key(self, role: str, scarcity_level: int) -> Tuple[str, int]:
+        return (role, scarcity_level)
+
+    def _ensure_posteriors(self, role: str, n_slots_remaining: int):
+        scarcity = max(1, n_slots_remaining)
+        key = self._context_key(role, scarcity)
+        if key not in self.posteriors:
+            self.posteriors[key] = [
+                BanditArmState(alpha=1.0, beta=1.0, trials=0, wins=0.0)
+                for _ in range(self.n_arms)
+            ]
+            logger.debug(f"Initialized bandit posteriors for role={role}, scarcity={scarcity}")
+
+    def _total_obs(self, key: Tuple[str, int]) -> int:
+        arms = self.posteriors.get(key, [])
+        return sum(a.trials for a in arms)
+
+    def compute_context(
+        self,
+        player,
+        auction_state,
+        pool_stats: Optional[dict] = None,
+    ) -> PlayerContext:
+        """Build context vector for a player given current auction state.
+
+        Args:
+            player: dict or object with name, role, projected_points attributes.
+            auction_state: current AuctionState (or dict with same keys).
+            pool_stats: optional dict with role-level pool means and stds.
+
+        Returns:
+            PlayerContext dataclass with all context features.
+        """
+        role = _get_attr(player, "role")
+        points = float(_get_attr(player, "projected_points", 6.5))
+        name = _get_attr(player, "name", "unknown")
+
+        total_slots = _get_attr(auction_state, "slots_total", {}).get(role, 1)
+        filled = _get_attr(auction_state, "slots_filled", {}).get(role, 0)
+        slots_remaining = max(total_slots - filled, 0)
+        scarcity = max(slots_remaining / max(total_slots, 1), 0.05)
+
+        budget_remaining = _get_attr(auction_state, "budget_remaining", 500.0)
+        total_budget_attr = _get_attr(auction_state, "total_budget", 500.0)
+        budget_fraction = budget_remaining / max(total_budget_attr, 1)
+
+        opponent_budget_avg = 0.0
+        opponent_budgets = _get_attr(auction_state, "opponent_budgets", [])
+        if opponent_budgets:
+            opponent_budget_avg = float(np.mean(opponent_budgets))
+        elif total_budget_attr > 0:
+            opponent_budget_avg = total_budget_attr * 0.6
+
+        if pool_stats and role in pool_stats:
+            role_mean = pool_stats[role].get("mean", 0.0)
+            role_std = pool_stats[role].get("std", 1.0)
+            points_z = (points - role_mean) / max(role_std, 0.01) if role_std > 0 else 0.0
+        else:
+            points_z = points / 15.0
+
+        vor = max(points - 6.5, 0.0)
+
+        return PlayerContext(
+            name=name,
+            role=role,
+            projected_points=points,
+            role_scarcity=scarcity,
+            value_over_replacement=vor,
+            budget_remaining_fraction=budget_fraction,
+            slots_remaining_in_role=slots_remaining,
+            round_number=_get_attr(auction_state, "round_number", 1),
+            opponent_budget_avg=opponent_budget_avg,
+        )
+
+    def select_bid(
+        self,
+        player,
+        auction_state,
+        pool_stats: Optional[dict] = None,
+    ) -> Tuple[int, float]:
+        """Select bid arm using Thompson Sampling (or UCB fallback).
+
+        Returns:
+            (arm_index, bid_amount_in_credits)
+        """
+        ctx = self.compute_context(player, auction_state, pool_stats)
+        role = ctx.role
+        n_slots = ctx.slots_remaining_in_role
+
+        self._ensure_posteriors(role, n_slots)
+        key = self._context_key(role, n_slots)
+        arms = self.posteriors[key]
+        total_obs = sum(a.trials for a in arms)
+
+        if total_obs >= self.min_obs_for_ts:
+            samples = [float(np.random.beta(a.alpha, max(a.beta, 0.01))) for a in arms]
+            arm_idx = int(np.argmax(samples))
+            logger.debug(
+                f"Thompson Sampling: role={role}, arms_sampled={samples[:5]}..., "
+                f"selected_arm={arm_idx}"
+            )
+        else:
+            values = []
+            for _i, arm in enumerate(arms):
+                if arm.trials == 0:
+                    values.append(float("inf"))
+                else:
+                    mean = arm.wins / arm.trials
+                    bonus = self.ucb_exploration * np.sqrt(
+                        np.log(max(total_obs, 1)) / arm.trials
+                    )
+                    values.append(mean + bonus)
+            arm_idx = int(np.argmax(values))
+            logger.debug(
+                f"UCB fallback: role={role}, total_obs={total_obs}, "
+                f"selected_arm={arm_idx}"
+            )
+
+        bid_amount = round(self.bid_arms[arm_idx] * self.total_budget)
+
+        budget_rem = _get_attr(auction_state, "budget_remaining", self.total_budget)
+        if bid_amount > budget_rem:
+            bid_amount = budget_rem
+            arm_idx = int(np.argmin(np.abs(self.bid_arms * self.total_budget - bid_amount)))
+
+        return arm_idx, bid_amount
+
+    def update(self, arm_idx: int, reward: float, player_role: str):
+        """Update Beta posterior for the selected arm.
+
+        Reward should be a normalized value: (player_season_value - cost) scaled.
+
+        Args:
+            arm_idx: index of the selected arm.
+            reward: normalized reward signal (higher = better purchase).
+            player_role: role of the purchased player.
+        """
+        reward_clipped = max(0.0, min(1.0, reward))
+        win = 1.0 if reward > 0 else 0.0
+
+        for key, arms in self.posteriors.items():
+            role, _scarcity = key
+            if role == player_role:
+                if arm_idx < len(arms):
+                    arm = arms[arm_idx]
+                    arm.trials += 1
+                    arm.wins += win
+                    arm.alpha += reward_clipped
+                    arm.beta += (1.0 - reward_clipped)
+                    logger.debug(
+                        f"Updated arm {arm_idx} for {player_role}: "
+                        f"trials={arm.trials}, alpha={arm.alpha:.2f}, beta={arm.beta:.2f}"
+                    )
+        for key, arms in self.posteriors.items():
+            role, _scarcity = key
+            if role == player_role:
+                for i, arm in enumerate(arms):
+                    if i == arm_idx:
+                        continue
+                    arm.beta = max(arm.beta, 1.001)
+
+    def get_arm_stats(self) -> Dict[str, dict]:
+        """Return arm statistics for analysis.
+
+        Returns:
+            dict mapping "role/scarcity/arm_idx" -> stats dict.
+        """
+        stats = {}
+        agg = {}
+        for (role, scarcity), arms in self.posteriors.items():
+            for idx, arm in enumerate(arms):
+                label = f"{role}/scarcity={scarcity}/arm={idx}"
+                stats[label] = {
+                    "trials": arm.trials,
+                    "wins": arm.wins,
+                    "alpha": arm.alpha,
+                    "beta": arm.beta,
+                    "win_rate": arm.wins / max(arm.trials, 1),
+                    "bid_fraction": float(self.bid_arms[idx]),
+                    "bid_amount": float(self.bid_arms[idx] * self.total_budget),
+                }
+                if idx not in agg:
+                    agg[idx] = {"trials": 0, "wins": 0.0, "alpha": 0.0, "beta": 0.0}
+                agg[idx]["trials"] += arm.trials
+                agg[idx]["wins"] += arm.wins
+                agg[idx]["alpha"] += arm.alpha
+                agg[idx]["beta"] += arm.beta
+        for idx, a in agg.items():
+            stats[idx] = {
+                "trials": a["trials"],
+                "wins": a["wins"],
+                "alpha": a["alpha"],
+                "beta": a["beta"],
+                "win_rate": a["wins"] / max(a["trials"], 1),
+                "bid_fraction": float(self.bid_arms[idx]),
+                "bid_amount": float(self.bid_arms[idx] * self.total_budget),
+            }
+        return stats
+
+    def exploration_bonus(self, player_role: str, player_points: float = 0.0) -> float:
+        """Compute exploration bonus for new/unknown player types.
+
+        Higher bonus when the bandit has limited experience with a role.
+        Encourages exploring sleeper players.
+
+        Returns:
+            Recommended extra bid amount in credits.
+        """
+        total_obs = 0
+        for (role, _scarcity), arms in self.posteriors.items():
+            if role == player_role:
+                total_obs += sum(a.trials for a in arms)
+
+        if total_obs == 0:
+            bonus = self.total_budget * 0.04
+        elif total_obs < 20:
+            bonus = self.total_budget * 0.025
+        elif total_obs < 50:
+            bonus = self.total_budget * 0.01
+        else:
+            bonus = 0.0
+
+        if bonus > 0:
+            base = max(player_points * 0.5, 0)
+            bonus += base * (1.0 / max(total_obs, 1)) * 20
+
+        logger.debug(f"Exploration bonus for {player_role}: {bonus:.1f} (obs={total_obs})")
+        return bonus
+
+    def recommend_bid_summary(
+        self,
+        player,
+        auction_state,
+        pool_stats: Optional[dict] = None,
+    ) -> dict:
+        """Full bidding recommendation for a player.
+
+        Returns a dict with arm index, bid amount, context, exploration bonus,
+        and total recommended bid.
+        """
+        arm_idx, bid_amount = self.select_bid(player, auction_state, pool_stats)
+        ctx = self.compute_context(player, auction_state, pool_stats)
+        bonus = self.exploration_bonus(ctx.role, ctx.projected_points)
+
+        return {
+            "player": _get_attr(player, "name", "unknown"),
+            "role": ctx.role,
+            "arm_index": arm_idx,
+            "base_bid": bid_amount,
+            "exploration_bonus": round(bonus),
+            "total_bid": round(bid_amount + bonus),
+            "recommended_bid": round(bid_amount + bonus),
+            "context": {
+                "points_zscore": round(
+                    (ctx.projected_points - 6.5) / 2.0, 2
+                ),
+                "role_scarcity": round(ctx.role_scarcity, 3),
+                "budget_remaining_frac": round(ctx.budget_remaining_fraction, 3),
+                "slots_remaining": ctx.slots_remaining_in_role,
+                "round": ctx.round_number,
+                "vor": round(ctx.value_over_replacement, 1),
+            },
+        }

src/optimization/budget_optimizer.py

@@ -0,0 +1,538 @@
+"""Bayesian Optimization for role-level budget allocation.
+
+Uses Gaussian Process regression to find optimal budget distribution
+across roles (P, D, C, A) that maximizes total projected team value.
+"""
+
+import logging
+from dataclasses import dataclass, field
+from typing import Dict, List, Optional, Tuple
+
+import numpy as np
+import pandas as pd
+from scipy.optimize import minimize
+from scipy.special import softmax
+
+from src.optimization.auction_solver import AuctionConfig
+
+logger = logging.getLogger(__name__)
+
+DEFAULT_ROSTER_QUOTAS = {"P": 3, "D": 8, "C": 8, "A": 6}
+
+
+@dataclass
+class RoleBudgetResult:
+    allocation: Dict[str, float]
+    total_value: float
+    role_values: Dict[str, float]
+    value_curves: Dict[str, Tuple[np.ndarray, np.ndarray]]
+
+
+class BudgetOptimizer:
+    """Bayesian Optimization for budget allocation across roster roles.
+
+    Finds the split of total_budget across P/D/C/A that yields the
+    highest possible team points via greedy fill within each role's budget.
+    Supports mid-auction adaptive rebalancing.
+    """
+
+    def __init__(
+        self,
+        total_budget: float = 500.0,
+        roster_quotas: Optional[Dict[str, int]] = None,
+        auction_config: Optional[AuctionConfig] = None,
+        random_state: int = 42,
+    ):
+        self.total_budget = total_budget
+        self.roster_quotas = roster_quotas or dict(DEFAULT_ROSTER_QUOTAS)
+        self.auction_config = auction_config or AuctionConfig()
+        self.random_state = random_state
+        self.rng = np.random.RandomState(random_state)
+        self._last_allocation: Optional[Dict[str, float]] = None
+        self._last_value: float = 0.0
+        self._optimization_history: List[dict] = []
+        self._gk_available = False
+
+    # ------------------------------------------------------------------
+    # Core optimization
+    # ------------------------------------------------------------------
+
+    def optimize(
+        self,
+        player_pool_df: pd.DataFrame,
+        n_calls: int = 50,
+        use_skopt: bool = True,
+    ) -> Dict[str, float]:
+        """Optimize budget allocation across roles.
+
+        Uses scikit-optimize GaussianProcessRegressor if available,
+        otherwise simplex-based local search.
+
+        Args:
+            player_pool_df: DataFrame with columns [name, role, projected_points].
+            n_calls: number of GP evaluations.
+            use_skopt: attempt Gaussian Process optimization.
+
+        Returns:
+            Dict mapping role -> recommended budget amount.
+        """
+        if "role" not in player_pool_df.columns or "projected_points" not in player_pool_df.columns:
+            raise ValueError("player_pool_df must have 'role' and 'projected_points' columns")
+
+        roles = list(self.roster_quotas.keys())
+        n_roles = len(roles)
+
+        if use_skopt:
+            try:
+                return self._optimize_gp(player_pool_df, roles, n_calls)
+            except ImportError:
+                logger.info("scikit-optimize not installed. Using local search.")
+            except Exception as exc:
+                logger.warning(f"GP optimization failed: {exc}. Using local search.")
+
+        return self._optimize_local(player_pool_df, roles, n_calls)
+
+    def _optimize_gp(
+        self, player_pool_df: pd.DataFrame, roles: list, n_calls: int
+    ) -> Dict[str, float]:
+        """Gaussian Process-based budget optimization."""
+        from skopt import gp_minimize
+        from skopt.space import Space
+        from skopt.learning import GaussianProcessRegressor
+
+        n_roles = len(roles)
+        space = Space([(0.01, 0.70) for _ in range(n_roles)])
+
+        def objective_wrapper(fractions):
+            fractions = np.array(fractions, dtype=float)
+            fractions = self._normalize_fractions(fractions)
+            value = self._objective(fractions, roles, player_pool_df)
+            self._optimization_history.append({
+                "fractions": fractions.tolist(),
+                "value": value,
+            })
+            return -value
+
+        def params_to_fractions(params):
+            return self._normalize_fractions(np.array(params, dtype=float))
+
+        result = gp_minimize(
+            objective_wrapper,
+            space,
+            n_calls=n_calls,
+            random_state=self.random_state,
+            n_initial_points=max(10, n_calls // 5),
+            verbose=False,
+            n_jobs=-1,
+        )
+
+        best_fractions = params_to_fractions(result.x)
+        self._last_allocation = self._fractions_to_allocation(best_fractions, roles)
+        self._last_value = -result.fun
+
+        logger.info(
+            f"GP Budget optimization complete: {self._last_allocation} "
+            f"=> value={self._last_value:.1f}"
+        )
+
+        return self._last_allocation
+
+    def _optimize_local(
+        self, player_pool_df: pd.DataFrame, roles: list, n_calls: int
+    ) -> Dict[str, float]:
+        """Local search optimization using Nelder-Mead simplex."""
+        n_roles = len(roles)
+
+        best_allocation = None
+        best_value = -float("inf")
+
+        for restart in range(max(n_calls // 10, 1)):
+            x0 = self._random_allocation(n_roles)
+
+            def simplex_objective(fractions):
+                fractions = self._normalize_fractions(np.array(fractions, dtype=float))
+                value = self._objective(fractions, roles, player_pool_df)
+                self._optimization_history.append({
+                    "fractions": fractions.tolist(),
+                    "value": value,
+                })
+                return -value
+
+            res = minimize(
+                simplex_objective,
+                x0,
+                method="Nelder-Mead",
+                options={"maxiter": max(n_calls // 3, 20), "xatol": 1e-3, "fatol": 1e-3},
+            )
+
+            fractions = self._normalize_fractions(np.array(res.x, dtype=float))
+            value = -res.fun
+
+            if value > best_value:
+                best_value = value
+                best_allocation = fractions
+
+        if best_allocation is None:
+            best_allocation = self._proportional_allocation(player_pool_df, roles)
+
+        self._last_allocation = self._fractions_to_allocation(best_allocation, roles)
+        self._last_value = best_value
+
+        logger.info(
+            f"Local optimization complete: {self._last_allocation} "
+            f"=> value={self._last_value:.1f}"
+        )
+
+        return self._last_allocation
+
+    # ------------------------------------------------------------------
+    # Adaptive rebalancing
+    # ------------------------------------------------------------------
+
+    def optimize_adaptive(
+        self,
+        player_pool_df: pd.DataFrame,
+        remaining_slots: Dict[str, int],
+        spent_per_role: Dict[str, float],
+        n_calls: int = 30,
+    ) -> Dict[str, float]:
+        """Mid-auction rebalancing — optimize remaining budget for unfilled slots.
+
+        Args:
+            player_pool_df: remaining available players pool.
+            remaining_slots: dict of role -> slots still needed.
+            spent_per_role: dict of role -> credits already spent.
+            n_calls: GP evaluation budget.
+
+        Returns:
+            Dict mapping role -> recommended budget for remaining slots.
+        """
+        remaining_budget = self.total_budget - sum(spent_per_role.values())
+        remaining_budget = max(1.0, remaining_budget)
+
+        if sum(remaining_slots.values()) == 0:
+            logger.info("All slots filled. No budget to allocate.")
+            return {r: 0.0 for r in self.roster_quotas}
+
+        saved_quotas = self.roster_quotas
+        saved_total = self.total_budget
+        self.roster_quotas = dict(remaining_slots)
+        self.total_budget = remaining_budget
+
+        available = player_pool_df[
+            player_pool_df["role"].isin(
+                [r for r, s in remaining_slots.items() if s > 0]
+            )
+        ]
+
+        if len(available) == 0:
+            logger.warning("No available players for remaining slots.")
+            self.roster_quotas = saved_quotas
+            self.total_budget = saved_total
+            return {r: 0.0 for r in saved_quotas}
+
+        allocation = self.optimize(
+            available,
+            n_calls=n_calls,
+            use_skopt=True,
+        )
+
+        self.roster_quotas = saved_quotas
+        self.total_budget = saved_total
+
+        result = {}
+        for role in saved_quotas:
+            result[role] = allocation.get(role, 0.0)
+        return result
+
+    # ------------------------------------------------------------------
+    # Objective function
+    # ------------------------------------------------------------------
+
+    def _objective(
+        self,
+        fractions: np.ndarray,
+        roles: list,
+        player_pool_df: pd.DataFrame,
+    ) -> float:
+        """Simulate greedy fill within role budgets; return total projected points.
+
+        For each role, pick the best players by projected_points until the
+        role budget or slot quota is exhausted.
+        """
+        total_value = 0.0
+
+        for i, role in enumerate(roles):
+            budget = fractions[i] * self.total_budget
+            slots = self.roster_quotas.get(role, 0)
+            role_players = player_pool_df[player_pool_df["role"] == role].copy()
+
+            if len(role_players) == 0 or slots == 0:
+                continue
+
+            role_players = role_players.sort_values(
+                "projected_points", ascending=False
+            )
+
+            total_cost = 0.0
+            filled = 0
+
+            for _, player in role_players.iterrows():
+                points = float(player["projected_points"])
+                estimated_price = self._estimate_price_simple(points, budget, role)
+
+                if total_cost + estimated_price > budget:
+                    continue
+
+                total_cost += estimated_price
+                total_value += points
+                filled += 1
+
+                if filled >= slots:
+                    break
+
+        return total_value
+
+    def _estimate_price_simple(
+        self, projected_points: float, role_budget: float, role: str
+    ) -> float:
+        """Simple price estimate: points * role_factor clamped within budget."""
+        role_factor = {"P": 4.0, "D": 2.5, "C": 3.0, "A": 4.5}.get(role, 3.0)
+        price = projected_points * role_factor
+        max_price = role_budget * 0.50
+        return min(price, max_price, role_budget)
+
+    # ------------------------------------------------------------------
+    # Allocation access and visualization data
+    # ------------------------------------------------------------------
+
+    def get_allocation(self) -> Dict[str, float]:
+        """Return last computed allocation (role -> budget amount)."""
+        if self._last_allocation is None:
+            return {r: self.total_budget / len(self.roster_quotas) for r in self.roster_quotas}
+        return dict(self._last_allocation)
+
+    def get_role_value_curves(
+        self,
+        player_pool_df: pd.DataFrame,
+        n_points: int = 20,
+    ) -> Dict[str, Tuple[np.ndarray, np.ndarray]]:
+        """Compute diminishing returns curves: budget vs. expected points per role.
+
+        Returns:
+            dict role -> (budget_array, value_array).
+        """
+        curves = {}
+        budget_step = self.total_budget / n_points
+
+        for role, slots in self.roster_quotas.items():
+            role_players = player_pool_df[player_pool_df["role"] == role].sort_values(
+                "projected_points", ascending=False
+            )
+
+            budgets = np.linspace(0, self.total_budget, n_points)
+            values = np.zeros(n_points)
+
+            for i, budget_limit in enumerate(budgets):
+                total_cost = 0.0
+                total_value = 0.0
+                filled = 0
+
+                for _, player in role_players.iterrows():
+                    points = float(player["projected_points"])
+                    price = self._estimate_price_simple(points, budget_limit, role)
+
+                    if total_cost + price > budget_limit:
+                        continue
+
+                    total_cost += price
+                    total_value += points
+                    filled += 1
+
+                    if filled >= slots:
+                        break
+
+                values[i] = total_value
+
+            curves[role] = (budgets.copy(), values.copy())
+
+        return curves
+
+    # ------------------------------------------------------------------
+    # Fallback: proportional allocation
+    # ------------------------------------------------------------------
+
+    def _proportional_allocation(
+        self, player_pool_df: pd.DataFrame, roles: list
+    ) -> np.ndarray:
+        """Allocate budget proportional to (points_variance * slots) per role."""
+        weights = np.zeros(len(roles))
+
+        for i, role in enumerate(roles):
+            role_players = player_pool_df[player_pool_df["role"] == role]
+            if len(role_players) > 1:
+                variance = role_players["projected_points"].var()
+            else:
+                variance = 1.0
+            slots = self.roster_quotas.get(role, 1)
+            weights[i] = variance * slots
+
+        weight_sum = weights.sum()
+        if weight_sum <= 0:
+            return np.ones(len(roles)) / len(roles)
+
+        fractions = weights / weight_sum
+        fractions = np.clip(fractions, 0.02, 0.70)
+        fractions = fractions / fractions.sum()
+
+        logger.info(
+            f"Proportional allocation (fallback): "
+            f"{dict(zip(roles, fractions.round(3)))}"
+        )
+        return fractions
+
+    # ------------------------------------------------------------------
+    # Utilities
+    # ------------------------------------------------------------------
+
+    def _normalize_fractions(self, fractions: np.ndarray) -> np.ndarray:
+        """Normalize fractions to sum to 1.0 with minimum per role."""
+        fractions = np.clip(fractions, 0.01, 0.70)
+        total = fractions.sum()
+        if total <= 0:
+            return np.ones_like(fractions) / len(fractions)
+        return fractions / total
+
+    def _fractions_to_allocation(
+        self, fractions: np.ndarray, roles: list
+    ) -> Dict[str, float]:
+        """Convert fractions to absolute budget per role."""
+        return {
+            role: round(float(fractions[i] * self.total_budget), 1)
+            for i, role in enumerate(roles)
+        }
+
+    def _random_allocation(self, n_roles: int) -> np.ndarray:
+        """Generate a random allocation via Dirichlet."""
+        alpha = np.ones(n_roles) * 2.0
+        return self.rng.dirichlet(alpha)
+
+    def get_optimization_trace(self) -> pd.DataFrame:
+        """Return DataFrame of all evaluated allocations during optimization."""
+        if not self._optimization_history:
+            return pd.DataFrame()
+        return pd.DataFrame(self._optimization_history)
+
+    def estimate_team_composition(
+        self,
+        player_pool_df: pd.DataFrame,
+        allocation: Optional[Dict[str, float]] = None,
+    ) -> pd.DataFrame:
+        """Given final allocation, return the recommended player selections.
+
+        Returns:
+            DataFrame with selected players, their estimated costs, and value.
+        """
+        roles = list(self.roster_quotas.keys())
+        alloc = allocation or self._last_allocation
+        if alloc is None:
+            alloc = self._proportional_allocation_as_dict(roles)
+
+        selections = []
+
+        for role in roles:
+            budget = alloc.get(role, 0.0)
+            slots = self.roster_quotas.get(role, 0)
+            role_players = player_pool_df[player_pool_df["role"] == role].sort_values(
+                "projected_points", ascending=False
+            )
+
+            total_cost = 0.0
+            filled = 0
+
+            for _, player in role_players.iterrows():
+                points = float(player["projected_points"])
+                price = self._estimate_price_simple(points, budget, role)
+
+                if total_cost + price > budget:
+                    continue
+
+                total_cost += price
+                name = player.get("name", player.get("player_name", "unknown"))
+                selections.append({
+                    "player": name,
+                    "role": role,
+                    "projected_points": points,
+                    "estimated_cost": round(price, 1),
+                    "value_ratio": round(points / max(price, 1), 3),
+                })
+                filled += 1
+
+                if filled >= slots:
+                    break
+
+        return pd.DataFrame(selections)
+
+    def _proportional_allocation_as_dict(self, roles: list) -> Dict[str, float]:
+        fractions = self._proportional_allocation(
+            pd.DataFrame(columns=["role", "projected_points"]), roles
+        )
+        return {roles[i]: round(float(fractions[i] * self.total_budget), 1) for i in range(len(roles))}
+
+    # ------------------------------------------------------------------
+    # Sensitivity analysis
+    # ------------------------------------------------------------------
+
+    def sensitivity_analysis(
+        self,
+        player_pool_df: pd.DataFrame,
+        role: Optional[str] = None,
+        delta_pct: float = 0.05,
+        n_steps: int = 11,
+        n_calls: int = 50,
+    ) -> dict:
+        """Test how shifting budget into/out of one role affects total value.
+
+        Args:
+            player_pool_df: current player pool.
+            role: role to perturb. If None, tests all roles.
+            delta_pct: fractional step size.
+            n_steps: number of steps in each direction.
+            n_calls: number of optimization calls (for compatibility).
+
+        Returns:
+            Dict with at least a 'per_role' key containing per-role analysis.
+        """
+        base_allocation = self.get_allocation()
+        roles_to_test = [role] if role else list(base_allocation.keys())
+        per_role = {}
+
+        for test_role in roles_to_test:
+            results = []
+            shifts = np.linspace(-delta_pct * n_steps, delta_pct * n_steps, 2 * n_steps + 1)
+
+            for shift in shifts:
+                adjusted = {}
+                for r, val in base_allocation.items():
+                    adjusted[r] = val * (1.0 + (shift if r == test_role else -shift / 3.0))
+
+                total_adj = sum(adjusted.values())
+                for r in adjusted:
+                    adjusted[r] = adjusted[r] / total_adj * self.total_budget
+
+                fractions = np.array([adjusted[r] / self.total_budget for r in base_allocation.keys()])
+                value = self._objective(
+                    fractions,
+                    list(base_allocation.keys()),
+                    player_pool_df,
+                )
+
+                results.append({
+                    "shift_pct": round(shift * 100, 1),
+                    "allocation": {r: round(v, 1) for r, v in adjusted.items()},
+                    "total_value": round(value, 1),
+                })
+
+            per_role[test_role] = pd.DataFrame(results)
+
+        return {"per_role": per_role, "base_allocation": base_allocation}

src/optimization/opponent_bidding_model.py

@@ -0,0 +1,412 @@
+"""Opponent bidding behavior modeling using LightGBM.
+
+Predicts what competitors will bid for each player in a live auction round.
+Supports Monte Carlo simulation of auction outcomes and win probability estimates.
+"""
+
+import logging
+from dataclasses import dataclass, field
+from typing import Dict, List, Optional, Tuple
+
+import numpy as np
+import pandas as pd
+
+logger = logging.getLogger(__name__)
+
+
+@dataclass
+class OpponentState:
+    budget_remaining: float = 500.0
+    initial_budget: float = 500.0
+    slots_filled: Dict[str, int] = field(default_factory=lambda: {"P": 0, "D": 0, "C": 0, "A": 0})
+    slots_total: Dict[str, int] = field(default_factory=lambda: {"P": 3, "D": 8, "C": 8, "A": 6})
+    round_number: int = 1
+    aggression_factor: float = 1.0
+
+
+def _get_attr(obj, key, default=None):
+    if isinstance(obj, dict):
+        return obj.get(key, default)
+    return getattr(obj, key, default)
+
+
+class OpponentBidModel:
+    """Predicts opponent bids using LightGBM with heuristic fallback.
+
+    Trains on historical auction logs and outputs estimated max opponent
+    bid, win probability per player, and Monte Carlo round simulations.
+    """
+
+    def __init__(self, random_state: int = 42):
+        self.random_state = random_state
+        self.model = None
+        self.fitted = False
+        self.feature_names: list = []
+        self.rng = np.random.RandomState(random_state)
+        self._role_scarcity_cache: Dict[str, float] = {}
+
+    # ------------------------------------------------------------------
+    # Feature engineering
+    # ------------------------------------------------------------------
+
+    def _extract_features(
+        self,
+        players_df: pd.DataFrame,
+        opponent_state: OpponentState,
+    ) -> pd.DataFrame:
+        """Build feature matrix for LightGBM prediction.
+
+        Args:
+            players_df: DataFrame with columns [player_name, player_role,
+                        player_projected_points, ...].
+            opponent_state: OpponentState describing current opponent.
+
+        Returns:
+            Feature DataFrame ready for model input.
+        """
+        df = players_df.copy()
+
+        df["role_P"] = (df["player_role"] == "P").astype(float)
+        df["role_D"] = (df["player_role"] == "D").astype(float)
+        df["role_C"] = (df["player_role"] == "C").astype(float)
+        df["role_A"] = (df["player_role"] == "A").astype(float)
+
+        df["budget_remaining_frac"] = (
+            opponent_state.budget_remaining / opponent_state.initial_budget
+        )
+
+        for role in ["P", "D", "C", "A"]:
+            slots_total = opponent_state.slots_total.get(role, 1)
+            slots_filled = opponent_state.slots_filled.get(role, 0)
+            scarcity = (slots_total - slots_filled) / slots_total
+            df[f"scarcity_{role}"] = scarcity
+
+        role_map = {"P": 0, "D": 1, "C": 2, "A": 3}
+        df["role_code"] = df["player_role"].map(role_map)
+
+        if "player_projected_points" in df.columns:
+            df["points_sq"] = df["player_projected_points"] ** 2
+            df["points_log"] = np.log1p(df["player_projected_points"].clip(lower=0))
+
+        df["slots_needed_total"] = sum(
+            opponent_state.slots_total.get(r, 0) - opponent_state.slots_filled.get(r, 0)
+            for r in ["P", "D", "C", "A"]
+        )
+        df["slots_needed_total"] = df["slots_needed_total"].clip(lower=1)
+
+        df["round_number"] = opponent_state.round_number
+
+        df["urgency"] = 1.0 - (opponent_state.budget_remaining / opponent_state.initial_budget)
+        df["aggression"] = opponent_state.aggression_factor
+
+        self.feature_names = [
+            "player_projected_points",
+            "role_P",
+            "role_D",
+            "role_C",
+            "role_A",
+            "role_code",
+            "budget_remaining_frac",
+            "scarcity_P",
+            "scarcity_D",
+            "scarcity_C",
+            "scarcity_A",
+            "points_sq",
+            "points_log",
+            "slots_needed_total",
+            "round_number",
+            "urgency",
+            "aggression",
+        ]
+
+        for col in self.feature_names:
+            if col not in df.columns:
+                df[col] = 0.0
+
+        return df[self.feature_names]
+
+    # ------------------------------------------------------------------
+    # Training
+    # ------------------------------------------------------------------
+
+    def fit(self, auction_logs: pd.DataFrame):
+        """Train LightGBM regressor on historical auction logs.
+
+        Args:
+            auction_logs: DataFrame with columns [player_name, player_role,
+                player_projected_points, opponent_budget_remaining,
+                opponent_slots_remaining, role_needed_count,
+                round_number, winning_bid].
+        """
+        required_cols = [
+            "player_role", "player_projected_points",
+            "winning_bid",
+        ]
+        for col in required_cols:
+            if col not in auction_logs.columns:
+                raise ValueError(f"Missing required column: '{col}' in auction_logs")
+
+        df = auction_logs.dropna(subset=required_cols).copy()
+
+        if len(df) < 20:
+            logger.warning(
+                f"Only {len(df)} auction records. Not enough to fit LightGBM. "
+                "Using heuristic fallback."
+            )
+            self.fitted = False
+            return
+
+        try:
+            import lightgbm as lgb
+        except ImportError:
+            logger.warning("LightGBM not available. Using heuristic fallback.")
+            self.fitted = False
+            return
+
+        dummy_state = OpponentState(
+            budget_remaining=df.get("opponent_budget_remaining", 500),
+            initial_budget=500,
+            slots_filled={"P": 0, "D": 0, "C": 0, "A": 0},
+            slots_total={"P": 3, "D": 8, "C": 8, "A": 6},
+            round_number=1,
+            aggression_factor=1.0,
+        )
+
+        df = df.rename(columns={
+            "player_role": "player_role",
+            "player_projected_points": "player_projected_points",
+        })
+
+        dummy_df = df[["player_role", "player_projected_points"]].copy()
+        dummy_df.columns = ["player_role", "player_projected_points"]
+
+        X = self._extract_features(dummy_df, dummy_state)
+
+        y = df["winning_bid"].astype(float)
+        y_min, y_max = y.min(), y.max()
+
+        self.model = lgb.LGBMRegressor(
+            n_estimators=100,
+            max_depth=6,
+            learning_rate=0.05,
+            num_leaves=31,
+            min_child_samples=10,
+            subsample=0.8,
+            colsample_bytree=0.8,
+            random_state=self.random_state,
+            verbose=-1,
+        )
+        self.model.fit(X, y)
+        self.fitted = True
+        self._y_min = y_min
+        self._y_max = y_max
+
+        logger.info(
+            f"OpponentBidModel trained on {len(df)} records. "
+            f"Target range: [{y_min:.0f}, {y_max:.0f}]"
+        )
+
+    # ------------------------------------------------------------------
+    # Prediction
+    # ------------------------------------------------------------------
+
+    def predict_opponent_bids(
+        self,
+        players_df: pd.DataFrame,
+        opponent_state: OpponentState,
+    ) -> pd.Series:
+        """Predict max opponent bid for each player.
+
+        Args:
+            players_df: DataFrame with player info.
+            opponent_state: current opponent state.
+
+        Returns:
+            Series of predicted max opponent bids (index = player index).
+        """
+        if not self.fitted or self.model is None:
+            bids = self._heuristic_bid(players_df, opponent_state)
+            return pd.Series(bids, index=players_df.index)
+
+        X = self._extract_features(players_df, opponent_state)
+        predictions = self.model.predict(X)
+        predictions = np.clip(predictions, 1, _get_attr(opponent_state, "budget_remaining", 500))
+
+        return pd.Series(predictions, index=players_df.index)
+
+    def predict_p_acquire(
+        self,
+        players_df: pd.DataFrame,
+        my_bids: np.ndarray,
+        opponent_state: OpponentState,
+        temperature: float = 0.1,
+    ) -> np.ndarray:
+        """Probability I acquire each player given my bids vs opponent.
+
+        Uses sigmoid: P = 1 / (1 + exp(-(my_bid - opp_bid) / temperature)).
+
+        Args:
+            players_df: DataFrame with player info.
+            my_bids: array of my bid amounts per player.
+            opponent_state: current opponent state.
+            temperature: softmax temperature (lower = sharper).
+
+        Returns:
+            Array of acquisition probabilities per player.
+        """
+        opp_bids = self.predict_opponent_bids(players_df, opponent_state).values
+        margin = np.array(my_bids, dtype=float) - opp_bids
+        scaled_temp = max(temperature * max(opp_bids.max(), 1), 0.01)
+        probabilities = 1.0 / (1.0 + np.exp(-margin / scaled_temp))
+        return np.clip(probabilities, 0.01, 0.99)
+
+    # ------------------------------------------------------------------
+    # Monte Carlo round simulation
+    # ------------------------------------------------------------------
+
+    def simulate_live_round(
+        self,
+        available_players: pd.DataFrame,
+        my_budget: float,
+        opponent_state: OpponentState,
+        n_sims: int = 1000,
+    ) -> dict:
+        """Monte Carlo simulation of a live auction round.
+
+        Args:
+            available_players: DataFrame of players up for bidding this round.
+            my_budget: my remaining budget.
+            opponent_state: opponent's current state.
+            n_sims: number of simulation runs.
+
+        Returns:
+            dict with expected_players_acquired, expected_cost, value_matrix.
+        """
+        opp_bids = self.predict_opponent_bids(available_players, opponent_state).values
+
+        n_players = len(available_players)
+        players_acquired = np.zeros(n_sims, dtype=int)
+        total_cost = np.zeros(n_sims, dtype=float)
+        value_matrix = np.zeros((n_sims, n_players), dtype=float)
+
+        for sim_idx in range(n_sims):
+            my_budget_left = my_budget
+            acquired = 0
+            cost = 0.0
+
+            for p_idx in range(n_players):
+                opp_bid = opp_bids[p_idx] + self.rng.normal(0, max(opp_bids[p_idx] * 0.15, 1))
+                opp_bid = max(opp_bid, 1)
+
+                my_bid = self._heuristic_bid_single(
+                    available_players.iloc[p_idx], opponent_state
+                )
+
+                my_bid = min(my_bid, my_budget_left)
+                if my_bid > opp_bid:
+                    acquired += 1
+                    cost += my_bid
+                    my_budget_left -= my_bid
+                    value_matrix[sim_idx, p_idx] = 1.0
+
+            players_acquired[sim_idx] = acquired
+            total_cost[sim_idx] = cost
+
+        return {
+            "expected_players_acquired": float(np.mean(players_acquired)),
+            "expected_cost": float(np.mean(total_cost)),
+            "cost_std": float(np.std(total_cost)),
+            "acquired_std": float(np.std(players_acquired)),
+            "cost_percentile_25": float(np.percentile(total_cost, 25)),
+            "cost_percentile_50": float(np.percentile(total_cost, 50)),
+            "cost_percentile_75": float(np.percentile(total_cost, 75)),
+            "acquisition_rate": float(players_acquired.mean() / n_players),
+            "value_matrix": value_matrix,
+            "n_sims": n_sims,
+        }
+
+    # ------------------------------------------------------------------
+    # Heuristic fallback
+    # ------------------------------------------------------------------
+
+    def _heuristic_bid_single(
+        self, player_row: pd.Series, opponent_state: OpponentState
+    ) -> float:
+        """Heuristic bid for a single player (scalar version)."""
+        role = player_row.get("player_role", None)
+        points = float(player_row.get("player_projected_points", 6.5))
+
+        if isinstance(role, pd.Series):
+            role = role.iloc[0]
+
+        scarcity = self._compute_role_scarcity(role, opponent_state)
+        aggression = _get_attr(opponent_state, "aggression_factor", 1.0)
+        budget_rem = _get_attr(opponent_state, "budget_remaining", 500.0)
+        budget_init = _get_attr(opponent_state, "initial_budget", 500.0) or _get_attr(opponent_state, "total_budget", 500.0)
+        base_bid = 0.4 * points * scarcity * aggression
+        bid = base_bid * (budget_rem / budget_init)
+        return max(bid, 1.0)
+
+    def _heuristic_bid(
+        self, players_df: pd.DataFrame, opponent_state: OpponentState
+    ) -> np.ndarray:
+        """Heuristic bid array for all players."""
+        bids = []
+        for _, row in players_df.iterrows():
+            bids.append(self._heuristic_bid_single(row, opponent_state))
+        return np.array(bids, dtype=float)
+
+    def _compute_role_scarcity(
+        self, role: Optional[str], opponent_state
+    ) -> float:
+        """Compute how scarce a role is for the opponent."""
+        slots_total_dict = _get_attr(opponent_state, "slots_total", {"P": 3, "D": 8, "C": 8, "A": 6})
+        slots_filled_dict = _get_attr(opponent_state, "slots_filled", {"P": 0, "D": 0, "C": 0, "A": 0})
+        if role is None or role not in slots_total_dict:
+            return 1.0
+        slots_total = slots_total_dict[role]
+        filled = slots_filled_dict.get(role, 0)
+        remaining = max(slots_total - filled, 1)
+        return slots_total / remaining
+
+    # ------------------------------------------------------------------
+    # Batch simulation with opponent model integration
+    # ------------------------------------------------------------------
+
+    def run_auction_simulation(
+        self,
+        player_groups: List[pd.DataFrame],
+        initial_budgets: List[float],
+        opponent_states: List[OpponentState],
+        n_sims: int = 500,
+    ) -> dict:
+        """Simulate multi-round auction against multiple opponents.
+
+        Args:
+            player_groups: list of DataFrames, one per round, with available players.
+            initial_budgets: my starting budget per round/concept.
+            opponent_states: OpponentState for each round.
+            n_sims: number of Monte Carlo runs.
+
+        Returns:
+            dict with aggregated simulation results.
+        """
+        all_results = []
+        for idx, (players, budget, opp_state) in enumerate(
+            zip(player_groups, initial_budgets, opponent_states)
+        ):
+            result = self.simulate_live_round(
+                players, budget, opp_state, n_sims=n_sims
+            )
+            result["round"] = idx
+            all_results.append(result)
+
+        total_acquired = sum(r["expected_players_acquired"] for r in all_results)
+        total_cost = sum(r["expected_cost"] for r in all_results)
+
+        return {
+            "rounds": all_results,
+            "total_expected_acquired": total_acquired,
+            "total_expected_cost": total_cost,
+            "n_sims": n_sims,
+        }