242 lines
9.4 KiB
Python
242 lines
9.4 KiB
Python
"""Stage 19 - the objection that the repairs compared against are from 1973.
|
|
|
|
Claim A as it stood was: the credit projection repairs regula falsi better
|
|
than Illinois (1971), Pegasus (1972) and Anderson-Bjorck (1973). A reader
|
|
answers that the field has moved on, and is right to: TOMS 748 (1995),
|
|
Chandrupatla (1997) and ITP (2021) all beat credit-repaired regula falsi.
|
|
Repairing a three-century-old interpolation better than the 1970s did is not
|
|
a modern result.
|
|
|
|
The credit projection is not tied to regula falsi, though. It is a safeguard
|
|
for ANY step rule. So the strong form of the claim is not "our repair of
|
|
their old rule beats their old repair" but
|
|
|
|
take the best interpolation in the field and replace its own acceptance
|
|
test with the credit projection.
|
|
|
|
Chandrupatla's method is exactly that pairing: an inverse quadratic
|
|
interpolation plus an acceptance test derived from where that interpolation
|
|
stays inside the bracket. Swapping the test for the credit projection, and
|
|
changing nothing else, is the sharpest test available of whether the
|
|
safeguard is worth anything against the modern state of the art.
|
|
|
|
Two further modern repairs of the false-position family are included so that
|
|
the 1970s are not the only comparison: Ford's Illinois-type variant (1995)
|
|
and Galdino's (2011).
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
|
|
from scipy.optimize._tstutils import (_APS_TESTS_DICTS,
|
|
_CHANDRUPATLA_TESTS_DICTS)
|
|
|
|
from simroot.baselines import Counter, _Budget, chandrupatla
|
|
from simroot.pure import _Safeguard
|
|
from stage18_bulletproof import falsi_driver, load_roots
|
|
|
|
XTOL, MAXFEV, ACC = 1e-12, 200, 1e-6
|
|
|
|
|
|
def chandrupatla_guard(f, a, b, xtol, maxfev, guard="own"):
|
|
"""Chandrupatla's method with its safeguard pluggable.
|
|
|
|
guard="own" the published acceptance test on xi and phi
|
|
guard="credit" always interpolate, and project the result onto the
|
|
interval the earned bisection budget allows
|
|
guard="none" always interpolate, clamp to the bracket only
|
|
|
|
Everything else - the interpolation, the bracket bookkeeping, the
|
|
termination test - is identical, so a difference can only be the guard.
|
|
"""
|
|
g = Counter(f, maxfev)
|
|
x1, x2 = b, a
|
|
f1 = f2 = 0.0
|
|
try:
|
|
f1, f2 = g(x1), g(x2)
|
|
if f1 == 0.0:
|
|
return x1, g.n, True
|
|
if f2 == 0.0:
|
|
return x2, g.n, True
|
|
x3, f3 = x2, f2
|
|
t = 0.5
|
|
sg = _Safeguard("credit", abs(b - a), xtol)
|
|
while True:
|
|
x = x1 + t * (x2 - x1)
|
|
w_before = abs(x2 - x1)
|
|
fx = g(x)
|
|
if (fx > 0.0) == (f1 > 0.0):
|
|
x3, f3 = x1, f1
|
|
x1, f1 = x, fx
|
|
else:
|
|
x3, f3 = x2, f2
|
|
x2, f2 = x1, f1
|
|
x1, f1 = x, fx
|
|
if guard == "credit":
|
|
sg.note_shrink(w_before, abs(x2 - x1))
|
|
xm, fm = (x1, f1) if abs(f1) < abs(f2) else (x2, f2)
|
|
if fm == 0.0:
|
|
return xm, g.n, True
|
|
tol = 2.0 * 2.22e-16 * abs(xm) + xtol
|
|
span = abs(x2 - x1)
|
|
if span == 0.0:
|
|
return xm, g.n, True
|
|
tl = tol / span
|
|
if tl > 0.5:
|
|
return xm, g.n, True
|
|
|
|
t = 0.5
|
|
if x3 != x2 and f3 != f2:
|
|
xi = (x1 - x2) / (x3 - x2)
|
|
ph = (f1 - f2) / (f3 - f2)
|
|
d1, d2, d3 = f1 - f2, f1 - f3, f2 - f3
|
|
usable = d1 != 0.0 and d2 != 0.0 and d3 != 0.0
|
|
if usable:
|
|
t_iqi = ((f1 / d1) * (f3 / (-d3))
|
|
+ ((x3 - x1) / (x2 - x1)) * (f1 / d2) * (f2 / d3))
|
|
if guard == "own":
|
|
if (0.0 < xi < 1.0
|
|
and 1.0 - math.sqrt(1.0 - xi) < ph
|
|
< math.sqrt(xi)):
|
|
t = t_iqi
|
|
else:
|
|
t = t_iqi
|
|
if guard == "credit" and math.isfinite(t):
|
|
lo = min(x1, x2)
|
|
hi = max(x1, x2)
|
|
cand = x1 + t * (x2 - x1)
|
|
proj = sg.project_credit(cand, lo, hi)
|
|
if x2 != x1:
|
|
t = (proj - x1) / (x2 - x1)
|
|
if not math.isfinite(t):
|
|
t = 0.5
|
|
t = min(max(t, tl), 1.0 - tl)
|
|
except _Budget:
|
|
return (x1 if abs(f1) < abs(f2) else x2), g.n, False
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# two modern repairs of the false-position family, so the 1970s are not the
|
|
# only comparison
|
|
# --------------------------------------------------------------------------
|
|
def falsi_modern(f, a, b, xtol, maxfev, variant):
|
|
"""variant="ford" Ford 1995, Illinois-type with the Anderson-Bjorck
|
|
factor bounded away from zero
|
|
variant="galdino" Galdino 2011, the factor f_moved/(f_moved + fx)
|
|
applied to both ordinates
|
|
"""
|
|
g = Counter(f, maxfev)
|
|
x = 0.5 * (a + b)
|
|
try:
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
if fa * fb > 0.0:
|
|
raise ValueError("no sign change")
|
|
while abs(b - a) > xtol:
|
|
lo, hi = (a, b) if a < b else (b, a)
|
|
denom = fb - fa
|
|
x_new = 0.5 * (a + b) if denom == 0.0 else (a * fb - b * fa) / denom
|
|
if not math.isfinite(x_new) or not (lo < x_new < hi):
|
|
x_new = 0.5 * (a + b)
|
|
if x_new in (a, b):
|
|
x_new = 0.5 * (a + b)
|
|
if x_new in (a, b):
|
|
return x_new, g.n, abs(b - a) <= xtol
|
|
fx = g(x_new)
|
|
if fx == 0.0:
|
|
return x_new, g.n, True
|
|
|
|
def factor(f_stale, f_moved):
|
|
if variant == "ford":
|
|
m = 1.0 - fx / f_moved if f_moved != 0.0 else 0.5
|
|
return m if m > 1e-3 else 0.5
|
|
m = f_moved / (f_moved + fx) if (f_moved + fx) != 0.0 else 0.5
|
|
return m if 0.0 < m < 1.0 else 0.5
|
|
|
|
if fa * fx < 0.0:
|
|
fa *= factor(fa, fb)
|
|
b, fb = x_new, fx
|
|
else:
|
|
fb *= factor(fb, fa)
|
|
a, fa = x_new, fx
|
|
x = x_new
|
|
return 0.5 * (a + b), g.n, True
|
|
except _Budget:
|
|
return x, g.n, False
|
|
|
|
|
|
def score(method, probs):
|
|
total = solved = worst = 0
|
|
for f, a, b, root in probs:
|
|
try:
|
|
x, n, conv = method(f, a, b, XTOL, MAXFEV)
|
|
except Exception:
|
|
x, n, conv = float("nan"), MAXFEV, False
|
|
ok = bool(conv and math.isfinite(x) and abs(x - root) <= ACC)
|
|
total += n if ok else MAXFEV
|
|
solved += ok
|
|
if ok:
|
|
worst = max(worst, n)
|
|
return total, solved, worst
|
|
|
|
|
|
def main():
|
|
probs = load_roots()
|
|
print("=" * 82)
|
|
print("A. the false-position family, 1971 to 2011, plus the credit "
|
|
"projection")
|
|
print("=" * 82)
|
|
print(f"{'repair':22s}{'year':>7s}{'evals':>10s}{'solved':>10s}"
|
|
f"{'worst':>8s}")
|
|
print("-" * 82)
|
|
rows = []
|
|
for label, year, m in (
|
|
("Illinois", "1971",
|
|
lambda f, a, b, xt, mx: falsi_driver(f, a, b, xt, mx, "illinois")),
|
|
("Pegasus", "1972",
|
|
lambda f, a, b, xt, mx: falsi_driver(f, a, b, xt, mx, "pegasus")),
|
|
("Anderson-Bjorck", "1973",
|
|
lambda f, a, b, xt, mx: falsi_driver(f, a, b, xt, mx, "anderson")),
|
|
("Ford (Illinois-type)", "1995",
|
|
lambda f, a, b, xt, mx: falsi_modern(f, a, b, xt, mx, "ford")),
|
|
("Galdino", "2011",
|
|
lambda f, a, b, xt, mx: falsi_modern(f, a, b, xt, mx, "galdino")),
|
|
("credit projection", "here",
|
|
lambda f, a, b, xt, mx: falsi_driver(f, a, b, xt, mx, "credit"))):
|
|
t, s, w = score(m, probs)
|
|
rows.append((label, t, s))
|
|
print(f"{label:22s}{year:>7s}{t:>10d}{s:>7d}/199{w:>8d}", flush=True)
|
|
|
|
print("\n" + "=" * 82)
|
|
print("B. the sharp test: Chandrupatla's interpolation, its guard swapped")
|
|
print("=" * 82)
|
|
print("same interpolation, same bookkeeping, same termination; only the "
|
|
"guard differs\n")
|
|
print(f"{'guard':30s}{'evals':>10s}{'solved':>10s}{'worst':>8s}")
|
|
print("-" * 82)
|
|
res = {}
|
|
for label, guard in (("Chandrupatla's own test", "own"),
|
|
("none (clamp to bracket only)", "none"),
|
|
("the credit projection", "credit")):
|
|
t, s, w = score(lambda f, a, b, xt, mx, _g=guard:
|
|
chandrupatla_guard(f, a, b, xt, mx, _g), probs)
|
|
res[guard] = (t, s, w)
|
|
print(f"{label:30s}{t:>10d}{s:>7d}/199{w:>8d}", flush=True)
|
|
t_pub, s_pub, _w = score(lambda f, a, b, xt, mx: chandrupatla(f, a, b, xt,
|
|
mx), probs)
|
|
print(f"{'(published implementation)':30s}{t_pub:>10d}{s_pub:>7d}/199")
|
|
drift = abs(res['own'][0] - t_pub) / t_pub
|
|
print(f"\n control drift {drift * 100:.2f}% - the pluggable version "
|
|
f"must reproduce the published one")
|
|
better = res["credit"][0] < res["own"][0]
|
|
print(f" credit {res['credit'][0]} against its own guard "
|
|
f"{res['own'][0]} => "
|
|
f"{'the safeguard beats the state of the art' if better else 'IT DOES NOT'}")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|