"""Reviewer sensitivity, not replacement trial evidence.
CR2 under a working independent diagonal covariance (v+tau2), conditional on fitted tau2.
Satterthwaite df from the Gaussian working-model quadratic form. Also cluster deletion.
Original author helpers are reused for extraction and fitted REML weights only.
"""
import json
from collections import defaultdict
import numpy as np
from scipy import stats
import meta

def cr2(y,X,labels,v,tau,L):
    y=np.asarray(y,float);X=np.asarray(X,float);v=np.asarray(v,float);L=np.asarray(L,float)
    w=1/(v+tau);Xs=X*np.sqrt(w[:,None]);ys=y*np.sqrt(w)
    B=np.linalg.inv(Xs.T@Xs);b=B@Xs.T@ys;r=ys-Xs@b;H=Xs@B@Xs.T;R=np.eye(len(y))-H
    meat=np.zeros_like(B);gs=[];minimum_eigenvalue=1.0
    for s in sorted(set(labels)):
        ii=np.flatnonzero(np.array(labels)==s);vals,vecs=np.linalg.eigh(np.eye(len(ii))-H[np.ix_(ii,ii)])
        minimum_eigenvalue=min(minimum_eigenvalue,float(vals.min()));assert vals.min()>1e-10
        A=(vecs*(1/np.sqrt(vals)))@vecs.T
        u=Xs[ii].T@A@r[ii];meat+=np.outer(u,u)
        gs.append(R[:,ii]@A@Xs[ii]@B@L)
    V=B@meat@B;G=np.array(gs);M=G@G.T;df=np.trace(M)**2/np.sum(M*M)
    est=float(L@b);se=float(np.sqrt(L@V@L));t=est/se
    return {'estimate':est,'se':se,'working_independence_satterthwaite_df':float(df),'ci95':[float(est-stats.t.ppf(.975,df)*se),float(est+stats.t.ppf(.975,df)*se)],'one_sided_positive_p':float(stats.t.sf(t,df)),'minimum_adjustment_eigenvalue':minimum_eigenvalue}

# Analytic fixture: independent observations, intercept only, CR2 equals usual mean SE and df n-1.
z=np.array([1.,2.,4.,5.,7.]);f=cr2(z,np.ones((5,1)),list('abcde'),np.ones(5),0,[1]);assert abs(f['se']-np.std(z,ddof=1)/np.sqrt(5))<1e-12;assert abs(f['working_independence_satterthwaite_df']-4)<1e-12
per=meta.primary_effects('data/extraction.csv');rows=meta.by_zone(per)+meta.secondary_effects('data/secondary.csv')
y=np.array([r['g'] for r in rows]);v=np.array([r['v'] for r in rows]);labels=[r['study'] for r in rows]
z=np.array([r['zone']=='beyond' for r in rows],float);X=np.column_stack([1-z,z]);tau=meta.reml_tau2(y,v,X)
primary=cr2(y,X,labels,v,tau,[-1,1]);assert abs(primary['estimate']-0.4836)<.00005
agg=defaultdict(list)
for e in per:agg[e['study'],e['interval_s']].append(e)
pi=[{'study':s,'interval_s':t,'g':float(np.mean([e['g'] for e in es])),'v':float(np.mean([e['v'] for e in es]))} for (s,t),es in agg.items()]
x=np.log10([r['interval_s'] for r in pi]);X=np.column_stack([np.ones(len(pi)),x]);y=np.array([r['g'] for r in pi]);v=np.array([r['v'] for r in pi]);labels=[r['study'] for r in pi];tau=meta.reml_tau2(y,v,X)
slope=cr2(y,X,labels,v,tau,[0,1]);assert abs(slope['estimate']-.6302)<.00005
loo={s:meta.slope_log_interval([r for r in pi if r['study']!=s]) for s in sorted(set(labels))}
within={s:len({r['interval_s'] for r in pi if r['study']==s}) for s in sorted(set(labels))}
results={'primary_working_independence_CR2':primary,'slope_working_independence_CR2':slope,'slope_cluster_deletion_author_CR1':loo,'distinct_targets_per_primary_study':within,'fixture_passed':True,'note':'Post hoc reviewer sensitivity, conditional on existing extraction and diagonal working covariance; not a unique correct covariance model and not new primary data.'}
with open('results/reviewer-diagnostics.json','w') as f:json.dump(results,f,indent=2)
print(json.dumps(results,indent=2))