Monte Carlo Simulation: How to go broke eight times faster with 8 Eurojackpot lines

"""
Eurojackpot Monte Carlo Simulation
===================================
Simulates realistic lottery participation over multiple years
using exact combinatorial probabilities and historical payout data.
"""

from __future__ import annotations

import math
import logging
from dataclasses import dataclass, field
from pathlib import Path
from typing import TYPE_CHECKING

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import numpy as np
import polars as pl

if TYPE_CHECKING:
    pass

# ---------------------------------------------------------------------------
# Logging
# ---------------------------------------------------------------------------
logging.basicConfig(level=logging.INFO, format="%(message)s")
log = logging.getLogger(__name__)


# ---------------------------------------------------------------------------
# Constants
# ---------------------------------------------------------------------------
MAIN_POOL = 50
MAIN_PICK = 5
EURO_POOL = 12
EURO_PICK = 2
TICKET_COST = 2.0
DRAWS_PER_YEAR = 104  # 2× per week

PRIZE_TIERS: dict[str, tuple[int, int]] = {
    "5 + 2": (5, 2), "5 + 1": (5, 1), "5 + 0": (5, 0),
    "4 + 2": (4, 2), "4 + 1": (4, 1), "4 + 0": (4, 0),
    "3 + 2": (3, 2), "2 + 2": (2, 2), "3 + 1": (3, 1),
    "3 + 0": (3, 0), "1 + 2": (1, 2), "2 + 1": (2, 1),
}


# ---------------------------------------------------------------------------
# Config dataclass
# ---------------------------------------------------------------------------
@dataclass
class SimConfig:
    players: int = 1_000
    years: int = 10
    fields_per_draw: int = 8
    output_dir: Path = field(default_factory=lambda: Path("."))

    @property
    def total_draws(self) -> int:
        return self.years * DRAWS_PER_YEAR

    @property
    def cost_per_draw(self) -> float:
        return self.fields_per_draw * TICKET_COST

    @property
    def total_spent(self) -> float:
        return self.total_draws * self.cost_per_draw


# ---------------------------------------------------------------------------
# Probability calculation
# ---------------------------------------------------------------------------
def calculate_exact_probabilities() -> dict[str, float]:
    """Exact Eurojackpot winning probabilities via combinatorics."""
    total_main = math.comb(MAIN_POOL, MAIN_PICK)
    total_euro = math.comb(EURO_POOL, EURO_PICK)

    probabilities: dict[str, float] = {}
    prob_sum = 0.0

    for tier, (match_m, match_e) in PRIZE_TIERS.items():
        ways_main = math.comb(MAIN_PICK, match_m) * math.comb(
            MAIN_POOL - MAIN_PICK, MAIN_PICK - match_m
        )
        ways_euro = math.comb(EURO_PICK, match_e) * math.comb(
            EURO_POOL - EURO_PICK, EURO_PICK - match_e
        )
        prob = (ways_main / total_main) * (ways_euro / total_euro)
        probabilities[tier] = prob
        prob_sum += prob

    probabilities["No Win"] = 1.0 - prob_sum
    return probabilities


# ---------------------------------------------------------------------------
# Data loading
# ---------------------------------------------------------------------------
def load_empirical_payouts(
    file_path: str | Path,
    probabilities: dict[str, float],
) -> tuple[np.ndarray, np.ndarray]:
    """Load historical payout data and align with probability tiers."""
    log.info("[1/3] Loading historical payouts…")

    df = pl.read_json(str(file_path))
    payouts_df = df.select(pl.col("prizeDistribution")).unnest("prizeDistribution")
    avg_payouts: dict[str, float] = payouts_df.mean().to_dicts()[0]

    tiers = list(probabilities.keys())
    payouts = np.array(
        [0.0 if t == "No Win" else float(avg_payouts.get(t) or 0.0) for t in tiers],
        dtype=np.float64,
    )
    probs = np.fromiter(
        (probabilities[t] for t in tiers), dtype=np.float64, count=len(tiers)
    )
    probs /= probs.sum()  # normalise to guard against floating-point drift

    return payouts, probs


