Contractive Floquet return map / floquet_toy.py
Mechanism confirmed, baseline not beaten
1#!/usr/bin/env python3
2"""Toy verification of a contractive Floquet return map.
3
4State is (theta, r), where theta is phase and r is transverse deviation from
5an ideal periodic orbit. One return period applies
6 r' = q*r + delta
7 theta' = theta + shear*r' (mod 2*pi)
8The orbit distance is |r|; pointwise Euclidean distance on the unit circle is
9computed against the phase-aligned target point.
10"""
11import json
12import math
13from pathlib import Path
14import numpy as np
15
16SEED = 2094
17rng = np.random.default_rng(SEED)
18
19
20def step(z, q, defect=0.0, shear=0.0):
21 theta, r = z
22 rn = q * r + defect
23 return np.array([theta + shear * rn, rn], dtype=float)
24
25
26def simulate(q, r0, n, defect=0.0, shear=0.0):
27 z = np.array([0.0, r0], dtype=float)
28 zs = [z.copy()]
29 for _ in range(n):
30 z = step(z, q, defect, shear)
31 zs.append(z.copy())
32 return np.asarray(zs)
33
34
35def orbital_error(z):
36 return np.abs(z[:, 1])
37
38
39def pointwise_error(z):
40 # Target orbit Gamma(theta)=(cos(theta), sin(theta)); learned point uses
41 # the model phase and radius 1+r. This intentionally retains phase drift.
42 theta, r = z[:, 0], z[:, 1]
43 target = np.stack([np.cos(theta), np.sin(theta)], axis=1)
44 model = np.stack([(1.0 + r) * np.cos(theta), (1.0 + r) * np.sin(theta)], axis=1)
45 # Compare to phase-zero target: pointwise error exposes phase mismatch.
46 target0 = np.array([1.0, 0.0])
47 return np.linalg.norm(model - target0[None, :], axis=1)
48
49
50def predicted_r(q, r0, n, delta):
51 if abs(1.0 - q) < 1e-12:
52 return r0 + n * delta
53 return q**n * r0 + delta * (1.0 - q**n) / (1.0 - q)
54
55
56def empirical_q(q, n_pairs=1000):
57 # For the affine map, this should equal |q| exactly, independently of
58 # phase shear when transverse distance is measured in the return coordinate.
59 rs1 = rng.normal(0, 0.2, n_pairs)
60 rs2 = rng.normal(0, 0.2, n_pairs)
61 d0 = np.abs(rs1 - rs2)
62 d1 = np.abs((q * rs1) - (q * rs2))
63 return float(np.max(d1 / d0))
64
65
66def main():
67 report = {"seed": SEED, "predictions": {}, "comparison": {}}
68
69 # Prediction 1: transition at q=1; contraction/expansion rate is log(q).
70 qs = [0.80, 0.95, 1.00, 1.05, 1.20]
71 n = 30
72 r0 = 0.1
73 threshold_rows = []
74 for q in qs:
75 z = simulate(q, r0, n)
76 observed = abs(z[-1, 1] / r0)
77 expected = abs(q) ** n
78 slope = math.log(observed) / n if observed > 0 else float("-inf")
79 threshold_rows.append({"q": q, "observed_gain": observed,
80 "predicted_gain": expected,
81 "observed_log_rate": slope,
82 "predicted_log_rate": math.log(abs(q))})
83 report["predictions"]["threshold_q_equals_1"] = threshold_rows
84
85 # Prediction 2: constant defect gives floor delta/(1-q), linearly in delta.
86 q = 0.8
87 deltas = [0.001, 0.003, 0.01, 0.03]
88 n = 200
89 defect_rows = []
90 for delta in deltas:
91 z = simulate(q, 0.0, n, defect=delta)
92 observed = abs(z[-1, 1])
93 expected = delta / (1 - q)
94 defect_rows.append({"delta": delta, "observed_floor": observed,
95 "predicted_floor": expected,
96 "relative_error": abs(observed-expected)/expected})
97 report["predictions"]["defect_floor"] = defect_rows
98
99 # Prediction 3: phase shear makes pointwise error grow, while orbital error
100 # follows q^n and remains bounded. Compare shear=0 and a sheared orbit.
101 q, r0, shear, n = 0.85, 0.15, 2.5, 80
102 z = simulate(q, r0, n, shear=shear)
103 orb = orbital_error(z)
104 point = pointwise_error(z)
105 report["predictions"]["phase_invariant_vs_pointwise"] = {
106 "q": q, "shear": shear, "initial_transverse": r0,
107 "orbital_error_at_80": float(orb[-1]),
108 "orbital_predicted_at_80": float(abs(q**n * r0)),
109 "pointwise_error_at_80": float(point[-1]),
110 "pointwise_peak": float(np.max(point)),
111 "orbital_max": float(np.max(orb)),
112 "interpretation": "phase shear leaves transverse contraction intact but increases phase-sensitive error"
113 }
114
115 # Secondary mini-experiment: unconstrained return map versus contractive map.
116 # Small random per-cycle model defect makes this a deployment-like test.
117 cycles = 1000
118 trials = 200
119 sigma = 0.004
120 configs = {"baseline_unconstrained": 1.02, "contractive_idea": 0.80}
121 for name, q in configs.items():
122 orbital_end, point_end, qhat = [], [], []
123 for _ in range(trials):
124 z = np.array([0.0, 0.10])
125 for _ in range(cycles):
126 d = float(rng.normal(0.0, sigma))
127 z = step(z, q, defect=d, shear=1.5)
128 orbital_end.append(abs(z[1]))
129 point_end.append(float(pointwise_error(z[None, :])[0]))
130 qhat.append(empirical_q(q, 100))
131 report["comparison"][name] = {
132 "q": q, "cycles": cycles, "defect_sigma": sigma,
133 "median_orbital_error": float(np.median(orbital_end)),
134 "p95_orbital_error": float(np.quantile(orbital_end, .95)),
135 "median_pointwise_error": float(np.median(point_end)),
136 "fraction_orbital_inside_0.05": float(np.mean(np.asarray(orbital_end) < .05)),
137 "qhat_max_over_trials": float(np.max(qhat))
138 }
139
140 # Overall numerical checks used to decide whether the mechanism manifested.
141 checks = {
142 "threshold_relative_errors": [abs(x["observed_gain"]-x["predicted_gain"])/max(x["predicted_gain"],1e-12) for x in threshold_rows],
143 "max_defect_relative_error": max(x["relative_error"] for x in defect_rows),
144 "empirical_q_contractive": report["comparison"]["contractive_idea"]["qhat_max_over_trials"] < 0.8000001,
145 "empirical_q_baseline_noncontractive": report["comparison"]["baseline_unconstrained"]["qhat_max_over_trials"] > 1.0,
146 }
147 report["checks"] = checks
148 out = Path("results.json")
149 out.write_text(json.dumps(report, indent=2))
150 print(json.dumps(report, indent=2))
151
152
153if __name__ == "__main__":
154 main()