"""Synthetic selection bias, not an RSI benchmark replication.
Run: python experiment.py --output ./results
Requires numpy and matplotlib. License: MIT.
"""
import argparse,csv,json,math,platform
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
parser=argparse.ArgumentParser();parser.add_argument('--output',default='results');args=parser.parse_args()
out=Path(args.output);out.mkdir(parents=True,exist_ok=True)
seed=20261003;repeats=1000;p=.6;n_dev=60;n_test=500;sizes=[1,8,32,128]
rng=np.random.default_rng(seed)
# Pre-generated, exchangeable candidates. Selection never uses test scores.
dev=rng.binomial(n_dev,p,size=(repeats,max(sizes)))
test=rng.binomial(n_test,p,size=(repeats,max(sizes)))
pmf=np.array([math.comb(n_dev,k)*p**k*(1-p)**(n_dev-k) for k in range(n_dev+1)])
tails=np.array([sum(pmf[k:]) for k in range(n_dev+1)])
def stats(a):return {'mean':float(a.mean()),'mc_se':float(a.std(ddof=1)/math.sqrt(repeats))}
def wilson(k,n,z=1.959963984540054):
 q=k/n;center=(q+z*z/(2*n))/(1+z*z/n);half=z*math.sqrt(q*(1-q)/n+z*z/(4*n*n))/(1+z*z/n)
 return [max(0,center-half),min(1,center+half)]
rows=[];raw=[]
for n in sizes:
 idx=np.argmax(dev[:,:n],axis=1);d=dev[np.arange(repeats),idx]/n_dev;t=test[np.arange(repeats),idx]/n_test
 gates={'fixed_075':dev[:,:n]>=45,'uncorrected_005':tails[dev[:,:n]]<=.05,'bonferroni':tails[dev[:,:n]]<=.05/n}
 row={'candidates':n,'dev':stats(d),'test':stats(t),'gap':stats(d-t),'gates':{}}
 for label,passing in gates.items():
  event=passing.any(axis=1);count=int(event.sum());q=(sum(pmf[45:]) if label=='fixed_075' else float(pmf[tails<=(.05/n if label=='bonferroni' else .05)].sum()))
  theory=1-(1-q)**n
  row['gates'][label]={'false_adoption_rate':float(event.mean()),'mc_se':float(math.sqrt(event.mean()*(1-event.mean())/repeats)),'wilson95':wilson(count,repeats),'theory':float(theory)}
  # Independent candidate evaluations imply the exact family-level calculation.
  assert abs(event.mean()-theory)<max(.015,5*math.sqrt(theory*(1-theory)/repeats))
 assert abs(t.mean()-p)<5*math.sqrt(p*(1-p)/(n_test*repeats))
 assert row['gates']['bonferroni']['theory']<=.05+1e-10
 for i in range(repeats):raw.append([n,i,int(idx[i]),float(d[i]),float(t[i])])
 rows.append(row)
assert all(rows[i]['dev']['mean']<=rows[i+1]['dev']['mean'] for i in range(3))
result={'protocol':{'seed':seed,'repeats':repeats,'p':p,'dev_tasks':n_dev,'test_tasks':n_test,'candidate_sizes':sizes,'selection':'dev argmax; first tie; candidates fixed before evaluation','scope':'independent synthetic noise, no real ability gain, no RSI replication'},'environment':{'python':platform.python_version(),'numpy':np.__version__,'matplotlib':matplotlib.__version__},'results':rows}
(out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
with (out/'trials.csv').open('w') as f:
 w=csv.writer(f);w.writerow(['candidates','replicate','selected_index','dev_score','test_score']);w.writerows(raw)
plt.rcParams.update({'font.size':11,'axes.spines.top':False,'axes.spines.right':False,'svg.fonttype':'none'})
fig,axes=plt.subplots(1,2,figsize=(11,4.3),layout='constrained');fig.set_facecolor('#faf9f6')
x=np.arange(len(sizes))
for ax in axes:ax.set_facecolor('#faf9f6');ax.set_xticks(x,labels=sizes);ax.set_xlabel('Number of candidates');ax.grid(axis='y',alpha=.18)
for key,color,label in [('dev','#81543d','Selected development score'),('test','#42676a','Independent test score')]:
 axes[0].errorbar(x,[r[key]['mean']*100 for r in rows],yerr=[r[key]['mc_se']*1.96*100 for r in rows],marker='o',color=color,label=label,capsize=3)
axes[0].axhline(60,color='#777',ls='--',lw=1);axes[0].set_ylim(55,85);axes[0].set_ylabel('Mean score (%)');axes[0].set_title('Same ability; a better-looking winner',loc='left',fontsize=12);axes[0].legend(frameon=False,fontsize=9)
for key,color,label in [('fixed_075','#81543d','Fixed score >= 75%'),('uncorrected_005','#a28755','Uncorrected p <= .05'),('bonferroni','#42676a','Bonferroni .05 / N')]:
 axes[1].plot(x,[r['gates'][key]['false_adoption_rate']*100 for r in rows],marker='o',color=color,label=label)
axes[1].set_ylim(-2,103);axes[1].set_ylabel('False adoption per search (%)');axes[1].set_title('More candidates; more false discoveries',loc='left',fontsize=12);axes[1].legend(frameon=False,fontsize=9)
fig.savefig(out/'selection-bias.svg');fig.savefig(out/'selection-bias.png',dpi=180);plt.close(fig)
print(json.dumps(result,indent=2))