import json
import numpy as np

SEED = 2060
rng = np.random.default_rng(SEED)
N = 16
# Random weighted undirected bipartite graph: its largest singular value is
# almost surely simple, so the negative flip eigenvalue is non-degenerate.
M = rng.normal(size=(8, 8))
S = np.block([[np.zeros((8, 8)), M], [M.T, np.zeros((8, 8))]])
S /= np.linalg.norm(M, 2)
w, V = np.linalg.eigh(S)
# J=-gamma*S has a unique dangerous flip multiplier -gamma, associated with S=+1.
flip_i = np.argmax(w)
vflip = V[:, flip_i]

def C_at(a):
    C = np.zeros((N, N)); C[a, a] = 1.0
    return C

def rho(A):
    return float(np.max(np.abs(np.linalg.eigvals(A))))

def flip_eigenvalue(A, reference):
    ev = np.linalg.eigvals(A)
    return float(ev[np.argmin(np.abs(ev - reference))].real)

def trajectory(A, x0, steps=80):
    xs = [x0.copy()]
    for _ in range(steps): xs.append(A @ xs[-1])
    return np.asarray(xs)

# For this symmetric J, u=v and u_a v_a=v_a^2. The claimed first-order
# formula predicts d lambda/d kappa = -v_a^2.
gamma = 0.92
J = -gamma * S
anchor = int(np.argmax(np.abs(vflip)))
random_anchor = int((anchor + 5) % N)
uv = float(vflip[anchor] ** 2)
small_ks = np.array([0.0, .001, .002, .005, .01])
lams = [flip_eigenvalue(J - k*C_at(anchor), -gamma) for k in small_ks]
observed_slope = float(np.polyfit(small_ks, lams, 1)[0])
predicted_slope = -uv

# For an unstable flip mode, the stated negative feedback has the wrong sign:
# lambda(k) ~= -gamma - k*uv, so rho should increase, not contract.
gamma_u = 1.08
Ju = -gamma_u * S
k_grid = np.linspace(0, 0.8, 161)
flip_selected = [flip_eigenvalue(Ju-k*C_at(anchor), -gamma_u) for k in k_grid]
rhos_selected = [rho(Ju-k*C_at(anchor)) for k in k_grid]
rhos_random = [rho(Ju-k*C_at(random_anchor)) for k in k_grid]
# Fit the unstable-mode slope at small gains and compare to prediction.
unstable_slope = float(np.polyfit(k_grid[:6], flip_selected[:6], 1)[0])
# Scan the whole range for a stable point (none is the expected result).
stable_idx = [i for i, r in enumerate(rhos_selected) if r < 1.0]

x0 = rng.normal(size=N); x0 /= np.linalg.norm(x0)
def case(kind, kappa):
    a = None if kind == 'none' else (random_anchor if kind == 'random' else anchor)
    A = Ju if a is None else Ju - kappa*C_at(a)
    x = trajectory(A, x0, 45)
    return {'rho': rho(A), 'A_K': float(np.linalg.norm(x[-1]-x[-3])),
            'norm_K': float(np.linalg.norm(x[-1]))}

out = {
 'setup': {'N': N, 'gamma_stable': gamma, 'gamma_unstable': gamma_u,
           'anchor_selected': anchor, 'anchor_random': random_anchor,
           'flip_coordinate_weight': uv},
 'prediction_1_first_order_slope': {
   'predicted_dlambda_dkappa': predicted_slope,
   'observed_dlambda_dkappa': observed_slope,
   'relative_error': abs(observed_slope-predicted_slope)/abs(predicted_slope),
   'samples': [[float(k), float(l)] for k,l in zip(small_ks,lams)]},
 'prediction_2_sign_for_unstable_flip': {
   'predicted_dlambda_dkappa': -uv,
   'observed_dlambda_dkappa': unstable_slope,
   'rho_kappa_0': rhos_selected[0],
   'rho_kappa_0.8': rhos_selected[-1],
   'first_stable_kappa': None if not stable_idx else float(k_grid[stable_idx[0]])},
 'prediction_3_anchor_selection': {
   'kappa': 0.2,
   'rho_random': rho(Ju-0.2*C_at(random_anchor)),
   'rho_selected': rho(Ju-0.2*C_at(anchor)),
   'selected_mode_weight': float(vflip[anchor]**2),
   'random_mode_weight': float(vflip[random_anchor]**2)},
 'mini_experiment': {
   'none': case('none', 0),
   'random': case('random', 0.2),
   'selected': case('selected', 0.2)}}
print(json.dumps(out, indent=2))