142 lines
6.1 KiB
Python
142 lines
6.1 KiB
Python
"""Stage 2b - is my ITP broken, or is that the method's constants?
|
|
|
|
ITP scored 182 against Brent's 82 on the stage-0 set, which is close to
|
|
bisection. Before ITP is used as the state-of-the-art competitor it has to be
|
|
established which of the two explanations holds, because publishing a
|
|
benchmark against a crippled competitor would be worthless.
|
|
|
|
Trace of the mechanism on sin(x) - 0.5 over [0, 1.5]:
|
|
|
|
x_f = regula falsi point 0.7519
|
|
x_half = 0.75
|
|
delta = kappa1 * (b-a)^kappa2 = (0.2/1.5) * 1.5^2 = 0.300
|
|
|x_half - x_f| = 0.0019 < delta => x_t := x_half
|
|
|
|
so the TRUNCATION step discards the interpolation and bisects. Both delta and
|
|
|x_half - x_f| are O((b-a)^2), so the comparison does not resolve as the
|
|
bracket shrinks: whichever constant is larger wins for the whole run. For a
|
|
nearly symmetric f, kappa1 wins and ITP is bisection.
|
|
|
|
That is the published rule, not a coding error - but it means the authors'
|
|
recommended kappa1 is not neutral across problems. The fair thing to do in
|
|
the benchmark is to give ITP its best kappa1, not its default one. This file
|
|
measures the sensitivity so the benchmark can be honest about it.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
|
|
from simroot.baselines import bisection, brent, itp
|
|
|
|
XTOL = 1e-12
|
|
MAXFEV = 300
|
|
|
|
PROBLEMS = [
|
|
("x^2 - 2", lambda x: x * x - 2.0, 0.0, 2.0, math.sqrt(2.0)),
|
|
("cos(x) - x", lambda x: math.cos(x) - x, 0.0, 1.0, 0.7390851332151607),
|
|
("exp(x) - 2", lambda x: math.exp(x) - 2.0, 0.0, 2.0, math.log(2.0)),
|
|
("x^3 - x - 2", lambda x: x ** 3 - x - 2.0, 1.0, 2.0, 1.5213797068045676),
|
|
("sin(x) - 0.5", lambda x: math.sin(x) - 0.5, 0.0, 1.5, math.asin(0.5)),
|
|
("x^10 - 1", lambda x: x ** 10 - 1.0, 0.0, 1.3, 1.0),
|
|
("exp(-x)(x-1)", lambda x: math.exp(-x) * (x - 1.0), -1.0, 40.0, 1.0),
|
|
("atan(x)", lambda x: math.atan(x), -1.0, 3.0, 0.0),
|
|
]
|
|
|
|
# kappa1 as a multiple of the authors' recommended 0.2/(b-a)
|
|
FACTORS = (1.0, 0.1, 0.01, 1e-3, 1e-4, 1e-6)
|
|
|
|
|
|
def main():
|
|
print("ITP sensitivity to kappa1 (kappa2 = 2, n0 = 1)\n")
|
|
head = (f"{'kappa1':>26s}" + "".join(f"{p[0][:12]:>14s}" for p in PROBLEMS)
|
|
+ f"{'TOTAL':>8s}")
|
|
print(head)
|
|
print("-" * len(head))
|
|
best_total = None
|
|
for fac in FACTORS:
|
|
cells, total = [], 0
|
|
for _, f, a, b, root in PROBLEMS:
|
|
k1 = fac * 0.2 / (b - a)
|
|
x, n, conv = itp(f, a, b, XTOL, MAXFEV, kappa1=k1)
|
|
ok = conv and abs(x - root) <= 1e-8
|
|
cells.append(f"{n:>14d}" if ok else f"{'FAIL(' + str(n) + ')':>14s}")
|
|
total += n
|
|
label = f"{fac:g} x 0.2/(b-a)"
|
|
print(f"{label:>26s}" + "".join(cells) + f"{total:>8d}")
|
|
best_total = total if best_total is None else min(best_total, total)
|
|
|
|
# per-problem best over the sweep: the most generous possible ITP
|
|
per_problem_best, cells = 0, []
|
|
for _, f, a, b, root in PROBLEMS:
|
|
best = None
|
|
for fac in FACTORS:
|
|
x, n, conv = itp(f, a, b, XTOL, MAXFEV, kappa1=fac * 0.2 / (b - a))
|
|
if conv and abs(x - root) <= 1e-8:
|
|
best = n if best is None else min(best, n)
|
|
cells.append(f"{best:>14d}" if best else f"{'FAIL':>14s}")
|
|
per_problem_best += best or MAXFEV
|
|
print("-" * len(head))
|
|
print(f"{'ITP oracle (best k1 per problem)':>26s}"[:26] + "".join(cells) +
|
|
f"{per_problem_best:>8d}")
|
|
|
|
ref = []
|
|
for name, m in (("bisection", bisection), ("brent", brent)):
|
|
cells, total = [], 0
|
|
for _, f, a, b, root in PROBLEMS:
|
|
x, n, conv = m(f, a, b, XTOL, MAXFEV)
|
|
cells.append(f"{n:>14d}")
|
|
total += n
|
|
ref.append((name, cells, total))
|
|
for name, cells, total in ref:
|
|
print(f"{name:>26s}" + "".join(cells) + f"{total:>8d}")
|
|
|
|
# ------------------------------------------------------------------
|
|
# n0 controls the slack in the projection step: r_j = 2^(n_max-j)*eps
|
|
# - (b-a)/2 with n_max = n_half + n0. A larger n0 lets the
|
|
# interpolation roam further before the minmax projection pulls it back
|
|
# to the midpoint. Both problems where the kappa1 sweep flatlined at the
|
|
# bisection count are exactly the ones where interpolation needs room.
|
|
# ------------------------------------------------------------------
|
|
print("\nITP sensitivity to n0 (kappa1 = 0.1 x 0.2/(b-a), kappa2 = 2)\n")
|
|
print(head)
|
|
print("-" * len(head))
|
|
for n0 in (0, 1, 2, 5, 10, 20):
|
|
cells, total = [], 0
|
|
for _, f, a, b, root in PROBLEMS:
|
|
x, n, conv = itp(f, a, b, XTOL, MAXFEV,
|
|
kappa1=0.1 * 0.2 / (b - a), n0=n0)
|
|
ok = conv and abs(x - root) <= 1e-8
|
|
cells.append(f"{n:>14d}" if ok else f"{'FAIL(' + str(n) + ')':>14s}")
|
|
total += n
|
|
print(f"{('n0 = ' + str(n0)):>26s}" + "".join(cells) + f"{total:>8d}")
|
|
|
|
# the most generous ITP possible: best over the whole (kappa1, n0) grid,
|
|
# per problem. If SIM still matches this, the comparison is beyond doubt.
|
|
cells, oracle_total = [], 0
|
|
for _, f, a, b, root in PROBLEMS:
|
|
best = None
|
|
for fac in FACTORS:
|
|
for n0 in (0, 1, 2, 5, 10, 20):
|
|
x, n, conv = itp(f, a, b, XTOL, MAXFEV,
|
|
kappa1=fac * 0.2 / (b - a), n0=n0)
|
|
if conv and abs(x - root) <= 1e-8:
|
|
best = n if best is None else min(best, n)
|
|
cells.append(f"{best:>14d}" if best else f"{'FAIL':>14s}")
|
|
oracle_total += best or MAXFEV
|
|
print("-" * len(head))
|
|
print(f"{'ITP oracle over k1 AND n0':>26s}" + "".join(cells) +
|
|
f"{oracle_total:>8d}")
|
|
print(f"{'brent (no tuning)':>26s}" +
|
|
"".join(f"{n:>14d}" for n in (9, 8, 10, 9, 10, 10, 17, 9)) +
|
|
f"{82:>8d}")
|
|
|
|
print("\nInterpretation: if the oracle row reaches brent, the "
|
|
"implementation is sound\nand ITP's performance is a function of "
|
|
"its three constants; the benchmark must\nthen report ITP at its "
|
|
"best, not at its defaults. If the oracle stays at the\nbisection "
|
|
"count, the implementation is wrong and must be fixed before use.")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|