#!/usr/bin/env python3
"""Compare Lyapunov candidates for x'=-x+4y, y'=-2y.

Python 3, standard library only. State and time are normalized.
Run `python3 lyapunov.py` for the worked example, or add `--check`.
The sampled checks are regressions; the completed-square identity is the proof.
"""

import argparse
import math


INITIAL = (1.0, 1.0)
DT = 0.05
STEPS = 120


def validate_state(state):
    if not isinstance(state, (list, tuple)) or len(state) != 2:
        raise ValueError("Use two finite state coordinates.")
    if not all(isinstance(value, (int, float)) and math.isfinite(value) for value in state):
        raise ValueError("Use two finite state coordinates.")


def validate_candidate(candidate):
    if candidate not in ("euclidean", "weighted"):
        raise ValueError("Choose the euclidean or weighted candidate.")


def finite_result(value):
    if not math.isfinite(value):
        raise ValueError("The calculation exceeds the finite numeric range.")
    return 0.0 if value == 0 else value


def state_at(initial, time):
    """Exact forward solution, with no numerical integration step error."""
    validate_state(initial)
    if not isinstance(time, (int, float)) or not math.isfinite(time) or time < 0:
        raise ValueError("Use a finite, nonnegative time.")
    if time == 0:
        return tuple(initial)
    decay = math.exp(-time)
    coupling = 4 * decay * -math.expm1(-time)
    return (finite_result(initial[0] * decay + initial[1] * coupling),
            finite_result(initial[1] * decay * decay))


def value(state, candidate):
    validate_state(state)
    validate_candidate(candidate)
    x, y = state
    return finite_result(x * x + (4 if candidate == "weighted" else 1) * y * y)


def derivative(state, candidate):
    """Analytic derivative grad(V) dot f at a state."""
    validate_state(state)
    validate_candidate(candidate)
    x, y = state
    # Multiplication avoids Python's exponent overflow exception before validation.
    residual = x - 2 * y
    return finite_result(-2 * residual * residual - 8 * y * y if candidate == "weighted"
                         else -2 * x * x + 8 * x * y - 4 * y * y)


def trajectory(initial=INITIAL):
    validate_state(initial)
    rows = []
    for index in range(STEPS + 1):
        time = index * DT
        state = state_at(initial, time)
        rows.append({"time": time, "state": state,
                     "valueEuclidean": value(state, "euclidean"),
                     "valueWeighted": value(state, "weighted"),
                     "derivativeEuclidean": derivative(state, "euclidean"),
                     "derivativeWeighted": derivative(state, "weighted")})
    return rows


def check():
    """Check worked values, the ODE, and the algebra on representative states."""
    assert value(INITIAL, "euclidean") == 2
    assert derivative(INITIAL, "euclidean") == 2
    assert value(INITIAL, "weighted") == 5
    assert derivative(INITIAL, "weighted") == -10
    assert state_at(INITIAL, 0) == INITIAL
    assert state_at((0, 0), 6) == (0, 0)
    for x in (-2, -0.3, 0, 0.4, 2):
        for y in (-2, -0.3, 0, 0.4, 2):
            state = (x, y)
            norm_squared = x * x + y * y
            weighted = value(state, "weighted")
            assert norm_squared <= weighted + 1e-14 <= 4 * norm_squared + 2e-14
            direct = 2 * x * (-x + 4 * y) + 8 * y * (-2 * y)
            assert math.isclose(derivative(state, "weighted"), direct, rel_tol=1e-12, abs_tol=1e-12)
            assert derivative(state, "weighted") < 0 if norm_squared > 0 else derivative(state, "weighted") == 0
            for time in (0.1, 1, 3):
                now = state_at(state, time)
                h = 1e-5
                before, after = state_at(state, time - h), state_at(state, time + h)
                rates = [(after[i] - before[i]) / (2 * h) for i in range(2)]
                assert math.isclose(rates[0], -now[0] + 4 * now[1], rel_tol=1e-7, abs_tol=1e-8)
                assert math.isclose(rates[1], -2 * now[1], rel_tol=1e-7, abs_tol=1e-8)
            rows = trajectory(state)
            assert all(b["valueWeighted"] <= a["valueWeighted"] + 1e-12 for a, b in zip(rows, rows[1:]))
    assert trajectory()[10]["valueEuclidean"] > value(INITIAL, "euclidean")


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--check", action="store_true", help="run numerical regression checks")
    args = parser.parse_args()
    if args.check:
        check()
        print("Checks passed: worked values, exact ODE, candidate derivatives, and weighted decrease.")
        return
    print("Normalized system: x' = -x + 4y, y' = -2y; initial state = (1, 1)")
    print(" t        x        y       Ve      dVe       Vw      dVw")
    for time in (0, 0.5, 1, 2, 6):
        state = state_at(INITIAL, time)
        print(f"{time:3.1f} {state[0]:8.6f} {state[1]:8.6f} "
              f"{value(state, 'euclidean'):8.6f} {derivative(state, 'euclidean'):8.6f} "
              f"{value(state, 'weighted'):8.6f} {derivative(state, 'weighted'):8.6f}")
    print("Ve = x^2 + y^2 fails the decrease test at (1, 1): dVe = +2.")
    print("Vw = x^2 + 4y^2 succeeds: dVw = -2(x - 2y)^2 - 8y^2 < 0 away from zero.")
    print("A failed candidate does not prove instability. Both exact state coordinates approach zero.")


if __name__ == "__main__":
    main()