Files
sim-numerics/diag_fun9_well.py

68 lines
2.6 KiB
Python

"""Why the fun9 well is closed only partly (Valeri, 2026-09-21).
fun9 = exp(x) - 2 - 0.01/x^2 + 2e-6/x^3, brackets [2e-4, 1 / 3 / 9 / 27 / 81], root 0.7032.
Every evaluation is recorded for the library Chandrupatla and for our variants. A run is
split at its first evaluation with f > 0 inside the bracket (the left end is negative, the
right end positive): before it the method is still looking for the root's side of the well,
after it it converges. Our stats give how many steps the decision/contra/regime refused
(pre_refused) and how many the safeguard forced. Nothing in the method changes.
python -u diag_fun9_well.py counts and phases
python -u diag_fun9_well.py fun9.4 plus the point sequences for one problem
"""
from __future__ import annotations
import sys
from simroot.baselines import chandrupatla_lib
from simroot.pure import pure_root
from sweep_full import MAXFEV, XTOL, load_problems
VARIANTS = ("none", "contra", "defect", "contra+flatanalogmempt", "defect+flatanalogmempt")
def recorded(f):
pts = []
def g(x):
v = f(x)
pts.append((float(x), float(v)))
return v
return g, pts
def phases(pts, a, b):
inner = [i for i, (x, v) in enumerate(pts) if x not in (a, b)]
first_pos = next((i for i in inner if pts[i][1] > 0.0), None)
n = len(pts)
return (first_pos + 1 if first_pos is not None else n), (n - first_pos - 1 if first_pos is not None else 0)
def main(show=None):
print(f"{'problem':<8} {'method':<26} {'n':>4} {'to first f>0':>13} {'after':>6} "
f"{'refused':>8} {'forced':>7}")
for s, pid, f, a, b, root in load_problems():
if not pid.startswith("fun9"):
continue
g, pts = recorded(f)
chandrupatla_lib(g, a, b, XTOL, MAXFEV)
p1, p2 = phases(pts, a, b)
print(f"{pid:<8} {'Chandrupatla':<26} {len(pts):4d} {p1:13d} {p2:6d}")
if pid == show:
print(" Chandrupatla: " + " ".join(f"{x:.4g}:{v:+.2g}" for x, v in pts))
for pre in VARIANTS:
g, pts = recorded(f)
st = {}
pure_root(g, a, b, XTOL, MAXFEV, similarity="ratio", scale="derived3",
guard="project", prestep=pre, stats=st)
p1, p2 = phases(pts, a, b)
print(f"{'':<8} {pre:<26} {len(pts):4d} {p1:13d} {p2:6d} "
f"{st.get('pre_refused', 0):8d} {st.get('forced', 0):7d}")
if pid == show:
print(f" {pre}: " + " ".join(f"{x:.4g}:{v:+.2g}" for x, v in pts))
print("diag done")
if __name__ == "__main__":
main(*sys.argv[1:2])