"""Reviewer's independent checks (scipy Fisher, Monte Carlo, grid breakdown)."""
import json, numpy as np
from scipy.stats import fisher_exact, betabinom
from fractions import Fraction
rows=json.load(open('/b/results/exact-grid.json'))
F=lambda d:Fraction(int(d['numerator']),int(d['denominator']))
out={}
# 1. Exact selected case with scipy fisher_exact + scipy betabinom pmf (float)
m,k,a=5,20,5; n=m*k
q=betabinom.pmf(np.arange(k+1),k,a,a)
arm=np.array([1.0])
for _ in range(m): arm=np.convolve(arm,q)
cache={}
rej=0.0
for x in range(n+1):
  for y in range(n+1):
    p=fisher_exact([[x,n-x],[y,n-y]])[1]
    if p<=0.05+1e-12: rej+=arm[x]*arm[y]
out['scipy_selected_fisher']=rej
# 2. Monte Carlo
rng=np.random.default_rng(20261007); R=200000
P=rng.beta(a,a,size=(R,2,m)); X=rng.binomial(k,P).sum(axis=2)
pv={}
hits=0
for A,B in X:
  key=(int(A),int(B))
  if key not in pv: pv[key]=fisher_exact([[A,n-A],[B,n-B]])[1]
  hits+=pv[key]<=0.05+1e-12
phat=hits/R; out['mc_selected_fisher']=phat; out['mc_se']=(phat*(1-phat)/R)**0.5; out['mc_reps']=R
# 3. grid breakdown
exc=[r for r in rows if F(r['fisher_null_rejection'])>Fraction(1,20)]
from collections import Counter
out['exceed_by_regime']=dict(Counter(r['regime'] for r in exc))
out['exceed_beta_only']=sum(r['regime'] in ('50','5','1') for r in exc)
out['beta_rows']=sum(r['regime'] in ('50','5','1') for r in rows)
mx=max(rows,key=lambda r:F(r['fisher_null_rejection']))
out['max_row']={k2:mx[k2] for k2 in ('donors_per_arm','cells_per_donor','regime')}
mxb=max([r for r in rows if r['regime'] in ('50','5','1')],key=lambda r:F(r['fisher_null_rejection']))
out['max_beta_row']={k2:mxb[k2] for k2 in ('donors_per_arm','cells_per_donor','regime')}; out['max_beta_value']=float(F(mxb['fisher_null_rejection']))
k1=[r for r in rows if r['cells_per_donor']==1]
out['k1_distinct_fisher_values_per_m']={m2:len({str(F(r['fisher_null_rejection'])) for r in k1 if r['donors_per_arm']==m2}) for m2 in (3,5,10)}
out['oracle_min_max']=[float(min(F(r['oracle_null_rejection']) for r in rows)),float(max(F(r['oracle_null_rejection']) for r in rows))]
out['oracle_below_0.02_count']=sum(F(r['oracle_null_rejection'])<Fraction(1,50) for r in rows)
indep=[float(F(r['fisher_null_rejection'])) for r in rows if r['regime']=='independent']
out['independent_fisher_range']=[min(indep),max(indep)]
# monotonic in rho for fixed m,k>1?
nonmono=[]
for m2 in (3,5,10):
  for k2 in (2,5,10,20):
    seq=[F(next(r for r in rows if r['donors_per_arm']==m2 and r['cells_per_donor']==k2 and r['regime']==g)['fisher_null_rejection']) for g in ('independent','50','5','1','perfect_copy')]
    if any(seq[i]>seq[i+1] for i in range(4)): nonmono.append([m2,k2,[round(float(s),4) for s in seq]])
out['non_monotone_in_rho']=nonmono
print(json.dumps(out,indent=1))