perf: release the GIL during s2 boolean operations

[thodson-usgs]

⚠️ I am not a CPython expert. Claude identified this performance improvement. I iterated with it for an hour to develop this PR, but before proceeding further, I'm seeking feedback from someone with more experience.

Summary

Wrap s2geog::s2_boolean_operation inside BooleanOp::operator() in a py::gil_scoped_release. Callers can then parallelize boolean ops across Python threads; index extraction stays under the GIL and only the pure-C++ s2 work runs with the GIL released.

 PyObjectGeography operator()(PyObjectGeography a, PyObjectGeography b) const {
     const auto& a_index = a.as_geog_ptr()->geog_index();
     const auto& b_index = b.as_geog_ptr()->geog_index();
-    std::unique_ptr<s2geog::Geography> geog_out =
-        s2geog::s2_boolean_operation(a_index, b_index, m_op_type, m_options);
+    std::unique_ptr<s2geog::Geography> geog_out;
+    {
+        py::gil_scoped_release release;
+        geog_out = s2geog::s2_boolean_operation(a_index, b_index, m_op_type, m_options);
+    }
     return make_py_geography(std::move(geog_out));
 }

Why this is safe. The released-GIL block only reads the two const S2ShapeIndex& references, and MutableS2ShapeIndex documents this pattern as thread-safe in s2/mutable_s2shape_index.h:

For efficiency, updates are batched together and applied lazily on the first subsequent query. Locking is used to ensure that MutableS2ShapeIndex has the same thread-safety properties as "vector": const methods are thread-safe, while non-const methods are not thread-safe.

That covers both the pre-built case and the lazy-first-build case; in the latter, s2's internal wait-mutex serializes the one-time build across threads.

Motivation

In a downstream gridding experiment, spherely.intersection over a large object array dominated the total wall time. py::vectorize holds the GIL throughout the batched call, so chunking the array across a ThreadPoolExecutor gave no speedup on its own. Releasing the GIL around the s2 work lets a batched call fan out across threads — in that experiment, end-to-end setup time dropped by about half.

Measured speedups

macOS arm64, 12-core M-series (6 P-cores + 6 E-cores), Python 3.13. Median of 3 runs of the benchmark script below.

threads	main	patched	speedup
1	995 ms	1032 ms	—
2	1038 ms	575 ms	1.80×
4	997 ms	395 ms	2.61×
8	1017 ms	628 ms	1.64×

Unpatched main does not scale with threads, as expected; single-threaded runtime is unchanged within noise.

Benchmark

Paste into a file and run against the installed spherely. No dependencies beyond spherely + numpy.

"""Minimal benchmark for PR #115 — threaded scaling of spherely.intersection."""

import time
from concurrent.futures import ThreadPoolExecutor

import numpy as np
import spherely


def grid_cells(ny, nx, lat=(-80.0, 80.0), lon=(-175.0, 175.0)):
    lats = np.linspace(lat[0], lat[1], ny + 1)
    lons = np.linspace(lon[0], lon[1], nx + 1)
    out = np.empty(ny * nx, dtype=object)
    for j in range(ny):
        for i in range(nx):
            shell = [
                (lons[i], lats[j]),
                (lons[i + 1], lats[j]),
                (lons[i + 1], lats[j + 1]),
                (lons[i], lats[j + 1]),
            ]
            out[j * nx + i] = spherely.create_polygon(shell, oriented=True)
    return out


# Two same-resolution grids, slightly offset so every element-wise pair
# overlaps. Replicated 25x for ~1 s serial runtime.
tgt = grid_cells(50, 100)
src = grid_cells(50, 100, lat=(-78.5, 81.5), lon=(-173.0, 177.0))
dst = np.tile(tgt, 25)
src_pairs = np.tile(src, 25)
n = len(dst)


def serial():
    spherely.intersection(dst, src_pairs)


def threaded(workers):
    splits = np.array_split(np.arange(n), workers)
    with ThreadPoolExecutor(max_workers=workers) as pool:
        list(pool.map(lambda s: spherely.intersection(dst[s], src_pairs[s]), splits))


