diff --git a/.github/workflows/build-wheels.yml b/.github/workflows/build-wheels.yml index 517cd55..125c28c 100644 --- a/.github/workflows/build-wheels.yml +++ b/.github/workflows/build-wheels.yml @@ -52,7 +52,7 @@ jobs: uses: pypa/cibuildwheel@v4.3.1 env: # All cibuildwheel options are here (nothing "hidden" in pyproject.toml) - CIBW_SKIP: "*-win32 *i686 cp3??t-* cp310-win_arm64" + CIBW_SKIP: "*-win32 *i686 cp313t-* cp310-win_arm64" CIBW_BUILD_FRONTEND: build CIBW_BUILD_VERBOSITY: 3 CIBW_TEST_REQUIRES: pytest scipy diff --git a/.github/workflows/test-linux.yml b/.github/workflows/test-linux.yml index 0cb36c9..b7de300 100644 --- a/.github/workflows/test-linux.yml +++ b/.github/workflows/test-linux.yml @@ -13,7 +13,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.10", "3.11", "3.12", "3.13", "3.14", "3.15"] + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14", "3.14t", "3.15"] resolution: ["lowest-direct", "highest"] exclude: - python-version: "3.12" @@ -22,6 +22,8 @@ jobs: resolution: "lowest-direct" - python-version: "3.14" resolution: "lowest-direct" + - python-version: "3.14t" + resolution: "lowest-direct" - python-version: "3.15" resolution: "lowest-direct" steps: @@ -43,8 +45,6 @@ jobs: - name: Install Ubuntu-build dependencies run: | - # sudo apt-get update - # sudo apt-get install -y libopenblas-dev gfortran --version gcc --version diff --git a/README.md b/README.md index 30c160c..f3369f2 100644 --- a/README.md +++ b/README.md @@ -14,9 +14,7 @@ of the original Fortran implementation, while working within the modern Python e ODR, also known as [errors-in-variables regression], is designed primarily for instances when both the explanatory and response variables have significant errors. -

- Deming regression; special case of ODR. -

