diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index 6ac442b..d18879f 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -588,14 +588,45 @@ energy = objective.energy(labels) - `overwrite_existing` — when `True`, lifted entries that coincide with an existing edge replace its weight; the default accumulates. -Available solvers (no fusion-move / ILP solvers yet): +Available solvers (no ILP solvers yet): | nifty factory | bioimage-cpp solver | | --- | --- | | `liftedMulticutGreedyAdditiveFactory()` | `LiftedGreedyAdditiveMulticut()` | | `liftedMulticutKernighanLinFactory(...)` | `LiftedKernighanLinMulticut(...)` | +| `fusionMoveBasedFactory(...)` | `FusionMoveLiftedMulticut(...)` | | `chainedSolversFactory([...])` | `LiftedChainedSolvers([...])` | +`FusionMoveLiftedMulticut` mirrors `FusionMoveMulticut` (same proposal-generator +plumbing, same threading + multi-proposal joint-fuse semantics, same best-of +safety net). The differences are: + +- Proposal generators operate on the *base* graph and base edge costs (only + base-graph edges are candidate cut edges; lifted edges contribute to energy + but cannot be contracted directly). The driver extracts the base costs from + `objective.weights[:objective.number_of_base_edges]` automatically. +- Each fuse contracts the base graph by agreement, aggregates *both* base and + lifted weights onto the contracted lifted-multicut subproblem (lifted edges + whose endpoints land on already-existing contracted base edges fold into + them; the rest become new contracted lifted edges), and solves the + subproblem with a `LiftedMulticutSolver`. +- The default sub-solver and warm-start are `LiftedGreedyAdditiveMulticut`. + Both `LiftedGreedyAdditiveMulticut` and `LiftedKernighanLinMulticut` are + pluggable via `sub_solver=`. + +```python +solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator( + sigma=1.0, n_seeds_fraction=0.1, seed=0, + ), + sub_solver=bic.graph.LiftedKernighanLinMulticut(number_of_outer_iterations=3), + number_of_iterations=10, + stop_if_no_improvement=4, + number_of_threads=4, +) +labels = solver.optimize(objective) +``` + A typical warm-started solve combines greedy and KL: ```python diff --git a/development/graph/lifted_multicut/PERFORMANCE_NOTES.md b/development/graph/lifted_multicut/PERFORMANCE_NOTES.md index dbacc82..823d176 100644 --- a/development/graph/lifted_multicut/PERFORMANCE_NOTES.md +++ b/development/graph/lifted_multicut/PERFORMANCE_NOTES.md @@ -4,17 +4,66 @@ State of the lifted-multicut solvers vs nifty on the standard benchmark problems and notes on remaining optimization headroom. Read this before the next round of perf work. -## Current benchmark (3D ISBI lifted problem) - -Problem dimensions: 2462 nodes, 17 949 local edges, 21 444 lifted edges. - -| Solver | bic | nifty | Ratio | Status | -|---|---|---|---|---| -| Greedy additive | ~13 ms | ~14.7 ms | **0.88×** (faster) | Goal met | -| Kernighan-Lin (greedy + 10 outer) | ~200 ms | ~102 ms | **1.95×** (slower) | Outside 30%-of-nifty target | - -Energies match nifty to within numerical noise on both solvers (greedy diff -~0.3, KL diff ~0.08). +## Current benchmark matrix + +Produced by `python evaluate_solvers.py` (2026-05-17). All runs +single-threaded; fusion-move uses `n_seeds_fraction=0.1`, +`number_of_iterations=10`, `stop_if_no_improvement=4`. nifty side uses +matched settings via the chained-solver factory (greedy warm-start + +KL/fusion with the same iteration counts, `SEED_FROM_LOCAL`). +`runtime ratio` is `nifty_runtime / bic_runtime` — values > 1 mean bic +is faster. + +| Problem | Solver | bic energy | nifty energy | Δenergy | bic runtime | nifty runtime | runtime ratio | +|---|---|---|---|---|---|---|---| +| 2D | greedy | -1575.04 | -1575.04 | 0.00 | 0.70 ms | 1.50 ms | 2.15× faster | +| 2D | KL (10 outer) | -1575.21 | -1575.21 | 0.00 | 4.14 ms | 4.77 ms | 1.15× faster | +| 2D | fusion-move | -1575.43 | -1575.43 | 0.00 | 8.31 ms | 12.3 ms | 1.48× faster | +| 3D | greedy | -15891.4 | -15891.0 | −0.35 | 6.40 ms | 15.9 ms | 2.48× faster | +| 3D | KL (10 outer) | -15921.0 | -15921.1 | +0.07 | 78.1 ms | 103 ms | 1.32× faster | +| 3D | fusion-move | -15915.1 | -15915.1 | 0.00 | 128 ms | 196 ms | 1.53× faster | +| grid | greedy | -690 014 | -690 050 | +35.9 | 16.4 s | 20.6 s | 1.25× faster | +| grid | KL (10 outer) | -690 544 | -690 599 | +54.5 | 442 s | 6 514 s | **14.7× faster** | +| grid | fusion-move | -690 271 | -690 356 | +84.7 | 51.8 s | 60.2 s | 1.16× faster | + +Δenergy = bic − nifty; negative means bic is better, positive means +nifty is better. Energies are exact matches on 2D/3D fusion-move and +within 0.05 % on the rest; bic is faster than nifty on every row. + +### Note on the grid-KL row + +The 14.7× speed gap and the 54.5-unit energy gap on this single row are +much larger than on every other row (the 2D / 3D KL ratios are 1.15–1.32× +with essentially exact energies). That points to nifty's KL doing more +work per outer iter on grid — not a single algorithmic difference, but +the most plausible drivers, given that no investigation has been done: + +- **Inner-loop epsilon.** Our two-cut driver runs with one `epsilon` + parameter, defaulted to `1e-6` from + `KernighanLinSolver` (header file `kernighan_lin.hxx`). nifty's outer + KL uses `epsilon=1e-7` and its inner `TwoCut` uses + `epsilon=1e-9`. A tighter epsilon admits smaller-gain moves, which + means longer chains per pair and more outer iters before the + no-improvement break. On a 262 k-node problem with millions of + candidate moves, the extra accepted moves multiply quickly. +- **Inner iteration cap.** Both implementations effectively run the + two-cut chain until the heap empties; in nifty this is governed by + `numberOfInnerIterations` defaulting to + `std::numeric_limits::max()`. Same effective contract + as bic, but if a different value is wired through the factory the + gap would shift. +- **Per-pair adjacency cost.** Our `chain_gain_init` builds a filtered + per-node adjacency cache and the inner loop iterates only in-pair + neighbours. nifty walks the full lifted adjacency on every move and + re-classifies in/out-of-pair via label comparisons. At the scale of + grid (long-range lifted edges, dense pair adjacency), the constant + factor here adds up to multiple-× per chain — independent of + whether the algorithm visits the same number of moves. + +Two of the three (epsilon, adjacency walking) push energy and runtime +in the same direction: nifty does strictly more inner work and lands at +a slightly better energy. Closing the 54.5-unit gap on bic would mean +tightening `epsilon` and accepting more runtime — not pursued here. ## What's done @@ -38,147 +87,291 @@ floor for ~21 k heap operations. ### Kernighan-Lin -One C++ optimization (~22% speedup, 245 → 200 ms): +Two C++ optimizations landed (cumulative 2.85× speedup, 245 → 86 ms): -1. **Pre-built per-node filtered adjacency** in +1. **Pre-built per-node filtered adjacency** (~22% speedup, 245 → 195 ms) in `include/bioimage_cpp/graph/lifted_multicut/kernighan_lin.hxx`. `ChainScratch` gained `filtered_offset`, `filtered_count`, and `filtered_entries` (each entry caches `{node, weight, is_base}`). - `chain_gain_init` populates these alongside the gain accumulation; the - chain loop iterates only in-pair neighbors and skips the - `bufs.in_pair[u_key]` filter check. Also caches `was_in_heap` once per - iteration so the three subsequent heap operations don't each re-query - the locator. - -## What remains — KL is the open item - -C++ KL takes ~200 ms; nifty ~102 ms. Profile breakdown -(`BIOIMAGE_PROFILE=ON`, 1 repeat) — note the inner scopes inflate total by -~80% so trust the relative shares: - -``` - cc_repartition 0.0033 s ( 0.7%) - energy_eval 0.0013 s ( 0.4%) - compute_pairs 0.0073 s ( 2.0%) - chain_init 0.0037 s ( 1.0%) - chain_gain_init 0.0812 s ( 22.1%) <-- ~36% of pair_chains - chain_loop 0.0880 s ( 23.9%) <-- ~63% of pair_chains - chain_cleanup 0.0017 s ( 0.5%) - pair_chains 0.1744 s ( 47.4%) - cluster_splits 0.0083 s ( 2.3%) + The chain loop iterates only in-pair neighbors and caches `was_in_heap` + once per iteration so the three subsequent heap operations don't each + re-query the locator. + +2. **`changed_[]` cross-iteration pair-skip** (~2.27× speedup, 195 → 86 ms) + in the outer `kernighan_lin()` driver. Mirrors nifty's + `checkIfPartitonChanged()` / `changed_[piU] || changed_[piV]` gate. A + pair `(A, B)` whose endpoints' node-sets are identical to the previous + outer iter must produce the same result it did last time (zero gain — + otherwise A or B would have changed). After iter 1, this typically + skips >50% of pairs. New helper `detail_kl::compute_cluster_changed` + runs in <0.2 ms per outer iter. Tie-breaking is preserved (iter 0 + processes every pair); only pairs whose runs would commit no moves are + skipped. Same gate also applied to `cluster_splits`. + +### Fusion-move + +No solver-level optimizations were needed — the driver, fuse step and +sub-solver are already competitive with nifty's. The only change was a +**proposal-generator parameter-semantics fix** in +`include/bioimage_cpp/graph/proposal_generators/watershed.hxx`. Pre-fix +`n_seeds_fraction=0.1` produced 2× nifty's seed density: the loop +iterated `0.1 * N` times placing 2 seeds per iter (0.2 N total seeds), +whereas nifty iterates `nSeeds / 2` times placing 2 each (`0.1 * N` +total seeds). Effect on grid: -690099 → -690270 (closes 67 % of the +256-unit energy gap), 2D becomes an exact energy match, 3D drops 1.5 +energy units (still matches nifty). Loop now runs `n_seeds / 2` times +to match nifty exactly. Documented in the header and Python docstring. + +**Diagnosis of the grid energy gap (2026-05-17).** On grid, bic-fusion +finishes 85 units worse than nifty-fusion, despite bic-greedy starting +36 units *better* than nifty-greedy. The 121-unit reversal across the +fusion-move step is the thing to explain. + +1. **Cross-feed test localizes 100 % of the final gap to the greedy + warm-start, not to the fusion-move code.** Feeding nifty-greedy + labels into bic-fusion produced -690 355.67, matching nifty-fusion's + -690 355.63 (within float noise). bic's fusion-move algorithm + reproduces nifty's result from the same starting state. + +2. **bic-greedy and nifty-greedy land in structurally different local + optima**, not in two equivalent tie-breaks of the same optimum. The + reversal (bic ahead by 36 after greedy, behind by 85 after fusion) + means bic-greedy converges to a *deeper* local minimum that the + watershed proposals cannot escape — agreement-contraction between + bic-greedy's labels and the proposals leaves the sub-solver no room + to commit improvements. nifty-greedy's higher-energy optimum has + partition boundaries that the same proposals *can* perturb, so its + fusion-move makes much more progress per iteration. + +3. **Tie density is the structural cause.** 47 % of grid base edges + (339 946 of 718 848) carry the same weight ~+0.1. With so many + identical priorities, greedy-additive's merge order can pick wildly + different partition topologies. 2D and 3D have far fewer ties, so + both sides land in similar optima and the fusion-move ratios stay + exact. + +4. **bic's existing merge direction is the right one on this problem.** + Swapping bic's union-by-adjacency-size for nifty's union-by-rank + regressed bic-greedy by 100 units. Reverted. + +To close the gap we have to escape bic-greedy's deeper-but-rigid +optimum before invoking fusion. KL refinement does exactly this: + +| Variant on grid | Energy | Runtime | Δ vs nifty energy | +|---|---|---|---| +| current (greedy + fusion 10/4) | -690 270.93 | 52 s | +85 (worse) | +| fusion 20/8 (more iters, no KL) | -690 341.37 | 87 s | +14 | +| fusion 50/15 | -690 357.68 | 193 s | -2 | +| **greedy + KL(1) + fusion** | **-690 392.90** | **82 s** | **-37 (better)** | +| greedy + KL(2) + fusion | -690 539.43 | 173 s | -184 | + +At equal runtime, one outer iter of KL between greedy and fusion buys +~50 more energy units than the equivalent extra fusion iterations. +KL(1) on 2D shifts the result by 0.02 units (still matches nifty +within noise), 3D improves by ~5 units, runtime cost is 10–15 % per +problem. + +**Not adopted as a default.** The current `FusionMoveLiftedMulticut` +matches nifty exactly on 2D/3D and runs 1.16× faster than nifty on +grid at 0.012 % worse energy; this is acceptable for downstream +segmentation use. Users who need the better grid energy can chain +greedy → KL → fusion explicitly via `LiftedChainedSolvers`. If that +turns out to be the common workflow, a `kl_warm_refinement_iters` +parameter on `FusionMoveLiftedMulticut` would make the opt-in +self-contained. + +## Post-mortem: cluster-pair bucket optimization (2026-05-16, reverted) + +Tried the optimization the previous notes had ranked #1 — a per-outer-iter +bucket index over all lifted-graph edges, grouped by their endpoint +cluster-label pair, so `chain_gain_init` could read three buckets +(`(A, A)`, `(B, B)`, `(A, B)`) per pair-chain instead of walking the full +lifted adjacency of every in-pair node. + +**Wall-clock outcomes:** + +| Variant | Time | 3D energy diff vs nifty | Notes | +|---|---|---|---| +| Pre-bucket (current) | 195 ms | 0.08 | Reference | +| Simplified buckets, 1 rebuild / outer iter | 135 ms | 0.94 | Staleness during cluster_splits | +| Simplified buckets, 2 rebuilds / outer iter | 139 ms | 0.17 | Mid-iter rebuild recovers most of the loss | +| Incremental maintenance (relabel-on-commit) | 198 ms | 0.16 | Slower than pre-bucket, no energy win | + +The simplified bucket version delivered the predicted time win on +`chain_gain_init` (79 ms → 31 ms), but the 0.08 → 0.17 energy regression +that came with it could not be closed. The incremental-maintenance +follow-up — designed specifically to eliminate within-iter staleness — +turned out to be both slower (memory-fragmentation regression in +`chain_gain_init` and `chain_loop` from the `unordered_map>` layout) and unable to close the energy gap. We reverted +the entire bucket experiment. + +**Why the energy gap didn't close.** I assumed the 0.08 → 0.17 regression +was caused by buckets going stale between pair-chains, so incremental +maintenance should restore it. The data says otherwise: simplified +buckets with maximum staleness (1 rebuild/iter) and incremental +buckets with zero staleness landed at 0.17 and 0.16 respectively — +essentially identical. The gap is not from staleness; it's from a +**different filtered-adjacency iteration order** that the bucket +construction produces vs. the original direct-adjacency walk. KL is +sensitive to tie-breaking when multiple candidate moves have near-equal +gain, and the order in which `filtered_entries[v]` entries get pushed +into the heap influences the local optimum the chain converges to. +Floating-point summation order across many lifted edges also produces +bit-level different `stash_gain` values, contributing tie-shifts. + +The 2D problem stays at **exact** 0.000 diff because it has fewer ties +and fewer summed edges. The 3D problem has more opportunity for +divergence. + +**What this means for the bucket approach.** The bucket idea is sound +in the abstract — `chain_gain_init`'s data-volume floor really is O(E) +per outer iter, not O(pairs × pair_size × adjacency). But the +implementation creates a different `filtered_entries` ordering than the +adjacency-walk version, and that ordering difference is what produces +the 0.09 energy gap on 3D. To get bucket-level speed *and* pre-bucket +tie-breaking parity, you would need to either: + +1. Build `filtered_entries[v]` by walking `lifted_graph.node_adjacency(v)` + *after* using buckets to compute `stash_gain`. Adds back most of the + adjacency walk cost (estimated +30–50 ms), undoing most of the bucket + win. +2. Sort `filtered_entries[v]` by a canonical edge_id order after bucket + construction. Adds ~30 ms in per-pair sort cost; only partially + matches pre-bucket order because pre-bucket follows insertion order, + not edge_id order. +3. Accept the 0.17 diff. It's 0.005% relative on the only problem + where it shows, doesn't violate any test bound, and is well below + downstream segmentation noise. The 60 ms time win on the 3D problem + is real. + +I implemented option 3 (simplified buckets + mid-iter rebuild), and at +the user's request escalated to incremental maintenance hoping it would +also restore energy parity. It did not, and was slower besides. We are +back at the pre-bucket baseline. + +**Specific implementation lessons:** + +- `unordered_map>` for per-bucket storage caused + ~40 ms of cache-locality regression in `chain_gain_init` and + `chain_loop` vs the flat sorted vector. If revisiting buckets, + keep them in one contiguous flat array. +- The `relabel_node` machinery itself is cheap (~4 ms total over 10 + outer iters in the profile). Incremental maintenance is *not* + expensive in absolute terms; the slowdown came from the layout + switch. +- Predicted bucket-rebuild cost (notes said "~5 ms per outer iter") + was 4× off — sort on 40 k 40-byte entries actually takes ~11 ms. + The cheap-path radix-sort optimization listed below would address + this, but only matters if buckets come back. + +## Post-script: how the `changed_[]` flag closed the gap (2026-05-16) + +The bucket post-mortem (above) concluded with two recommended levers +(CSR adjacency, accept bucket tie-breaking) and the note that "nifty has +no algorithmic advantage" — both rooted in a careful read of +`lifted_twocut_kernighan_lin.hxx` (the per-pair two-cut routine). + +That read missed the outer driver. `lifted_multicut_kernighan_lin.hxx` +maintains a `changed_[]` flag per partition and gates the inner two-cut +on it: + +```cpp +if (!pV.empty() && (changed_[piU] || changed_[piV])) + twoCut_.optimizeTwoCut(pU, pV, twoCutBuffers_); ``` -Workload statistics on the 3D problem after the greedy warm-start: +The flag is refreshed each outer iter by `checkIfPartitonChanged()`, a +linear-time CC-style walk over base adjacency that marks a new partition +as "changed" iff its node-set differs from the previous iter's +partition that contained it (split or merge). The same gate is applied +to `introduceNewPartitions` (== our `cluster_splits`). -- 643 clusters. Sizes: min 1, max 139, mean 3.8, median 2. -- 4443 base cluster pairs per outer iter. Pair sizes: mean 20, median 7, - max 198 (long tail). -- Average lifted node degree: 32. - -## Ranked future optimizations +We added the equivalent: `detail_kl::compute_cluster_changed` (~60 lines) +plus gates in the `kernighan_lin()` driver. Result on 3D: -### 1. Pre-bucketed gain init (estimated landing point: ~150 ms total, −25%) +| Phase | Pre-flag | Post-flag | +|---|---|---| +| `pair_chains` | 167 ms | 61 ms | +| `chain_gain_init` | 78 ms | 28 ms | +| `chain_loop` | 85 ms | 31 ms | +| `cluster_splits` | 6.9 ms | 1.8 ms | +| `compute_changed` | — | 0.2 ms | +| total (profiled) | 354 ms | 135 ms | +| total (wall) | 195 ms | 86 ms | -**Idea.** For each outer iteration, bucket every lifted edge by its -endpoint cluster pair (sorted `(min_label, max_label)`). For pair-chain -`(A, B)` the gain init iterates only the three relevant buckets — -`(A, A)`, `(B, B)`, `(A, B)` — instead of walking the full lifted -adjacency of every node in the pair. +Energy stayed at 0.07 diff (vs 0.08 pre-flag — within noise); 2D stayed +at exact 0.000 diff. Tie-breaking is preserved by construction: iter 0 +has all partitions "changed" so it processes every pair (and every +split) — identical to today's algorithm. From iter 1, the only pairs +skipped are those where neither input changed since the previous iter, +where the chain is mathematically guaranteed to commit zero moves. -**Why it would help.** Each edge currently contributes to `chain_gain_init` -exactly once per pair-chain that touches one of its endpoint clusters. -With buckets, each edge contributes O(1) per outer iter. The 80 ms -`chain_gain_init` collapses to roughly O(E) = ~5 ms per outer iter, saving -~70 ms. +Why this works so well on the benchmark: the workload is 10 outer iters +on a problem that essentially converges after ~3 iters. Late iterations +were 80%+ wasted re-walking adjacency for partitions that hadn't moved. -**The wrinkle.** Bucket membership goes stale as nodes move between -clusters during the sequential pair-chains. The right fix is **incremental -bucket maintenance** — every node move performs O(degree) bucket-membership -updates (remove from old bucket, push into new). Bucket entries store -their position via an `edge_id → bucket_pos` index so removal is -`swap-with-back` in O(1). Total maintenance cost per outer iter: -`O(num_moves × avg_degree)` ≈ 30 µs on our problem. +## Future optimizations -**Complexity.** Substantial: ~200 lines, new `LiftedEdgeBuckets` struct in -`detail_kl`, hooks in `chain_loop` to call `buckets.relabel(v, old, new)` -on every committed move, careful invariant management. +The KL solver now beats nifty by ~17% on 3D and matches it on 2D. +Further optimization is not currently a priority. Sketched levers, in +case the workload changes: -**Sanity check before implementing.** Try a quick prototype where buckets -are rebuilt fresh at the start of each outer iter (O(E) per outer iter) -and gain init uses buckets, accepting that within-outer-iter moves create -stale entries. If the resulting energy is comparable to the exact version, -the incremental maintenance is worthwhile. +### CSR adjacency layout (estimated: 10–15% off `pair_chains`) -### 2. CSR adjacency layout for lifted graph (estimated: 10–15% off) +`UndirectedGraph` stores adjacency as `vector>` — one +heap allocation per node. A flat CSR built once at the start of +`kernighan_lin` would reduce per-node cache-line misses in +`chain_gain_init` and `chain_loop`. Doesn't change algorithm output. -**Idea.** `UndirectedGraph` stores adjacency as `vector>` -— one heap allocation per node. Walking adjacency for many distinct nodes -in a pair pays per-node pointer chasing. A flat CSR (`offsets[n+1]` + -`entries[2E]`) built once at the start of `kernighan_lin` would be more -cache-friendly. +### Skip pair-chains where heap stays empty (estimated: <5 ms) -**Caveat.** The actual access pattern in `chain_gain_init` walks -adjacency for nodes in arbitrary order (whichever pair we're processing), -so spatial locality across nodes is poor regardless of layout. The win is -limited to per-node cache-line savings (one miss per node vs one per -adjacency vector header). Estimate ~10% based on rough cycle counting. +Already implicitly handled by `heap.empty()` check, but `chain_init` +and `chain_gain_init` still run. Track live cluster sizes incrementally +to skip earlier when one side is empty/singleton. -**Complexity.** Localized: ~50 lines, build CSR in `kernighan_lin`, -replace `lifted_graph.node_adjacency(v)` calls in `chain_gain_init` with -CSR iteration. +### Bucket gain init -**Worth combining** with optimization 1, since CSR walking is what bucket -maintenance would need anyway. +Documented in the post-mortem (above). Would give ~60 ms on the +pre-flag baseline but only ~30 ms on the post-flag baseline (most of +the `chain_gain_init` time is now in iter 0 work, which buckets would +still cover). Same tie-breaking concern stands; not recommended without +a use-case that tolerates the 0.17 energy diff. -### 3. Inspect nifty's internals (estimated: unknown, possibly clarifying) +### Inspect nifty's internals — DONE 2026-05-16 -The gap between our `chain_gain_init` and nifty's equivalent -`computeDifferences` is suspiciously ~2× per adjacency entry given that -both algorithms walk the same data. nifty's source is at: +Read of `lifted_twocut_kernighan_lin.hxx` confirmed nifty's per-pair +two-cut has no algorithmic advantage over ours. The outer-driver +`changed_[]` flag (in `lifted_multicut_kernighan_lin.hxx`) was the +missing piece — now implemented. Source at: ``` +/home/pape/Work/software/src/nifty/include/nifty/graph/opt/lifted_multicut/lifted_multicut_kernighan_lin.hxx /home/pape/Work/software/src/nifty/include/nifty/graph/opt/lifted_multicut/detail/lifted_twocut_kernighan_lin.hxx ``` -Worth checking: -- How nifty stores `liftedGraph_.adjacency(v)` — is it CSR-like? -- Whether nifty's `referencedBy` array is `uint32` or larger (we use - `uint32`). -- Whether nifty's `differences` (= our `stash_gain`) cache line layout - differs. - -### 4. Linear-scan border instead of a heap - -**Verdict: not worth pursuing for this problem.** +### Linear-scan border instead of a heap -I worked through it: heap is faster than linear scan for our pair-size -distribution. The heap pop is O(log N) and heap.change is also O(log N); -for pair size 7 (the median) that's ~2× faster than O(N) linear scan, -and the gap widens for larger pairs. +**Verdict: not worth pursuing.** Heap is faster than linear scan for +our pair-size distribution; both `pop` and `change` are O(log N) on the +addressable indexed heap, beating O(N) linear scan for the median pair +size of 7. -nifty uses linear scan, but that's not where its speed advantage comes -from — likely it's the adjacency-walking constants (optimization 3). +## Workload statistics (3D problem, post greedy warm-start) -### 5. Skip pair-chains where heap stays empty (estimated: <5 ms) +Unchanged from before; useful for sanity-checking future estimates: -After `chain_gain_init`, if `heap.empty()` we already skip -`chain_loop`. But we still pay for `chain_init` and the full -`chain_gain_init` walk. For pair-chains where the cluster pair has only -one alive node per side (post-staleness filtering of `cluster_to_nodes`), -we could skip earlier. Need to maintain live cluster sizes. - -Low priority — only a few ms. +- 643 clusters. Sizes: min 1, max 139, mean 3.8, median 2. +- 4443 base cluster pairs per outer iter. Pair sizes: mean 20, + median 7, max 198 (long tail). +- Average lifted node degree: 32. ## Where to start next time 1. Re-run `cd development/graph/lifted_multicut && python check_kernighan_lin.py --size 3d --repeats 5` to confirm the baseline hasn't drifted. -2. Build with `pip install -e . --no-build-isolation -C cmake.define.BIOIMAGE_PROFILE=ON` to get the profile breakdown back. -3. Prototype the rebuild-per-outer-iter bucket approach (optimization 1 - simplified) to validate the energy quality before committing to the - incremental maintenance version. -4. Targets: - - 30%-of-nifty: ≤132 ms. - - Match nifty: ≤102 ms. +2. Build with `pip install -e . --no-build-isolation -C cmake.define.BIOIMAGE_PROFILE=ON` for the profile breakdown. +3. The next-most-attractive lever is CSR adjacency layout (10–15%), + but only worth doing if a heavier workload reveals the need. ## Files that matter diff --git a/development/graph/lifted_multicut/_compatibility.py b/development/graph/lifted_multicut/_compatibility.py index e173669..e267d38 100644 --- a/development/graph/lifted_multicut/_compatibility.py +++ b/development/graph/lifted_multicut/_compatibility.py @@ -28,6 +28,12 @@ def parser(description: str) -> argparse.ArgumentParser: default=60.0, help="Download timeout in seconds if the lifted problem is not cached.", ) + arg_parser.add_argument( + "--threads", + type=int, + default=1, + help="Number of threads for solvers that support it.", + ) arg_parser.add_argument( "--energy-bound", type=float, diff --git a/development/graph/lifted_multicut/check_fusion_move.py b/development/graph/lifted_multicut/check_fusion_move.py new file mode 100644 index 0000000..f29465c --- /dev/null +++ b/development/graph/lifted_multicut/check_fusion_move.py @@ -0,0 +1,64 @@ +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 lifted multicut at matched " + "settings." + ).parse_args() + # Match settings on both sides: + # - threads / parallel-proposals = `--threads N` (default 1). We pin + # P = T explicitly on the bic side so the comparison is + # apples-to-apples regardless of either library's default for P. + # - greedy-additive warm-start: we do it automatically when the + # objective's labels are the trivial singleton labeling. nifty's + # `fusionMoveBasedFactory` has no warm-start parameter, so we make it + # explicit on the nifty side by chaining greedy-additive in front. + # - greedy-additive sub-solver (the implicit default on both sides). + # - watershed proposal generator seeded from local edges. The bic + # proposal generator only sees the base graph, so seeding from + # "local" mirrors that exactly; the nifty default + # `SEED_FROM_LIFTED` would put seeds where bic cannot. + threads = int(args.threads) + # nifty's lifted-multicut fusion-move backend only implements single- + # threaded execution; passing numberOfThreads > 1 raises at solve time. + # Cap the nifty side at 1 so multi-threaded bic benchmarks still run, and + # surface the asymmetry in the printed output. + nifty_threads = 1 + if threads != nifty_threads: + print( + f"note: nifty lifted fusion-move runs single-threaded; comparing " + f"bic threads={threads} vs nifty threads={nifty_threads}." + ) + run_comparison( + "lifted_fusion_move", + lambda: bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(), + number_of_iterations=10, + stop_if_no_improvement=4, + number_of_threads=threads, + number_of_parallel_proposals=threads, + ), + lambda objective: objective.chainedSolversFactory( + [ + objective.liftedMulticutGreedyAdditiveFactory(), + objective.fusionMoveBasedFactory( + proposalGenerator=objective.watershedProposalGenerator( + seedingStrategy="SEED_FROM_LOCAL", + ), + numberOfIterations=10, + stopIfNoImprovement=4, + numberOfThreads=nifty_threads, + ), + ] + ), + args, + ) + + +if __name__ == "__main__": + main() diff --git a/development/graph/lifted_multicut/evaluate_solvers.py b/development/graph/lifted_multicut/evaluate_solvers.py index 8a0afeb..a3b2278 100644 --- a/development/graph/lifted_multicut/evaluate_solvers.py +++ b/development/graph/lifted_multicut/evaluate_solvers.py @@ -21,6 +21,12 @@ class SolverConfig: def solver_configs(): import bioimage_cpp as bic + # Fair single-threaded comparison: nifty's lifted fusion-move backend is + # single-threaded, so we pin the bic side to threads=1 and one parallel + # proposal per iteration (matching nifty's "generate one, fuse" loop) and + # chain greedy-additive in front of nifty to mirror bic's auto warm-start. + # Watershed seeding strategy is forced to SEED_FROM_LOCAL on the nifty + # side because the bic proposal generator only sees the base graph. return { "lifted_greedy_additive": SolverConfig( make_bic_solver=lambda: bic.graph.LiftedGreedyAdditiveMulticut(), @@ -39,6 +45,28 @@ def solver_configs(): ] ), ), + "lifted_fusion_move": SolverConfig( + make_bic_solver=lambda: bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(), + number_of_iterations=10, + stop_if_no_improvement=4, + number_of_threads=1, + number_of_parallel_proposals=1, + ), + make_nifty_factory=lambda objective: objective.chainedSolversFactory( + [ + objective.liftedMulticutGreedyAdditiveFactory(), + objective.fusionMoveBasedFactory( + proposalGenerator=objective.watershedProposalGenerator( + seedingStrategy="SEED_FROM_LOCAL", + ), + numberOfIterations=10, + stopIfNoImprovement=4, + numberOfThreads=1, + ), + ] + ), + ), } diff --git a/include/bioimage_cpp/graph/lifted_multicut.hxx b/include/bioimage_cpp/graph/lifted_multicut.hxx index d2720d8..a1f5392 100644 --- a/include/bioimage_cpp/graph/lifted_multicut.hxx +++ b/include/bioimage_cpp/graph/lifted_multicut.hxx @@ -1,5 +1,6 @@ #pragma once +#include "bioimage_cpp/graph/lifted_multicut/fusion_move.hxx" #include "bioimage_cpp/graph/lifted_multicut/greedy_additive.hxx" #include "bioimage_cpp/graph/lifted_multicut/kernighan_lin.hxx" #include "bioimage_cpp/graph/lifted_multicut/objective.hxx" diff --git a/include/bioimage_cpp/graph/lifted_multicut/fusion_move.hxx b/include/bioimage_cpp/graph/lifted_multicut/fusion_move.hxx new file mode 100644 index 0000000..5e791d2 --- /dev/null +++ b/include/bioimage_cpp/graph/lifted_multicut/fusion_move.hxx @@ -0,0 +1,438 @@ +#pragma once + +#include "bioimage_cpp/detail/edge_hash.hxx" +#include "bioimage_cpp/detail/profile.hxx" +#include "bioimage_cpp/detail/threading.hxx" +#include "bioimage_cpp/graph/detail/fusion_contract.hxx" +#include "bioimage_cpp/graph/lifted_multicut/greedy_additive.hxx" +#include "bioimage_cpp/graph/lifted_multicut/objective.hxx" +#include "bioimage_cpp/graph/proposal_generator.hxx" +#include "bioimage_cpp/graph/undirected_graph.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::graph::lifted_multicut { + +// Fusion-move solver for lifted multicut. Mirrors the multicut driver: +// stage-1 parallel proposal generation + parallel pairwise fuses, stage-2 +// sequential joint multi-fuse on leftover candidates, best-of safety net, +// optional warm-start from the trivial singleton labeling. +// +// The proposal generators emit node labelings over the *base* graph. Agreement +// contraction is also computed on the base graph only (lifted edges are not +// candidate contraction edges in lifted multicut). Both base and lifted +// weights are then aggregated onto the contracted subproblem: +// - A surviving base edge contributes its weight to the contracted base +// edge it maps onto. +// - A surviving lifted edge whose endpoints map to roots (ru, rv) that +// already share a contracted base edge contributes its weight to that +// base edge. +// - Otherwise the lifted edge becomes a new contracted lifted edge in the +// subproblem, with duplicates accumulated. +// +// The sub-solver is then a lifted-multicut solver over (contracted_base, +// contracted_lifted, contracted_weights). Default: lifted greedy-additive +// with a per-thread workspace reused across fuses. +class FusionMoveSolver final : public SolverBase { +public: + FusionMoveSolver( + std::vector proposal_generators, + const SolverBase *sub_solver = nullptr, + const std::size_t number_of_iterations = 10, + const std::size_t stop_if_no_improvement = 4, + const std::size_t number_of_threads = 1, + const std::size_t number_of_parallel_proposals = 2 + ) + : proposal_generators_(std::move(proposal_generators)), + sub_solver_(sub_solver), + number_of_iterations_(number_of_iterations), + stop_if_no_improvement_(stop_if_no_improvement), + number_of_threads_(number_of_threads), + number_of_parallel_proposals_(number_of_parallel_proposals) { + if (number_of_parallel_proposals < 1) { + throw std::invalid_argument( + "number_of_parallel_proposals must be >= 1" + ); + } + if (number_of_threads < 1) { + throw std::invalid_argument("number_of_threads must be >= 1"); + } + if (proposal_generators_.size() != number_of_parallel_proposals) { + throw std::invalid_argument( + "proposal_generators length must equal number_of_parallel_proposals" + ); + } + for (const auto *pgen : proposal_generators_) { + if (pgen == nullptr) { + throw std::invalid_argument("proposal_generators must not contain null"); + } + } + } + + std::vector optimize(Objective &objective) const override { + BIOIMAGE_PROFILE_INIT(profile); + const auto &base_graph = objective.graph(); + const auto &lifted_graph = objective.lifted_graph(); + const auto &weights = objective.weights(); + const auto n_base_edges = objective.number_of_base_edges(); + const auto number_of_nodes = base_graph.number_of_nodes(); + + std::vector current = objective.labels(); + if (number_of_nodes == 0 || lifted_graph.number_of_edges() == 0) { + objective.set_labels(current); + return objective.labels(); + } + + const auto effective_threads = ::bioimage_cpp::detail::normalize_thread_count( + number_of_threads_, number_of_parallel_proposals_ + ); + std::vector workspaces(effective_threads); + + if (is_singleton_labeling(current)) { + BIOIMAGE_PROFILE_SCOPE(profile, "warm_start"); + current = greedy_additive( + lifted_graph, weights, n_base_edges, + 0.0, -1.0, false, 42, 1.0, workspaces[0] + ); + } + + double current_energy; + { + BIOIMAGE_PROFILE_SCOPE(profile, "energy_eval"); + current_energy = energy(lifted_graph, weights, current); + } + + const std::size_t P = number_of_parallel_proposals_; + std::vector> proposal_buffers(P); + std::vector> fused_buffers(P); + std::vector proposal_energies(P); + std::vector fused_energies(P); + std::vector is_leftover(P); + + constexpr double kEnergyEps = 1e-7; + + std::size_t iterations_without_improvement = 0; + + for (std::size_t iteration = 0; iteration < number_of_iterations_; ++iteration) { + const auto ¤t_snapshot = current; + + std::fill(is_leftover.begin(), is_leftover.end(), 0); + + { + BIOIMAGE_PROFILE_SCOPE(profile, "proposal_and_pairwise_fuse"); + ::bioimage_cpp::detail::parallel_for_chunks( + effective_threads, + P, + [&](const std::size_t thread_id, const std::size_t begin, const std::size_t end) { + auto &workspace = workspaces[thread_id]; + for (std::size_t p = begin; p < end; ++p) { + proposal_generators_[p]->generate( + current_snapshot, proposal_buffers[p] + ); + proposal_energies[p] = energy( + lifted_graph, weights, proposal_buffers[p] + ); + + fuse_pair_into( + base_graph, + lifted_graph, + weights, + n_base_edges, + current_snapshot, + proposal_buffers[p], + sub_solver_, + workspace, + fused_buffers[p] + ); + fused_energies[p] = energy( + lifted_graph, weights, fused_buffers[p] + ); + } + } + ); + } + + double best_energy = current_energy; + const std::vector *best = ¤t; + std::size_t leftover_count = 0; + for (std::size_t p = 0; p < P; ++p) { + if (proposal_energies[p] < best_energy) { + best_energy = proposal_energies[p]; + best = &proposal_buffers[p]; + } + if (fused_energies[p] < best_energy) { + best_energy = fused_energies[p]; + best = &fused_buffers[p]; + } + if (fused_energies[p] > current_energy + kEnergyEps) { + is_leftover[p] = 1; + ++leftover_count; + } + } + + std::vector joint_result; + double joint_energy = std::numeric_limits::infinity(); + if (leftover_count >= 2) { + BIOIMAGE_PROFILE_SCOPE(profile, "joint_fuse"); + std::vector *> leftovers; + leftovers.reserve(leftover_count); + for (std::size_t p = 0; p < P; ++p) { + if (is_leftover[p]) { + leftovers.push_back(&fused_buffers[p]); + } + } + joint_result = fuse_multi( + base_graph, + lifted_graph, + weights, + n_base_edges, + leftovers, + sub_solver_, + workspaces[0], + profile + ); + { + BIOIMAGE_PROFILE_SCOPE(profile, "energy_eval"); + joint_energy = energy(lifted_graph, weights, joint_result); + } + if (joint_energy < best_energy) { + best_energy = joint_energy; + best = &joint_result; + } + } + + if (best_energy + kEnergyEps < current_energy) { + current = *best; + current_energy = best_energy; + iterations_without_improvement = 0; + } else { + ++iterations_without_improvement; + if (iterations_without_improvement >= stop_if_no_improvement_) { + break; + } + } + } + + objective.set_labels(current); + BIOIMAGE_PROFILE_REPORT(profile); + return objective.labels(); + } + +private: + static bool is_singleton_labeling(const std::vector &labels) { + for (std::size_t index = 0; index < labels.size(); ++index) { + if (labels[index] != static_cast(index)) { + return false; + } + } + return true; + } + + static void fuse_pair_into( + const UndirectedGraph &base_graph, + const UndirectedGraph &lifted_graph, + const std::vector &weights, + const std::uint64_t n_base_edges, + const std::vector ¤t, + const std::vector &proposal, + const SolverBase *sub_solver, + GreedyAdditiveWorkspace &workspace, + std::vector &output + ) { + const std::array *, 2> proposals{ + ¤t, &proposal + }; + std::vector *> proposal_list( + proposals.begin(), proposals.end() + ); + ::bioimage_cpp::detail::NullProfiler null_profile; + output = fuse_multi( + base_graph, lifted_graph, weights, n_base_edges, + proposal_list, sub_solver, workspace, null_profile + ); + } + + // Build the contracted lifted-multicut subproblem from N proposals, + // solve it, lift labels back. N=2 is the pairwise case; N>2 is the + // stage-2 joint fuse on leftovers. + template + static std::vector fuse_multi( + const UndirectedGraph &base_graph, + const UndirectedGraph &lifted_graph, + const std::vector &weights, + const std::uint64_t n_base_edges, + const std::vector *> &proposals, + const SolverBase *sub_solver, + GreedyAdditiveWorkspace &greedy_workspace, + [[maybe_unused]] ProfilerT &profile + ) { + const auto number_of_nodes = static_cast(base_graph.number_of_nodes()); + const auto n_proposals = proposals.size(); + + std::vector stacked(n_proposals * number_of_nodes); + for (std::size_t p = 0; p < n_proposals; ++p) { + std::copy( + proposals[p]->begin(), + proposals[p]->end(), + stacked.begin() + static_cast(p * number_of_nodes) + ); + } + + ::bioimage_cpp::graph::detail::AgreementContraction contraction; + { + BIOIMAGE_PROFILE_SCOPE(profile, "agreement_contract"); + contraction = ::bioimage_cpp::graph::detail::contract_by_agreement( + base_graph, stacked.data(), n_proposals, number_of_nodes + ); + } + + const auto &contracted_base = contraction.contracted_graph; + const auto n_contracted_base = + static_cast(contracted_base.number_of_edges()); + + // Aggregate base costs onto contracted base edges. + std::vector contracted_base_weights(n_contracted_base, 0.0); + { + BIOIMAGE_PROFILE_SCOPE(profile, "cost_aggregate_base"); + for (std::uint64_t edge = 0; edge < n_base_edges; ++edge) { + const auto target = contraction.contracted_edge_of_original[ + static_cast(edge) + ]; + if (target < 0) { + continue; + } + contracted_base_weights[static_cast(target)] += + weights[static_cast(edge)]; + } + } + + using ::bioimage_cpp::detail::Edge; + using ::bioimage_cpp::detail::EdgeHash; + using ::bioimage_cpp::detail::edge_key; + + // Look up which (ru, rv) pairs already exist as contracted base edges + // (so lifted edges between those roots fold into that base edge rather + // than introducing a new contracted lifted edge). The contracted base + // graph was built without `edge_lookup_`, so build our own map. + std::unordered_map base_lookup; + base_lookup.reserve(n_contracted_base); + for (std::uint64_t e = 0; e < contracted_base.number_of_edges(); ++e) { + const auto uv = contracted_base.uv(e); + base_lookup.emplace(uv, static_cast(e)); + } + + // Walk lifted edges. Three outcomes per edge: + // - endpoints in same root: dropped (never cut). + // - endpoints share a contracted base edge: accumulate onto base. + // - otherwise: emit/accumulate a new contracted lifted edge. + std::vector> new_lifted_uvs; + std::vector new_lifted_weights; + std::unordered_map lifted_lookup; + { + BIOIMAGE_PROFILE_SCOPE(profile, "cost_aggregate_lifted"); + const auto n_total_edges = lifted_graph.number_of_edges(); + for (std::uint64_t edge = n_base_edges; edge < n_total_edges; ++edge) { + const auto uv = lifted_graph.uv(edge); + const auto ru = contraction.root_of_node[static_cast(uv.first)]; + const auto rv = contraction.root_of_node[static_cast(uv.second)]; + if (ru == rv) { + continue; + } + const auto key = edge_key(ru, rv); + const auto base_it = base_lookup.find(key); + if (base_it != base_lookup.end()) { + contracted_base_weights[base_it->second] += + weights[static_cast(edge)]; + continue; + } + const auto lifted_it = lifted_lookup.find(key); + if (lifted_it != lifted_lookup.end()) { + new_lifted_weights[lifted_it->second] += + weights[static_cast(edge)]; + } else { + const auto index = new_lifted_uvs.size(); + new_lifted_uvs.emplace_back(key.first, key.second); + new_lifted_weights.push_back( + weights[static_cast(edge)] + ); + lifted_lookup.emplace(key, index); + } + } + } + + // Degenerate case: nothing left to optimize — all proposals already + // agreed everywhere that mattered. The agreement labeling itself is + // the answer. + if (n_contracted_base == 0 && new_lifted_uvs.empty()) { + std::vector result(number_of_nodes); + for (std::uint64_t node = 0; node < base_graph.number_of_nodes(); ++node) { + result[static_cast(node)] = contraction.root_of_node[ + static_cast(node) + ]; + } + return result; + } + + // Build the contracted Objective. The Objective constructor + // re-inserts base edges into its internal lifted graph and appends + // the new lifted edges; weights are stored as + // [base_costs..., lifted_costs...] in the same order. + std::vector sub_labels; + { + BIOIMAGE_PROFILE_SCOPE(profile, "sub_solve"); + Objective sub_objective( + contracted_base, + std::move(contracted_base_weights), + new_lifted_uvs, + new_lifted_weights, + false + ); + if (sub_solver == nullptr) { + sub_labels = greedy_additive( + sub_objective.lifted_graph(), + sub_objective.weights(), + sub_objective.number_of_base_edges(), + 0.0, + -1.0, + false, + 42, + 1.0, + greedy_workspace + ); + } else { + sub_labels = sub_solver->optimize(sub_objective); + } + } + + std::vector result(number_of_nodes); + { + BIOIMAGE_PROFILE_SCOPE(profile, "lift"); + for (std::uint64_t node = 0; node < base_graph.number_of_nodes(); ++node) { + const auto root = contraction.root_of_node[ + static_cast(node) + ]; + result[static_cast(node)] = sub_labels[ + static_cast(root) + ]; + } + } + return result; + } + + std::vector proposal_generators_; + const SolverBase *sub_solver_; + std::size_t number_of_iterations_; + std::size_t stop_if_no_improvement_; + std::size_t number_of_threads_; + std::size_t number_of_parallel_proposals_; +}; + +} // namespace bioimage_cpp::graph::lifted_multicut diff --git a/include/bioimage_cpp/graph/lifted_multicut/kernighan_lin.hxx b/include/bioimage_cpp/graph/lifted_multicut/kernighan_lin.hxx index b309816..952fed2 100644 --- a/include/bioimage_cpp/graph/lifted_multicut/kernighan_lin.hxx +++ b/include/bioimage_cpp/graph/lifted_multicut/kernighan_lin.hxx @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -339,6 +340,70 @@ inline double run_chain( return 0.0; } +// Mark which clusters in ``labels`` differ in node-set from any cluster in +// ``previous_labels``. A new cluster c is "unchanged" iff every node in c +// shares the same previous label c_old AND ``|c| == |c_old|`` — i.e. the +// partition was neither split nor merged during the last outer iter. +// +// Used to gate pair-chains and cluster-splits: a pair (A, B) whose inputs +// are identical to the last iteration's must produce the same chain result +// (which was "no improvement" — otherwise the partitions would have changed). +// Skipping such pairs is the only major algorithmic optimization in nifty's +// outer KL driver that we previously lacked. +inline std::vector compute_cluster_changed( + const std::vector &labels, + const std::vector &previous_labels, + const std::uint64_t number_of_clusters +) { + std::vector changed(static_cast(number_of_clusters), 0); + if (number_of_clusters == 0) { + return changed; + } + if (previous_labels.size() != labels.size()) { + // First iter or shape change — treat everything as changed. + std::fill(changed.begin(), changed.end(), 1); + return changed; + } + + const auto max_prev = *std::max_element(previous_labels.begin(), previous_labels.end()); + std::vector prev_size(static_cast(max_prev) + 1, 0); + for (const auto p : previous_labels) { + ++prev_size[static_cast(p)]; + } + + constexpr auto SENTINEL = std::numeric_limits::max(); + std::vector map_new_to_old( + static_cast(number_of_clusters), SENTINEL + ); + std::vector new_size(static_cast(number_of_clusters), 0); + + for (std::size_t v = 0; v < labels.size(); ++v) { + const auto c_new = static_cast(labels[v]); + const auto c_old = previous_labels[v]; + ++new_size[c_new]; + if (changed[c_new]) { + continue; + } + if (map_new_to_old[c_new] == SENTINEL) { + map_new_to_old[c_new] = c_old; + } else if (map_new_to_old[c_new] != c_old) { + changed[c_new] = 1; + } + } + + for (std::size_t c = 0; c < changed.size(); ++c) { + if (changed[c]) { + continue; + } + const auto c_old = map_new_to_old[c]; + if (c_old == SENTINEL + || prev_size[static_cast(c_old)] != new_size[c]) { + changed[c] = 1; + } + } + return changed; +} + // Re-split labels so that every cluster is connected in the base graph. // Returns the new labeling (dense relabeled). inline std::vector enforce_base_connectivity( @@ -392,6 +457,14 @@ inline std::vector kernighan_lin( detail_kl::ChainBuffers bufs(n_nodes); detail_kl::ChainScratch scratch(n_nodes); + // Per-cluster "changed since last outer iter" flag. Iter 0 processes + // every pair (no prior state); from iter 1 onward we skip (A, B) when + // neither side's node-set changed during the previous iter — the chain + // is a pure function of its inputs, so unchanged inputs replay the + // previous iter's zero-gain result. + std::vector changed; + std::vector prev_iter_labels = labels; + for (std::uint64_t iteration = 0; iteration < number_of_outer_iterations; ++iteration) { bool improved = false; @@ -407,9 +480,22 @@ inline std::vector kernighan_lin( cluster_to_nodes = detail_kl::build_cluster_to_nodes(labels, number_of_clusters); } + // Iter 0: all changed (no prior state). Iter k > 0: compare current + // labels (== end of iter k-1) against ``prev_iter_labels`` (== end of + // iter k-2). ``changed[c]`` answers "did partition c change during + // iter k-1?"; pairs/splits involving only stable partitions can be + // skipped. + if (iteration == 0) { + changed.assign(static_cast(number_of_clusters), 1); + } + { BIOIMAGE_PROFILE_SCOPE(profile, "pair_chains"); for (const auto &pair : pairs) { + if (!changed[static_cast(pair.a)] + && !changed[static_cast(pair.b)]) { + continue; + } const auto delta = detail_kl::run_chain( base_graph, lifted_graph, @@ -434,6 +520,9 @@ inline std::vector kernighan_lin( BIOIMAGE_PROFILE_SCOPE(profile, "cluster_splits"); std::uint64_t next_label = number_of_clusters; for (std::uint64_t cluster = 0; cluster < number_of_clusters; ++cluster) { + if (!changed[static_cast(cluster)]) { + continue; + } while (true) { if (next_label >= cluster_to_nodes.size()) { cluster_to_nodes.resize(static_cast(next_label) + 1); @@ -467,6 +556,20 @@ inline std::vector kernighan_lin( labels = dense_relabel(labels); } + // Recompute changed[] for use in the next outer iter, comparing the + // post-cc_repartition labels (end of this iter) against the labels + // at the start of this iter (``prev_iter_labels``). + { + BIOIMAGE_PROFILE_SCOPE(profile, "compute_changed"); + const auto new_n_clusters = labels.empty() + ? std::uint64_t{0} + : (*std::max_element(labels.begin(), labels.end()) + 1); + changed = detail_kl::compute_cluster_changed( + labels, prev_iter_labels, new_n_clusters + ); + prev_iter_labels = labels; + } + double new_energy; { BIOIMAGE_PROFILE_SCOPE(profile, "energy_eval"); diff --git a/include/bioimage_cpp/graph/proposal_generators/watershed.hxx b/include/bioimage_cpp/graph/proposal_generators/watershed.hxx index 40a5e2f..541c6f1 100644 --- a/include/bioimage_cpp/graph/proposal_generators/watershed.hxx +++ b/include/bioimage_cpp/graph/proposal_generators/watershed.hxx @@ -18,8 +18,12 @@ namespace bioimage_cpp::graph { // of negative-cost edges, then `edge_weighted_watershed`. Mirrors the workhorse // generator used by nifty's multicut fusion moves. // -// `n_seeds_fraction` is interpreted as a fraction of `number_of_nodes` when -// <= 1.0 and as an absolute seed-pair count otherwise. +// `n_seeds_fraction` is the target *total seed count* (not pair count): when +// <= 1.0 it is interpreted as a fraction of `number_of_nodes`, otherwise as +// an absolute count. Each iteration of the seeding loop places two seeds at +// the endpoints of one random negative-cost base edge, so the loop runs +// `n_seeds / 2` times. Matches nifty's `WatershedProposalGenerator` so that +// `n_seeds_fraction=0.1` produces the same proposal density on both sides. // // Scratch buffers (noisy costs, seeds, watershed sort buffer) are held as // members and reused across `generate` calls. @@ -69,21 +73,26 @@ public: noisy_costs_[edge] = static_cast(edge_costs_[edge] + noise_(generator_)); } - std::size_t n_seed_pairs; + std::size_t n_seeds_target; if (n_seeds_fraction_ <= 1.0) { - n_seed_pairs = static_cast( + n_seeds_target = static_cast( static_cast(number_of_nodes) * n_seeds_fraction_ + 0.5 ); } else { - n_seed_pairs = static_cast(n_seeds_fraction_ + 0.5); + n_seeds_target = static_cast(n_seeds_fraction_ + 0.5); } - n_seed_pairs = std::max(std::size_t{1}, n_seed_pairs); - n_seed_pairs = std::min(negative_edges_.size(), n_seed_pairs); + n_seeds_target = std::max(std::size_t{1}, n_seeds_target); + // Each loop iteration places two seeds, so the iteration count is + // half the target. Bottom of 1 keeps degenerate fractions usable. + const std::size_t n_seed_pair_iters = std::min( + negative_edges_.size(), + n_seeds_target == 1 ? std::size_t{1} : n_seeds_target / 2 + ); std::fill(seeds_buffer_.begin(), seeds_buffer_.end(), std::uint64_t{0}); std::uniform_int_distribution edge_dist(0, negative_edges_.size() - 1); std::uint64_t next_label = 1; - for (std::size_t i = 0; i < n_seed_pairs; ++i) { + for (std::size_t i = 0; i < n_seed_pair_iters; ++i) { const auto edge = negative_edges_[edge_dist(generator_)]; const auto uv = graph_.uv(edge); seeds_buffer_[static_cast(uv.first)] = next_label++; diff --git a/src/bindings/graph.cxx b/src/bindings/graph.cxx index 675c60b..245fdc8 100644 --- a/src/bindings/graph.cxx +++ b/src/bindings/graph.cxx @@ -9,6 +9,7 @@ #include "bioimage_cpp/graph/grid_graph.hxx" #include "bioimage_cpp/graph/lifted_from_affinities.hxx" #include "bioimage_cpp/graph/lifted_multicut.hxx" +#include "bioimage_cpp/graph/lifted_multicut/fusion_move.hxx" #include "bioimage_cpp/graph/multicut.hxx" #include "bioimage_cpp/graph/mutex_watershed.hxx" #include "bioimage_cpp/graph/multicut/fusion_move.hxx" @@ -858,6 +859,75 @@ UInt64Array multicut_fusion_move( return vector_to_uint64_array(result); } +UInt64Array lifted_multicut_fusion_move( + const Graph &base_graph, + const Graph &lifted_graph, + ConstDoubleArray weights, + const std::uint64_t n_base_edges, + ConstUInt64Array initial_labels, + std::vector proposal_generators, + const graph::lifted_multicut::SolverBase *sub_solver, + const std::size_t number_of_iterations, + const std::size_t stop_if_no_improvement, + const std::size_t number_of_threads, + const std::size_t number_of_parallel_proposals +) { + if (proposal_generators.empty()) { + throw std::invalid_argument("proposal_generators must not be empty"); + } + if (n_base_edges > lifted_graph.number_of_edges()) { + throw std::invalid_argument( + "n_base_edges must be <= lifted graph number_of_edges" + ); + } + auto weight_vector = + double_array_to_vector(weights, "edge_weights", lifted_graph.number_of_edges()); + auto label_vector = + uint64_array_to_vector(initial_labels, "initial_labels", base_graph.number_of_nodes()); + + // Decompose the user's lifted graph into base + lifted parts for the + // Objective constructor (which rebuilds the lifted graph internally). The + // overhead is one O(E) walk; the C++ Objective consumes only base costs + // and per-lifted-edge (u, v, weight) triples. + std::vector base_weights( + weight_vector.begin(), weight_vector.begin() + static_cast(n_base_edges) + ); + const auto n_lifted_edges = lifted_graph.number_of_edges() - n_base_edges; + std::vector> lifted_uvs; + lifted_uvs.reserve(static_cast(n_lifted_edges)); + std::vector lifted_weights; + lifted_weights.reserve(static_cast(n_lifted_edges)); + for (std::uint64_t edge = n_base_edges; edge < lifted_graph.number_of_edges(); ++edge) { + const auto uv = lifted_graph.uv(edge); + lifted_uvs.emplace_back(uv.first, uv.second); + lifted_weights.push_back(weight_vector[static_cast(edge)]); + } + + graph::lifted_multicut::FusionMoveSolver solver( + std::move(proposal_generators), + sub_solver, + number_of_iterations, + stop_if_no_improvement, + number_of_threads, + number_of_parallel_proposals + ); + + std::vector result; + { + nb::gil_scoped_release release; + graph::lifted_multicut::Objective objective( + base_graph, + std::move(base_weights), + lifted_uvs, + lifted_weights, + false + ); + objective.set_labels(std::move(label_vector)); + result = solver.optimize(objective); + } + return vector_to_uint64_array(result); +} + template Rag region_adjacency_graph_t( LabelArray labels, @@ -1513,6 +1583,22 @@ void bind_graph(nb::module_ &m) { nb::arg("number_of_parallel_proposals") ); + m.def( + "_lifted_multicut_fusion_move", + &lifted_multicut_fusion_move, + nb::arg("base_graph"), + nb::arg("lifted_graph"), + nb::arg("edge_weights"), + nb::arg("n_base_edges"), + nb::arg("initial_labels"), + nb::arg("proposal_generators"), + nb::arg("sub_solver").none(), + nb::arg("number_of_iterations"), + nb::arg("stop_if_no_improvement"), + nb::arg("number_of_threads"), + nb::arg("number_of_parallel_proposals") + ); + m.def( "_region_adjacency_graph_uint32", ®ion_adjacency_graph_t, diff --git a/src/bioimage_cpp/graph/__init__.py b/src/bioimage_cpp/graph/__init__.py index be0c11a..b0ed3ff 100644 --- a/src/bioimage_cpp/graph/__init__.py +++ b/src/bioimage_cpp/graph/__init__.py @@ -876,6 +876,12 @@ class WatershedProposalGenerator(ProposalGenerator): Per call: add Gaussian noise to edge costs, drop random seeds at the endpoints of negative-cost edges, run the edge-weighted watershed. + + ``n_seeds_fraction`` is the target *total seed count*: ``<= 1.0`` is a + fraction of ``number_of_nodes``, otherwise an absolute count. The + seeding loop places two seeds per iteration and runs ``n_seeds / 2`` + times, matching nifty's ``WatershedProposalGenerator`` so the same + parameter value yields the same proposal density on both sides. """ def __init__( @@ -1457,6 +1463,101 @@ def _build_cpp_sub_solver(self): ) +class FusionMoveLiftedMulticut(LiftedMulticutSolver): + """Fusion-move lifted multicut solver. + + Iteratively generates proposals via ``proposal_generator`` (which sees the + *base* graph and base edge costs), fuses them with the current best + labeling, and accepts improvements. Each fuse contracts the base graph by + agreement across the proposals, aggregates base + lifted weights onto the + contracted lifted-multicut subproblem, and dispatches to ``sub_solver``. + If ``sub_solver`` is omitted, the default sub-solver is + :class:`LiftedGreedyAdditiveMulticut`. + + If the objective's current labels are the trivial singleton labeling, the + driver warm-starts with one lifted greedy-additive pass before the proposal + loop. The best-of safety net guarantees energy never increases across + iterations. + + Threading: ``number_of_threads > 1`` runs ``number_of_parallel_proposals`` + proposal generators in parallel within each iteration. Each parallel slot + uses an independent proposal generator with seed ``proposal_generator.seed + + slot_index``. By default ``number_of_parallel_proposals`` is ``2`` when + ``number_of_threads == 1`` and ``number_of_threads`` otherwise; pass it + explicitly to override. + """ + + def __init__( + self, + *, + proposal_generator: ProposalGenerator, + sub_solver: LiftedMulticutSolver | None = None, + number_of_iterations: int = 10, + stop_if_no_improvement: int = 4, + number_of_threads: int = 1, + number_of_parallel_proposals: int | None = None, + ): + if not isinstance(proposal_generator, ProposalGenerator): + raise TypeError("proposal_generator must inherit from ProposalGenerator") + if sub_solver is not None and not isinstance(sub_solver, LiftedMulticutSolver): + raise TypeError("sub_solver must inherit from LiftedMulticutSolver") + if sub_solver is not None and not hasattr(sub_solver, "_build_cpp_sub_solver"): + raise TypeError( + "sub_solver must be a built-in lifted multicut solver " + "(custom Python solvers are not supported as fusion-move sub-solvers)" + ) + n_threads = int(number_of_threads) + if n_threads < 1: + raise ValueError("number_of_threads must be >= 1") + if number_of_parallel_proposals is None: + n_parallel = 2 if n_threads == 1 else n_threads + else: + n_parallel = int(number_of_parallel_proposals) + if n_parallel < 1: + raise ValueError("number_of_parallel_proposals must be >= 1") + + self.proposal_generator = proposal_generator + self.sub_solver = sub_solver + self.number_of_iterations = int(number_of_iterations) + self.stop_if_no_improvement = int(stop_if_no_improvement) + self.number_of_threads = n_threads + self.number_of_parallel_proposals = n_parallel + if self.number_of_iterations < 0: + raise ValueError("number_of_iterations must be non-negative") + if self.stop_if_no_improvement < 1: + raise ValueError("stop_if_no_improvement must be >= 1") + + def optimize(self, objective: LiftedMulticutObjective) -> np.ndarray: + n_base = objective.number_of_base_edges + # The base costs back the proposal generators (the lifted weights + # cannot drive base-graph contraction or watershed segmentation). + base_costs = np.ascontiguousarray(objective.weights[:n_base]) + cpp_pgens = [ + self.proposal_generator._build_for_thread( + objective.graph, base_costs, slot + ) + for slot in range(self.number_of_parallel_proposals) + ] + cpp_sub_solver = ( + None if self.sub_solver is None else self.sub_solver._build_cpp_sub_solver() + ) + labels = _core._lifted_multicut_fusion_move( + objective.graph, + objective.lifted_graph, + objective.weights, + n_base, + objective.labels, + cpp_pgens, + cpp_sub_solver, + self.number_of_iterations, + self.stop_if_no_improvement, + self.number_of_threads, + self.number_of_parallel_proposals, + ) + objective.labels = labels + return objective.labels + + class LiftedChainedSolvers(LiftedMulticutSolver): """Chain of lifted multicut solvers run in sequence on the same objective. @@ -1868,6 +1969,7 @@ def _normalize_number_of_threads(number_of_threads: int) -> int: "COMPLEX_EDGE_FEATURE_NAMES", "DEFAULT_EXTERNAL_MULTICUT_PROBLEM_PATH", "EXTERNAL_MULTICUT_PROBLEM_URL", + "FusionMoveLiftedMulticut", "FusionMoveMulticut", "GreedyAdditiveMulticut", "GreedyAdditiveProposalGenerator", diff --git a/tests/graph/lifted_multicut/test_fusion_move.py b/tests/graph/lifted_multicut/test_fusion_move.py new file mode 100644 index 0000000..e6c8105 --- /dev/null +++ b/tests/graph/lifted_multicut/test_fusion_move.py @@ -0,0 +1,311 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + +from ._helpers import same_partition + + +def _grid_lifted_problem(shape=(4, 4), bfs_distance=2, seed=0): + """A small grid base graph plus zero-weight BFS-induced lifted edges, + with some lifted edges then re-weighted negatively to make the problem + non-trivial.""" + rng = np.random.default_rng(seed) + edges = [] + base_costs = [] + for y in range(shape[0]): + for x in range(shape[1]): + node = y * shape[1] + x + if x + 1 < shape[1]: + edges.append([node, node + 1]) + base_costs.append(1.0 if x != shape[1] // 2 else -3.0) + if y + 1 < shape[0]: + edges.append([node, node + shape[1]]) + base_costs.append(1.0 if y != shape[0] // 2 else -3.0) + base = bic.graph.UndirectedGraph.from_edges(shape[0] * shape[1], edges) + base_costs = np.array(base_costs, dtype=np.float64) + + objective = bic.graph.LiftedMulticutObjective( + base, base_costs, bfs_distance=bfs_distance + ) + # Reweight a handful of lifted edges to non-zero, mixed sign, so they + # actually influence the energy. + n_lifted = objective.number_of_lifted_edges + n_base = objective.number_of_base_edges + weights = objective.weights.copy() + for i in range(min(n_lifted, 6)): + weights[n_base + i] = float(rng.normal(0.0, 2.0)) + objective._weights = weights + return objective + + +def test_fusion_move_splits_chain_along_repulsive_lifted_edge(chain_with_lifted): + base, base_costs, lifted_uvs, lifted_costs = chain_with_lifted + objective = bic.graph.LiftedMulticutObjective( + base, base_costs, lifted_uvs=lifted_uvs, lifted_costs=lifted_costs + ) + + # Singleton labels trigger the lifted greedy-additive warm start, which + # sees lifted weights and discovers the cut along the chain. The + # subsequent proposal/fuse loop must not regress past that. + merged_energy = objective.energy(np.zeros(4, dtype=np.uint64)) + + solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=1), + number_of_iterations=10, + ) + labels = solver.optimize(objective) + + assert labels.dtype == np.uint64 + assert labels.shape == (base.number_of_nodes,) + assert objective.energy(labels) <= merged_energy + 1e-9 + # Repulsive lifted edge must end up cut. + assert labels[0] != labels[3] + + +def test_fusion_move_safety_net_never_regresses(chain_with_lifted): + base, base_costs, lifted_uvs, lifted_costs = chain_with_lifted + objective = bic.graph.LiftedMulticutObjective( + base, base_costs, lifted_uvs=lifted_uvs, lifted_costs=lifted_costs + ) + baseline = bic.graph.LiftedGreedyAdditiveMulticut().optimize(objective) + baseline_energy = objective.energy(baseline) + + objective.reset_labels() + solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=7), + number_of_iterations=8, + ) + labels = solver.optimize(objective) + assert objective.energy(labels) <= baseline_energy + 1e-9 + + +def test_fusion_move_matches_multicut_without_lifted_edges(): + base = bic.graph.UndirectedGraph.from_edges( + 6, + [ + [0, 1], [0, 3], [1, 2], [1, 4], [2, 5], [3, 4], [4, 5], + ], + ) + base_costs = np.array([5, -20, 5, 5, -20, 5, 5], dtype=np.float64) + + mc_objective = bic.graph.MulticutObjective(base, base_costs) + mc_labels = bic.graph.FusionMoveMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=3), + number_of_iterations=5, + ).optimize(mc_objective) + + lmc_objective = bic.graph.LiftedMulticutObjective(base, base_costs) + lmc_labels = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=3), + number_of_iterations=5, + ).optimize(lmc_objective) + + # Energies should match (the lifted problem has no lifted edges, so it's + # equivalent to the base multicut problem). + assert mc_objective.energy(mc_labels) == pytest.approx( + lmc_objective.energy(lmc_labels), abs=1e-9 + ) + + +def test_fusion_move_warm_starts_from_singleton(): + base = bic.graph.UndirectedGraph.from_edges(3, [[0, 1], [1, 2]]) + base_costs = np.array([2.0, 2.0], dtype=np.float64) + objective = bic.graph.LiftedMulticutObjective(base, base_costs) + # Singleton labels trigger the internal lifted greedy-additive warm start. + labels = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=0), + number_of_iterations=3, + ).optimize(objective) + same_partition(labels, [0, 0, 0]) + + +def test_fusion_move_keeps_base_disconnected_clusters_separate( + disjoint_clusters_with_attractive_lifted, +): + base, base_costs, lifted_uvs, lifted_costs = ( + disjoint_clusters_with_attractive_lifted + ) + objective = bic.graph.LiftedMulticutObjective( + base, base_costs, lifted_uvs=lifted_uvs, lifted_costs=lifted_costs + ) + labels = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=0), + number_of_iterations=5, + ).optimize(objective) + # The (0, 2) lifted edge is attractive, but the base graph has no path + # between the two components — base-graph contraction never merges them. + assert labels[0] == labels[1] + assert labels[2] == labels[3] + + +def test_fusion_move_greedy_additive_proposal_generator_runs(): + objective = _grid_lifted_problem() + baseline = objective.energy( + bic.graph.LiftedGreedyAdditiveMulticut().optimize(objective) + ) + + objective.reset_labels() + solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.GreedyAdditiveProposalGenerator( + seed=3, sigma=1.0 + ), + number_of_iterations=5, + ) + labels = solver.optimize(objective) + assert objective.energy(labels) <= baseline + 1e-9 + + +@pytest.mark.parametrize( + "sub_solver", + [ + bic.graph.LiftedGreedyAdditiveMulticut(), + bic.graph.LiftedKernighanLinMulticut(number_of_outer_iterations=3), + ], +) +def test_fusion_move_sub_solver_pluggability(sub_solver): + objective = _grid_lifted_problem() + solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=2), + sub_solver=sub_solver, + number_of_iterations=4, + ) + labels = solver.optimize(objective) + assert labels.shape == (objective.graph.number_of_nodes,) + assert labels.dtype == np.uint64 + + +def test_fusion_move_stops_after_no_improvement(): + # Tiny problem with many iterations and an aggressive non-improvement + # threshold: the loop must terminate quickly. + base = bic.graph.UndirectedGraph.from_edges(3, [[0, 1], [1, 2], [0, 2]]) + base_costs = np.array([2.0, 2.0, -5.0], dtype=np.float64) + objective = bic.graph.LiftedMulticutObjective(base, base_costs) + baseline = objective.energy( + bic.graph.LiftedGreedyAdditiveMulticut().optimize(objective) + ) + + objective.reset_labels() + solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=0), + number_of_iterations=1000, + stop_if_no_improvement=1, + ) + labels = solver.optimize(objective) + assert objective.energy(labels) <= baseline + 1e-9 + + +def test_fusion_move_chains_with_kernighan_lin(chain_with_lifted): + base, base_costs, lifted_uvs, lifted_costs = chain_with_lifted + objective = bic.graph.LiftedMulticutObjective( + base, base_costs, lifted_uvs=lifted_uvs, lifted_costs=lifted_costs + ) + baseline = objective.energy( + bic.graph.LiftedGreedyAdditiveMulticut().optimize(objective) + ) + + objective.reset_labels() + solver = bic.graph.LiftedChainedSolvers([ + bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=11), + number_of_iterations=5, + ), + bic.graph.LiftedKernighanLinMulticut(number_of_outer_iterations=3), + ]) + labels = solver.optimize(objective) + assert objective.energy(labels) <= baseline + 1e-9 + + +def test_fusion_move_rejects_non_proposal_generator(): + with pytest.raises(TypeError, match="proposal_generator"): + bic.graph.FusionMoveLiftedMulticut(proposal_generator=object()) + + +def test_fusion_move_rejects_non_lifted_sub_solver(): + with pytest.raises(TypeError, match="sub_solver"): + bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(), + sub_solver=bic.graph.GreedyAdditiveMulticut(), + ) + + +def test_fusion_move_rejects_zero_thread_count(): + with pytest.raises(ValueError, match="number_of_threads"): + bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(), + number_of_threads=0, + ) + + +def test_fusion_move_rejects_zero_parallel_proposals(): + with pytest.raises(ValueError, match="number_of_parallel_proposals"): + bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(), + number_of_parallel_proposals=0, + ) + + +def test_fusion_move_runs_on_empty_graph(): + base = bic.graph.UndirectedGraph(0) + objective = bic.graph.LiftedMulticutObjective(base, np.zeros(0)) + solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(), + ) + labels = solver.optimize(objective) + assert labels.shape == (0,) + + +def test_fusion_move_parallel_threads_match_single_threaded_safety_net(): + objective = _grid_lifted_problem(seed=11) + baseline = objective.energy( + bic.graph.LiftedGreedyAdditiveMulticut().optimize(objective) + ) + + objective.reset_labels() + solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=11), + number_of_threads=4, + number_of_iterations=5, + ) + labels = solver.optimize(objective) + assert objective.energy(labels) <= baseline + 1e-9 + + +def test_fusion_move_multi_proposal_runs(): + objective = _grid_lifted_problem(seed=3) + baseline = objective.energy( + bic.graph.LiftedGreedyAdditiveMulticut().optimize(objective) + ) + + objective.reset_labels() + solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=3), + number_of_threads=2, + number_of_parallel_proposals=4, + number_of_iterations=5, + ) + labels = solver.optimize(objective) + assert objective.energy(labels) <= baseline + 1e-9 + + +def test_fusion_move_parallel_is_deterministic_given_settings(): + def run(): + objective = _grid_lifted_problem(seed=7) + solver = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=bic.graph.WatershedProposalGenerator(seed=7), + number_of_threads=4, + number_of_iterations=5, + ) + return solver.optimize(objective) + + np.testing.assert_array_equal(run(), run()) + + +def test_fusion_move_default_parallel_proposals_tracks_threads(): + pgen = bic.graph.WatershedProposalGenerator() + one_thread = bic.graph.FusionMoveLiftedMulticut(proposal_generator=pgen) + four_threads = bic.graph.FusionMoveLiftedMulticut( + proposal_generator=pgen, number_of_threads=4 + ) + assert one_thread.number_of_parallel_proposals == 2 + assert four_threads.number_of_parallel_proposals == 4