import json,csv,itertools,math
from pathlib import Path
H=Path(__file__).resolve().parent
raw=json.loads((H/"raw_worldbank.json").read_text())
# rebuild indicator->code map from screen2 definitions is gone; use data identity to drop exact duplicates
sig={k:tuple(sorted((i,round(v[1],12)) for i,v in raw[k].items())) for k in raw}
uniq={}
for k,s in sig.items():
    uniq.setdefault(s,k)
keep=sorted(uniq.values())
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 pval(r,n):
    if n<4 or abs(r)>=1: return 0.0
    t=r*math.sqrt((n-2)/(1-r*r)); return 2*(1-0.5*(1+math.erf(abs(t)/math.sqrt(2))))
res=[]
for a,b in itertools.combinations(keep,2):
    A,B=raw[a],raw[b]; ks=sorted(set(A)&set(B))
    if len(ks)<50: continue
    r=pear(rank([A[k][1] for k in ks]),rank([B[k][1] for k in ks]))
    if abs(r)>=0.95: continue   # near-duplicate measurements
    res.append({"var_a":a,"var_b":b,"n":len(ks),"spearman":round(r,3),"p":pval(r,len(ks))})
m=len(res); o=sorted(range(m),key=lambda i:res[i]["p"]); prev=1.0
for ri in range(m-1,-1,-1):
    i=o[ri]; prev=min(prev,res[i]["p"]*m/(ri+1)); res[i]["p_bh"]=prev
res.sort(key=lambda r_:abs(r_["spearman"]),reverse=True)
surv=sum(1 for r_ in res if r_["p_bh"]<0.05); bonf=sum(1 for r_ in res if r_["p"]<0.05/m)
with (H/"results_deduped.csv").open("w",newline="") as f:
    w=csv.DictWriter(f,fieldnames=list(res[0])); w.writeheader(); w.writerows(res)
(L:=H/"FINDINGS_DEDUPED.md").write_text("\n".join([
"# Association screen — deduplicated","",
f"After removing exact-duplicate indicators and near-duplicate pairs (|ρ| ≥ 0.95): **{m} tests**.","",
f"- Survive BH: **{surv} / {m}** ({100*surv/m:.1f}%)",
f"- Survive Bonferroni: **{bonf} / {m}**","",
"## Strongest non-duplicate associations","",
"| A | B | n | ρ | p (BH) |","|---|---|---:|---:|---:|"] +
[f"| {r_['var_a']} | {r_['var_b']} | {r_['n']} | {r_['spearman']:+.3f} | {r_['p_bh']:.1e} |" for r_ in res[:15]] +
["","## Why the survival rate is still high","",
"Even after deduplication, a large share of tests survive because the indicators are highly",
"correlated by construction: income, health, education and infrastructure are not separate",
"measurements of independent things. With ~200 countries and a dominant latent factor, the",
"effective number of independent tests is a small fraction of the nominal count.","",
"The honest interpretation: this screen maps the correlation structure of a statistical",
"catalogue, not of the world. Its outputs are hypotheses for designs that can separate",
"an indicator from national income.",""]))
print(f"unique indicators={len(keep)} tests={m} survive_BH={surv} bonferroni={bonf}")
print("top:", [(r_['var_a'],r_['var_b'],r_['spearman']) for r_ in res[:6]])