Files
sim-numerics/stage3_bound_check.py

67 lines
2.6 KiB
Python

"""Stage 3 - the worst-case bound, checked against the runs that violate it.
Claim made for guard="bracket": one bisection is forced for every similarity
step that failed to halve the bracket, so at least one evaluation in two is a
bisection, so the bracket halves every two evaluations, so
nfev <= 2*ceil(log2(w/xtol)) + 4.
The Chandrupatla run reports 19 SOLVED violations of exactly that bound, so
either the reasoning or the implementation is wrong. This prints the actual
numbers for the worst offenders so the answer comes from the run, not from
re-reading the proof.
For orientation: Brent is NOT claimed to satisfy a 2x bound - his guarantee is
quadratic in the bisection count - so his violations are expected and are
printed only as a scale reference.
"""
from __future__ import annotations
import math
from scipy.optimize._tstutils import _CHANDRUPATLA_TESTS_DICTS
from simroot.baselines import bisection, brent
from simroot.pure import pure_root
XTOL = 1e-12
MAXFEV = 200
ACC = 1e-6
def main():
rows = []
for d in _CHANDRUPATLA_TESTS_DICTS:
a, b = float(d["bracket"][0]), float(d["bracket"][1])
f, root, pid = d["f"], float(d["root"]), d["ID"]
w = abs(b - a)
bound = 2 * math.ceil(math.log2(w / XTOL)) + 4
x_g, n_g, c_g = pure_root(f, a, b, XTOL, MAXFEV, similarity="ratio",
scale="derived", guard="bracket")
x_bi, n_bi, c_bi = bisection(f, a, b, XTOL, MAXFEV)
x_br, n_br, c_br = brent(f, a, b, XTOL, MAXFEV)
ok_g = c_g and abs(x_g - root) <= ACC
if ok_g and n_g > bound:
rows.append((pid, w, bound, n_g, n_bi if c_bi else -1,
n_br if c_br else -1))
rows.sort(key=lambda r: r[3] - r[2], reverse=True)
print(f"{'problem':10s}{'width':>12s}{'2x bound':>10s}"
f"{'bracket guard':>15s}{'bisection':>11s}{'brent':>8s}{'excess':>8s}")
print("-" * 74)
for pid, w, bound, n_g, n_bi, n_br in rows:
print(f"{pid:10s}{w:12.3g}{bound:10d}{n_g:15d}{n_bi:11d}"
f"{n_br:8d}{n_g - bound:8d}")
print(f"\n{len(rows)} solved runs exceed the bound.")
if rows:
worst = rows[0]
print(f"worst: {worst[0]}, {worst[3]} evaluations against a bound of "
f"{worst[2]}, while plain bisection took {worst[4]}.")
print("\nIf bisection itself needs MORE than ceil(log2(w/xtol)) "
"evaluations, the\nbound is mis-stated rather than violated: the "
"bisection count includes the\ntwo initial endpoint "
"evaluations, and the guarded method pays them too.")
if __name__ == "__main__":
main()