import csv, json, random, math
from pathlib import Path
H=Path(__file__).resolve().parent; ROOT=H.parents[1]
random.seed(9); N=20000
genes=list(csv.DictReader(open(ROOT/'data/processed/gene_reference.csv')))
pub=json.loads((H/'raw_pubmed.json').read_text())['queries']
for g in genes:
    g['pub']=float(pub[g['gene']]['pgx']); g['all']=float(pub[g['gene']]['all'])
    g['cov']=float(g['ancestry_groups_with_data']); g['pairs']=float(g['n_actionable_pairs'])
def rank(v):
    o=sorted(range(len(v)),key=lambda i:v[i]); r=[0.0]*len(v); i=0
    while i<len(o):
        j=i
        while j+1<len(o) and v[o[j+1]]==v[o[i]]: j+=1
        a=(i+j)/2+1
        for x in range(i,j+1): r[o[x]]=a
        i=j+1
    return r
def pear(x,y):
    n=len(x); mx=sum(x)/n; my=sum(y)/n
    num=sum((a-mx)*(b-my) for a,b in zip(x,y))
    dx=sum((a-mx)**2 for a in x)**.5; dy=sum((b-my)**2 for b in y)**.5
    return num/(dx*dy) if dx and dy else 0.0
def spear(x,y): return pear(rank(x),rank(y))
def pperm(x,y,obs,n=N):
    ys=list(y); h=0
    for _ in range(n):
        random.shuffle(ys)
        if abs(spear(x,ys))>=abs(obs)-1e-12: h+=1
    return (h+1)/(n+1)
def test(sub,label):
    x=[g['pub'] for g in sub]; y=[g['cov'] for g in sub]
    r=spear(x,y); p=pperm(x,y,r)
    return {"subset":label,"n":len(sub),"spearman":round(r,3),"p_perm":round(p,5)}
T=[]
T.append(test(genes,"all genes"))
T.append(test([g for g in genes if not g['gene'].startswith('CYP')],"excluding CYP family"))
T.append(test([g for g in genes if not g['gene'].startswith('HLA')],"excluding HLA"))
T.append(test([g for g in genes if g['gene'].startswith('CYP')],"CYP family only"))
# residuals from a log-linear fit on ranks
rx,ry=rank([g['pub'] for g in genes]),rank([g['cov'] for g in genes])
n=len(rx); mx,my=sum(rx)/n,sum(ry)/n
b=sum((a-mx)*(c-my) for a,c in zip(rx,ry))/sum((a-mx)**2 for a in rx)
for g,a,c in zip(genes,rx,ry): g['resid']=round(c-(my+b*(a-mx)),2)
genes.sort(key=lambda g:g['resid'])
neg=[(g['gene'],g['pub'],g['cov'],g['resid']) for g in genes[:5]]
pos=[(g['gene'],g['pub'],g['cov'],g['resid']) for g in genes[-5:]]
L=["# Robustness of the attention–coverage finding","",
"Does the association survive removing the family that supplies most of the literature?","",
"| Subset | n | ρ | p |","|---|---:|---:|---:|"]
for t in T: L.append(f"| {t['subset']} | {t['n']} | {t['spearman']:+.3f} | {t['p_perm']:.4f} |")
L+=["","## Genes that deviate most from the attention model","",
"Residual = coverage rank − predicted coverage rank from publication rank.","",
"**Coverage lower than attention predicts** (candidate gaps):","",
"| Gene | PGx publications | Groups | Residual |","|---|---:|---:|---:|"]
for a,b_,c,d in neg: L.append(f"| {a} | {int(b_)} | {int(c)}/9 | {d:+.1f} |")
L+=["","**Coverage higher than attention predicts:**","",
"| Gene | PGx publications | Groups | Residual |","|---|---:|---:|---:|"]
for a,b_,c,d in pos: L.append(f"| {a} | {int(b_)} | {int(c)}/9 | {d:+.1f} |")
L+=["","## Verdict","",
"The finding survives removing the CYP family only if the second row stays positive and material. "
"If it collapses, the association is a family effect, not a general one — and the headline claim "
"should be withdrawn rather than qualified.",""]
(H/'ROBUSTNESS.md').write_text("\n".join(L))
json.dump({"tests":T,"low_residual":neg,"high_residual":pos},open(H/'robustness.json','w'),indent=2)
for t in T: print(f"{t['subset']:26s} n={t['n']:2d} rho={t['spearman']:+.3f} p={t['p_perm']:.4f}")
print("gaps:",neg[:3])