parallel processing of disconnected graphs

[DirkEilander]

Kind of request

None

Enhancement Description

Currently all links are saved in an single flat array. Investigate ways to distinguish between disconnected graphs (i.e. different basins within the same data) or stream order which can be processed in parallel.

Use case

No response

Additional Context

No response

[hellolulujiang]

The two halves of this issue need different structures, and both fall out of the
cell ordering rather than out of any change to the flat array. I have been
working on exactly this; implementation and numbers are on
https://github.com/hellolulujiang/pyflwdir/tree/orderings, see also #114. I read
"stream order which can be processed in parallel" as an ordering of the cells
that splits across workers, not as Strahler order.

A note on the numbers before anything else. They are one 4000 x 4000 window of
MERIT Hydro at 30 arcsec, 16 million cells with 9,755 pits, timed on my personal
MacBook, which keeps running other work. The cores share memory bandwidth with
whatever else is active, so the timings, the threaded ones especially, are
indicative: lean on the comparisons, not the absolute times. Kernel time
excludes building the ordering, which is built once and cached.

Disconnected graphs: a depth-first ordering separates them by itself

A depth-first traversal from the pits finishes one pit's whole tree before
starting the next, so each basin is one contiguous run of the sequence and the
runs start at the pits:

s = seq_dfs                                    # preorder, pits first
start = np.append(np.flatnonzero(idxs_ds[s] == s), s.size)
# basin b occupies s[start[b]:start[b+1]]

No labelling pass and no extra array; basins() holds the same information but
itself needs a sequence first. On the window above, the depth-first sequence has
exactly one contiguous run per basin, all 9,755 of them, where the breadth-first
sequence interleaves the basins level by level.

The runs share nothing, so they go to threads with no atomics and no locks.
Cells in loops stay outside the sequence as they do today. A raster that is a single
basin gets nothing from this.

Within one basin: layers work, with one extra rule

When the raster is one basin, or a few large basins leave the other cores idle,
the split has to go inside a basin. The rule that is easy to miss: grouping
cells so that none drains into another in its own group is not enough, because
two cells can drain into the same cell, and then two threads write it at
once. So a cell joins the earliest layer that its upstream cells allow and
in which no other cell already writes to its downstream cell. A layer then runs
with no atomics, the result is identical for any thread count, and the threaded
path can be tested against the serial one exactly. The constraint is nearly
free, 3 extra layers out of 4,180 here.

What each is worth

Splitting by basin sped accuflux up about three and a half times on eight
cores, part of which is the depth-first ordering itself rather than the
threading, and the result matches the serial sweep to rounding. It cannot reach
the core count, for two reasons that are worth naming: the kernel is limited by
memory traffic, which the cores share, and the basins are unequal, so the
largest one bounds how evenly the work can spread. Which of the two dominates I
have not measured.

Splitting by layer never got past about one and a half times and fell back
beyond four threads. That is structural rather than bad luck: every layer is a
barrier, and near the headwaters a layer holds only a handful of cells, so the
scheme pays one synchronisation per layer while giving the threads almost
nothing to do between barriers. On a single-basin raster the size of the Rhine
example it loses at any thread count.

So: basin-level first, it is one flatnonzero over a sequence worth having
anyway. The layering matters as the correctness rule if these kernels are ever
threaded, more than as a speedup. examples/bench_parallel.py on the branch
holds both implementations and reruns both experiments on any d8 raster.

Happy to open a pull request for the basin split once #115 and the depth-first
ordering are settled; it is a few lines on top of them.