156 lines
6.1 KiB
Python
156 lines
6.1 KiB
Python
"""Stage 18 - attack the two claims that look like wins.
|
|
|
|
A claim is worth publishing only if it survives the objection a reader will
|
|
raise first. For each of the two candidate wins the objection is obvious, and
|
|
this file runs it rather than answering it in prose.
|
|
|
|
CLAIM A. The credit safeguard rescues unsafe interpolation rules: regula
|
|
falsi goes from 133 solved and 18799 evaluations to 193 and 8056.
|
|
OBJECTION: nobody uses raw regula falsi. It has been fixed since 1971 -
|
|
Illinois, Pegasus, Anderson-Bjorck - and those fixes are the real
|
|
comparison. Beating the unfixed version proves nothing.
|
|
TEST: the same rule, the same driver, the same everything, with the
|
|
stale-endpoint correction swapped for each classical fix in turn and for
|
|
the credit projection. Only the repair differs.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
from simroot.baselines import _opposite, _same_sign # 2026-09-21 fix
|
|
|
|
import math
|
|
|
|
import numpy as np
|
|
from scipy.optimize._tstutils import (_APS_TESTS_DICTS,
|
|
_CHANDRUPATLA_TESTS_DICTS)
|
|
|
|
from simroot.baselines import Counter, _Budget
|
|
from simroot.pure import _Safeguard
|
|
|
|
XTOL, MAXFEV, ACC = 1e-12, 200, 1e-6
|
|
|
|
|
|
# ==========================================================================
|
|
# CLAIM A: the credit projection against the classical repairs of regula falsi
|
|
# ==========================================================================
|
|
def falsi_driver(f, a, b, xtol, maxfev, repair):
|
|
"""Regula falsi with a pluggable repair for the stale endpoint.
|
|
|
|
Everything is shared: the interpolation, the bracket update, the budget,
|
|
the duplicate handling. `repair` is the only difference, so a difference
|
|
in the result cannot come from anywhere else.
|
|
|
|
"none" the unfixed rule
|
|
"illinois" halve the stale ordinate (Hafner/Dowell 1971)
|
|
"pegasus" f * fx / (fx + f_new) (Dowell/Jarratt 1972)
|
|
"anderson" f * (1 - f_new/f_moved) (Anderson/Bjorck 1973)
|
|
"credit" leave the ordinates alone; project the step onto the
|
|
interval the earned bisection budget allows
|
|
"""
|
|
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 _same_sign(fa, fb):
|
|
raise ValueError("no sign change")
|
|
sg = _Safeguard("credit", abs(b - a), xtol)
|
|
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 repair == "credit":
|
|
x_new = sg.project_credit(x_new, lo, hi)
|
|
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
|
|
w_before = abs(b - a)
|
|
fx = g(x_new)
|
|
if fx == 0.0:
|
|
return x_new, g.n, True
|
|
if _opposite(fa, fx):
|
|
if repair == "illinois":
|
|
fa *= 0.5
|
|
elif repair == "pegasus":
|
|
fa = fa * fb / (fb + fx)
|
|
elif repair == "anderson":
|
|
m = 1.0 - fx / fb
|
|
fa = fa * (m if m > 0.0 else 0.5)
|
|
b, fb = x_new, fx
|
|
else:
|
|
if repair == "illinois":
|
|
fb *= 0.5
|
|
elif repair == "pegasus":
|
|
fb = fb * fa / (fa + fx)
|
|
elif repair == "anderson":
|
|
m = 1.0 - fx / fa
|
|
fb = fb * (m if m > 0.0 else 0.5)
|
|
a, fa = x_new, fx
|
|
if repair == "credit":
|
|
sg.note_shrink(w_before, abs(b - a))
|
|
x = x_new
|
|
return 0.5 * (a + b), g.n, True
|
|
except _Budget:
|
|
return x, g.n, False
|
|
|
|
|
|
def load_roots():
|
|
out = []
|
|
for d in _APS_TESTS_DICTS:
|
|
a, b = d["bracket"]
|
|
args = d.get("args", ())
|
|
out.append(((lambda ff, aa: lambda t: ff(t, *aa))(d["f"], args),
|
|
float(a), float(b), float(d["root"])))
|
|
for d in _CHANDRUPATLA_TESTS_DICTS:
|
|
a, b = d["bracket"]
|
|
out.append(((lambda ff: lambda t: ff(t))(d["f"]),
|
|
float(a), float(b), float(d["root"])))
|
|
return out
|
|
|
|
|
|
def claim_a():
|
|
probs = load_roots()
|
|
print("=" * 84)
|
|
print("CLAIM A - the credit projection against the classical repairs")
|
|
print("=" * 84)
|
|
print("same interpolation, same driver, same budget; only the repair "
|
|
"differs\n")
|
|
print(f"{'repair':16s}{'year':>7s}{'evals':>10s}{'solved':>10s}"
|
|
f"{'worst':>8s}{'over bisect+3':>15s}")
|
|
print("-" * 84)
|
|
rows = {}
|
|
for repair, year in (("none", "300 BC"), ("illinois", "1971"),
|
|
("pegasus", "1972"), ("anderson", "1973"),
|
|
("credit", "here")):
|
|
total = solved = worst = over = 0
|
|
for f, a, b, root in probs:
|
|
try:
|
|
x, n, conv = falsi_driver(f, a, b, XTOL, MAXFEV, repair)
|
|
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)
|
|
if n > math.ceil(math.log2(abs(b - a) / XTOL)) + 4:
|
|
over += 1
|
|
rows[repair] = (total, solved, worst, over)
|
|
print(f"{repair:16s}{year:>7s}{total:>10d}{solved:>7d}/199"
|
|
f"{worst:>8d}{over:>15d}", flush=True)
|
|
best_classical = min(rows[r][0] for r in ("illinois", "pegasus",
|
|
"anderson"))
|
|
c = rows["credit"][0]
|
|
print(f"\n best classical repair {best_classical}, credit {c} "
|
|
f"=> {'SURVIVES' if c < best_classical else 'FAILS'} the objection")
|
|
return rows
|
|
|
|
|
|
if __name__ == "__main__":
|
|
claim_a()
|