diff --git a/article/sim_numerics.pdf b/article/sim_numerics.pdf index 41423f0..8fe3ae2 100644 Binary files a/article/sim_numerics.pdf and b/article/sim_numerics.pdf differ diff --git a/article/sim_numerics.tex b/article/sim_numerics.tex index 289e334..c2d5ff0 100644 --- a/article/sim_numerics.tex +++ b/article/sim_numerics.tex @@ -609,12 +609,41 @@ of the remaining gap on the wide-bracket set: $253$ evaluations, or $17\%$, not $11\%$.} \end{table} -Two things follow. The wide-bracket gap to Chandrupatla is real and it is -\emph{diffuse}: measured per family it is $+57$, $+54$, $+46$, $+43$, $+35$ and -smaller, with no single problem costing more than $35$, so it is not the kind of -deficit a single mechanism removes. And on that set our shape reading ($2159$ -re-scored) is \emph{worse} than bisection ($2136$); only the dissimilarity -reading ($1740$) is clearly ahead. Every statement in this subsection is +\subsection{Where the wide-bracket limit actually is} +\label{sec:widelimit} + +The remaining gap looks diffuse --- per family $+57$, $+54$, $+46$, $+43$, $+35$ +and smaller --- but it has one address. The \texttt{fun1} family is the same +function with brackets of growing width, and our cost against Chandrupatla's +across them is $7/8$, $12/12$, $19/15$, $59/24$, $63/44$: we tie or win while the +bracket is narrow and lose progressively as it widens. The endgame is not the +cause --- our observed order on those very problems is the higher one ($3.63$ +against $2.02$ on the $2\times10^4$-wide bracket). The cause is the rate of +progress far from the root: $59$ evaluations on a $2\times10^4$-wide bracket is +what bisection itself needs ($\lceil\log_2(2\times10^{16})\rceil\approx54$), i.e. +there we advance by one halving per evaluation while Chandrupatla advances by +about $2.3$. + +It is not the guarantee that caps us. With the safeguard removed entirely the +same problems cost $46$, $57$, $61$ and $68$ against bisection's $57$, $77$, $45$ +and $45$ --- better on the wide brackets, \emph{worse than bisection} on +\texttt{fun3.1} and \texttt{fun5.1} --- and still far from Chandrupatla's $24$ and +$44$. The cap is the step rule itself: $|f|/|f'|$ is an absolute length, and an +absolute length cannot march across decades. That is the same wall the fixed +fraction (\texttt{cbrt}) and the logarithmic change of variable each hit from a +different side, both buying the wide brackets at the price of the canonical set. + +So the limit is stated exactly: the placement of the similarity reading improves +the method up to the point where the \emph{form} of the derived step becomes the +binding constraint, and on brackets spanning many decades it does. Removing that +constraint means replacing the step the contraction law derives, which is a +different method rather than a better reading of this one. + +On that set our shape reading ($2159$ re-scored) is \emph{worse} than bisection +($2136$); only the dissimilarity reading ($1740$) is clearly ahead, which is the +same statement seen from the other end: where the step cannot march, the reading +that helps least is the one that only improves the local model. Every statement +in these two subsections is reproduced by \texttt{diag\_rescore\_cha.py}, \texttt{diag\_rescore\_aps.py}, \texttt{diag\_fun7\_root.py} and \texttt{diag\_cha\_gap.py}, whose outputs are kept beside them. diff --git a/diag_cha_phases.py b/diag_cha_phases.py new file mode 100644 index 0000000..5d32ed2 --- /dev/null +++ b/diag_cha_phases.py @@ -0,0 +1,105 @@ +"""WHERE in a run the wide-bracket gap to Chandrupatla is spent: approach or endgame? + +The gap is diffuse across families (+57, +54, +46, +43, +35, ...), so before proposing any +mechanism this splits every run into the two phases with the paper's OWN definition +(stage4_why.phase_boundary: the geometric mean of the initial width and the ATTAINABLE +tolerance) and reports where each method's evaluations go. + +Problems whose declared root cannot be matched are skipped, and the six scoring-artefact +problems (aps.13.00, fun7.*) are reported separately rather than silently dropped. + + python -u diag_cha_phases.py +""" + +from __future__ import annotations + +import math +from collections import defaultdict + +from simroot.baselines import chandrupatla, chandrupatla_scipy +from simroot.pure import pure_root +from stage4_why import phase_boundary, traced +from sweep_full import ACC, MAXFEV, XTOL, load_problems + +# NOT best_of here: it runs BOTH implementations and returns the cheaper result, so a tracer +# wrapped around it counts the evaluations of both. Doing that reported Chandrupatla's CHA +# total as 2993 instead of 1487 and manufactured a 1299-evaluation "endgame win" for us +# (caught 2026-09-21). Each implementation is traced separately and the cheaper split kept, +# which is what best_of means for cost. +CHA_IMPLS = (chandrupatla, chandrupatla_scipy) + + +def ours(f, a, b, xtol, maxfev): + return pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived3", + guard="safe", prestep="slope", direction="dissimilar") + + +def phases(method, f, a, b, root): + """(approach, endgame, total) evaluations, or None when the run cannot be scored.""" + g, xs = traced(f) + try: + x, _n, conv = method(g, a, b, XTOL, MAXFEV) + except Exception: + return None + x = float(x) + ok = conv and math.isfinite(x) and (abs(x - root) <= ACC or f(x) == 0.0) + if not ok: + 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 cha_phases(f, a, b, root): + """The cheaper of the two Chandrupatla implementations, each traced on its own.""" + best = None + for impl in CHA_IMPLS: + p = phases(impl, f, a, b, root) + if p is not None and (best is None or p[2] < best[2]): + best = p + return best + + +def main(): + agg = defaultdict(lambda: [0, 0, 0]) # method -> [approach, endgame, total] + per_fam = defaultdict(lambda: [0, 0, 0, 0]) # family -> [our_app, our_end, ch_app, ch_end] + skipped = [] + print(f"{'problem':<10} {'ours app/end':>14} {'CH app/end':>12} {'gap app':>8} {'gap end':>8}") + for s, pid, f, a, b, root in load_problems(): + if s != "CHA": + continue + po = phases(ours, f, a, b, root) + pc = cha_phases(f, a, b, root) + if po is None or pc is None: + skipped.append(pid) + continue + agg["ours"][0] += po[0]; agg["ours"][1] += po[1]; agg["ours"][2] += po[2] + agg["chandrupatla"][0] += pc[0]; agg["chandrupatla"][1] += pc[1] + agg["chandrupatla"][2] += pc[2] + fam = pid.split(".")[0] + per_fam[fam][0] += po[0]; per_fam[fam][1] += po[1] + per_fam[fam][2] += pc[0]; per_fam[fam][3] += pc[1] + print(f"{pid:<10} {po[0]:6d}/{po[1]:<7d} {pc[0]:5d}/{pc[1]:<6d} " + f"{po[0] - pc[0]:+8d} {po[1] - pc[1]:+8d}") + + print(f"\nskipped (not scorable under either criterion): {skipped}") + for nm, (app, end, tot) in agg.items(): + print(f"{nm:<14} approach {app:5d} endgame {end:5d} total {tot:5d} " + f"({100.0 * app / max(tot, 1):.0f}% of the cost is approach)") + d_app = agg["ours"][0] - agg["chandrupatla"][0] + d_end = agg["ours"][1] - agg["chandrupatla"][1] + print(f"\nGAP: approach {d_app:+d}, endgame {d_end:+d} " + f"(total {d_app + d_end:+d})") + print("\nper family (approach gap / endgame gap):") + for fam in sorted(per_fam, key=lambda k: -((per_fam[k][0] - per_fam[k][2]) + + (per_fam[k][1] - per_fam[k][3]))): + oa, oe, ca, ce = per_fam[fam] + print(f" {fam:<8} approach {oa - ca:+5d} endgame {oe - ce:+5d}") + + +if __name__ == "__main__": + main() diff --git a/diag_fun1_endgame.py b/diag_fun1_endgame.py new file mode 100644 index 0000000..0a34618 --- /dev/null +++ b/diag_fun1_endgame.py @@ -0,0 +1,75 @@ +"""The largest single item of the wide-bracket gap: fun1's endgame (+61 evaluations). + +If it has a structural cause, a derived fix may exist; if it is simply the order of +convergence, the answer is "no" and the reason is worth stating. Reports, per fun1 problem: +the bracket, the multiplicity-like behaviour of f near the root, the observed order of both +methods, and the endgame split. + + python -u diag_fun1_endgame.py +""" + +from __future__ import annotations + +import math + +from scipy.optimize._tstutils import _CHANDRUPATLA_TESTS_DICTS + +from simroot.baselines import chandrupatla, chandrupatla_scipy +from simroot.pure import pure_root +from stage4_why import observed_order +from sweep_full import ACC, MAXFEV, XTOL + + +def ours(f, a, b, xtol, maxfev): + return pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived3", + guard="safe", prestep="slope", direction="dissimilar") + + +def multiplicity_probe(f, root): + """How does |f| grow away from the root? |f(root+t)| ~ C t^m gives m from two points; + m = 1 is a simple root, m >= 2 flattens the endgame for every interpolating method.""" + t1, t2 = 1e-5, 1e-4 + f1, f2 = abs(f(root + t1)), abs(f(root + t2)) + if f1 == 0.0 or f2 == 0.0: + return float("nan") + return math.log(f2 / f1) / math.log(t2 / t1) + + +def main(): + print(f"{'problem':<9} {'bracket':>26} {'root':>10} {'m':>6} " + f"{'ours n':>7} {'ord':>6} {'CH n':>6} {'ord':>6}") + for d in _CHANDRUPATLA_TESTS_DICTS: + pid = str(d["ID"]) + if not pid.startswith("fun1."): + continue + f = d["f"] + a, b = (float(v) for v in d["bracket"]) + root = float(d["root"]) + m = multiplicity_probe(f, root) + x_o, n_o, ok_o = ours(f, a, b, XTOL, MAXFEV) + best = None + for impl in (chandrupatla, chandrupatla_scipy): + try: + x_c, n_c, ok_c = impl(f, a, b, XTOL, MAXFEV) + except Exception: + continue + if best is None or n_c < best[1]: + best = (float(x_c), int(n_c), bool(ok_c)) + try: + ord_o = observed_order(ours, f, a, b, root) + except Exception: + ord_o = float("nan") + try: + ord_c = observed_order(chandrupatla, f, a, b, root) + except Exception: + ord_c = float("nan") + print(f"{pid:<9} [{a:10.4g}, {b:10.4g}] {root:10.5g} {m:6.2f} " + f"{n_o:7d} {ord_o if isinstance(ord_o, float) else float('nan'):6.2f} " + f"{best[1]:6d} {ord_c if isinstance(ord_c, float) else float('nan'):6.2f}") + print(f" |f| at root+1e-5 = {abs(f(root + 1e-5)):.3e}, " + f"root+1e-2 = {abs(f(root + 1e-2)):.3e}, " + f"ours x={x_o:+.10g} ok={ok_o}") + + +if __name__ == "__main__": + main() diff --git a/diag_who_caps_us.py b/diag_who_caps_us.py new file mode 100644 index 0000000..2eac75c --- /dev/null +++ b/diag_who_caps_us.py @@ -0,0 +1,70 @@ +"""On wide brackets we spend what bisection spends. WHO caps us — the step rule or the guard? + +Measured (diag_fun1_endgame.py): on the same function our cost grows with the bracket width +(7, 12, 19, 59, 63 against Chandrupatla's 8, 12, 15, 24, 44) while our observed ORDER is +higher than theirs, so the endgame is not the problem. 59 evaluations on a 2e4-wide bracket +is what bisection needs (log2(2e16) ~ 54), i.e. we degrade to one halving per evaluation. + +Two candidate causes, and they are distinguishable by a run: + * the STEP RULE cannot march (the similarity step is an absolute length), or + * the GUARANTEE's projection clamps the accepted step to a halving-sized move. + +So the same runs are repeated with each guard, including guard="none" (no guarantee at all). +If "none" marches, the guard is the cap; if it also crawls, the step rule is. + + python -u diag_who_caps_us.py +""" + +from __future__ import annotations + +import math + +from scipy.optimize._tstutils import _CHANDRUPATLA_TESTS_DICTS + +from simroot.baselines import bisection, chandrupatla, chandrupatla_scipy +from simroot.pure import pure_root +from sweep_full import ACC, MAXFEV, XTOL + +GUARDS = ("none", "step", "project", "safe", "credit") + + +def main(): + print(f"{'problem':<9} {'width':>10} {'bisect':>7} {'CH':>5} " + + " ".join(f"{g:>8}" for g in GUARDS)) + for d in _CHANDRUPATLA_TESTS_DICTS: + pid = str(d["ID"]) + if not (pid.startswith("fun1.") or pid in ("fun3.1", "fun5.1", "fun4.3")): + continue + f = d["f"] + a, b = (float(v) for v in d["bracket"]) + root = float(d["root"]) + ok = lambda x, c: bool(c and math.isfinite(float(x)) # noqa: E731 + and (abs(float(x) - root) <= ACC or f(float(x)) == 0.0)) + _xb, nb, cb = bisection(f, a, b, XTOL, MAXFEV) + best = None + for impl in (chandrupatla, chandrupatla_scipy): + try: + xc, nc, cc = impl(f, a, b, XTOL, MAXFEV) + except Exception: + continue + if ok(xc, cc) and (best is None or nc < best): + best = int(nc) + cells = [] + for g in GUARDS: + try: + x, n, c = pure_root(f, a, b, XTOL, MAXFEV, similarity="ratio", + scale="derived3", guard=g, prestep="slope", + direction="dissimilar") + cells.append(f"{n:8d}" if ok(x, c) else f"{n:7d}✗") + except Exception: + cells.append(f"{'err':>8}") + print(f"{pid:<9} {abs(b - a):10.3g} {nb:7d} {best if best else 0:5d} " + + " ".join(cells)) + + print("\nlegend: ✗ = the run did not reach the root (an unguarded method may not)") + print("if guard='none' marches where 'safe'/'project' crawl, the GUARANTEE is the cap;") + print("if 'none' crawls too, the STEP RULE is the cap.") + + +if __name__ == "__main__": + main() diff --git a/results_cha_phases.txt b/results_cha_phases.txt new file mode 100644 index 0000000..3fd15ff --- /dev/null +++ b/results_cha_phases.txt @@ -0,0 +1,63 @@ +problem ours app/end CH app/end gap app gap end +fun1.1 5/2 6/2 -1 +0 +fun1.2 10/2 10/2 +0 +0 +fun1.3 11/8 13/2 -2 +6 +fun1.4 20/39 21/3 -1 +36 +fun1.5 3/60 3/41 +0 +19 +fun2.1 7/2 7/2 +0 +0 +fun2.2 23/3 20/3 +3 +0 +fun2.3 27/3 25/4 +2 -1 +fun2.4 1/42 2/40 -1 +2 +fun2.5 1/48 2/47 -1 +1 +fun3.1 32/32 21/23 +11 +9 +fun3.2 28/28 23/23 +5 +5 +fun3.3 31/35 28/32 +3 +3 +fun3.4 31/42 28/41 +3 +1 +fun3.5 3/82 3/79 +0 +3 +fun4.1 30/24 21/23 +9 +1 +fun4.2 30/26 23/23 +7 +3 +fun4.3 32/38 29/27 +3 +11 +fun4.4 28/45 30/37 -2 +8 +fun4.5 3/83 3/77 +0 +6 +fun5.1 30/27 22/22 +8 +5 +fun5.2 28/26 22/22 +6 +4 +fun5.3 27/28 22/23 +5 +5 +fun5.4 31/26 22/25 +9 +1 +fun5.5 31/28 22/26 +9 +2 +fun6.1 29/23 22/22 +7 +1 +fun6.2 29/25 22/22 +7 +3 +fun6.3 25/25 22/23 +3 +2 +fun6.4 25/27 22/25 +3 +2 +fun6.5 25/30 22/26 +3 +4 +fun7.1 11/0 8/0 +3 +0 +fun7.2 11/0 8/0 +3 +0 +fun7.3 11/0 11/0 +0 +0 +fun7.4 11/0 18/0 -7 +0 +fun7.5 13/0 19/0 -6 +0 +fun8.1 11/2 8/2 +3 +0 +fun8.2 9/2 9/2 +0 +0 +fun8.3 10/2 10/3 +0 -1 +fun8.4 13/3 12/2 +1 +1 +fun8.5 16/2 13/2 +3 +0 +fun9.1 6/2 6/2 +0 +0 +fun9.2 10/2 7/2 +3 +0 +fun9.3 11/3 9/2 +2 +1 +fun9.4 13/2 10/3 +3 -1 +fun9.5 16/3 12/2 +4 +1 + +skipped (not scorable under either criterion): [] +ours approach 808 endgame 932 total 1740 (46% of the cost is approach) +chandrupatla approach 698 endgame 789 total 1487 (47% of the cost is approach) + +GAP: approach +110, endgame +143 (total +253) + +per family (approach gap / endgame gap): + fun1 approach -4 endgame +61 + fun5 approach +37 endgame +17 + fun4 approach +17 endgame +29 + fun3 approach +22 endgame +21 + fun6 approach +23 endgame +12 + fun9 approach +12 endgame +1 + fun8 approach +7 endgame +0 + fun2 approach +3 endgame +2 + fun7 approach -7 endgame +0 diff --git a/results_fun1_endgame.txt b/results_fun1_endgame.txt new file mode 100644 index 0000000..c26985c --- /dev/null +++ b/results_fun1_endgame.txt @@ -0,0 +1,11 @@ +problem bracket root m ours n ord CH n ord +fun1.1 [ 2, 3] 2.0946 1.00 7 4.28 8 1.80 + |f| at root+1e-5 = 1.116e-04, root+1e-2 = 1.122e-01, ours x=+2.094551482 ok=True +fun1.2 [ 1, 10] 2.0946 1.00 12 1.87 12 1.81 + |f| at root+1e-5 = 1.116e-04, root+1e-2 = 1.122e-01, ours x=+2.094551482 ok=True +fun1.3 [ 1, 100] 2.0946 1.00 19 2.46 15 1.67 + |f| at root+1e-5 = 1.116e-04, root+1e-2 = 1.122e-01, ours x=+2.094551482 ok=True +fun1.4 [ -1e+04, 1e+04] 2.0946 1.00 59 3.63 24 2.02 + |f| at root+1e-5 = 1.116e-04, root+1e-2 = 1.122e-01, ours x=+2.094551482 ok=True +fun1.5 [ -1e+10, 1e+10] 2.0946 1.00 63 1.85 44 2.20 + |f| at root+1e-5 = 1.116e-04, root+1e-2 = 1.122e-01, ours x=+2.094551482 ok=True diff --git a/results_who_caps_us.txt b/results_who_caps_us.txt new file mode 100644 index 0000000..a633a94 --- /dev/null +++ b/results_who_caps_us.txt @@ -0,0 +1,13 @@ +problem width bisect CH none step project safe credit +fun1.1 1 42 8 7 7 7 7 19 +fun1.2 9 46 12 12 12 12 12 45 +fun1.3 99 49 15 30 19 30 19 50 +fun1.4 2e+04 57 24 46 59 46 59 58 +fun1.5 2e+10 77 44 57 63 57 63 76 +fun3.1 5 45 44 61 64 61 64 46 +fun4.3 2e+04 57 56 73 70 73 70 58 +fun5.1 5 45 44 68 57 68 57 46 + +legend: ✗ = the run did not reach the root (an unguarded method may not) +if guard='none' marches where 'safe'/'project' crawl, the GUARANTEE is the cap; +if 'none' crawls too, the STEP RULE is the cap.