1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6
  7def rotation(phi):
  8    return np.array([[math.cos(phi), -math.sin(phi)],
  9                     [math.sin(phi), math.cos(phi)]], dtype=float)
 10
 11
 12def certificate_scalar_rotation(a, phi, theta):
 13    """Exact min_gamma ||gamma exp(-j theta) a R(phi)-I||_2.
 14
 15    A planar rotation is a complex scalar. For phase error d in (-pi,pi],
 16    the positive-gamma optimum is gamma=cos(d)/a when cos(d)>0 and the
 17    optimum boundary is gamma -> 0 otherwise.
 18    """
 19    d = (phi - theta + math.pi) % (2 * math.pi) - math.pi
 20    c = math.cos(d)
 21    if c > 0:
 22        gamma = c / a
 23        z = abs(math.sin(d))
 24    else:
 25        gamma = 0.0
 26        z = 1.0
 27    # The certificate is only meaningful for z < 1.
 28    Gamma = math.asin(min(1.0, z))
 29    return gamma, z, Gamma, d
 30
 31
 32def run():
 33    rng = np.random.default_rng(2193)
 34    rows = []
 35
 36    # Prediction 1: Gamma_theta(a R(phi)) = |wrapped(phi-theta)|
 37    phase_errors = np.linspace(-0.9, 0.9, 19)
 38    err = []
 39    for d in phase_errors:
 40        _, z, g, wrapped = certificate_scalar_rotation(1.7, d, 0.0)
 41        err.append(abs(g - abs(d)))
 42    pred1 = float(max(err))
 43
 44    # Prediction 2: phase composition is additive for products.
 45    n = 8
 46    a = 0.91
 47    phis = rng.uniform(-0.28, 0.28, size=n)
 48    P = np.eye(2)
 49    for phi in phis:
 50        P = (a * rotation(phi)) @ P
 51    direct_angle = math.atan2(P[1, 0], P[0, 0])
 52    expected_angle = ((float(np.sum(phis)) + math.pi) % (2 * math.pi)) - math.pi
 53    pred2 = abs(((direct_angle - expected_angle + math.pi) % (2 * math.pi)) - math.pi)
 54
 55    # Prediction 3: for phase sum pi, I+P loses invertibility precisely at a^N=1.
 56    # The predicted distance of the exact product SRG point from -1 equals |1-a^N|.
 57    gains = np.linspace(0.82, 1.18, 73)
 58    sigma = []
 59    margin = []
 60    for gain in gains:
 61        P = (gain * rotation(math.pi / n))
 62        Pn = np.linalg.matrix_power(P, n)
 63        sigma.append(float(np.linalg.svd(np.eye(2) + Pn, compute_uv=False)[-1]))
 64        margin.append(abs(1.0 - gain ** n))
 65    sigma = np.asarray(sigma)
 66    margin = np.asarray(margin)
 67    pred3 = float(np.max(np.abs(sigma - margin)))
 68    observed_cross = float(gains[np.argmin(sigma)])
 69    predicted_cross = 1.0
 70
 71    # Prediction 4 / sharp boundary: a product with phase pi has bounded
 72    # powers below gain 1 and amplification above gain 1.
 73    horizons = np.arange(1, 13)
 74    below, above = 0.97, 1.03
 75    below_norm = np.array([below ** k for k in horizons])
 76    above_norm = np.array([above ** k for k in horizons])
 77    # Fit log growth rates; predicted slopes are log(gain).
 78    slope_below = float(np.polyfit(horizons, np.log(below_norm), 1)[0])
 79    slope_above = float(np.polyfit(horizons, np.log(above_norm), 1)[0])
 80
 81    # Small optimization demo: phase-margin penalty moves phase away from the
 82    # dangerous product angle pi, while a vanilla objective leaves it there.
 83    # Parameters are block phases; target is deliberately unstable at pi.
 84    N = 10
 85    target = np.full(N, math.pi / N)
 86    def optimize(use_penalty, steps=250, lr=0.08, lam=3.0):
 87        x = target.copy()
 88        for _ in range(steps):
 89            # synthetic fitting pressure toward the dangerous target
 90            grad = 2.0 * (x - target)
 91            total = float(np.sum(x))
 92            # keep total phase at least 0.25 rad from pi (a differentiable
 93            # hinge proxy for the estimated -1 exclusion margin)
 94            if use_penalty:
 95                dist = abs(((total - math.pi + math.pi) % (2 * math.pi)) - math.pi)
 96                if dist < 0.25:
 97                    sign = 1.0 if total >= math.pi else -1.0
 98                    grad += lam * 2.0 * (0.25 - dist) * sign
 99            x -= lr * grad
100        P = np.linalg.matrix_power(rotation(float(np.mean(x))), N)
101        return x, float(np.linalg.svd(np.eye(2) + P, compute_uv=False)[-1]), float(abs(1 + np.exp(1j * np.sum(x))))
102    x0, s0, m0 = optimize(False)
103    x1, s1, m1 = optimize(True)
104
105    result = {
106        "certificate_phase_max_abs_error": pred1,
107        "composition_angle_abs_error": pred2,
108        "invertibility_margin_max_abs_error": pred3,
109        "invertibility_crossing_observed": observed_cross,
110        "invertibility_crossing_predicted": predicted_cross,
111        "power_log_slope_below_observed_predicted": [slope_below, math.log(below)],
112        "power_log_slope_above_observed_predicted": [slope_above, math.log(above)],
113        "optimization_baseline_sigma_min": s0,
114        "optimization_phase_margin_sigma_min": s1,
115        "optimization_baseline_margin": m0,
116        "optimization_phase_margin_margin": m1,
117        "n_blocks": n,
118        "gain_sweep": {"gains": gains.tolist(), "sigma_min": sigma.tolist(), "exact_margin": margin.tolist()},
119        "phase_sweep": {"errors": phase_errors.tolist()}
120    }
121    Path("results.json").write_text(json.dumps(result, indent=2))
122    print(json.dumps(result, indent=2))
123
124
125if __name__ == "__main__":
126    run()