Implement parallelization for fusion movs

Author: constantinpape

AGENTS.md

@@ -33,6 +33,57 @@ Before introducing union-finds, priority queues, edge hashing, stride math, thre
 
 If a needed helper does not exist but is generally useful, add it to `detail/` as a header-only utility with a focused API rather than inlining it into one algorithm. Keep the contract small and well-documented; over time `detail/` is what lets multiple modules stay small and consistent.
 
+## Reusable algorithmic infrastructure
+
+Some larger pieces of infrastructure sit above `detail/` but are still
+intended to be reused across objective types. Check these before starting a
+new fusion-move / proposal-based / contraction-based solver:
+
+- `include/bioimage_cpp/graph/detail/fusion_contract.hxx` — objective-agnostic
+  agreement-projection primitive. `contract_by_agreement(graph, proposals,
+  n_proposals, ...)` returns `{contracted_graph, contracted_edge_of_original,
+  root_of_node}`. Already supports N ≥ 1 proposals; reused by both pairwise
+  and joint multi-proposal fuses.
+- `include/bioimage_cpp/graph/proposal_generator.hxx` — `ProposalGeneratorBase`
+  abstract class. Concrete generators (Watershed, GreedyAdditiveMulticut) live
+  in `proposal_generators/` and depend only on `(graph, edge_costs)` plus an
+  RNG seed. They emit `std::vector<std::uint64_t>` node labelings and are
+  therefore reusable across multicut, lifted multicut, mincut, etc.
+- `include/bioimage_cpp/graph/multicut/greedy_additive.hxx` — exposes
+  `GreedyAdditiveWorkspace` so multiple invocations on different graphs
+  share scratch buffers. Use this pattern (workspace + `reset(graph)`) when
+  a fusion-move driver calls a sub-solver inside its iteration loop.
+- `UndirectedGraph::from_sorted_unique_edges(N, edges, populate_lookup=false)`
+  — bulk graph construction without the per-edge hash insertion in
+  `insert_edge`. Pair with `populate_lookup=false` when the consumer only
+  walks edges / adjacency (the multicut sub-solvers do).
+- `detail/threading.hxx::parallel_for_chunks` — the only threading primitive
+  we use. New parallel solvers should not introduce alternatives.
+
+When porting fusion moves to a new objective (e.g. lifted multicut):
+
+1. The driver loop in `multicut/fusion_move.hxx::FusionMoveSolver::optimize`
+   is short and dense. Duplicate it for the new objective rather than
+   abstracting it via a template/CRTP base — the moving parts (cost
+   aggregation, energy evaluator, sub-solver type) are objective-specific
+   and template gymnastics buy little.
+2. Reuse `contract_by_agreement` unchanged; it operates on the *base* graph
+   only.
+3. Write a new `fuse_multi(...)` that aggregates *both* base and lifted (or
+   other auxiliary) weights through `contraction.contracted_edge_of_original`
+   and `contraction.root_of_node`, calls the new objective's sub-solver, and
+   lifts labels back via `root_of_node`.
+4. Reuse the existing `WatershedProposalGenerator` and
+   `GreedyAdditiveMulticutProposalGenerator` verbatim; they only depend on
+   the base graph + base costs and emit node labelings. Add objective-specific
+   generators (e.g. `GreedyAdditiveLiftedMulticutProposalGenerator`) only if a
+   meaningful new proposal strategy emerges.
+5. Reuse the per-thread parallel pattern from `optimize`: stage-1 parallel
+   proposal generation + parallel pairwise fuse, stage-2 sequential joint
+   multi-fuse on leftover candidates. Per-thread `GreedyAdditiveWorkspace`
+   becomes per-thread `<NewObjective>Workspace` if the new sub-solver follows
+   the same pattern.
+
 ## Dependencies
 
 **Allowed**: C++20 stdlib, `nanobind`, `numpy`, `scikit-build-core`, `cmake`, `pytest`, `cibuildwheel` (CI only).

MIGRATION_GUIDE.md

@@ -466,12 +466,19 @@ Intentional differences vs. nifty:
   proposal-generator object is built lazily when `optimize` is called.
 - The driver warm-starts from the trivial singleton labeling by running the
   default greedy-additive sub-solver once before the proposal loop.
