1import json, sys
  2import numpy as np
  3import torch
  4
  5sys.path.insert(0, '/home/maxwelhelp/all/math2nn')
  6from bench import get_dataset, make_model, train_model, sweep_baseline, make_report
  7from bench.protocol import evaluate
  8
  9SEEDS = tuple(range(8))
 10SWEEP_SEEDS = (0, 1, 2, 3)
 11# Same union on both sides: baseline's standard knob is calibration quantile;
 12# training lr/epochs are also shared and all are evaluated on both sides.
 13GRID = [
 14    {'lr': 0.0015, 'epochs': 12, 'alpha': 0.05},
 15    {'lr': 0.0030, 'epochs': 12, 'alpha': 0.10},
 16    {'lr': 0.0060, 'epochs': 12, 'alpha': 0.15},
 17]
 18
 19
 20def fit_and_score(seed, cfg, return_sig=False):
 21    torch.manual_seed(seed); np.random.seed(seed)
 22    d = get_dataset('tabular', seed, n_train=400, n_test=400)
 23    # Fixed independent calibration split carved from training data; model fitting
 24    # uses only the first 300 points, calibration uses the remaining 100.
 25    fit = dict(d)
 26    fit['xtr'], fit['ytr'] = d['xtr'][:300], d['ytr'][:300]
 27    model = make_model('mlp_tiny', d['input_shape'], d['out_dim'])
 28    model, _, _ = train_model(model, fit, epochs=cfg['epochs'], lr=cfg['lr'], batch=128)
 29    model = model.cpu()
 30    model.eval()
 31    with torch.no_grad():
 32        pred_fit = model(d['xtr'][:300]).detach().cpu().numpy().reshape(-1)
 33        pred_cal = model(d['xtr'][300:]).detach().cpu().numpy().reshape(-1)
 34        pred_test = model(d['xte']).detach().cpu().numpy().reshape(-1)
 35    yfit = d['ytr'][:300].numpy().reshape(-1)
 36    ycal = d['ytr'][300:].numpy().reshape(-1)
 37    ytest = d['yte'].numpy().reshape(-1)
 38    mse = float(np.mean((pred_test-ytest)**2))
 39    # Absolute residual score; max/order-statistic boundary is projective.
 40    scores = np.abs(pred_cal-ycal)
 41    n = len(scores)
 42    # split conformal upper quantile with finite-sample conservative index
 43    rank = min(n-1, int(np.ceil((n+1)*(1-cfg['alpha'])))-1)
 44    threshold = float(np.partition(scores, rank)[rank])
 45    boundary = np.flatnonzero(scores >= threshold - 1e-12)
 46    # For the usual continuous case K=1 and R|K follows Beta(1,n).
 47    k = int(len(boundary))
 48    beta_q = float(1.0-cfg['alpha']**(1.0/n)) if k == 1 else float(threshold)
 49    test_risk = float(np.mean(np.abs(pred_test-ytest) > threshold))
 50    # deletion check on actual trained-model calibration residuals
 51    # (recompute order statistic, with fixed predictions; this is the rule's
 52    # explicit boundary map and does not use test labels).
 53    preserved=[]; equivalence=[]
 54    for i in range(n):
 55        rem=np.delete(scores,i); r2=min(n-2, int(np.ceil((len(rem)+1)*(1-cfg['alpha'])))-1)
 56        b2=np.flatnonzero(rem >= np.partition(rem,r2)[r2]-1e-12)
 57        b2full=np.array([j if j<i else j+1 for j in b2])
 58        same=np.array_equal(np.sort(b2full), np.sort(boundary))
 59        preserved.append(same); equivalence.append(bool(scores[i] <= threshold+1e-12)==same)
 60    sig={'threshold':threshold, 'boundary_size':k,
 61         'empirical_test_risk':test_risk,
 62         'beta_predicted_mean_risk':float(k/(n+1)) if k else None,
 63         'deletion_projectivity':float(np.mean(preserved)),
 64         'deletion_equivalence':float(np.mean(equivalence))}
 65    return (mse, sig) if return_sig else mse
 66
 67
 68def make_fn(cfg):
 69    return lambda seed: fit_and_score(seed, cfg)
 70
 71
 72def main():
 73    # Baseline is ordinary split-conformal selective regression: same trained
 74    # network, but no boundary-indexed beta diagnostic. Sweep all configs.
 75    base = sweep_baseline(make_fn, GRID, seeds=SWEEP_SEEDS)
 76    best = base['best_cfg']
 77    # evaluate baseline best on all paired seeds
 78    base_full = evaluate(make_fn(best), seeds=SEEDS)
 79    base = dict(base); base['full'] = base_full
 80    idea_cfgs = GRID  # same union, 3 settings; idea selects best on sweep seeds
 81    idea_trials=[]
 82    for cfg in idea_cfgs:
 83        r=evaluate(make_fn(cfg), seeds=SWEEP_SEEDS)
 84        idea_trials.append((r['mean'],cfg,r))
 85    _, idea_cfg, _ = min(idea_trials, key=lambda x:x[0])
 86    idea_full=evaluate(make_fn(idea_cfg), seeds=SEEDS)
 87    # independent mechanism signature from all trained models at selected cfg
 88    signatures=[fit_and_score(s, idea_cfg, True)[1] for s in SEEDS]
 89    keys=['boundary_size','empirical_test_risk','beta_predicted_mean_risk','deletion_projectivity','deletion_equivalence']
 90    avg={k:float(np.mean([x[k] for x in signatures if x[k] is not None])) for k in keys}
 91    # Quantitative prediction: continuous scalar residual boundary should be K=1,
 92    # projective, and observed risk near beta mean. Use honest tolerance.
 93    confirmed=(abs(avg['boundary_size']-1.0)<0.25 and avg['deletion_projectivity']>0.98
 94               and abs(avg['empirical_test_risk']-avg['beta_predicted_mean_risk'])<0.08)
 95    sig={'prediction':'trained scalar residual order-statistic has K=1 and projective deletion; risk approximately Beta(1,n) mean',
 96         'observed_mean':avg, 'per_seed':signatures, 'confirmed':bool(confirmed)}
 97    report=make_report('tabular','mlp_tiny',base,idea_full,{'mechanism_signature':sig,
 98        'selection':{'baseline_grid':GRID,'idea_grid':GRID,'baseline_best':best,'idea_best':idea_cfg},
 99        'structural_match':'Independent neural calibration residuals define a scalar acceptance set; explicit order-statistic boundary is the stated projective structure.'})
100    with open('bench_report.json','w') as f: json.dump(report,f,indent=2)
101    print(json.dumps(report,indent=2))
102
103if __name__=='__main__': main()