Sweeping a cube along its chunk lines

This API runs neighbourhood kernels over chunked cubes. It batches the cells owned by each storage chunk with the halo needed by their rings, reducing the repeated chunk decoding caused by scalar-at-a-time access.

The traversal follows the cube's existing chunk grid. Each callback receives an ordinary in-memory cube containing the owned cells and the halo cells reached by their rings. A foreign chunk can serve several halos and may be read again.

The API separates a plan from its runner. The plan records the read order, chunk partitioning, halo, and estimated work before data access begins.

A    = dggread("dem.zarr")[:elevation]
plan = chunkplan(A; halo = 1)          # no data is read here

slope(c, v, vs) = maximum(abs(v - u) for u in vs; init = 0.0)

out = zeros(Float64, size(A))
mapneighbors!(out, slope, A, plan)

The plan

chunkplan reads boundaries from the data's chunk grid, including irregular layouts such as one chunk per ancestor subtree. It finds each halo by walking its boundary through halo, so planning reads metadata and performs CPU work without loading data chunks.

split divides a plan into contiguous pieces with similar chunk counts. For weighted work, unequal worker capacities, or grouping chunks by shared inputs, use the partitioning API. Build the plan before assigning work so its region conversion is shared by in-process tasks.

DiscreteGlobalGrids.chunkplan Function

chunkplan(A::AbstractDimArray; halo = 1, chunks = :auto, spatialdim = nothing,
          connectivity = Vertex()) -> MapChunkPlan
chunkplan(lk::AbstractCellLookup, bounds; halo = 1, connectivity = Vertex())

Plan a chunk-following sweep of A: which cells each chunk owns, and which cells outside it its rings reach.

chunks = :auto takes the boundaries from the data's own chunk grid (DiskArrays.eachchunk), which is the whole point — reading along them is what makes each chunk decode once. An Integer fixes a uniform chunk length in cells instead, and an in-memory array, having no chunk grid, is one chunk.

halo = n makes n rings of input context available to a chunk callback. It does not change the neighborhood passed by mapneighbors!: that operation selects neighbors with neighborhood = Disc(k) or Ring(k). The default plan halo is 1; a sweep requires at least its selected radius. spatialdim and connectivity follow mapneighbors.

The second form plans against a cell axis alone, with bounds either a chunk length or the chunk ranges themselves — the form to reach for when the cube is not the thing being read.

Building the plan walks each chunk's boundary once. That is CPU, not IO: no chunk of the data is touched until foreachchunk runs the plan.

plan = chunkplan(A; halo = 1)
nchunks(plan)                             # what it will read
foreachchunk(A, plan) do cc               # ... and reading it
    r = mapneighbors(f, chunkcube(cc))
    out[ownedindices(cc)] = parent(r)[localindices(cc)]
end

source

DiscreteGlobalGrids.MapChunkPlan Type

MapChunkPlan

What a chunk-following sweep will read, in what order — a value, not a running traversal.

Build one with chunkplan and run it with foreachchunk. It is worth being a value because three separate decisions are made on it:

• order: plan[i] is the ith chunk to visit. Reorder by rebuilding the plan from a permutation of its chunks.

• parallelism: split(plan, n) cuts it into n plans over disjoint chunks. Running those on separate tasks is what parallelises a sweep; there is no threaded keyword here because the split IS the decision.

• cost: nchunks(plan) and the per-chunk chunkhalo lengths say how much will be read before anything is.

The halo is computed when the plan is built, which costs one region walk per chunk and no IO. Chunk boundaries come from the data's own chunk grid, so a store whose chunks are irregular — one per ancestor subtree, say — is planned on its real boundaries and not on an assumed uniform length.

source

DiscreteGlobalGrids.MapChunk Type

MapChunk

One chunk of a MapChunkPlan: the axis indices it owns, and the axis indices outside it that its cells' rings reach.

ownedindices(mc) is the first, chunkhalo(mc) the second; both are indices in the cube's own cell axis, and both ascend. The halo holds only indices the axis really has — a ring that leaves the axis altogether is clipped here exactly as neighbors clips it.

