"""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()