1import json, time
  2import numpy as np
  3
  4
  5def huber(r, delta):
  6    r = np.asarray(r)
  7    a = np.abs(r)
  8    return np.where(a <= delta, 0.5*r*r, delta*(a - 0.5*delta))
  9
 10
 11def huber_barycenter(Q, weights=None, delta=1.0, iterations=12, sort_output=True):
 12    """Safeguarded pointwise displacement-Huber barycenter for [objects, quantiles]."""
 13    Q = np.asarray(Q, dtype=float)
 14    n, m = Q.shape
 15    if weights is None:
 16        weights = np.ones(n) / n
 17    weights = np.asarray(weights, dtype=float)
 18    weights = weights / weights.sum()
 19    q = np.sum(weights[:, None] * Q, axis=0)
 20    lo, hi = Q.min(axis=0), Q.max(axis=0)
 21    # Newton is used where local curvature exists; otherwise bisection is safe.
 22    for _ in range(iterations):
 23        r = q[None, :] - Q
 24        score = np.sum(weights[:, None] * np.clip(r, -delta, delta), axis=0)
 25        active = (np.abs(r) < delta).astype(float)
 26        curvature = np.sum(weights[:, None] * active, axis=0)
 27        # Maintain a score bracket: score is monotone increasing in q.
 28        hi = np.where(score > 0, q, hi)
 29        lo = np.where(score < 0, q, lo)
 30        newton = q - score / np.maximum(curvature, 1e-12)
 31        bisect = 0.5 * (lo + hi)
 32        # Also reject Newton steps outside the valid bracket.
 33        q = np.where((curvature > 1e-10) & (newton >= lo) & (newton <= hi), newton, bisect)
 34    return np.sort(q) if sort_output else q
 35
 36
 37def objective(q, Q, weights, delta):
 38    return np.sum(weights[:, None] * huber(q[None, :] - Q, delta))
 39
 40
 41def mechanism_checks():
 42    delta = 1.0
 43    # Prediction 1: for [0,0,T]/3, q=T/3 until T=1.5, then q=delta/2.
 44    Ts = np.linspace(0, 8, 33)
 45    observed, predicted = [], []
 46    for T in Ts:
 47        Q = np.array([[0.0], [0.0], [T]])
 48        observed.append(huber_barycenter(Q, delta=delta, iterations=20, sort_output=False)[0])
 49        predicted.append(min(T/3.0, delta/2.0))
 50    observed, predicted = np.array(observed), np.array(predicted)
 51    transition_T = Ts[np.where(np.abs(observed - delta/2) < 2e-4)[0][0]]
 52    max_err = float(np.max(np.abs(observed-predicted)))
 53
 54    # Prediction 2: inliers within delta remain exactly quadratic/arithmetic.
 55    shifts = np.linspace(-0.9, 0.9, 19)
 56    small_err = max(abs(huber_barycenter(np.array([[-s], [0.0], [s]]), delta=delta,
 57                                         iterations=12, sort_output=False)[0]) for s in shifts)
 58
 59    # Prediction 3: one remote token's score contribution is capped at delta/n.
 60    scores = []
 61    for T in [2, 10, 100, 1000]:
 62        Q = np.array([[0.0], [0.0], [T]])
 63        q = huber_barycenter(Q, delta=delta, iterations=20, sort_output=False)[0]
 64        scores.append(abs(np.clip(q-T, -delta, delta)/3.0))
 65    return {
 66        'predicted_transition_T': 1.5, 'observed_transition_T_grid': float(transition_T),
 67        'transition_grid_step': 0.25, 'max_saturation_curve_error': max_err,
 68        'predicted_small_shift_error': 0.0, 'observed_small_shift_max_error': float(small_err),
 69        'predicted_outlier_score_limit': 1/3, 'observed_outlier_score_max': float(max(scores)),
 70        'checks_pass': bool(max_err < 0.03 and small_err < 1e-6 and max(scores) <= 1/3+1e-8)
 71    }
 72
 73
 74def classification_experiment(seed=7):
 75    rng = np.random.default_rng(seed)
 76    m, nobj, ntrain, ntest = 32, 12, 1200, 800
 77    shape = np.linspace(-0.7, 0.7, m)
 78    def make(count, corruption_rate):
 79        y = rng.integers(0, 2, count)
 80        center = (2*y-1)*0.85
 81        Q = center[:, None, None] + shape[None, None, :] + rng.normal(0, .16, (count,nobj,m))
 82        bad = rng.random((count,nobj)) < corruption_rate
 83        Q += bad[:, :, None] * rng.choice([-1, 1], (count,nobj,1)) * 5.0
 84        return Q, y
 85    trainQ, trainy = make(ntrain, 0.0)
 86    test_clean, testy = make(ntest, 0.0)
 87    test_bad, _ = make(ntest, 0.5)
 88    def pool(A, mode):
 89        out=[]; t=time.perf_counter()
 90        for x in A:
 91            out.append(x.mean(axis=0) if mode == 'mean' else huber_barycenter(x, delta=1.0, iterations=8))
 92        return np.asarray(out), time.perf_counter()-t
 93    results={}
 94    for mode in ['mean','huber']:
 95        Xtr,ttr=pool(trainQ,mode); Xc,tc=pool(test_clean,mode); Xb,tb=pool(test_bad,mode)
 96        cent=np.stack([Xtr[trainy==k].mean(0) for k in [0,1]])
 97        def acc(X,y): return float(np.mean(np.argmin(((X[:,None,:]-cent[None,:,:])**2).mean(2),1)==y))
 98        ac,ab=acc(Xc,testy),acc(Xb,testy)
 99        results[mode]={'clean_accuracy':ac,'corrupt_accuracy':ab,'degradation':ac-ab,
100                       'pool_ms_per_sample_clean':1000*tc/ntest,'pool_ms_per_sample_train':1000*ttr/ntrain}
101    results['corruption_rate']=0.5
102    results['huber_over_mean_cost_ratio']=results['huber']['pool_ms_per_sample_clean']/results['mean']['pool_ms_per_sample_clean']
103    return results
104
105
106if __name__ == '__main__':
107    report={'mechanism_checks':mechanism_checks(),'classification':classification_experiment()}
108    with open('results.json','w') as f: json.dump(report,f,indent=2)
109    print(json.dumps(report,indent=2))