haloexec

Cellular Automaton execution engine using Domain Decomposition with Halo Zones (Ghost Cell Pattern), integrated with real dissmodel via a pip dependency — without living inside the dissmodel core repository, which is currently under JOSS review (issue #10827).

Note

Full documentation: see the docs/ directory for the complete reference — Index, Theory & Core Concepts, Architecture & Design Patterns, API Reference, Tutorials & Recipes.

Motivation

The block+halo pattern is not specific to any single domain model or even to cellular automata in general — it is generic domain decomposition for any transition rule with a local neighborhood dependency, applied to a grid too large to fit entirely in RAM.

dissmodel core (0.6.3, PyPI) already ships RasterCellularAutomaton (dissmodel.geo.raster.cellular_automaton) with the rule(arrays: dict) -> dict contract, but it processes the whole grid at once — there is no chunking. haloexec fills that gap.

Design: same rule() contract

HaloChunkedRasterCellularAutomaton extends RasterCellularAutomaton and preserves exactly the same rule() contract as the base class. Any rule already written for dissmodel runs in blocks+halo by swapping only the base class — zero change to the model's scientific logic:

import numpy as np
from dissmodel.core import Environment
from dissmodel.geo.raster.backend import RasterBackend
from haloexec import HaloChunkedRasterCellularAutomaton
from dissmodel.visualization.raster_map import RasterMap

class GameOfLife(HaloChunkedRasterCellularAutomaton):
    def rule(self, arrays):
        state = arrays["state"]
        neighbors = self.backend.focal_sum_mask(state == 1)
        born = (state == 0) & (neighbors == 3)
        survive = (state == 1) & np.isin(neighbors, [2, 3])
        return {"state": np.where(born | survive, 1, 0).astype(np.uint8)}

backend = RasterBackend(shape=(200, 200))
backend.set("state", np.random.randint(0, 2, (200, 200)).astype(np.uint8))

env = Environment(start_time=1, end_time=50)
GameOfLife(backend=backend, block_h=50, block_w=50, halo=1)

RasterMap(
    backend=backend,
    band="state",
    color_map={0: "#ffffff", 1: "#2f8f6e"},
    labels={0: "morta", 1: "viva"},
    title=f"Padrões clássicos sobre fronteiras de bloco",
)
env.run()

The same GameOfLife class, inheriting from plain RasterCellularAutomaton (no block_h/block_w/halo), runs monolithically with no change to rule() — that property is exactly what the equivalence tests verify.

How it works internally

HaloChunkedRasterCellularAutomaton.execute():

1. Snapshots the global grid and applies np.pad (global halo at the domain's outer boundary).
2. For each block (Block/make_blocks, pure partitioning, no dissmodel dependency), builds a temporary RasterBackend holding only that block's sub-grid+halo.
3. Swaps self.backend for the block's temporary backend — this guarantees that internal rule calls such as self.backend.focal_sum_mask(...) operate on the correct local shape (focal_sum_mask uses self.shape of the active backend).
4. Calls self.rule(block_backend.snapshot()) — same signature as always.
5. Trims the halo off the result and writes it into the corresponding position of the new global grid.
6. Restores self.backend to the real global backend.

Theoretical foundation

• Engineering pattern: Kjolstad, F. B.; Snir, M. Ghost Cell Pattern. In: Proceedings of the 2010 Workshop on Parallel Programming Patterns (ParaPLoP '10). ACM, 2010. https://doi.org/10.1145/1953611.1953615
• Direct application to geospatial LULC-CA: Xia, W. et al. Dynamic Load Balancing Based on Hypergraph Partitioning for Parallel Geospatial Cellular Automata Models. ISPRS Int. J. Geo-Inf., 14(3):109, 2025. https://doi.org/10.3390/ijgi14030109

Installation (development)

pip install -e ".[dev]"
pytest

Depends on dissmodel>=0.6.3 (PyPI) as a real runtime dependency.

Finding: boundary_value=0 is not safe for every domain

Not every domain can safely default the halo's outer boundary fill to 0. If 0 is a valid class code for some array (say, a land-use class or a channel type — not a "no data" sentinel), filling the domain's outer edge with 0 creates phantom migration sources at that edge that don't exist in the monolithic run. This surfaced in testing even with a single block covering the entire grid — i.e., it was never a block-to-block boundary bug, it was an outer-domain-edge fill bug.

Fix: boundary_value at every layer now accepts a dict {name: value} in addition to a scalar, via engine.resolve_boundary_value(). Names ending in "_past" automatically fall back to the base name's value if they have no entry of their own — without that fallback, the bug reappeared in a subtler form (the "_past" array generated by synchronize() kept using 0 even with {"class_code": -1} configured, since the key being looked up was literally "class_code_past").

boundary_value = {"land_use": 0, "elevation": -9999.0, "soil_class": -1}
MyModelHalo(backend=backend, block_h=10, block_w=10, halo=1,
            boundary_value=boundary_value)

Whenever a domain uses 0 as a valid class code in some array, boundary_value must be a dict aligned to that array's real nodata value — never rely on the scalar default of 0.

Finding: the correct halo depth is the dependency chain's depth, not the rule's nominal shift radius

The most consequential lesson from validating this engine against a real (not synthetic) dataset: halo depth is not the same thing as neighborhood radius.

A rule that only reads a raw neighbor value (e.g., "is this cell adjacent to water?") has a 1-hop dependency — halo=1 is enough. But a rule that reads a derived quantity computed from neighbors (e.g., "how many of my neighbors have a lower elevation, and what do their neighbors look like") has a 2-hop dependency, not 1 — with halo=1, the derived quantity computed inside the block's own halo ring is already wrong (its own 2-hop neighbors were zero-padded by the local shift), and that error contaminates the block's core on write-back.

In a real hydrological flow model validated this way, this showed up as dozens of divergent cells in the first time step, growing into the thousands by step 19 — even after the boundary_value fix above was already applied. The root cause was exactly this: the model computed a "how many lower neighbors do I have" quantity and then used the neighbor's version of that same quantity to derive flow — a 2-hop read hiding behind what looked like a 1-hop rule.

Why synthetic tests can miss this entirely: if a synthetic test always places its source condition in a full column or row at the domain's outer edge, the zero-padding artifact from halo=1 already exists equally on both sides (monolithic and block+halo) at that edge — the bug only appears when a source condition sits near an internal block boundary, which typically only happens with real, irregularly-shaped domains.

Fix: there is no code fix for this — halo depth is the responsibility of whoever configures the model, not something the engine can infer on its own from the rule's source code. What changed was making this dependency explicit in documentation, plus a regression test that pins halo=1 as insufficient and halo=2 as correct for that specific model.

General takeaway: the correct halo is not the nominal radius of the shift operation used inside the rule — it is the real depth of the dependency chain. If a rule reads a quantity derived from neighbors (rather than a neighbor's raw value), the dependency is N+1 hops, not N. Synthetic tests with sources only at the domain's outer boundary can mask this; validation against real, irregular data is what tends to expose it.

Two input paths: GeoTIFF/VRT or Zarr

geotiff_io.py/mosaic_io.py (TIFF/VRT, via rasterio) and zarr_io.py (via zarr) are interchangeable input paths — both populate the same MemmapRasterWorkspace, block by block, without materializing the full array in RAM. Swap the loader; the rest of the pipeline (halo, disk, models) doesn't change.

• GeoTIFF/VRT: for data already in traditional raster form, or for tile mosaics (see "Tile mosaics" below).
• Zarr: designed to consume DerivedVariable directly from disscube — which already stores data natively in Zarr, already aligned to the master grid via its GridAligner (per-operator resampling and fine alignment for categorical variables — more rigorous than doing this by hand). Supports both a group with multiple variables (variable_map optional if names differ) and a single array, plus variables with a time dimension (time_index, disscube's "Temporal Backend").

from haloexec import MemmapRasterWorkspace, load_zarr_into_workspace

ws = MemmapRasterWorkspace.create(root="workspace", shape=(height, width),
                                   arrays={"land_use": np.int16, "elevation": np.float32},
                                   block_h=512, block_w=512, halo=1)
load_zarr_into_workspace(ws, "data/derived/BDC_100m/009002/abc123/",
                          variable_map={"land_use": "land_use", "elevation": "elevation"})

Requires the optional zarr extra (pip install -e ".[zarr]").

Tile mosaics (multiple satellite files)

Handled by a separate package: geomosaic — tile discovery, grid contract validation, and VRT construction, with no dependency on haloexec. Mosaicking strictly precedes any execution engine; haloexec doesn't know (and doesn't need to know) that a mosaic exists underneath — a VRT opens via rasterio.open() exactly like a plain GeoTIFF, including for blocks that straddle a tile boundary.

from geomosaic import discover_tiles, build_mosaic_contract, write_vrt
from haloexec import MemmapRasterWorkspace, load_geotiff_into_workspace

tiles = discover_tiles("data/mapbiomas_tiles/")
contract = build_mosaic_contract(tiles)
vrt_path = write_vrt(contract, "data/mosaic.vrt")

ws = MemmapRasterWorkspace.create(root="workspace", shape=(contract.mosaic_height, contract.mosaic_width),
                                   arrays={"land_use": np.int16}, block_h=512, block_w=512, halo=1)
load_geotiff_into_workspace(ws, vrt_path, [("land_use", "int16", 0)])

tests/test_geomosaic_integration.py proves this integration (requires the optional geomosaic extra, used only in the test — never imported by haloexec's runtime code).

Loading GeoTIFF straight to disk

geotiff_io.py (load_geotiff_into_workspace) loads a real GeoTIFF block by block directly into a MemmapRasterWorkspace, via rasterio.windows.Window — never materializing a whole band in RAM. It uses the band_spec convention already established by dissmodel.io.raster.load_geotiff (a list of (name, dtype, nodata)) instead of hardcoded band names.

from haloexec import MemmapRasterWorkspace, load_geotiff_into_workspace

ws = MemmapRasterWorkspace.create(
    root="/data/workspace", shape=(height, width),
    arrays={"land_use": np.int16, "elevation": np.float32, "soil_class": np.int16, "mask": np.uint8},
    block_h=512, block_w=512, halo=2,
)
load_geotiff_into_workspace(ws, "domain.tif", BAND_SPEC + [("mask", "uint8", 0)])

Requires the optional geotiff extra (pip install -e ".[geotiff]"). Validated (tests/test_geotiff_io_equivalence.py) with an exact round-trip against direct rasterio reads, including irregular blocks and a block larger than the grid.

Iterative convergence (unbounded spatial dependency)

convergence.py (sweep_until_convergence) solves a DIFFERENT problem from what a fixed-size halo solves: connectivity, flow routing, watershed delineation — where a cell's value can, in principle, depend on the entire domain, not just its immediate neighbors.

Generalized from a domain-specific student prototype (tidal connectivity via scipy.ndimage.binary_propagation), though that prototype's specific rule is not part of this primitive — only the orchestration pattern was extracted: a small halo plus repeated global sweeps until no block changes, instead of one large halo. Each sweep writes its result immediately back (Gauss-Seidel, via write_block_core_in_place — no ping-pong), so a block processed later in the same sweep already sees the update from a block processed earlier, accelerating convergence.

from haloexec import MemmapRasterWorkspace, sweep_until_convergence
from scipy.ndimage import binary_propagation

def connectivity_rule(window, halo=1):
    connected = window["connected"].astype(bool)
    permeable = window["permeable"].astype(bool)
    propagated = binary_propagation(connected, mask=permeable)
    return {"connected": propagated[halo:-halo, halo:-halo].astype(np.uint8)}

info = sweep_until_convergence(ws, connectivity_rule, boundary_value=0)
# info = {"sweeps": N, "blocks_changed_total": M, "converged": True}

Validated against what the original prototype lacked: tests/test_convergence.py proves the blocks+sweeps version converges to a result exactly identical to a monolithic binary_propagation run on the whole domain at once — across 4 block configurations (including a permeability maze forcing connectivity to cross several block boundaries before converging) plus 5 stress seeds and a non-convergence case (must raise RuntimeError, not hang silently). The original prototype never had this kind of proof.

Educational examples (examples/)

• examples/gol_patterns/gol_patterns_haloexec.py — Game of Life with classic patterns (glider, blinker, beacon, toad, block, pulsar) deliberately positioned across block boundaries, via dissmodel_ca (PATTERNS) and dissmodel.visualization.RasterMap. tests/test_gol_patterns_example.py proves monolithic-vs-blocks equivalence for the same scenario (found in the process: the original coordinates had beacon overlapping pulsar — place() silently overwrites, with no error; fixed).

Equivalence tests

tests/test_gameoflife_from_geotiff.py closes the loop mosaic→TIFF→disk→halo with a simple model: Game of Life loaded from a GeoTIFF (simulating an already-materialized mosaic) straight into MemmapRasterWorkspace, including a "large" case (2000×2000 = 4 million cells). Before this test, TIFF was only validated by round-trip (without running any model on top of it), and Game of Life was only tested with synthetic data in RAM/disk — never the two together.

tests/test_equivalence.py proves, using real dissmodel Environment/RasterBackend (not an isolated harness), that the result of GameOfLifeHalo (blocks+halo) is cell-for-cell identical to GameOfLifeMono (monolithic) after N time steps, across 9 distinct configurations (evenly divisible grid, grid with remainder, 1-row blocks, a block larger than the grid, and a stress case with 5 random seeds).

Disk layer (grids larger than RAM)

disk_backend.py (MemmapRasterWorkspace) and disk_sync_model.py (DiskChunkedSyncRasterModel) generalize the same domain decomposition for grids that don't fit in memory, using np.memmap with double-buffering and checkpointing. No dependency on dissmodel in disk_backend.py — reusable by any framework.

Extracted and generalized from a student prototype (a domain-specific preprocessing pipeline) that already had correct memmap+double-buffer +checkpoint mechanics, but tied to domain-specific state names. Here arrays are named generically (a name -> dtype dict), with no coupling to any particular domain's variable names.

from haloexec import DiskChunkedSyncRasterModel, MemmapRasterWorkspace, workspace_arrays_for_sync_model

class MyModelDiskHalo(DiskChunkedSyncRasterModel, MyModel):
    pass

arrays = workspace_arrays_for_sync_model(
    base={"land_use": np.int16, "elevation": np.float32},
    land_use_types=["land_use", "elevation"],
)
ws = MemmapRasterWorkspace.create(root="/data/workspace", shape=(50000, 50000),
                                   arrays=arrays, block_h=512, block_w=512, halo=1)
ws.fill("land_use", initial_land_use)
ws.fill("elevation", initial_elevation)

env = Environment(start_time=1, end_time=100)
MyModelDiskHalo(workspace=ws)
env.run()

Central technical point: SyncRasterModel.synchronize() (which generates the "<name>_past" arrays) does a .copy() of the whole array — applied naively to a memmap, that would materialize the whole grid in RAM. DiskChunkedSyncRasterModel replicates that logic locally (without importing or modifying installed dissmodel), copying block by block. This is the piece to reconcile when this package eventually migrates into core.

Sparsity and disk cost (memmap solves RAM, not disk)

Worth separating two costs that are easy to conflate:

• RAM — solved by np.memmap: the kernel brings in pages on demand, so iterating over the whole grid never materializes it all at once. Measured with real data (a 371M-cell grid, two float32 arrays): peak RSS of 3.0 GB for a 5.6 GB workspace.
• Disk — memmap doesn't help here, nor should it: a file is a file.

For disk, create() sizes the .dat files without pre-writing zeros, so they start out sparse: only what gets written occupies disk blocks. Reading a never-written region still returns zero — a POSIX guarantee, identical to what pre-writing zeros gave, so the semantics haven't changed. Measured: a 4000×4000 float64 array goes from 122 MB of real disk usage to 0 MB until something is written.

The limit of this saving, and why it isn't automatic: it only shows up if the loader leaves blocks unwritten. A loader that fills every block — including empty ones, with a sentinel like NaN — makes the file dense again. And leaving a region unwritten means that region reads as zero, not "absent". For domains where 0 is a valid value (a class code, sea-level elevation) the two cases become indistinguishable.

Proposed next step (not implemented): distinguishing "absent" from "valid zero" needs one more channel that this format doesn't have today. The natural design mirrors what GeoTIFF already does:

	GeoTIFF	equivalent here
"this block doesn't exist"	TileOffsets[i] == 0	index of present blocks
"this cell is invalid"	band's nodata value	per-array nodata in metadata.json

The cost is low — a typical workspace has ~1500 blocks, so the index is on the order of 1500 bits, versus the ~371 MB a per-cell mask would cost. read_block_* for an absent block would return the declared nodata without touching disk. The point of attention is disk_sync_model.py, which today writes every block unconditionally in write_block_core — the change would need to decide whether a fully-nodata block should be written back at all.

Finding: Zarr axis order is not guaranteed to be (y, x)

While investigating real (not synthetic) integration with disscube, CubeClient.load()/to_lucc_data() were found to do .transpose("y", "x") defensively before using any array — evidence that the axis order stored on disk isn't guaranteed. Reproduced with real xarray, writing exactly the way disscube's VariableWriter writes (da.to_dataset(...).to_zarr(...)): a square array with (x, y) axes instead of (y, x) has the same shape in both cases — the shape check in zarr_io.py didn't catch the inversion. Without a fix, this corrupted row/column silently, with no error at all.

Fix: Zarr v3 stores the real axis order in a native format field (arr.metadata.dimension_names, not attrs — unlike the older Zarr v2/xarray convention, _ARRAY_DIMENSIONS). load_zarr_into_workspace now reads that metadata and normalizes each block's read to (y, x)/(time, y, x) regardless of the physical on-disk order. Tested against three real cases (correct y,x, inverted x,y on a square array, inverted x,y,time with a time dimension) in tests/test_zarr_axis_order_regression.py, using real xarray to write — not plain zarr — to reproduce exactly what disscube produces.

Testing with large files

scripts/generate_and_benchmark.py generates synthetic data directly on disk, block by block (never materializing the whole grid in RAM to generate it — deterministic RNG per block position, hence reproducible) and runs Game of Life via MemmapRasterWorkspace, measuring wall time and real memory footprint.

python scripts/generate_and_benchmark.py \
    --shape 40000 40000 --block 512 512 --halo 1 \
    --generations 3 --density 0.35 --seed 42 \
    --root /tmp/haloexec_bench

Options: --shape HEIGHT WIDTH, --block BLOCK_H BLOCK_W, --halo, --generations, --density (fraction of initially alive cells), --seed, --root (workspace directory), --keep (don't delete the workspace at the end, to inspect the generated files).

Reported metric and why: the script separates RssAnon (real heap allocated by the process — the metric that proves whether the grid was materialized or not) from RssFile (page cache for touched mmap pages, reclaimable by the kernel, which grows with cumulative volume touched but doesn't represent "retained" memory). Using only ru_maxrss/VmRSS (RssAnon+RssFile combined) is misleading for mmap-based workflows — it grows even when the process never held the whole grid at once.

Reference result measured in an environment with ~3.9GB of RAM: a 40000×40000 grid (1.6 billion cells, ~1.5GB per array — more than 40% of the container's total RAM) ran with RssAnon around 130MB and a delta of ~4MB from the start, confirming the process footprint doesn't scale with grid size.

To test a different model instead of Game of Life, swap out the script's _game_of_life_rule function for any dict[str, np.ndarray] -> dict[str, np.ndarray] rule — the generation/benchmark mechanics don't change.

Migration path to dissmodel core

Once the JOSS review stabilizes, HaloChunkedRasterCellularAutomaton should migrate into dissmodel.geo.raster (or an equivalent module), keeping the same public API. Since the package already depends on dissmodel and reuses its real classes (not a reimplementation), the migration is literally moving the folder — no logic rewrite. tests/test_equivalence.py serves as the regression suite to confirm this.

No patch from this package should be applied to dissmodel core while it remains under review.