1"""Fading-memory habituation gate: math checks and a tiny controlled experiment."""
  2import json, math, random
  3from pathlib import Path
  4import numpy as np
  5
  6SEED = 1625
  7np.random.seed(SEED)
  8random.seed(SEED)
  9
 10
 11def simulate(rho, beta, stim, n=None, a0=0.0):
 12    """Update fading memory then apply nonlinear gain."""
 13    if n is not None:
 14        stim = [stim] * n
 15    a = float(a0)
 16    states, gains = [], []
 17    for s in stim:
 18        a = rho * a + (1.0 - rho) * float(s)
 19        states.append(a)
 20        gains.append(1.0 / (1.0 + beta * a))
 21    return np.asarray(states), np.asarray(gains)
 22
 23
 24def half_life(rho):
 25    return math.log(0.5) / math.log(rho)
 26
 27
 28def main():
 29    # Prediction 1: constant stimulation converges to a=s and gain 1/(1+beta*s).
 30    steady_rows = []
 31    for rho in [0.5, 0.8, 0.9, 0.97, 0.99]:
 32        # Enough steps for the slowest pole to settle accurately.
 33        n = 3000
 34        for beta in [0.5, 2.0, 5.0]:
 35            s = 1.7
 36            a, g = simulate(rho, beta, s, n=n)
 37            pred = 1.0 / (1.0 + beta * s)
 38            steady_rows.append({
 39                "rho": rho, "beta": beta,
 40                "a_error": abs(float(a[-1]) - s),
 41                "gain_observed": float(g[-1]),
 42                "gain_predicted": pred,
 43                "gain_abs_error": abs(float(g[-1]) - pred),
 44            })
 45
 46    # Prediction 2: after withdrawal a_k/a_0=rho^k; recovery half-life is ln(.5)/ln(rho).
 47    recovery_rows = []
 48    for rho in [0.5, 0.8, 0.9, 0.97, 0.99]:
 49        a0 = 2.3
 50        kmax = 3000
 51        states, _ = simulate(rho, 0.0, [0.0] * kmax, a0=a0)
 52        ratios = states / a0
 53        observed = (int(np.flatnonzero(ratios <= 0.5)[0]) + 1) if np.any(ratios <= 0.5) else None
 54        predicted = half_life(rho)
 55        predicted_index = int(math.ceil(predicted))
 56        # Fit log ratio versus k; exact prediction is slope log(rho).
 57        ks = np.arange(1, min(kmax, max(20, int(predicted * 8))))
 58        slope = float(np.polyfit(ks, np.log(ratios[ks - 1]), 1)[0])
 59        recovery_rows.append({
 60            "rho": rho, "half_life_predicted": predicted,
 61            "half_index_observed": observed,
 62            "half_index_predicted": predicted_index,
 63            "half_index_relative_error": abs(observed - predicted_index) / max(1, predicted_index),
 64            "log_slope_observed": slope,
 65            "log_slope_predicted": math.log(rho),
 66        })
 67
 68    # Prediction 3: rho<1 is stable/bounded for nonnegative bounded stimulation;
 69    # rho=1 freezes the state and rho>1 grows after withdrawal.
 70    stability_rows = []
 71    for rho in [0.5, 0.9, 0.99, 1.0, 1.01, 1.1]:
 72        states, _ = simulate(rho, 1.0, [0.0] * 100, a0=1.0)
 73        # Direct withdrawal exposes the pole: rho=1 does not decay and rho>1 grows.
 74        post = states
 75        stability_rows.append({
 76            "rho": rho,
 77            "constant_final_state": float(states[-1]),
 78            "withdrawal_final_state": float(post[-1]),
 79            "predicted_pole_behavior": "decay" if rho < 1 else ("constant" if rho == 1 else "growth"),
 80            "bounded_over_run": bool(np.all(np.isfinite(np.r_[states, post])) and np.max(np.abs(np.r_[states, post])) < 1e6),
 81        })
 82
 83    # Controlled mini-experiment: repeated familiar pulses versus novel pulses.
 84    # The baseline passes amplitude unchanged. The proposed gate suppresses repeated
 85    # stimulation while its state is reset between independent sequences.
 86    rho, beta, s = 0.9, 2.0, 1.0
 87    familiar = [s] * 20 + [0.0] * 20 + [s]
 88    novel = [0.0] * 20 + [s]
 89    a_f, g_f = simulate(rho, beta, familiar)
 90    a_n, g_n = simulate(rho, beta, novel)
 91    baseline_f = np.ones(len(familiar))
 92    baseline_n = np.ones(len(novel))
 93    # Matched linear-in-state gate: choose c so its steady gain equals the
 94    # nonlinear gate at s. This is not a true LTI input-output system because
 95    # it multiplies v by a state-dependent factor, but is a useful ablation.
 96    c = beta / (1.0 + beta * s)
 97    a_lin, _ = simulate(rho, 0.0, familiar)
 98    lin_g_f = 1.0 - c * a_lin
 99    a_lin_n, _ = simulate(rho, 0.0, novel)
100    lin_g_n = 1.0 - c * a_lin_n
101    comparison = {
102        "rho": rho, "beta": beta, "matched_linear_c": c,
103        "baseline_repeated_pulse_gain": float(baseline_f[19]),
104        "idea_repeated_pulse_gain": float(g_f[19]),
105        "linear_repeated_pulse_gain": float(lin_g_f[19]),
106        "baseline_post_withdrawal_gain": float(baseline_f[39]),
107        "idea_post_withdrawal_gain": float(g_f[39]),
108        "linear_post_withdrawal_gain": float(lin_g_f[39]),
109        "baseline_novel_gain": float(baseline_n[-1]),
110        "idea_novel_gain_after_reset": float(g_n[-1]),
111        "linear_novel_gain_after_reset": float(lin_g_n[-1]),
112        "repeated_suppression_ratio": float(g_f[19] / baseline_f[19]),
113        "novel_retention_ratio": float(g_n[-1] / baseline_n[-1]),
114    }
115
116    result = {
117        "seed": SEED,
118        "predictions": {
119            "steady_state": steady_rows,
120            "recovery": recovery_rows,
121            "stability": stability_rows,
122        },
123        "comparison": comparison,
124    }
125    out = Path("results.json")
126    out.write_text(json.dumps(result, indent=2))
127    print(json.dumps(result, indent=2))
128
129
130if __name__ == "__main__":
131    main()