1import json
 2import numpy as np
 3
 4SEED = 2060
 5rng = np.random.default_rng(SEED)
 6N = 16
 7# Random weighted undirected bipartite graph: its largest singular value is
 8# almost surely simple, so the negative flip eigenvalue is non-degenerate.
 9M = rng.normal(size=(8, 8))
10S = np.block([[np.zeros((8, 8)), M], [M.T, np.zeros((8, 8))]])
11S /= np.linalg.norm(M, 2)
12w, V = np.linalg.eigh(S)
13# J=-gamma*S has a unique dangerous flip multiplier -gamma, associated with S=+1.
14flip_i = np.argmax(w)
15vflip = V[:, flip_i]
16
17def C_at(a):
18    C = np.zeros((N, N)); C[a, a] = 1.0
19    return C
20
21def rho(A):
22    return float(np.max(np.abs(np.linalg.eigvals(A))))
23
24def flip_eigenvalue(A, reference):
25    ev = np.linalg.eigvals(A)
26    return float(ev[np.argmin(np.abs(ev - reference))].real)
27
28def trajectory(A, x0, steps=80):
29    xs = [x0.copy()]
30    for _ in range(steps): xs.append(A @ xs[-1])
31    return np.asarray(xs)
32
33# For this symmetric J, u=v and u_a v_a=v_a^2. The claimed first-order
34# formula predicts d lambda/d kappa = -v_a^2.
35gamma = 0.92
36J = -gamma * S
37anchor = int(np.argmax(np.abs(vflip)))
38random_anchor = int((anchor + 5) % N)
39uv = float(vflip[anchor] ** 2)
40small_ks = np.array([0.0, .001, .002, .005, .01])
41lams = [flip_eigenvalue(J - k*C_at(anchor), -gamma) for k in small_ks]
42observed_slope = float(np.polyfit(small_ks, lams, 1)[0])
43predicted_slope = -uv
44
45# For an unstable flip mode, the stated negative feedback has the wrong sign:
46# lambda(k) ~= -gamma - k*uv, so rho should increase, not contract.
47gamma_u = 1.08
48Ju = -gamma_u * S
49k_grid = np.linspace(0, 0.8, 161)
50flip_selected = [flip_eigenvalue(Ju-k*C_at(anchor), -gamma_u) for k in k_grid]
51rhos_selected = [rho(Ju-k*C_at(anchor)) for k in k_grid]
52rhos_random = [rho(Ju-k*C_at(random_anchor)) for k in k_grid]
53# Fit the unstable-mode slope at small gains and compare to prediction.
54unstable_slope = float(np.polyfit(k_grid[:6], flip_selected[:6], 1)[0])
55# Scan the whole range for a stable point (none is the expected result).
56stable_idx = [i for i, r in enumerate(rhos_selected) if r < 1.0]
57
58x0 = rng.normal(size=N); x0 /= np.linalg.norm(x0)
59def case(kind, kappa):
60    a = None if kind == 'none' else (random_anchor if kind == 'random' else anchor)
61    A = Ju if a is None else Ju - kappa*C_at(a)
62    x = trajectory(A, x0, 45)
63    return {'rho': rho(A), 'A_K': float(np.linalg.norm(x[-1]-x[-3])),
64            'norm_K': float(np.linalg.norm(x[-1]))}
65
66out = {
67 'setup': {'N': N, 'gamma_stable': gamma, 'gamma_unstable': gamma_u,
68           'anchor_selected': anchor, 'anchor_random': random_anchor,
69           'flip_coordinate_weight': uv},
70 'prediction_1_first_order_slope': {
71   'predicted_dlambda_dkappa': predicted_slope,
72   'observed_dlambda_dkappa': observed_slope,
73   'relative_error': abs(observed_slope-predicted_slope)/abs(predicted_slope),
74   'samples': [[float(k), float(l)] for k,l in zip(small_ks,lams)]},
75 'prediction_2_sign_for_unstable_flip': {
76   'predicted_dlambda_dkappa': -uv,
77   'observed_dlambda_dkappa': unstable_slope,
78   'rho_kappa_0': rhos_selected[0],
79   'rho_kappa_0.8': rhos_selected[-1],
80   'first_stable_kappa': None if not stable_idx else float(k_grid[stable_idx[0]])},
81 'prediction_3_anchor_selection': {
82   'kappa': 0.2,
83   'rho_random': rho(Ju-0.2*C_at(random_anchor)),
84   'rho_selected': rho(Ju-0.2*C_at(anchor)),
85   'selected_mode_weight': float(vflip[anchor]**2),
86   'random_mode_weight': float(vflip[random_anchor]**2)},
87 'mini_experiment': {
88   'none': case('none', 0),
89   'random': case('random', 0.2),
90   'selected': case('selected', 0.2)}}
91print(json.dumps(out, indent=2))