Implement fusion moves (#10)

* Implement edge-weighted graph watershed

* Add fusion moves implementation

* Add comparison to nifty

* Optimize fusion moves

* 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).
@@ -120,6 +171,25 @@ Tests run against the installed Python package. For each public function, cover:
 
 Correct first, fast second. Benchmark before adding complexity. Prefer algorithmic improvements over build-level tweaks. Release the GIL for expensive kernels. Add threading only after the single-threaded implementation is stable; it must be portable and user-controllable.
 
+### Profiling
+
+Measure before optimizing. The codebase carries a lightweight per-phase profiling utility for exactly this:
+
+- Header: `include/bioimage_cpp/detail/profile.hxx`. Macros: `BIOIMAGE_PROFILE_INIT(name)`, `BIOIMAGE_PROFILE_SCOPE(name, "label")`, `BIOIMAGE_PROFILE_REPORT(name)`.
+- Gated behind the `BIOIMAGE_PROFILE` compile-time flag. Outside of profile builds the macros expand to no-ops and a `NullProfiler` stub so the same code compiles unchanged.
+- Enable via CMake option: `pip install -e . --no-build-isolation -C cmake.define.BIOIMAGE_PROFILE=ON`. Rebuild without the flag for production work.
+- Reports per-phase wall-clock totals to stderr at the end of the instrumented scope (e.g., at the end of `optimize`). Same labels accumulate across multiple invocations of the scope (e.g., per-iteration phases).
+
+Workflow when adding or chasing a performance issue:
+
+1. **Compare standalone primitives first** (`development/.../check_*.py` scripts vs. nifty). If a primitive is already fast, the gap is elsewhere — don't optimize it speculatively.
+2. **Instrument the suspect function** by wrapping each logical phase in a `BIOIMAGE_PROFILE_SCOPE`. Pick labels that map to one operation each ("agreement_contract", "sub_solve", "energy_eval", ...), not full call paths.
+3. **Build with `BIOIMAGE_PROFILE=ON` and run a realistic problem** — typically the external multicut instance loaded by the comparison scripts. Run with `--repeats 1` so the report isn't drowned out.
+4. **Optimize the largest phase** (50% phase beats two 10% phases combined). Re-measure after each change; verify no other phase regressed.
+5. **Strip the instrumentation when done** only if it adds clutter; otherwise leave it in place — it's free when the flag is off.
+
+Don't add `std::chrono` snippets ad hoc; use the existing macros so future profiling sessions land in a consistent format.
+
 ## Documentation
 
 `README.md` covers: what `bioimage-cpp` is and isn't, install/build, minimal examples, design philosophy. Public functions need concise docstrings documenting input shapes, supported dtypes, output shapes/dtypes, copy behavior, background-label behavior (if relevant), and axis/coordinate conventions.

CMakeLists.txt

@@ -27,4 +27,9 @@ target_compile_options(_core PRIVATE
     $<$<NOT:$<CONFIG:Debug>>:-O3>
 )
 
+option(BIOIMAGE_PROFILE "Enable per-phase profiling instrumentation (development only)" OFF)
+if(BIOIMAGE_PROFILE)
+    target_compile_definitions(_core PRIVATE BIOIMAGE_PROFILE)
+endif()
+
 install(TARGETS _core LIBRARY DESTINATION bioimage_cpp)

MIGRATION_GUIDE.md

@@ -254,6 +254,50 @@ Notes:
 - `number_of_threads=0` uses the library default; pass a positive integer for a
   fixed thread count.
 
+## Edge-Weighted Watershed
+
+Nifty:
+
+```python
+import nifty.graph as ng
+
+labels = ng.edgeWeightedWatershedsSegmentation(graph, edge_weights, seeds)
+```
+
+bioimage-cpp:
+
+```python
+import bioimage_cpp as bic
+
+labels = bic.graph.edge_weighted_watershed(graph, edge_weights, seeds)
+```
+
+Notes:
+
+- Only the Kruskal variant of nifty's algorithm is provided. Edges are visited
+  in ascending weight order; two distinct components merge iff at least one is
+  unlabeled (seed label `0`), so seed boundaries are preserved.
+- `graph` may be an `UndirectedGraph` or a `RegionAdjacencyGraph`.
+- `edge_weights` must be 1D with length `graph.number_of_edges`. Supported
+  dtypes are `float32` and `float64`; other floating dtypes are promoted to
+  `float32` (matching nifty, whose Python binding is `float32`-only). Non-float
+  dtypes raise `TypeError`.
+- `seeds` must be 1D with length `graph.number_of_nodes`. Supported dtypes are
+  `uint32`, `uint64`, `int32`, `int64`. The value `0` marks unlabeled nodes;
+  non-zero ids are propagated along low-weight paths. Signed seed arrays must
+  not contain negative values.
+- The output is 1D with length `graph.number_of_nodes` and the same dtype as
+  `seeds`. Seed label values are preserved (no dense relabeling). Nodes that
+  no seed can reach remain `0`.
+
+Intentional differences vs. nifty:
+
+- No priority-queue variant — only the simpler sort + union-find Kruskal flow.
+  For the same input it matches nifty's default behavior (which also dispatches
+  to the Kruskal implementation).
+- No carving / background-bias variant. Build a carving prior into the edge
+  weights before calling the function if needed.
+
 ## Multicut
 
 Nifty exposes multicut through an objective + factory-style solver hierarchy.
