1"""Gram-Whitened Directional Pooling: numerical MVP and mechanism checks."""
  2import json
  3import numpy as np
  4
  5
  6def unit_rows(x):
  7    x = np.asarray(x, dtype=float)
  8    return x / np.maximum(np.linalg.norm(x, axis=1, keepdims=True), 1e-12)
  9
 10
 11def softmax_weights(theta, atoms, epsilon):
 12    theta, atoms = unit_rows(theta), unit_rows(atoms)
 13    logits = theta @ atoms.T / epsilon
 14    logits -= logits.max(axis=1, keepdims=True)
 15    e = np.exp(logits)
 16    return e / e.sum(axis=1, keepdims=True)
 17
 18
 19def gram_and_pool(theta, q, atoms, epsilon, features):
 20    """Finite quadrature implementation; features has shape [K,C]."""
 21    W = softmax_weights(theta, atoms, epsilon)
 22    q = np.asarray(q, dtype=float); q = q / q.sum()
 23    c = W.T @ (q[:, None] * features)
 24    G = W.T @ (q[:, None] * W)
 25    return W, G, c
 26
 27
 28def whiten(c, G, lam=0.0):
 29    A = (G + G.T) / 2 + lam * np.eye(G.shape[0])
 30    val, vec = np.linalg.eigh(A)
 31    inv = vec @ np.diag(1.0 / np.maximum(val, 1e-12)) @ vec.T
 32    return inv @ c
 33
 34
 35def projected_energy(c, G, lam=0.0):
 36    z = whiten(c, G, lam)
 37    return float(np.sum(c * z))
 38
 39
 40def circle(K=512):
 41    a = np.arange(K) * 2 * np.pi / K
 42    return np.stack([np.cos(a), np.sin(a)], axis=1)
 43
 44
 45def mechanism_checks(seed=7):
 46    rng = np.random.default_rng(seed)
 47    theta, q = circle(), np.ones(512) / 512
 48    angles = np.array([0.0, .18, 1.7, 3.2])
 49    atoms = np.stack([np.cos(angles), np.sin(angles)], axis=1)
 50    ang = np.arctan2(theta[:, 1], theta[:, 0])
 51    f = (1.2*np.cos(2*ang) + .5*np.sin(5*ang))[:, None]
 52    W, G, c = gram_and_pool(theta, q, atoms, .22, f)
 53    signal_energy = float(np.sum(q*f[:, 0]**2))
 54    # Prediction 1: PoU is exact up to floating point.
 55    pou_err = float(np.max(np.abs(W.sum(1)-1)))
 56    # Prediction 2: lambda=0 gives an orthogonal projection, hence energy <= ||f||^2.
 57    z = whiten(c, G)
 58    reconstruction = W @ z[:, 0]
 59    reconstruction_energy = float(np.sum(q*reconstruction**2))
 60    projection_identity_error = abs(reconstruction_energy-projected_energy(c,G))
 61    # Prediction 3: ridge energy is decreasing; at large lambda, lambda*E -> c^T c.
 62    lambdas = np.array([0, 1e-5, 1e-3, 1e-2, 1e-1, 1, 10, 100, 1000.])
 63    energies = np.array([projected_energy(c,G,x) for x in lambdas])
 64    asymptotic_target = float(np.sum(c*c))
 65    asymptotic_observed = float(lambdas[-1]*energies[-1])
 66    # Increasing a dense atom family eventually captures the smooth signal.
 67    counts, captured = [], []
 68    for n in [2,3,4,6,8,12,16,24,32]:
 69        a = np.linspace(0, 2*np.pi, n, endpoint=False)
 70        at = np.stack([np.cos(a),np.sin(a)],1)
 71        _, gg, cc = gram_and_pool(theta,q,at,.18,f)
 72        counts.append(n); captured.append(projected_energy(cc,gg))
 73    return {
 74        'prediction_pou_max_error': pou_err,
 75        'prediction_projection_energy_le_signal': bool(projected_energy(c,G) <= signal_energy + 2e-5),
 76        'signal_energy': signal_energy,
 77        'projection_energy': projected_energy(c,G),
 78        'projection_identity_abs_error': projection_identity_error,
 79        'lambdas': lambdas.tolist(), 'ridge_energies': energies.tolist(),
 80        'prediction_ridge_monotone': bool(np.all(np.diff(energies) <= 1e-8)),
 81        'ridge_asymptotic_predicted_lambdaE': asymptotic_target,
 82        'ridge_asymptotic_observed_lambdaE_at_1000': asymptotic_observed,
 83        'atom_counts': counts, 'captured_energy': captured,
 84        'gram_condition_number': float(np.linalg.cond(G)),
 85        'energy_fraction': projected_energy(c,G)/signal_energy,
 86    }
 87
 88
 89def classification_experiment(seed=11, N=1600):
 90    """Matched nearest-centroid comparison on rotated frequency signals."""
 91    rng=np.random.default_rng(seed); K=64; theta=circle(K); q=np.ones(K)/K
 92    atoms=np.stack([np.cos(np.arange(4)*np.pi/2),np.sin(np.arange(4)*np.pi/2)],1)
 93    raw=[]; white=[]; y=[]
 94    for label in range(2):
 95        for _ in range(N//2):
 96            phase=rng.uniform(0,2*np.pi); a=np.arctan2(theta[:,1],theta[:,0])+phase
 97            sig=(np.cos((2 if label==0 else 5)*a)+.15*rng.normal(size=K))[:,None]
 98            _,G,c=gram_and_pool(theta,q,atoms,.25,sig)
 99            raw.append(c[:,0]); white.append(whiten(c,G,.01)[:,0]); y.append(label)
100    y=np.asarray(y); raw=np.asarray(raw); white=np.asarray(white)
101    order=rng.permutation(len(y)); split=int(.7*len(y)); tr,te=order[:split],order[split:]
102    def acc(A):
103        mu=np.stack([A[tr][y[tr]==j].mean(0) for j in [0,1]])
104        pred=((A[te,None]-mu[None])**2).sum(2).argmin(1)
105        return float((pred==y[te]).mean())
106    return {'raw_soft_pool_accuracy':acc(raw), 'gram_whitened_accuracy':acc(white), 'n_test':len(te)}
107
108
109def main():
110    out={'mechanism':mechanism_checks(), 'classification':classification_experiment()}
111    with open('results.json','w') as f: json.dump(out,f,indent=2)
112    print(json.dumps(out,indent=2))
113
114if __name__=='__main__': main()