#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
fairness_audit.py — is the draw fair in SPACE (not just over time)?
===================================================================
Time-domain tests (recency, overlaps, diagonals) can all pass while a physical
machine still pairs certain numbers together. This audits that: for Joker (5/45)
and Eurojackpot (5/50) it checks whether pairs and triples of numbers co-occur
more (or less) often than a fair draw allows.

  * every pair's co-occurrence as a z-score vs its Binomial expectation,
  * the max |z| judged against a fair-null (random draw-histories), to beat the
    multiple-testing trap of ~1000 pairs,
  * triples vs Poisson,
  * the "consecutive numbers never come up" folklore.

Needs each game's hist_df.csv (st1..st5). Writes 2 PNGs to
static/img/blog/fairness-audit/ and prints the full summary.
"""
from __future__ import annotations
from math import comb
from itertools import combinations
from collections import Counter
from pathlib import Path

import numpy as np
import pandas as pd
from scipy import stats
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

ST = ["st1", "st2", "st3", "st4", "st5"]
GAMES = {"Joker (5/45)": ("games/joker/data/hist_df.csv", 45),
         "Eurojackpot (5/50)": ("games/eurojackpot/data/hist_df.csv", 50)}
COL = {"Joker (5/45)": "#e0a417", "Eurojackpot (5/50)": "#3ba7d6"}
IMG = Path("static/img/blog/fairness-audit"); IMG.mkdir(parents=True, exist_ok=True)
B_NULL = 300
plt.rcParams.update({"figure.dpi": 130, "font.size": 11, "axes.grid": True, "grid.alpha": .25})


def pair_z(draws, npool, n):
    p = comb(npool - 2, 3) / comb(npool, 5); exp = n * p
    C = Counter()
    for d in draws:
        for pr in combinations(sorted(int(x) for x in d), 2):
            C[pr] += 1
    allp = list(combinations(range(1, npool + 1), 2))
    obs = np.array([C.get(pr, 0) for pr in allp])
    z = (obs - exp) / np.sqrt(exp * (1 - p))
    return z, obs, exp, p, C


zfig, zax = plt.subplots(1, 2, figsize=(11, 4.3))
nfig, nax = plt.subplots(1, 2, figsize=(11, 4.3))
rng = np.random.default_rng(0)

for gi, (game, (path, npool)) in enumerate(GAMES.items()):
    h = pd.read_csv(path); h.columns = [str(c).strip() for c in h.columns]
    draws = h[ST].astype(int).to_numpy(); n = len(draws); c = COL[game]

    freq = np.bincount(draws.ravel(), minlength=npool + 1)[1:]
    _, ps = stats.chisquare(freq)
    z, obs, exp, p, C = pair_z(draws, npool, n)
    realmax = float(np.abs(z).max())

    # fair-null distribution of max|z|
    nullmax = []
    for _ in range(B_NULL):
        rd = np.array([rng.choice(npool, 5, replace=False) + 1 for _ in range(n)])
        zz, *_ = pair_z(rd, npool, n)
        nullmax.append(float(np.abs(zz).max()))
    nullmax = np.array(nullmax)

    # triples + consecutive
    p_tri = comb(npool - 3, 2) / comb(npool, 5); et = n * p_tri
    Ct = Counter()
    for d in draws:
        for tr in combinations(sorted(int(x) for x in d), 3):
            Ct[tr] += 1
    obs_cons = sum(v for (a, b), v in C.items() if b == a + 1)
    exp_cons = (npool - 1) * exp

    print(f"\n=== {game} | {n} draws ===")
    print(f"  singles uniformity chi2 p={ps:.3f} (freq {freq.min()}-{freq.max()}, exp {5*n/npool:.0f})")
    print(f"  pairs: exp/pair={exp:.2f}  z mean={z.mean():+.3f} sd={z.std():.3f}  |z|>2={np.mean(np.abs(z)>2):.4f} (null .0455)  |z|>3={np.mean(np.abs(z)>3):.4f} (null .0027)")
    print(f"  max|z| real={realmax:.2f}  fair-null mean={nullmax.mean():.2f} 95%<= {np.percentile(nullmax,95):.2f}  p={(nullmax>=realmax).mean():.3f}")
    print(f"  triples: exp/triple={et:.2f}  busiest={max(Ct.values())} (Poisson-consistent)")
    print(f"  consecutive (i,i+1): obs={obs_cons} exp={exp_cons:.0f} ratio={obs_cons/exp_cons:.3f}")

    ax = zax[gi]
    zs = np.sort(z)
    q = stats.norm.ppf((np.arange(1, len(zs) + 1) - 0.5) / len(zs))
    ax.plot([-4, 4], [-4, 4], "k--", lw=1.4, label="fair machine  N(0,1)")
    ax.scatter(q, zs, s=7, color=c, alpha=.55, label="observed pair z-scores")
    ax.set_xlim(-4, 4); ax.set_ylim(-4, 4)
    ax.set_title(f"{game}\npair z-scores vs normal (Q–Q)", fontsize=10)
    ax.set_xlabel("theoretical normal quantile"); ax.set_ylabel("observed pair z-score"); ax.legend(fontsize=9)

    ax = nax[gi]
    ax.hist(nullmax, bins=28, color="#999", alpha=.75, label="fair-null max|z|")
    ax.axvline(realmax, color=c, lw=2.5, label=f"real = {realmax:.2f}")
    ax.set_title(f"{game}\nmost-extreme pair vs null  (p={(nullmax>=realmax).mean():.2f})", fontsize=10)
    ax.set_xlabel("max |z| over all pairs"); ax.set_ylabel("null trials"); ax.legend(fontsize=9)

for fig, name in [(zfig, "img1"), (nfig, "img2")]:
    fig.tight_layout(); fig.savefig(IMG / f"{name}.png", bbox_inches="tight"); plt.close(fig)
print(f"\ncharts -> {IMG}")