"""Numerical MVP for the PIPO-PITO bounded recurrent gain idea.

The scalar positive plant is x' = (-d + w)x + p, y=x.  A fixed recurrent
input gain uses w=gamma.  The proposed controller uses w'=b-a*y*w.
This is intentionally a small mechanism test, not a neural benchmark.

Run:
  /home/maxwelhelp/main/bin/python3 pipo_pito_experiment.py
"""
import math
import json
import numpy as np


def controller_constant_y(a, b, V, w0, dt=1e-4):
    """Test w'=b-a*V*w against the PITO upper-bound prediction."""
    eps = 2.0*b/(a*V)
    predicted_T = max(0.0, math.log(2*w0/eps)/(a*V)) if w0 > eps else 0.0
    horizon = max(2*predicted_T, 4/(a*V), 0.1)
    n = int(math.ceil(horizon/dt))
    w = float(w0)
    crossing = 0.0 if w <= eps else None
    trace = []
    for i in range(1, n+1):
        w += dt*(b-a*V*w)
        if crossing is None and w <= eps:
            crossing = i*dt
        if i > 10 and w > eps/2 + 1e-8:
            trace.append((i*dt, math.log(w-eps/2)))
    slope = float('nan')
    if len(trace) > 3:
        xy = np.asarray(trace)
        slope = -float(np.polyfit(xy[:, 0], xy[:, 1], 1)[0])
    # Exact crossing is useful to distinguish the conservative theorem bound.
    exact = (0.0 if w0 <= eps else
             max(0., math.log((w0-eps/2)/(eps/2))/(a*V)))
    return {"epsilon": eps, "V": V, "V_times_epsilon": V*eps,
            "predicted_T_bound": predicted_T, "observed_T": crossing,
            "exact_T": exact, "rate_predicted": a*V,
            "rate_observed": slope, "rate_relative_error": abs(slope-a*V)/(a*V)}


def closed_loop(gamma, regulated, seed=0, dt=0.002, horizon=80.0):
    rng = np.random.default_rng(seed)
    d, p, a, b = 0.4, 0.02, 1.0, 0.2
    x = 0.1 + 0.01*rng.random()
    w = gamma if not regulated else 2.0
    max_x = x
    for i in range(int(horizon/dt)):
        if regulated:
            w = max(0.0, w + dt*(b-a*x*w))
        gain = w if regulated else gamma
        x = max(0.0, x + dt*((-d+gain)*x+p))
        max_x = max(max_x, x)
        if x > 1e100:
            break
    # For p>0, gamma<d converges, gamma=d grows linearly, gamma>d grows exponentially.
    unbounded = (not regulated and gamma >= d and x > 10.0)
    return {"final_x": float(x), "max_x": float(max_x), "final_w": float(w),
            "unbounded_by_80": bool(unbounded),
            "predicted_regime": ("bounded" if gamma < d else "unbounded (marginal at equality)")}


def main():
    a, b, w0 = 1.3, 0.7, 3.0
    epsilons = [0.2, 0.4, 0.8, 1.2]
    pito = [controller_constant_y(a, b, 2*b/(a*eps), w0) for eps in epsilons]
    # Fixed plant threshold: coefficient -d+gamma changes sign at gamma=d.
    gammas = [0.2, 0.39, 0.4, 0.6, 0.8]
    fixed = {str(g): closed_loop(g, False) for g in gammas}
    regulated = closed_loop(0.8, True)
    max_rate_error = max(r["rate_relative_error"] for r in pito)
    max_bound_slack = max((r["observed_T"] or 0)-r["predicted_T_bound"] for r in pito)
    result = {
      "parameters": {"a": a, "b": b, "d": 0.4, "p": 0.02,
                     "w0": w0, "dt": 0.002, "horizon": 80.0},
      "pito_sweep": pito,
      "fixed_gain_sweep": fixed,
      "regulated_at_nominal_gamma_0.8": regulated,
      "quantitative_checks": {
        "PITO_rate_prediction": {"predicted": "a*V", "max_relative_error": max_rate_error},
        "PITO_threshold_prediction": {"predicted": "V*epsilon=2*b/a",
                                       "observed_values": [r["V_times_epsilon"] for r in pito],
                                       "target": 2*b/a},
        "crossing_bound": {"predicted": "observed_T <= predicted_T_bound",
                            "maximum_observed_minus_predicted": max_bound_slack},
        "plant_boundary": {"predicted": "gamma=d=0.4; gamma<d bounded, gamma>=d unbounded",
                            "observed": {g: fixed[str(g)]["unbounded_by_80"] for g in gammas}},
        "regulated_boundedness": {"predicted": "feedback reduces gain and bounds x",
                                   "observed_max_x": regulated["max_x"],
                                   "observed_final_w": regulated["final_w"]}
      }
    }
    print(json.dumps(result, indent=2))


if __name__ == "__main__":
    main()