import json, math, os import numpy as np SEED = 2803 rng = np.random.default_rng(SEED) def linear_path(gains, choices): d = 1.0 out = [d] for a, i in zip(gains, choices): d *= a[i] out.append(d) return np.asarray(out) def mixed_boundary_sweep(): # Two scalar candidate maps: one contractive and one expansive. lo, hi = 0.60, 1.40 ps = np.linspace(0.0, 1.0, 21) n = 20000 horizon = 30 rows = [] for p in ps: # candidate 1 (hi) is selected with probability p choices = rng.random((n, horizon)) < p # Vectorized path distances; this is also a direct Monte Carlo estimate. logs = np.where(choices, math.log(hi), math.log(lo)) logd = np.cumsum(logs, axis=1) # log of E[d_t], stabilized, equals t*log((1-p)lo+p*hi) theoretically. t = np.arange(1, horizon + 1) log_mean_d = np.array([np.logaddexp.reduce(logd[:, j]) - math.log(n) for j in range(horizon)]) slope = np.polyfit(t[-15:], log_mean_d[-15:], 1)[0] predicted = math.log((1-p)*lo + p*hi) rows.append({"p_hi": float(p), "predicted_log_exponent": predicted, "observed_log_exponent": float(slope), "predicted_mean_gain": (1-p)*lo+p*hi}) # The expected-product boundary is mean gain = 1. predicted_p = (1-lo)/(hi-lo) observed_p = min(rows, key=lambda r: abs(r["observed_log_exponent"]))["p_hi"] return rows, predicted_p, observed_p def expected_product_check(): # Time-inhomogeneous independent routing: E product is product of expected gains. probs = np.array([.15, .75, .35, .60, .20, .80, .45, .55]) gains = np.array([[.55, 1.35], [.70, 1.20], [.50, 1.50], [.80, 1.10], [.60, 1.30], [.75, 1.25], [.55, 1.40], [.65, 1.15]]) n = 300000 choices = rng.random((n, len(probs))) < probs[None, :] paths = np.where(choices, gains[:, 1], gains[:, 0]) empirical = np.mean(np.prod(paths, axis=1)) predicted = float(np.prod((1-probs)*gains[:, 0] + probs*gains[:, 1])) relerr = abs(empirical-predicted)/predicted # Also verify the finite-horizon geometric rate for a stationary contractive mix. m = .94 horizon = 12 stationary_choices = rng.random((n, horizon)) < .4 stationary_paths = np.where(stationary_choices, 1.3, .7) empirical_h = np.mean(np.prod(stationary_paths, axis=1)) predicted_h = .94**horizon return {"predicted_product": predicted, "empirical_product": float(empirical), "relative_error": float(relerr), "stationary_horizon": horizon, "predicted_geometric_product": predicted_h, "empirical_geometric_product": float(empirical_h)} def varying_schedule_check(): # A deterministic time-varying schedule tests the product of p_t-weighted gains. p = np.array([.1,.7,.2,.8,.3,.6,.4,.5,.25,.75]) lo = np.array([.55,.65,.60,.70,.58,.62,.56,.68,.61,.59]) hi = np.array([1.20,1.10,1.30,1.05,1.25,1.15,1.18,1.08,1.22,1.12]) means = (1-p)*lo+p*hi predicted = np.cumprod(means) n = 250000 choices = rng.random((n, len(p))) < p products = np.cumprod(np.where(choices, hi, lo), axis=1) observed = np.mean(products, axis=0) relerr = np.max(np.abs(observed-predicted)/(predicted+1e-12)) return {"predicted_final_product": float(predicted[-1]), "observed_final_product": float(observed[-1]), "max_relative_error_all_times": float(relerr), "mean_log_gain": float(np.mean(np.log(means)))} def recurrent_stability_comparison(): # Same nonlinear tanh setup: unconstrained baseline versus random candidate maps. # Distances are measured on paired trajectories with identical inputs and routing. local = np.random.default_rng(SEED + 1) dim, T, trials = 24, 200, 300 x = local.normal(size=(T, dim)).astype(np.float32) def orthogonal(scale): q, _ = np.linalg.qr(local.normal(size=(dim, dim))) return (scale*q).astype(np.float32) W_base = orthogonal(1.18) W_lo, W_hi = orthogonal(.68), orthogonal(1.28) U = (local.normal(size=(dim, dim))*0.18).astype(np.float32) p = .35 base_slopes, idea_slopes = [], [] for _ in range(trials): h = local.normal(size=dim).astype(np.float32); hp = h.copy(); hp[0] += 1 hi = local.normal(size=dim).astype(np.float32); hip = hi.copy(); hip[0] += 1 db, di = [], [] for t in range(T): h, hp = np.tanh(W_base@h + U@x[t]), np.tanh(W_base@hp + U@x[t]) choose = local.random() < p W = W_hi if choose else W_lo hi, hip = np.tanh(W@hi + U@x[t]), np.tanh(W@hip + U@x[t]) db.append(np.linalg.norm(h-hp)+1e-12); di.append(np.linalg.norm(hi-hip)+1e-12) tt = np.arange(T) base_slopes.append(np.polyfit(tt[-80:], np.log(db[-80:]), 1)[0]) idea_slopes.append(np.polyfit(tt[-80:], np.log(di[-80:]), 1)[0]) expected_gain = (1-p)*.68+p*1.28 return {"baseline_predicted_upper_log_rate": math.log(1.18), "baseline_observed_mean_log_distance_slope": float(np.mean(base_slopes)), "idea_predicted_log_mean_gain": math.log(expected_gain), "idea_observed_mean_log_distance_slope": float(np.mean(idea_slopes)), "idea_mean_gain": expected_gain, "idea_fraction_decaying": float(np.mean(np.asarray(idea_slopes)<0)), "baseline_fraction_decaying": float(np.mean(np.asarray(base_slopes)<0))} def main(): product = expected_product_check() schedule = varying_schedule_check() boundary_rows, predicted_p, observed_p = mixed_boundary_sweep() stability = recurrent_stability_comparison() # Boundary tolerance is one grid step; product checks use a 2% MC tolerance. product_ok = product["relative_error"] < .02 and schedule["max_relative_error_all_times"] < .02 boundary_ok = abs(predicted_p-observed_p) <= .06 # Nonlinear hidden states should show the predicted ordering and contraction signal. stability_ok = (stability["idea_observed_mean_log_distance_slope"] < 0 and stability["baseline_observed_mean_log_distance_slope"] > stability["idea_observed_mean_log_distance_slope"]) result = {"seed": SEED, "expected_product_check": product, "time_varying_schedule_check": schedule, "boundary": {"predicted_p_hi": predicted_p, "observed_grid_p_hi": observed_p, "rows": boundary_rows}, "nonlinear_stability": stability, "predictions_confirmed": {"expected_product": product_ok, "mixed_gain_boundary": boundary_ok, "nonlinear_attraction": stability_ok}, "worked": bool(product_ok and boundary_ok and stability_ok)} print(json.dumps(result, indent=2)) if __name__ == "__main__": main()