1import json, math, os
  2import numpy as np
  3
  4SEED = 2803
  5rng = np.random.default_rng(SEED)
  6
  7
  8def linear_path(gains, choices):
  9    d = 1.0
 10    out = [d]
 11    for a, i in zip(gains, choices):
 12        d *= a[i]
 13        out.append(d)
 14    return np.asarray(out)
 15
 16
 17def mixed_boundary_sweep():
 18    # Two scalar candidate maps: one contractive and one expansive.
 19    lo, hi = 0.60, 1.40
 20    ps = np.linspace(0.0, 1.0, 21)
 21    n = 20000
 22    horizon = 30
 23    rows = []
 24    for p in ps:
 25        # candidate 1 (hi) is selected with probability p
 26        choices = rng.random((n, horizon)) < p
 27        # Vectorized path distances; this is also a direct Monte Carlo estimate.
 28        logs = np.where(choices, math.log(hi), math.log(lo))
 29        logd = np.cumsum(logs, axis=1)
 30        # log of E[d_t], stabilized, equals t*log((1-p)lo+p*hi) theoretically.
 31        t = np.arange(1, horizon + 1)
 32        log_mean_d = np.array([np.logaddexp.reduce(logd[:, j]) - math.log(n) for j in range(horizon)])
 33        slope = np.polyfit(t[-15:], log_mean_d[-15:], 1)[0]
 34        predicted = math.log((1-p)*lo + p*hi)
 35        rows.append({"p_hi": float(p), "predicted_log_exponent": predicted,
 36                     "observed_log_exponent": float(slope),
 37                     "predicted_mean_gain": (1-p)*lo+p*hi})
 38    # The expected-product boundary is mean gain = 1.
 39    predicted_p = (1-lo)/(hi-lo)
 40    observed_p = min(rows, key=lambda r: abs(r["observed_log_exponent"]))["p_hi"]
 41    return rows, predicted_p, observed_p
 42
 43
 44def expected_product_check():
 45    # Time-inhomogeneous independent routing: E product is product of expected gains.
 46    probs = np.array([.15, .75, .35, .60, .20, .80, .45, .55])
 47    gains = np.array([[.55, 1.35], [.70, 1.20], [.50, 1.50], [.80, 1.10],
 48                      [.60, 1.30], [.75, 1.25], [.55, 1.40], [.65, 1.15]])
 49    n = 300000
 50    choices = rng.random((n, len(probs))) < probs[None, :]
 51    paths = np.where(choices, gains[:, 1], gains[:, 0])
 52    empirical = np.mean(np.prod(paths, axis=1))
 53    predicted = float(np.prod((1-probs)*gains[:, 0] + probs*gains[:, 1]))
 54    relerr = abs(empirical-predicted)/predicted
 55    # Also verify the finite-horizon geometric rate for a stationary contractive mix.
 56    m = .94
 57    horizon = 12
 58    stationary_choices = rng.random((n, horizon)) < .4
 59    stationary_paths = np.where(stationary_choices, 1.3, .7)
 60    empirical_h = np.mean(np.prod(stationary_paths, axis=1))
 61    predicted_h = .94**horizon
 62    return {"predicted_product": predicted, "empirical_product": float(empirical),
 63            "relative_error": float(relerr), "stationary_horizon": horizon,
 64            "predicted_geometric_product": predicted_h,
 65            "empirical_geometric_product": float(empirical_h)}
 66
 67
 68def varying_schedule_check():
 69    # A deterministic time-varying schedule tests the product of p_t-weighted gains.
 70    p = np.array([.1,.7,.2,.8,.3,.6,.4,.5,.25,.75])
 71    lo = np.array([.55,.65,.60,.70,.58,.62,.56,.68,.61,.59])
 72    hi = np.array([1.20,1.10,1.30,1.05,1.25,1.15,1.18,1.08,1.22,1.12])
 73    means = (1-p)*lo+p*hi
 74    predicted = np.cumprod(means)
 75    n = 250000
 76    choices = rng.random((n, len(p))) < p
 77    products = np.cumprod(np.where(choices, hi, lo), axis=1)
 78    observed = np.mean(products, axis=0)
 79    relerr = np.max(np.abs(observed-predicted)/(predicted+1e-12))
 80    return {"predicted_final_product": float(predicted[-1]),
 81            "observed_final_product": float(observed[-1]),
 82            "max_relative_error_all_times": float(relerr),
 83            "mean_log_gain": float(np.mean(np.log(means)))}
 84
 85
 86def recurrent_stability_comparison():
 87    # Same nonlinear tanh setup: unconstrained baseline versus random candidate maps.
 88    # Distances are measured on paired trajectories with identical inputs and routing.
 89    local = np.random.default_rng(SEED + 1)
 90    dim, T, trials = 24, 200, 300
 91    x = local.normal(size=(T, dim)).astype(np.float32)
 92    def orthogonal(scale):
 93        q, _ = np.linalg.qr(local.normal(size=(dim, dim)))
 94        return (scale*q).astype(np.float32)
 95    W_base = orthogonal(1.18)
 96    W_lo, W_hi = orthogonal(.68), orthogonal(1.28)
 97    U = (local.normal(size=(dim, dim))*0.18).astype(np.float32)
 98    p = .35
 99    base_slopes, idea_slopes = [], []
