255 lines
9.5 KiB
Python
255 lines
9.5 KiB
Python
"""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 (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", best_of(chandrupatla, chandrupatla_scipy)),
|
|
("toms748 (scipy)", toms748),
|
|
("brent", best_of(brent, brent_scipy)),
|
|
("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):
|
|
"""Charged evaluation count per problem, and whether it was solved."""
|
|
out = {}
|
|
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)
|
|
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()
|