"""Stage 4 - WHERE the similarity method wins, and WHY. The totals say the parameter-free similarity rule loses to Brent and TOMS748 on Alefeld-Potra-Shi and beats them on Chandrupatla. A total does not explain anything, so this file takes the runs apart. Every root search has two phases with completely different economics: APPROACH getting from the initial bracket into the neighbourhood of the root. Cost here is governed by how boldly a method is willing to jump, and by how often its safeguard drags it back to the midpoint. Bracket width drives this phase. ENDGAME going from the neighbourhood down to the tolerance. Cost here is governed by the ORDER of convergence and nothing else: secant 1.618, Brent's inverse quadratic 1.839, TOMS748's cubic up to 2.7. Prediction, if the method is what the derivation says it is: - its step is |f|/|f'| with a secant slope, so its order is ~1.618 and it must LOSE the endgame to both Brent and TOMS748, by a roughly constant number of evaluations per problem; - its approach carries no cautious truncation, so it must WIN the approach, by an amount that GROWS with the width of the bracket. - therefore: narrow brackets -> the endgame dominates -> we lose; wide brackets -> the approach dominates -> we win. All three are measured below. If the width correlation is absent, the explanation is wrong and must not go into the paper. """ from __future__ import annotations import math from scipy.optimize._tstutils import (_APS_TESTS_DICTS, _CHANDRUPATLA_TESTS_DICTS) from simroot.baselines import bisection, brent, toms748 from simroot.pure import pure_root XTOL = 1e-12 MAXFEV = 200 ACC = 1e-6 EPS = 2.220446049250313e-16 def phase_boundary(width, scale=0.0, xtol=XTOL): """Where the approach ends and the endgame begins. A fixed relative distance does not work: at 1e-4 of a bracket of width 1e12 the boundary sits at 1e8, which the very first evaluation clears, so every wide problem would report its entire cost as endgame. The boundary is the geometric mean of the initial width and the tolerance that is actually ATTAINABLE. The attainable one is not xtol: near abscissae of magnitude 1e12 the spacing of floating point is about 2e-4, so no method can drive the bracket to 1e-12 there and the journey is shorter than xtol suggests. Using xtol regardless put the boundary too low and charged the whole tail to the approach - measured on a width of 1e12, bisection split 40 / 15 instead of evenly. `scale` is the largest |abscissa| in the bracket; pass it whenever the bracket is far from the origin.""" attainable = max(xtol, EPS * abs(scale)) return math.sqrt(width * attainable) def ours(f, a, b, xtol, maxfev): return pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived", guard="step") METHODS = [("OURS (similarity)", ours), ("brent", brent), ("toms748", toms748), ("bisection", bisection)] def traced(f): """Wrap f so the sequence of abscissae it is asked about is recorded.""" xs = [] def g(x): xs.append(x) return f(x) return g, xs def phases(method, f, a, b, root): """Split a run into approach and endgame evaluations.""" g, xs = traced(f) try: x, n, conv = method(g, a, b, XTOL, MAXFEV) except Exception: return None if not conv or abs(x - root) > ACC: return None target = phase_boundary(abs(b - a), max(abs(a), abs(b))) approach = len(xs) for i, xi in enumerate(xs): if abs(xi - root) <= target: approach = i + 1 break return approach, len(xs) - approach, len(xs) def observed_order(method, f, a, b, root): """Order p from the tail of the running-best error sequence: e_{n+1} ~ C e_n^p => p = log(e_{n+1}/e_n) / log(e_n/e_{n-1}).""" g, xs = traced(f) try: method(g, a, b, XTOL, MAXFEV) except Exception: return float("nan") best, errs = float("inf"), [] for xi in xs: e = abs(xi - root) if e < best: best = e errs.append(e) errs = [e for e in errs if e > 1e-14] if len(errs) < 4: return float("nan") ps = [] for i in range(1, len(errs) - 1): d1, d2 = errs[i] / errs[i - 1], errs[i + 1] / errs[i] if 0 < d1 < 1 and 0 < d2 < 1: ps.append(math.log(d2) / math.log(d1)) if not ps: return float("nan") return sorted(ps)[len(ps) // 2] def replay_bracket(f, a, b, xs): """Reconstruct the bracket a bracketing method maintained, from the sequence of abscissae it probed, and score every probe against it. Returns (positions, shrinks, skipped) where positions[i] where probe i fell inside the bracket, 0..1 (0.5 is the midpoint, i.e. a bisection step) shrinks[i] the fraction of the bracket that probe i removed skipped probes that fell outside the proven bracket and so cannot be scored against it. Counting them matters: dropping them silently would flatter any method that probes outside, and scipy's toms748 does. The reconstruction is exact for a bracketing method because the bracket is determined by the signs of f, which are recomputed here. """ lo, hi = (a, b) if a < b else (b, a) flo = f(lo) positions, shrinks, skipped = [], [], 0 for xi in xs: if not (lo < xi < hi): skipped += 1 continue w0 = hi - lo positions.append((xi - lo) / w0) fi = f(xi) if flo * fi < 0.0: hi = xi else: lo, flo = xi, fi shrinks.append(1.0 - (hi - lo) / w0) return positions, shrinks, skipped 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 main(): probs = load() # ------------------------------------------------------------------ # 1. approach vs endgame, summed over the problems every method solves # ------------------------------------------------------------------ print("=" * 78) print("1. WHERE THE EVALUATIONS GO (only problems all four methods solve)") print("=" * 78) for setname in ("APS", "CHA"): rows = {} common = [] for s, pid, f, a, b, root in probs: if s != setname: continue r = {nm: phases(m, f, a, b, root) for nm, m in METHODS} if all(v is not None for v in r.values()): common.append((pid, abs(b - a), r)) for nm, _m in METHODS: app = sum(c[2][nm][0] for c in common) end = sum(c[2][nm][1] for c in common) rows[nm] = (app, end, app + end) print(f"\n{setname}: {len(common)} problems solved by all four") print(f"{'method':22s}{'approach':>11s}{'endgame':>10s}{'total':>9s}" f"{'approach %':>12s}") print("-" * 64) for nm, _m in METHODS: app, end, tot = rows[nm] print(f"{nm:22s}{app:>11d}{end:>10d}{tot:>9d}" f"{100.0 * app / tot:>11.1f}%") # ------------------------------------------------------------------ # 2. observed order of convergence # ------------------------------------------------------------------ print("\n" + "=" * 78) print("2. OBSERVED ORDER OF CONVERGENCE (median over solved problems)") print("=" * 78) print(f"{'method':22s}{'median order':>14s}{'problems':>11s} theory") print("-" * 66) theory = {"OURS (similarity)": "1.618 (secant step |f|/|f'|)", "brent": "1.839 (inverse quadratic)", "toms748": "up to 2.7 (cubic)", "bisection": "1.0 (linear, rate 1/2)"} for nm, m in METHODS: ps = [observed_order(m, f, a, b, root) for _s, _pid, f, a, b, root in probs] ps = [p for p in ps if math.isfinite(p) and 0.5 < p < 5.0] med = sorted(ps)[len(ps) // 2] if ps else float("nan") print(f"{nm:22s}{med:>14.3f}{len(ps):>11d} {theory[nm]}") # ------------------------------------------------------------------ # 3. the prediction: our margin over brent grows with bracket width # ------------------------------------------------------------------ print("\n" + "=" * 78) print("3. MARGIN OVER BRENT BY BRACKET WIDTH (negative = we spend fewer)") print("=" * 78) buckets = [(0, 1e1), (1e1, 1e3), (1e3, 1e6), (1e6, 1e15)] labels = ["width < 10", "10 - 1e3", "1e3 - 1e6", "> 1e6"] print(f"{'bracket width':16s}{'problems':>10s}{'OURS':>9s}{'brent':>8s}" f"{'diff':>8s}{'diff/problem':>14s}") print("-" * 66) for (lo, hi), lab in zip(buckets, labels): n_ours = n_brent = cnt = 0 for _s, _pid, f, a, b, root in probs: w = abs(b - a) if not (lo <= w < hi): continue ro = phases(ours, f, a, b, root) rb = phases(brent, f, a, b, root) if ro is None or rb is None: continue cnt += 1 n_ours += ro[2] n_brent += rb[2] if cnt: d = n_ours - n_brent print(f"{lab:16s}{cnt:>10d}{n_ours:>9d}{n_brent:>8d}{d:>8d}" f"{d / cnt:>14.2f}") print("\nThe trend is not monotone, so bracket width is a correlate and " "not the cause.\nSection 4 tests the mechanism directly.") # ------------------------------------------------------------------ # 4. the mechanism: how often does each method's safeguard fire? # An interpolation step that gets thrown away and replaced by a # midpoint is a wasted evaluation. If our advantage is real, our # safeguard must fire LESS often than Brent's exactly on the set # where we win, and not much less where we lose. # ------------------------------------------------------------------ print("\n" + "=" * 78) print("4. HOW OFTEN THE SAFEGUARD THROWS AN INTERPOLATION AWAY") print("=" * 78) print(f"{'set':6s}{'method':22s}{'steps':>9s}{'forced to midpoint':>21s}" f"{'share':>9s}") print("-" * 68) for setname in ("APS", "CHA"): for nm in ("OURS (similarity)", "brent"): steps = forced = 0 for s, _pid, f, a, b, root in probs: if s != setname: continue st = {} if nm == "brent": x, n, conv = brent(f, a, b, XTOL, MAXFEV, stats=st) fired = st.get("bisect", 0) else: x, n, conv = pure_root(f, a, b, XTOL, MAXFEV, similarity="ratio", scale="derived", guard="step", stats=st) fired = st.get("forced", 0) if not conv or abs(x - root) > ACC: continue steps += n - 2 # the two endpoint evaluations are free forced += fired share = 100.0 * forced / steps if steps else float("nan") print(f"{setname:6s}{nm:22s}{steps:>9d}{forced:>21d}" f"{share:>8.1f}%") print("\nREFUTED: Brent throws away 53% of its interpolations on APS, " "where it WINS,\nand only 34% on Chandrupatla, where it LOSES. " "The safeguard rate is not the\nmechanism. Section 5 tests the " "next candidate.") # ------------------------------------------------------------------ # 5. how BOLD are the steps? # # Brent accepts an interpolation only if 2|p| < 3*xm*q, i.e. only if # it lands within three quarters of the way to the far end; otherwise # the midpoint is used. The similarity rule has no such cap - its step # is |f|/|f'| and may cross almost the whole bracket in one go. On a # wide bracket whose root sits far from the middle, that cap is the # difference between arriving in one step and bisecting towards it. # # Measured as the position t of each probe inside the bracket it was # proposed in, t = (x - lo)/(hi - lo). t near 0 or 1 is a bold step; # t = 0.5 is bisection. The bracket is reconstructed from the sign of # f at each probe, which is exact for a bracketing method. # ------------------------------------------------------------------ print("\n" + "=" * 78) print("5. HOW BOLD IS EACH STEP (position inside the bracket it was " "proposed in)") print("=" * 78) print(f"{'set':6s}{'method':22s}{'probes':>9s}{'bold (t<.1 or t>.9)':>22s}" f"{'share':>9s}") print("-" * 70) for setname in ("APS", "CHA"): for nm, m in (("OURS (similarity)", ours), ("brent", brent)): probes = bold = 0 for s, _pid, f, a, b, root in probs: if s != setname: continue g, xs = traced(f) try: x, _n, conv = m(g, a, b, XTOL, MAXFEV) except Exception: continue if not conv or abs(x - root) > ACC: continue positions, _shrinks, _skipped = replay_bracket(f, a, b, xs[2:]) probes += len(positions) bold += sum(1 for t in positions if t < 0.1 or t > 0.9) share = 100.0 * bold / probes if probes else float("nan") print(f"{setname:6s}{nm:22s}{probes:>9d}{bold:>22d}{share:>8.1f}%") print("\nREFUTED: on CHA Brent is BOLDER than we are (56.3% vs 41.8%) and " "still loses.\nBoldness is not the mechanism either. Section 6 " "measures step QUALITY rather\nthan step size.") # ------------------------------------------------------------------ # 6. step quality: how much of the bracket does one probe remove? # # Boldness is how far a probe goes; quality is whether it went the # right way. A probe at t = 0.95 that lands on the wrong side of the # root removes 5% of the bracket; the same probe on the right side # removes 95%. This is the one quantity that must differ if a method # with a LOWER order of convergence nevertheless needs fewer probes. # ------------------------------------------------------------------ print("\n" + "=" * 78) print("6. STEP QUALITY (fraction of the bracket removed per probe)") print("=" * 78) print(f"{'set':6s}{'method':22s}{'probes':>9s}{'median removed':>17s}" f"{'mean removed':>15s}{'probes skipped':>16s}") print("-" * 86) for setname in ("APS", "CHA"): for nm, m in (("OURS (similarity)", ours), ("brent", brent), ("toms748", toms748), ("bisection", bisection)): shrinks, skipped = [], 0 for s, _pid, f, a, b, root in probs: if s != setname: continue g, xs = traced(f) try: x, _n, conv = m(g, a, b, XTOL, MAXFEV) except Exception: continue if not conv or abs(x - root) > ACC: continue _pos, sh, sk = replay_bracket(f, a, b, xs[2:]) shrinks.extend(sh) skipped += sk if shrinks: srt = sorted(shrinks) med = srt[len(srt) // 2] print(f"{setname:6s}{nm:22s}{len(shrinks):>9d}{med:>17.3f}" f"{sum(shrinks) / len(shrinks):>15.3f}{skipped:>16d}") print("\nBisection removes exactly 0.500 by construction and is the " "reference row.") if __name__ == "__main__": main()