1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 17
  6rng = np.random.default_rng(SEED)
  7
  8# Analytic continuous-time latent flow: a circular orbit with fixed radius.
  9# The speed parameter makes clock-rate perturbations explicit.
 10def flow(z0, t, omega=1.0):
 11    z0 = np.asarray(z0, dtype=float)
 12    r = np.linalg.norm(z0)
 13    phase = math.atan2(z0[1], z0[0])
 14    a = phase + omega * np.asarray(t)
 15    return r * np.stack([np.cos(a), np.sin(a)], axis=-1)
 16
 17def discrepancy(x, y, t, alpha, omega_x=1.0, omega_y=1.0):
 18    return np.linalg.norm(flow(x, t, omega_x) - flow(y, alpha * t, omega_y), axis=1).max()
 19
 20def optimize_affine_warp(x, y, t, grid, omega_x=1.0, omega_y=1.0):
 21    vals = np.array([discrepancy(x, y, t, a, omega_x, omega_y) for a in grid])
 22    i = int(vals.argmin())
 23    return float(vals[i]), float(grid[i])
 24
 25def make_point(r, phase):
 26    return np.array([r * np.cos(phase), r * np.sin(phase)])
 27
 28def empirical_delta(ds, q=0.10):
 29    # Largest empirical threshold with approximately q false-closeness.
 30    return float(np.quantile(np.asarray(ds), q, method='linear'))
 31
 32def main():
 33    # A moderate horizon avoids trivial saturation of a max distance at 2r.
 34    t = np.linspace(0, 2.0, 81)
 35    warp_grid = np.linspace(0.55, 1.05, 1001)
 36    rows = {}
 37
 38    # Prediction 1: if y runs at omega_y, the optimal warp slope is
 39    # alpha*=omega_x/omega_y and removes clock mismatch when geometry agrees.
 40    mismatch = np.linspace(0.0, 0.35, 8)
 41    p1 = []
 42    x = make_point(1.0, 0.0)
 43    for dm in mismatch:
 44        omega_y = 1.0 + dm
 45        y = x.copy()
 46        ordinary = discrepancy(x, y, t, 1.0, omega_y=omega_y)
 47        warped, alpha = optimize_affine_warp(x, y, t, warp_grid, omega_y=omega_y)
 48        predicted_alpha = 1.0 / omega_y
 49        p1.append({'mismatch': float(dm), 'ordinary_D': float(ordinary),
 50                   'warped_D': warped, 'alpha_star': alpha,
 51                   'predicted_alpha': predicted_alpha})
 52
 53    # Prediction 2: distinct invariant radii cannot be removed by a time warp;
 54    # D >= |r_x-r_y| and should increase with the radial gap.
 55    gaps = np.linspace(0.02, 0.40, 8)
 56    p2 = []
 57    for g in gaps:
 58        x = make_point(1.0, 0.2)
 59        y = make_point(1.0 + g, 0.2)
 60        D, alpha = optimize_affine_warp(x, y, t, warp_grid)
 61        p2.append({'radius_gap': float(g), 'D': D, 'lower_bound': float(g), 'alpha_star': alpha})
 62
 63    # Prediction 3: L_exp=(max(0,m-D))^2 activates exactly for D<m;
 64    # therefore activation probability follows the empirical CDF and crosses
 65    # q near delta_emp(q).
 66    negatives = []
 67    for _ in range(300):
 68        x = make_point(rng.uniform(0.85, 1.15), rng.uniform(-math.pi, math.pi))
 69        y = make_point(rng.uniform(0.65, 1.35), rng.uniform(-math.pi, math.pi))
 70        negatives.append(optimize_affine_warp(x, y, t, warp_grid)[0])
 71    negatives = np.asarray(negatives)
 72    q = 0.10
 73    delta = empirical_delta(negatives, q)
 74    margins = np.array([0.05, 0.15, 0.25, 0.35, 0.50])
 75    activation = [float(np.mean(negatives < m)) for m in margins]
 76
 77    # Tiny controlled comparison on clock-perturbed copies: ordinary pointwise
 78    # separation treats speed change as a difference, whereas the warp-aware D
 79    # largely removes it. This is the intended benefit, not a trained model.
 80    ordinary_clock, warped_clock = [], []
 81    for _ in range(100):
 82        r = rng.uniform(.9, 1.1)
 83        phase = rng.uniform(-math.pi, math.pi)
 84        dm = rng.uniform(0.0, .35)
 85        x = make_point(r, phase)
 86        y = x.copy()
 87        ordinary_clock.append(discrepancy(x, y, t, 1.0, omega_y=1.0 + dm))
 88        warped_clock.append(optimize_affine_warp(x, y, t, warp_grid, omega_y=1.0 + dm)[0])
 89
 90    alpha_err = [abs(x['alpha_star'] - x['predicted_alpha']) for x in p1]
 91    result = {
 92      'seed': SEED,
 93      'prediction_1_clock_mismatch': p1,
 94      'prediction_2_radius_gap': p2,
 95      'prediction_3_margin_transition': {
 96          'delta_emp_q0.10': delta, 'margins': margins.tolist(),
 97          'activation_rates': activation, 'observed_rate_at_delta': float(np.mean(negatives < delta))},
 98      'mini_comparison': {
 99          'clock_pairs_ordinary_mean_D': float(np.mean(ordinary_clock)),
100          'clock_pairs_warped_mean_D': float(np.mean(warped_clock)),
101          'relative_reduction': float(1 - np.mean(warped_clock) / np.mean(ordinary_clock))},
102      'checks': {
103        'clock_warp_reduces_at_max_mismatch': bool(p1[-1]['warped_D'] < p1[-1]['ordinary_D']),
104        'alpha_matches_inverse_speed_max_error': float(max(alpha_err)),
105        'radius_lower_bound': bool(all(x['D'] + 1e-8 >= x['lower_bound'] for x in p2)),
106        'radius_monotone': bool(np.all(np.diff([x['D'] for x in p2]) > 0)),
107        'activation_monotone': bool(np.all(np.diff(activation) >= 0)),
108        'delta_is_quantile': bool(abs(np.mean(negatives < delta) - q) < 0.02)}
109    }
110    Path('results.json').write_text(json.dumps(result, indent=2))
111    print(json.dumps(result, indent=2))
112
113if __name__ == '__main__':
114    main()