Digital-Recurrence Lyapunov Monitor / digital_recurrence.py

Mechanism failed

Raw ⬇ ZIP
  1"""Digital-recurrence monitor and a small logistic-map validation experiment."""
  2import json
  3import math
  4import numpy as np
  5
  6
  7class RecurrenceMonitor:
  8    def __init__(self, quantization=0.0, max_states=1000000):
  9        self.quantization = float(quantization)
 10        self.max_states = int(max_states)
 11
 12    def _key(self, z):
 13        a = np.asarray(z)
 14        if self.quantization > 0:
 15            # The infinity-norm tolerance is represented by a deterministic grid key.
 16            a = np.round(a / self.quantization).astype(np.int64)
 17        return a.tobytes()
 18
 19    def first_recurrence(self, trajectory):
 20        """Return (tau, period), where tau is the first repeated state's time."""
 21        seen = {}
 22        for t, z in enumerate(trajectory):
 23            if len(seen) >= self.max_states:
 24                break
 25            key = self._key(z)
 26            if key in seen:
 27                s = seen[key]
 28                return {"tau": int(t), "start": int(s), "period": int(t - s)}
 29            seen[key] = t
 30        return {"tau": None, "start": None, "period": None}
 31
 32
 33def logistic_trajectory(x0, n, dtype):
 34    x = np.array(x0, dtype=dtype)
 35    out = np.empty(n + 1, dtype=dtype)
 36    out[0] = x
 37    one = dtype(1.0)
 38    four = dtype(4.0)
 39    for i in range(n):
 40        x = dtype(four * x * (one - x))
 41        out[i + 1] = x
 42    return out
 43
 44
 45def finite_difference_slope(x0, horizon, dtype, delta=1e-7):
 46    """Slope of log forecast separation, matching the stated lambda estimator."""
 47    # Keep the perturbation representable in the selected precision.
 48    a = np.array(x0, dtype=dtype)
 49    b = np.array(x0 + delta, dtype=dtype)
 50    xa = logistic_trajectory(a, horizon, dtype)
 51    xb = logistic_trajectory(b, horizon, dtype)
 52    d = np.abs(xb.astype(np.float64) - xa.astype(np.float64))
 53    valid = d[1:] > 0
 54    if valid.sum() < 3:
 55        return float("nan")
 56    y = np.log(np.maximum(d[1:][valid], np.finfo(np.float64).tiny))
 57    t = np.arange(1, horizon + 1, dtype=np.float64)[valid]
 58    return float(np.polyfit(t, y, 1)[0])
 59
 60
 61def restart_lambda(x0, horizon, dtype, repeats=32):
 62    vals = []
 63    # Independent starts are deterministic, but separated by a fixed irrational-ish offset.
 64    for r in range(repeats):
 65        start = (float(x0) + (r + 1) * 0.00123456789) % 0.999
 66        v = finite_difference_slope(start, horizon, dtype)
 67        if np.isfinite(v):
 68            vals.append(v)
 69    return float(np.mean(vals)) if vals else float("nan")
 70
 71
 72def run():
 73    seed = 2386
 74    np.random.seed(seed)
 75    x0 = 0.1234567
 76    eps = 1e-3
 77    results = []
 78    # Long trajectories locate the recurrence scale; horizons test the predicted kink.
 79    for name, dtype, max_n in [("float16", np.float16, 20000),
 80                                ("float32", np.float32, 12000),
 81                                ("float64", np.float64, 20000)]:
 82        tr = logistic_trajectory(x0, max_n, dtype)
 83        rec = RecurrenceMonitor(max_states=max_n + 2).first_recurrence(tr)
 84        row = {"dtype": name, "recurrence": rec, "horizons": []}
 85        for n in [32, 128, 512, 2048, 4096, 8192, 12000]:
 86            if n > max_n:
 87                continue
 88            one = finite_difference_slope(x0, n, dtype)
 89            # Restart segments are short relative to the long-rollout test.
 90            k = min(128, n)
 91            rest = restart_lambda(x0, k, dtype)
 92            collapse = abs(one - rest) / (abs(rest) + eps) if np.isfinite(one) and np.isfinite(rest) else float("nan")
 93            row["horizons"].append({"N": n, "lambda_one": one,
 94                                    "lambda_restart": rest, "collapse": collapse,
 95                                    "alarm": bool(rec["tau"] is not None and n >= rec["tau"] and collapse > 1.0)})
 96        results.append(row)
 97    return {"seed": seed, "results": results}
 98
 99
100def json_safe(value):
101    if isinstance(value, dict):
102        return {k: json_safe(v) for k, v in value.items()}
103    if isinstance(value, list):
104        return [json_safe(v) for v in value]
105    if isinstance(value, float) and not math.isfinite(value):
106        return None
107    return value
108
109
110if __name__ == "__main__":
111    print(json.dumps(json_safe(run()), indent=2, allow_nan=False))