"""Verifier-side aggregate checks using only newly generated model fits."""
import json
import math
from pathlib import Path
import numpy as np
from scipy.stats import beta, binom, norm, t, nct

JOB = Path('work/jobs/job-9687e2777a2fbfd255c4af888e10184f')
EVIDENCE = JOB/'evidence'
generated = EVIDENCE/'results'
if not (generated/'associations.json').exists():
    choices=list(JOB.glob('.run-*/results/associations.json'))
    if len(choices)!=1: raise RuntimeError('Fresh run outputs are not available')
    generated=choices[0].parent
records=json.loads((generated/'associations.json').read_text())
summary=json.loads((generated/'R1.json').read_text())
rows=list(records.values()) if isinstance(records,dict) else records

def bh(values):
    values=np.asarray(values)
    order=np.argsort(values,kind='stable')
    sorted_q=np.minimum.accumulate((values[order]*len(values)/np.arange(1,len(values)+1))[::-1])[::-1]
    output=np.empty(len(values));output[order]=np.minimum(1,sorted_q)
    return output

def share(flags):
    count=sum(bool(f) for f in flags);size=len(flags)
    low=0 if count==0 else beta.ppf(.025,count,size-count+1)
    high=1 if count==size else beta.ppf(.975,count+1,size-count)
    return {'k':count,'n':size,'share':round(count/size,3),'ci':[round(float(low),3),round(float(high),3)]}

def median(values):
    x=np.sort([v for v in values if v is not None and math.isfinite(v)])
    n=len(x);low=int(binom.ppf(.025,n,.5));high=n-low+1
    interval=[round(float(x[low-1]),3),round(float(x[high-1]),3)] if low>0 else [None,None]
    return {'n':n,'median':round(float(np.median(x)),3),'ci':interval,
            'quartiles':[round(float(v),3) for v in np.quantile(x,[.25,.75])]}

def close(actual,expected):
    if isinstance(expected,dict):
        for key,value in expected.items(): close(actual[key],value)
    elif isinstance(expected,list):
        assert len(actual)==len(expected)
        for a,b in zip(actual,expected):close(a,b)
    elif isinstance(expected,(float,int)):
        assert abs(actual-expected)<1e-9,(actual,expected)
    else:assert actual==expected,(actual,expected)

q=bh([r['replication']['p'] for r in rows])
q_heterogeneity=bh([r['heterogeneity_p'] for r in rows])
q_difference=bh([r['published_difference_p'] for r in rows])
replicated=[];informative=[]
for i,r in enumerate(rows):
    assert abs(float(r['q'])-q[i])<1e-10
    assert abs(float(r['heterogeneity_q'])-q_heterogeneity[i])<1e-10
    assert abs(float(r['published_difference_q'])-q_difference[i])<1e-10
    fit=r['replication'];published=r['published_b'][0]
    df=math.inf if fit['df'] is None else fit['df']
    if not math.isfinite(df):assert not fit['weighted']
    critical=t.ppf(.975,df);ncp=published/fit['se']
    power=float(norm.cdf(-critical-ncp)+norm.sf(critical-ncp)) if math.isinf(df) else float(nct.cdf(-critical,df,ncp)+nct.sf(critical,df,ncp))
    assert abs(power-r['power'])<1e-7,(r['id'],power,r['power'])
    informative.append(power>=.8)
    replicated.append(np.sign(fit['b'])==np.sign(published) and q[i]<.05)
    assert bool(r['informative'])==informative[-1]
    assert bool(r['replicated'])==replicated[-1]
    assert abs(fit['p']-2*t.sf(abs(fit['b']/fit['se']),df))<1e-9
    original=r['original']
    assert bool(r['reproduced']) == (r['published']['low']<=original['estimate']<=r['published']['high'])
    assert abs(r['ratio']-fit['b']/published)<1e-10
    assert abs(r['ratio_own']-fit['b']/r['harmonized']['b'])<1e-10

checks={
 'replicated':share(replicated),
 'informative':sum(informative),
 'replicated_informative':share([a for a,b in zip(replicated,informative) if b]),
 'ratio':median([r['ratio'] for r in rows]),
 'ratio_informative':median([r['ratio'] for r,b in zip(rows,informative) if b]),
 'ratio_own':median([r['ratio_own'] for r in rows]),
 'heterogeneous':share(q_heterogeneity<.05),
 'differs_from_published':share(q_difference<.05),
 'in_published_ci':share([r['published']['low']<=r['replication']['estimate']<=r['published']['high'] for r in rows]),
 'same_sign_p05':share([r['same_sign'] and r['replication']['p']<.05 for r in rows]),
 'reversed':share([not r['same_sign'] and q[i]<.05 for i,r in enumerate(rows)]),
 'replicated_reproduced':share([replicated[i] for i,r in enumerate(rows) if r['reproduced']]),
}
close(summary['replication'],checks)
close(summary['reproduction']['reproduced'],share([r['reproduced'] for r in rows]))
close(summary['reproduction']['ratio_original'],median([r['ratio_original'] for r in rows]))
evidence={'source':'fresh offline reproduction outputs, not declared results',
 'checks_passed':['Independent BH correction for replication, heterogeneity and published-difference tests',
 'Replication test p-values and noncentral-t power', 'Replication and informative classifications',
 'Exact binomial intervals for key shares','Order-statistic median intervals and quartiles',
 'Log-scale effect ratios and primary aggregate summaries'],
 'associations_checked':len(rows),'scope':'Does not independently authenticate every original paper extraction or judgment of a departure.'}
EVIDENCE.mkdir(exist_ok=True)
(EVIDENCE/'independent-aggregation.json').write_text(json.dumps(evidence,indent=2)+'\n')
print(json.dumps(evidence,indent=2))