Files
sim-numerics/audit_chandrupatla_library.py
vvsandClaude Opus 5 de90c7a089 Audit on Valeri's challenge: our Chandrupatla implementation stops below 2*xtol (tl = tol/span > 1/2) while every other method, SciPy's included, stops at xtol; the protocol corrected that factor for ITP and not here, and best_of let the looser run win. Corrected it costs Chandrupatla 2 evaluations on APS (2785 -> 2787) and 18 on CHA (2423 -> 2441): real, in the competitor's favour, and not what changes the picture. The public scipy.optimize.elementwise.find_root, the private _chandrupatla and ours corrected agree on all 199 problems exactly.
What the tables' totals hid: problem by problem the scheduled projection beats Chandrupatla on APS 87/24/42, exact sign test p = 9e-5, while trailing by 0.7% in total; on CHA Chandrupatla is genuinely better (26 to 9, p = 0.006). The earned projection ties him on APS problems (79/3/71, p = 0.57). Benchmark not yet switched to the library baseline - awaiting Valeri's decision together with the framing of the proposed method.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
2026-09-21 10:14:42 +01:00

62 lines
2.6 KiB
Python

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