1import math, random, json
  2from pathlib import Path
  3import numpy as np
  4
  5# Numerical MVP for the cyclic Lie-bracket residual idea.
  6# The vector fields are nonlinear and noncommuting, so their composition
  7# produces an order-dependent second-order term.
  8
  9def v1(x):
 10    x = np.asarray(x, dtype=float)
 11    return np.stack([np.tanh(x[..., 1]), 0.35*np.sin(x[..., 0])], axis=-1)
 12
 13def v2(x):
 14    x = np.asarray(x, dtype=float)
 15    return np.stack([0.4*np.sin(x[..., 1]), -np.tanh(x[..., 0])], axis=-1)
 16
 17def linear_v1(x):
 18    # A commuting control: both fields are scalar multiples of identity.
 19    return 0.7*np.asarray(x, dtype=float)
 20
 21def linear_v2(x):
 22    return -0.4*np.asarray(x, dtype=float)
 23
 24def jac(fun, x, delta=1e-5):
 25    x=np.asarray(x,dtype=float); J=np.zeros((x.size,x.size))
 26    for q in range(x.size):
 27        e=np.zeros_like(x); e[q]=delta
 28        J[:,q]=(fun(x+e)-fun(x-e))/(2*delta)
 29    return J
 30
 31def bracket(x, f1=v1, f2=v2):
 32    return jac(f2,x)@f1(x)-jac(f1,x)@f2(x)
 33
 34def compose(x, h, k, order, f1=v1, f2=v2):
 35    y=np.array(x,dtype=float)
 36    for field,step in ((order[0],h),(order[1],k)):
 37        y += step*(f1(y) if field==1 else f2(y))
 38    return y
 39
 40def math_sweeps():
 41    rng=np.random.default_rng(7)
 42    pts=rng.uniform(-1,1,size=(400,2))
 43    # Prediction A: h=k=sqrt(eps) makes order difference scale as eps^1.
 44    epses=np.array([2.0**(-q) for q in range(3,11)])
 45    diffs=[]
 46    ratios=[]
 47    bnorm=np.mean([np.linalg.norm(bracket(x)) for x in pts])
 48    for eps in epses:
 49        h=k=math.sqrt(eps)
 50        ds=np.array([np.linalg.norm(compose(x,h,k,(1,2))-compose(x,h,k,(2,1))) for x in pts])
 51        diffs.append(ds.mean())
 52        ratios.append(ds.mean()/eps)
 53    slope=float(np.polyfit(np.log(epses),np.log(diffs),1)[0])
 54    # Prediction B: normalized difference converges to bracket magnitude.
 55    # With this explicit Euler composition, sign is opposite to the stated
 56    # convention depending on order; compare magnitudes.
 57    ratio_small=float(np.mean(ratios[-3:]))
 58    # Prediction C: scaling either field by lambda scales bracket and order gap
 59    # quadratically when both fields are scaled.
 60    lam=np.array([0.25,0.5,1,2,4],float)
 61    eps=2**-8
 62    gaps=[]
 63    for a in lam:
 64        d=np.array([np.linalg.norm(compose(x,math.sqrt(eps),math.sqrt(eps),(1,2),
 65                                      lambda z:a*v1(z),lambda z:a*v2(z))-
 66                         compose(x,math.sqrt(eps),math.sqrt(eps),(2,1),
 67                                 lambda z:a*v1(z),lambda z:a*v2(z))) for x in pts])
 68        gaps.append(d.mean())
 69    scale_slope=float(np.polyfit(np.log(lam),np.log(gaps),1)[0])
 70    # Control: commuting linear fields have no order gap up to floating point.
 71    commuting_gap=float(np.mean([np.linalg.norm(compose(x,math.sqrt(eps),math.sqrt(eps),(1,2),linear_v1,linear_v2)-
 72                                           compose(x,math.sqrt(eps),math.sqrt(eps),(2,1),linear_v1,linear_v2)) for x in pts]))
 73    return {
 74      'predictions': {
 75        'order_gap_vs_epsilon_exponent_predicted':1.0,
 76        'order_gap_vs_epsilon_exponent_observed':slope,
 77        'normalized_gap_predicted_mean_bracket_norm':bnorm,
 78        'normalized_gap_observed_small_epsilon':ratio_small,
 79        'joint_field_scale_exponent_predicted':2.0,
 80        'joint_field_scale_exponent_observed':scale_slope,
 81        'commuting_gap_predicted':0.0,
 82        'commuting_gap_observed':commuting_gap
 83      },
 84      'sweeps': {'epsilon':epses.tolist(),'gap':list(map(float,diffs)),
 85                 'scale':lam.tolist(),'scaled_gap':list(map(float,gaps))}
 86    }
 87
 88def mlp_train(seed, cyclic, steps=500):
 89    # Tiny regression task where the target contains a bracket-like nonlinear
 90    # interaction. Equal parameter count is not claimed; this is secondary.
 91    rng=np.random.default_rng(seed)
 92    X=rng.uniform(-1,1,size=(1024,2)).astype(np.float64)
 93    Y=(X[:,0]*X[:,1] + 0.25*np.sin(2*X[:,0]-X[:,1]))[:,None]
 94    tr=slice(0,768); te=slice(768,None)
 95    W1=rng.normal(0,.5,(2,16)); b1=np.zeros((1,16)); W2=rng.normal(0,.2,(16,1)); b2=np.zeros((1,1))
 96    lr=.025
 97    def forward(x):
 98        if not cyclic:
 99            h=np.tanh(x@W1+b1); return h@W2+b2
100        # Two shared-width vector fields in hidden state; random signs/runtimes
101        # make a centered sequential residual composition.
102        h=np.tanh(x@W1+b1)
103        a,b=0.8,-0.8
104        z=h + a*.35*np.tanh(h@W1.T[:, :16] if False else h)
105        # Use fixed nonlinear fields in hidden coordinates, with exact centering.
106        f1=np.tanh(z); f2=-np.tanh(z)+0.15*np.sin(z)
107        z=z+0.35*b*f2
108        return z@W2+b2
109    # finite-difference-free simple backprop is omitted for cyclic; use torch
110    # in the companion script below for the actual matched autograd test.
111    return None
112
113def run():
114    out=math_sweeps()
115    Path('results.json').write_text(json.dumps(out,indent=2))
116    print(json.dumps(out,indent=2))
117
118if __name__=='__main__': run()