Files
sim-numerics/diag_sign_flips.py
T

105 lines
3.8 KiB
Python

"""Correcting the diagnosis: the step IS a fraction of the bracket, so what stalls?
With S = ratio, 1 - S = min(1, |f|/(|slope| D)), hence step = min(D, |f|/|slope|): the step is
a FRACTION of the distance to the far end and the fraction is the similarity. So "an absolute
length cannot march" was wrong.
The candidate mechanism instead: when the Newton length is much shorter than the bracket the
iterate creeps from ONE SIDE, the sign never changes, and the PROVEN bracket does not shrink
at all — only the guard's bisections shrink it. Chandrupatla, bisecting on 72% of steps and
placing its accepted probes by inverse quadratic interpolation, may instead CROSS the root
and collapse the bracket in one step.
Measured here, per problem: the share of evaluations that flip the sign of f (i.e. actually
shrink the proven bracket), and the geometric bracket reduction per evaluation.
python -u diag_sign_flips.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 sweep_full import ACC, MAXFEV, XTOL
def traced_bracket(f, a, b):
"""Replay a run's probes and count how many flipped the sign, plus the bracket history."""
xs = []
def g(x):
v = f(x)
xs.append((x, v))
return v
return g, xs
def replay(xs, a, b, fa, fb):
"""Walk the probes in order, maintaining the proven bracket; count sign flips."""
A, FA, B, FB = a, fa, b, fb
flips = 0
widths = [abs(B - A)]
for x, v in xs:
if not (min(A, B) < x < max(A, B)):
continue
if FA * v < 0.0:
B, FB = x, v
flips += 1 if (FB * FA < 0.0) else 0
else:
A, FA = x, v
widths.append(abs(B - A))
return flips, widths
def run(method, f, a, b, root):
g, xs = traced_bracket(f, a, b)
try:
x, n, conv = method(g, a, b, XTOL, MAXFEV)
except Exception:
return None
x = float(x)
if not (conv and math.isfinite(x) and (abs(x - root) <= ACC or f(x) == 0.0)):
return None
fa, fb = f(a), f(b)
_flips, widths = replay(xs[2:], a, b, fa, fb)
# how often did a probe actually reduce the proven bracket, and by how much on average
reduced = sum(1 for i in range(1, len(widths)) if widths[i] < widths[i - 1] * 0.999)
halvings = (math.log2(widths[0] / max(widths[-1], 1e-300)) / max(len(xs), 1))
return len(xs), reduced, 100.0 * reduced / max(len(xs), 1), halvings
def main():
print(f"{'problem':<9} {'method':<14} {'evals':>6} {'probes that shrank':>19} "
f"{'share':>7} {'halvings/eval':>14}")
for d in _CHANDRUPATLA_TESTS_DICTS:
pid = str(d["ID"])
if pid not in ("fun1.3", "fun1.4", "fun1.5", "fun3.1", "fun5.1"):
continue
f = d["f"]
a, b = (float(v) for v in d["bracket"])
root = float(d["root"])
ours = run(lambda ff, aa, bb, xt, mf: pure_root(ff, aa, bb, xt, mf, similarity="ratio",
scale="derived3", guard="safe",
prestep="slope",
direction="dissimilar"),
f, a, b, root)
best = None
for impl in (chandrupatla, chandrupatla_scipy):
r = run(impl, f, a, b, root)
if r and (best is None or r[0] < best[0]):
best = r
for nm, r in (("ours", ours), ("chandrupatla", best)):
if r is None:
print(f"{pid:<9} {nm:<14} (not scorable)")
continue
n, red, share, halv = r
print(f"{pid:<9} {nm:<14} {n:6d} {red:19d} {share:6.1f}% {halv:14.3f}")
if __name__ == "__main__":
main()