source

DiscreteGlobalGrids.ownedindices Function

ownedindices(mc::MapChunk) -> UnitRange{Int}
ownedindices(cc::ChunkCube) -> UnitRange{Int}

The indices in the CELL AXIS OF THE CUBE THE CALLER PASSED that this chunk owns — the cells it produces results for. Contiguous, because a chunk is a contiguous block of that axis.

This is the axis index, not the complete level's numbering: on a cube over a subset of a level the two differ, and it is the caller's axis a result is written back to. axisindices names every cell of a loaded chunk in the same space, halo included.

source

DiscreteGlobalGrids.chunkhalo Function

chunkhalo(mc::MapChunk) -> Vector{Int}

The axis indices outside ownedindices that this chunk's cells reach, ascending. These are the cells a sweep over the chunk has to read but does not produce a result for.

source

chunkhalo(cc::ChunkCube) -> Vector{Int}

The indices in chunkcube that are HALO — read to complete the owned cells' rings, and not results to keep. The complement of localindices over the cube, so the two together are all of it.

These are indices in the block, not in the axis; chunkhalo of the plan's MapChunk gives the axis indices they came from.

source

DiscreteGlobalGrids.ChunkedLookups.nchunks Method

nchunks(plan::MapChunkPlan) -> Int

How many chunks the plan visits.

source

DiscreteGlobalGrids.halowidth Function

halowidth(plan::MapChunkPlan) -> Int

How many rings of context each chunk's halo carries. 1 is what a one-ring stencil needs and is the default.

source

Base.split Method

split(plan::MapChunkPlan, n::Integer) -> Vector{MapChunkPlan}

Cut the plan into at most n plans over disjoint chunks, in order.

This is how a chunked sweep is parallelised: run the pieces on separate tasks. Each chunk writes only the indices it owns, so pieces running at once cannot collide however the cut falls — and the result does not depend on how many pieces there were.

source

Running it

foreachchunk hands the callback a ChunkCube containing the chunk's cells and halo in memory, represented as an ordinary cube over a CellLookup. Package operations work on it as they do on any one-level cell cube; ownership metadata remains alongside the cube.

A chunk is a partial grid, so its indices are chunk-local. The owned cells form a contiguous range in the block, making localindices a range; ownedindices maps that range to the caller's cell axis. axisindices names every cell of the block, halo included, in that same axis, so axisindices(cc)[localindices(cc)] == ownedindices(cc). It is what lets a sweep over a chunk report numbers the caller can use, and it is how a field request is translated onto a chunk.

DiscreteGlobalGrids.foreachchunk Function

foreachchunk(f, A::AbstractDimArray, plan = chunkplan(A; kw...); kw...)

Run plan over A, calling f(cc::ChunkCube) once per chunk.

Each call gets the chunk's cells and its halo as an in-memory cube — one contiguous read for the chunk itself, and one read per foreign chunk the halo touches, never more. f's return value is discarded; this is the primitive the result-producing forms are built on, and the one to reach for when what a pass produces is not one value per cell.

Chunks are visited in the plan's order, on the calling task. To use more than one, split the plan and run the pieces — see MapChunkPlan.

plan = chunkplan(A; halo = 1)
@sync for p in Base.split(plan, Threads.nthreads())
    Threads.@spawn foreachchunk(A, p) do cc
        write!(out, ownedindices(cc), mine(mapneighbors(f, chunkcube(cc)), cc))
    end
end

source

DiscreteGlobalGrids.ChunkCube Type

ChunkCube

One chunk of a cube, read: its own cells and its halo, in memory, as an ordinary DimArray over a CellLookup.

This is what foreachchunk hands its callback, and the reason it is worth handing over rather than hiding: every verb in this package already works on a cube, so a chunk-following pass is written in the same vocabulary as the whole-cube pass it replaces.

A chunk is a partial grid, so an index into its cube is CHUNK-LOCAL and means nothing to the caller. Its accessors name cells in the two index spaces:

