1import json, math
  2from pathlib import Path
  3import numpy as np
  4from scipy.optimize import nnls
  5
  6SEED = 2184
  7rng = np.random.default_rng(SEED)
  8
  9
 10def lasso_cd(A, b, lam=0.04, iters=1200):
 11    # objective 0.5||Ac-b||^2 + lam||c||_1, deterministic coordinate descent
 12    A = np.asarray(A, float); b = np.asarray(b, float)
 13    L = A.shape[1]; c = np.zeros(L)
 14    col = (A*A).sum(0) + 1e-12
 15    for _ in range(iters):
 16        old = c.copy()
 17        for k in range(L):
 18            rho = A[:, k] @ (b - A @ c + A[:, k] * c[k])
 19            z = abs(rho) - lam
 20            c[k] = np.sign(rho) * max(z, 0.) / col[k]
 21        if np.max(abs(c-old)) < 1e-9: break
 22    return c
 23
 24
 25def weak_patch_integral(n, sigma=1.0, repeats=300):
 26    # Integral of a smooth function with iid observation noise; Monte Carlo weak integration.
 27    vals=[]
 28    for _ in range(repeats):
 29        x = rng.random(n)
 30        phi = 6*x*(1-x) # compact-like bump on patch, normalized below
 31        truth = np.sin(2*np.pi*x) + .35*x
 32        vals.append(np.mean(phi*(truth + sigma*rng.normal(size=n))) / np.mean(phi))
 33    return np.std(vals, ddof=1)
 34
 35
 36def consensus_error(p, m, trials=30000):
 37    # binary local support decisions, modal support with random tie breaking
 38    x = rng.random((trials, m)) < p
 39    votes = x.sum(1)
 40    return np.mean(votes < (m/2)) if m % 2 else np.mean(votes <= (m//2-1))
 41
 42
 43def patch_regression_demo(npatch=80, q=24, noise=.22, lam=.045):
 44    # Two candidate weak operators; true support is term 0. Each local patch gets q
 45    # quadrature samples and noisy b. LASSO support is then pooled by modal consensus.
 46    local = []
 47    for j in range(npatch):
 48        # normalized test-function integrations create mildly varying local designs
 49        A = rng.normal(size=(3,2)); A[:,0] += 1.4
 50        ctrue = np.array([1.0, 0.0])
 51        b = A @ ctrue + noise*rng.normal(size=3)/math.sqrt(q)
 52        c = lasso_cd(A, b, lam)
 53        local.append(c)
 54    local=np.array(local)
 55    supports=np.abs(local)>0.16
 56    true_support=np.array([True, False])
 57    local_exact=np.all(supports == true_support[None, :], axis=1)
 58    modal_support=np.sum(supports, axis=0) >= np.ceil(npatch/2)
 59    modal_exact=bool(np.all(modal_support == true_support))
 60    return float(np.mean(local_exact)), modal_exact, local
 61
 62
 63def main():
 64    global rng
 65    # Prediction 1: weak integration averages independent noise, std proportional n^-1/2.
 66    ns=np.array([8,16,32,64,128,256])
 67    stds=np.array([weak_patch_integral(int(n)) for n in ns])
 68    slope=float(np.polyfit(np.log(ns), np.log(stds), 1)[0])
 69
 70    # Prediction 2: local decisions only improve consensus when p > 1/2.
 71    ps=[.40,.55,.70,.85]
 72    ms=[3,7,15,31]
 73    rows=[]
 74    for p in ps:
 75        for m in ms:
 76            obs=consensus_error(p,m)
 77            bound=math.exp(-2*m*(p-.5)**2)
 78            rows.append({'p':p,'m':m,'observed_error':float(obs),'hoeffding_bound':float(bound)})
 79
 80    # Prediction 3: with a fixed local accuracy p=.7, more patches reduce errors;
 81    # check monotonicity and compare the observed trend to the theoretical bound.
 82    fixed=[r for r in rows if r['p']==.70]
 83    monotonic=all(fixed[i]['observed_error'] >= fixed[i+1]['observed_error']-0.004 for i in range(len(fixed)-1))
 84    below_bound=all(r['observed_error'] <= r['hoeffding_bound']+0.015 for r in rows if r['p']>.5)
 85
 86    # Secondary MVP comparison: local sparse fits versus modal support on noisy patches.
 87    rng=np.random.default_rng(SEED+9)
 88    reps=[]
 89    for _ in range(100):
 90        frac, modal, local=patch_regression_demo()
 91        reps.append((frac,modal))
 92    local_support=float(np.mean([x[0] for x in reps]))
 93    consensus_support_accuracy=float(np.mean([x[1] for x in reps]))
 94
 95    result={
 96      'seed':SEED,
 97      'prediction_checks':{
 98        'weak_noise_scaling':{'predicted_loglog_slope':-0.5,'observed_loglog_slope':slope,'tolerance':0.12,'passed':abs(slope+0.5)<0.12,'n':ns.tolist(),'std':stds.tolist()},
 99        'consensus_transition':{'prediction':'majority improves for p>0.5 and worsens for p<0.5','rows':rows,'passed':(rows[0]['observed_error']>0.5 and rows[4]['observed_error']<0.5 and rows[8]['observed_error']<0.5)},
100        'patch_count_scaling':{'prediction':'error decreases with m for p=.70 and stays below Hoeffding upper bound','rows':fixed,'monotonic_with_sampling_tolerance':monotonic,'below_bound_with_sampling_tolerance':below_bound,'passed':monotonic and below_bound}
101      },
102      'secondary_mvp':{'local_exact_support_accuracy_mean':local_support,'modal_consensus_exact_accuracy':consensus_support_accuracy,'interpretation':'consensus converts noisy local support votes into a region-level decision'},
103      'all_mechanism_checks_passed':bool(abs(slope+0.5)<0.12 and rows[0]['observed_error']>0.5 and rows[4]['observed_error']<0.5 and monotonic and below_bound)
104    }
105    Path('results.json').write_text(json.dumps(result,indent=2))
106    print(json.dumps(result,indent=2))
107
108if __name__=='__main__': main()