import csv, json, io, random, urllib.request
from pathlib import Path
H=Path(__file__).resolve().parent; ROOT=H.parents[1]
random.seed(4); N=20000
rows=[]
for off in range(0,6000,1000):
    u=f"https://api.cpicpgx.org/v1/change_log?limit=1000&offset={off}"
    r=urllib.request.Request(u,headers={"Accept":"text/csv","User-Agent":"incoprea-research/1.0"})
    body=urllib.request.urlopen(r,timeout=90).read().decode()
    chunk=list(csv.DictReader(io.StringIO(body)))
    rows+=chunk
    if len(chunk)<1000: break
print("change_log rows:",len(rows))
freq={}; first={}
for r in rows:
    if r.get("type")!="FREQUENCY": continue
    g=r.get("entityid") or ""
    freq[g]=freq.get(g,0)+1
    d=r.get("date") or ""
    if d and (g not in first or d<first[g]): first[g]=d
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['cov']=float(g['ancestry_groups_with_data'])
    g['ev']=float(freq.get(g['gene'],0)); g['first']=first.get(g['gene'],"")
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 sp(x,y): return pear(rank(x),rank(y))
def pp(x,y,obs,n=N):
    ys=list(y); h=0
    for _ in range(n):
        random.shuffle(ys)
        if abs(sp(x,ys))>=abs(obs)-1e-12: h+=1
    return (h+1)/(n+1)
def partial(x,y,z):
    rx,ry,rz=rank(x),rank(y),rank(z); n=len(rx)
    mx,my,mz=sum(rx)/n,sum(ry)/n,sum(rz)/n
    szz=sum((a-mz)**2 for a in rz)
    bx=sum((a-mz)*(b-mx) for a,b in zip(rz,rx))/szz
    by=sum((a-mz)*(b-my) for a,b in zip(rz,ry))/szz
    ex=[a-mx-bx*(c-mz) for a,c in zip(rx,rz)]; ey=[a-my-by*(c-mz) for a,c in zip(ry,rz)]
    return pear(ex,ey)
x=[g['pub'] for g in genes]; y=[g['cov'] for g in genes]; z=[g['ev'] for g in genes]
r1=sp(x,y); r2=sp(z,y); r3=sp(x,z); rp=partial(x,y,z)
res={"n":len(genes),"rows_scanned":len(rows),
 "publications_to_coverage":{"rho":round(r1,3),"p":round(pp(x,y,r1),5)},
 "curation_to_coverage":{"rho":round(r2,3),"p":round(pp(z,y,r2),5)},
 "publications_to_curation":{"rho":round(r3,3),"p":round(pp(x,z,r3),5)},
 "publications_to_coverage_given_curation":{"rho":round(rp,3)},
 "with_events":sum(1 for g in genes if g['ev']>0),"genes":len(genes)}
(H/'mediation.json').write_text(json.dumps(res,indent=2))
L=["# Mediation test: is coverage driven by curation or by literature?","",
f"Scanned {len(rows)} CPIC change-log rows; {res['with_events']} of {len(genes)} genes have at least one FREQUENCY event.","",
"| Path | ρ | p |","|---|---:|---:|",
f"| publications → coverage | {r1:+.3f} | {res['publications_to_coverage']['p']} |",
f"| curation events → coverage | {r2:+.3f} | {res['curation_to_coverage']['p']} |",
f"| publications → curation events | {r3:+.3f} | {res['publications_to_curation']['p']} |",
f"| **publications → coverage, controlling for curation** | **{rp:+.3f}** | partial |","",
"## Reading","",
f"If curation fully mediated the association, the partial correlation would fall toward zero. It is {rp:+.3f}. "
"That means curation activity and publications are not interchangeable explanations in these data; the "
"coverage–attention link is not simply 'the gene that was edited most often'.","",
"## Limits","",
"- n = 21; partial correlations at this size are unstable.",
"- Change-log events count edits of any kind within a frequency table, not additions of ancestry groups.",
"- Event counts are not dated coverage levels; this is mediation-style evidence, not a longitudinal test.",
"- The true longitudinal test still requires two dated CPIC snapshots.",""]
(H/'MEDIATION.md').write_text("\n".join(L))
print(json.dumps(res,indent=1))