Debugging Non-Deterministic Output in Parallel Grid Generation

When a grid generator that hashes byte-identical on a single worker starts emitting a different raster on every parallel run, the cause is almost never the algorithm — it is that concurrency has exposed an ordering assumption the serial code silently satisfied.

Part of Async Execution for Large Grids: that workflow covers distributing tile generation across a worker pool; this page isolates the specific regression where the same inputs and the same seed no longer reproduce the same output once the work is spread over more than one worker, and walks through the four mechanisms that cause it and the reproducibility gate that catches all of them.

Root Cause: Concurrency Exposes Four Hidden Ordering Assumptions

Determinism in a generation pipeline means one thing operationally: fixed inputs plus a fixed seed produce a bit-for-bit identical artifact, so a content hash is a stable identity you can gate on. Serial code gets this almost for free because there is exactly one execution order. Parallelizing the same code does not add randomness — it removes the guarantee that operations run in the order the author assumed. Four distinct mechanisms surface once that guarantee is gone.

Per-worker unseeded RNG. The most common cause. A module-level np.random call, or a default_rng() seeded once in the parent process, does not survive a fork or a spawn into worker processes the way the author expects. With spawn (the default on macOS and Windows, and increasingly on Linux), each worker re-imports the module and re-initializes its generator from OS entropy. With fork, workers inherit an identical generator state, so several tiles draw the same stream — deterministic per run but wrong, and different the moment worker count or chunking changes.

Race-order reductions. When per-tile results are combined — summing densities into a global normalizer, concatenating point batches, merging a histogram — the combine step often runs in completion order, not tile order. as_completed, dict insertion from callbacks, and unordered map results all deliver tiles in whatever order the scheduler finished them, which varies run to run.

Floating-point non-associativity. Even when the reduction consumes tiles in a fixed order, floating-point addition is not associative: $(a + b) + c \neq a + (b + c)$ in general, because each intermediate result is rounded. A tree reduction over 64 chunks and a sequential reduction over the same 64 chunks produce sums that differ in the last few ULPs, which is enough to flip a content hash.

$$\left| \sum_{i=1}^{n} x_i - \mathrm{fl}\!\left(\sum_{i=1}^{n} x_i\right) \right| \le \frac{(n-1)\,\varepsilon}{1 - (n-1)\varepsilon} \sum_{i=1}^{n} |x_i|$$

The bound scales with the number of chunks $n$, so re-chunking alone changes the rounding envelope.

Unordered container iteration. Iterating a set of cell ids, or a dict built from concurrent inserts, and feeding that order into any order-sensitive step (seed derivation, point emission, edge assignment) injects run-to-run variation that has nothing to do with the RNG at all.

Minimal Reproducer: Make the Non-Determinism Fail on Demand

Reproduce the defect deterministically before touching it. Pin the toolchain so the scheduler and array semantics match across machines.

Copy
numpy==1.26.*
dask==2024.*        # major pin; distributed scheduler

pythonCopy
import hashlib
import numpy as np
from concurrent.futures import ProcessPoolExecutor

def make_tile(tile_id: int, size: int = 256) -> np.ndarray:
    # BUG: each worker seeds from OS entropy on spawn, so output varies per run.
    rng = np.random.default_rng()
    return rng.random((size, size), dtype=np.float64)

def generate(n_tiles: int = 16, workers: int = 4) -> str:
    with ProcessPoolExecutor(max_workers=workers) as ex:
        # BUG: results arrive in completion order, not tile order.
        tiles = list(ex.map(make_tile, range(n_tiles)))
    grid = np.concatenate(tiles, axis=0)
    return hashlib.sha256(grid.tobytes()).hexdigest()

if __name__ == "__main__":
    print(generate(), generate())   # two different hashes -> non-deterministic

Two runs print two different digests. Drop workers=1 and the RNG bug still fires (each map call still re-seeds), which confirms the seeding defect is independent of the ordering defect — you must fix both.

Fix: Derive Per-Tile Seeds and Reduce in a Fixed Order

The fix has three moving parts: give every tile a reproducible, independent seed derived from a single base seed and the tile id; force any cross-tile reduction to run in tile-id order; and make the reduction itself order-insensitive where you can. Deriving seeds by tile id rather than worker id is what keeps output invariant to worker count and to the chunk size you tune for spatial tiling.

pythonCopy
import hashlib
import numpy as np
from concurrent.futures import ProcessPoolExecutor

def tile_rng(base_seed: int, tile_id: int) -> np.random.Generator:
    """Independent, reproducible stream per tile — invariant to worker count.
    Spawning with the tile id as the spawn key gives a stream that depends only
    on (base_seed, tile_id), never on which worker or in what order it runs."""
    child = np.random.SeedSequence(entropy=base_seed, spawn_key=(tile_id,))
    return np.random.default_rng(child)