def median_ms(fn, trials=3):
    fn()  # warm
    runs = []
    for _ in range(trials):
        t0 = time.perf_counter()
        fn()
        runs.append(time.perf_counter() - t0)
    return 1000 * sorted(runs)[len(runs) // 2]


print(f"1 thread  : {median_ms(serial):6.0f} ms")
for w in (2, 4, 8):
    print(f"{w} threads : {median_ms(lambda w=w: threaded(w)):6.0f} ms")

Tests

tests/test_boolean_operations_concurrency.py adds two tests (5 items with parametrization) that target the failure modes specific to this change:

• test_concurrent_shared_inputs_match_serial, parametrized over the four boolean ops, exercises concurrent reads on shared index state. Each op runs on a 4-thread ThreadPoolExecutor behind a threading.Barrier so all threads enter the released-GIL window simultaneously; outputs are compared bit-for-bit against a serial reference.
• test_concurrent_lazy_index_race does the same but on freshly-constructed Geographies, to surface any first-access materialization inside the index object that runs with the GIL released.

Open questions for the maintainer

1. TSAN. I built spherely with -fsanitize=thread against the conda-forge s2geography / s2geometry and ran the concurrency tests; no warnings. Caveat: TSAN only instruments what was compiled with the flag, and conda-forge's s2geography is not — so the wrapper is covered but the index internals aren't. Wiring TSAN into CI would need either a TSAN-built s2geography variant (doesn't exist on conda-forge) or building s2geography from source in the CI job. Left out of CI for now; happy to add it gated on merge if you'd like.
2. Scope. Only BooleanOp is patched. Predicates (intersects, contains, within, touches) go through the same index-based code path and would benefit identically; happy to expand in this PR or leave for a follow-up.

Appendix: performance canary

Not included in the PR — this test ran a 4000-pair overlapping-polygon workload across a ThreadPoolExecutor to directly confirm the GIL was released. It's CPU-intensive and timing-sensitive (reliable on workstation hardware, flaky on shared CI runners). Reproduced here for anyone who wants to verify the release locally.

@pytest.mark.slow
def test_gil_release_actually_enables_parallelism():
    """Directly confirms the GIL was released: if it wasn't, the threaded
    path serializes on the GIL and runs no faster than a single thread."""
    if N_THREADS < 2:
        pytest.skip("need at least 2 cores for a meaningful speedup measurement")

    big1 = spherely.from_wkt("POLYGON ((-80 -40, 80 -40, 80 40, -80 40, -80 -40))")
    big2 = spherely.from_wkt("POLYGON ((-40 -80, 40 -80, 40 80, -40 80, -40 -80))")
    n = 4000
    dst = np.array([big1] * n, dtype=object)
    src = np.array([big2] * n, dtype=object)

    def serial():
        spherely.intersection(dst, src)

    def threaded():
        splits = np.array_split(np.arange(n), N_THREADS)

        def worker(idx):
            spherely.intersection(dst[idx], src[idx])

        with ThreadPoolExecutor(max_workers=N_THREADS) as pool:
            list(pool.map(worker, splits))

    def median_of(fn, trials=3):
        serial()  # warm any caches
        samples = []
        for _ in range(trials):
            t0 = time.perf_counter()
            fn()
            samples.append(time.perf_counter() - t0)
        return sorted(samples)[len(samples) // 2]

    t_serial = median_of(serial)
    t_threaded = median_of(threaded)
    # Expect at least 1.3× on 4 cores. Well below the theoretical 4× to
    # tolerate bench noise, but any number <1.1 means the release almost
    # certainly did not land.
    assert t_threaded < t_serial / 1.3, (
        f"threaded={t_threaded * 1000:.1f}ms was not meaningfully faster than "
        f"serial={t_serial * 1000:.1f}ms — GIL release may not be in effect"
    )

On the 12-core M-series machine used for the speedup table above it passes with threaded ≈ 2.6× serial; on unpatched main it fails with threaded ≈ serial within noise.