@@ -364,6 +408,84 @@ Intentional differences vs. nifty:
   `GreedyAdditiveMulticut` and no fallthrough is given — the greedy solver
   already operates on each connected component internally.
 
+## Fusion Moves (Multicut)
+
+Nifty exposes the fusion-move multicut solver via the factory hierarchy with a
+chosen proposal generator and sub-solver factory.
+
+Nifty:
+
+```python
+import nifty.graph.opt.multicut as nmc
+
+objective = nmc.multicutObjective(graph, edge_costs)
+pgen = nmc.watershedProposals(sigma=1.0, numberOfSeeds=0.1)
+factory = nmc.fusionMoveBasedFactory(
+    proposalGenerator=pgen,
+    fusionMove=nmc.fusionMoveSettings(
+        mcFactory=nmc.greedyAdditiveFactory(),
+    ),
+    numberOfIterations=10,
+    stopIfNoImprovement=4,
+)
+labels = factory.create(objective).optimize()
+```
+
+bioimage-cpp:
+
+```python
+import bioimage_cpp as bic
+
+objective = bic.graph.MulticutObjective(graph, edge_costs)
+solver = bic.graph.FusionMoveMulticut(
+    proposal_generator=bic.graph.WatershedProposalGenerator(
+        sigma=1.0, n_seeds_fraction=0.1, seed=0,
+    ),
+    sub_solver=bic.graph.GreedyAdditiveMulticut(),
+    number_of_iterations=10,
+    stop_if_no_improvement=4,
+)
+labels = solver.optimize(objective)
+```
+
+Proposal generators:
+
+| nifty proposal generator | bioimage-cpp proposal generator |
+| --- | --- |
+| `watershedProposals(sigma=..., numberOfSeeds=...)` | `WatershedProposalGenerator(sigma=..., n_seeds_fraction=..., seed=...)` |
+| `greedyAdditiveProposals(sigma=..., weightStopCond=..., nodeNumStopCond=...)` | `GreedyAdditiveProposalGenerator(sigma=..., weight_stop=..., node_num_stop=..., seed=...)` |
+
+Sub-solvers: any built-in multicut solver (`GreedyAdditiveMulticut`,
+`GreedyFixationMulticut`, `KernighanLinMulticut`). If `sub_solver` is omitted,
+the default is `GreedyAdditiveMulticut` constructed with no-noise defaults.
+
+Intentional differences vs. nifty:
+
+- Single object construction: no separate factory / solver step.
+- Proposal generators are Python classes carrying their settings; the C++
+  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 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:
+
+- Custom Python proposal generators are not yet supported; subclass
+  `ProposalGenerator` and provide your own `_build` returning a C++
+  proposal-generator object if you need to extend the set.
+
 ## Segmentation Overlaps
 
 Nifty:

development/graph/multicut/check_fusion_move.py

@@ -0,0 +1,43 @@
+from __future__ import annotations
+
+import bioimage_cpp as bic
+
+from _compatibility import parser, run_comparison
+
+
+def main() -> None:
+    args = parser(
+        "Compare bioimage-cpp and nifty fusion-move multicut at matched settings."
+    ).parse_args()
+    # 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(),
+            fusionMove=objective.fusionMoveSettings(
+                mcFactory=objective.greedyAdditiveFactory(),
+            ),
+            numberOfIterations=10,
+            stopIfNoImprovement=4,
+            numberOfThreads=threads,
+            warmStartGreedy=True,
+        ),
+        args,
+    )
+
+
+if __name__ == "__main__":
+    main()

include/bioimage_cpp/detail/indexed_heap.hxx

@@ -133,6 +133,28 @@ public:
         sift_up(pos);
     }
 
+    // Bulk-load the heap from a list of entries in O(n) via Floyd's heapify.
+    // Replaces any current heap contents. Preconditions:
+    //   - the heap is empty before this call (call `clear()`/`reset_capacity`
+    //     first if needed);
+    //   - every key in `entries` is unique.
+    // Compared to N successive `push` calls (O(n log n)) this matters when
+    // initializing a heap from a large edge set in one shot.
+    void build_heap(std::vector<Entry> entries) {
+        heap_ = std::move(entries);
+        for (std::size_t pos = 0; pos < heap_.size(); ++pos) {
+            locator_.set(heap_[pos].key, pos);
+        }
+        if (heap_.size() < 2) {
+            return;
+        }
+        // Sift down from the last internal node back to the root. After each
+        // sift_down, the subtree rooted at `pos` is heap-ordered.
+        for (std::size_t pos = heap_.size() / 2; pos-- > 0;) {
+            sift_down(pos);
+        }
+    }
+
     // Precondition: `key` is currently in the heap.
     void change(const KeyT &key, PriorityT priority) {
         const auto pos = locator_.at(key);