#!/usr/bin/env python3
"""Monte Carlo: how long does $100 actually last on slots?

Isolates VOLATILITY from RTP. All three slot models below are scaled in code to
exactly 96.00% RTP, so any difference in survival is variance, not payback.

Metrics per configuration:
  median spins until the bankroll cannot cover the next bet
  share still playing at 100 / 500 / 1000 spins
  share that ever reached 2x the starting bankroll at any point
  share ending a 500-spin session ahead, and the average balance then

Seed 20260804.
"""
import random
from statistics import mean, median

random.seed(20260804)

RTP = 0.96

def build(table):
    """Scale the paying buckets so RTP is exactly 0.96; 0x absorbs the rest."""
    paid = sum(m * p for m, p in table)
    scale = RTP / paid
    scaled = [(m, p * scale) for m, p in table]
    p_pay = sum(p for _, p in scaled)
    assert p_pay < 1, "paying probability exceeded 1"
    dist = [(0.0, 1 - p_pay)] + scaled
    assert abs(sum(m * p for m, p in dist) - RTP) < 1e-12
    return dist

# same RTP, deliberately different shapes
LOW  = build([(0.5, 0.30), (1.0, 0.20), (2.0, 0.09), (5.0, 0.012)])
MED  = build([(0.5, 0.18), (1.0, 0.09), (2.0, 0.055), (5.0, 0.022), (20.0, 0.0065), (100.0, 0.0012)])
HIGH = build([(0.5, 0.06), (2.0, 0.035), (10.0, 0.010), (50.0, 0.0030), (500.0, 0.00055), (5000.0, 0.00004)])

def variance(dist):
    m1 = sum(m * p for m, p in dist)
    m2 = sum(m * m * p for m, p in dist)
    return m2 - m1 * m1

def hit_rate(dist):
    return sum(p for m, p in dist if m > 0)

def run(dist, bankroll, bet, n, max_spins=5000):
    spins_lasted, doubled, alive = [], 0, {100: 0, 500: 0, 1000: 0}
    ahead500, bal500 = 0, []
    for _ in range(n):
        bal = float(bankroll)
        peak_hit = False
        s = 0
        while bal >= bet and s < max_spins:
            bal += spin(dist, random.random()) * bet - bet
            s += 1
            if not peak_hit and bal >= bankroll * 2:
                peak_hit = True
            for k in alive:
                if s == k:
                    alive[k] += 1
            if s == 500:
                bal500.append(bal)
                if bal > bankroll:
                    ahead500 += 1
        if s < 500:           # busted before 500 spins
            bal500.append(0.0)
        spins_lasted.append(s)
        doubled += peak_hit
    return {
        'median_spins': median(spins_lasted),
        'mean_spins': mean(spins_lasted),
        'alive': {k: v / n for k, v in alive.items()},
        'doubled': doubled / n,
        'ahead500': ahead500 / n,
        'avg_bal500': mean(bal500),
    }

def spin(dist, r):
    acc = 0.0
    for m, p in dist:
        acc += p
        if r < acc:
            return m
    return dist[-1][0]

N = 50_000
BANKROLL = 100

print(f"All models: RTP exactly {RTP*100:.2f}%. Bankroll ${BANKROLL}. {N:,} sessions each.\n")
for name, dist in (('LOW VOL', LOW), ('MEDIUM VOL', MED), ('HIGH VOL', HIGH)):
    print(f"{name:11} hit rate {hit_rate(dist)*100:5.2f}%   variance {variance(dist):9.2f}   "
          f"top multiplier x{max(m for m,_ in dist):g}")
print()

for bet in (0.20, 1.00, 2.00):
    print(f"================ ${bet:.2f} per spin ({BANKROLL/bet:.0f} spins of bankroll) ================")
    for name, dist in (('LOW VOL', LOW), ('MEDIUM VOL', MED), ('HIGH VOL', HIGH)):
        r = run(dist, BANKROLL, bet, N)
        print(f"--- {name}")
        print(f"    median spins survived   : {r['median_spins']:,.0f}   (mean {r['mean_spins']:,.0f})")
        print(f"    still playing @100 spins: {r['alive'][100]*100:5.1f}%")
        print(f"    still playing @500 spins: {r['alive'][500]*100:5.1f}%")
        print(f"    still playing @1000     : {r['alive'][1000]*100:5.1f}%")
        print(f"    ever doubled to ${BANKROLL*2}   : {r['doubled']*100:5.2f}%")
        print(f"    ahead after 500 spins   : {r['ahead500']*100:5.1f}%   avg balance ${r['avg_bal500']:.2f}")
        print()