"""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()