100    for _ in range(trials):
101        h = local.normal(size=dim).astype(np.float32); hp = h.copy(); hp[0] += 1
102        hi = local.normal(size=dim).astype(np.float32); hip = hi.copy(); hip[0] += 1
103        db, di = [], []
104        for t in range(T):
105            h, hp = np.tanh(W_base@h + U@x[t]), np.tanh(W_base@hp + U@x[t])
106            choose = local.random() < p
107            W = W_hi if choose else W_lo
108            hi, hip = np.tanh(W@hi + U@x[t]), np.tanh(W@hip + U@x[t])
109            db.append(np.linalg.norm(h-hp)+1e-12); di.append(np.linalg.norm(hi-hip)+1e-12)
110        tt = np.arange(T)
111        base_slopes.append(np.polyfit(tt[-80:], np.log(db[-80:]), 1)[0])
112        idea_slopes.append(np.polyfit(tt[-80:], np.log(di[-80:]), 1)[0])
113    expected_gain = (1-p)*.68+p*1.28
114    return {"baseline_predicted_upper_log_rate": math.log(1.18),
115            "baseline_observed_mean_log_distance_slope": float(np.mean(base_slopes)),
116            "idea_predicted_log_mean_gain": math.log(expected_gain),
117            "idea_observed_mean_log_distance_slope": float(np.mean(idea_slopes)),
118            "idea_mean_gain": expected_gain,
119            "idea_fraction_decaying": float(np.mean(np.asarray(idea_slopes)<0)),
120            "baseline_fraction_decaying": float(np.mean(np.asarray(base_slopes)<0))}
121
122
123def main():
124    product = expected_product_check()
125    schedule = varying_schedule_check()
126    boundary_rows, predicted_p, observed_p = mixed_boundary_sweep()
127    stability = recurrent_stability_comparison()
128    # Boundary tolerance is one grid step; product checks use a 2% MC tolerance.
129    product_ok = product["relative_error"] < .02 and schedule["max_relative_error_all_times"] < .02
130    boundary_ok = abs(predicted_p-observed_p) <= .06
131    # Nonlinear hidden states should show the predicted ordering and contraction signal.
132    stability_ok = (stability["idea_observed_mean_log_distance_slope"] < 0 and
133                    stability["baseline_observed_mean_log_distance_slope"] > stability["idea_observed_mean_log_distance_slope"])
134    result = {"seed": SEED, "expected_product_check": product,
135              "time_varying_schedule_check": schedule,
136              "boundary": {"predicted_p_hi": predicted_p, "observed_grid_p_hi": observed_p,
137                           "rows": boundary_rows},
138              "nonlinear_stability": stability,
139              "predictions_confirmed": {"expected_product": product_ok,
140                                        "mixed_gain_boundary": boundary_ok,
141                                        "nonlinear_attraction": stability_ok},
142              "worked": bool(product_ok and boundary_ok and stability_ok)}
143    print(json.dumps(result, indent=2))
144
145if __name__ == "__main__":
146    main()