1import json
  2import math
  3import random
  4import numpy as np
  5
  6# Barrier-controlled basin switching MVP.
  7# Potential U(x)=(x^2-1)^2/4 has minima +/-1 and saddle 0.
  8# Delta V = U(0)-U(+/-1)=1/4; U''(min)=2, U''(saddle)=-1.
  9# For dX=-U'(X)dt+sqrt(2 eps)dW, Kramers predicts
 10# k(eps) ~= sqrt(2)/(2*pi)*exp(-0.25/eps).
 11
 12SEED = 1306
 13DT = 0.01
 14N_PATHS = 220
 15T = 1200.0
 16STEPS = int(T / DT)
 17BARRIER = 0.25
 18PREF = math.sqrt(2.0) / (2.0 * math.pi)
 19
 20
 21def potential(x):
 22    return 0.25 * (x*x - 1.0)**2
 23
 24
 25def grad(x):
 26    return x*x*x - x
 27
 28
 29def simulate_rate(eps, seed, n=N_PATHS, t=T):
 30    """Estimate committed basin switching rate using +/-0.5 hysteresis."""
 31    rng = np.random.default_rng(seed)
 32    x = np.full(n, -1.0)
 33    state = np.full(n, -1, dtype=np.int8)
 34    switches = 0
 35    # Count only transitions after reaching the opposite committed region.
 36    for _ in range(int(t / DT)):
 37        x += -DT * grad(x) + math.sqrt(2.0 * eps * DT) * rng.standard_normal(n)
 38        new = state.copy()
 39        new[(state == -1) & (x > 0.5)] = 1
 40        new[(state == 1) & (x < -0.5)] = -1
 41        switches += int(np.count_nonzero(new != state))
 42        state = new
 43    # Each path has an opportunity to switch in either direction; this is a
 44    # conservative event rate and is comparable across noise settings.
 45    return switches / (n * t), switches
 46
 47
 48def feedback_rate(q, seed, n=N_PATHS, t=T):
 49    """Use eps=DeltaV/q while in the escape-to-target mode."""
 50    eps = np.clip(BARRIER / q, 1e-4, 0.20)
 51    return simulate_rate(float(eps), seed, n, t)
 52
 53
 54def linear_fit(x, y):
 55    return float(np.polyfit(np.asarray(x), np.asarray(y), 1)[0])
 56
 57
 58def main():
 59    random.seed(SEED)
 60    np.random.seed(SEED)
 61
 62    # Prediction 1: log transition rate against 1/eps has slope -DeltaV.
 63    eps_values = np.array([0.045, 0.055, 0.070, 0.090, 0.120, 0.160])
 64    raw = []
 65    for i, eps in enumerate(eps_values):
 66        rate, events = simulate_rate(float(eps), SEED + i)
 67        raw.append({"eps": float(eps), "rate": rate, "events": events})
 68    usable = [r for r in raw if r["events"] >= 3]
 69    slope = linear_fit([1/r["eps"] for r in usable], [math.log(r["rate"]) for r in usable])
 70    slope_error = abs(slope + BARRIER) / BARRIER
 71
 72    # Prediction 2: feedback makes the exponential part exp(-q), so
 73    # log(rate/PREF) versus q should have slope -1.
 74    q_values = np.array([1.6, 2.0, 2.5, 3.0, 3.5, 4.0])
 75    fb = []
 76    for i, q in enumerate(q_values):
 77        rate, events = feedback_rate(float(q), SEED + 100 + i)
 78        fb.append({"q": float(q), "eps": BARRIER/float(q), "rate": rate, "events": events})
 79    usable_fb = [r for r in fb if r["events"] >= 3]
 80    feedback_slope = linear_fit([r["q"] for r in usable_fb],
 81                                [math.log(r["rate"]/PREF) for r in usable_fb])
 82    feedback_error = abs(feedback_slope + 1.0)
 83
 84    # Prediction 3: changing q by delta q changes rate ratio by exp(-delta q).
 85    q_lo, q_hi = 2.0, 3.5
 86    lo = next(r for r in fb if r["q"] == q_lo)
 87    hi = next(r for r in fb if r["q"] == q_hi)
 88    observed_ratio = hi["rate"] / lo["rate"]
 89    predicted_ratio = math.exp(-(q_hi-q_lo))
 90    ratio_factor = observed_ratio / predicted_ratio if predicted_ratio else float("inf")
 91
 92    # Small optimizer comparison: equal Euler updates, beginning in the left
 93    # basin and targeting the right basin. Fixed-noise SGD is compared with
 94    # barrier feedback (escape q=2 while left, safe q=8 after reaching right).
 95    def train_controller(mode, seed, steps=120000):
 96        rng = np.random.default_rng(seed)
 97        x = -1.0
 98        reached = None
 99        switches = 0
100        for k in range(steps):
101            if mode == "fixed":
102                eps = BARRIER / 3.0
103            else:
104                # Known toy barrier estimate; the current basin is identified
105                # by sign and target is the lower-loss right basin.
106                eps = BARRIER / (2.0 if x < 0 else 8.0)
107            x += -DT * grad(x) + math.sqrt(2*eps*DT) * rng.standard_normal()
108            if reached is None and x > 0.5:
109                reached = k + 1
110            if x > 0.5 and x < 0.0:
111                switches += 1
112        return reached, x, switches
113    baseline = train_controller("fixed", SEED + 500)
114    idea = train_controller("controlled", SEED + 501)
115
116    result = {
117        "known_barrier": BARRIER,
118        "kramers_prefactor": PREF,
119        "prediction_1_rate_slope": {"predicted": -BARRIER, "observed": slope,
120                                    "relative_error": slope_error, "data": raw},
121        "prediction_2_feedback_slope": {"predicted": -1.0, "observed": feedback_slope,
122                                         "absolute_error": feedback_error, "data": fb},
123        "prediction_3_ratio": {"q_low": q_lo, "q_high": q_hi,
124                               "predicted_rate_ratio": predicted_ratio,
125                               "observed_rate_ratio": observed_ratio,
126                               "observed_over_predicted": ratio_factor},
127        "mini_optimizer": {"baseline_fixed_noise": baseline,
128                            "barrier_controlled": idea,
129                            "steps": 120000, "dt": DT},
130        "criterion": {"slope_tolerance_fraction": 0.20,
131                      "feedback_slope_tolerance": 0.20,
132                      "ratio_factor_tolerance": [0.5, 2.0]},
133        "worked": bool(slope_error <= 0.20 and feedback_error <= 0.20 and
134                        0.5 <= ratio_factor <= 2.0)
135    }
136    print(json.dumps(result, indent=2))
137
138if __name__ == "__main__":
139    main()