# ---------------------------------------------------------------------------
# Monte Carlo simulation
# ---------------------------------------------------------------------------
def run_monte_carlo(
    payouts: np.ndarray,
    probs: np.ndarray,
    cfg: SimConfig,
    *,
    seed: int | None = None,
) -> pl.DataFrame:
    """
    Vectorised Monte Carlo simulation.

    Memory layout: (total_draws, players, fields_per_draw)
    All draws are sampled in a single call for maximum throughput.
    """
    rng = np.random.default_rng(seed)

    log.info("[2/3] Simulating %d players over %d draws…", cfg.players, cfg.total_draws)
    log.info(
        "      Fields/draw: %d  |  Cost/draw: €%.2f  |  Total per player: €%.2f",
        cfg.fields_per_draw, cfg.cost_per_draw, cfg.total_spent,
    )

    # --- sample outcomes for every (draw, player, field) triple at once ---
    outcome_indices = rng.choice(
        len(payouts),
        size=(cfg.total_draws, cfg.players, cfg.fields_per_draw),
        p=probs,
    )
    raw_winnings: np.ndarray = payouts[outcome_indices]

    # sum over fields, subtract per-draw cost, cumsum over time
    net_per_draw = raw_winnings.sum(axis=2) - cfg.cost_per_draw  # (draws, players)
    bankroll = np.cumsum(net_per_draw, axis=0)                   # (draws, players)

    return pl.DataFrame(
        bankroll,
        schema=[f"Player_{i}" for i in range(cfg.players)],
    )


# ---------------------------------------------------------------------------
# Visualisation helpers
# ---------------------------------------------------------------------------
_PALETTE = {
    "jackpot":   "#00E5FF",
    "major":     "#FFD700",
    "minor":     "#FF9800",
    "profit":    "#00ff88",
    "loss":      "#ff4d4d",
    "worst":     "#ffffff",
    "breakeven": "#00bfff",
    "bg_fig":    "#0f1117",
    "bg_ax":     "#111827",
    "bg_legend": "#1f2937",
    "bg_text":   "#1f2937",
    "edge_text": "#ff4d4d",
}


def _profit_color(balance: float) -> str:
    if balance >= 1_000_000:
        return _PALETTE["jackpot"]
    if balance >= 100_000:
        return _PALETTE["major"]
    if balance >= 10_000:
        return _PALETTE["minor"]
    if balance > 0:
        return _PALETTE["profit"]
    return _PALETTE["loss"]


def _apply_dark_style(fig: plt.Figure, ax: plt.Axes) -> None:
    fig.patch.set_facecolor(_PALETTE["bg_fig"])
    ax.set_facecolor(_PALETTE["bg_ax"])


def _add_stats_box(
    ax: plt.Axes,
    avg: float,
    winners: int,
    players: int,
    fields: int,
) -> None:
    ax.text(
        0.02, 0.04,
        f"Stake per draw: €{fields * TICKET_COST:.2f}\n"
        f"Average result: €{avg:,.0f}\n"
        f"Profitable players: {winners}/{players}",
        transform=ax.transAxes,
        fontsize=12,
        bbox=dict(
            boxstyle="round,pad=0.6",
            facecolor=_PALETTE["bg_text"],
            edgecolor=_PALETTE["edge_text"],
            alpha=0.9,
        ),
    )


def _euro_formatter(x: float, _pos: int) -> str:
    return f"€ {int(x):,}"


