278 lines
11 KiB
Python
278 lines
11 KiB
Python
"""Stage 16 - similarity used to IMPROVE an existing method, not to replace it.
|
|
|
|
Everything that worked in this project worked that way. The credit safeguard
|
|
did not produce a new root finder that beats the field; it took regula falsi,
|
|
the secant rule and inverse quadratic interpolation and made each of them
|
|
solve 193 of 199 with a near-bisection worst case. Building a method from
|
|
similarity landed second behind Chandrupatla 1997; grafting a measured
|
|
quantity onto existing methods was the part that paid.
|
|
|
|
This stage tries the one graft that comes straight from the similarity
|
|
literature and has not been tried here: **similarity replaces a hard
|
|
partition by a graded weight**.
|
|
|
|
Every strong bracketing method contains such a partition:
|
|
|
|
Brent three empirical inequalities decide interpolation OR bisection
|
|
Chandrupatla the criterion 1-sqrt(1-xi) < phi < sqrt(xi) decides the same
|
|
|
|
In both, the quantity being tested is continuous and the decision taken on it
|
|
is a step function. The similarity move is to keep the quantity and drop the
|
|
step:
|
|
|
|
r = (how far outside the safe region the interpolation falls), 0 = safe
|
|
w = S(r), S the tent similarity, S(0) = 1, no constant
|
|
d = w * d_interpolation + (1 - w) * d_bisection
|
|
|
|
At r = 0 this is exactly the original method. At r >= 1 it is exactly
|
|
bisection. In between, where the original flips discontinuously, it blends.
|
|
|
|
The comparison is against the SAME implementation with the blend switched
|
|
off, so the difference can only be the blend. Our Brent is separately
|
|
validated against scipy.optimize.brentq to 0.2% of evaluations.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
|
|
from scipy.optimize._tstutils import (_APS_TESTS_DICTS,
|
|
_CHANDRUPATLA_TESTS_DICTS)
|
|
|
|
from simroot.baselines import Counter, _Budget, brent, chandrupatla
|
|
|
|
XTOL, MAXFEV, ACC = 1e-12, 200, 1e-6
|
|
|
|
# the tent similarity of the paper: S(0) = 1, S(1) = 0, no constant
|
|
TENT = lambda r: max(0.0, 1.0 - r)
|
|
|
|
|
|
def brent_blend(f, a, b, xtol, maxfev, mode="blend"):
|
|
"""Brent's method with its accept/reject replaced by a similarity blend.
|
|
|
|
The only change from simroot.baselines.brent is the four lines marked
|
|
below; with mode="off" the two are the same algorithm. mode is a name, not a
|
|
boolean carrying three states."""
|
|
g = Counter(f, maxfev)
|
|
b_ = b
|
|
try:
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
c, fc = a, fa
|
|
d = e = b - a
|
|
while True:
|
|
if fb * fc > 0.0:
|
|
c, fc = a, fa
|
|
d = e = b - a
|
|
if abs(fc) < abs(fb):
|
|
a, fa = b, fb
|
|
b, fb = c, fc
|
|
c, fc = a, fa
|
|
tol1 = 2.0 * 2.22e-16 * abs(b) + 0.5 * xtol
|
|
xm = 0.5 * (c - b)
|
|
if abs(xm) <= tol1 or fb == 0.0:
|
|
return b, g.n, True
|
|
if abs(e) >= tol1 and abs(fa) > abs(fb):
|
|
s = fb / fa
|
|
if a == c:
|
|
p = 2.0 * xm * s
|
|
q = 1.0 - s
|
|
else:
|
|
q = fa / fc
|
|
r = fb / fc
|
|
p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0))
|
|
q = (q - 1.0) * (r - 1.0) * (s - 1.0)
|
|
if p > 0.0:
|
|
q = -q
|
|
p = abs(p)
|
|
bound = min(3.0 * xm * q - abs(tol1 * q), abs(e * q))
|
|
if mode == "project":
|
|
# The lesson from the root-finding safeguard, applied to
|
|
# somebody else's method: when the interpolation violates
|
|
# the safety condition, do not discard it and do not
|
|
# average it with the midpoint - keep as much of it as the
|
|
# condition allows, by clamping it to the boundary.
|
|
if 2.0 * p < bound:
|
|
e, d = d, p / q
|
|
elif bound > 0.0 and q != 0.0:
|
|
e, d = d, math.copysign(0.5 * bound / abs(q),
|
|
p / q if q else xm)
|
|
else:
|
|
d = e = xm
|
|
elif mode == "off":
|
|
if 2.0 * p < bound:
|
|
e, d = d, p / q
|
|
else:
|
|
d = e = xm
|
|
else:
|
|
# ---- the graft -------------------------------------
|
|
# Brent's own test is 2p < bound. The ratio of the two
|
|
# is continuous; the decision on it is a step. Keep the
|
|
# ratio, weight by similarity, blend the two candidate
|
|
# steps. r = 0 reproduces Brent exactly, r >= 1 gives
|
|
# bisection exactly.
|
|
r_trust = (2.0 * p / bound) if bound > 0.0 else math.inf
|
|
w = TENT(r_trust) if math.isfinite(r_trust) else 0.0
|
|
d_int = p / q if q != 0.0 else xm
|
|
e, d = d, w * d_int + (1.0 - w) * xm
|
|
# ----------------------------------------------------
|
|
else:
|
|
d = e = xm
|
|
a, fa = b, fb
|
|
b = b + (d if abs(d) > tol1 else math.copysign(tol1, xm))
|
|
b_ = b
|
|
fb = g(b)
|
|
except _Budget:
|
|
return b_, g.n, False
|
|
|
|
|
|
def chandrupatla_blend(f, a, b, xtol, maxfev, mode="blend"):
|
|
"""Chandrupatla's method with its acceptance region softened the same way.
|
|
|
|
Its test is 1 - sqrt(1-xi) < phi < sqrt(xi). The distance outside that
|
|
interval, scaled by its own width, is the continuous quantity; the
|
|
decision on it is a step."""
|
|
g = Counter(f, maxfev)
|
|
x1, x2 = b, a
|
|
f1 = f2 = 0.0
|
|
try:
|
|
f1, f2 = g(x1), g(x2)
|
|
if f1 == 0.0:
|
|
return x1, g.n, True
|
|
if f2 == 0.0:
|
|
return x2, g.n, True
|
|
x3, f3 = x2, f2
|
|
t = 0.5
|
|
while True:
|
|
x = x1 + t * (x2 - x1)
|
|
fx = g(x)
|
|
if (fx > 0.0) == (f1 > 0.0):
|
|
x3, f3 = x1, f1
|
|
x1, f1 = x, fx
|
|
else:
|
|
x3, f3 = x2, f2
|
|
x2, f2 = x1, f1
|
|
x1, f1 = x, fx
|
|
xm, fm = (x1, f1) if abs(f1) < abs(f2) else (x2, f2)
|
|
if fm == 0.0:
|
|
return xm, g.n, True
|
|
tol = 2.0 * 2.22e-16 * abs(xm) + xtol
|
|
span = abs(x2 - x1)
|
|
if span == 0.0:
|
|
return xm, g.n, True
|
|
tl = tol / span
|
|
if tl > 0.5:
|
|
return xm, g.n, True
|
|
t = 0.5
|
|
if x3 != x2 and f3 != f2:
|
|
xi = (x1 - x2) / (x3 - x2)
|
|
ph = (f1 - f2) / (f3 - f2)
|
|
if 0.0 < xi < 1.0:
|
|
lo = 1.0 - math.sqrt(1.0 - xi)
|
|
hi = math.sqrt(xi)
|
|
d1, d2, d3 = f1 - f2, f1 - f3, f2 - f3
|
|
inside = lo < ph < hi
|
|
if d1 != 0.0 and d2 != 0.0 and d3 != 0.0:
|
|
t_iqi = ((f1 / d1) * (f3 / (-d3))
|
|
+ ((x3 - x1) / (x2 - x1)) * (f1 / d2)
|
|
* (f2 / d3))
|
|
if mode == "off":
|
|
if inside:
|
|
t = t_iqi
|
|
else:
|
|
# distance outside the safe interval, in units of
|
|
# the interval's own width
|
|
width = hi - lo
|
|
if width > 0.0:
|
|
out = max(0.0, lo - ph, ph - hi) / width
|
|
w = TENT(out)
|
|
t = w * t_iqi + (1.0 - w) * 0.5
|
|
elif inside:
|
|
t = t_iqi
|
|
if not math.isfinite(t):
|
|
t = 0.5
|
|
t = min(max(t, tl), 1.0 - tl)
|
|
except _Budget:
|
|
return (x1 if abs(f1) < abs(f2) else x2), g.n, False
|
|
|
|
|
|
def load():
|
|
out = []
|
|
for d in _APS_TESTS_DICTS:
|
|
a, b = d["bracket"]
|
|
args = d.get("args", ())
|
|
out.append(("APS", (lambda ff, aa: lambda t: ff(t, *aa))(d["f"], args),
|
|
float(a), float(b), float(d["root"])))
|
|
for d in _CHANDRUPATLA_TESTS_DICTS:
|
|
a, b = d["bracket"]
|
|
out.append(("CHA", (lambda ff: lambda t: ff(t))(d["f"]),
|
|
float(a), float(b), float(d["root"])))
|
|
return out
|
|
|
|
|
|
def measure(method, probs):
|
|
t = {"APS": 0, "CHA": 0}
|
|
ok = {"APS": 0, "CHA": 0}
|
|
worst = 0
|
|
for s, 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
|
|
good = bool(conv and math.isfinite(x) and abs(x - root) <= ACC)
|
|
t[s] += n if good else MAXFEV
|
|
ok[s] += good
|
|
if good:
|
|
worst = max(worst, n)
|
|
return t, ok, worst
|
|
|
|
|
|
def main():
|
|
probs = load()
|
|
rows = [
|
|
("brent, as published", lambda f, a, b, xt, mx: brent(f, a, b, xt, mx)),
|
|
("brent, blend off (same code)",
|
|
lambda f, a, b, xt, mx: brent_blend(f, a, b, xt, mx, mode="off")),
|
|
("brent + similarity blend",
|
|
lambda f, a, b, xt, mx: brent_blend(f, a, b, xt, mx, mode="blend")),
|
|
("brent + projection onto the bound",
|
|
lambda f, a, b, xt, mx: brent_blend(f, a, b, xt, mx,
|
|
mode="project")),
|
|
("chandrupatla, as published",
|
|
lambda f, a, b, xt, mx: chandrupatla(f, a, b, xt, mx)),
|
|
("chandrupatla, blend off",
|
|
lambda f, a, b, xt, mx: chandrupatla_blend(f, a, b, xt, mx,
|
|
mode="off")),
|
|
("chandrupatla + similarity blend",
|
|
lambda f, a, b, xt, mx: chandrupatla_blend(f, a, b, xt, mx,
|
|
mode="blend")),
|
|
]
|
|
print(f"{'method':34s}{'APS':>8s}{'CHA':>8s}{'both':>8s}"
|
|
f"{'solved':>10s}{'worst':>8s}")
|
|
print("-" * 78)
|
|
base = {}
|
|
for name, m in rows:
|
|
t, ok, worst = measure(m, probs)
|
|
both = t["APS"] + t["CHA"]
|
|
base[name] = both
|
|
print(f"{name:34s}{t['APS']:>8d}{t['CHA']:>8d}{both:>8d}"
|
|
f"{ok['APS'] + ok['CHA']:>7d}/199{worst:>8d}", flush=True)
|
|
|
|
print("\nThe two 'blend off' rows must reproduce the published methods -"
|
|
"\nif they do not, the blend is being compared against a different"
|
|
"\nalgorithm and the experiment is void.")
|
|
for pub, off, on in (("brent, as published", "brent, blend off (same code)",
|
|
"brent + similarity blend"),
|
|
("chandrupatla, as published", "chandrupatla, blend off",
|
|
"chandrupatla + similarity blend")):
|
|
drift = abs(base[off] - base[pub]) / base[pub]
|
|
gain = (base[pub] - base[on]) / base[pub]
|
|
print(f"\n {pub.split(',')[0]:14s} control drift {drift * 100:5.2f}%"
|
|
f" blend changes the total by {gain * 100:+6.2f}%")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|