+![Deming regression; special case of ODR.](https://github.com/HugoMVale/odrpack-python/blob/main/docs/odr.png) [errors-in-variables regression]: https://en.wikipedia.org/wiki/Errors-in-variables_models [odrpack95]: https://github.com/HugoMVale/odrpack95 diff --git a/meson.build b/meson.build index 2010ede..8e66718 100644 --- a/meson.build +++ b/meson.build @@ -1,7 +1,7 @@ project( 'odrpack', ['cpp', 'fortran'], - version : '0.7.0', + version : '0.7.1', meson_version: '>=1.5.0', default_options: [ 'buildtype=release', @@ -15,6 +15,15 @@ project( python = import('python').find_installation(pure: false) python_dep = python.dependency() +# Free-threaded CPython (3.13t, 3.14t, ...): nanobind -- both the static library +# and the extension -- must be compiled with NB_FREE_THREADED. Otherwise the +# module does not declare Py_MOD_GIL_NOT_USED (CPython silently re-enables the +# GIL on import) and nanobind's internals are not thread-safe. Global arguments +# are used on purpose, because they also reach the nanobind subproject. +if python.get_variable('Py_GIL_DISABLED', 0) == 1 + add_global_arguments('-DNB_FREE_THREADED', language: 'cpp') +endif + # nanobind (from subprojects) nanobind_dep = dependency('nanobind', static: true) diff --git a/src/odrpack/__odrpack.py.cpp b/src/odrpack/__odrpack.py.cpp index 3e143ba..6bde26d 100644 --- a/src/odrpack/__odrpack.py.cpp +++ b/src/odrpack/__odrpack.py.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -24,13 +25,23 @@ struct Context { nb::callable fcn_f; nb::callable fcn_fjacb; nb::callable fcn_fjacd; + nb::object odr_stop; + std::exception_ptr error; }; /* Callback function invoked from `odrpack` to evaluate the model function and, optionally, its -Jacobians. The actual python functions are stored in the `Context` struct, which is passed through -the void pointer argument `data`. This is used to call the appropriate functions without -creating closures. +Jacobians. The actual python functions are stored in the `Context` struct, which is passed +through the void pointer argument `data`. This is used to call the appropriate functions +without creating closures. + +Exceptions must never propagate out of this function: it is called from Fortran, and unwinding +a C++ exception through Fortran frames is not guaranteed to work. Therefore: +- `OdrStop` (or a subclass): the regression is stopped gracefully (`istop = 1`) and a solution + is returned. +- Any other exception (Python or C++): it is stored in the context, the regression is aborted + (`istop = -1`) and the exception is rethrown by the caller after the Fortran routine has + returned normally. */ void fcn(const int *n_ptr, const int *m_ptr, const int *q_ptr, const int *npar_ptr, const int *ldifx_ptr, const double beta[], const double xplusd[], const int ifixb[], const int ifixx[], @@ -50,7 +61,13 @@ void fcn(const int *n_ptr, const int *m_ptr, const int *q_ptr, const int *npar_p (m == 1) ? std::initializer_list{n} : std::initializer_list{m, n}); // Retrieve the model context - const auto *context = static_cast(data); + auto *context = static_cast(data); + + // A previous callback already failed: do not call back into Python, just keep aborting. + if (context->error) { + *istop = -1; + return; + } *istop = 0; try { @@ -76,14 +93,18 @@ void fcn(const int *n_ptr, const int *m_ptr, const int *q_ptr, const int *npar_p } } catch (const nb::python_error &e) { - // temporary solution: need to figure out how to do this the right way - std::string ewhat = e.what(); - if (ewhat.find("OdrStop") != std::string::npos) { - std::cerr << ewhat << std::endl; + if (e.matches(context->odr_stop)) { + // `OdrStop` (or a subclass) raised by the user: stop gracefully + std::cerr << e.what() << std::endl; *istop = 1; } else { - throw; + context->error = std::current_exception(); + *istop = -1; } + } catch (...) { + // Any other exception + context->error = std::current_exception(); + *istop = -1; } }; @@ -175,8 +196,9 @@ int odr_wrapper( if (rwork) lrwork = rwork.value().size(); if (iwork) liwork = iwork.value().size(); - // Define the context for the user-supplied model function and its Jacobians. - Context context = {fcn_f, fcn_fjacb, fcn_fjacd}; + // Define the context for the user-supplied model function and its Jacobians + nb::object odr_stop = nb::module_::import_("odrpack.exceptions").attr("OdrStop"); + Context context = {fcn_f, fcn_fjacb, fcn_fjacd, odr_stop}; // Open files int lunrpt = 6; @@ -221,6 +243,11 @@ int odr_wrapper( if (ierr != 0) throw std::runtime_error("Error closing error file."); } + // Rethrow the exception (if any) that aborted the regression in a callback. + // This is done only now, after the Fortran routine has returned normally and the files + // have been closed. + if (context.error) std::rethrow_exception(context.error); + return info; } diff --git a/subprojects/nanobind.wrap b/subprojects/nanobind.wrap index 2c8768f..75b2ffd 100644 --- a/subprojects/nanobind.wrap +++ b/subprojects/nanobind.wrap @@ -1,14 +1,14 @@ [wrap-file] -directory = nanobind-2.12.0 -source_url = https://github.com/wjakob/nanobind/archive/refs/tags/v2.12.0.tar.gz -source_filename = nanobind-2.12.0.tar.gz -source_hash = 01f1f0cd0398743c18f33d07ae36ad410bd7f4a1e90683b508504de897d6e629 -source_fallback_url = https://wrapdb.mesonbuild.com/v2/nanobind_2.12.0-1/get_source/nanobind-2.12.0.tar.gz -patch_filename = nanobind_2.12.0-1_patch.zip -patch_url = https://wrapdb.mesonbuild.com/v2/nanobind_2.12.0-1/get_patch -patch_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/nanobind_2.12.0-1/nanobind_2.12.0-1_patch.zip -patch_hash = 4b8158bdb359218bcfb7a5b4459fbfa5919171d9ea2400de3a84cf9ab7c7308b -wrapdb_version = 2.12.0-1 +directory = nanobind-2.15.0 +source_url = https://github.com/wjakob/nanobind/archive/refs/tags/v2.15.0.tar.gz +source_filename = nanobind-2.15.0.tar.gz +source_hash = 36c8760b3acb25643cd89d549782d8c67bd82ce54c6238608787c22a34fd490f +source_fallback_url = https://wrapdb.mesonbuild.com/v2/nanobind_2.15.0-1/get_source/nanobind-2.15.0.tar.gz +patch_filename = nanobind_2.15.0-1_patch.zip +patch_url = https://wrapdb.mesonbuild.com/v2/nanobind_2.15.0-1/get_patch +patch_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/nanobind_2.15.0-1/nanobind_2.15.0-1_patch.zip +patch_hash = 4a3956c67792830e0c85dff3cf88b30a5547bd087c33842fa5f1596ca86ccabf +wrapdb_version = 2.15.0-1 [provide] nanobind = nanobind_dep diff --git a/subprojects/robin-map.wrap b/subprojects/robin-map.wrap index 0bea538..38e327e 100644 --- a/subprojects/robin-map.wrap +++ b/subprojects/robin-map.wrap @@ -1,13 +1,14 @@ [wrap-file] -directory = robin-map-1.4.0 -source_url = https://github.com/Tessil/robin-map/archive/refs/tags/v1.4.0.tar.gz -source_filename = robin-map-1.4.0.tar.gz -source_hash = 7930dbf9634acfc02686d87f615c0f4f33135948130b8922331c16d90a03250c -patch_filename = robin-map_1.4.0-1_patch.zip -patch_url = https://wrapdb.mesonbuild.com/v2/robin-map_1.4.0-1/get_patch -patch_hash = feb14b6752b7d439fb2f3ee968e595a9a3de00ef8cb029488af2ceb4f504b95d -source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/robin-map_1.4.0-1/robin-map-1.4.0.tar.gz -wrapdb_version = 1.4.0-1 +directory = robin-map-1.4.1 +source_url = https://github.com/Tessil/robin-map/archive/refs/tags/v1.4.1.tar.gz +source_filename = robin-map-1.4.1.tar.gz +source_hash = 0e3f53a377fdcdc5f9fed7a4c0d4f99e82bbb64175233bd13427fef9a771f4a1 +source_fallback_url = https://wrapdb.mesonbuild.com/v2/robin-map_1.4.1-1/get_source/robin-map-1.4.1.tar.gz +patch_filename = robin-map_1.4.1-1_patch.zip +patch_url = https://wrapdb.mesonbuild.com/v2/robin-map_1.4.1-1/get_patch +patch_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/robin-map_1.4.1-1/robin-map_1.4.1-1_patch.zip +patch_hash = 8137d5283004db445d744d603e600119047fed865cbd740f8fec4421dc1a45b6 +wrapdb_version = 1.4.1-1 [provide] -robin-map = robin_map_dep +dependency_names = tsl-robin-map diff --git a/tests/test_callback_exceptions.py b/tests/test_callback_exceptions.py new file mode 100644 index 0000000..960ca84 --- /dev/null +++ b/tests/test_callback_exceptions.py @@ -0,0 +1,155 @@ +"""Exceptions raised in the user-supplied callbacks (model function and Jacobians). + +The callbacks are invoked from Fortran. A C++ exception must never be unwound through +Fortran frames (not guaranteed to work, and fatal without unwind tables), so the extension +stores the exception, aborts the regression through the `istop` flag, and rethrows after +the Fortran routine has returned normally. These tests check the observable consequences. +""" +import traceback + +import numpy as np +import pytest + +from odrpack import OdrStop, odr_fit + +XDATA = np.linspace(0, 2, 20) +YDATA = 1.5*np.exp(0.7*XDATA) +BETA0 = np.array([1., 1.]) + + +def f_ok(x, beta): + return beta[0]*np.exp(beta[1]*x) + + +def jb_ok(x, beta): + return np.stack([np.exp(beta[1]*x), beta[0]*x*np.exp(beta[1]*x)])[:, None, :] + + +def jx_ok(x, beta): + return (beta[0]*beta[1]*np.exp(beta[1]*x))[None, None, :] + + +class Counter: + """Callable wrapper that counts calls and fails from the `after`-th one on. + + The first calls are done by `odr_fit` itself (argument checks), before entering + Fortran, so `after` must be large enough for the failure to happen inside + the iterations. + """ + + def __init__(self, fn, exc=None, after=6, value=None): + self.fn, self.exc, self.after, self.value = fn, exc, after, value + self.calls = 0 + + def __call__(self, x, beta): + self.calls += 1 + if self.calls > self.after: + if self.exc is not None: + raise self.exc + return self.value + return self.fn(x, beta) + + +class CustomError(Exception): + pass + + +def _fit(f=f_ok, **kwargs): + return odr_fit(f, XDATA, YDATA, BETA0, **kwargs) + + +def test_exception_type_message_and_traceback_are_preserved(): + def failing_model(x, beta): + return model(x, beta) + + model = Counter(f_ok, CustomError("custom message")) + with pytest.raises(CustomError, match="custom message") as excinfo: + _fit(failing_model) + tb = "".join(traceback.format_tb(excinfo.value.__traceback__)) + assert "failing_model" in tb # the user's frame is in the traceback + + +@pytest.mark.parametrize("exc", [ValueError("v"), ZeroDivisionError("z"), KeyError("k"), + KeyboardInterrupt(), CustomError("c")]) +def test_various_exception_types_propagate(exc): + with pytest.raises(type(exc)): + _fit(Counter(f_ok, exc)) + + +def test_no_python_callback_after_the_first_failure(): + """The optimiser must not keep calling the user's function once it has failed.""" + model = Counter(f_ok, CustomError("boom")) + with pytest.raises(CustomError): + _fit(model) + assert model.calls == model.after + 1 + + +def test_exception_in_jacobian_beta(): + jb = Counter(jb_ok, CustomError("jb"), after=3) + with pytest.raises(CustomError, match="jb"): + _fit(f_ok, jac_beta=jb, jac_x=jx_ok) + assert jb.calls == jb.after + 1 + + +def test_exception_in_jacobian_x(): + jx = Counter(jx_ok, CustomError("jx"), after=3) + with pytest.raises(CustomError, match="jx"): + _fit(f_ok, jac_beta=jb_ok, jac_x=jx) + assert jx.calls == jx.after + 1 + + +def test_wrong_return_type_in_callback_raises_instead_of_crashing(): + """A conversion error (a C++ exception that is not a python_error) in the callback.""" + with pytest.raises(Exception) as excinfo: + _fit(Counter(f_ok, value="not an array")) + assert not isinstance(excinfo.value, OdrStop) + + +def test_files_are_closed_after_an_exception(tmp_path): + """The Fortran units opened for report/error files must be released. + + Reusing the same file name would fail if the previous unit were still open. + """ + rpt, err = str(tmp_path / "x.rpt"), str(tmp_path / "x.err") + with pytest.raises(CustomError): + _fit(Counter(f_ok, CustomError("boom")), report="short", rptfile=rpt, errfile=err) + sol = _fit(report="short", rptfile=rpt, errfile=err) + assert sol.success + (tmp_path / "x.rpt").unlink() # not locked (matters on Windows) + + +def test_library_is_healthy_after_an_exception(): + ref = _fit() + for exc in (CustomError("a"), ValueError("b")): + with pytest.raises(type(exc)): + _fit(Counter(f_ok, exc)) + sol = _fit() + assert np.array_equal(sol.beta, ref.beta) + assert sol.info == ref.info and sol.nfev == ref.nfev + + +def test_odrstop_is_still_graceful(capfd): + model = Counter(f_ok, OdrStop("stop")) + sol = _fit(model) + assert not sol.success + assert sol.info // 10000 == 5 + + +def test_odrstop_then_other_exception_in_same_fit_is_not_swallowed(): + """`OdrStop` first, anything else later: a later, different exception must still propagate + (if the optimiser calls the function again after a graceful stop).""" + state = {"n": 0} + + def model(x, beta): + state["n"] += 1 + if state["n"] == 7: + raise OdrStop("stop") + if state["n"] > 7: + raise CustomError("later") + return f_ok(x, beta) + + try: + sol = _fit(model) + except CustomError: + return + assert not sol.success diff --git a/tests/test_multiprocessing.py b/tests/test_multiprocessing.py index 071af3c..eab8904 100644 --- a/tests/test_multiprocessing.py +++ b/tests/test_multiprocessing.py @@ -1,7 +1,8 @@ +import multiprocessing import time -from multiprocessing.pool import Pool import numpy as np +import pytest from odrpack import odr_fit @@ -48,7 +49,19 @@ def f3(x: np.ndarray, beta: np.ndarray) -> np.ndarray: case3 = (f3, x3, y3, np.ones_like(beta3)) -def test_multiple_processes(): +@pytest.mark.parametrize("method", ["fork", "spawn", "forkserver"]) +def test_multiple_processes(method): + """Fits run in worker processes must match the serial ones. + + All start methods available on the platform are tested explicitly, because + the default differs between platforms and Python versions (`fork`, `spawn`, + `forkserver`) and each one stresses something different: `fork` copies the + parent's state, `spawn`/`forkserver` import the extension from scratch in + each worker and require everything to be picklable. + """ + if method not in multiprocessing.get_all_start_methods(): + pytest.skip(f"start method '{method}' not available on this platform") + ctx = multiprocessing.get_context(method) # ref solutions sol1 = odr_fit(*case1) @@ -56,12 +69,14 @@ def test_multiple_processes(): sol3 = odr_fit(*case3) # multiple processes - pool = Pool() + pool = ctx.Pool() num_jobs = 10 cases = [case1, case2, case3] - solutions = pool.starmap(odr_fit, cases*num_jobs) - pool.close() - pool.join() + try: + solutions = pool.starmap(odr_fit, cases*num_jobs) + finally: + pool.close() + pool.join() for i in range(0, len(solutions), len(cases)): assert np.allclose(solutions[i].beta, sol1.beta) diff --git a/tests/test_odr_fit.py b/tests/test_odr_fit.py index 5203087..d416285 100644 --- a/tests/test_odr_fit.py +++ b/tests/test_odr_fit.py @@ -658,23 +658,6 @@ def test_OLS(case1): assert np.allclose(sol1.delta, np.zeros_like(sol1.delta)) -def test_exception_odrstop(): - - xdata = np.array([1.0, 2.0, 3.0, 4.0]) - ydata = np.array([1.0, 2.0, 3.0, 4.0]) - beta0 = np.array([1.0, 1.0]) - - def f(x: np.ndarray, beta: np.ndarray) -> np.ndarray: - if beta[0] > 0: - raise OdrStop("Oops!") - return beta[0] * np.exp(beta[1]*x) - - with pytest.raises(OdrStop): - sol = odr_fit(f, xdata, ydata, beta0) - assert not sol.success - assert sol.info == 51000 - - @pytest.mark.xfail(reason="Fails in conda-forge::linux_aarch64") def test_compare_scipy(case1, case2, case3): @@ -709,3 +692,44 @@ def test_compare_scipy(case1, case2, case3): assert np.all(np.max(wd*abs(sol1.delta - sol2[3]['delta']), -1) / (np.max(case['xdata'], -1) - np.min(case['xdata'], -1)) < 1e-4) + + +def _stopping_model(exc, after=5): + """Model function that raises `exc` after `after` calls (the first calls are + the checks done by `odr_fit` itself, before entering the Fortran code).""" + calls = {"n": 0} + + def f(x: np.ndarray, beta: np.ndarray) -> np.ndarray: + calls["n"] += 1 + if calls["n"] > after: + raise exc + return beta[0] * np.exp(beta[1]*x) + + return f + + +def test_odrstop_inside_iterations_stops_gracefully(): + xdata = np.linspace(0, 2, 20) + ydata = 1.5*np.exp(0.7*xdata) + sol = odr_fit(_stopping_model(OdrStop("stop")), xdata, ydata, np.array([1., 1.])) + assert not sol.success + assert sol.info // 10000 == 5 # stopped by user from the model function + + +def test_odrstop_subclass_stops_gracefully(): + class MyStop(OdrStop): + pass + + xdata = np.linspace(0, 2, 20) + ydata = 1.5*np.exp(0.7*xdata) + sol = odr_fit(_stopping_model(MyStop("stop")), xdata, ydata, np.array([1., 1.])) + assert not sol.success + + +def test_other_exceptions_inside_iterations_propagate(): + xdata = np.linspace(0, 2, 20) + ydata = 1.5*np.exp(0.7*xdata) + # The message mentions "OdrStop", but the exception type is not `OdrStop`. + with pytest.raises(ValueError, match="OdrStop"): + odr_fit(_stopping_model(ValueError("this is not an OdrStop")), + xdata, ydata, np.array([1., 1.])) diff --git a/tests/test_threading.py b/tests/test_threading.py new file mode 100644 index 0000000..d4945cb --- /dev/null +++ b/tests/test_threading.py @@ -0,0 +1,429 @@ +"""Thread-safety tests for odrpack. + +These tests are meant to detect problems that only show up when several threads +use the library at the same time, in particular on free-threaded CPython builds +(3.13t, 3.14t, ...). They are also meaningful on regular (GIL) builds: the +model callbacks are Python code (often sleeping), so the GIL is released while +a fit is in progress and other threads can enter the Fortran core. Hidden +global/static state in the Fortran code, or shared mutable state in the wrapper, +would therefore corrupt results there as well. + +What is checked: + +1. The extension module declares itself free-threading-safe, i.e. importing it + does not silently re-enable the GIL (free-threaded builds only). +2. Concurrent fits (heterogeneous problems, sleeping and non-sleeping + callbacks) give bit-identical results to the same fits run serially. +3. Input arrays shared between threads are neither modified nor a source of + cross-talk. +4. Re-entrancy: a fit can be started from inside the model function of another + fit (this fails if the Fortran core keeps any static work storage). +5. Exceptions raised in the callbacks stay confined to the thread that raised + them. +6. Concurrent report/error files (Fortran I/O units) do not get mixed up. +7. The small helper functions of the extension module are thread-safe. +""" +import os +import subprocess +import sys +import sysconfig +import threading +import time +from concurrent.futures import ThreadPoolExecutor + +import numpy as np +import pytest + +from odrpack import OdrStop, odr_fit +from odrpack.__odrpack import loc_iwork, loc_rwork, stop_message, workspace_dimensions + +N_THREADS = 8 +N_ROUNDS = 3 +DELAY = 0.0002 # s; makes the model function release the GIL on regular builds + +FREE_THREADED_BUILD = bool(sysconfig.get_config_var("Py_GIL_DISABLED")) + + +# %% Problem definitions (module-level, deterministic) + +def _make_cases(): + cases = {} + + # m=1, q=1, analytic Jacobians, with weights + def f1(x, beta): + time.sleep(DELAY) + return beta[0] + beta[1] * x + beta[2] * x**2 + beta[3] * x**3 + + def jb1(x, beta): + return np.stack([np.ones_like(x), x, x**2, x**3])[:, None, :] + + def jx1(x, beta): + return (beta[1] + 2 * beta[2] * x + 3 * beta[3] * x**2)[None, None, :] + + x = np.linspace(-10., 10., 21) + beta_true = np.array([1., -2., 0.1, -0.1]) + rng = np.random.default_rng(1) + y = f1(x, beta_true) + rng.normal(0, 0.05, x.shape) + cases["m1q1_jac"] = dict( + f=f1, xdata=x, ydata=y, beta0=np.ones(4), jac_beta=jb1, jac_x=jx1, + weight_x=0.5, weight_y=2.0) + + # m=2, q=1, finite differences, central scheme, bounds + def f2(x, beta): + time.sleep(DELAY) + return (beta[0] * x[0, :])**3 + x[1, :]**beta[1] + + x = np.linspace(-10., 10., 41) + x = np.vstack((x, 10 + x / 2)) + y = f2(x, np.array([2., 2.])) + cases["m2q1_bounds"] = dict( + f=f2, xdata=x, ydata=y, beta0=np.array([1.5, 1.5]), + bounds=(np.array([0.5, 0.5]), np.array([3., 3.])), + diff_scheme="central") + + # m=3, q=2, OLS-free ODR with fixed parameter + def f3(x, beta): + time.sleep(DELAY) + y = np.zeros((2, x.shape[-1])) + y[0, :] = (beta[0] * x[0, :])**3 + x[1, :]**beta[1] + np.exp(x[2, :] / 2) + y[1, :] = (beta[2] * x[0, :])**2 + x[1, :]**beta[1] + return y + + x = np.linspace(-1., 1., 31) + x = np.vstack((x, np.exp(x), x**2)) + y = f3(x, np.array([1., 2., 3.])) + cases["m3q2_fix"] = dict( + f=f3, xdata=x, ydata=y, beta0=np.ones(3), + fix_beta=np.array([False, True, False])) + + # OLS, no sleep (pure compute, maximises true parallelism) + def f4(x, beta): + return beta[0] * np.exp(beta[1] * x) + + x = np.linspace(0., 2., 200) + rng = np.random.default_rng(2) + y = f4(x, np.array([1.5, 0.7])) * (1 + 0.01 * rng.uniform(-1, 1, x.shape)) + cases["ols_nosleep"] = dict( + f=f4, xdata=x, ydata=y, beta0=np.array([1., 1.]), task="OLS") + + # implicit ODR (ellipse), no sleep + def f5(x, beta): + return beta[0] * x[0]**2 + 2 * beta[1] * x[0] * x[1] + beta[2] * x[1]**2 - 1 + + t = np.linspace(0, 2 * np.pi, 40, endpoint=False) + rng = np.random.default_rng(3) + x = np.vstack((2 * np.cos(t) + rng.normal(0, 0.01, t.shape), + 1 * np.sin(t) + rng.normal(0, 0.01, t.shape))) + cases["implicit"] = dict( + f=f5, xdata=x, ydata=np.zeros(t.size), beta0=np.array([0.3, 0.0, 0.8]), + task="implicit-ODR") + + return cases + + +CASES = _make_cases() +RESULT_FIELDS = ("beta", "delta", "eps", "sd_beta", "cov_beta", "xplusd", "yest") +SCALAR_FIELDS = ("info", "nfev", "njev", "niter", "irank", "res_var", + "sum_square", "sum_square_delta", "sum_square_eps") + + +def _run(name, **overrides): + kwargs = {**CASES[name], **overrides} + return odr_fit(**kwargs) + + +def _assert_identical(sol, ref, label=""): + for field in RESULT_FIELDS: + a, b = getattr(sol, field), getattr(ref, field) + assert np.array_equal(a, b, equal_nan=True), f"{label}: `{field}` differs" + for field in SCALAR_FIELDS: + a, b = getattr(sol, field), getattr(ref, field) + assert a == b, f"{label}: `{field}` differs ({a} != {b})" + + +@pytest.fixture(scope="module") +def reference(): + """Serial reference solutions.""" + return {name: _run(name) for name in CASES} + + +def _run_in_threads(tasks, n_threads=N_THREADS): + """Run callables in a pool of `n_threads` threads. + + The first `n_threads` tasks are released simultaneously by a barrier, to + maximise the overlap between fits; the remaining tasks just queue up. + """ + n_start = min(n_threads, len(tasks)) + barrier = threading.Barrier(n_start) + lock = threading.Lock() + counter = [0] + + def wrapper(task): + with lock: + first_batch = counter[0] < n_start + counter[0] += 1 + if first_batch: + barrier.wait(timeout=60) + return task() + + with ThreadPoolExecutor(max_workers=n_threads) as pool: + futures = [pool.submit(wrapper, t) for t in tasks] + return [f.result() for f in futures] + + +# %% 1. Module declares free-threading support + +@pytest.mark.skipif(not FREE_THREADED_BUILD, reason="needs a free-threaded build of CPython") +def test_import_does_not_enable_gil(): + """Importing odrpack must not silently re-enable the GIL. + + If the extension is not compiled as free-threading-compatible (nanobind needs + `NB_FREE_THREADED`), CPython re-enables the GIL at import time and emits a + RuntimeWarning. Everything would still 'work', but not free-threaded, and all + the other tests in this module would be testing nothing. + + Runs in a fresh interpreter, in the default configuration (no `-X gil=0`, + no PYTHON_GIL), which is what users get. + """ + code = ( + "import sys\n" + "assert not sys._is_gil_enabled()\n" + "import odrpack, odrpack.__odrpack\n" + "print('GIL_ENABLED=%s' % sys._is_gil_enabled())\n" + ) + env = {k: v for k, v in os.environ.items() if k != "PYTHON_GIL"} + proc = subprocess.run([sys.executable, "-W", "always", "-c", code], + capture_output=True, text=True, env=env) + assert "GIL_ENABLED=False" in proc.stdout, ( + "odrpack re-enabled the GIL on import.\n" + proc.stdout + proc.stderr) + assert "has not declared that it can run safely without the GIL" not in proc.stderr + + +@pytest.mark.skipif(not FREE_THREADED_BUILD, reason="needs a free-threaded build of CPython") +def test_gil_still_disabled_in_this_process(): + """The GIL must not have been re-enabled by importing/using odrpack here.""" + _run("ols_nosleep") + assert not sys._is_gil_enabled() + + +# %% 2. Concurrent fits == serial fits + +@pytest.mark.parametrize("name", list(CASES)) +def test_same_problem_in_many_threads(name, reference): + """Many threads solving the *same* problem at the same moment.""" + results = _run_in_threads([lambda: _run(name)] * (N_THREADS * 2)) + for i, sol in enumerate(results): + _assert_identical(sol, reference[name], f"{name}[{i}]") + + +def test_different_problems_in_many_threads(reference): + """Heterogeneous problems (different n, m, q, npar, options) at the same time. + + Catches static/global work storage in the Fortran core, which would be + clobbered by problems of different size. + """ + names = list(CASES) * N_ROUNDS + tasks = [(lambda n=n: (n, _run(n))) for n in names] + for name, sol in _run_in_threads(tasks): + _assert_identical(sol, reference[name], name) + + +def test_threads_with_shared_inputs(reference): + """All threads share the very same input arrays (and function objects).""" + before = {k: {a: b.copy() for a, b in c.items() if isinstance(b, np.ndarray)} + for k, c in CASES.items()} + tasks = [(lambda n=n: (n, _run(n))) for n in list(CASES) * N_ROUNDS] + for name, sol in _run_in_threads(tasks): + _assert_identical(sol, reference[name], name) + for k, arrays in before.items(): + for a, b in arrays.items(): + assert np.array_equal(CASES[k][a], b), f"input `{a}` of case `{k}` was modified" + + +def test_results_are_independent_objects(): + """Result arrays returned to different threads must not alias each other.""" + results = _run_in_threads([lambda: _run("ols_nosleep")] * N_THREADS) + for i, a in enumerate(results): + for b in results[i + 1:]: + for field in ("beta", "delta", "eps", "rwork", "iwork"): + assert not np.shares_memory(getattr(a, field), getattr(b, field)) + # mutate one, the others must be unaffected + ref = results[1].beta.copy() + results[0].beta[:] = -1 + assert np.array_equal(results[1].beta, ref) + + +# %% 3. Re-entrancy + +def test_nested_fit_inside_model_function(reference): + """A fit started from inside the model function of another fit. + + Deterministic, single-threaded detector of static storage in the Fortran core. + """ + inner_ref = reference["ols_nosleep"] + base = CASES["m2q1_bounds"] + inner_results = [] + + def f_outer(x, beta): + inner_results.append(_run("ols_nosleep")) + return base["f"](x, beta) + + sol = odr_fit(**{**base, "f": f_outer}) + _assert_identical(sol, reference["m2q1_bounds"], "outer") + assert len(inner_results) > 0 + for s in inner_results: + _assert_identical(s, inner_ref, "inner") + + +def test_nested_fit_inside_model_function_threaded(reference): + def task(): + base = CASES["m3q2_fix"] + inner = [] + + def f_outer(x, beta): + inner.append(_run("m1q1_jac")) + return base["f"](x, beta) + + outer = odr_fit(**{**base, "f": f_outer}) + return outer, inner + + for outer, inner in _run_in_threads([task] * N_THREADS): + _assert_identical(outer, reference["m3q2_fix"], "outer") + for s in inner: + _assert_identical(s, reference["m1q1_jac"], "inner") + + +# %% 4. Exceptions stay in their thread + +class _Boom(Exception): + pass + + +def test_exceptions_do_not_leak_between_threads(reference): + base = CASES["m2q1_bounds"] + + def make_failing(exc, after): + calls = {"n": 0} + + def f(x, beta): + calls["n"] += 1 + if calls["n"] > after: + raise exc + return base["f"](x, beta) + return f + + def task_ok(): + return "ok", _run("m2q1_bounds") + + def task_boom(): + try: + odr_fit(**{**base, "f": make_failing(_Boom("boom"), 5)}) + except _Boom: + return "boom", None + return "no-exception", None + + def task_stop(): + # `OdrStop` raised from within the iteration is handled by the library + try: + sol = odr_fit(**{**base, "f": make_failing(OdrStop("stop"), 5)}) + except OdrStop: + return "odrstop-raised", None + return "odrstop-handled", sol + + kinds = ["ok", "boom", "stop"] + tasks = {"ok": task_ok, "boom": task_boom, "stop": task_stop} + + # serial expectation for the "stop" kind + expected_stop = tasks["stop"]()[0] + + order = kinds * (N_THREADS * 2) + results = _run_in_threads([tasks[k] for k in order]) + for k, (label, sol) in zip(order, results): + if k == "ok": + assert label == "ok" + _assert_identical(sol, reference["m2q1_bounds"], "ok-thread") + elif k == "boom": + assert label == "boom" + else: + assert label == expected_stop + + +# %% 5. Fortran I/O units + +def test_report_and_error_files_do_not_mix(tmp_path): + """Every thread writes its own report/error file: contents must equal serial ones.""" + names = ["m1q1_jac", "m2q1_bounds", "m3q2_fix", "ols_nosleep"] + + def run_with_files(name, tag): + rpt = tmp_path / f"{tag}_{name}.rpt" + err = tmp_path / f"{tag}_{name}.err" + sol = _run(name, report="long", rptfile=str(rpt), errfile=str(err)) + return rpt.read_text(), err.read_text(), sol + + ref = {n: run_with_files(n, "ref") for n in names} + + tasks = [(lambda n=n, i=i: (n, run_with_files(n, f"t{i}"))) + for i, n in enumerate(names * N_THREADS)] + for name, (rpt, err, sol) in _run_in_threads(tasks): + assert rpt == ref[name][0], f"report of `{name}` differs from serial one" + assert err == ref[name][1], f"error file of `{name}` differs from serial one" + assert len(rpt) > 0 + _assert_identical(sol, ref[name][2], name) + + +def test_shared_report_file_name_fails_cleanly(tmp_path): + """Several threads use the same report file name (a user error). + + The Fortran runtime refuses to open one file in two units at once, which + surfaces as a RuntimeError. The only requirements are: no crash/hang, no + other exception type, and the library keeps working afterwards. + """ + path = str(tmp_path / "shared.rpt") + + def task(): + try: + return _run("ols_nosleep", report="short", rptfile=path) + except RuntimeError: + return None + + results = _run_in_threads([task] * N_THREADS) + assert any(r is not None for r in results) + assert all(r is None or r.success for r in results) + # still healthy afterwards + assert _run("ols_nosleep", report="short", rptfile=path).success + + +def test_report_to_stdout_from_many_threads(capfd, reference): + """Concurrent use of the default Fortran output unit must not crash or corrupt results.""" + results = _run_in_threads([lambda: _run("ols_nosleep", report="short")] * N_THREADS) + capfd.readouterr() + for sol in results: + _assert_identical(sol, reference["ols_nosleep"], "stdout-report") + + +# %% 6. Helper functions + +def test_helper_functions_in_threads(): + ref_dims = workspace_dimensions(50, 2, 3, 4, True) + ref_iw = loc_iwork(2, 3, 4) + ref_rw = loc_rwork(50, 2, 3, 4, 1, 1, True) + ref_msg = [stop_message(i) for i in (1, 2, 3, 4, 1000, 10000, 20000, 51000)] + + def task(): + for _ in range(200): + assert workspace_dimensions(50, 2, 3, 4, True) == ref_dims + assert loc_iwork(2, 3, 4) == ref_iw + assert loc_rwork(50, 2, 3, 4, 1, 1, True) == ref_rw + assert [stop_message(i) for i in (1, 2, 3, 4, 1000, 10000, 20000, 51000)] == ref_msg + return True + + assert all(_run_in_threads([task] * N_THREADS)) + + +# %% 7. Long-ish stress (many fits, more threads than cores) + +def test_stress(reference): + names = list(CASES) + tasks = [(lambda n=names[i % len(names)]: (n, _run(n))) for i in range(64)] + for name, sol in _run_in_threads(tasks, n_threads=16): + _assert_identical(sol, reference[name], name)