334 lines
13 KiB
Python
334 lines
13 KiB
Python
"""Stage 9 - wall-clock time, and what it is worth.
|
|
|
|
Evaluations of f are the right cost model for a root finder because in any
|
|
real problem f is the expensive part. Wall-clock time measures something
|
|
different and narrower: the method's OWN arithmetic per step. It matters only
|
|
when f is cheap, and it is implementation- and language-bound, so a timing
|
|
comparison is evidence about this code rather than about the algorithms.
|
|
|
|
Three measurements, in increasing order of usefulness to the paper:
|
|
|
|
1. TOTAL TIME over the 199 benchmark problems, whose f are all cheap. This
|
|
is the pessimal case for a method that does more arithmetic per step.
|
|
2. OVERHEAD PER EVALUATION, i.e. total time divided by evaluations, with the
|
|
cost of f itself subtracted out. This is the machine-comparable number.
|
|
3. BREAK-EVEN COST OF f: how expensive f must be before the difference in
|
|
evaluations outweighs the difference in overhead. This is the number a
|
|
practitioner actually needs, and it is the one worth publishing.
|
|
|
|
Timing method: each configuration runs the whole set R times and the MINIMUM
|
|
is kept. The minimum is the right estimator here - noise from other processes
|
|
can only ever add time, never remove it - and it is not a summary statistic
|
|
standing in for a decision.
|
|
|
|
Only pure-Python implementations are compared, so the arithmetic is
|
|
commensurable. scipy's toms748 is pure Python too and is included; a compiled
|
|
brentq would not be comparable and is excluded.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
import time
|
|
|
|
from scipy.optimize._tstutils import (_APS_TESTS_DICTS,
|
|
_CHANDRUPATLA_TESTS_DICTS)
|
|
|
|
from simroot.baselines import bisection, brent, itp, ridders, toms748
|
|
from simroot.pure import pure_root
|
|
|
|
XTOL = 1e-12
|
|
MAXFEV = 200
|
|
ACC = 1e-6
|
|
REPEATS = 5
|
|
|
|
|
|
def ours(f, a, b, xtol, maxfev):
|
|
return pure_root(f, a, b, xtol, maxfev, similarity="ratio",
|
|
scale="derived3", guard="project")
|
|
|
|
|
|
def ours_inlined(f, a, b, xtol, maxfev):
|
|
"""The SAME algorithm as pure_root(ratio, derived3, project), written as
|
|
one flat loop instead of a similarity callable plus three helper
|
|
functions.
|
|
|
|
This exists only to answer one question: is our per-step cost algorithmic,
|
|
or is it Python call overhead in a modular implementation? It must return
|
|
exactly the same iterates and the same evaluation count as pure_root - the
|
|
caller checks that - so any time it saves is overhead, not arithmetic.
|
|
"""
|
|
n = 0
|
|
x = 0.5 * (a + b)
|
|
try:
|
|
fa = f(a); fb = f(b); n = 2
|
|
if fa == 0.0:
|
|
return a, n, True
|
|
if fb == 0.0:
|
|
return b, n, True
|
|
if fa * fb > 0.0:
|
|
raise ValueError("f(a) and f(b) must straddle zero")
|
|
A, FA, B = a, fa, b
|
|
if abs(fa) <= abs(fb):
|
|
x, fx, xp, fp = a, fa, b, fb
|
|
else:
|
|
x, fx, xp, fp = b, fb, a, fa
|
|
xp2 = fp2 = None
|
|
eps = 0.5 * xtol
|
|
w0 = abs(B - A)
|
|
n_max = 2 * (math.ceil(math.log2(w0 / xtol)) if w0 > xtol else 0)
|
|
allow = 2.0 ** n_max * eps
|
|
s2 = s1 = w0
|
|
while abs(B - A) > xtol:
|
|
lo = A if A < B else B
|
|
hi = B if A < B else A
|
|
far = B if abs(B - x) >= abs(A - x) else A
|
|
# three-point divided difference, inline
|
|
if x == xp:
|
|
slope = 0.0
|
|
else:
|
|
d12 = (fx - fp) / (x - xp)
|
|
if not math.isfinite(d12):
|
|
slope = 0.0 # as _slope3 does: no usable slope
|
|
elif xp2 is None or xp2 == xp or xp2 == x:
|
|
slope = d12
|
|
else:
|
|
d01 = (fp - fp2) / (xp - xp2)
|
|
slope = d12 + (d12 - d01) / (x - xp2) * (x - xp)
|
|
if not math.isfinite(slope):
|
|
slope = d12
|
|
sc = abs(slope) * abs(far - x)
|
|
if sc == 0.0 or not math.isfinite(sc):
|
|
xn = 0.5 * (A + B)
|
|
else:
|
|
z = abs(fx) / sc
|
|
step = (z if z < 1.0 else 1.0) * (far - x) # 1 - max(0,1-z)
|
|
# NB the proposed configuration is guard="project", which has
|
|
# NO step test - applying one here made this a different
|
|
# algorithm on 132 of the 199 problems, and the equality
|
|
# check above caught it.
|
|
xn = x + step
|
|
if allow == allow and allow != math.inf: # finite
|
|
left = hi - allow
|
|
right = lo + allow
|
|
if left > right:
|
|
xn = 0.5 * (lo + hi)
|
|
elif xn < left:
|
|
xn = left
|
|
elif xn > right:
|
|
xn = right
|
|
step = xn - x
|
|
s2 = s1
|
|
s1 = abs(step)
|
|
if not math.isfinite(xn) or not (lo < xn < hi):
|
|
xn = 0.5 * (A + B)
|
|
if xn == x or xn == xp:
|
|
xn = 0.5 * (A + B)
|
|
if xn == x or xn == xp:
|
|
return xn, n, abs(B - A) <= xtol
|
|
if n >= maxfev:
|
|
return x, n, False
|
|
fn = f(xn); n += 1
|
|
allow *= 0.5
|
|
xp2, fp2, xp, fp, x, fx = xp, fp, x, fx, xn, fn
|
|
if fn == 0.0:
|
|
return xn, n, True
|
|
if FA * fn < 0.0:
|
|
B = xn
|
|
else:
|
|
A, FA = xn, fn
|
|
return 0.5 * (A + B), n, True
|
|
except ValueError:
|
|
raise
|
|
except Exception:
|
|
return x, n, False
|
|
|
|
|
|
METHODS = [
|
|
("OURS (similarity)", ours),
|
|
("OURS, inlined", ours_inlined),
|
|
("brent", brent),
|
|
("toms748 (scipy)", toms748),
|
|
("itp", itp),
|
|
("ridders", ridders),
|
|
("bisection", bisection),
|
|
]
|
|
|
|
|
|
def load():
|
|
out = []
|
|
for d in _APS_TESTS_DICTS:
|
|
a, b = d["bracket"]
|
|
args = d.get("args", ())
|
|
out.append(((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(((lambda ff: lambda x: ff(x))(d["f"]),
|
|
float(a), float(b), float(d["root"])))
|
|
return out
|
|
|
|
|
|
def sweep_once(method, probs):
|
|
"""One pass over every problem; returns (seconds, evaluations, solved)."""
|
|
t0 = time.perf_counter()
|
|
evals = solved = 0
|
|
for 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
|
|
evals += n
|
|
solved += bool(conv and math.isfinite(x) and abs(x - root) <= ACC)
|
|
return time.perf_counter() - t0, evals, solved
|
|
|
|
|
|
def make_costly(f, work):
|
|
"""f with a controlled amount of extra arithmetic per call, standing in
|
|
for an f that is expensive because it runs a simulation or reads a table."""
|
|
if work <= 0:
|
|
return f
|
|
|
|
def g(x):
|
|
s = 0.0
|
|
for i in range(work):
|
|
s += (i + x) * 1.0000001
|
|
return f(x) + s * 0.0
|
|
|
|
return g
|
|
|
|
|
|
def calibrate(probs, work):
|
|
"""The wall-clock cost of one call at this work level, measured the same
|
|
way the solvers will pay it."""
|
|
# calibrating on one problem's f measured that f's quirks, not the work
|
|
# level: the result came out non-monotonic in the work, which is proof
|
|
# the estimate was noise. Sample across the set instead.
|
|
sample = probs[::10]
|
|
pairs = []
|
|
for f, a, b, _r in sample:
|
|
g = make_costly(f, work)
|
|
for k in range(10):
|
|
pairs.append((g, a + (b - a) * (k + 0.5) / 10.0))
|
|
best = math.inf
|
|
for _ in range(REPEATS):
|
|
t0 = time.perf_counter()
|
|
for g, p in pairs:
|
|
g(p)
|
|
best = min(best, (time.perf_counter() - t0) / len(pairs))
|
|
return best
|
|
|
|
|
|
def main():
|
|
probs = load()
|
|
print(f"{len(probs)} problems, {REPEATS} repeats, minimum kept\n")
|
|
|
|
|
|
# the inlined variant is only admissible evidence if it is the same
|
|
# algorithm; anything else and its timing means nothing
|
|
mism = []
|
|
for f, a, b, root in probs:
|
|
try:
|
|
r1 = ours(f, a, b, XTOL, MAXFEV)
|
|
except Exception:
|
|
r1 = None
|
|
try:
|
|
r2 = ours_inlined(f, a, b, XTOL, MAXFEV)
|
|
except Exception:
|
|
r2 = None
|
|
if (r1 is None) != (r2 is None) or (r1 and r2 and r1[1] != r2[1]):
|
|
mism.append((r1, r2))
|
|
print(f"inlined variant vs pure_root: {len(mism)} of {len(probs)} "
|
|
f"problems where the evaluation count differs")
|
|
methods = list(METHODS)
|
|
if mism:
|
|
print(" NOT the same algorithm, so it is DROPPED from the tables")
|
|
print(" rather than shown with a caveat nobody reads. The question")
|
|
print(" it was meant to settle - whether our per-step cost is")
|
|
print(" arithmetic or Python call overhead - therefore stays open.")
|
|
print(f" first divergence: {mism[0]}")
|
|
methods = [(n, m) for n, m in methods if "inlined" not in n]
|
|
print()
|
|
|
|
rows = []
|
|
for name, m in methods:
|
|
secs, evals, solved = min(sweep_once(m, probs) for _ in range(REPEATS))
|
|
rows.append((name, secs, evals, solved))
|
|
|
|
print("1. TOTAL TIME on the benchmark as it stands - every f here is")
|
|
print(" cheap, which is the worst case for a method that does more")
|
|
print(" arithmetic per step.")
|
|
print(f"{'method':22s}{'ms':>10s}{'evals':>9s}{'solved':>9s}"
|
|
f"{'vs ours':>10s}")
|
|
print("-" * 62)
|
|
base = next(r[1] for r in rows if r[0].startswith("OURS"))
|
|
for name, secs, evals, solved in rows:
|
|
print(f"{name:22s}{secs * 1e3:>10.1f}{evals:>9d}{solved:>7d}/199"
|
|
f"{(secs / base - 1) * 100:>9.1f}%")
|
|
|
|
# ------------------------------------------------------------------
|
|
# 2. The crossover, measured rather than modelled.
|
|
#
|
|
# An earlier version estimated a per-evaluation overhead by measuring
|
|
# the cost of f separately and subtracting evals * t_f. It produced
|
|
# NEGATIVE overheads for four of six methods, which is impossible and
|
|
# proved the estimate wrong: sampling f across the whole bracket
|
|
# charges it for expensive regions the solvers never visit.
|
|
#
|
|
# Nothing is subtracted here. f is made progressively more expensive
|
|
# by a controlled amount of extra arithmetic, and the totals are timed
|
|
# directly. The crossover is then read off, not derived.
|
|
# ------------------------------------------------------------------
|
|
print("\n2. WHERE THE CROSSOVER IS (f made progressively more expensive,")
|
|
print(" nothing subtracted; total ms over all 199 problems)")
|
|
works = [0, 3, 10, 30, 100, 300]
|
|
costs = [calibrate(probs, w) for w in works]
|
|
if any(costs[i + 1] <= costs[i] for i in range(len(costs) - 1)):
|
|
print(" WARNING: the calibration is not monotonic in the work "
|
|
"level, so the\n machine is too noisy for this measurement "
|
|
"and the crossover below\n must not be quoted.")
|
|
head = (f"{'cost of one f call':22s}"
|
|
+ "".join(f"{c * 1e6:>11.2f}us" for c in costs))
|
|
print(head)
|
|
print("-" * len(head))
|
|
trend = {}
|
|
for name, m in methods:
|
|
cells = []
|
|
for w in works:
|
|
pr = [(make_costly(f, w), a, b, r) for f, a, b, r in probs]
|
|
secs = min(sweep_once(m, pr)[0] for _ in range(2))
|
|
cells.append(secs)
|
|
trend[name] = cells
|
|
print(f"{name:22s}" + "".join(f"{c * 1e3:>13.1f}" for c in cells))
|
|
|
|
print("\n crossover against ours:")
|
|
ours_t = trend[next(n for n, _ in methods if n.startswith('OURS'))]
|
|
for name, _m in methods:
|
|
if name.startswith("OURS"):
|
|
continue
|
|
# A crossing must PERSIST: taking the first one reported that we
|
|
# overtake Brent at 2.3 us, which the very next measurement
|
|
# contradicted. Require us to be ahead at that level and at every
|
|
# level above it, and say plainly when no such point exists.
|
|
cross = None
|
|
for i, c in enumerate(costs):
|
|
if all(trend[name][k] > ours_t[k] for k in range(i, len(costs))):
|
|
cross = c
|
|
break
|
|
flips = sum(1 for k in range(len(costs) - 1)
|
|
if (trend[name][k] > ours_t[k]) !=
|
|
(trend[name][k + 1] > ours_t[k + 1]))
|
|
note = f" [{flips} sign changes - noisy]" if flips > 1 else ""
|
|
if cross is None:
|
|
print(f" {name:20s} faster than us at every level measured, "
|
|
f"up to f = {costs[-1] * 1e6:.1f} us{note}")
|
|
elif cross == costs[0]:
|
|
print(f" {name:20s} we are faster at every level, including "
|
|
f"the raw benchmark{note}")
|
|
else:
|
|
print(f" {name:20s} we are faster from f = "
|
|
f"{cross * 1e6:.1f} us per call upwards{note}")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|