analyze.py

1#!/usr/bin/env python3 2"""Aggregate fixed benchmark outcomes and audit persisted numerical artifacts.""" 3from __future__ import annotations 4import json, math 5from pathlib import Path 6import numpy as np 7import pandas as pd 8ROOT=Path(__file__).resolve().parents[1] 9records=[json.loads(p.read_text()) for p in sorted((ROOT/'results').glob('n*.json'))] 10if len(records)!=60:raise RuntimeError(f'Expected 60 instances, found {len(records)}') 11rows=[];trials=[] 12for r in records: 13 flat={k:v for k,v in r.items() if not isinstance(v,(dict,list))} 14 for alg in ('queyranne','qaoa_p1','qaoa_p2'): 15 for k,v in r[alg].items(): 16 if isinstance(v,(int,float,str,bool)):flat[f'{alg}_{k}']=v 17 for alg in ('SA','Uniform','QAOA-p1','QAOA-p2'): 18 ar=[a for a in r['trials'] if a['algorithm']==alg] 19 flat[f'{alg}_hit_rate']=float(np.mean([x['hit'] for x in ar])) 20 flat[f'{alg}_mean_gap_pct']=float(np.mean([(x['J']-r['Jstar'])/r['Jstar']*100 for x in ar])) 21 flat[f'{alg}_mean_sec']=float(np.mean([x['sec'] for x in ar])) 22 for x in ar: 23 trials.append({'instance':r['instance'],'n':r['n'],'family':r['family'], 24 'rep':r['rep'],'Jstar':r['Jstar'], 25 'gap_pct':(x['J']-r['Jstar'])/r['Jstar']*100,**x}) 26 flat['queyranne_gap_pct']=(r['queyranne']['J']-r['Jstar'])/r['Jstar']*100 27 flat['uniform_theoretical_hit256']=float(-np.expm1(256*np.log1p(-r['uniform_pstar']))) 28 rows.append(flat) 29df=pd.DataFrame(rows);tdf=pd.DataFrame(trials) 30df.to_csv(ROOT/'results/instances.csv',index=False,float_format='%.17g') 31tdf.to_csv(ROOT/'results/trials.csv',index=False,float_format='%.17g') 32 33# Stratified cluster bootstrap: networks are the unit, readout batches are not 34# independent networks. Hold 4 sizes x 3 families fixed; resample 5 instances/cell. 35rng=np.random.default_rng(8890021) 36strata=[g.index.to_numpy() for _,g in df.groupby(['n','family'],sort=True)] 37indices=np.concatenate([rng.choice(ix,size=(10000,len(ix)),replace=True) for ix in strata],axis=1) 38def ci(values): 39 v=np.asarray(values,float);dist=v[indices].mean(axis=1) 40 return list(map(float,np.quantile(dist,[.025,.975]))) 41summary={} 42for alg in ('Queyranne','SA','Uniform','QAOA-p1','QAOA-p2'): 43 if alg=='Queyranne': 44 rates=df.queyranne_hit.astype(float).to_numpy();gaps=df.queyranne_gap_pct.to_numpy() 45 hits=int(rates.sum());denom=len(rates) 46 medsec=float(df.queyranne_sec.median()) 47 else: 48 rates=df[f'{alg}_hit_rate'].to_numpy();sub=tdf[tdf.algorithm==alg] 49 gaps=sub.gap_pct.to_numpy();hits=int(sub.hit.sum());denom=len(sub) 50 medsec=float(sub.sec.median()) 51 summary[alg]={'hits':hits,'attempts':denom,'hit_rate':float(rates.mean()), 52 'hit_rate_95_stratified_cluster_bootstrap':ci(rates), 53 'mean_gap_pct':float(gaps.mean()),'median_gap_pct':float(np.median(gaps)), 54 'max_gap_pct':float(gaps.max()),'median_recorded_sec':medsec} 55 56byn=[] 57for n,g in df.groupby('n',sort=True): 58 ns={'n':int(n),'instances':len(g),'partitions_each':int(g.partition_count.iloc[0]), 59 'Jstar_min':float(g.Jstar.min()),'Jstar_max':float(g.Jstar.max()), 60 'Jstar_median':float(g.Jstar.median()),'queyranne_hit_rate':float(g.queyranne_hit.mean()), 61 'queyranne_hits':int(g.queyranne_hit.sum()), 62 'exact_median_sec':float(g.exact_sec.median()), 63 'queyranne_median_sec':float(g.queyranne_sec.median()), 64 'SA_median_sec':float(tdf[(tdf.n==n)&(tdf.algorithm=='SA')].sec.median()), 65 'qaoa_p2_training_median_sec':float(g.qaoa_p2_training_sec.median()), 66 'qaoa_p2_total_median_sec':float((g.exact_sec+g.qaoa_p2_training_sec).median()), 67 'qaoa_p2_mean_pstar':float(g.qaoa_p2_pstar.mean()), 68 'qaoa_p2_mean_lift':float((g.qaoa_p2_pstar/g.uniform_pstar).mean()), 69 'uniform_theoretical_hit256':float(g.uniform_theoretical_hit256.mean()), 70 'qaoa_p2_theoretical_hit256':float(g.qaoa_p2_theoretical_hit256.mean()), 71 'queyranne_median_calls':float(g.queyranne_calls.median()), 72 'queyranne_median_unique':float(g.queyranne_unique_computations.median())} 73 for alg in ('SA','Uniform','QAOA-p1','QAOA-p2'): 74 ns[f'{alg}_hit_rate']=float(g[f'{alg}_hit_rate'].mean()) 75 ns[f'{alg}_hits']=int(tdf[(tdf.n==n)&(tdf.algorithm==alg)].hit.sum()) 76 byn.append(ns) 77 78# Post-run mathematical diagnostic, not an outcome-dependent benchmark change: 79# exhaustively scan every elementary submodularity inequality for all instances. 80witnesses=[];audit=[] 81for r in records: 82 dat=np.load(ROOT/'data'/f'{r["instance"]}.npz',allow_pickle=False) 83 n=r['n'];D=1<<(n-1);full=(1<<n)-1;J=dat['J'];star=float(J.min()) 84 assert len(J)==D-1 85 assert np.all(np.isfinite(J)) and np.min(J)>=0 86 assert abs(star-r['Jstar'])<1e-14 87 assert int(np.argmin(J))==r['exact_z'] 88 assert abs(J[r['queyranne']['z']]-r['queyranne']['J'])<1e-12 89 # f(empty)=f(full)=0, symmetric lift of the canonical table. 90 masks=np.arange(1<<n,dtype=np.int64) 91 canonical=np.where((masks&1)==1,masks,full^masks) 92 vals=np.zeros(1<<n);ok=canonical!=full;vals[ok]=J[canonical[ok]>>1] 93 lowest=0.;wit=None;tested=0 94 for i in range(n): 95 for j in range(i+1,n): 96 bases=masks[(masks&((1<<i)|(1<<j)))==0] 97 slack=vals[bases|(1<<i)]+vals[bases|(1<<j)]-vals[bases]-vals[bases|(1<<i)|(1<<j)] 98 ix=int(np.argmin(slack));tested+=len(slack) 99 if float(slack[ix])<lowest: 100 lowest=float(slack[ix]);s=int(bases[ix]) 101 wit={'S':s,'i':i,'j':j, 102 'J_S_plus_i':float(vals[s|(1<<i)]),'J_S_plus_j':float(vals[s|(1<<j)]), 103 'J_S':float(vals[s]),'J_S_plus_i_j':float(vals[s|(1<<i)|(1<<j)])} 104 witnesses.append({'instance':r['instance'],'minimum_submodular_slack':lowest, 105 'inequalities_checked':tested,'witness':wit}) 106 # Audit every saved trial and every saved readout (no additional experiments). 107 max_norm_err=0.;hits={} 108 for p in (1,2): 109 probs=dat[f'qaoa_p{p}_probabilities'] 110 max_norm_err=max(max_norm_err,abs(float(probs.sum())-1)) 111 assert np.min(probs)>=0 112 ps=float(probs[dat['optimal_z']].sum()) 113 assert abs(ps-r[f'qaoa_p{p}']['pstar'])<1e-14 114 for t in [x for x in r['trials'] if x['algorithm']==f'QAOA-p{p}']: 115 shots=dat[f'qaoa_p{p}_shots'][t['trial']];valid=shots[shots<len(J)] 116 assert len(shots)==256 117 assert abs(float(J[valid].min())-t['J'])<1e-14 118 assert bool(np.isin(shots,dat['optimal_z']).any())==t['hit'] 119 for t in r['trials']: 120 assert t['z']>=0 121 assert abs(J[t['z']]-t['J'])<1e-12 122 assert bool(J[t['z']]<=star+1e-10+1e-8*abs(star))==t['hit'] 123 audit.append({'instance':r['instance'],'pass':True,'probability_norm_abs_error':max_norm_err}) 124 125contrast={ 126 'qaoa_p2_minus_uniform_observed_hit_rate':float((df['QAOA-p2_hit_rate']-df['Uniform_hit_rate']).mean()), 127 'qaoa_p2_minus_uniform_observed_95_ci':ci(df['QAOA-p2_hit_rate']-df['Uniform_hit_rate']), 128 'SA_minus_qaoa_p2_observed_hit_rate':float((df['SA_hit_rate']-df['QAOA-p2_hit_rate']).mean()), 129 'SA_minus_qaoa_p2_observed_95_ci':ci(df['SA_hit_rate']-df['QAOA-p2_hit_rate']), 130 'qaoa_p2_mean_pstar':float(df.qaoa_p2_pstar.mean()), 131 'qaoa_p2_median_pstar':float(df.qaoa_p2_pstar.median()), 132 'qaoa_p2_median_pstar_lift_over_valid_uniform':float((df.qaoa_p2_pstar/df.uniform_pstar).median()), 133 'qaoa_p2_lift_min':float((df.qaoa_p2_pstar/df.uniform_pstar).min()), 134 'qaoa_p2_lift_max':float((df.qaoa_p2_pstar/df.uniform_pstar).max()), 135 'qaoa_p2_mean_invalid_probability':float(df.qaoa_p2_p_invalid.mean()), 136 'qaoa_p2_mean_theoretical_hit256':float(df.qaoa_p2_theoretical_hit256.mean()), 137 'uniform_mean_theoretical_hit256':float(df.uniform_theoretical_hit256.mean()), 138 'qaoa_p2_nfev_min':int(df.qaoa_p2_nfev.min()), 139 'qaoa_p2_nfev_max':int(df.qaoa_p2_nfev.max()), 140 'qaoa_p2_nfev_sum':int(df.qaoa_p2_nfev.sum()), 141 'qaoa_p1_nfev_sum':int(df.qaoa_p1_nfev.sum()), 142 'non_submodular_instances':sum(x['minimum_submodular_slack']< -1e-10 for x in witnesses), 143 'elementary_submodularity_inequalities_checked':sum(x['inequalities_checked'] for x in witnesses), 144 'max_covariance_relative_residual':float(df.covariance_relative_residual.max()), 145 'max_conditional_vs_joint_phi_difference':float(df.SI_crosscheck_max_abs_error.max()), 146 'exact_min_smaller_part_size':int(df.exact_smaller_part.min()), 147 'exact_max_smaller_part_size':int(df.exact_smaller_part.max()), 148 'ties_instances':int((df.n_tied_optima>1).sum()), 149} 150out={'summary':summary,'by_size':byn,'comparisons':contrast, 151 'totals':{'networks':60,'partitions':int(df.partition_count.sum()), 152 'stochastic_trials':len(tdf),'readout_shots_per_qaoa_depth':60*10*256}, 153 'bootstrap':{'replicates':10000,'unit':'network','strata':'node count x family','seed':8890021}} 154(ROOT/'results/summary.json').write_text(json.dumps(out,indent=2)+'\n') 155(ROOT/'validation/submodularity_witnesses.json').write_text(json.dumps(witnesses,indent=2)+'\n') 156(ROOT/'validation/persisted_artifact_audit.json').write_text(json.dumps({'passed':True,'instances':audit},indent=2)+'\n') 157print(json.dumps(out,indent=2)) 158