1"""Small reproducible MVP for a truncated causal Volterra stabilizer."""
  2import json
  3from itertools import combinations_with_replacement
  4from pathlib import Path
  5import numpy as np
  6
  7SEED = 264
  8
  9
 10def lag_tuples(L, order):
 11    """Ordered simplex: 0 <= lag_n <= ... <= lag_1 < L."""
 12    return list(combinations_with_replacement(range(L), order))
 13
 14
 15class TruncatedVolterra:
 16    """Scalar finite-order causal compensator summing kernels of orders 2..N."""
 17    def __init__(self, L=1, order=2, weights=None):
 18        self.L, self.order = L, order
 19        self.weights = {}
 20        weights = {} if weights is None else weights
 21        for n in range(2, order + 1):
 22            self.weights[n] = np.asarray(weights.get(n, np.zeros(len(lag_tuples(L, n)))))
 23
 24    def __call__(self, history):
 25        total = 0.0
 26        for n in range(2, self.order + 1):
 27            tuples = lag_tuples(self.L, n)
 28            phi = np.array([np.prod([history[i] for i in t]) for t in tuples])
 29            total += float(phi @ self.weights[n])
 30        return total
 31
 32
 33def verify_math():
 34    q, c = 0.12, 0.08
 35    amplitudes = np.geomspace(1e-3, 0.4, 40)
 36    ratio = abs(c * amplitudes**3) / abs(q * amplitudes**2)
 37    slope = np.polyfit(np.log(amplitudes), np.log(ratio), 1)[0]
 38    return {
 39        "ordered_tuple_counts_L4": {str(n): len(lag_tuples(4, n)) for n in (2, 3)},
 40        "remainder_over_quadratic_loglog_slope": float(slope),
 41        "ratio_small_amplitude": float(ratio[0]),
 42        "ratio_large_amplitude": float(ratio[-1]),
 43        "stable_linear_transition": bool(abs(0.72) < 1),
 44    }
 45
 46
 47def fit_coefficients(seed=SEED):
 48    rng = np.random.default_rng(seed)
 49    A, q, c = 0.72, 0.55, 0.22
 50    x = rng.uniform(-0.35, 0.35, 500)
 51    y = A*x + q*x*x + c*x**3 + rng.normal(0, 0.001, len(x))
 52    residual = y - A*x
 53    q2 = float(np.linalg.lstsq((x*x)[:, None], residual, rcond=None)[0][0])
 54    q3, c3 = np.linalg.lstsq(np.column_stack([x*x, x**3]), residual, rcond=None)[0]
 55    return A, q, c, q2, float(q3), float(c3)
 56
 57
 58def rollout(A, q, c, compensator, x0, steps=1000):
 59    history = [float(x0)] * compensator.L
 60    trajectory = []
 61    for _ in range(steps):
 62        h = history[0]
 63        nxt = A*h + q*h*h + c*h**3 - compensator(history)
 64        history = [nxt] + history[:-1]
 65        trajectory.append(nxt)
 66        if not np.isfinite(nxt) or abs(nxt) > 1e12:
 67            return np.asarray(trajectory + [np.inf] * (steps-len(trajectory)))
 68    return np.asarray(trajectory)
 69
 70
 71def sweep(A, q, c, q2, q3, c3, amplitudes):
 72    controllers = {
 73        "baseline": TruncatedVolterra(L=1, order=2, weights={2: [0.0]}),
 74        "quadratic": TruncatedVolterra(L=1, order=2, weights={2: [q2]}),
 75        "cubic": TruncatedVolterra(L=1, order=3, weights={2: [q3], 3: [c3]}),
 76    }
 77    result = {}
 78    for name, ctrl in controllers.items():
 79        final, maxima, divergent = [], [], 0
 80        for x0 in amplitudes:
 81            z = rollout(A, q, c, ctrl, x0)
 82            bad = (not np.all(np.isfinite(z))) or np.max(np.abs(z)) > 1e6
 83            divergent += int(bad)
 84            final.append(np.inf if bad else abs(float(z[-1])))
 85            maxima.append(np.inf if bad else float(np.max(np.abs(z))))
 86        finite = np.isfinite(final)
 87        result[name] = {
 88            "mean_final_abs_state": float(np.mean(np.asarray(final)[finite])) if finite.any() else float("inf"),
 89            "max_abs_state": float(np.max(np.asarray(maxima)[finite])) if finite.any() else float("inf"),
 90            "divergent_rollouts": divergent,
 91            "bounded_initial_amplitude_max": float(amplitudes[np.where(finite)[0][-1]]) if finite.any() else 0.0,
 92            "parameter_count": 0 if name == "baseline" else (1 if name == "quadratic" else 2),
 93        }
 94    return result
 95
 96
 97def run():
 98    A, q, c, q2, q3, c3 = fit_coefficients()
 99    amplitudes = np.geomspace(0.02, 3.0, 60)
100    return {
101        "seed": SEED,
102        "plant": {"A": A, "quadratic": q, "cubic": c},
103        "learned_coefficients": {"qhat_quadratic": q2, "qhat_cubic": q3, "chat_cubic": c3},
104        "math_check": verify_math(),
105        "rollout_steps": 1000,
106        "amplitude_sweep": [float(x) for x in amplitudes],
107        "rollout_metrics": sweep(A, q, c, q2, q3, c3, amplitudes),
108    }
109
110
111if __name__ == "__main__":
112    result = run()
113    Path("results.json").write_text(json.dumps(result, indent=2))
114    print(json.dumps(result, indent=2))