import json, math, time
import numpy as np


def rot(theta):
    return np.array([[math.cos(theta), -math.sin(theta)],
                     [math.sin(theta), math.cos(theta)]], dtype=float)


def J_from_R(R):
    return (np.eye(R.shape[0]) + R) / 2.0


def fejer(J, y0, K):
    z = y0.copy()
    a = y0.copy()
    for _ in range(K):
        z = 2.0 * (J @ z) - z
        a += z
    return a / (K + 1)


def ppm(J, y0, n):
    y = y0.copy()
    for _ in range(n):
        y = J @ y
    return y


def residual(J, y):
    return np.linalg.norm(J @ y - y)


def run():
    y0 = np.array([1.0, 0.0])
    rows = []
    # Prediction 1: worst-case over rotation angle equals 1/(K+1).
    for K in [1, 3, 7, 15, 31]:
        th = np.linspace(0, 2*np.pi, 200001)
        # exact closed form for the residual of the averaged reflection
        vals = np.abs(np.sin((K + 1) * th / 2.0)) / (K + 1)
        observed = float(vals.max())
        predicted = 1.0 / (K + 1)
        rows.append({'test':'worst_case_bound','K':K,'observed':observed,
                     'predicted':predicted,'relative_error':abs(observed-predicted)/predicted})
    # Prediction 2: transition at (K+1) theta/2 = pi/2, i.e. theta*=pi/(K+1).
    # Sweep first crossing of 90% of worst case for K=15.
    K = 15
    angles = np.linspace(1e-7, np.pi, 200000)
    vals = np.abs(np.sin((K + 1) * angles / 2.0)) / (K + 1)
    target = .9 / (K + 1)
    idx = np.where(vals >= target)[0][0]
    observed_theta = float(angles[idx])
    predicted_theta = math.asin(.9) * 2.0 / (K + 1)
    rows.append({'test':'transition_90pct','K':K,'observed_theta':observed_theta,
                 'predicted_theta':predicted_theta,'observed_Ktheta':observed_theta*K,
                 'predicted_Ktheta':predicted_theta*K,
                 'relative_error':abs(observed_theta-predicted_theta)/predicted_theta})
    # Prediction 3: low-frequency residual is approximately theta/2 and does not
    # improve with K while K theta << 1; report ratios over a parameter sweep.
    low_rows = []
    for K in [3, 7, 15, 31]:
        theta = 0.01 / (K + 1)  # safely in the low-frequency regime
        R = rot(theta); J = J_from_R(R)
        got = residual(J, fejer(J, y0, K))
        predicted = theta / 2.0
        low_rows.append({'K':K,'theta':theta,'observed':got,'predicted_theta_over_2':predicted,
                         'ratio':got/predicted})
    # Direct numerical identity check on several angles and K.
    identity_err = 0.0
    for theta in np.linspace(0.03, 5.9, 12):
        R = rot(float(theta)); J = J_from_R(R)
        for K in [0,1,2,5,11]:
            yh = fejer(J, y0, K)
            rhs = (np.linalg.matrix_power(R, K+1) @ y0 - y0) / (2*(K+1))
            identity_err = max(identity_err, np.linalg.norm((J @ yh-yh)-rhs))
    # Equal oracle evaluations: Fejer K uses K J calls; PPM uses K calls.
    comparison = []
    theta = 0.35
    R = rot(theta); J = J_from_R(R)
    for K in [3,7,15,31]:
        yf = fejer(J,y0,K); yp = ppm(J,y0,K)
        comparison.append({'evaluations':K,'fejer_residual':residual(J,yf),
                           'ppm_residual':residual(J,yp),
                           'ratio_fejer_over_ppm':residual(J,yf)/residual(J,yp)})
    out = {'identity_max_abs_error':identity_err,'predictions':rows,
           'low_frequency':low_rows,'equal_eval_rotation_theta':theta,
           'comparison':comparison,
           'notes':'R is a rotation, hence nonexpansive; J=(I+R)/2 is an exact resolvent of a linear maximal monotone skew operator when -1 is absent.'}
    print(json.dumps(out, indent=2))

if __name__ == '__main__':
    run()