"""
Which de-vig method gives the best-calibrated 'true' probability from Pinnacle's
1X2 closing odds (PSH/PSD/PSA)?  Read-only calibration study, no betting.

The live pipeline uses plain MULTIPLICATIVE de-vig (live_monitor.py:442-443,627-630):
    p_i = (1/o_i) / sum_j(1/o_j)
which is known to over-state longshot probabilities (it spreads the book's vig
proportionally, but real vig sits disproportionately on longshots). Our CLV edge
lives in the longshot tail, so a biased fair there would distort every CLV number.

Compare three de-vig methods against ACTUAL results:
  multiplicative : p_i = r_i / sum(r),                r_i = 1/o_i
  power          : p_i = r_i^k,  k chosen so sum(p)=1 (shrinks longshots more)
  shin           : insider-trader model, solve z so sum(p)=1 (favorite-longshot fix)

Scoring (lower=better): log-loss & Brier vs realized FTR (which outcome happened).
Plus a calibration table BY PREDICTED-PROB BIN, with the longshot bins (p<0.25)
broken out -- that is the slice our strategy actually bets. If multiplicative
over-states longshots, its predicted prob in those bins will exceed the empirical
win frequency, and Shin/power should be closer.
"""
import csv, glob, os, math
from collections import defaultdict

DIR = r'E:\soccer\league_data'

def _f(x):
    try:
        v = float(x)
        return v if v > 1.0 else None     # valid decimal odds are > 1.0
    except (TypeError, ValueError):
        return None

# ---- de-vig methods: take (oh,od,oa) -> dict {'H','D','A'} summing to 1 ----
def dv_mult(o):
    r = [1.0/x for x in o]; s = sum(r)
    return [ri/s for ri in r]

def dv_power(o):
    r = [1.0/x for x in o]
    # find k so sum(r_i^k) = 1.  r_i<1, sum(r)>1 -> need k>1. bisect.
    lo, hi = 0.5, 10.0
    for _ in range(60):
        k = (lo+hi)/2
        s = sum(ri**k for ri in r)
        if s > 1.0: lo = k
        else:       hi = k
    k = (lo+hi)/2
    p = [ri**k for ri in r]; s = sum(p)
    return [pi/s for pi in p]                  # tiny renorm for residual

def dv_shin(o):
    r = [1.0/x for x in o]; B = sum(r)
    def psum(z):
        return sum((math.sqrt(z*z + 4*(1-z)*ri*ri/B) - z)/(2*(1-z)) for ri in r)
    # at z=0, sum = sqrt(B) > 1; increasing z lowers the sum. bisect for sum=1.
    lo, hi = 0.0, 0.99
    for _ in range(60):
        z = (lo+hi)/2
        if psum(z) > 1.0: lo = z
        else:             hi = z
    z = (lo+hi)/2
    p = [(math.sqrt(z*z + 4*(1-z)*ri*ri/B) - z)/(2*(1-z)) for ri in r]
    s = sum(p)
    return [pi/s for pi in p]

METHODS = [('multiplicative', dv_mult), ('power', dv_power), ('shin', dv_shin)]
IDX = {'H': 0, 'D': 1, 'A': 2}

def load():
    rows = []
    for path in glob.glob(os.path.join(DIR, '*.csv')):
        with open(path, newline='', encoding='latin-1') as fh:
            for r in csv.DictReader(fh):
                ftr = (r.get('FTR') or '').strip().upper()
                if ftr not in ('H','D','A'): continue
                oh, od, oa = _f(r.get('PSH')), _f(r.get('PSD')), _f(r.get('PSA'))
                if None in (oh, od, oa): continue
                rows.append((oh, od, oa, ftr))
    return rows

def main():
    rows = load()
    print(f'Loaded {len(rows)} matches with Pinnacle closing 1X2 odds\n')
    booksum = sum(1/oh + 1/od + 1/oa for oh,od,oa,_ in rows) / len(rows)
    print(f'mean Pinnacle overround = {booksum:.4f}  ({(booksum-1)*100:.2f}% vig)\n')

    # overall log-loss / brier
    print('OVERALL fit vs realized outcome (lower = better):')
    stats = {}
    for name, fn in METHODS:
        ll = br = 0.0
        for oh,od,oa,ftr in rows:
            p = fn((oh,od,oa)); y = IDX[ftr]
            ll += -math.log(max(p[y], 1e-12))
            br += sum((p[k] - (1.0 if k==y else 0.0))**2 for k in range(3))
        stats[name] = (ll/len(rows), br/len(rows))
        print(f'   {name:14s} logloss={ll/len(rows):.5f}  brier={br/len(rows):.5f}')

    # calibration by predicted-prob bin (per-outcome predictions pooled H/D/A)
    edges = [0.0,0.05,0.10,0.15,0.20,0.25,0.35,0.50,0.70,1.01]
    print('\nCALIBRATION by predicted-prob bin (pooled over H/D/A outcomes):')
    print('   bin            n      mult_pred  mult_emp |  shin_pred  shin_emp |  pow_pred  pow_emp')
    # precompute per-method predictions
    preds = {name: [] for name,_ in METHODS}   # list of (p_outcome, hit 0/1)
    for oh,od,oa,ftr in rows:
        y = IDX[ftr]
        for name, fn in METHODS:
            p = fn((oh,od,oa))
            for k in range(3):
                preds[name].append((p[k], 1.0 if k==y else 0.0))
    for i in range(len(edges)-1):
        lo, hi = edges[i], edges[i+1]
        cells = {}
        for name,_ in METHODS:
            sel = [(pp,hh) for pp,hh in preds[name] if lo <= pp < hi]
            n = len(sel)
            pm = sum(pp for pp,_ in sel)/n if n else 0
            em = sum(hh for _,hh in sel)/n if n else 0
            cells[name] = (n, pm, em)
        n = cells['multiplicative'][0]
        tag = ' <- LONGSHOT' if hi <= 0.25 else ''
        print(f'   [{lo:.2f},{hi:.2f})  {n:6d}   {cells["multiplicative"][1]:7.3f}  {cells["multiplicative"][2]:7.3f}  | '
              f'{cells["shin"][1]:7.3f}  {cells["shin"][2]:7.3f}  | '
              f'{cells["power"][1]:7.3f}  {cells["power"][2]:7.3f}{tag}')

    # longshot-only summary: mean predicted vs empirical for p<0.25
    print('\nLONGSHOT TAIL (predicted p<0.25, the slice the strategy bets):')
    for name,_ in METHODS:
        sel = [(pp,hh) for pp,hh in preds[name] if pp < 0.25]
        n = len(sel)
        pm = sum(pp for pp,_ in sel)/n
        em = sum(hh for _,hh in sel)/n
        bias = pm - em
        print(f'   {name:14s} n={n:6d}  mean_pred={pm:.4f}  emp_freq={em:.4f}  '
              f'bias(pred-emp)={bias:+.4f}  ({"OVER" if bias>0 else "under"}-states)')

if __name__ == '__main__':
    main()