vvsandClaude Opus 5
1b6c626f67
(2) The proposed method is again the scheduled projection (2806/3220/6026), its budget multiplier DECLARED as the one free constant: the canonical total is the same for every multiplier 1.25-4, the wide-bracket one is not. The earned projection is the constant-free variant (3568/3020/6588). What the old tables hid: problem by problem the proposed rule beats Chandrupatla on the canonical set 87/24/42, exact sign test p = 9e-5, while trailing him by 0.7% in total; on the wide-bracket set he is genuinely stronger (26 to 9, p = 0.006). Gap decomposition recomputed with the library: 244 (was 262), 160 of it on the flat families (two thirds), fun1 58 (a quarter, not a third as the abstract said). Also corrected on the way: the 'subsumed by Chandrupatla' pin claimed he is better on the worst single problem - true only with the looser tolerance (84 against the credit rule's 83 once fixed); §8.3 still called zero slack degenerate, which Table 1a refutes; §5.2's Brent figure is our pure-Python implementation and now says so. verify_article.py 286/286, 802 tests pass twice. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
269 lines
9.4 KiB
Python
269 lines
9.4 KiB
Python
"""Stage 13 - the whole thing recomputed, asking one question:
|
|
|
|
is there ANY axis on which this method is ahead?
|
|
|
|
The ranking tables became unreadable because 216 of the 228 configurations in
|
|
them are variants of our own idea, so a place number says almost nothing.
|
|
Here the field is the proposed rule and the real competitors only, and every
|
|
axis ends in a plain verdict: AHEAD, LEVEL or BEHIND.
|
|
|
|
"LEVEL" is not a courtesy. It means the gap is inside the spread measured by
|
|
perturbing every evaluation of f by one ulp - if arithmetic noise can move the
|
|
total by more than the gap, the gap is not a result.
|
|
|
|
Axes:
|
|
A evaluations of f, each set and both
|
|
B problems solved
|
|
C head to head, problem by problem, with an exact sign test
|
|
D worst single problem
|
|
E robustness to +-1 ulp
|
|
F worst-case guarantee
|
|
G free constants, counted honestly for every method including ours
|
|
H wall-clock time
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
import random
|
|
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 (best_of, bisection, brent, brent_scipy,
|
|
chandrupatla, chandrupatla_scipy, itp, ridders,
|
|
toms748)
|
|
from simroot.pure import pure_root
|
|
|
|
XTOL, MAXFEV, ACC = 1e-12, 200, 1e-6
|
|
SEED = 20260824
|
|
REAL = 40
|
|
|
|
|
|
def ours(f, a, b, xtol, maxfev):
|
|
return pure_root(f, a, b, xtol, maxfev, similarity="ratio",
|
|
scale="derived3", guard="credit")
|
|
|
|
|
|
OURS = "OURS"
|
|
FIELD = [
|
|
(OURS, ours),
|
|
("chandrupatla", CHANDRUPATLA_BASELINE),
|
|
("toms748", toms748),
|
|
("brent", BRENT_BASELINE),
|
|
("itp", itp),
|
|
("ridders", ridders),
|
|
("bisection", bisection),
|
|
]
|
|
|
|
# Counted by reading each method's definition, and counted the same way for
|
|
# ours. A constant is free if changing it leaves a working method whose
|
|
# performance differs - i.e. if it had to be picked.
|
|
CONSTANTS = {
|
|
OURS: (0, "the credit endowment is one unconstrained step, not a "
|
|
"number: every value in (0,1] promises the same bisection+1 "
|
|
"bound and 1 is the cheapest of them, while 0 locks the "
|
|
"mechanism into bisection for ever"),
|
|
"chandrupatla": (0, "the acceptance test is derived from where inverse "
|
|
"quadratic interpolation stays inside the bracket"),
|
|
"toms748": (2, "the interpolation order schedule and mu"),
|
|
"brent": (3, "three empirical acceptance inequalities"),
|
|
"itp": (3, "kappa1, kappa2, n0"),
|
|
"ridders": (0, "none"),
|
|
"bisection": (0, "none"),
|
|
}
|
|
GUARANTEE = {
|
|
OURS: "bisection + 1 step",
|
|
"chandrupatla": "bisection-bounded by construction",
|
|
"toms748": "bisection-bounded",
|
|
"brent": "quadratic in the bisection count",
|
|
"itp": "bisection + n0, optimal minmax",
|
|
"ridders": "bisection-bounded",
|
|
"bisection": "exactly the bisection bound",
|
|
}
|
|
|
|
|
|
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 run(method, probs, maxfev=MAXFEV):
|
|
per = {}
|
|
for _s, pid, f, a, b, root in probs:
|
|
try:
|
|
x, n, conv = method(f, a, b, XTOL, maxfev)
|
|
except Exception:
|
|
x, n, conv = float("nan"), maxfev, False
|
|
ok = bool(conv and math.isfinite(x) and abs(x - root) <= ACC)
|
|
per[pid] = (n if ok else maxfev, ok, n)
|
|
return per
|
|
|
|
|
|
def sums(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 sign_p(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 verdict(ours_v, other_v, noise, lower_is_better=True):
|
|
d = ours_v - other_v
|
|
if abs(d) <= noise:
|
|
return "LEVEL"
|
|
if (d < 0) == lower_is_better:
|
|
return "AHEAD"
|
|
return "BEHIND"
|
|
|
|
|
|
def main():
|
|
probs = load()
|
|
t0 = time.time()
|
|
print(f"{len(probs)} problems; competitors only; seed {SEED}\n", flush=True)
|
|
|
|
per = {n: run(m, probs) for n, m in FIELD}
|
|
tot = {n: sums(per[n], probs) for n, _ in FIELD}
|
|
|
|
# ---- E first: the noise floor, because every verdict needs it ------
|
|
print("Measuring the noise floor (+-1 ulp, "
|
|
f"{REAL} realisations)...", flush=True)
|
|
spread = {}
|
|
for name, m in FIELD:
|
|
vals = []
|
|
for k in range(REAL):
|
|
r = random.Random(SEED + k)
|
|
|
|
def perturb(f, rr=r):
|
|
seen = {}
|
|
|
|
def g(x):
|
|
if x not in seen:
|
|
v = f(x)
|
|
if v != 0.0 and math.isfinite(v):
|
|
v = math.nextafter(
|
|
v, math.inf if rr.random() < 0.5 else -math.inf)
|
|
seen[x] = v
|
|
return seen[x]
|
|
return g
|
|
|
|
pp = [(s, pid, perturb(f), a, b, root)
|
|
for s, pid, f, a, b, root in probs]
|
|
t, _ = sums(run(m, pp), pp)
|
|
vals.append(t["APS"] + t["CHA"])
|
|
spread[name] = max(vals) - min(vals)
|
|
print(f" {name:14s} spread {spread[name]:>5d} "
|
|
f"({time.time() - t0:.0f}s)", flush=True)
|
|
noise = spread[OURS]
|
|
print(f"\n our own spread is {noise} evaluations; any gap smaller than "
|
|
f"that is\n not a result, and is called LEVEL below.\n", flush=True)
|
|
|
|
ot, ook = tot[OURS]
|
|
ours_both = ot["APS"] + ot["CHA"]
|
|
|
|
print("=" * 92, flush=True)
|
|
print("THE VERDICT, AXIS BY AXIS", flush=True)
|
|
print("=" * 92, flush=True)
|
|
print(f"{'competitor':14s}{'A evals both':>14s}{'verdict':>9s}"
|
|
f"{'A APS':>8s}{'verdict':>9s}{'A CHA':>8s}{'verdict':>9s}"
|
|
f"{'B solved':>10s}{'G consts':>10s}", flush=True)
|
|
print("-" * 92, flush=True)
|
|
print(f"{OURS:14s}{ours_both:>14d}{'-':>9s}{ot['APS']:>8d}{'-':>9s}"
|
|
f"{ot['CHA']:>8d}{'-':>9s}"
|
|
f"{ook['APS'] + ook['CHA']:>7d}/199{CONSTANTS[OURS][0]:>10d}",
|
|
flush=True)
|
|
for name, _m in FIELD:
|
|
if name == OURS:
|
|
continue
|
|
t, ok = tot[name]
|
|
both = t["APS"] + t["CHA"]
|
|
print(f"{name:14s}{both:>14d}"
|
|
f"{verdict(ours_both, both, noise):>9s}"
|
|
f"{t['APS']:>8d}{verdict(ot['APS'], t['APS'], noise):>9s}"
|
|
f"{t['CHA']:>8d}{verdict(ot['CHA'], t['CHA'], noise):>9s}"
|
|
f"{ok['APS'] + ok['CHA']:>7d}/199{CONSTANTS[name][0]:>10d}",
|
|
flush=True)
|
|
|
|
print("\n" + "=" * 92, flush=True)
|
|
print("C HEAD TO HEAD, and D WORST SINGLE PROBLEM", flush=True)
|
|
print("=" * 92, flush=True)
|
|
# The worst single problem must be taken over the problems BOTH solve.
|
|
# Over all problems it is the evaluation cap for every method, so the axis
|
|
# would report the budget rather than the methods - it did, 200 against
|
|
# 200 for all six.
|
|
print(f"{'competitor':14s}{'W/T/L':>14s}{'sign p':>11s}{'verdict':>10s}"
|
|
f"{'their worst':>13s}{'our worst':>11s}{'verdict':>10s}", flush=True)
|
|
print("-" * 92, flush=True)
|
|
for name, _m in FIELD:
|
|
if name == OURS:
|
|
continue
|
|
w = t_ = l = 0
|
|
ow = tw = 0
|
|
for _s, pid, *_ in probs:
|
|
(no, oko, _), (nc, okc, _) = per[OURS][pid], per[name][pid]
|
|
if not (oko and okc):
|
|
continue
|
|
ow = max(ow, no)
|
|
tw = max(tw, nc)
|
|
if no < nc:
|
|
w += 1
|
|
elif no > nc:
|
|
l += 1
|
|
else:
|
|
t_ += 1
|
|
p = sign_p(w, l)
|
|
v = ("LEVEL" if p > 0.05 else ("AHEAD" if w > l else "BEHIND"))
|
|
print(f"{name:14s}{f'{w}/{t_}/{l}':>14s}{p:>11.2e}{v:>10s}"
|
|
f"{tw:>13d}{ow:>11d}"
|
|
f"{verdict(ow, tw, 0):>10s}", flush=True)
|
|
|
|
print("\n" + "=" * 92, flush=True)
|
|
print("E ROBUSTNESS, F GUARANTEE, G WHAT EACH METHOD HAD TO CHOOSE",
|
|
flush=True)
|
|
print("=" * 92, flush=True)
|
|
for name, _m in FIELD:
|
|
c, why = CONSTANTS[name]
|
|
print(f" {name:14s} spread {spread[name]:>5d} "
|
|
f"constants {c} guarantee: {GUARANTEE[name]}", flush=True)
|
|
print(f" {'':14s} {why}", flush=True)
|
|
|
|
print("\n" + "=" * 92, flush=True)
|
|
print("SUMMARY", flush=True)
|
|
print("=" * 92, flush=True)
|
|
ahead, level, behind = [], [], []
|
|
for name, _m in FIELD:
|
|
if name == OURS:
|
|
continue
|
|
both = tot[name][0]["APS"] + tot[name][0]["CHA"]
|
|
v = verdict(ours_both, both, noise)
|
|
(ahead if v == "AHEAD" else level if v == "LEVEL" else behind).append(name)
|
|
print(f" by total evaluations over both sets:", flush=True)
|
|
print(f" ahead of : {', '.join(ahead) or 'nobody'}", flush=True)
|
|
print(f" level with: {', '.join(level) or 'nobody'}", flush=True)
|
|
print(f" behind : {', '.join(behind) or 'nobody'}", flush=True)
|
|
print(f"\nfinished in {time.time() - t0:.0f}s", flush=True)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|