"""The full benchmark: similarity root finding against the state of the art. Problem sets, both taken from scipy so that nothing is transcribed by hand: APS 154 problems - Alefeld, Potra & Shi, "Algorithm 748: Enclosing Zeros of Continuous Functions", ACM TOMS 21(3), 1995. The canonical set for bracketed root finding. Chandrupatla 45 problems - Chandrupatla, "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives", Adv. Eng. Software 28, 1997. Cost model: number of evaluations of f, which is the only cost that matters for a bracketed solver. A run counts as SOLVED when it reports convergence AND lands within `acc` of the published root. Reporting rules, chosen so the numbers cannot flatter the new method: * failures are charged at the full evaluation budget, so a method cannot buy a low total by giving up early; * the total is reported over ALL problems, not over the ones a method happened to solve; * ITP is reported twice - at the authors' recommended constants, and as an ORACLE that is given the best of a (kappa1, n0) grid per problem, i.e. tuned individually to every problem it is asked to solve. The oracle is not an achievable method; it is there so that no one can say the state of the art was handicapped; * every method that carries tuning constants is marked in the table. """ from __future__ import annotations import math import sys from scipy.optimize._tstutils import _APS_TESTS_DICTS, _CHANDRUPATLA_TESTS_DICTS from simroot.baselines import (anderson_bjorck, bisection, brent, illinois, itp, pegasus, regula_falsi, ridders, toms748) from simroot.pure import pure_root from simroot.sim import sim_root XTOL = 1e-12 MAXFEV = 200 ACC = 1e-6 # a run must land this close to the published root def load(which): out = [] if which in ("aps", "all"): for d in _APS_TESTS_DICTS: a, b = d["bracket"] args = d.get("args", ()) f = d["f"] out.append((d["ID"], (lambda ff, aa: lambda x: ff(x, *aa))(f, args), float(a), float(b), float(d["root"]))) if which in ("chandrupatla", "all"): for d in _CHANDRUPATLA_TESTS_DICTS: a, b = d["bracket"] f = d["f"] out.append((d["ID"], (lambda ff: lambda x: ff(x))(f), float(a), float(b), float(d["root"]))) return out ITP_GRID = [(fac, n0) for fac in (1.0, 0.1, 0.01, 1e-4) for n0 in (0, 1, 2, 5, 10)] def itp_oracle(f, a, b, xtol, maxfev): """The most generous ITP possible: the best of the grid, per problem.""" best = None for fac, n0 in ITP_GRID: try: x, n, conv = itp(f, a, b, xtol, maxfev, kappa1=fac * 0.2 / (b - a), n0=n0) except ValueError: continue if conv and (best is None or n < best[1]): best = (x, n, True) return best if best is not None else (0.5 * (a + b), maxfev, False) # name -> (callable, has_tuning_constants) METHODS = { "bisection": (bisection, False), "regula_falsi": (regula_falsi, False), "illinois": (illinois, True), "pegasus": (pegasus, True), "anderson_bjorck": (anderson_bjorck, True), "ridders": (ridders, False), "brent": (brent, True), "toms748": (toms748, True), "itp (recommended)": (itp, True), "itp (per-problem oracle)": (itp_oracle, True), "PURE.ratio.step": (lambda f, a, b, xt, mx: pure_root( f, a, b, xt, mx, similarity="ratio", scale="derived", guard="step"), False), "PURE.exp_abs.step": (lambda f, a, b, xt, mx: pure_root( f, a, b, xt, mx, similarity="exp_abs", scale="derived", guard="step"), False), "PURE.rational.step": (lambda f, a, b, xt, mx: pure_root( f, a, b, xt, mx, similarity="rational", scale="derived", guard="step"), False), "PURE.ratio.both": (lambda f, a, b, xt, mx: pure_root( f, a, b, xt, mx, similarity="ratio", scale="derived", guard="both"), False), "PURE.ratio.bracket": (lambda f, a, b, xt, mx: pure_root( f, a, b, xt, mx, similarity="ratio", scale="derived", guard="bracket"), False), "PURE.ratio.none": (lambda f, a, b, xt, mx: pure_root( f, a, b, xt, mx, similarity="ratio", scale="derived", guard="none"), False), "SIM.pair2d.inv": (lambda f, a, b, xt, mx: sim_root( f, a, b, xt, mx, similarity="inv", prototypes="pair2d", guard=True), False), "SIM.pair2d.neg_sq": (lambda f, a, b, xt, mx: sim_root( f, a, b, xt, mx, similarity="neg_sq", prototypes="pair2d", guard=True), False), } def run(problems, label): results = {} for name, (m, _tuned) in METHODS.items(): total, solved, fails, per, over = 0, 0, [], {}, [] for pid, f, a, b, root in problems: note = "" try: x, n, conv = m(f, a, b, XTOL, MAXFEV) except Exception as exc: # a crash is a failure, loudly x, n, conv = float("nan"), MAXFEV, False note = repr(exc)[:60] ok = bool(conv) and math.isfinite(x) and abs(x - root) <= ACC # failures are charged the full budget total += n if ok else MAXFEV per[pid] = (n, ok, note) if ok: solved += 1 else: fails.append(pid) # Does the run respect the "no worse than twice bisection" bound? # Only SOLVED runs can be judged against it: a run that failed sat # at the evaluation cap, which says nothing about the bound and # would silently fill this column with budget exhaustion. bound = 2 * math.ceil(math.log2(abs(b - a) / XTOL)) + 4 if ok and n > bound: over.append(pid) results[name] = dict(total=total, solved=solved, n=len(problems), fails=fails, per=per, over=over) print(f"\n{'=' * 96}\n{label} ({len(problems)} problems, " f"xtol={XTOL:g}, budget={MAXFEV} evaluations, accepted within " f"{ACC:g} of the published root)\n{'=' * 96}") print(f"{'method':28s}{'tuned':>7s}{'solved':>9s}{'failed':>8s}" f"{'total fev':>11s}{'mean/solved':>13s}{'over 2x bisect':>16s}") print("-" * 96) order = sorted(results, key=lambda k: (len(results[k]["fails"]), results[k]["total"])) for name in order: r = results[name] tuned = "yes" if METHODS[name][1] else "no" mean = (sum(v[0] for v in r["per"].values() if v[1]) / r["solved"] if r["solved"] else float("nan")) print(f"{name:28s}{tuned:>7s}{r['solved']:>6d}/{r['n']:<3d}" f"{len(r['fails']):>8d}{r['total']:>11d}{mean:>13.2f}" f"{len(r['over']):>16d}") return results def head_to_head(results, a_name, b_name, problems): ra, rb = results[a_name], results[b_name] win = loss = tie = 0 both = 0 for pid, _f, _a, _b, _root in problems: na, oka, _ = ra["per"][pid] nb, okb, _ = rb["per"][pid] if not (oka and okb): continue both += 1 if na < nb: win += 1 elif na > nb: loss += 1 else: tie += 1 print(f"\n{a_name} vs {b_name}, on the {both} problems both solve: " f"{win} wins / {tie} ties / {loss} losses") def main(): which = sys.argv[1] if len(sys.argv) > 1 else "all" aps = load("aps") cha = load("chandrupatla") if which in ("aps", "all"): r = run(aps, "Alefeld-Potra-Shi 1995") head_to_head(r, "PURE.ratio.step", "brent", aps) head_to_head(r, "PURE.ratio.step", "toms748", aps) head_to_head(r, "PURE.ratio.step", "itp (per-problem oracle)", aps) head_to_head(r, "SIM.pair2d.inv", "brent", aps) for nm in ("PURE.ratio.step", "brent", "toms748"): if r[nm]["fails"]: print(f" {nm} failed on: {', '.join(r[nm]['fails'][:12])}" + (" ..." if len(r[nm]["fails"]) > 12 else "")) if which in ("chandrupatla", "all"): r2 = run(cha, "Chandrupatla 1997") head_to_head(r2, "PURE.ratio.step", "brent", cha) if __name__ == "__main__": main()