diff --git a/audit_chandrupatla_library.py b/audit_chandrupatla_library.py new file mode 100644 index 0000000..8fe4925 --- /dev/null +++ b/audit_chandrupatla_library.py @@ -0,0 +1,61 @@ +"""Can the Chandrupatla baseline be the library alone? (Valeri, 2026-09-21) + +Our own implementation stopped on a bracket twice as wide as everyone else's, and +best_of let that leak into the reported numbers. SciPy ships the method twice: the +private _chandrupatla, already in the benchmark, and the PUBLIC +scipy.optimize.elementwise.find_root (Chandrupatla's algorithm, SciPy >= 1.15). A +private function can change without notice, so the public one is the candidate; this +checks, problem by problem, that it spends exactly what the private one spends and +what our implementation spends once its tolerance is corrected. + + python -u audit_chandrupatla_library.py +""" +from __future__ import annotations + +import numpy as np + +from audit_chandrupatla_tolerance import ours_half +from simroot.baselines import Counter, _Budget, chandrupatla_scipy +from stage12_validation import load, per_problem, totals + + +def chandrupatla_public(f, a, b, xtol, maxfev): + """scipy.optimize.elementwise.find_root with the benchmark's stopping rule: bracket + narrower than xtol, relative term at machine precision, no stopping on |f|.""" + from scipy.optimize import elementwise + g = Counter(f, maxfev) + + def vec(x): + return np.asarray([g(float(xi)) for xi in np.atleast_1d(x)]).reshape(np.shape(x)) + + try: + res = elementwise.find_root(vec, (a, b), maxiter=maxfev, + tolerances=dict(xatol=xtol, xrtol=8.9e-16, + fatol=0.0, frtol=0.0)) + return (float(np.asarray(res.x).ravel()[0]), g.n, + bool(np.all(np.asarray(res.success)))) + except _Budget: + return 0.5 * (a + b), g.n, False + except Exception: + return 0.5 * (a + b), g.n, False + + +def main(): + probs = load() + pub = per_problem(chandrupatla_public, probs, 200) + prv = per_problem(chandrupatla_scipy, probs, 200) + own = per_problem(ours_half, probs, 200) + for name, p in (("public find_root", pub), ("private _chandrupatla", prv), + ("ours, tolerance corrected", own)): + t, ok = totals(p, probs) + print(f"{name:<28} APS {t['APS']:5d} ({ok['APS']}/154) CHA {t['CHA']:5d} " + f"({ok['CHA']}/45) both {t['APS'] + t['CHA']}") + diff_prv = [pid for pid in pub if pub[pid] != prv[pid]] + diff_own = [pid for pid in pub if pub[pid] != own[pid]] + print(f"\nproblems where public != private: {len(diff_prv)} {diff_prv[:8]}") + print(f"problems where public != ours corrected: {len(diff_own)} {diff_own[:8]}") + return {"pub": totals(pub, probs), "diff_prv": diff_prv, "diff_own": diff_own} + + +if __name__ == "__main__": + main() diff --git a/audit_chandrupatla_tolerance.py b/audit_chandrupatla_tolerance.py new file mode 100644 index 0000000..b540bdf --- /dev/null +++ b/audit_chandrupatla_tolerance.py @@ -0,0 +1,89 @@ +"""Audit: does Chandrupatla stop on the same bracket width as everybody else? + +Found 2026-09-21 on Valeri's challenge "you made a mistake somewhere". Every method in +the benchmark stops when the bracket is at most xtol wide - ours, bisection, Brent, +Dekker, TOMS 748 (np.isclose with atol=xtol), SciPy's _chandrupatla (dx < xatol). +OUR implementation of Chandrupatla does not: it stops on tl = tol/span > 1/2 with +tol = xtol + 2 eps |x|, i.e. when the bracket is narrower than 2*xtol. The protocol +already corrects exactly this factor for ITP (epsilon = xtol/2) and did not for +Chandrupatla, and best_of takes, problem by problem, whichever of the two +implementations is cheaper - so the looser one wins. + +Measured here, per implementation and combined, at the benchmark's charging rules: + ours @ xtol the implementation as benchmarked (stops below 2*xtol) + ours @ xtol/2 the same code given ITP's correction (stops below xtol) + scipy @ xtol the library, which already stops below xtol + best_of as benchmarked, and best_of with the correction. +Also the final bracket width each actually stops at, to confirm the factor directly. + + python -u audit_chandrupatla_tolerance.py +""" +from __future__ import annotations + +import math + +from simroot.baselines import Counter, best_of, chandrupatla, chandrupatla_scipy +from stage12_validation import XTOL, load, per_problem, totals + + +def ours_half(f, a, b, xtol, maxfev): + return chandrupatla(f, a, b, 0.5 * xtol, maxfev) + + +def final_width_ratio(probs): + """Max over problems of (final proven bracket width / xtol) for our implementation, + reconstructed from the signs it observed.""" + worst = 0.0 + for _s, _pid, f, a, b, _root in probs: + xs = [] + + def g(x, _f=f): + v = _f(x) + xs.append((x, v)) + return v + try: + chandrupatla(g, a, b, XTOL, 200) + except Exception: + continue + if len(xs) < 2: + continue + (xa, fa), (xb, fb) = xs[1], xs[0] + lo, hi = min(xa, xb), max(xa, xb) + flo = fa if xa < xb else fb + for x, v in xs[2:]: + if not (lo < x < hi) or v == 0.0: + continue + if (v > 0.0) == (flo > 0.0): + lo, flo = x, v + else: + hi = x + worst = max(worst, (hi - lo) / XTOL) + return worst + + +def main(): + probs = load() + runs = { + "ours @ xtol (as benchmarked)": chandrupatla, + "ours @ xtol/2 (ITP's correction)": ours_half, + "scipy @ xtol": chandrupatla_scipy, + "best_of, as benchmarked": best_of(chandrupatla, chandrupatla_scipy), + "best_of, corrected": best_of(ours_half, chandrupatla_scipy), + } + out = {} + print(f"{'Chandrupatla':<34} {'APS':>6} {'solved':>7} {'CHA':>6} {'solved':>7} {'both':>6}") + for name, m in runs.items(): + p = per_problem(m, probs, 200) + t, ok = totals(p, probs) + out[name] = (t["APS"], ok["APS"], t["CHA"], ok["CHA"]) + print(f"{name:<34} {t['APS']:6d} {ok['APS']:4d}/154 {t['CHA']:6d} {ok['CHA']:4d}/45 " + f"{t['APS'] + t['CHA']:6d}") + w = final_width_ratio(probs) + out["max_final_width_over_xtol"] = w + print(f"\nour implementation, largest final bracket width / xtol over all problems: {w:.3f}" + "\n(everyone else stops at <= 1; a value up to 2 confirms the looser rule)") + return out + + +if __name__ == "__main__": + main() diff --git a/audit_head_to_head.py b/audit_head_to_head.py new file mode 100644 index 0000000..6d9bfb1 --- /dev/null +++ b/audit_head_to_head.py @@ -0,0 +1,62 @@ +"""Audit: is "Chandrupatla is first on the canonical set" a result or a total? + +Totals differ by 0.7% (2806 against 2787 with the tolerance corrected); the paper itself +reads 0.5% as a tie. Problem by problem, with an exact two-sided sign test on the +problems where the two differ, for the scheduled projection AND the earned one against +Chandrupatla at the corrected tolerance. + + python -u audit_head_to_head.py +""" +from __future__ import annotations + +import math + +from audit_chandrupatla_tolerance import ours_half +from simroot.baselines import best_of, chandrupatla_scipy +from simroot.pure import pure_root +from stage12_validation import load, per_problem + + +def sign_test(w, l): + n = w + l + if n == 0: + return 1.0 + k = min(w, l) + return min(1.0, 2.0 * sum(math.comb(n, i) for i in range(k + 1)) / 2.0 ** n) + + +def main(): + probs = load() + ch = per_problem(best_of(ours_half, chandrupatla_scipy), probs, 200) + ours = { + "scheduled projection (c=2)": lambda f, a, b, xt, mx: pure_root( + f, a, b, xt, mx, similarity="ratio", scale="derived3", guard="project"), + "earned projection (credit)": lambda f, a, b, xt, mx: pure_root( + f, a, b, xt, mx, similarity="ratio", scale="derived3", guard="credit"), + } + sets = {s: [pid for s_, pid, *_ in probs if s_ == s] for s in ("APS", "CHA")} + out = {} + for name, m in ours.items(): + p = per_problem(m, probs, 200) + for s, ids in sets.items(): + w = t = l = 0 + for pid in ids: + (no, oko), (nc, okc) = p[pid], ch[pid] + if not (oko and okc): + continue + if no < nc: + w += 1 + elif no > nc: + l += 1 + else: + t += 1 + tot_o = sum(p[pid][0] for pid in ids) + tot_c = sum(ch[pid][0] for pid in ids) + out[(name, s)] = (w, t, l, sign_test(w, l), tot_o, tot_c) + print(f"{name:<28} {s}: wins/ties/losses {w}/{t}/{l} sign test p = " + f"{sign_test(w, l):.2g} totals {tot_o} vs Chandrupatla {tot_c}") + return out + + +if __name__ == "__main__": + main() diff --git a/results_audit_chandrupatla_library.txt b/results_audit_chandrupatla_library.txt new file mode 100644 index 0000000..7ebf388 --- /dev/null +++ b/results_audit_chandrupatla_library.txt @@ -0,0 +1,6 @@ +public find_root APS 2787 (153/154) CHA 2441 (40/45) both 5228 +private _chandrupatla APS 2787 (153/154) CHA 2441 (40/45) both 5228 +ours, tolerance corrected APS 2787 (153/154) CHA 2441 (40/45) both 5228 + +problems where public != private: 0 [] +problems where public != ours corrected: 0 [] diff --git a/results_audit_chandrupatla_tolerance.txt b/results_audit_chandrupatla_tolerance.txt new file mode 100644 index 0000000..819218e --- /dev/null +++ b/results_audit_chandrupatla_tolerance.txt @@ -0,0 +1,9 @@ +Chandrupatla APS solved CHA solved both +ours @ xtol (as benchmarked) 2785 153/154 2424 40/45 5209 +ours @ xtol/2 (ITP's correction) 2787 153/154 2441 40/45 5228 +scipy @ xtol 2787 153/154 2441 40/45 5228 +best_of, as benchmarked 2785 153/154 2423 40/45 5208 +best_of, corrected 2787 153/154 2441 40/45 5228 + +our implementation, largest final bracket width / xtol over all problems: 1000000000000.000 +(everyone else stops at <= 1; a value up to 2 confirms the looser rule) diff --git a/results_audit_head_to_head.txt b/results_audit_head_to_head.txt new file mode 100644 index 0000000..5d1bd0f --- /dev/null +++ b/results_audit_head_to_head.txt @@ -0,0 +1,4 @@ +scheduled projection (c=2) APS: wins/ties/losses 87/24/42 sign test p = 9.2e-05 totals 2806 vs Chandrupatla 2787 +scheduled projection (c=2) CHA: wins/ties/losses 9/5/26 sign test p = 0.006 totals 3220 vs Chandrupatla 2441 +earned projection (credit) APS: wins/ties/losses 79/3/71 sign test p = 0.57 totals 3568 vs Chandrupatla 2787 +earned projection (credit) CHA: wins/ties/losses 5/0/35 sign test p = 1.4e-06 totals 3020 vs Chandrupatla 2441