202 lines
8.4 KiB
Python
202 lines
8.4 KiB
Python
"""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()
|