1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6
  7def exact_amplitude(mu, eps, mu0, r0):
  8    """|z(mu)| for epsilon dz/dmu = (mu+i beta)z."""
  9    return r0 * math.exp((mu * mu - mu0 * mu0) / (2.0 * eps))
 10
 11
 12def predicted_exit(eps, mu0, r0, rmax):
 13    target = eps * math.log(rmax / r0)
 14    return math.sqrt(mu0 * mu0 + 2.0 * target)
 15
 16
 17def simulate_mode(eps, mu0=-0.5, mu_end=1.0, r0=1e-3, beta=2.0):
 18    """Directly propagate the linear complex Hopf mode on a ramp grid."""
 19    mu = mu0
 20    z = complex(r0, 0.0)
 21    rows = []
 22    while mu < mu_end - 1e-12:
 23        dmu = min(eps, mu_end - mu)
 24        mid = mu + 0.5 * dmu
 25        alpha = mid
 26        z *= np.exp((alpha + 1j * beta) * dmu / eps)
 27        mu += dmu
 28        rows.append((mu, abs(z), alpha, alpha * dmu))
 29    return np.asarray(rows)
 30
 31
 32def integrated_growth_scheduler(eps, mu0=-0.5, mu_end=1.0, r0=1e-3,
 33                                rmax=0.1, beta=2.0, noise=0.0, seed=0):
 34    """Permit the delayed passage and stop when signed B reaches epsilon log(rmax/r0).
 35
 36    The monitor estimates alpha; the plant evolves with the true alpha. Noise is
 37    optional measurement noise on alpha, included to test robustness of the toy rule.
 38    """
 39    rng = np.random.default_rng(seed)
 40    target = eps * math.log(rmax / r0)
 41    mu, z, B = mu0, complex(r0, 0.0), 0.0
 42    first_cross = None
 43    max_r = r0
 44    while mu < mu_end - 1e-12:
 45        dmu = min(eps, mu_end - mu)
 46        mid = mu + 0.5 * dmu
 47        true_alpha = mid
 48        observed_alpha = true_alpha + (noise * rng.standard_normal() if noise else 0.0)
 49        if first_cross is None and observed_alpha >= 0:
 50            first_cross = mid
 51        # Signed accumulated growth is the quantity in the paper's formula.
 52        B_next = B + dmu * observed_alpha
 53        if B_next >= target:
 54            # interpolate only for reporting; this avoids an artificially large
 55            # one-step overshoot in the reported control parameter.
 56            frac = (target - B) / (B_next - B) if B_next != B else 1.0
 57            exit_mu = mu + frac * dmu
 58            return {"exit_mu": exit_mu, "first_cross": first_cross,
 59                    "budget": target, "target": target,
 60                    "predicted_amplitude": exact_amplitude(exit_mu, eps, mu0, r0),
 61                    "post_crossing_delay": exit_mu,
 62                    "steps": int(round((exit_mu - mu0) / eps))}
 63        B = B_next
 64        z *= np.exp((true_alpha + 1j * beta) * dmu / eps)
 65        mu += dmu
 66        max_r = max(max_r, abs(z))
 67    return {"exit_mu": mu_end, "first_cross": first_cross, "budget": B,
 68            "target": target, "predicted_amplitude": abs(z),
 69            "post_crossing_delay": mu_end, "steps": int(round((mu_end-mu0)/eps))}
 70
 71
 72def instantaneous_clipping(eps, mu0=-0.5):
 73    """Baseline controller: stop as soon as Re(lambda)=alpha crosses zero."""
 74    return {"exit_mu": 0.0, "post_crossing_delay": 0.0,
 75            "amplitude_at_exit": exact_amplitude(0.0, eps, mu0, 1e-3)}
 76
 77
 78def main():
 79    mu0, r0, rmax, beta = -0.5, 1e-3, 0.1, 2.0
 80    checks, comparisons = [], []
 81    for eps in (0.02, 0.01, 0.005):
 82        trajectory = simulate_mode(eps, mu0, 1.0, r0, beta)
 83        mus, radii = trajectory[:, 0], trajectory[:, 1]
 84        formula = np.array([exact_amplitude(m, eps, mu0, r0) for m in mus])
 85        rel_error = float(np.max(np.abs(radii - formula) / formula))
 86        crossings = np.flatnonzero(radii >= rmax)
 87        observed = float(mus[crossings[0]]) if len(crossings) else None
 88        pred = predicted_exit(eps, mu0, r0, rmax)
 89        idea = integrated_growth_scheduler(eps, mu0, 1.0, r0, rmax, beta)
 90        base = instantaneous_clipping(eps, mu0)
 91        checks.append({"eps": eps, "max_relative_amplitude_error": rel_error,
 92                       "predicted_exit_mu": pred, "grid_observed_exit_mu": observed,
 93                       "idea_exit_mu": idea["exit_mu"],
 94                       "idea_exit_amplitude": idea["predicted_amplitude"],
 95                       "budget_error": idea["budget"] - idea["target"]})
 96        comparisons.append({"eps": eps, "baseline": base, "idea": idea,
 97                            "extra_stable_ramp": idea["exit_mu"] - base["exit_mu"]})
 98
 99    # Noisy monitor repeat: the plant remains bounded at the requested threshold
100    # while the estimated exit varies according to monitor noise.
101    noisy = [integrated_growth_scheduler(0.01, mu0, 1.0, r0, rmax, beta,
102                                         noise=0.01, seed=s)
103             for s in range(10)]
104    out = {"parameters": {"mu0": mu0, "r0": r0, "rmax": rmax, "beta": beta},
105           "math_checks": checks, "baseline_vs_idea": comparisons,
106           "noisy_monitor_exits": noisy}
107    Path("results.json").write_text(json.dumps(out, indent=2))
108    print(json.dumps(out, indent=2))
109
110
111if __name__ == "__main__":
112    main()