From b2b878ba5452619a7de4b880f6293ac0376fdb79 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 6 Sep 2026 15:25:52 -0500 Subject: [PATCH] Add OpenMP to build and enable 5-D multilinear interpolation --- .github/workflows/deploy.yaml | 1 + AGENTS.md | 25 ++++- geotiepoints/multilinear.py | 3 +- geotiepoints/multilinear_cython.pyx | 8 +- geotiepoints/tests/test_multilinear.py | 40 +++++++ setup.py | 149 ++++++++++++++++++++++++- 6 files changed, 218 insertions(+), 8 deletions(-) diff --git a/.github/workflows/deploy.yaml b/.github/workflows/deploy.yaml index 107f4a5..bc687d5 100644 --- a/.github/workflows/deploy.yaml +++ b/.github/workflows/deploy.yaml @@ -57,6 +57,7 @@ jobs: CIBW_SKIP: "cp39-* cp310-* *-manylinux_i686 *-musllinux_i686 *-musllinux_aarch64 *-win32" CIBW_ARCHS: "${{ matrix.arch }}" CIBW_TEST_SKIP: "*_arm64" + CIBW_BEFORE_ALL_MACOS: brew install libomp - uses: actions/upload-artifact@v7 with: diff --git a/AGENTS.md b/AGENTS.md index db1a9c2..609582f 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -115,8 +115,15 @@ graph without computing. `_simple_modis_interpolator` re-acquires the GIL to call `scipy.ndimage.map_coordinates`. - `_modis_utils.xyz2lonlat` deliberately upcasts to float64 internally even for float32 input ("64-bit precision matters apparently"). +- **`multilinear_cython.pyx`'s `prange`/`parallel()` loops are live**, one per kernel (1-D through + 5-D). `setup.py` compiles that extension — and only that extension — with OpenMP; see "Build, test, + lint". Every temporary is assigned inside the `parallel()` block so Cython privatises it per thread, + and the only write is `output[i]` at the loop index, so results are thread-count independent. + Anything assigned inside a `parallel()` block becomes thread-private, which is easy to trip over + when adding a variable that is meant to be shared. - Free-threading support is intentional: `freethreading_compatible=True` in `setup.py`, cp314t wheels, - and a `Free Threading :: 1 - Unstable` classifier. + and a `Free Threading :: 1 - Unstable` classifier. OpenMP is independent of the GIL, so libgomp and + `freethreading_compatible=True` coexist fine. - `_simple_modis_interpolator.pyx` redundantly re-applies the file-level directives as per-function decorators; `_modis_interpolator.pyx` does not. @@ -133,6 +140,17 @@ make -C doc doctest `--cython-coverage` is a **custom `setup.py` flag** that adds the `linetrace`/`profile` directives and `CYTHON_TRACE` macros; without it Cython line coverage is empty. +`multilinear_cython` is built with OpenMP whenever it is available. The `USE_OMP` environment variable +controls this and defaults to `probe`, which picks `/openmp` for MSVC, `-fopenmp`/`-lgomp` for +gcc/conda, and `-Xpreprocessor -fopenmp`/`-lomp` for macOS clang after locating libomp via +`brew ls --verbose libomp` or `port contents libomp`. It also accepts `gcc`/`clang`/`msvc`, and +`USE_OMP=0` forces a serial build — useful for A/B testing and as an escape hatch when a platform +misbehaves. A failed probe degrades to a serial build rather than failing, so watch for the +`Will use ... for OpenMP.` line in build logs. The other three extensions are **deliberately** built +without OpenMP so they don't gain a needless libgomp dependency; `ldd` on the built `.so` files is the +quickest way to confirm this. The macOS wheel job installs libomp via `CIBW_BEFORE_ALL_MACOS` and +cibuildwheel's `delocate` bundles `libomp.dylib` into the wheel. + - Tests load HDF5 fixtures by path relative to the test file (`../../testdata/`), so they only work from a source checkout, never from an installed wheel. - There is **no `conftest.py` and no pytest configuration at all** — no markers, no ini options. @@ -159,9 +177,8 @@ Verified as of this writing; fix them only when the task calls for it. - `_modis_utils.pyx:178` — the error message hardcodes "(10 rows per scan)" regardless of resolution. - `_modis_interpolator.pyx:4` imports `scanline_mapblocks` from `.simple_modis_interpolator` rather than from `._modis_utils` where it is defined, coupling the two MODIS front-ends for no reason. -- `multilinear_cython.pyx` uses `prange`, but `setup.py` compiles with only `-O3` and no - `-fopenmp`/`/openmp`, so those loops are serial. `multilinear_interpolation_5d` is unreachable — the - dispatcher raises for `d > 4`. +- The 5-D multilinear path (`multilinear_interpolation_5d`) was unreachable until recently and is + **experimental**: it now has scipy-comparison coverage in `test_multilinear.py` but no production use. - **Three Earth radii**: `6370997.0` (`__init__.py`, `geointerpolator.py`, `_modis_utils.pyx`), `6370.997` km (`_modis_interpolator.pyx`), `6371008.7714` (`viiinterpolator.py`). - `AbstractMultipleInterpolator.interpolate` (inherited by both `Multiple*Interpolator` classes) diff --git a/geotiepoints/multilinear.py b/geotiepoints/multilinear.py index aedad4e..4db35e8 100644 --- a/geotiepoints/multilinear.py +++ b/geotiepoints/multilinear.py @@ -7,8 +7,7 @@ def mlinspace(smin, smax, orders): if len(orders) == 1: - res = np.atleast_2d( - np.linspace(np.array(smin), np.array(smax), np.array(orders))) + res = np.atleast_2d(np.linspace(smin[0], smax[0], orders[0])) return res.copy() # workaround for strange bug else: meshes = np.meshgrid( diff --git a/geotiepoints/multilinear_cython.pyx b/geotiepoints/multilinear_cython.pyx index c69c624..95a90ed 100644 --- a/geotiepoints/multilinear_cython.pyx +++ b/geotiepoints/multilinear_cython.pyx @@ -26,8 +26,10 @@ def multilinear_interpolation(floating[:] smin, floating[:] smax, long[:] orders cdef floating[:] vals cdef floating[:] res - if d > 4: + if d > 5: raise Exception("Can't interpolate in dimension strictly greater than 5") + if d < 1: + raise Exception("Can't interpolate in dimension less than 1") with nogil: for i in range(n_v): @@ -41,6 +43,10 @@ def multilinear_interpolation(floating[:] smin, floating[:] smax, long[:] orders multilinear_interpolation_3d(smin, smax, orders, vals, n_s, s, res) elif d == 4: multilinear_interpolation_4d(smin, smax, orders, vals, n_s, s, res) + elif d == 5: + # EXPERIMENTAL: this kernel was unreachable until now (the dispatcher + # raised for d > 4) and has had no coverage. See test_multilinear.py. + multilinear_interpolation_5d(smin, smax, orders, vals, n_s, s, res) return result_arr diff --git a/geotiepoints/tests/test_multilinear.py b/geotiepoints/tests/test_multilinear.py index f8b0a8f..ffc8483 100644 --- a/geotiepoints/tests/test_multilinear.py +++ b/geotiepoints/tests/test_multilinear.py @@ -2,6 +2,8 @@ import unittest import numpy as np +import pytest +from scipy.interpolate import RegularGridInterpolator from geotiepoints.multilinear import MultilinearInterpolator @@ -72,6 +74,44 @@ def tearDown(self): return +def _reference_function(x): + """Smooth function of the ``d`` coordinates in the rows of ``x``.""" + return np.sum(np.sin(x + np.arange(x.shape[0])[:, None]), axis=0) + + +@pytest.mark.parametrize("d", [1, 2, 3, 4, 5]) +def test_multilinear_matches_scipy(d): + """Compare each dimensionality of the Cython kernels against scipy.""" + smin = [-1.0] * d + smax = [1.0] * d + orders = [4] * d + + interp = MultilinearInterpolator(smin, smax, orders) + interp.set_values(np.atleast_2d(_reference_function(interp.grid))) + + axes = [np.linspace(smin[i], smax[i], orders[i]) for i in range(d)] + scipy_interp = RegularGridInterpolator( + axes, interp.values[0].reshape(orders), method="linear") + + rng = np.random.default_rng(1234) + # in-bounds points only: outside the grid the Cython kernels extrapolate + # linearly, which is an intentional difference from scipy's handling + points = rng.uniform(-1.0, 1.0, (d, 500)) + + result = interp(points) + expected = scipy_interp(points.T) + + np.testing.assert_allclose(result[0], expected, rtol=1e-12, atol=1e-12) + + +def test_multilinear_too_many_dimensions(): + """A 6-D input is rejected with a message that matches the guard.""" + interp = MultilinearInterpolator([-1.0] * 6, [1.0] * 6, [4] * 6) + interp.set_values(np.zeros((1, 4 ** 6))) + with pytest.raises(Exception, match="strictly greater than 5"): + interp(np.zeros((6, 3))) + + def suite(): """The suite for Multilinear Interpolator""" loader = unittest.TestLoader() diff --git a/setup.py b/setup.py index f965f26..aa459ba 100644 --- a/setup.py +++ b/setup.py @@ -1,5 +1,7 @@ """Setting up the geo_interpolator project.""" +import os +import re import sys from setuptools import setup, find_packages @@ -16,6 +18,33 @@ else: extra_compile_args = ["-O3"] +# Only this extension uses ``prange``; the MODIS extensions are deliberately +# built without OpenMP so they don't gain a needless libgomp dependency. +OMP_EXTENSION = "geotiepoints.multilinear_cython" + +OMP_SETTING_TABLE = { + '1': 'probe', + '0': None, + 'gcc': 'gomp', + 'gomp': 'gomp', + 'clang': 'omp', + 'omp': 'omp', + 'msvc': 'msvc', + 'probe': 'probe', +} + +OMP_COMPILE_ARGS = { + 'gomp': ['-fopenmp'], + 'omp': ['-Xpreprocessor', '-fopenmp'], + 'msvc': ['/openmp'], +} + +OMP_LINK_ARGS = { + 'gomp': ['-lgomp'], + 'omp': ['-lomp'], + 'msvc': [], +} + EXTENSIONS = [ Extension( 'geotiepoints.multilinear_cython', @@ -70,7 +99,125 @@ ext.define_macros = define_macros ext.cython_directives.update(cython_directives) -cmdclass = versioneer.get_cmdclass(cmdclass={"build_ext": build_ext}) + +class build_ext_subclass(build_ext): + """Add OpenMP flags to the one extension that uses ``prange``.""" + + def build_extensions(self): + omp_compile_args, omp_link_args = _omp_compile_link_args(self.compiler.compiler_type) + for ext in self.extensions: + if ext.name != OMP_EXTENSION: + continue + ext.extra_compile_args = list(ext.extra_compile_args or []) + omp_compile_args + ext.extra_link_args = list(ext.extra_link_args or []) + omp_link_args + build_ext.build_extensions(self) + + +def _omp_compile_link_args(compiler): + """Get the OpenMP compile and link arguments for this compiler and platform.""" + try: + use_omp = OMP_SETTING_TABLE[os.environ.get('USE_OMP', 'probe')] + except KeyError: + raise ValueError("Unknown USE_OMP value %r, expected one of: %s" + % (os.environ.get('USE_OMP'), ", ".join(sorted(OMP_SETTING_TABLE)))) + + compile_args = [] + link_args = [] + if use_omp == "probe": + use_omp, compile_args, link_args = _probe_omp_for_compiler_and_platform(compiler) + + print(f"Will use {use_omp} for OpenMP." if use_omp else "OpenMP support not available.") + compile_args = compile_args + OMP_COMPILE_ARGS.get(use_omp, []) + link_args = link_args + OMP_LINK_ARGS.get(use_omp, []) + print(f"Compiler: {compiler} / OpenMP: {use_omp} / " + f"OpenMP compile args: {compile_args} / OpenMP link args: {link_args}") + return compile_args, link_args + + +def _probe_omp_for_compiler_and_platform(compiler): + compile_args = [] + link_args = [] + if compiler == "msvc": + use_omp = "msvc" + elif _is_conda_interpreter(): + # Conda provides its own compiler which does support openmp + use_omp = "gomp" + elif _is_macOS(): + # OpenMP is not supported with system clang but homebrew and macports have libomp packages + compile_args, link_args = _macOS_omp_options_from_probe() + if not (compile_args or link_args): + print("Probe for libomp failed, skipping use of OpenMP with clang.") + print("It may be possible to build with OpenMP using USE_OMP=clang with CFLAGS and " + "LDFLAGS explicit settings to use libomp.") + use_omp = None + else: + use_omp = "omp" + else: + use_omp = "gomp" + return use_omp, compile_args, link_args + + +def _is_conda_interpreter(): + """Is the running interpreter from Anaconda, miniconda, or conda-forge? + + Modern conda-forge builds don't always mention conda in ``sys.version``, so + the environment variable is checked first. + + """ + if os.environ.get("CONDA_PREFIX"): + return True + return 'conda' in sys.version or 'Continuum' in sys.version + + +def _is_macOS(): + return 'darwin' in sys.platform + + +def _macOS_omp_options_from_probe(): + """Get common include and library paths for libomp installation on macOS. + + For example ``(['-I/opt/local/include/libomp'], ['-L/opt/local/lib/libomp'])``. + + """ + for cmd in ["brew ls --verbose libomp", "port contents libomp"]: + inc, lib = _compile_link_paths_from_manifest(cmd) + if inc and lib: + return [f"-I{inc}"], [f"-L{lib}"] + return [], [] + + +def _compile_link_paths_from_manifest(cmd): + """Parse include and library paths from macOS package managers. + + Example executions:: + + # Homebrew + $ brew ls --verbose libomp + /opt/homebrew/Cellar/libomp/15.0.7/include/omp.h + /opt/homebrew/Cellar/libomp/15.0.7/lib/libomp.dylib + + # MacPorts + $ port contents libomp + Port libomp contains: + /opt/local/include/libomp/omp.h + /opt/local/lib/libomp/libomp.dylib + + """ + from subprocess import run + query = run(cmd, shell=True, check=False, capture_output=True) + if query.returncode != 0: + return None, None + manifest = query.stdout.decode("UTF-8") + # find all the unique directories mentioned in the manifest + dirs = set(os.path.split(filename)[0] for filename in re.findall(r'^\s*(/.*?)\s*$', manifest, re.MULTILINE)) + # find a unique libdir and incdir + inc = tuple(d for d in dirs if re.search(r'/include(\W|$)', d)) + lib = tuple(d for d in dirs if re.search(r'/lib(\W|$)', d)) + # only return success if there's no ambiguity + return (inc + lib) if len(inc) == 1 and len(lib) == 1 else (None, None) + + +cmdclass = versioneer.get_cmdclass(cmdclass={"build_ext": build_ext_subclass}) with open('README.md', 'r') as readme: README = readme.read()