import json, math
import numpy as np

# Exact conformal two-state cocycle: A_i = r_i R(theta_i).
# Every product is scalar times orthogonal, hence lambda_+=lambda_- exactly.
def stationary(P):
    w, v = np.linalg.eig(P.T)
    q = np.real(v[:, np.argmin(abs(w - 1))])
    return q / q.sum()

def R(t):
    return np.array([[math.cos(t), -math.sin(t)], [math.sin(t), math.cos(t)]])

P = np.array([[.82, .18], [.27, .73]])
p = stationary(P)
r = np.array([.985, 1.015])
theta = np.array([.31, -.47])
A = np.array([r[i] * R(theta[i]) for i in range(2)])
lam = float(p @ np.log(r))
rows = []
for d in np.logspace(-8, -1, 15):
    # Radial perturbation of mode zero; matrix spectral distance is exactly d.
    new_r = np.array([r[0] + d, r[1]])
    exact = float(p @ np.log(new_r))
    delta = abs(exact - lam)
    rows.append({'d': float(d), 'observed_delta': delta,
                 'exact_delta': float(p[0] * math.log1p(d / r[0])),
                 'invlog': 1 / abs(math.log(d)),
                 'ratio_to_invlog': float(delta * abs(math.log(d)))})

# Predictions: exact zero gap; locally linear response; inverse-log is a loose envelope.
local = rows[:7]
slope = float(np.polyfit(np.log([x['d'] for x in local]),
                          np.log([x['observed_delta'] for x in local]), 1)[0])
C = float(max(x['ratio_to_invlog'] for x in rows))
coverage = float(np.mean([x['observed_delta'] <= C * x['invlog'] * (1 + 1e-12) for x in rows]))
Ccal = max(x['ratio_to_invlog'] for x in rows[:8])
heldout = rows[8:]
heldout_coverage = float(np.mean([x['observed_delta'] <= Ccal * x['invlog'] for x in heldout]))
out = {
    'stationary': p.tolist(), 'exact_lambda_plus_minus': lam,
    'max_exact_gap': 0.0, 'local_loglog_slope': slope,
    'expected_local_slope': 1.0, 'all_envelope_coverage': coverage,
    'calibration_C': Ccal, 'heldout_envelope_coverage': heldout_coverage,
    'rows': rows
}
with open('verification.json', 'w') as f:
    json.dump(out, f, indent=2)
print(json.dumps(out, indent=2))