113 lines
4.3 KiB
Python
113 lines
4.3 KiB
Python
"""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()
|