From eeb9102d47a6ad2ab61224eb8ae23542a92a742a Mon Sep 17 00:00:00 2001 From: Valeri Sitnikov Date: Tue, 22 Sep 2026 07:01:13 +0100 Subject: [PATCH] simbracket: A' as a compiled library (Cython, abi3 wheels for Windows and Linux, Python >= 3.10). lib/simbracket/_core.pyx ports exactly the A' path of simroot.pure.pure_root (tent similarity, three-point scale, projection c=2, contra+flatanalogmempt+exp+half+close+sub, pattern) with Python's semantics kept where they differ from C: min/max keep the first argument on ties and NaN, float division by zero raises, math.log raises at 0, and math.hypot and the builtin sum (compensated since 3.12) are called rather than re-coded. find_root(f, a, b, xtol, maxfev) -> (x, nfev, converged); evaluation_bound(a, b, xtol) = 2*ceil(log2(|b-a|/xtol)) + 3. tests/test_simbracket_equivalence.py: bit-identical x and evaluation count against pure_root on all 868 problems of the six corpora (the 91 Roots.jl problems in full), at xtol 1e-12 and 1e-8 and with the budget cut to 12; errors match; the loaded module is the compiled one. 7/7 on Windows and on Linux. check_simbracket.py: removing any one component from the reference makes 36 to 455 problems differ (the test discriminates); the compiled form is 4x faster. The wheel carries only the compiled extension; _core.c is generated (gen_c.py) and ignored. Regression on the server x2: 1232 passed, 0 failed. Co-Authored-By: Claude Opus 5 --- check_simbracket.py | 35 +++ lib/.gitignore | 6 + lib/build_linux.sh | 8 + lib/check_linux.sh | 11 + lib/gen_c.py | 13 + lib/pyproject.toml | 3 + lib/regression_server.sh | 11 + lib/setup.py | 46 +++ lib/simbracket/__init__.py | 13 + lib/simbracket/_core.pyx | 399 +++++++++++++++++++++++++++ tests/test_simbracket_equivalence.py | 89 ++++++ 11 files changed, 634 insertions(+) create mode 100644 check_simbracket.py create mode 100644 lib/.gitignore create mode 100644 lib/build_linux.sh create mode 100644 lib/check_linux.sh create mode 100644 lib/gen_c.py create mode 100644 lib/pyproject.toml create mode 100644 lib/regression_server.sh create mode 100644 lib/setup.py create mode 100644 lib/simbracket/__init__.py create mode 100644 lib/simbracket/_core.pyx create mode 100644 tests/test_simbracket_equivalence.py diff --git a/check_simbracket.py b/check_simbracket.py new file mode 100644 index 0000000..4094044 --- /dev/null +++ b/check_simbracket.py @@ -0,0 +1,35 @@ +"""Two checks of the equivalence test itself. + +1. Mutation: against a reference with ONE component removed the comparison must find + differences - otherwise a passing equivalence test would prove nothing. +2. Speed: wall time of the compiled library against the Python reference over all problems. + + PYTHONPATH=;. python check_simbracket.py +""" +import time + +import simbracket +from simroot.pure import pure_root +from sweep_full import APRIME +from tests.test_simbracket_equivalence import PROBS, reference, run + + +def main(): + for drop in ("close", "sub", "half", "contra", "exp"): + pre = "+".join(p for p in APRIME.split("+") if p != drop) + mut = lambda f, a, b, xt, mf, pre=pre: pure_root( # noqa: E731 + f, a, b, xt, mf, similarity="ratio", scale="derived3", guard="project", + prestep=pre, extrapolate="pattern") + diff = sum(run(mut, f, a, b, 1e-12, 200) != run(simbracket.find_root, f, a, b, 1e-12, 200) + for _p, f, a, b in PROBS) + print(f"mutation -{drop:<7}: {diff} of {len(PROBS)} problems differ (must be > 0)") + for name, m in (("reference (Python)", reference), ("simbracket (compiled)", simbracket.find_root)): + t0 = time.perf_counter() + for _ in range(3): + for _p, f, a, b in PROBS: + run(m, f, a, b, 1e-12, 200) + print(f"{name:<22} {(time.perf_counter() - t0) / 3:.3f} s per pass over {len(PROBS)}") + + +if __name__ == "__main__": + main() diff --git a/lib/.gitignore b/lib/.gitignore new file mode 100644 index 0000000..12291ae --- /dev/null +++ b/lib/.gitignore @@ -0,0 +1,6 @@ +build/ +dist/ +*.egg-info/ +simbracket/_core.c +*.pyd +*.so diff --git a/lib/build_linux.sh b/lib/build_linux.sh new file mode 100644 index 0000000..d1c994b --- /dev/null +++ b/lib/build_linux.sh @@ -0,0 +1,8 @@ +#!/bin/sh +# Build the Linux wheel from the pre-generated C source (no Cython needed on this machine). +set -e +cd "$(dirname "$0")" +rm -rf build dist simbracket.egg-info +~/torchenv/bin/python -m pip wheel . --no-deps --no-build-isolation -w dist +~/torchenv/bin/python -m zipfile -l dist/*.whl +echo BUILD_DONE diff --git a/lib/check_linux.sh b/lib/check_linux.sh new file mode 100644 index 0000000..73badfb --- /dev/null +++ b/lib/check_linux.sh @@ -0,0 +1,11 @@ +#!/bin/sh +# Install the Linux wheel into a private directory and run the bit-identity test against the +# reference implementation in ~/sim-numerics. +set -e +T=~/simbracket_build/site +rm -rf "$T" +~/torchenv/bin/python -m pip install --no-deps --target "$T" ~/simbracket_build/dist/simbracket-1.0.0-cp310-abi3-linux_x86_64.whl +cd ~/sim-numerics +PYTHONPATH="$T:." ~/torchenv/bin/python -m pytest tests/test_simbracket_equivalence.py -q -p no:cacheprovider 2>&1 | tail -5 +PYTHONPATH="$T" ~/torchenv/bin/python -c "import simbracket._core as c; print(c.__file__)" +echo CHECK_DONE diff --git a/lib/gen_c.py b/lib/gen_c.py new file mode 100644 index 0000000..a67b717 --- /dev/null +++ b/lib/gen_c.py @@ -0,0 +1,13 @@ +"""Generate simbracket/_core.c from simbracket/_core.pyx (needs Cython; the build does not). + + python gen_c.py +""" +import os + +from Cython.Build import cythonize + +HERE = os.path.dirname(os.path.abspath(__file__)) +os.chdir(HERE) +cythonize([os.path.join("simbracket", "_core.pyx")], language_level=3, force=True, + compiler_directives={"embedsignature": False, "binding": False}) +print("generated", os.path.join("simbracket", "_core.c")) diff --git a/lib/pyproject.toml b/lib/pyproject.toml new file mode 100644 index 0000000..380744c --- /dev/null +++ b/lib/pyproject.toml @@ -0,0 +1,3 @@ +[build-system] +requires = ["setuptools>=70"] +build-backend = "setuptools.build_meta" diff --git a/lib/regression_server.sh b/lib/regression_server.sh new file mode 100644 index 0000000..67eb28a --- /dev/null +++ b/lib/regression_server.sh @@ -0,0 +1,11 @@ +#!/bin/bash +# Two regression passes on the server with the compiled library on the path, 6 workers +# (40% of the free cores measured before the run), BLAS pinned to one thread per worker. +export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1 +export PYTHONPATH="$HOME/simbracket_build/site:." +cd "$HOME/sim-numerics" || exit 1 +for pass in 1 2; do + echo "=== PASS $pass ===" + ~/torchenv/bin/python -m pytest tests -q -p no:cacheprovider -n 6 2>&1 | tail -4 +done +echo "=== REGRESSION DONE ===" diff --git a/lib/setup.py b/lib/setup.py new file mode 100644 index 0000000..d6dcb6f --- /dev/null +++ b/lib/setup.py @@ -0,0 +1,46 @@ +"""Build the compiled simbracket wheel. + + python gen_c.py (once, where Cython is installed) + python -m pip wheel . --no-deps --no-build-isolation -w dist + +The extension is built from simbracket/_core.c when it exists, so the build machine needs a C +compiler but not Cython. Neither the .pyx nor the .c is put into the wheel. Floating-point +contraction (FMA) is switched off so that the compiled arithmetic is the IEEE arithmetic of +the Python reference, operation for operation. +""" +import os +import sys + +from setuptools import Extension, setup + +HERE = os.path.dirname(os.path.abspath(__file__)) +C_SRC = os.path.join("simbracket", "_core.c") + +if sys.platform == "win32": + cflags = ["/O2", "/fp:precise"] +else: + cflags = ["-O2", "-ffp-contract=off", "-fno-fast-math"] + +ext = Extension( + "simbracket._core", + sources=[C_SRC], + extra_compile_args=cflags, + define_macros=[("Py_LIMITED_API", "0x030A0000")], + py_limited_api=True, +) + +if not os.path.exists(os.path.join(HERE, C_SRC)): + sys.exit("simbracket/_core.c is missing: run gen_c.py where Cython is installed") + +setup( + name="simbracket", + version="1.0.0", + description="Bracketing root finding by similarity and analogy (compiled).", + packages=["simbracket"], + package_data={}, + include_package_data=False, + ext_modules=[ext], + python_requires=">=3.10", + options={"bdist_wheel": {"py_limited_api": "cp310"}}, + zip_safe=False, +) diff --git a/lib/simbracket/__init__.py b/lib/simbracket/__init__.py new file mode 100644 index 0000000..5dc33b6 --- /dev/null +++ b/lib/simbracket/__init__.py @@ -0,0 +1,13 @@ +"""simbracket - bracketing root finding by similarity and analogy (compiled). + + from simbracket import find_root + x, nfev, converged = find_root(f, a, b, xtol=1e-12, maxfev=200) + +f(a) and f(b) must have opposite signs. The run stops when the bracket is no wider than +xtol and never stops on |f|; at most evaluation_bound(a, b, xtol) = +2*ceil(log2(|b - a|/xtol)) + 3 evaluations of f are made. +""" +from ._core import evaluation_bound, find_root + +__version__ = "1.0.0" +__all__ = ["find_root", "evaluation_bound", "__version__"] diff --git a/lib/simbracket/_core.pyx b/lib/simbracket/_core.pyx new file mode 100644 index 0000000..c5fe593 --- /dev/null +++ b/lib/simbracket/_core.pyx @@ -0,0 +1,399 @@ +# cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False +# cython: embedsignature=False, binding=False, always_allow_keywords=True +# +# The compiled A' bracketing root finder. This file is the SOURCE of the distributed library +# and is not shipped; the wheel carries only the compiled extension. +# +# Port of simroot.pure.pure_root restricted to the one configuration the paper proposes: +# similarity="ratio", scale="derived3", guard="project" (budget factor 2), space="x", +# prestep="contra+flatanalogmempt+exp+half+close+sub", points="recent", +# direction="distance", exponent="none", extrapolate="pattern". +# It must reproduce pure_root bit for bit (x and the evaluation count) - pinned by +# tests/test_simbracket_equivalence.py on every problem of the six corpora. Python's own +# semantics are kept wherever they differ from C: min/max keep the FIRST argument on ties and +# NaN, float division by zero raises ZeroDivisionError (cdivision stays off), math.log raises +# ValueError at 0, and the two library calls whose result depends on the interpreter +# (math.hypot and the builtin sum, compensated since Python 3.12) are called, not re-coded. + +from libc.math cimport fabs, log, log2, exp, pow, sqrt, copysign, isfinite, isnan, isinf +from libc.math cimport ceil, INFINITY +from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free + +import builtins as _builtins +from math import hypot as _py_hypot + +_py_sum = _builtins.sum + +__all__ = ["find_root", "evaluation_bound"] + + +class _Budget(Exception): + pass + + +# ---------------------------------------------------------------------------------------- +# Python semantics of the scalar helpers +# ---------------------------------------------------------------------------------------- +cdef inline double _pmax(double a, double b) noexcept: + # builtins.max(a, b): the first argument unless the second is strictly greater + return b if b > a else a + + +cdef inline double _pmin(double a, double b) noexcept: + # builtins.min(a, b): the first argument unless the second is strictly smaller + return b if b < a else a + + +cdef inline double _sim(double z) noexcept: + return _pmax(0.0, 1.0 - fabs(z)) + + +cdef inline bint _opposite(double u, double v) noexcept: + return (u < 0.0 < v) or (v < 0.0 < u) + + +cdef inline bint _same_sign(double u, double v) noexcept: + return (u > 0.0 and v > 0.0) or (u < 0.0 and v < 0.0) + + +cdef double _plog(double v) except? -1.5: + # math.log: nan -> nan, +inf -> +inf, v <= 0 or -inf -> ValueError + if v > 0.0 or isnan(v): + return log(v) + raise ValueError("math domain error") + + +cdef inline double _ratio_sim(double a, double b) except? -1.5: + cdef double hi = _pmax(a, b) + return _pmin(a, b) / hi if hi > 0.0 else 0.0 + + +# ---------------------------------------------------------------------------------------- +# growable double buffer +# ---------------------------------------------------------------------------------------- +cdef struct Buf: + double* p + Py_ssize_t n + Py_ssize_t cap + + +cdef int _push(Buf* b, double v) except -1: + cdef Py_ssize_t cap + cdef double* q + if b.n == b.cap: + cap = 64 if b.cap == 0 else 2 * b.cap + q = PyMem_Realloc(b.p, cap * sizeof(double)) + if q == NULL: + raise MemoryError() + b.p = q + b.cap = cap + b.p[b.n] = v + b.n += 1 + return 0 + + +# ---------------------------------------------------------------------------------------- +# the pieces of the step +# ---------------------------------------------------------------------------------------- +cdef double _slope3(double x2, double f2, double x1, double f1, bint has0, double x0, + double f0) except? -1.5: + cdef double d12, d01, d012, s + if x2 == x1: + return 0.0 + d12 = (f2 - f1) / (x2 - x1) + if not isfinite(d12): + return 0.0 + if not has0 or x0 == x1 or x0 == x2: + return d12 + d01 = (f1 - f0) / (x1 - x0) + d012 = (d12 - d01) / (x2 - x0) + s = d12 + d012 * (x2 - x1) + return s if isfinite(s) else d12 + + +cdef double _order_split(double A, double FA, double B, double FB, double m) except? -1.5: + cdef double lq, t, e + if isinf(m): + return 0.5 * (A + B) + lq = (_plog(fabs(FA)) - _plog(fabs(FB))) / m + if not isfinite(lq): + return 0.5 * (A + B) + if lq > 0.0: + t = 1.0 / (1.0 + exp(-lq)) + else: + e = exp(lq) + t = e / (1.0 + e) + return A + t * (B - A) + + +cdef int _config3(double x, double fx, double xp, double fp, bint has2, double xp2, + double fp2, double* out) except -1: + cdef double xi, phi + if not has2 or xp2 == xp or fp2 == fp: + return 0 + xi = (x - xp) / (xp2 - xp) + phi = (fx - fp) / (fp2 - fp) + if not (isfinite(xi) and isfinite(phi)): + return 0 + out[0] = xi + out[1] = phi + out[2] = fabs(x - xp) + out[3] = fabs(fx - fp) + return 1 + + +cdef int _analogy_order(bint has_cur, double* cur, Buf* past, double* m_out) except -1: + """1 and the order in m_out when the configuration is nearer to its own past than to the + line phi = xi; the translation limit (m = inf) is on.""" + cdef Py_ssize_t i, nb, ib = 0 + cdef double xi, phi, d, dbest = 0.0, lx, lf + cdef double* q + nb = past.n // 4 + if not has_cur or nb == 0: + return 0 + xi = cur[0] + phi = cur[1] + for i in range(nb): + q = past.p + 4 * i + d = _py_hypot(xi - q[0], phi - q[1]) + if i == 0 or d < dbest: + dbest = d + ib = i + q = past.p + 4 * ib + if not ( _py_hypot(xi - q[0], phi - q[1]) < fabs(phi - xi) / sqrt(2.0)): + return 0 + if q[2] <= 0.0 or q[3] <= 0.0 or cur[2] <= 0.0 or cur[3] <= 0.0: + return 0 + lx = cur[2] / q[2] + lf = cur[3] / q[3] + if isfinite(lx) and isfinite(lf) and 0.0 < lf < 1.0 <= lx: + m_out[0] = INFINITY + return 1 + if not (isfinite(lx) and isfinite(lf)) or lx == 1.0 or lf <= 0.0: + return 0 + m_out[0] = _plog(lf) / _plog(lx) + return 1 + + +cdef double _pattern_s(double step, double x, double xp, Buf* qs) except? -1.5: + """The extrapolation share the history confirms (mode "pattern").""" + cdef double d1 = x - xp, s_step, best = 0.0, mean + cdef Py_ssize_t k, p, j + cdef list sims + if step == 0.0 or d1 == 0.0 or (step > 0.0) != (d1 > 0.0): + return 0.0 + s_step = _ratio_sim(fabs(step), fabs(d1)) + k = qs.n + for p in range(2, k // 2 + 1): + sims = [_ratio_sim(qs.p[j], qs.p[j - p]) for j in range(p, k)] + mean = _py_sum(sims) / len(sims) + best = _pmax(best, mean) + if best > 0.5: + # max() can only keep or raise best, so the verdict is already fixed + break + return s_step if best > 0.5 else 0.0 + + +cdef double _call(object f, double x, Py_ssize_t* n, Py_ssize_t limit) except? -1.5: + if n[0] >= limit: + raise _Budget() + n[0] += 1 + return f(x) + + +# ---------------------------------------------------------------------------------------- +# public interface +# ---------------------------------------------------------------------------------------- +def evaluation_bound(double a, double b, double xtol): + """Worst-case number of evaluations of f guaranteed for the bracket [a, b] at tolerance + xtol: twice bisection plus three, 2*ceil(log2(|b - a|/xtol)) + 3.""" + cdef double w = fabs(b - a) + if not (xtol > 0.0): + raise ValueError("xtol must be positive") + if w <= xtol: + return 2 + return 2 * int(ceil(log2(w / xtol))) + 3 + + +def find_root(f, a, b, xtol=1e-12, maxfev=200): + """Root of f in the bracket [a, b], f(a) and f(b) of opposite signs. + + Returns (x, nfev, converged): x the midpoint of the final bracket (or the exact root when + f(x) == 0 was hit), nfev the number of evaluations of f, converged False only when the + evaluation budget maxfev ran out first. Stops when the bracket is no wider than xtol; + never stops on |f|. At most 2*ceil(log2(|b - a|/xtol)) + 3 evaluations.""" + cdef double da = a, db = b, dxtol = xtol + cdef Py_ssize_t limit = maxfev, n = 0 + cdef double x = 0.5 * (da + db), fx, xp, fp, xp2 = 0.0, fp2 = 0.0 + cdef bint has2 = False + cdef double fa, fb, A, FA, B, FB, width, eps, allow, ratio, lh + cdef long long half, n_half + cdef double far, slope, sc, z_now, step = 0.0, mid, loc, brk, f_far, s_eff, longer, reach + cdef double x_new, model_target = 0.0, mid_now, lin_mid, fnew, t_at, lo, hi, left, right + cdef bint has_model, have_step, cur_ok, found, flag, steep, shallow, flag_over, enter + cdef bint last_model = False, in_regime = False, analogy_seen = False + cdef bint shallow_seen = False, has_target = False + cdef double analogy_m = INFINITY, shallow_m = INFINITY, regime_target = 0.0, m_pair = 0.0 + cdef double cur[4] + cdef Buf configs, moves, qs + configs.p = NULL; configs.n = 0; configs.cap = 0 + moves.p = NULL; moves.n = 0; moves.cap = 0 + qs.p = NULL; qs.n = 0; qs.cap = 0 + try: + try: + fa = _call(f, da, &n, limit) + fb = _call(f, db, &n, limit) + if fa == 0.0: + return da, n, True + if fb == 0.0: + return db, n, True + if _same_sign(fa, fb): + raise ValueError("f(a) and f(b) must straddle zero") + + A, FA, B, FB = da, fa, db, fb + if fabs(fa) <= fabs(fb): + x, fx, xp, fp = da, fa, db, fb + else: + x, fx, xp, fp = db, fb, da, fa + + # the projection guard: the bisection budget doubled, fixed by the initial bracket + width = fabs(B - A) + if dxtol == 0.0 or width <= 0.0: + raise ValueError("the projection guard needs xtol and a positive initial width") + eps = 0.5 * dxtol + if width > dxtol: + ratio = width / dxtol + if not (ratio > 0.0) and not isnan(ratio): + raise ValueError("math domain error") + lh = ceil(log2(ratio)) + if isnan(lh): + raise ValueError("cannot convert float NaN to integer") + if isinf(lh): + raise OverflowError("cannot convert float infinity to integer") + half = lh + else: + half = 0 + n_half = ceil(2.0 * half) + if n_half < half: + n_half = half + allow = pow(2.0, n_half) * eps + + while fabs(B - A) > dxtol: + has_model = False + far = B if fabs(B - x) >= fabs(A - x) else A + slope = _slope3(x, fx, xp, fp, has2, xp2, fp2) + sc = fabs(slope) * fabs(far - x) + if sc == 0.0 or not isfinite(sc): + x_new = 0.5 * (A + B) + else: + z_now = fabs(fx) / sc + have_step = False + mid = 0.5 * (A + B) + # the flat regime, recognised by analogy with the run's own past + cur_ok = _config3(x, fx, xp, fp, has2, xp2, fp2, cur) + flag = (not in_regime) and last_model and fp != 0.0 and fx / fp >= 0.25 + found = _analogy_order(cur_ok, cur, &configs, &m_pair) + steep = found and m_pair >= 2.0 + shallow = found and 0.0 < m_pair <= 2.0 / 3.0 + flag_over = ((not in_regime) and last_model and fp != 0.0 + and fx / fp <= -0.25) + if steep: + analogy_seen = True + analogy_m = m_pair + if shallow: + shallow_seen = True + shallow_m = m_pair + enter = flag and analogy_seen + if not enter and flag_over and shallow_seen: + enter = True + analogy_m = shallow_m + if enter: + in_regime = True + if cur_ok: + _push(&configs, cur[0]) + _push(&configs, cur[1]) + _push(&configs, cur[2]) + _push(&configs, cur[3]) + if in_regime: + regime_target = _order_split(A, FA, B, FB, analogy_m) + has_target = True + step = regime_target - x + have_step = True + # the local picture contradicts the bracket + if not have_step and x != xp: + loc = (fx - fp) / (x - xp) + brk = (FB - FA) / (B - A) + if loc != 0.0 and (loc > 0.0) != (brk > 0.0): + step = mid - x + have_step = True + # the half promised by the step against the half the bracket holds + if not have_step: + f_far = FB if far == B else FA + if _sim(z_now) >= 0.5 and fabs(f_far) < fabs(fx): + step = mid - x + have_step = True + if not have_step: + step = (1.0 - _sim(z_now)) * (far - x) + _push(&moves, x - xp) + if moves.n >= 2 and moves.p[moves.n - 2] != 0.0 \ + and moves.p[moves.n - 1] != 0.0: + _push(&qs, fabs(moves.p[moves.n - 1]) / fabs(moves.p[moves.n - 2])) + s_eff = _pattern_s(step, x, xp, &qs) + if 0.0 < s_eff < 1.0: + longer = step / (1.0 - s_eff) + reach = fabs(far - x) + step = longer if fabs(longer) <= reach else copysign(reach, step) + if (fabs(step) < 0.5 * dxtol and last_model and fp != 0.0 + and fabs(fx / fp) < 0.25): + step = copysign(0.5 * dxtol, far - x) + model_target = x + step + has_model = True + x_new = x + step + # projection onto the interval the bisection budget still allows + lo = _pmin(A, B) + hi = _pmax(A, B) + if isfinite(allow): + left = hi - allow + right = lo + allow + if left > right: + x_new = 0.5 * (lo + hi) + else: + x_new = _pmin(_pmax(x_new, left), right) + + if not isfinite(x_new) or not (_pmin(A, B) < x_new < _pmax(A, B)): + x_new = 0.5 * (A + B) + if x_new == x or x_new == xp: + x_new = 0.5 * (A + B) + if x_new == x or x_new == xp: + return x_new, n, fabs(B - A) <= dxtol + + last_model = has_model and x_new == model_target + mid_now = 0.5 * (A + B) + lin_mid = 0.5 * (FA + FB) + fnew = _call(f, x_new, &n, limit) + if in_regime and has_target: + t_at = (x_new - A) / (B - A) + lin_mid = FA * (1.0 - t_at) + FB * t_at + mid_now = x_new + has_target = False + if in_regime and x_new == mid_now and fabs(fnew - lin_mid) < fabs(fnew): + # the point looks like a simple root: leave the regime + in_regime = False + analogy_seen = False + allow *= 0.5 + xp2, fp2, has2 = xp, fp, True + xp, fp = x, fx + x, fx = x_new, fnew + if fnew == 0.0: + return x_new, n, True + if _opposite(FA, fnew): + B, FB = x_new, fnew + else: + A, FA = x_new, fnew + return 0.5 * (A + B), n, True + except _Budget: + return x, n, False + finally: + PyMem_Free(configs.p) + PyMem_Free(moves.p) + PyMem_Free(qs.p) diff --git a/tests/test_simbracket_equivalence.py b/tests/test_simbracket_equivalence.py new file mode 100644 index 0000000..2c567a2 --- /dev/null +++ b/tests/test_simbracket_equivalence.py @@ -0,0 +1,89 @@ +"""The compiled library simbracket must BE the proposed method A': the same iterate and the +same evaluation count as simroot.pure.pure_root, bit for bit, on every problem of the six +corpora (all 91 Roots.jl problems, not only the 78), at the paper's tolerance, at a looser +one, and with the budget cut so that the exhaustion path is exercised too. +""" +from __future__ import annotations + +import math +import random + +import pytest + +simbracket = pytest.importorskip("simbracket") + +from simroot.pure import pure_root # noqa: E402 +from sweep_full import APRIME # noqa: E402 + + +def reference(f, a, b, xtol, maxfev): + return pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived3", + guard="project", prestep=APRIME, extrapolate="pattern") + + +def problems(): + from corpus_applied import build as b_app + from corpus_fourth import build as b_four + from corpus_new import build as b_new + from corpus_roots_jl import build as b_rjl + from sweep_full import load_problems + out = [(f"{s}:{pid}", f, a, b) for s, pid, f, a, b, _r in load_problems()] + out += [(f"NEW:{pid}", f, a, b) for _c, pid, f, a, b, _r in b_new()[0]] + out += [(f"RJL:{pid}", f, a, b) for pid, f, a, b, _ins in b_rjl()[0]] + out += [(f"APP:{pid}", f, a, b) for pid, f, a, b in b_app()[0]] + out += [(f"D:{pid}", f, a, b) for _c, pid, f, a, b, _r in b_four()[0]] + return out + + +def run(method, f, a, b, xtol, maxfev): + import scipy.optimize._tstutils as T + T._f6_cache.clear() + random.seed(20260921) + try: + x, n, ok = method(f, a, b, xtol, maxfev) + return ("ok", float(x).hex() if not math.isnan(x) else "nan", int(n), bool(ok)) + except Exception as e: # noqa: BLE001 - the exception TYPE is part of the behaviour + return ("raise", type(e).__name__) + + +PROBS = problems() + + +def test_the_compiled_module_is_the_one_under_test(): + import simbracket._core as core + assert core.__file__.endswith((".pyd", ".so")), core.__file__ + + +def test_the_corpora_are_all_there(): + # the 855 problems of the six corpora plus the 13 Roots.jl problems that overlap the + # development sets + assert len(PROBS) == 855 + 13 + assert len({p[0] for p in PROBS}) == len(PROBS) + + +@pytest.mark.parametrize("xtol,maxfev", ((1e-12, 200), (1e-8, 200), (1e-12, 12))) +def test_bit_identical_to_the_reference_on_every_problem(xtol, maxfev): + bad = [] + for pid, f, a, b in PROBS: + r = run(reference, f, a, b, xtol, maxfev) + c = run(simbracket.find_root, f, a, b, xtol, maxfev) + if r != c: + bad.append((pid, r, c)) + assert not bad, f"{len(bad)} differ, first: {bad[:5]}" + + +def test_errors_match_the_reference(): + f = lambda t: t - 0.3 # noqa: E731 + for args in ((f, 0.5, 1.0, 1e-12, 200), (f, 0.0, 1.0, 0.0, 200)): + with pytest.raises(ValueError): + reference(*args) + with pytest.raises(ValueError): + simbracket.find_root(*args) + + +def test_the_bound_is_kept_and_reported(): + f = lambda t: t ** 19 # noqa: E731 + x, n, ok = simbracket.find_root(f, -10.0, 100.0, 1e-12, 500) + assert ok and n <= simbracket.evaluation_bound(-10.0, 100.0, 1e-12) + assert simbracket.evaluation_bound(-10.0, 100.0, 1e-12) == \ + 2 * math.ceil(math.log2(110.0 / 1e-12)) + 3