343 lines
14 KiB
Python
343 lines
14 KiB
Python
"""Full audit of the algorithm: functional, logical, algorithmic.
|
|
|
|
Every finding below is DEMONSTRATED by a probe rather than asserted from
|
|
reading, because reading is how the last four defects got in. Each probe
|
|
prints PASS or the failure it provoked.
|
|
|
|
Sections
|
|
1 contracts and degenerate input
|
|
2 the bracketing invariant
|
|
3 the convergence claim: what "converged" is allowed to mean
|
|
4 the credit safeguard's arithmetic over its full domain
|
|
5 wasted evaluations: does the method ever re-ask a point it knows?
|
|
6 the competitors, since every conclusion rests on them
|
|
7 determinism
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
|
|
from simroot.baselines import (BASELINES, chandrupatla, chandrupatla_scipy,
|
|
dekker)
|
|
from simroot.pure import GUARDS, PURE_SIMS, _Safeguard, pure_root
|
|
|
|
XTOL, MAXFEV = 1e-12, 200
|
|
FINDINGS = []
|
|
|
|
|
|
def probe(name, fn):
|
|
try:
|
|
ok, detail = fn()
|
|
except Exception as exc:
|
|
ok, detail = False, f"raised {type(exc).__name__}: {exc}"
|
|
print(f" {'PASS' if ok else 'FAIL'} {name:56s} {detail}")
|
|
if not ok:
|
|
FINDINGS.append((name, detail))
|
|
|
|
|
|
def ours(f, a, b, guard="credit", **kw):
|
|
return pure_root(f, a, b, XTOL, MAXFEV, similarity="ratio",
|
|
scale="derived3", guard=guard, **kw)
|
|
|
|
|
|
# ==========================================================================
|
|
def sec1():
|
|
print("\n1. CONTRACTS AND DEGENERATE INPUT")
|
|
|
|
def reversed_bracket():
|
|
# a > b is not a bracket in the sense the docstring implies; what
|
|
# happens?
|
|
x, n, conv = ours(lambda t: t * t - 2.0, 2.0, 0.0)
|
|
inside = min(0.0, 2.0) <= x <= max(0.0, 2.0)
|
|
return inside, f"returned x={x:.6g} conv={conv} inside={inside}"
|
|
|
|
def no_sign_change():
|
|
try:
|
|
ours(lambda t: t * t + 1.0, 0.0, 2.0)
|
|
return False, "accepted a bracket with no sign change"
|
|
except ValueError:
|
|
return True, "ValueError, as it should"
|
|
|
|
def zero_tolerance():
|
|
x, n, conv = pure_root(lambda t: t * t - 2.0, 0.0, 2.0, 0.0, MAXFEV,
|
|
similarity="ratio", scale="derived3",
|
|
guard="credit")
|
|
return n <= MAXFEV, f"xtol=0 -> n={n} conv={conv} (must not hang)"
|
|
|
|
def negative_tolerance():
|
|
x, n, conv = pure_root(lambda t: t * t - 2.0, 0.0, 2.0, -1.0, MAXFEV,
|
|
similarity="ratio", scale="derived3",
|
|
guard="credit")
|
|
return n <= MAXFEV, f"xtol<0 -> n={n} conv={conv}"
|
|
|
|
def nan_function():
|
|
x, n, conv = ours(lambda t: float("nan"))if False else ours(
|
|
lambda t: float("nan") if 0.3 < t < 0.7 else t * t - 2.0, 0.0, 2.0)
|
|
return (n <= MAXFEV and not math.isnan(x)), \
|
|
f"f returns NaN on part of the bracket -> x={x} n={n} conv={conv}"
|
|
|
|
def infinite_function():
|
|
# a genuine pole INSIDE a genuine bracket: f(-1) < 0 < f(1) with a
|
|
# singularity at 0. The first version of this probe had both ends
|
|
# negative, so it measured the author's inattention rather than the
|
|
# method - the signature refused it, correctly.
|
|
def f(t):
|
|
if t == 0.0:
|
|
return math.inf
|
|
return 1.0 / t
|
|
x, n, conv = ours(f, -1.0, 1.0)
|
|
return n <= MAXFEV, (f"pole at 0 inside [-1,1] -> x={x:.6g} n={n} "
|
|
f"conv={conv} (no root exists; must not hang "
|
|
f"or crash)")
|
|
|
|
probe("reversed bracket stays inside", reversed_bracket)
|
|
probe("no sign change is refused", no_sign_change)
|
|
probe("xtol = 0 terminates", zero_tolerance)
|
|
probe("xtol < 0 terminates", negative_tolerance)
|
|
probe("NaN inside the bracket", nan_function)
|
|
probe("pole inside the bracket", infinite_function)
|
|
|
|
|
|
# ==========================================================================
|
|
def sec2():
|
|
print("\n2. THE BRACKETING INVARIANT")
|
|
|
|
fs = [(lambda t: t * t - 2.0, 0.0, 2.0),
|
|
(lambda t: math.exp(-t) * (t - 1.0), -1.0, 40.0),
|
|
(lambda t: t ** 10 - 1.0, 0.0, 1.3),
|
|
(lambda t: math.atan(t), -1.0, 3.0)]
|
|
|
|
def never_leaves():
|
|
bad = []
|
|
for guard in GUARDS:
|
|
for f, a, b in fs:
|
|
probes = []
|
|
def g(t, _p=probes, _f=f):
|
|
_p.append(t)
|
|
return _f(t)
|
|
try:
|
|
ours(g, a, b, guard=guard)
|
|
except Exception as exc:
|
|
bad.append((guard, repr(exc)[:40]))
|
|
continue
|
|
lo, hi = min(a, b), max(a, b)
|
|
out = [t for t in probes if not (lo <= t <= hi)]
|
|
if out:
|
|
bad.append((guard, f"probed outside: {out[:3]}"))
|
|
return not bad, ("no probe left the initial bracket"
|
|
if not bad else str(bad[:3]))
|
|
|
|
def sign_change_preserved():
|
|
"""Every returned interval must still straddle the root."""
|
|
bad = []
|
|
for guard in GUARDS:
|
|
for f, a, b in fs:
|
|
x, n, conv = ours(f, a, b, guard=guard)
|
|
if conv and abs(f(x)) > 1e-6 and abs(x - a) > 1e-9:
|
|
# allow the tolerance itself, but a converged answer whose
|
|
# |f| is large is a false claim
|
|
bad.append((guard, x, f(x)))
|
|
return not bad, ("converged answers all have small |f|"
|
|
if not bad else str(bad[:3]))
|
|
|
|
probe("no probe ever leaves the initial bracket", never_leaves)
|
|
probe("a converged answer really is near a root", sign_change_preserved)
|
|
|
|
|
|
# ==========================================================================
|
|
def sec3():
|
|
print("\n3. WHAT 'CONVERGED' IS ALLOWED TO MEAN")
|
|
|
|
def converged_implies_tolerance():
|
|
"""conv=True must imply the final bracket is within xtol - otherwise
|
|
the flag is an opinion."""
|
|
bad = []
|
|
for guard in GUARDS:
|
|
f = lambda t: math.exp(-t) * (t - 1.0)
|
|
x, n, conv = pure_root(f, -1.0, 40.0, 1e-300, 400,
|
|
similarity="ratio", scale="derived3",
|
|
guard=guard)
|
|
if conv and abs(x - 1.0) > 1e-8:
|
|
bad.append((guard, x, n))
|
|
return not bad, ("an impossible tolerance is reported as failure"
|
|
if not bad else str(bad))
|
|
|
|
def budget_respected():
|
|
bad = []
|
|
for guard in GUARDS:
|
|
for budget in (0, 1, 2, 3, 7, 25):
|
|
_x, n, _c = pure_root(lambda t: t * t - 2.0, 0.0, 2.0, XTOL,
|
|
budget, similarity="ratio",
|
|
scale="derived3", guard=guard)
|
|
if n > budget:
|
|
bad.append((guard, budget, n))
|
|
return not bad, "never spends more than the budget" if not bad else str(bad)
|
|
|
|
probe("conv=True implies the tolerance was met", converged_implies_tolerance)
|
|
probe("the evaluation budget is never exceeded", budget_respected)
|
|
|
|
|
|
# ==========================================================================
|
|
def sec4():
|
|
print("\n4. THE CREDIT SAFEGUARD OVER ITS FULL DOMAIN")
|
|
|
|
def credit_overflow():
|
|
"""credit accumulates log2 of every shrink. 2**(C-1) overflows a
|
|
double at C ~ 1025. Is C bounded in the implementation?"""
|
|
sg = _Safeguard("credit", 1.0, XTOL)
|
|
sg.credit = 1100.0
|
|
try:
|
|
sg.project_credit(0.9, 0.0, 1.0)
|
|
return True, "no overflow at C=1100"
|
|
except OverflowError as exc:
|
|
return False, f"OverflowError at C=1100 ({exc}); C is unbounded"
|
|
|
|
def credit_reachable():
|
|
"""A large credit is REACHABLE - sixty nine-decade steps reach 1012 -
|
|
so what matters is that it is harmless, not that it is small. An
|
|
earlier version of this probe asserted C < 1000, which after the fix
|
|
reports a defect that no longer exists."""
|
|
sg = _Safeguard("credit", 1.0, XTOL)
|
|
sg.credit = 1.0
|
|
w = 1.0
|
|
for _ in range(60):
|
|
sg.note_shrink(w, w * 1e-9) # nine decades per step
|
|
w *= 1e-9
|
|
big = sg.credit
|
|
x = sg.project_credit(0.9, 0.0, 1.0)
|
|
return x == 0.9, (f"C reaches {big:.0f} and a step is then "
|
|
f"unrestricted (x={x}), which is the correct "
|
|
f"behaviour: credit that large restricts nothing")
|
|
|
|
def wide_dynamic_range():
|
|
"""An end-to-end problem with a huge bracket and a tight tolerance."""
|
|
f = lambda t: t - 1.0
|
|
x, n, conv = pure_root(f, -1e150, 1e150, 1e-12, 4000,
|
|
similarity="ratio", scale="derived3",
|
|
guard="credit")
|
|
return (conv and abs(x - 1.0) < 1e-6), \
|
|
f"bracket 1e150 wide -> x={x:.6g} n={n} conv={conv}"
|
|
|
|
probe("credit cannot overflow the exponent", credit_overflow)
|
|
probe("credit stays in a reachable range", credit_reachable)
|
|
probe("a 1e150-wide bracket is handled", wide_dynamic_range)
|
|
|
|
|
|
# ==========================================================================
|
|
def sec5():
|
|
print("\n5. WASTED EVALUATIONS")
|
|
|
|
def repeats():
|
|
"""The duplicate check compares only against the last two points. Is
|
|
a third-back point ever re-evaluated?"""
|
|
worst = 0
|
|
where = None
|
|
from scipy.optimize._tstutils import _APS_TESTS_DICTS, \
|
|
_CHANDRUPATLA_TESTS_DICTS
|
|
probs = [((lambda ff, aa: lambda t: ff(t, *aa))(d["f"],
|
|
d.get("args", ())),
|
|
float(d["bracket"][0]), float(d["bracket"][1]), d["ID"])
|
|
for d in _APS_TESTS_DICTS]
|
|
probs += [((lambda ff: lambda t: ff(t))(d["f"]),
|
|
float(d["bracket"][0]), float(d["bracket"][1]), d["ID"])
|
|
for d in _CHANDRUPATLA_TESTS_DICTS]
|
|
for f, a, b, pid in probs:
|
|
seen = []
|
|
def g(t, _s=seen, _f=f):
|
|
_s.append(t)
|
|
return _f(t)
|
|
try:
|
|
ours(g, a, b)
|
|
except Exception:
|
|
continue
|
|
dup = len(seen) - len(set(seen))
|
|
if dup > worst:
|
|
worst, where = dup, pid
|
|
return worst == 0, (f"no repeated evaluation" if worst == 0
|
|
else f"{worst} repeats, worst on {where}")
|
|
|
|
probe("no point is evaluated twice", repeats)
|
|
|
|
|
|
# ==========================================================================
|
|
def sec6():
|
|
print("\n6. THE COMPETITORS")
|
|
|
|
def dekker_sane():
|
|
"""Dekker solves 24/45 on Chandrupatla, which is either its known slow
|
|
convergence or a bug in our transcription. It has no library version,
|
|
so the check is against its defining property: it must never be worse
|
|
than its own secant/bisection choice on an easy problem."""
|
|
x, n, conv = dekker(lambda t: t * t - 2.0, 0.0, 2.0, XTOL, MAXFEV)
|
|
return (conv and abs(x - math.sqrt(2.0)) < 1e-8 and n < 30), \
|
|
f"x^2-2 -> n={n} conv={conv} x={x:.12f}"
|
|
|
|
def chandrupatla_agrees_with_scipy():
|
|
from scipy.optimize._tstutils import _APS_TESTS_DICTS
|
|
diff = 0
|
|
for d in _APS_TESTS_DICTS:
|
|
a, b = float(d["bracket"][0]), float(d["bracket"][1])
|
|
args = d.get("args", ())
|
|
f = (lambda ff, aa: lambda t: ff(t, *aa))(d["f"], args)
|
|
x1, n1, c1 = chandrupatla(f, a, b, XTOL, MAXFEV)
|
|
x2, n2, c2 = chandrupatla_scipy(f, a, b, XTOL, MAXFEV)
|
|
if c1 and c2 and abs(x1 - x2) > 1e-6:
|
|
diff += 1
|
|
return diff == 0, f"{diff} problems where the two disagree on the root"
|
|
|
|
def exact_zero_handling():
|
|
"""f(x) == 0 exactly: every method must return that x, not keep going."""
|
|
bad = []
|
|
for name, m in BASELINES.items():
|
|
x, n, conv = m(lambda t: t - 1.0, 1.0, 5.0, XTOL, MAXFEV)
|
|
if not (conv and abs(x - 1.0) < 1e-12 and n <= 3):
|
|
bad.append((name, n, conv, x))
|
|
return not bad, "all return the exact root at once" if not bad else str(bad)
|
|
|
|
probe("dekker behaves on an easy problem", dekker_sane)
|
|
probe("our chandrupatla agrees with scipy's root", chandrupatla_agrees_with_scipy)
|
|
probe("an exact zero at an endpoint costs <= 3 evaluations", exact_zero_handling)
|
|
|
|
|
|
# ==========================================================================
|
|
def sec7():
|
|
print("\n7. DETERMINISM")
|
|
|
|
def repeatable():
|
|
f = lambda t: math.cos(t) - t
|
|
runs = {ours(f, 0.0, 1.0) for _ in range(5)}
|
|
return len(runs) == 1, f"{len(runs)} distinct results over 5 runs"
|
|
|
|
def similarity_independence():
|
|
"""Changing the similarity must change only the step, never the
|
|
contract: all nine must stay inside the bracket and respect the
|
|
budget."""
|
|
bad = []
|
|
for sim in PURE_SIMS:
|
|
x, n, conv = pure_root(lambda t: t * t - 2.0, 0.0, 2.0, XTOL,
|
|
MAXFEV, similarity=sim, scale="derived3",
|
|
guard="credit")
|
|
if not (0.0 <= x <= 2.0 and n <= MAXFEV):
|
|
bad.append((sim, x, n))
|
|
return not bad, "all nine respect the contract" if not bad else str(bad)
|
|
|
|
probe("the same call gives the same answer", repeatable)
|
|
probe("every similarity respects the contract", similarity_independence)
|
|
|
|
|
|
def main():
|
|
print("=" * 78)
|
|
print("ALGORITHM AUDIT")
|
|
print("=" * 78)
|
|
for s in (sec1, sec2, sec3, sec4, sec5, sec6, sec7):
|
|
s()
|
|
print("\n" + "=" * 78)
|
|
print(f"FINDINGS: {len(FINDINGS)}")
|
|
print("=" * 78)
|
|
for name, detail in FINDINGS:
|
|
print(f" - {name}: {detail}")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|