Batched BFS distance maps

cuplan.distance_maps(grid, sources, backend="auto") computes the exact 4-connected distance from each of n sources to every cell — n flood fills as one batch.

Why it is the workhorse

Every serious MAPF implementation runs on exact goal-distance tables: they are the only admissible heuristic worth having on a map with walls, and PIBT queries them millions of times. pymapf builds them with one backward Dijkstra per agent, serially. On a unit-cost grid Dijkstra is breadth-first search, and BFS is a wavefront — every frontier cell can be expanded simultaneously.

Parallelization

The batch is a (n, height, width) volume advanced one wave per sweep:

• CPU reference — four array shifts OR-ed together per wave, vectorized over the whole batch with NumPy.
• CUDA — a gather ("pull") kernel: at wave t every unlabelled free cell checks whether any neighbour was labelled t−1 (Merrill, Garland and Grimshaw 2012, Scalable GPU graph traversal, PPoPP). Pull needs no atomics and stays coalesced; the cost of touching settled cells each wave is amortized across the batch.

The number of kernel launches equals the graph diameter, not the number of sources — batching 512 sources costs barely more than 16. Measured scaling is in Benchmarks.

Semantics

• Returns int32 distances, -1 for unreachable cells (blocked cells included).
• Matches pymapf.algorithms.search.distance_table values exactly (tested against a reference Python BFS, and CPU == CUDA is asserted on GPU machines).

Batched grid BFS: one flood-fill distance map per source.

This is the workhorse primitive of the library. MAPF solvers consume exact goal-distance tables — pymapf computes one backward Dijkstra per agent, serially. On a 4-connected grid with unit edge costs Dijkstra is BFS, and BFS is a wavefront: every cell on the frontier can be expanded simultaneously. Batching N sources into one (N, H, W) volume turns the whole heuristic-table build into makespan fully parallel sweeps.

Reference: Merrill, D.; Garland, M.; and Grimshaw, A. 2012. Scalable GPU graph traversal. PPoPP 2012: 117-128 (frontier-parallel BFS).

distance_maps

distance_maps(grid, sources, backend='auto', timings=None)

Return exact 4-connected distances from each source to every cell.

Parameters:

Name	Type	Description	Default
grid	Grid	the occupancy grid.	required
sources	ndarray	(n, 2) array of (row, col) source cells. Each must be a free cell.	required
backend	Backend	"auto", "cpu" or "cuda".	'auto'
timings	dict[str, float] | None	optional dict the CUDA backend fills with per-phase wall times in seconds — h2d (upload + device allocation), kernel (the wave loop, including the per-wave termination check, which is a 4-byte device read), and d2h (copying the finished maps back). The CPU backend leaves it untouched. Phase boundaries are device synchronization points, so profiling adds a small cost; benchmark totals should come from an unprofiled run.	None

Returns:

Type	Description
ndarray	(n, height, width) int32 array; entry [i, r, c] is the
ndarray	length of a shortest path from sources[i] to (r, c), or
ndarray	-1 where unreachable (including blocked cells).

Source code in cuplan/bfs.py

def distance_maps(
    grid: Grid,
    sources: np.ndarray,
    backend: Backend = "auto",
    timings: dict[str, float] | None = None,
) -> np.ndarray:
    """Return exact 4-connected distances from each source to every cell.

    Args:
        grid: the occupancy grid.
        sources: ``(n, 2)`` array of ``(row, col)`` source cells. Each
            must be a free cell.
        backend: ``"auto"``, ``"cpu"`` or ``"cuda"``.
        timings: optional dict the CUDA backend fills with per-phase
            wall times in seconds — ``h2d`` (upload + device
            allocation), ``kernel`` (the wave loop, including the
            per-wave termination check, which is a 4-byte device read),
            and ``d2h`` (copying the finished maps back). The CPU
            backend leaves it untouched. Phase boundaries are device
            synchronization points, so profiling adds a small cost;
            benchmark totals should come from an unprofiled run.

    Returns:
        ``(n, height, width)`` int32 array; entry ``[i, r, c]`` is the
        length of a shortest path from ``sources[i]`` to ``(r, c)``, or
        ``-1`` where unreachable (including blocked cells).
    """
    sources = np.atleast_2d(np.asarray(sources, dtype=np.int32))
    if sources.ndim != 2 or sources.shape[1] != 2:
        raise ValueError("sources must have shape (n, 2)")
    for r, c in sources:
        if not grid.is_free((int(r), int(c))):
            raise ValueError(f"source ({r}, {c}) is blocked or out of bounds")
    if resolve_backend(backend) == "cuda":
        return _distance_maps_cuda(grid, sources, timings)
    return _distance_maps_cpu(grid, sources)