"""Stage 5 - does a three-point scale lift the order from 1.618 to 1.84?

Stage 4 measured the entire deficit against Brent as the ENDGAME, and the
endgame as nothing but the order of convergence: ours 1.584, Brent 1.809,
TOMS748 2.654.  The order follows from how the local scale of the value axis
is read: a two-point divided difference has slope error O(h), which gives the
secant order 1.618.  A three-point divided difference has error O(h^2), which
should give ~1.84 - Brent's order - with no new constant, because a divided
difference is a computation, not a fitted coefficient.

If that is right, `derived3` must (a) measure at order ~1.84, and (b) close
most of the APS gap.  If the order does not move, the explanation from
stage 4 was wrong and the paper must not carry it.

This is the SMOKE run: eight problems, before anything goes to the server.
"""
from __future__ import annotations

import math

from simroot.baselines import bisection, brent, toms748
from simroot.pure import pure_root

XTOL = 1e-12
MAXFEV = 200

PROBLEMS = [
    ("x^2 - 2", lambda x: x * x - 2.0, 0.0, 2.0, math.sqrt(2.0)),
    ("cos(x) - x", lambda x: math.cos(x) - x, 0.0, 1.0, 0.7390851332151607),
    ("exp(x) - 2", lambda x: math.exp(x) - 2.0, 0.0, 2.0, math.log(2.0)),
    ("x^3 - x - 2", lambda x: x ** 3 - x - 2.0, 1.0, 2.0, 1.5213797068045676),
    ("sin(x) - 0.5", lambda x: math.sin(x) - 0.5, 0.0, 1.5, math.asin(0.5)),
    ("x^10 - 1  STALL", lambda x: x ** 10 - 1.0, 0.0, 1.3, 1.0),
    ("exp(-x)(x-1) FLAT", lambda x: math.exp(-x) * (x - 1.0), -1.0, 40.0, 1.0),
    ("atan(x)", lambda x: math.atan(x), -1.0, 3.0, 0.0),
]


def order_of(method, f, a, b, root):
    """Observed order from the running-best error sequence."""
    xs = []

    def g(x):
        xs.append(x)
        return f(x)

    try:
        method(g, a, b, XTOL, MAXFEV)
    except Exception:
        return float("nan")
    best, errs = float("inf"), []
    for xi in xs:
        e = abs(xi - root)
        if e < best:
            best = e
            errs.append(e)
    errs = [e for e in errs if e > 1e-14]
    if len(errs) < 4:
        return float("nan")
    ps = []
    for i in range(1, len(errs) - 1):
        d1, d2 = errs[i] / errs[i - 1], errs[i + 1] / errs[i]
        if 0 < d1 < 1 and 0 < d2 < 1:
            ps.append(math.log(d2) / math.log(d1))
    return sorted(ps)[len(ps) // 2] if ps else float("nan")


def main():
    methods = [
        ("bisection", lambda f, a, b: bisection(f, a, b, XTOL, MAXFEV)),
        ("brent", lambda f, a, b: brent(f, a, b, XTOL, MAXFEV)),
        ("toms748", lambda f, a, b: toms748(f, a, b, XTOL, MAXFEV)),
        ("OURS 2-point .step", lambda f, a, b: pure_root(
            f, a, b, XTOL, MAXFEV, scale="derived", guard="step")),
        ("OURS 3-point .step", lambda f, a, b: pure_root(
            f, a, b, XTOL, MAXFEV, scale="derived3", guard="step")),
        ("OURS 3-point .both", lambda f, a, b: pure_root(
            f, a, b, XTOL, MAXFEV, scale="derived3", guard="both")),
    ]
    width = max(len(n) for n, _ in methods) + 1
    header = (f"{'method':{width}s}" +
              "".join(f"{p[0][:13]:>15s}" for p in PROBLEMS) +
              f"{'TOTAL':>8s}{'FAIL':>6s}")
    print(header)
    print("-" * len(header))
    for name, fn in methods:
        cells, total, fails = [], 0, 0
        for _, f, a, b, root in PROBLEMS:
            x, n, conv = fn(f, a, b)
            ok = conv and abs(x - root) <= 1e-8
            cells.append(f"{n:>15d}" if ok else f"{'FAIL(' + str(n) + ')':>15s}")
            total += n
            fails += 0 if ok else 1
        print(f"{name:{width}s}" + "".join(cells) + f"{total:>8d}{fails:>6d}")

    print("\nobserved order (median over the eight problems)")
    print(f"{'method':22s}{'order':>9s}   expected")
    print("-" * 56)
    exp = {"bisection": "1.0", "brent": "1.839", "toms748": "up to 2.7",
           "OURS 2-point .step": "1.618", "OURS 3-point .step": "~1.84",
           "OURS 3-point .both": "~1.84"}
    for name, fn in methods:
        ps = [order_of(lambda f_, a_, b_, xt, mx, _fn=fn: _fn(f_, a_, b_),
                       f, a, b, root)
              for _, f, a, b, root in PROBLEMS]
        ps = [p for p in ps if math.isfinite(p) and 0.5 < p < 5.0]
        med = sorted(ps)[len(ps) // 2] if ps else float("nan")
        print(f"{name:22s}{med:>9.3f}   {exp[name]}")


if __name__ == "__main__":
    main()