accessor	what it gives
chunkcube	the in-memory cube: the chunk's cells AND its halo
localindices	the indices IN THAT CUBE the chunk owns
ownedindices	those same owned cells' indices in the CALLER'S CELL AXIS
axisindices	EVERY cube cell's index in the caller's cell axis, halo included

The owned pair is one set of cells in the two numberings, so axisindices(cc)[localindices(cc)] == ownedindices(cc) always.

A result computed on the cube is meaningful for the owned indices only: a halo cell is there to complete its neighbours' rings, and its own ring is missing whatever fell outside the block. r[localindices(cc)] selects the part to keep and ownedindices(cc) says where it belongs.

Because the halo holds every axis neighbour of every owned cell, a stencil that reaches no further than the plan's halowidth computes exactly what the whole-axis sweep computes for those cells.

source

DiscreteGlobalGrids.chunkcube Function

chunkcube(cc::ChunkCube) -> AbstractDimArray

The chunk's cells and its halo, in memory, over a CellLookup.

source

DiscreteGlobalGrids.localindices Function

localindices(cc::ChunkCube) -> UnitRange{Int}

The indices in chunkcube that the chunk OWNS — the ones a result is kept for. Contiguous: a chunk is a run of the axis, and every halo cell is outside that run, so the owned cells stay together when the two are sorted into one axis.

source

DiscreteGlobalGrids.axisindices Function

axisindices(cc::ChunkCube) -> Vector{Int}

Where every cell of chunkcube came from in the CALLER'S CELL AXIS, in the cube's own cell order: entry k is the axis index of the cube's cell k, halo cells included. Ascending, since the halo below the chunk, the chunk itself and the halo above it are stored in that order.

This is the translation between the two numberings, and it is what lets a sweep over a chunk report indices the caller can use: Index(Local()) on a chunk is Value(axisindices(cc)), so a request answered chunk by chunk names the same cells the whole-axis sweep would have.

source

DiscreteGlobalGrids.globalindices Function

globalindices(mc::MapChunk) -> UnitRange{Int}
globalindices(cc::ChunkCube) -> UnitRange{Int}

Deprecated. Use ownedindices, which this forwards to, so existing calls keep their old behaviour exactly.

source

The sweeps built on it

Halo width describes available input data. The neighborhood selector chooses which cells the callback receives. Repeating a one-ring sweep is a different operation from one sweep over a wider neighborhood.

mapneighbors! keeps that invariant itself. Its halo defaults to the radius of the neighborhood selector, k for both Disc(k) and Ring(k) and one ring for Disc(0), and a supplied plan must carry at least that many rings: halowidth(plan) narrower than the radius throws an ArgumentError naming both widths. A wider halo than the radius is allowed. With A a stored cube, out a destination of the same shape and kernel a (cell, value, values) function:

plan = chunkplan(A; halo = 3)
mapneighbors!(out, kernel, A, plan; neighborhood = Ring(3))   # halowidth(plan) ≥ 3
mapneighbors!(out, kernel, A; neighborhood = Disc(3))         # plans halo = 3 itself

mapneighbors! is the streaming form: results are written into dest a chunk at a time, so neither the input nor the output has to fit in memory. mapneighbors with pass = Values() takes the same route by itself whenever the cube's data is chunked, and collects the results.

A field request takes it too. mapneighbors!(dest, f, A, plan; needs = (Value(dem), Centroid())) states the request once, about the cube that was passed, and the route translates it onto each chunk: Index(Local()) keeps answering the caller's cell-axis index, never a chunk-local one, and each stored Value is read along its own chunk grid the way the swept data is. mapneighbors and foreachneighbors reach for it by themselves under the rule Values() uses — a chunked parent and order = StorageOrder().

plan = chunkplan(A; halo = 1)
out  = zeros(Float64, size(A))
mapneighbors!(out, steepest, A, plan; needs = (Value(A), Centroid()))

Neighbors uses a different callback contract: it supplies cell handles and the callback reads values from the original array. Use Values() when a chunked sweep should stream the fields through the traversal.

The neighbourhood API documents these kernels and callback forms.

Index