scripts/post_analysis.py

import json,collections,numpy as np,pandas as pd
from pathlib import Path
rows=json.load(open('reports/test_subject_results.json'));primary=json.load(open('reports/demo_subject_models_test.json'));rng=np.random.default_rng(2026)
p=np.array([r['metrics']['accuracy'] for r in sorted(primary,key=lambda r:r['subject'])]);out={}
for model in ['precue','late_matched','permuted','loso_bandpower','loso_svm']:
 r=sorted([r for r in rows if r['model']==model],key=lambda r:r['subject'])
 if len(r)!=15:continue
 a=np.array([q['metrics']['accuracy'] for q in r]);d=p-a;boot=rng.choice(d,(10000,15),replace=True).mean(1);out['primary_minus_'+model]={'mean_accuracy_difference':float(d.mean()),'paired_subject_bootstrap_ci95':np.quantile(boot,[.025,.975]).tolist(),'control_accuracy':float(a.mean())}
# Matched duration comparison explicitly separates late imagination from pre-onset fixation.
pre=sorted([r for r in rows if r['model']=='precue'],key=lambda r:r['subject']);late=sorted([r for r in rows if r['model']=='late_matched'],key=lambda r:r['subject'])
if len(pre)==len(late)==15:
 d=np.array([b['metrics']['accuracy']-a['metrics']['accuracy'] for a,b in zip(pre,late)]);out['matched_late_minus_precue']={'mean_accuracy_difference':float(d.mean()),'paired_subject_bootstrap_ci95':np.quantile(rng.choice(d,(10000,15),replace=True).mean(1),[.025,.975]).tolist()}
json.dump(out,open('reports/control_comparisons.json','w'),indent=2)
# Per-class recall and calibration for primary, preserving the classifier's actual predicted label.
ys=[];ps=[];cs=[]
for r in primary:
 z=np.load('reports/predictions/'+r['name']+'.npz');ys.extend(z['y']);ps.extend(z['pred']);cs.extend(z['prob'][np.arange(len(z['pred'])),z['pred']])
y=np.array(ys);p=np.array(ps);c=np.array(cs);bins=[];ece=0
for low in np.arange(0,1,.1):
 mask=(c>=low)&(c<low+.1 if low<.9 else c<=1)
 if mask.any():
  accuracy=float((y[mask]==p[mask]).mean());confidence=float(c[mask].mean());count=int(mask.sum());ece+=count/len(y)*abs(accuracy-confidence);bins.append({'lower':float(low),'n':count,'accuracy':accuracy,'confidence':confidence})
json.dump({'expected_calibration_error_10_bins':ece,'bins':bins,'per_class_recall':[float(np.mean(p[y==k]==k)) for k in range(5)],'note':'Descriptive final-test calibration; no calibration was fitted on test.'},open('reports/primary_calibration.json','w'),indent=2)
print(out)

channel_effects=[]
for model in ['bandpower','svm','mdm','ts','fbts','fbcsp','eegnet','cbramod','cbramod_partial','cbramod_full','reve_probe','reve_partial','reve_full','csbrain_probe','csbrain_partial','csbrain_full']:
 base=sorted([r for r in rows if r['model']==model and r['channels']==64 and r['montage']=='distributed' and r['seed'] in ['deterministic',2026]],key=lambda r:r['subject'])
 if len(base)!=15:continue
 for n in [32,16,8]:
  reduced=sorted([r for r in rows if r['model']==model and r['channels']==n and r['montage']=='distributed' and r['seed'] in ['deterministic',2026]],key=lambda r:r['subject'])
  if len(reduced)!=15:continue
  loss=np.array([a['metrics']['accuracy']-b['metrics']['accuracy'] for a,b in zip(base,reduced)])
  channel_effects.append({'model':model,'channels':n,'accuracy_loss_vs_64':float(loss.mean()),'paired_subject_bootstrap_ci95':np.quantile(rng.choice(loss,(10000,15),replace=True).mean(1),[.025,.975]).tolist()})
json.dump(channel_effects,open('reports/channel_effects.json','w'),indent=2)

# Amplitude strata are defined by each person's training data, never final-test accuracy.
from bci_data import cached
import joblib,torch
quality=[]
for r in primary:
 subject=r['subject'];path=Path(r['checkpoint']);state=torch.load(path,map_location='cpu',weights_only=False) if path.suffix=='.pt' else joblib.load(path);idx=state['indices']
 train,_,_=cached(subject,'train');test,y,_=cached(subject,'test',evaluation=True)
 def amplitude(x):
  x=x[:,idx,128:640];x=x-x.mean(1,keepdims=True);return np.ptp(x,axis=-1).max(1)
 threshold=float(np.quantile(amplitude(train),.95));high=amplitude(test)>threshold;pred=np.load('reports/predictions/'+r['name']+'.npz')['pred'];correct=pred==y
 quality.append({'subject':subject,'train_95th_percentile_peak_to_peak_microvolt':threshold,'high_amplitude_trials':int(high.sum()),'high_amplitude_correct':int(correct[high].sum()),'other_trials':int((~high).sum()),'other_correct':int(correct[~high].sum())})
high_n=sum(r['high_amplitude_trials'] for r in quality);other_n=sum(r['other_trials'] for r in quality)
json.dump({'definition':'Maximum peak-to-peak amplitude over selected channels after CAR, relative to within-subject training 95th percentile. Descriptive amplitude flag, not a validated artifact detector. No trials excluded from primary results.','subjects':quality,'high_amplitude_accuracy':sum(r['high_amplitude_correct'] for r in quality)/high_n if high_n else None,'high_amplitude_n':high_n,'other_accuracy':sum(r['other_correct'] for r in quality)/other_n,'other_n':other_n},open('reports/amplitude_quality_analysis.json','w'),indent=2)