Symmetry-Resolved Fourier Bifurcation Monitor / symmetry_monitor.py
Failed on benchmark
1"""Symmetry-resolved Fourier bifurcation monitor and toy experiment."""
2import json
3import numpy as np
4
5
6def cyclic_shift(n):
7 S = np.zeros((n, n), dtype=complex)
8 for i in range(n):
9 S[(i + 1) % n, i] = 1.0
10 return S
11
12
13def spatial_fourier(n):
14 j = np.arange(n)
15 # columns are eigenvectors of the cyclic shift
16 return np.exp(2j * np.pi * np.outer(j, j) / n) / np.sqrt(n)
17
18
19def circulant_from_modes(modes):
20 n = len(modes)
21 U = spatial_fourier(n)
22 return U @ np.diag(modes) @ U.conj().T
23
24
25def sector_projectors(n):
26 U = spatial_fourier(n)
27 return [np.outer(U[:, k], U[:, k].conj()) for k in range(n)]
28
29
30def restricted_spectrum(M, n=None):
31 """Return each cyclic sector's scalar block and margin |1-lambda|."""
32 if n is None:
33 n = M.shape[0]
34 U = spatial_fourier(n)
35 blocks = np.diag(U.conj().T @ M @ U)
36 margins = np.abs(1.0 - blocks)
37 return blocks, margins
38
39
40def temporal_fft_symbols(period_map):
41 """Temporal Fourier coefficients of a sequence of one-step matrices."""
42 period_map = np.asarray(period_map)
43 return np.fft.fft(period_map, axis=0) / len(period_map)
44
45
46def winding_number(characteristic_values):
47 """Discrete phase winding around a closed contour (integer after rounding)."""
48 z = np.asarray(characteristic_values)
49 phase = np.unwrap(np.angle(z))
50 return int(np.rint((phase[-1] - phase[0]) / (2 * np.pi)))
51
52
53def circle_winding(eigenvalue, radius=0.15, points=512):
54 """Winding of chi(mu)=mu-lambda as mu traverses a circle centered at 1."""
55 theta = np.linspace(0, 2*np.pi, points + 1)
56 mu = 1 + radius * np.exp(1j * theta)
57 return winding_number(mu - eigenvalue)
58
59
60def run():
61 rng = np.random.default_rng(7)
62 n = 8
63 # A nontrivial equivariant ring: each spatial Fourier sector has a known gain.
64 # Sector 3 is deliberately the first one to hit the discrete-time unit circle.
65 base = np.array([0.62, 0.71, 0.80, 0.93, 0.77, 0.80, 0.71, 0.62])
66 boundary_rows = []
67 for gamma in np.linspace(0.70, 1.30, 25):
68 M = circulant_from_modes(gamma * base)
69 lam, margins = restricted_spectrum(M, n)
70 boundary_rows.append((gamma, float(np.max(np.abs(lam))), int(np.argmin(np.abs(1-lam)))))
71 # Predicted first unit-circle crossing is gamma=1/max(base)=1/base[3].
72 predicted_boundary = 1.0 / np.max(base)
73 observed_boundary = min(boundary_rows, key=lambda x: abs(x[1]-1))[0]
74
75 # Margin law: for the critical real sector, m=|1-gamma*a|, slope a.
76 gammas = np.linspace(0.55, 1.02, 48)
77 margins = []
78 for gamma in gammas:
79 lam, mm = restricted_spectrum(circulant_from_modes(gamma * base), n)
80 margins.append(mm[3])
81 margins = np.asarray(margins)
82 fit = np.polyfit(gammas[gammas < predicted_boundary], margins[gammas < predicted_boundary], 1)
83 predicted_slope = -base[3]
84 observed_slope = fit[0]
85
86 # Critical-sector prediction: add a small perturbation and measure the dominant FFT mode.
87 gamma = 1.03
88 M = circulant_from_modes(gamma * base)
89 x = rng.normal(size=n)
90 amplitudes = []
91 for _ in range(20):
92 x = M.real @ x
93 amplitudes.append(np.abs(np.fft.fft(x)))
94 dominant_sector = int(np.argmax(np.mean(amplitudes[-8:], axis=0)))
95
96 # Winding transition: circle centered at 1 encloses lambda iff its distance is < radius.
97 w_inside = circle_winding(1.0 + 0.05j, radius=0.15)
98 w_outside = circle_winding(1.0 + 0.30j, radius=0.15)
99
100 # Mini comparison: optimize a scalar gain per sector toward target gains.
101 # Baseline uses one global LR multiplier from the worst full-Jacobian margin;
102 # blockwise monitor damps only the near-critical sector.
103 target = np.array([0.45, 0.52, 0.58, 0.985, 0.50, 0.58, 0.52, 0.45])
104 init = np.array([0.55, 0.60, 0.66, 1.08, 0.63, 0.66, 0.60, 0.55])
105 def optimize(blockwise):
106 p = init.copy()
107 losses = []
108 for _ in range(80):
109 loss = float(np.mean((p-target)**2))
110 losses.append(loss)
111 grad = 2*(p-target)/n
112 # monitor return-map margins; near one gets damping, but no clipping of params
113 mm = np.abs(1-p)
114 if blockwise:
115 scale = np.minimum(1.0, mm/0.20)
116 # The critical sector is slowed, other sectors retain the base step.
117 step = 0.22 * scale
118 else:
119 scale = min(1.0, float(np.min(mm))/0.20)
120 step = np.full(n, 0.22 * scale)
121 p -= step * grad
122 return losses, p
123 loss_global, final_global = optimize(False)
124 loss_block, final_block = optimize(True)
125
126 result = {
127 "predictions": {
128 "unit_circle_boundary": {"predicted": predicted_boundary, "observed_grid": observed_boundary, "relative_error": abs(observed_boundary-predicted_boundary)/predicted_boundary},
129 "margin_slope": {"predicted": predicted_slope, "observed_fit": float(observed_slope), "relative_error": abs(observed_slope-predicted_slope)/abs(predicted_slope)},
130 "critical_spatial_sector": {"predicted": 3, "observed_dominant_fft": dominant_sector},
131 "winding": {"inside_circle": w_inside, "outside_circle": w_outside}
132 },
133 "comparison": {
134 "global_monitor_final_loss": loss_global[-1],
135 "blockwise_monitor_final_loss": loss_block[-1],
136 "global_monitor_steps_to_1e-4": next((i for i,x in enumerate(loss_global) if x < 1e-4), None),
137 "blockwise_monitor_steps_to_1e-4": next((i for i,x in enumerate(loss_block) if x < 1e-4), None),
138 "global_final": final_global.tolist(), "blockwise_final": final_block.tolist()
139 },
140 "checks": {"boundary_ok_5pct": bool(abs(observed_boundary-predicted_boundary)/predicted_boundary < .05),
141 "slope_ok_10pct": bool(abs(observed_slope-predicted_slope)/abs(predicted_slope) < .10),
142 "sector_ok": bool(dominant_sector == 3),
143 "winding_ok": bool(w_inside == 1 and w_outside == 0)}
144 }
145 print(json.dumps(result, indent=2))
146
147if __name__ == "__main__":
148 run()