"""Stage 12 - the validation runs added on revision. Five independent questions, each printed as it completes so the log can be read at any moment: 1. RE-RANKING with the baselines that were missing - Chandrupatla's own method and Dekker - against the proposed rule. 2. STABILITY to arithmetic noise. The margin over TOMS 748 is 0.5% and the paper itself reports that libm moves counts by a unit between platforms. Every f is perturbed by +-1 ulp at random, and the distribution of the total is reported for each method. Without this, "a tie" is rhetoric. 3. THE BUDGET MULTIPLIER. The method fixes it at 2 and argues that 1 is degenerate. Showing 1 fails is not showing 2 is right; the sweep over 1, 1.25, 1.5, 2, 3, 4 answers whether the result is flat in it. 4. A LARGER EVALUATION BUDGET. On Chandrupatla the guarantee alone can reach 160 evaluations against a budget of 200, so some "failures" may be the budget rather than the method. Re-run at 400. 5. STATISTICS. A sum of 2806 against 2820 is 14 evaluations over 154 problems. The head-to-head count gets an exact sign test and the total gets a bootstrap interval, so that "tie" and "win" are separated by evidence rather than by wording. Seeded; the seed is printed. Nothing here selects a configuration. """ from __future__ import annotations import math import random import sys import time from scipy.optimize._tstutils import (_APS_TESTS_DICTS, _CHANDRUPATLA_TESTS_DICTS) from simroot.baselines import BRENT_BASELINE, CHANDRUPATLA_BASELINE from simroot.baselines import (anderson_bjorck, best_of, bisection, brent, brent_scipy, chandrupatla, chandrupatla_scipy, dekker, itp, ridders, toms748) from simroot.pure import pure_root XTOL = 1e-12 ACC = 1e-6 SEED = 20260824 def ours(f, a, b, xtol, maxfev, budget_factor=2.0): return pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived3", guard="project", budget_factor=budget_factor) # Competitors are taken at the better of SciPy's implementation and ours, # problem by problem, so no result can be attributed to a weak reimplementation. # ITP has no library version, which is why it is validated against its own # theorem instead; the false-position family and bisection are not in SciPy. METHODS = [ ("OURS (similarity)", ours), ("chandrupatla", CHANDRUPATLA_BASELINE), ("toms748 (scipy)", toms748), ("brent", BRENT_BASELINE), ("itp", itp), ("ridders", ridders), ("anderson_bjorck", anderson_bjorck), ("dekker", dekker), ("bisection", bisection), ] def load(): out = [] for d in _APS_TESTS_DICTS: a, b = d["bracket"] args = d.get("args", ()) out.append(("APS", d["ID"], (lambda ff, aa: lambda x: ff(x, *aa))(d["f"], args), float(a), float(b), float(d["root"]))) for d in _CHANDRUPATLA_TESTS_DICTS: a, b = d["bracket"] out.append(("CHA", d["ID"], (lambda ff: lambda x: ff(x))(d["f"]), float(a), float(b), float(d["root"]))) return out def per_problem(method, probs, maxfev, criterion="root"): """Charged evaluation count per problem, and whether it was solved. criterion="root": within ACC of the published root (this stage's original scoring); "sign": the paper's single criterion since 2026-09-22, a sign change within 1e-6 or f(x) == 0.""" from corpus_roots_jl import solved out = {} for _s, pid, f, a, b, root in probs: try: x, n, conv = method(f, a, b, XTOL, maxfev) x = float(x) except Exception: x, n, conv = float("nan"), maxfev, False if criterion == "sign": ok = bool(conv and solved(f, x)) else: ok = bool(conv and math.isfinite(x) and abs(x - root) <= ACC) out[pid] = (n if ok else maxfev, ok) return out def totals(per, probs): t = {"APS": 0, "CHA": 0} ok = {"APS": 0, "CHA": 0} for s, pid, *_ in probs: t[s] += per[pid][0] ok[s] += per[pid][1] return t, ok def section(n, title): print("\n" + "=" * 78, flush=True) print(f"{n}. {title}", flush=True) print("=" * 78, flush=True) def main(): probs = load() t0 = time.time() print(f"seed {SEED}; {len(probs)} problems; started", flush=True) # ------------------------------------------------------------------ section(1, "RE-RANKING WITH THE MISSING BASELINES") base = {} for name, m in METHODS: p = per_problem(m, probs, 200) base[name] = p t, ok = totals(p, probs) print(f" {name:22s} APS {t['APS']:>6d} {ok['APS']:>3d}/154 " f"CHA {t['CHA']:>6d} {ok['CHA']:>3d}/45 " f"both {t['APS'] + t['CHA']:>6d}", flush=True) # ------------------------------------------------------------------ section(2, "STABILITY TO +-1 ULP PERTURBATION OF f") rng = random.Random(SEED) REAL = 40 print(f" {REAL} realisations, every evaluation of f nudged one ulp " f"either way\n", flush=True) def perturb(f, r): # The nudge is CACHED per abscissa. Drawing it afresh on every call # would model a noisy function rather than one computed to within an # ulp, and would penalise any method that asks about a point twice. seen = {} def g(x): if x in seen: return seen[x] v = f(x) if v != 0.0 and math.isfinite(v): v = math.nextafter(v, math.inf if r.random() < 0.5 else -math.inf) seen[x] = v return v return g dist = {name: [] for name, _ in METHODS} for k in range(REAL): r = random.Random(SEED + k) pp = [(s, pid, perturb(f, r), a, b, root) for s, pid, f, a, b, root in probs] for name, m in METHODS: t, _ok = totals(per_problem(m, pp, 200), pp) dist[name].append(t["APS"] + t["CHA"]) if (k + 1) % 5 == 0: print(f" ... {k + 1}/{REAL} realisations, " f"{time.time() - t0:.0f}s", flush=True) print(f"\n {'method':22s}{'min':>8s}{'median':>9s}{'max':>8s}" f"{'spread':>9s}{'unperturbed':>13s}", flush=True) print(" " + "-" * 70, flush=True) for name, _m in METHODS: v = sorted(dist[name]) t, _ = totals(base[name], probs) print(f" {name:22s}{v[0]:>8d}{v[len(v) // 2]:>9d}{v[-1]:>8d}" f"{v[-1] - v[0]:>9d}{t['APS'] + t['CHA']:>13d}", flush=True) print("\n If the spread of a method is comparable to the gap between two" "\n methods, that gap is not a result.", flush=True) # ------------------------------------------------------------------ section(3, "IS THE BUDGET MULTIPLIER 2 LOAD-BEARING?") print(f" {'factor':>8s}{'APS':>9s}{'APS ok':>9s}{'CHA':>9s}" f"{'CHA ok':>9s}{'both':>9s}", flush=True) print(" " + "-" * 54, flush=True) for fac in (1.0, 1.25, 1.5, 2.0, 3.0, 4.0): p = per_problem(lambda f, a, b, xt, mx, _f=fac: ours(f, a, b, xt, mx, _f), probs, 200) t, ok = totals(p, probs) mark = " <- the method" if fac == 2.0 else "" print(f" {fac:>8.2f}{t['APS']:>9d}{ok['APS']:>6d}/154" f"{t['CHA']:>9d}{ok['CHA']:>6d}/45" f"{t['APS'] + t['CHA']:>9d}{mark}", flush=True) # ------------------------------------------------------------------ section(4, "A LARGER EVALUATION BUDGET (400 INSTEAD OF 200)") print(f" {'method':22s}{'solved@200':>12s}{'solved@400':>12s}" f"{'evals@400':>12s}", flush=True) print(" " + "-" * 60, flush=True) for name, m in METHODS: p4 = per_problem(m, probs, 400) t4, ok4 = totals(p4, probs) _t2, ok2 = totals(base[name], probs) print(f" {name:22s}{ok2['APS'] + ok2['CHA']:>9d}/199" f"{ok4['APS'] + ok4['CHA']:>9d}/199" f"{t4['APS'] + t4['CHA']:>12d}", flush=True) # ------------------------------------------------------------------ section(5, "STATISTICS ON THE MARGINS") ourp = base["OURS (similarity)"] def sign_test(w, l): n = w + l if n == 0: return 1.0 k = min(w, l) p = sum(math.comb(n, i) for i in range(k + 1)) / 2.0 ** n return min(1.0, 2.0 * p) print(f" {'competitor':18s}{'W/T/L':>14s}{'sign test p':>14s}" f"{'total diff':>13s}{'bootstrap 95%':>22s}", flush=True) print(" " + "-" * 82, flush=True) for name, _m in METHODS: if name.startswith("OURS"): continue opp = base[name] w = t_ = l = 0 diffs = [] for _s, pid, *_ in probs: (no, oko), (nc, okc) = ourp[pid], opp[pid] if not (oko and okc): continue diffs.append(no - nc) if no < nc: w += 1 elif no > nc: l += 1 else: t_ += 1 if not diffs: continue p = sign_test(w, l) tot = sum(diffs) r = random.Random(SEED) boots = [] for _ in range(4000): boots.append(sum(diffs[r.randrange(len(diffs))] for _ in range(len(diffs)))) boots.sort() lo, hi = boots[int(0.025 * len(boots))], boots[int(0.975 * len(boots))] print(f" {name:18s}{f'{w}/{t_}/{l}':>14s}{p:>14.2e}" f"{tot:>13d}{f'[{lo}, {hi}]':>22s}", flush=True) print("\n Negative total difference means we spend fewer evaluations.", flush=True) print(" An interval straddling zero means the totals are indistinguishable.", flush=True) print(f"\nfinished in {time.time() - t0:.0f}s", flush=True) if __name__ == "__main__": main()