"""
Blend modello Dixon-Coles + mercato: il primo vero test di valore aggiunto.

Domanda: combinare le probabilità del modello con quelle del mercato produce
un log loss migliore del mercato da solo? Se sì, il modello contiene
informazione che il mercato non ha già incorporato.

Protocollo (senza look-ahead):
  - il peso del blend si impara sulla stagione di TRAIN (default 2024/25)
  - la valutazione avviene solo sulla stagione di TEST (default 2025/26)
  - bootstrap a coppie (1000 ricampionamenti) per la significatività

Due metodi di blend, il migliore viene scelto in base al TRAIN:
  - lineare:    p = w * modello + (1-w) * mercato
  - geometrico: p ∝ modello^w * mercato^(1-w)   (pooling in log-spazio)

Uso:
    python scripts/blend.py
    python scripts/blend.py --train-seasons 2024/25 --test-seasons 2025/26
"""

import argparse
import sys
from pathlib import Path

import numpy as np
import pandas as pd
from scipy.optimize import minimize_scalar

PROJECT_ROOT = Path(__file__).resolve().parent.parent
BT_PATH = PROJECT_ROOT / "data" / "processed" / "backtest_base.csv"
OUT_PATH = PROJECT_ROOT / "data" / "processed" / "backtest_blend.csv"

RESULT_INDEX = {"H": 0, "D": 1, "A": 2}


def normalize(p: np.ndarray) -> np.ndarray:
    p = np.clip(p, 1e-12, 1)
    return p / p.sum(axis=1, keepdims=True)


def linear_blend(p_model: np.ndarray, p_mkt: np.ndarray, w: float) -> np.ndarray:
    return normalize(w * p_model + (1 - w) * p_mkt)


def geometric_blend(p_model: np.ndarray, p_mkt: np.ndarray, w: float) -> np.ndarray:
    return normalize(np.exp(w * np.log(np.clip(p_model, 1e-12, 1))
                            + (1 - w) * np.log(np.clip(p_mkt, 1e-12, 1))))


def logloss(p: np.ndarray, outcome: np.ndarray) -> float:
    p = normalize(p)
    return float(-np.mean(np.log(p[np.arange(len(outcome)), outcome])))


def fit_weight(blend_fn, p_model, p_mkt, outcome) -> float:
    res = minimize_scalar(lambda w: logloss(blend_fn(p_model, p_mkt, w), outcome),
                          bounds=(0.0, 1.0), method="bounded")
    return float(res.x)


def bootstrap_diff(p_a: np.ndarray, p_b: np.ndarray, outcome: np.ndarray,
                   n_boot: int = 1000, seed: int = 42) -> tuple[float, float, float]:
    """
    Differenza di log loss (A - B) con intervallo di confidenza al 95%
    tramite bootstrap a coppie sulle stesse partite.
    """
    rng = np.random.default_rng(seed)
    n = len(outcome)
    ll_a = -np.log(normalize(p_a)[np.arange(n), outcome])
    ll_b = -np.log(normalize(p_b)[np.arange(n), outcome])
    diff = ll_a - ll_b
    boots = np.array([diff[rng.integers(0, n, n)].mean() for _ in range(n_boot)])
    return float(diff.mean()), float(np.quantile(boots, 0.025)), float(np.quantile(boots, 0.975))


def main() -> int:
    parser = argparse.ArgumentParser(description="Blend modello+mercato con pesi appresi")
    parser.add_argument("--train-seasons", nargs="+", default=["2024/25"])
    parser.add_argument("--test-seasons", nargs="+", default=["2025/26"])
    parser.add_argument("--input", default=str(BT_PATH),
                        help="CSV di backtest da valutare (default: backtest_base.csv)")
    args = parser.parse_args()

    in_path = Path(args.input)
    if not in_path.exists():
        print(f"Manca {in_path}: esegui prima scripts/backtest.py")
        return 1

    bt = pd.read_csv(in_path).dropna(subset=["mkt_h", "mkt_d", "mkt_a"])
    train = bt[bt["season"].isin(args.train_seasons)]
    test = bt[bt["season"].isin(args.test_seasons)]
    if train.empty or test.empty:
        print("Stagioni di train o test non trovate nel backtest.")
        return 1

    def unpack(df):
        return (df[["model_h", "model_d", "model_a"]].to_numpy(),
                df[["mkt_h", "mkt_d", "mkt_a"]].to_numpy(),
                df["result"].map(RESULT_INDEX).to_numpy())

    pm_tr, pk_tr, y_tr = unpack(train)
    pm_te, pk_te, y_te = unpack(test)
    print(f"Train: {len(train)} partite {args.train_seasons} | "
          f"Test: {len(test)} partite {args.test_seasons}")

    # Apprendimento pesi sul train e scelta del metodo migliore (sempre sul train)
    methods = {}
    for name, fn in [("lineare", linear_blend), ("geometrico", geometric_blend)]:
        w = fit_weight(fn, pm_tr, pk_tr, y_tr)
        methods[name] = (fn, w, logloss(fn(pm_tr, pk_tr, w), y_tr))
        print(f"  blend {name}: peso modello w = {w:.3f} "
              f"(log loss train {methods[name][2]:.4f})")
    best_name = min(methods, key=lambda k: methods[k][2])
    best_fn, best_w, _ = methods[best_name]
    print(f"  scelto (dal train): {best_name}")

    # Valutazione out-of-sample sul test
    p_blend = best_fn(pm_te, pk_te, best_w)
    print(f"\n=== Test {args.test_seasons} ({len(test)} partite) ===")
    print(f"{'':<22}{'log loss':>10}")
    print(f"{'Mercato (closing)':<22}{logloss(pk_te, y_te):>10.4f}")
    print(f"{'Modello DC':<22}{logloss(pm_te, y_te):>10.4f}")
    print(f"{'Blend ' + best_name:<22}{logloss(p_blend, y_te):>10.4f}")

    mean_diff, lo, hi = bootstrap_diff(p_blend, pk_te, y_te)
    verdict = "il blend BATTE il mercato" if hi < 0 else (
        "il mercato batte il blend" if lo > 0 else
        "differenza NON significativa (compatibile con il rumore)")
    print(f"\nDelta log loss blend - mercato: {mean_diff:+.4f} "
          f"[IC 95%: {lo:+.4f}, {hi:+.4f}]")
    print(f"Verdetto: {verdict}")

    print("\nPer lega (test):")
    test = test.copy()
    test[["blend_h", "blend_d", "blend_a"]] = p_blend
    for lg, grp in test.groupby("league"):
        y = grp["result"].map(RESULT_INDEX).to_numpy()
        ll_b = logloss(grp[["blend_h", "blend_d", "blend_a"]].to_numpy(), y)
        ll_k = logloss(grp[["mkt_h", "mkt_d", "mkt_a"]].to_numpy(), y)
        print(f"  {lg:<5} blend {ll_b:.4f}   mercato {ll_k:.4f}   gap {ll_b - ll_k:+.4f}")

    test.to_csv(OUT_PATH, index=False)
    print(f"\nPredizioni blend salvate in: {OUT_PATH}")
    return 0


if __name__ == "__main__":
    sys.exit(main())