def make_tile(args) -> tuple[int, np.ndarray]:
    tile_id, base_seed, size = args
    rng = tile_rng(base_seed, tile_id)
    return tile_id, rng.random((size, size), dtype=np.float64)   # id travels with data

def generate(base_seed: int, n_tiles: int = 16, workers: int = 4, size: int = 256) -> str:
    work = [(t, base_seed, size) for t in range(n_tiles)]
    with ProcessPoolExecutor(max_workers=workers) as ex:
        results = list(ex.map(make_tile, work))
    # Fixed-order combine: sort by tile id, NOT completion order.
    results.sort(key=lambda pair: pair[0])
    grid = np.concatenate([arr for _, arr in results], axis=0)
    return hashlib.sha256(grid.tobytes()).hexdigest()

SeedSequence.spawn is the load-bearing call: it hashes the base entropy with a per-child spawn key, so children are independent regardless of how many workers draw them or in what order. Never seed a worker from its process id, its rank, or the wall clock — all three re-introduce the defect.

For the reduction itself, pin the order and, where the quantity is a count or a sum of small integers, accumulate in integer space so non-associativity cannot bite:

pythonCopy
import numpy as np

def combine_counts(tiles: dict[int, np.ndarray]) -> np.ndarray:
    """Deterministic global histogram: fixed key order + integer accumulation."""
    total = None
    for tile_id in sorted(tiles):                 # sorted -> fixed reduction order
        counts = tiles[tile_id].astype(np.int64)  # integers are associative
        total = counts if total is None else total + counts
    return total

def combine_float_mean(tiles: dict[int, np.ndarray]) -> np.ndarray:
    """When floats are unavoidable, use math.fsum-style pairwise sum in fixed order."""
    ordered = [tiles[t] for t in sorted(tiles)]
    return np.add.reduce(np.stack(ordered), axis=0) / len(ordered)  # deterministic given order

Under Dask, the same principle applies: set optimize_graph aside and control determinism at the reduction. Use da.reduction or .sum(split_every=...) with a fixed split_every so the tree shape is stable across runs, and always name tiles by their spatial index rather than relying on partition arrival order.

Verification Step: Gate the Content Hash

Assert reproducibility across two independent knobs — repeated runs at fixed concurrency, and runs at different worker counts — because each catches a different one of the four mechanisms. Wire this into the same CI flow described in CI/CD Integration for Spatial Data so a regression cannot merge.

pythonCopy
import pytest

def test_parallel_output_is_deterministic():
    seed = 20260623

    # 1. Repeat at fixed concurrency: catches unseeded RNG and unordered iteration.
    a = generate(seed, workers=4)
    b = generate(seed, workers=4)
    assert a == b, "non-deterministic at fixed worker count"

    # 2. Vary concurrency: catches race-order reduction and float non-associativity.
    for w in (1, 2, 8):
        assert generate(seed, workers=w) == a, f"output depends on worker count {w}"

    # 3. A different seed MUST change the output (guards against a frozen constant).
    assert generate(seed + 1, workers=4) != a, "seed has no effect — RNG not wired in"

The third assertion is easy to omit and expensive to lose: without it, a bug that silently produces a constant grid passes the first two checks trivially. All three together certify that output is a pure function of (inputs, base_seed) and nothing else.

Edge Cases & Gotchas

Float reduction across an antimeridian split. When tiles are indexed by longitude and a global normalizer sums per-tile densities, tiles straddling ±180° are often re-ordered by a downstream sort that treats +179° and −179° as far apart. If the reduction order then depends on that sort, the float sum shifts. Index tiles by an integer tile id derived from the grid origin, never by raw longitude, so the antimeridian never participates in reduction ordering.

Fork vs. spawn seed inheritance. On Linux with the default fork start method, workers inherit the parent’s generator state, so unfixed code is deterministic within a run but silently correlated across tiles — the same random field repeats. Switching CI to spawn (or a different container base image) then flips it to per-run randomness, so a pipeline that “passed” locally breaks in CI. Always derive per-tile seeds explicitly; never rely on inherited state under either start method.

GEOS and GDAL thread pools. Geometry operations and raster warps can spin their own thread pools whose scheduling is non-deterministic. A parallel shapely set operation feeding polygon tessellation algorithms may return vertices in thread-completion order. Sort output geometries by a stable key (centroid, then id) before hashing, and pin OMP_NUM_THREADS/GDAL_NUM_THREADS in CI so the thread count itself is reproducible.

Frequently Asked Questions

Related

• Async Execution for Large Grids — the parent workflow for distributing tile generation that this debugging guide troubleshoots.
• Tuning Dask Chunk Size for Spatial Tiling — chunk choices that must not change output once seeds are derived by tile id.
• Polygon Tessellation Algorithms for Synthetic Spatial Data — a downstream stage where thread-order geometry output introduces the same class of non-determinism.
• CI/CD Integration for Spatial Data — where the content-hash reproducibility gate belongs in the pipeline.