-- A best-of-three safety net (current, proposal, fused) keeps the running
-  energy monotonically non-increasing across iterations.
-- The current implementation is single-threaded and uses pairwise
-  (proposal, current) fuses. `number_of_threads` and
-  `number_of_parallel_proposals` are kept on the API for forward compatibility
-  but must be left at their defaults (`1` and `2` respectively).
+- A best-of safety net keeps the running energy monotonically non-increasing
+  across iterations (compared against current, proposals, fused, and the
+  stage-2 joint fuse).
+- Parallel proposal generation and a multi-proposal joint fuse are supported:
+  `number_of_threads=T` runs `number_of_parallel_proposals=P` proposal
+  generators in parallel within each iteration. By default `P=2` when `T=1`
+  and `P=T` when `T>1`; pass an explicit `number_of_parallel_proposals` to
+  override. Each parallel slot uses an independent proposal generator with
+  seed `proposal_generator.seed + slot_index` so the result is deterministic
+  for a given `(seed, T, P)`. When at least two parallel pairwise fuses fail
+  to improve on the current best, a joint multi-proposal fuse runs over the
+  surviving fused candidates (matches nifty's `ccFusionMoveBased` stage-2
+  behaviour).
 
 Notes:
 

development/graph/multicut/check_fusion_move.py

@@ -7,18 +7,23 @@
 
 def main() -> None:
     args = parser(
-        "Compare bioimage-cpp and nifty fusion-move multicut at default settings."
+        "Compare bioimage-cpp and nifty fusion-move multicut at matched settings."
     ).parse_args()
-    # bioimage-cpp warm-starts from the trivial singleton labeling with a
-    # greedy-additive pass before the proposal loop. Nifty's ccFusionMoveBased
-    # exposes the same behaviour via the `warmStartGreedy=True` flag (which
-    # internally chains a greedyAdditiveFactory in front of the fusion-move
-    # factory). Both sides therefore enter the proposal loop from the same
-    # starting point.
+    # Match settings on both sides:
+    #   - threads / parallel-proposals = `--threads N` (default 1). We set
+    #     P = T explicitly on both sides so the comparison is apples-to-apples
+    #     regardless of either library's API default for P.
+    #   - greedy-additive warm-start (we do it automatically when initial
+    #     labels are the trivial singleton; nifty exposes it as
+    #     `warmStartGreedy=True`).
+    #   - greedy-additive sub-solver.
+    threads = int(args.threads)
     run_comparison(
         "fusion_move",
         lambda: bic.graph.FusionMoveMulticut(
             proposal_generator=bic.graph.WatershedProposalGenerator(),
+            number_of_threads=threads,
+            number_of_parallel_proposals=threads,
         ),
         lambda objective: objective.ccFusionMoveBasedFactory(
             proposalGenerator=objective.watershedCcProposals(),
@@ -27,7 +32,7 @@ def main() -> None:
             ),
             numberOfIterations=10,
             stopIfNoImprovement=4,
-            numberOfThreads=1,
+            numberOfThreads=threads,
             warmStartGreedy=True,
         ),
         args,

include/bioimage_cpp/detail/profile.hxx

@@ -5,6 +5,12 @@
 // and adds no overhead. Use it in development/profile-mode builds to find
 // hotspots; do not enable it in production wheels.
 
+// `NullProfiler` is always available; helper templates can request "no
+// profiling here" via the same type regardless of build mode.
+namespace bioimage_cpp::detail {
+struct NullProfiler {};
+} // namespace bioimage_cpp::detail
+
 #ifdef BIOIMAGE_PROFILE
 #include <chrono>
 #include <cstdio>
@@ -63,18 +69,31 @@ private:
     Profiler::Clock::time_point start_;
 };
 
+// Overload that accepts a NullProfiler and does nothing. Lets the same
+// macros expand cleanly when a profiled translation unit calls into a helper
+// that explicitly does not want to participate in profiling (e.g. work inside
+// parallel workers, where we measure wall-clock at the dispatch level).
+class ProfileTimerNull {
+public:
+    ProfileTimerNull(NullProfiler &, const char *) {}
+};
+
+inline ProfileTimer make_profile_timer(Profiler &profiler, const char *name) {
+    return ProfileTimer(profiler, name);
+}
+
+inline ProfileTimerNull make_profile_timer(NullProfiler &profiler, const char *name) {
+    return ProfileTimerNull(profiler, name);
+}
+
 } // namespace bioimage_cpp::detail
 
 #define BIOIMAGE_PROFILE_INIT(var) ::bioimage_cpp::detail::Profiler var;
-#define BIOIMAGE_PROFILE_SCOPE(var, name) ::bioimage_cpp::detail::ProfileTimer _bp_##__LINE__(var, name);
+#define BIOIMAGE_PROFILE_SCOPE(var, name) auto _bp_##__LINE__ = ::bioimage_cpp::detail::make_profile_timer(var, name);
 #define BIOIMAGE_PROFILE_REPORT(var) (var).report();
 
 #else
 
-namespace bioimage_cpp::detail {
-struct NullProfiler {};
-} // namespace bioimage_cpp::detail
-
 #define BIOIMAGE_PROFILE_INIT(var) ::bioimage_cpp::detail::NullProfiler var;
 #define BIOIMAGE_PROFILE_SCOPE(var, name) (void)var;
 #define BIOIMAGE_PROFILE_REPORT(var) (void)var;