Add OpenMP to build and enable 5-D multilinear interpolation

.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:

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)

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(

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
 

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()

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()