#!/usr/bin/env python3
"""Exact-arithmetic checks for the evidence/discovery unit. Python standard library only.
Run: python foundations-evidence-discovery-check.py
This checks examples and finite algorithm equivalences, not general probability theorems.
Verification gates execute in both normal and optimized (-O) Python.
"""
from fractions import Fraction as F
from itertools import product
import json

def check(condition, message):
    """Verification gates remain active with python -O."""
    if not condition:
        raise AssertionError(message)

def ebh(values,q):
    m=len(values)
    order=sorted(range(m),key=lambda i:(-values[i],i))
    k=max([0]+[j for j,i in enumerate(order,1) if q*j*values[i]>=m])
    return frozenset(order[:k])

def bh(values,q):
    m=len(values)
    order=sorted(range(m),key=lambda i:(values[i],i))
    k=max([0]+[j for j,i in enumerate(order,1) if m*values[i]<=q*j])
    return frozenset(order[:k])

def by(values,q):
    h=sum((F(1,j) for j in range(1,len(values)+1)),F(0))
    return bh(values,q/h)

def closure(local,m):
    c=[False]*(1<<m)
    for mask in sorted(range(1,1<<m),key=int.bit_count,reverse=True):
        c[mask]=local[mask] and all(c[mask|1<<j] for j in range(m) if not(mask>>j&1))
    return c

def fdp_bounds(c,m):
    a=[mask.bit_count() if not c[mask] else 0 for mask in range(1<<m)]
    for j in range(m):
        for mask in range(1<<m):
            if mask>>j&1:a[mask]=max(a[mask],a[mask^(1<<j)])
    return a

def fdp(S,I0):return F(len(S&I0),max(len(S),1))

def stair(p,m,q):
    h=sum((F(1,j) for j in range(1,m+1)),F(0));c=q/(m*h)
    if p==0:return F(m,1)/q
    for j in range(1,m+1):
        if p<=j*c:return F(m,1)/(q*j)
    return F(0)

def main():
    result={}
    e=[F(30),F(18),F(17),F(2),F(0)]
    check(ebh(e,F(1,10))=={0,1,2}, 'e-BH exact rejection set')
    check(not ebh([F(30),F(18),F(16),F(2),F(0)],F(1,10)), 'e-BH threshold below boundary')
    result['e_bh_example']={'input':list(map(str,e)),'rejected_1_based':[1,2,3]}
    p=list(map(F,['.01','.02','.04','.20']))
    check(by(p,F(1,10))=={0,1}, 'BY exact rejection set');check(bh(p,F(1,10))=={0,1,2}, 'BH exact rejection set')
    result['by_example']={'H4':'25/12','c':'3/250','rejected_1_based':[1,2]}
    # Arbitrary-dependence counterexample; verify each discrete marginal exactly.
    atoms=[(F(1,10),(F(1,10),F(1))),(F(1,10),(F(1),F(1,10))),
           (F(1,10),(F(1,5),F(1,5))),(F(7,10),(F(1),F(1)))]
    for i in [0,1]:
        for t in sorted({v[i] for _,v in atoms}):check(sum(w for w,v in atoms if v[i]<=t)<=t, 'null marginal superuniformity')
    ordinary=sum(w*fdp(bh(v,F(1,5)),{0,1}) for w,v in atoms)
    corrected=sum(w*fdp(by(v,F(1,5)),{0,1}) for w,v in atoms)
    check(ordinary==F(3,10) and corrected==0, 'dependent-p FDR values')
    result['dependent_p_counterexample']={'BH_FDR':str(ordinary),'BY_FDR':str(corrected)}
    # Three-variable filtering counterexample, including one nonnull.
    atoms=[(F(1,15),(F(15),F(0),F(30))),(F(1,15),(F(0),F(15),F(30))),
           (F(13,15),(F(0),F(0),F(30)))]
    check(all(sum(w*v[i] for w,v in atoms)==1 for i in [0,1]), 'null e-value expectation')
    orig=sum(w*fdp(ebh(v,F(1,10)),{0,1}) for w,v in atoms)
    filt=sum(w*fdp(ebh(v,F(1,10))&{0,1},{0,1}) for w,v in atoms)
    check((orig,filt)==(F(1,15),F(2,15)), 'filtering FDR values')
    result['filtering_counterexample']={'original_FDR':str(orig),'filtered_FDR':str(filt)}
    # Full three-node Bonferroni closure.
    p=list(map(F,['.01','.04','.20']));local=[False]*8;local_p={}
    for mask in range(1,8):
        v=min(F(1),mask.bit_count()*min(p[j] for j in range(3) if mask>>j&1))
        local_p[str(mask)]=str(v);local[mask]=(v<=F(1,20))
    c=closure(local,3);check([i for i in range(3) if c[1<<i]]==[0], 'Bonferroni closure singleton')
    result['bonferroni_closure']={'local_p_by_bitmask':local_p,'closed_rejected_bitmasks':[i for i in range(1,8) if c[i]]}
    local=[False]+[6**i.bit_count()>=20 for i in range(1,8)];c=closure(local,3);a=fdp_bounds(c,3)
    check(all(a[i]==1 for i in range(1,8)) and a[0]==0, 'product-e simultaneous FDP bound')
    result['product_e_fdp_bounds']={'bounds_by_bitmask':a,'all_three_FDP_upper':'1/3'}
    # Test Boolean closure and subset-max transform on every possible local result table for m=3.
    for bits in range(128):
        local=[False]+[bool(bits>>(i-1)&1) for i in range(1,8)]
        c=closure(local,3)
        check(all(c[i]==all(local[j] for j in range(1,8) if j&i==i) for i in range(1,8)), 'full local-table closure equivalence')
        a=fdp_bounds(c,3)
        check(all(a[r]==max([0]+[i.bit_count() for i in range(1,8) if i&r==i and not c[i]]) for r in range(8)), 'full subset-max FDP equivalence')
    result['all_three_node_local_tables_checked']=128
    # Exact BY/e-BH staircase equivalence at every discontinuity and nearby points.
    count=0;m=3;q=F(11,100);h=F(11,6);c=q/(m*h)
    grid=sorted(set([F(0),F(1)]+[j*c for j in range(1,4)]+[j*c+F(1,10000) for j in range(1,4)]))
    for p in product(grid,repeat=m):
        check(by(p,q)==ebh([stair(x,m,q) for x in p],q), 'BY staircase e-BH equivalence')
        count+=1
    result['staircase_equivalence_cases']=count
    # Static likelihood-ratio examples and regression correction.
    check((F(3,4)**4)/(F(1,2)**4)==F(81,16), 'all-success split likelihood ratio')
    check((F(3,4)**2*F(1,4)**2)/(F(1,2)**4)==F(9,16), 'mixed split likelihood ratio')
    result['split_likelihood_ratios']=['81/16','9/16']
    ys=[(F(1,4),F(-1)),(F(1,2),F(0)),(F(1,4),F(1))]
    risk=lambda a:sum(w*(y-a)**2 for w,y in ys)
    check(risk(F(0))==F(1,2) and risk(F(1,2))-risk(F(0))==F(1,4), 'regression excess risk')
    result['regression_realized_excess']=['0','1/4']
    result['status']='PASS'
    result['scope']='Exact finite examples and algorithm equivalences only; not a substitute for proof review.'
    print(json.dumps(result,ensure_ascii=False,indent=2))

if __name__=='__main__':main()