# ---------------------------------------------------------------------------
# Plot generation
# ---------------------------------------------------------------------------
def analyze_and_plot(df: pl.DataFrame, cfg: SimConfig) -> None:
    """Generate log-scale (illusion) and linear-scale (reality) charts."""
    log.info("[3/3] Creating visualisations…")

    players = len(df.columns)
    final_balances = np.array(df.row(-1))

    avg_balance   = float(final_balances.mean())
    best_balance  = float(final_balances.max())
    worst_balance = float(final_balances.min())
    best_idx  = int(final_balances.argmax())
    worst_idx = int(final_balances.argmin())

    best_col  = df.columns[best_idx]
    worst_col = df.columns[worst_idx]

    winners_mask = final_balances > 0
    winners_count = int(winners_mask.sum())

    winner_cols = [df.columns[i] for i, w in enumerate(winners_mask) if w and df.columns[i] != best_col]
    loser_cols  = [df.columns[i] for i, w in enumerate(winners_mask) if not w]

    # cap losers shown to 500 for readability
    sample_losers = loser_cols[:500]

    plt.style.use("dark_background")

    chart_configs: list[tuple[str, str, bool, float | None]] = [
        (
            "eurojackpot_1_illusion_log.png",
            f"The Illusion: {cfg.years} Years Eurojackpot ({cfg.fields_per_draw} Fields)",
            True, None,
        ),
        (
            "eurojackpot_2_reality_linear.png",
            f"The Reality: {cfg.years} Years Eurojackpot ({cfg.fields_per_draw} Fields)",
            False, 5_000.0,
        ),
    ]

    for filename, title, use_log, ylim_max in chart_configs:
        fig, ax = plt.subplots(figsize=(12, 7))
        _apply_dark_style(fig, ax)

        # losers (faint red)
        for col in sample_losers:
            ax.plot(df[col].to_numpy(), color=_PALETTE["loss"],
                    alpha=0.15, linewidth=0.5, zorder=1)

        # non-best winners
        for col in winner_cols:
            series = df[col].to_numpy()
            ax.plot(series, color=_profit_color(series[-1]),
                    alpha=0.7, linewidth=0.8, zorder=3)

        # best & worst highlighted
        ax.plot(
            df[best_col].to_numpy(),
            color=_profit_color(best_balance),
            linewidth=1.5,
            label=f"Best: €{best_balance:,.0f}",
            zorder=6,
        )
        ax.plot(
            df[worst_col].to_numpy(),
            color=_PALETTE["worst"],
            linestyle="--",
            linewidth=1.5,
            label=f"Worst: €{worst_balance:,.0f}",
            zorder=6,
        )

        ax.axhline(0, color=_PALETTE["breakeven"], linestyle=":",
                   linewidth=1.5, label="Breakeven (€0)")

        if use_log:
            ax.set_yscale("symlog", linthresh=1_000)
            ax.set_ylabel("Net Bankroll (Log Scale)")
        else:
            ax.set_ylim(worst_balance - 500, ylim_max)
            ax.set_ylabel("Net Bankroll (Linear Scale)")

        ax.set_title(title, fontsize=18, fontweight="bold")
        ax.set_xlabel("Number of Draws (Tuesday & Friday)")
        ax.grid(True, axis="y", linestyle="--", alpha=0.15)
        ax.yaxis.set_major_formatter(ticker.FuncFormatter(_euro_formatter))

        legend = ax.legend(loc="upper right", fontsize=11)
        legend.get_frame().set_facecolor(_PALETTE["bg_legend"])
        legend.get_frame().set_edgecolor("#374151")

        _add_stats_box(ax, avg_balance, winners_count, players, cfg.fields_per_draw)

        plt.tight_layout()
        out_path = cfg.output_dir / filename
        fig.savefig(out_path, dpi=300, bbox_inches="tight")
        plt.close(fig)
        log.info("[+] Saved: %s", out_path)

    log.info("\nAll charts generated successfully.")


# ---------------------------------------------------------------------------
# Entry point
# ---------------------------------------------------------------------------
def main() -> None:
    data_file = Path("eurojackpot-history.json")

    if not data_file.exists():
        log.error("Error: '%s' not found.", data_file)
        return

    cfg = SimConfig(players=10_000, years=10, fields_per_draw=8)

    probabilities      = calculate_exact_probabilities()
    payouts, probs     = load_empirical_payouts(data_file, probabilities)
    df_simulation      = run_monte_carlo(payouts, probs, cfg, seed=42)
    analyze_and_plot(df_simulation, cfg)


if __name__ == "__main__":
    main()