1"""Representation-Invariant Authority Demand: small numerical verification.
  2
  3Toy system: xdot = alpha*x + u, x in R^2, ||u||_inf <= rho.
  4Safe set: h(x)=1-||x||^2 >= 0.  On the boundary, the exact IAD is alpha,
  5because r=2 alpha and a=2 rho (|x_1|+|x_2|), whose worst ratio is alpha/rho.
  6The exact boundary feasibility threshold is therefore rho=alpha.
  7"""
  8import json
  9import numpy as np
 10
 11SEED = 3140
 12rng = np.random.default_rng(SEED)
 13alpha = 0.73
 14N = 200_000
 15# Uniform angular boundary sample, plus axes so the supremum is represented.
 16theta = rng.uniform(0, 2*np.pi, N)
 17theta = np.concatenate([theta, [0, np.pi/2, np.pi, 3*np.pi/2]])
 18x = np.column_stack([np.cos(theta), np.sin(theta)])
 19
 20def base_h(z):
 21    return 1.0 - np.sum(z*z, axis=-1)
 22
 23def grad_base_h(z):
 24    return -2.0*z
 25
 26def transformed(z, name):
 27    h = base_h(z)
 28    if name == 'h':
 29        q, qp = h, np.ones_like(h)
 30    elif name == '2h':
 31        q, qp = 2*h, 2*np.ones_like(h)
 32    elif name == 'h+h^3':
 33        q, qp = h + h**3, 1 + 3*h*h
 34    elif name == 'log(1+h)':
 35        q, qp = np.log1p(h), 1/(1+h)
 36    elif name == 'h^3 (invalid)':
 37        q, qp = h**3, 3*h*h
 38    else:
 39        raise ValueError(name)
 40    # On/near the boundary grad(q(h)) = q'(h) grad(h).
 41    grad = qp[:, None] * grad_base_h(z)
 42    return q, grad
 43
 44def iad_for(rho, name):
 45    _, grad = transformed(x, name)
 46    f = alpha*x
 47    r = -np.sum(grad*f, axis=1)
 48    # support of [-rho,rho]^2: rho * L1 norm of grad
 49    a = rho*np.sum(np.abs(grad), axis=1)
 50    d = np.zeros_like(r)
 51    valid = a > 1e-14
 52    d[valid] = np.maximum(0, r[valid]) / a[valid]
 53    d[~valid & (r > 0)] = np.inf
 54    return float(np.max(d)), float(np.min(-r+a)), float(np.mean(r))
 55
 56def softsup(values, tau=0.02):
 57    m = np.max(values)
 58    return float(m + tau*np.log(np.mean(np.exp((values-m)/tau))))
 59
 60def run():
 61    names = ['h', '2h', 'h+h^3', 'log(1+h)', 'h^3 (invalid)']
 62    rows = {}
 63    for name in names:
 64        D, _, mean_r = iad_for(1.0, name)
 65        # raw boundary outward penalty changes with barrier representation;
 66        # IAD is a ratio and should not (for valid q'(0)>0).
 67        rows[name] = {'D_at_rho_1': D, 'raw_mean_r': mean_r}
 68
 69    # Transition: exact sampled min controlled derivative at each actuator scale.
 70    scales = np.array([0.50, 0.65, 0.72, 0.73, 0.74, 0.85, 1.00])
 71    transition = []
 72    for rho in scales:
 73        D, min_deriv, _ = iad_for(rho, 'h')
 74        transition.append({'rho': float(rho), 'D': D,
 75                           'min_boundary_derivative': min_deriv,
 76                           'feasible_sampled': bool(min_deriv >= -1e-10)})
 77
 78    # Empirical soft supremum converges to the exact supremum as sampling grows.
 79    # Compute pointwise demand for the base representation for a stable diagnostic.
 80    _, grad = transformed(x, 'h')
 81    r = -np.sum(grad*(alpha*x), axis=1)
 82    a = np.sum(np.abs(grad), axis=1)  # rho=1
 83    d = np.maximum(0, r)/a
 84    soft = softsup(d)
 85    # Finite-difference check of the analytic gradient on several boundary points.
 86    eps = 1e-6
 87    check = x[:32]
 88    numeric = np.empty_like(check)
 89    for j in range(2):
 90        step = np.zeros_like(check)
 91        step[:, j] = eps
 92        numeric[:, j] = (base_h(check + step) - base_h(check - step)) / (2 * eps)
 93    analytic = grad_base_h(check)
 94    grad_err = float(np.max(np.abs(numeric - analytic)))
 95
 96    result = {
 97        'seed': SEED, 'alpha_exact_threshold': alpha,
 98        'representation_checks': rows,
 99        'transition': transition,
100        'softsup_tau_0.02_rho_1': soft,
101        'max_boundary_gradient_finite_difference_error': grad_err,
102        'valid_invariance_relative_spread': float((max(rows[n]['D_at_rho_1'] for n in names[:4]) - min(rows[n]['D_at_rho_1'] for n in names[:4]))/alpha),
103        'invalid_h3_D': rows['h^3 (invalid)']['D_at_rho_1'],
104        'notes': 'For h^3, q_prime(0)=0, so the defining-function regularity assumption fails and the boundary gradient vanishes.'
105    }
106    with open('results.json', 'w') as f:
107        json.dump(result, f, indent=2)
108    print(json.dumps(result, indent=2))
109
110if __name__ == '__main__':
111    run()