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