1import json
  2import numpy as np
  3
  4
  5def build_q(h, q, theta):
  6    h = max(float(h), 1e-6)
  7    theta = max(float(theta), 1e-6)
  8    c = np.sqrt(h * theta)
  9    Q = np.array([
 10        [0, 0, -q/(2*theta), c/theta, -c/theta],
 11        [0, 0, 0, 1, 1],
 12        [0, 1, 0, 0, 0],
 13        [0, 0, 1, 0, 0],
 14        [1, 0, 0, 0, 0]], dtype=float)
 15    # This is the algebraic inverse of the displayed Q.  Solving the
 16    # first/fourth rows gives a PLUS q/(4c) term in rows 4 and 5.
 17    Qi = np.array([
 18        [0, 0, 0, 0, 1],
 19        [0, 0, 1, 0, 0],
 20        [0, 0, 0, 1, 0],
 21        [theta/(2*c), .5, 0, q/(4*c), 0],
 22        [-theta/(2*c), .5, 0, -q/(4*c), 0]], dtype=float)
 23    return Q, Qi
 24
 25
 26def displayed_inverse(h, q, theta):
 27    h = max(float(h), 1e-6); theta = max(float(theta), 1e-6)
 28    c = np.sqrt(h * theta)
 29    return np.array([
 30        [0, 0, 0, 0, 1], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0],
 31        [theta/(2*c), .5, 0, -q/(4*c), 0],
 32        [-theta/(2*c), .5, 0, -q/(4*c), 0]], dtype=float)
 33
 34
 35def characteristic_gate(a, alpha=2.0, eps=1e-8):
 36    left, right = a[:-1], a[1:]
 37    s = np.abs(right-left) / (np.abs(right)+np.abs(left)+eps)
 38    g = 1.0 / (1.0 + alpha*s)
 39    return left * g, g
 40
 41
 42def scalar_tv_gate(a, alpha=2.0, eps=1e-8):
 43    left, right = a[:-1], a[1:]
 44    num = np.linalg.norm(right-left, axis=1)
 45    den = np.linalg.norm(right, axis=1)+np.linalg.norm(left, axis=1)+eps
 46    g = 1.0/(1.0+alpha*num/den)
 47    return left*g[:, None], np.broadcast_to(g[:, None], left.shape)
 48
 49
 50def metrics(pred, gates, target, base):
 51    acoustic_rms = float(np.sqrt(np.mean(pred[:,3:]**2)))
 52    base_acoustic = float(np.sqrt(np.mean(base[:,3:]**2)))
 53    return {
 54        'acoustic_rms': acoustic_rms,
 55        'acoustic_suppression_fraction': 1-acoustic_rms/base_acoustic,
 56        'equilibrium_mse_to_clean': float(np.mean((pred[:,:3]-target[:,:3])**2)),
 57        'total_mse_to_clean': float(np.mean((pred-target)**2)),
 58        'mean_gate_equilibrium': float(np.mean(gates[:,:3])),
 59        'mean_gate_acoustic': float(np.mean(gates[:,3:])),
 60        'min_gate': float(np.min(gates)), 'max_gate': float(np.max(gates)),
 61    }
 62
 63
 64def main():
 65    rng = np.random.default_rng(631)
 66    identity_errors = []
 67    literal_errors = []
 68    for _ in range(1000):
 69        h = 10**rng.uniform(-6, 1); theta = 10**rng.uniform(-6, 1)
 70        q = rng.normal() * np.sqrt(h*theta)
 71        Q, Qi = build_q(h, q, theta)
 72        identity_errors.append(max(np.max(np.abs(Q@Qi-np.eye(5))), np.max(np.abs(Qi@Q-np.eye(5)))))
 73        Qil = displayed_inverse(h, q, theta)
 74        literal_errors.append(max(np.max(np.abs(Q@Qil-np.eye(5))), np.max(np.abs(Qil@Q-np.eye(5)))))
 75
 76    n = 256; x = np.linspace(-1, 1, n)
 77    shock = .5*(1+np.tanh(x/.045))
 78    a = np.zeros((n,5))
 79    a[:,0] = 1.0 + .4*shock
 80    a[:,1] = .15*np.sin(2*np.pi*x)
 81    a[:,2] = .08*np.cos(np.pi*x)
 82    envelope = np.exp(-(x/.22)**2)
 83    ringing = envelope*((-1.0)**np.arange(n))
 84    a[:,3] = .18*ringing; a[:,4] = -.15*ringing
 85    a += rng.normal(0, .002, a.shape)
 86    clean = a.copy(); clean[:,3:] = 0.0
 87    base, base_g = a[:-1], np.ones((n-1,5)); target = clean[:-1]
 88    gated, g_char = characteristic_gate(a)
 89    scalar, g_scalar = scalar_tv_gate(a)
 90
 91    Q, Qi = build_q(2.0, .7, 1.3); r = rng.normal(size=5)
 92    result = {
 93        'math': {
 94            'corrected_inverse_max_identity_error': float(max(identity_errors)),
 95            'corrected_inverse_median_identity_error': float(np.median(identity_errors)),
 96            'literal_displayed_inverse_max_identity_error': float(max(literal_errors)),
 97            'literal_displayed_inverse_median_identity_error': float(np.median(literal_errors)),
 98            'representative_roundtrip_error': float(np.max(np.abs(Q @ (Qi @ r) - r))),
 99            'gate_min': float(np.min(g_char)), 'gate_max': float(np.max(g_char)),
100            'all_gates_bounded': bool(np.all((g_char > 0) & (g_char <= 1))),
101        },
102        'baseline_no_gate': metrics(base, base_g, target, base),
103        'characteristic_gate': metrics(gated, g_char, target, base),
104        'scalar_tv_gate': metrics(scalar, g_scalar, target, base),
105        'setup': {'n': n, 'alpha': 2.0, 'seed': 631,
106                  'description': 'smooth equilibrium jump plus localized alternating acoustic ringing'}
107    }
108    print(json.dumps(result, indent=2))
109
110if __name__ == '__main__':
111    main()