diff --git a/scripts/run_benchmark/param_sweep/baysor_params.yaml b/scripts/run_benchmark/param_sweep/baysor_params.yaml new file mode 100644 index 000000000..56088f369 --- /dev/null +++ b/scripts/run_benchmark/param_sweep/baysor_params.yaml @@ -0,0 +1,64 @@ +# Parameter sweep for the baysor transcript-assignment method. +# Committed source of truth, read at runtime from GitHub via a raw URL by +# run_test_baysor_nebius.sh (the Nebius compute env pulls the repo but cannot see the +# launch host's local files). +# +# Consumed by the run_benchmark workflow via the `method_parameters_yaml` setting +# (src/workflows/run_benchmark/main.nf). For every method the workflow builds: +# * one "default" variant using the `default:` args below, and +# * one extra variant per value in each `sweep:` list, with that ONE arg overridden. +# The benchmark varies a SINGLE parameter at a time (a "star" around the default, not +# a full grid), so total baysor variants = 1 default + sum(sweep list lengths) = 11. +# +# See src/methods_transcript_assignment/baysor/NOTES.md ("Optimization / tuning") for the +# tiers and rationale. +# +# WHAT BAYSOR ACTUALLY VARIES HERE: baysor runs Baysor v0.7.1 through sopa +# (sopa.segmentation.baysor). It re-segments transcripts starting from the upstream +# segmentation, fed in as a per-transcript prior. Coordinates are in ~1 um/px global space, so +# `scale` (Baysor's typical cell radius) is effectively in um. `scale = -1.0` is the sentinel +# for "estimate the radius from min_molecules_per_cell". +# +# NON-DEFAULT AUDIT (config default -> Baysor/sopa true default): +# * prior_segmentation_confidence: config 0.8 DEVIATES from Baysor default 0.2 (Baysor only +# recommends >0.7 when the prior is highly reliable) -> forced onto a sweep axis, walking +# back toward the tool default (0.2, and sopa's example 0.5). +# * min_molecules_per_cell: config 50 DEVIATES from sopa's example 10 (Baysor has no hard +# default) -> forced onto a sweep axis toward the smaller value. +# * scale: config -1.0 (auto-estimate) is a deliberate mode, not a tool default; swept to +# explicit um radii straddling sopa's example 6.25. +# * scale_std ("25%") and n_clusters (4) MATCH the Baysor defaults; probed lightly as Tier-2 +# companion axes (default value omitted from each list, as required). +# * force_2d (true) matches sopa's default and the data is 2D -> left fixed (flipping to 3D is +# meaningless here). transcripts_key / coordinate_system are structural, not swept. +# Perf/resource knobs (patch_width/patch_overlap, JULIA_NUM_THREADS) are hardcoded, not swept. +# +# SUBMITTABLE-NOW: the Nebius launch runs `--revision build/main`, and every swept arg below is +# ALREADY exposed in src/methods_transcript_assignment/baysor/config.vsh.yaml today, so the +# build/main container accepts them with NO rebuild. +parameters: + baysor: + # Baseline == the component's shipped defaults. + default: + force_2d: true + min_molecules_per_cell: 50 + scale: -1.0 + scale_std: "25%" + n_clusters: 4 + prior_segmentation_confidence: 0.8 + sweep: + # ---- Tier 1: cell-size / recall levers (sharpest) ---- + # scale: target cell radius in um (coords ~1 um/px). Default -1.0 = auto-estimate. + # Explicit radii straddling sopa's example 6.25 um; grounded in typical cell radii ~4-10 um. + scale: [5.0, 7.5, 10.0] + # min_molecules_per_cell: min molecules to call a cell; also drives the auto scale estimate. + # DEVIATES (config 50 vs sopa 10) -> walk toward the smaller, more permissive value. + min_molecules_per_cell: [10, 25] + # prior_segmentation_confidence: trust the prior vs re-segment (in [0,1]). + # DEVIATES (config 0.8 vs Baysor 0.2) -> walk back toward the tool default and sopa's 0.5. + prior_segmentation_confidence: [0.2, 0.5] + # ---- Tier 2: companion knobs at the Baysor default ---- + # scale_std: prior std of cell radius, relative to scale. Default "25%" -> wider tolerance. + scale_std: ["50%"] + # n_clusters: number of molecule clusters (~ major cell types). Default 4; Baysor rec. 3-15. + n_clusters: [6, 8] diff --git a/scripts/run_benchmark/param_sweep/clustermap_params.yaml b/scripts/run_benchmark/param_sweep/clustermap_params.yaml new file mode 100644 index 000000000..d434c00b1 --- /dev/null +++ b/scripts/run_benchmark/param_sweep/clustermap_params.yaml @@ -0,0 +1,75 @@ +# Parameter sweep for the clustermap (ClusterMap, He et al. 2021) transcript-assignment method. +# Committed source of truth, read at runtime from GitHub via a raw URL by +# run_test_clustermap_nebius.sh (the Nebius compute env pulls the repo but cannot see the +# launch host's local files). +# +# Consumed by the run_benchmark workflow via the `method_parameters_yaml` setting +# (src/workflows/run_benchmark/main.nf). For every method the workflow builds: +# * one "default" variant using the `default:` args below, and +# * one extra variant per value in each `sweep:` list, with that ONE arg overridden. +# The benchmark varies a SINGLE parameter at a time (a "star" around the default, not a full +# grid), so total clustermap variants = 1 default + sum(sweep list lengths) = 14. +# +# See src/methods_transcript_assignment/clustermap/NOTES.md ("Optimization / tuning") for the +# tiers, the ClusterMap-source audit, and the inert-arg findings. +# +# WHAT CLUSTERMAP ACTUALLY VARIES HERE: ClusterMap is segmentation-free — it density-peak- +# clusters the RNA spots (+ injected DAPI points) into cells. The swept knobs are the levers +# that reach the tool AND change the output, confirmed by reading script.py and the upstream +# ClusterMap source: +# xy_radius (spot-neighbourhood radius, px), cell_num_threshold (DPC cell-count threshold), +# dapi_grid_interval (DAPI-point density), min_spot_per_cell (small-cell filter), +# pct_filter (low-density spot removal), gauss_blur (DAPI blur on/off). +# +# NON-DEFAULT AUDIT vs ClusterMap true defaults: three config defaults DEVIATE from the tool +# and are therefore forced onto sweep axes (each list walks that arg back TOWARD its tool +# default, since the default variant already covers the shipped value): +# gauss_blur True -> tool False (sweep [false]) +# pct_filter 0.0 -> tool 0.1 (sweep toward 0.1) +# cell_num_threshold 0.1 -> tool 0.01 (sweep toward 0.01) +# A fourth deviation, `contamination` (0 -> tool 0.1), is DEAD CODE while LOF=False (the shipped +# default) and sklearn's LocalOutlierFactor rejects contamination=0, so LOF/contamination form a +# coupled, unsafe pair in the star model and are NOT swept (see NOTES.md). `window_size` (700) is +# a memory/perf tiling knob, not a ClusterMap parameter, so it is left fixed. +# +# INERT config args EXCLUDED from the sweep (they never reach ClusterMap): `use_dapi` (never read +# by script.py), `add_dapi` and `use_genedis` (both overwritten by par["xy_radius"] via a script +# bug). Sweeping any of them would just produce duplicate variants. Documented in NOTES.md. +# +# SUBMITTABLE-NOW: the Nebius launch runs `--revision build/main`, so the sweep touches ONLY args +# already exposed + correctly forwarded in the build/main container -- all six below exist in +# src/methods_transcript_assignment/clustermap/config.vsh.yaml today, so no rebuild is needed. +parameters: + clustermap: + # Baseline == the component's shipped defaults (the "default variant"). + default: + window_size: 700 + xy_radius: 40 + z_radius: 0 + fast_preprocess: false + gauss_blur: true + sigma: 1.0 + pct_filter: 0.0 + LOF: false + contamination: 0 + min_spot_per_cell: 5 + dapi_grid_interval: 5 + cell_num_threshold: 0.1 + sweep: + # ---- Tier 1: sharpest quality levers ---- + # xy_radius: spot-neighbourhood radius in PIXELS (~1 um/px grid). Default 40 is large vs + # typical cell radii; straddle tighter and looser neighbourhoods. + xy_radius: [20, 30, 60] + # cell_num_threshold: DPC cell-count threshold, "larger -> more cells". DEVIATION (0.1 vs + # tool 0.01); 0.01 = tool default (fewer/larger cells), 0.2 = more/smaller cells. + cell_num_threshold: [0.01, 0.05, 0.2] + # dapi_grid_interval: density of injected DAPI points; smaller = denser/slower. + dapi_grid_interval: [3, 10] + # min_spot_per_cell: min transcripts to keep a cell; higher deletes small clusters. + min_spot_per_cell: [10, 20] + # ---- Tier 2: quality knobs / must-sweep deviations ---- + # pct_filter: low-density spot removal. DEVIATION (0.0 removes nothing vs tool 0.1); walk up. + pct_filter: [0.05, 0.1] + # gauss_blur: DAPI blur before binarising. DEVIATION (True vs tool False); boolean -> the + # single non-default value is the tool default False. + gauss_blur: [false] diff --git a/scripts/run_benchmark/param_sweep/comseg_params.yaml b/scripts/run_benchmark/param_sweep/comseg_params.yaml new file mode 100644 index 000000000..2c655cec2 --- /dev/null +++ b/scripts/run_benchmark/param_sweep/comseg_params.yaml @@ -0,0 +1,74 @@ +# Parameter sweep for the comseg (ComSeg, Python/sopa) transcript-assignment method. +# Committed source of truth, read at runtime from GitHub via a raw URL by +# run_test_comseg_nebius.sh (the Nebius compute env pulls the repo but cannot see the +# launch host's local files). +# +# Consumed by the run_benchmark workflow via the `method_parameters_yaml` setting +# (src/workflows/run_benchmark/main.nf). For every method the workflow builds: +# * one "default" variant using the `default:` args below, and +# * one extra variant per value in each `sweep:` list, with that ONE arg overridden. +# The benchmark varies a SINGLE parameter at a time (a "star" around the default, not +# a full grid), so total comseg variants = 1 default + sum(sweep list lengths) = 13. +# +# See src/methods_transcript_assignment/comseg/NOTES.md ("Optimization / tuning") for the +# tiers and rationale. +# +# WHAT COMSEG ACTUALLY VARIES HERE: comseg runs sopa's ComSeg wrapper, which hands a config +# dict to the ComSeg point-cloud clustering algorithm (Defard et al., Comms Biol 2024). The +# component exposes exactly six SCIENCE knobs of that dict (mean_cell_diameter, max_cell_radius, +# alpha, min_rna_per_cell, norm_vector, allow_disconnected_polygon); the other config knobs +# (patch_width/overlap/transcript_patch_width, n_workers, worker_memory_limit) are memory/perf +# and resource knobs, and gene_column/transcripts_key/coordinate_system are structural, so none +# of those is swept. +# +# NON-DEFAULT AUDIT (vs sopa's own _get_default_config, the tool baseline): TWO component +# defaults deviate from the tool default and neither is a perf/resource knob, so each is FORCED +# onto a sweep axis -- +# * allow_disconnected_polygon: component TRUE vs sopa False -> sweep the flip (false). +# * min_rna_per_cell: component 5 vs sopa 20 -> sweep walks 5 back toward 20. +# The other four (mean_cell_diameter=15, max_cell_radius=25, alpha=0.5, norm_vector=false) sit at +# or near sopa's defaults but are the algorithm's core science levers, so they are swept too. +# +# SUBMITTABLE-NOW: the Nebius launch runs `--revision build/main`, so the sweep touches ONLY +# args already exposed in the build/main container -- all six knobs below already exist in +# src/methods_transcript_assignment/comseg/config.vsh.yaml today, so NO rebuild is needed. +# ComSeg's deeper clustering knobs (co-expression graph n_neighbors, k_nearest_neighbors, +# leiden resolution, min_nb_rna_patch) are NOT plumbed through the component config dict at all +# (ComSeg/sopa use their own defaults); exposing them would need a config arg + container +# rebuild, so they are documented in NOTES.md as future work and NOT here. +parameters: + comseg: + # Baseline == the component's shipped defaults (the six exposed ComSeg science knobs). + default: + mean_cell_diameter: 15.0 + max_cell_radius: 25.0 + alpha: 0.5 + min_rna_per_cell: 5 + norm_vector: false + allow_disconnected_polygon: true + sweep: + # ---- Tier 1: cell-size geometry (µm), the sharpest levers on what gets grouped ---- + # mean_cell_diameter: expected mean cell diameter in µm; sets the graph scale ComSeg + # builds neighborhoods on. Shipped 15µm (a typical mammalian cell). Straddle smaller / + # larger tissue cells. (sopa's true default DERIVES this from the prior-segmentation + # areas; the component instead pins 15µm.) + mean_cell_diameter: [10.0, 20.0] + # max_cell_radius: max centroid->RNA distance (µm) for association; caps how far a cell + # can reach out to claim transcripts. Shipped 25µm; sopa's rule-of-thumb is + # mean_cell_diameter*1.75 (~26µm at 15µm), the ComSeg docs example uses 50µm. Walk from + # tighter (20) up toward the docs value (50). + max_cell_radius: [20.0, 37.5, 50.0] + # alpha: alphashape parameter (0..1) for the cell polygon; 1 == convex hull, lower == + # tighter/more concave boundary hugging the points. Shipped 0.5. Sweep tighter, looser, + # and full convex hull. + alpha: [0.25, 0.75, 1.0] + # min_rna_per_cell: minimum transcripts for a cell to be kept. Shipped 5 (relaxed); + # sopa's default is 20. NON-DEFAULT AUDIT -> walk 5 back toward the sopa default 20 + # (stricter = fewer, more confident cells). + min_rna_per_cell: [10, 20] + # norm_vector (bool): normalize the per-cell expression vectors before the co-expression + # graph / clustering. Shipped false (== sopa default). Boolean -> one flip. + norm_vector: [true] + # allow_disconnected_polygon (bool): allow a cell's boundary to be a multi-part polygon. + # Shipped true; sopa's default is false. NON-DEFAULT AUDIT -> sweep the flip (false). + allow_disconnected_polygon: [false] diff --git a/scripts/run_benchmark/param_sweep/fastreseg_params.yaml b/scripts/run_benchmark/param_sweep/fastreseg_params.yaml new file mode 100644 index 000000000..119741064 --- /dev/null +++ b/scripts/run_benchmark/param_sweep/fastreseg_params.yaml @@ -0,0 +1,65 @@ +# Parameter sweep for the fastreseg (FastReseg, R) transcript-assignment method. +# Committed source of truth, read at runtime from GitHub via a raw URL by +# run_test_fastreseg_nebius.sh (the Nebius compute env pulls the repo but cannot see the +# launch host's local files). +# +# Consumed by the run_benchmark workflow via the `method_parameters_yaml` setting +# (src/workflows/run_benchmark/main.nf). For every method the workflow builds: +# * one "default" variant using the `default:` args below, and +# * one extra variant per value in each `sweep:` list, with that ONE arg overridden. +# The benchmark varies a SINGLE parameter at a time (a "star" around the default, not +# a full grid), so total fastreseg variants = 1 default + sum(sweep list lengths) = 11. +# +# See src/methods_transcript_assignment/fastreseg/NOTES.md ("Optimization / tuning") for the +# tiers and rationale. +# +# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +# NEEDS REBUILD (unlike the other six sibling sweeps, this one is NOT submittable as-is): +# these four knobs did NOT exist as Viash args until this change — they were HARDCODED inside +# script.R's fastReseg_full_pipeline(...) call. Exposing them adds four `--...` args to +# config.vsh.yaml + wires them through orchestrator.sh -> script.R. The `build/main` container +# run by run_test_fastreseg_nebius.sh's `--revision build/main` will REJECT these args until +# the component is rebuilt. Before launching you MUST: +# 1. commit + push this file AND the config/orchestrator/script.R edits, +# 2. `viash ns build` (regenerate target/) and rebuild the fastreseg container on ghcr +# with the build_main tag (see the check-component skill), then +# 3. launch. Otherwise every fastreseg variant fails with an "unknown option" error. +# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +# +# NON-DEFAULT AUDIT: every hardcoded value in script.R exactly matched FastReseg's TRUE package +# default for fastReseg_full_pipeline() (pixel_size 0.18, zstep_size 0.8, +# molecular_distance_cutoff 2.7, flagCell_lrtest_cutoff 5, svmClass_score_cutoff -2, +# groupTranscripts_method "dbscan", spatialMergeCheck_method "leidenCut", cutoff_spatialMerge +# 0.5, all NULL cutoffs). So NONE was a pre-tuned deviation demanding a sweep axis; the config +# defaults below are set to those same package defaults and the sweep walks each knob AWAY from +# its default in both directions. +# +# NOT SWEPT (documented in NOTES.md): pixel_size / zstep_size are Tier-0 dataset-physical +# geometry constants (still hardcoded, not exposed); groupTranscripts_method / +# spatialMergeCheck_method are Tier-2/3 categorical toggles (hardcoded). The four swept below +# are the Tier-1 error-detection / correction levers. +parameters: + fastreseg: + # Baseline == FastReseg's own package defaults (previously hardcoded in script.R). + default: + molecular_distance_cutoff: 2.7 + flagCell_lrtest_cutoff: 5 + svmClass_score_cutoff: -2 + cutoff_spatialMerge: 0.5 + sweep: + # ---- Tier 1: molecular neighborhood distance (SHARPEST lever) ---- + # Max molecule-to-molecule distance (um) for grouping transcripts into candidate cells. + # Default 2.7; smaller => tighter/more fragmented groups, larger => looser groups that + # may merge neighbouring cells. Widest sweep (4 values straddling 2.7). + molecular_distance_cutoff: [1.5, 2.0, 3.5, 4.0] + # ---- Tier 1: cell-flagging aggressiveness ---- + # lrtest_nlog10P cutoff to flag putative mis-segmented cells. Default 5; LOWER flags more + # cells for re-segmentation (aggressive correction), HIGHER is conservative. + flagCell_lrtest_cutoff: [3, 8] + # ---- Tier 1: transcript high/low-score class boundary (SVM) ---- + # Transcript-score cutoff separating high vs low score classes. Default -2; straddle it. + svmClass_score_cutoff: [-3, -1] + # ---- Tier 1: spatial-merge constraint ---- + # Fractional (0-1) spatial constraint for accepting a group-merge event. Default 0.5; + # lower is stricter (fewer merges), higher is more permissive. + cutoff_spatialMerge: [0.3, 0.7] diff --git a/scripts/run_benchmark/param_sweep/pciseq_params.yaml b/scripts/run_benchmark/param_sweep/pciseq_params.yaml new file mode 100644 index 000000000..cd04d5799 --- /dev/null +++ b/scripts/run_benchmark/param_sweep/pciseq_params.yaml @@ -0,0 +1,66 @@ +# Parameter sweep for the pciseq (pciSeq, Qian et al. 2020) transcript-assignment method. +# Committed source of truth, read at runtime from GitHub via a raw URL by +# run_test_pciseq_nebius.sh (the Nebius compute env pulls the repo but cannot see the +# launch host's local files). +# +# Consumed by the run_benchmark workflow via the `method_parameters_yaml` setting +# (src/workflows/run_benchmark/main.nf). For every method the workflow builds: +# * one "default" variant using the `default:` args below, and +# * one extra variant per value in each `sweep:` list, with that ONE arg overridden. +# The benchmark varies a SINGLE parameter at a time (a "star" around the default, not +# a full grid), so total pciseq variants = 1 default + sum(sweep list lengths) = 14. +# +# See src/methods_transcript_assignment/pciseq/NOTES.md ("Optimization / tuning") for the +# tiers and rationale. +# +# WHAT PCISEQ ACTUALLY VARIES HERE: pciSeq is a Variational-Bayes probabilistic model that +# jointly cell-types cells and assigns each transcript ("spot") to a candidate parent cell, +# driven by an scRNA-seq reference. The component forwards nine of pciSeq's `opts` knobs +# (script.py L176-188); this sweep varies the five with the highest impact on WHICH spots go +# to WHICH cell (the recall/precision + prior levers). save_data is an IO/debug flag (left +# fixed). transcripts_key/coordinate_system are structural (not swept). +# +# NON-DEFAULT AUDIT: every one of the nine exposed `opts` defaults in config.vsh.yaml +# (max_iter, CellCallTolerance, rGene, Inefficiency, InsideCellBonus, MisreadDensity, +# SpotReg, nNeighbors, rSpot, save_data) MATCHES pciSeq's own config.py default EXACTLY +# (verified against acycliq/pciSeq@issue_14 pciSeq/config.py). So NO knob is forced onto a +# sweep axis by the "deviation is evidence" rule — the ranges below are chosen on first +# principles, each straddling pciSeq's default (which the "default" variant already covers, +# hence the default value is omitted from every sweep list). +# +# SUBMITTABLE AS-IS: the Nebius launch runs `--revision build/main`, and every arg below is +# ALREADY exposed in src/methods_transcript_assignment/pciseq/config.vsh.yaml today, so no +# `viash ns build` / container rebuild is needed. Tier-2 knobs that ARE exposed but deferred +# from this first (light) test run: SpotReg, rSpot, CellCallTolerance, max_iter. Tier-3 pciSeq +# priors NOT surfaced by the component (cell_radius, cell_type_prior, mean_gene_counts_per_class, +# mean_gene_counts_per_cell) would need a new config arg + rebuild — see NOTES.md. +parameters: + pciseq: + # Baseline == the component's shipped defaults, which equal pciSeq's own config.py defaults. + default: + InsideCellBonus: 2 + MisreadDensity: 0.00001 + nNeighbors: 3 + Inefficiency: 0.2 + rGene: 20 + sweep: + # ---- Tier 1: what gets assigned to which cell (sharpest recall/precision levers) ---- + # InsideCellBonus: log-likelihood boost for a spot INSIDE a cell boundary. Higher => trust + # the segmentation more (spots pulled into the enclosing cell); 0/low => trust it less + # (spots freer to reassign by expression). Straddle default 2 both ways. + InsideCellBonus: [0, 1, 4, 6] + # MisreadDensity: uniform background misread rate (per pixel). Higher => more spots called + # background/misread (sparser, higher precision); lower => almost everything assigned to a + # cell. Swept on a log scale around the 1e-5 default. + MisreadDensity: [0.000001, 0.0001, 0.001] + # nNeighbors: number of candidate parent cells considered per spot (a "super-neighbor" + # background is always added on top). Fewer => stricter local assignment; more => a spot can + # reach farther cells. Straddle default 3. + nNeighbors: [2, 6] + # ---- Tier 1/2: per-gene efficiency Gamma prior (in-situ capture efficiency) ---- + # Inefficiency: mean of the per-gene efficiency Gamma prior (observed/theoretical count + # ratio). Lower/higher shifts how many counts the model expects per cell. Straddle 0.2. + Inefficiency: [0.1, 0.5] + # rGene: controls the VARIANCE of that efficiency Gamma prior (higher rGene => tighter prior + # around Inefficiency; lower => more per-gene flexibility). Straddle default 20. + rGene: [10, 40] diff --git a/scripts/run_benchmark/param_sweep/proseg_params.yaml b/scripts/run_benchmark/param_sweep/proseg_params.yaml new file mode 100644 index 000000000..66893c116 --- /dev/null +++ b/scripts/run_benchmark/param_sweep/proseg_params.yaml @@ -0,0 +1,59 @@ +# Parameter sweep for the proseg (Proseg, Rust) transcript-assignment method. +# Committed source of truth, read at runtime from GitHub via a raw URL by +# run_test_proseg_nebius.sh (the Nebius compute env pulls the repo but cannot see the +# launch host's local files). +# +# Consumed by the run_benchmark workflow via the `method_parameters_yaml` setting +# (src/workflows/run_benchmark/main.nf). For every method the workflow builds: +# * one "default" variant using the `default:` args below, and +# * one extra variant per value in each `sweep:` list, with that ONE arg overridden. +# The benchmark varies a SINGLE parameter at a time (a "star" around the default, not +# a full grid), so total proseg variants = 1 default + sum(sweep list lengths) = 10. +# +# See src/methods_transcript_assignment/proseg/NOTES.md ("Optimization / tuning") for the +# tiers and rationale. +# +# !!! NEEDS-REBUILD !!! This sweep touches NEWLY EXPOSED config arguments. proseg's model +# knobs (cell_compactness, nuclear_reassignment_prob, diffusion_probability, ncomponents, +# diffusion_sigma_far/near) were HARDCODED in script.py and are exposed as config arguments +# in this branch. The Nebius launch runs `--revision build/main`, whose container does NOT +# yet know these args, so this sweep WILL FAIL until the component is regenerated +# (`viash ns build`) and its container rebuilt (see the check-component skill). Commit + +# push the config/script changes AND rebuild the build_main image before launching. +# +# WHAT PROSEG ACTUALLY VARIES HERE: proseg runs a probabilistic (Cellular-Potts-like) sampler +# over a voxel lattice, using the prior nuclear segmentation, and reassigns transcripts to +# cells while regularising cell shape. The sweep walks the strongest model knobs. +# +# NON-DEFAULT AUDIT: the five previously-exposed args (voxel_layers, samples, burnin_samples, +# voxel_size, burnin_voxel_size) all MATCH proseg's own CLI defaults -> none is forced onto a +# sweep axis. Of the six previously-hardcoded knobs, only cell_compactness DEVIATED from +# proseg's true default (script.py hardcoded 0.03 vs proseg 0.04). Per the audit rule the +# deviating value must be covered by a sweep: the config default is reset to proseg's 0.04 +# (the "default variant") and cell_compactness is swept over values straddling the old 0.03 +# and both directions. --nthreads is a pure performance knob (computed from meta['cpus']) and +# is left fixed, not swept. +parameters: + proseg: + # Baseline == proseg's own CLI defaults for every knob (note cell_compactness=0.04, + # proseg's true default, NOT the previously hardcoded 0.03). + default: + cell_compactness: 0.04 + nuclear_reassignment_prob: 0.2 + diffusion_probability: 0.2 + ncomponents: 10 + sweep: + # ---- Tier 1: transcript-reassignment / cell-shape levers (strongest) ---- + # cell_compactness: exponential prior on cell compactness; SMALLER = more compact cells. + # Straddles the old hardcoded 0.03 (more compact) and adds a looser 0.06. proseg default 0.04. + cell_compactness: [0.02, 0.03, 0.06] + # nuclear_reassignment_prob: prior prob. of moving a transcript off its nuclear-prior label. + # Lower trusts the prior segmentation more; higher reassigns more freely. proseg default 0.2. + nuclear_reassignment_prob: [0.05, 0.5] + # diffusion_probability: prior prob. a transcript diffused from its cell of origin; higher + # treats more transcripts as diffusion noise -> more aggressive repositioning. proseg default 0.2. + diffusion_probability: [0.1, 0.35] + # ---- Tier 2: expression-model capacity ---- + # ncomponents: mixture components for cellular gene expression; more = finer cell-state + # modelling at higher cost / overfitting risk. proseg default 10. + ncomponents: [5, 15] diff --git a/scripts/run_benchmark/param_sweep/segger_params.yaml b/scripts/run_benchmark/param_sweep/segger_params.yaml new file mode 100644 index 000000000..b42ecf950 --- /dev/null +++ b/scripts/run_benchmark/param_sweep/segger_params.yaml @@ -0,0 +1,53 @@ +# Parameter sweep for the segger (GPU, GNN) transcript-assignment method. +# Committed source of truth, read at runtime from GitHub via a raw URL by +# run_test_segger_nebius.sh (the Nebius compute env pulls the repo but cannot see the +# launch host's local files). +# +# Consumed by the run_benchmark workflow via the `method_parameters_yaml` setting +# (src/workflows/run_benchmark/main.nf). For every method the workflow builds: +# * one "default" variant using the `default:` args below, and +# * one extra variant per value in each `sweep:` list, with that ONE arg overridden. +# The benchmark varies a SINGLE parameter at a time (a "star" around the default, not +# a full grid), so total segger variants = 1 default + sum(sweep list lengths) = 8. +# The workflow also allows at most ONE non-default method across the pipeline at a time; +# segger is that one non-default method here (every other stage stays on its single default). +# +# See src/methods_transcript_assignment/segger/NOTES.md ("Optimization / tuning") for the +# tiers and rationale. +# +# NON-DEFAULT AUDIT: every one of segger's exposed knobs already ships at segger 0.1.0's OWN +# default (n_epochs=20, prediction_graph_buffer_ratio=0.05, prediction_mode=cell, +# node_representation_dim=128 -- verified against dpeerlab/segger@0233cf62). So NOTHING is a +# pre-existing deviation forced onto an axis; this sweep is genuine exploration AROUND segger's +# defaults. (The 0.5 in the NOTES arg-table was the OLD --prediction-expansion-ratio default; +# 0.1.0 renamed the flag to --prediction-graph-buffer-ratio and lowered its default to 0.05 == +# what we ship.) +# +# node_representation_dim is deliberately NOT swept: it is a PCA-ceiling upper bound / capacity +# knob (script.py auto-retries at segger's reported ceiling on small panels), not a quality axis. +# +# SUBMITTABLE-NOW: the Nebius launch runs `--revision build/main`, and every arg swept below is +# already exposed in src/methods_transcript_assignment/segger/config.vsh.yaml, so NO container +# rebuild is needed (segger's GPU/RAPIDS image is expensive to rebuild -- avoided on purpose). +# segger is GPU-only; scheduling is handled by the gpuh100 label + src/base/labels_nebius_test.config. +parameters: + segger: + # Baseline == the component's shipped defaults (== segger 0.1.0's own defaults). + default: + n_epochs: 20 + prediction_graph_buffer_ratio: 0.05 + prediction_mode: cell + node_representation_dim: 128 + sweep: + # ---- Tier 1: which polygon set drives the prediction graph ---- + # prediction_mode: default `cell`; sweep the other two of segger's three choices. + prediction_mode: [nucleus, uniform] + # ---- Tier 1: polygon buffer/expansion = the recall<->precision dial ---- + # prediction_graph_buffer_ratio: expand each cell polygon by this fraction of its + # equivalent radius. Larger = capture more surrounding transcripts (more recall, less + # precision). Default 0.05 (tight); walk UP. 0.5 == the old --prediction-expansion-ratio + # default, a natural upper anchor. + prediction_graph_buffer_ratio: [0.1, 0.25, 0.5] + # ---- Tier 2: GNN training epochs (quality vs. speed) ---- + # More epochs = better convergence at a linear GPU-time cost (60 ~ 3x default train time). + n_epochs: [40, 60] diff --git a/scripts/run_benchmark/param_sweep/singler_params.yaml b/scripts/run_benchmark/param_sweep/singler_params.yaml new file mode 100644 index 000000000..50649e37d --- /dev/null +++ b/scripts/run_benchmark/param_sweep/singler_params.yaml @@ -0,0 +1,43 @@ +# Parameter sweep for the singler (SingleR, singler-py reference-correlation +# labeller) cell-type-annotation method. Committed source of truth for +# run_test_singler_nebius.sh (read from GitHub via a raw URL, since the Nebius +# compute env pulls the repo but cannot see the launch host's local files). +# +# Consumed by the run_benchmark workflow via the `method_parameters_yaml` setting +# (src/workflows/run_benchmark/main.nf). For every method the workflow builds: +# * one "default" variant using the `default:` args below, and +# * one extra variant per value in each `sweep:` list, with that ONE arg overridden. +# The benchmark varies a SINGLE parameter at a time (a "star" around the default, not +# a full grid), so total singler variants = 1 default + sum(sweep list lengths) = 4. +# +# See src/methods_cell_type_annotation/singler/NOTES.md ("Optimization / tuning") for +# the tiers and rationale. +# +# WHY celltype_key is the axis: it is the ONLY exposed method-relevant argument, and the +# scRNA-seq reference is annotated at nested granularities (cell_type -> cell_type_level2 +# -> level3 -> level4). Selecting the granularity the classifier trains on is the single +# meaningful lever for a reference labeller (Tier 1). The real singler-py algorithm knobs +# (marker_method, num_de, quantile, use_fine_tune/fine_tune_threshold, aggregate) are NOT +# exposed by the component and are documented in NOTES.md as future work (they would need +# a config expose + script wiring + container rebuild before they could be swept). +# +# !!! SUBMITTABLE-BUT-NO-OP CAVEAT !!! +# celltype_key already exists in build/main's config, so this sweep LAUNCHES without a +# rebuild. BUT build/main's baked-in script.py hardcodes the reference label column +# (ref_labels = ...column("cell_type")) and never reads par['celltype_key'] -- so on the +# deployed build/main container all four variants below produce IDENTICAL annotations +# (a no-op). The sweep only measures anything after a one-line wiring fix +# (...column(par["celltype_key"])) AND a container rebuild. See NOTES.md. +parameters: + singler: + # Baseline == the component's shipped default (the primary annotation column). + default: + celltype_key: cell_type + sweep: + # celltype_key: which reference obs column supplies the training labels, i.e. the + # annotation GRANULARITY the SingleR classifier is trained/scored against. The + # reference carries cell_type (default, covered by the default variant) plus three + # progressively finer levels. Finer levels add classes (higher resolution) but make + # the per-label correlation harder to separate on a small iST panel -- the knee + # between resolution and accuracy is what this axis is meant to find (once wired). + celltype_key: [cell_type_level2, cell_type_level3, cell_type_level4] diff --git a/src/base/labels_nebius_test.config b/src/base/labels_nebius_test.config new file mode 100644 index 000000000..f2684ed0e --- /dev/null +++ b/src/base/labels_nebius_test.config @@ -0,0 +1,188 @@ +// labels_nebius_test.config — TEST-CROP variant of labels_nebius.config. +// +// For runs over the small test crops (resources_test), every resource tier is downgraded: +// small CPU / memory / time requests, a low memory cap, and NO hard node selectors on the +// CPU/mem tiers (a small request fits any node group, so the scheduler packs many test +// tasks per node instead of pinning each to a specific large group). GPU tiers keep their +// node selectors (GPU nodes are a dedicated group) but request less. Retry/errorStrategy, +// get_memory() clamp and tracing are identical to the production config. +// +// Use with: tw launch ... --config src/base/labels_nebius_test.config +// Do NOT use for full/production data — the memory tiers are far too small. + +def exitStrat(task, max_attempts = 3) { + println "Determining exit strategy for task (attempt '${task.attempt}', exit status '${task.exitStatus}')" + + // if the component failed 3 times, ignore the error so the workflow can continue + if (task.attempt >= 3) { + return 'ignore' + } + // when an aws spot instance is reclaimed, nextflow seems to use exit code 2147483647 + if (task.exitStatus == null || task.exitStatus <= -1 || task.exitStatus > 2100000000 || !(task.exitStatus.toString().isNumber())) { + return 'retry' + } + // if component failed, retry once + if (task.exitStatus == 1 && task.attempt < 2) { + return 'retry' + } + // if component ran out of memory, retry with more memory and disk + if (task.exitStatus in [137, 139] && task.attempt < max_attempts) { + return 'retry' + } + return 'ignore' +} + + +process { + + // Default disk space (test crops are small) + disk = 30.GB + + // Always pull the latest image digest so nodes never serve a stale cached image. + pod = [[imagePullPolicy: 'Always']] + + errorStrategy = { exitStrat(task) } + maxRetries = 3 + // Low cap for test: even the escalated attempt-3 stays small and schedulable anywhere. + maxMemory = 120.GB + + // ---- CPU (downgraded) ---- + withLabel: lowcpu { cpus = 2 } + withLabel: midcpu { cpus = 4 } + withLabel: highcpu { cpus = 8 } + + // ---- Memory (downgraded; NO nodeSelector so small requests schedule on any node group) ---- + withLabel: lowmem { + memory = { get_memory( 8.GB * task.attempt ) } + disk = 30.GB + } + withLabel: midmem { + memory = { get_memory( 16.GB * task.attempt ) } + disk = 40.GB + } + withLabel: highmem { + memory = { get_memory( 32.GB * task.attempt ) } + disk = 50.GB + } + withLabel: veryhighmem { + memory = { get_memory( 64.GB * task.attempt ) } + disk = 60.GB + } + + withLabel: lowsharedmem { + containerOptions = { workflow.containerEngine != 'singularity' ? "--shm-size ${String.format("%.0f",task.memory.mega * 0.05)}" : ""} + } + withLabel: midsharedmem { + containerOptions = { workflow.containerEngine != 'singularity' ? "--shm-size ${String.format("%.0f",task.memory.mega * 0.1)}" : ""} + } + withLabel: highsharedmem { + containerOptions = { workflow.containerEngine != 'singularity' ? "--shm-size ${String.format("%.0f",task.memory.mega * 0.25)}" : ""} + } + + // ---- GPU (keep node selectors — GPU is a dedicated group — but request less) ---- + withLabel: gpu { + cpus = 4 + accelerator = 1 + memory = 20.GB + disk = 50.GB + // runAsUser: 0 — see labels_nebius.config (rapidsai/base non-root user can't write the + // root-owned Nextflow task scratch dir otherwise). Harmless for already-root GPU images. + pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r'], [imagePullPolicy: 'Always'], [runAsUser: 0]] + containerOptions = { workflow.containerEngine == "singularity" ? '--nv': + ( workflow.containerEngine == "docker" ? '--gpus all': null ) } + } + withLabel: midgpu { + cpus = 4 + accelerator = 1 + memory = 20.GB + disk = 50.GB + pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r'], [imagePullPolicy: 'Always']] + containerOptions = { workflow.containerEngine == "singularity" ? '--nv': + ( workflow.containerEngine == "docker" ? '--gpus all': null ) } + } + withLabel: highgpu { + cpus = 4 + accelerator = 1 + memory = 20.GB + disk = 50.GB + pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r'], [imagePullPolicy: 'Always']] + containerOptions = { workflow.containerEngine == "singularity" ? '--nv': + ( workflow.containerEngine == "docker" ? '--gpus all': null ) } + } + withLabel: biggpu { + cpus = 4 + accelerator = 1 + memory = 20.GB + disk = 50.GB + pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00t775jb99svb7k5r'], [imagePullPolicy: 'Always']] + containerOptions = { workflow.containerEngine == "singularity" ? '--nv': + ( workflow.containerEngine == "docker" ? '--gpus all': null ) } + } + withLabel: gpuhighmem { + cpus = 8 + accelerator = 1 + memory = { [ 40.GB * task.attempt, 100.GB ].min() } + disk = 60.GB + pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00dhcgx1xqjskycvc'], [imagePullPolicy: 'Always'], [runAsUser: 0]] + containerOptions = { workflow.containerEngine == "singularity" ? '--nv': + ( workflow.containerEngine == "docker" ? '--gpus all': null ) } + } + withLabel: gpuh100 { + cpus = 8 + accelerator = 1 + memory = { [ 40.GB * task.attempt, 120.GB ].min() } + disk = 60.GB + // runAsUser: 0 + /dev/shm emptyDir — required by segger (see labels_nebius.config). + pod = [[nodeSelector: 'nebius.com/node-group-id=mk8snodegroup-e00jp7hyqr094tmy35'], [imagePullPolicy: 'Always'], [runAsUser: 0], [emptyDir: [medium: 'Memory', sizeLimit: '16Gi'], mountPath: '/dev/shm']] + containerOptions = { workflow.containerEngine == "singularity" ? '--nv': + ( workflow.containerEngine == "docker" ? '--gpus all': null ) } + } + + // ---- Time (downgraded — test crops finish fast) ---- + withLabel: hightime { time = 2.h } + withLabel: veryhightime { time = 3.h } + withLabel: veryveryhightime { time = 4.h } + + // publishStates needs a little disk + memory + withName:'.*publishStatesProc' { + memory = '8GB' + disk = '50GB' + } + + // similarity metric (downgraded) + withName: '.*similarity_process' { + memory = '24.GB' + disk = '50.GB' + } +} + +def get_memory(to_compare) { + if (!process.containsKey("maxMemory") || !process.maxMemory) { + return to_compare + } + + try { + if (process.containsKey("maxRetries") && process.maxRetries && task.attempt == (process.maxRetries as int)) { + return process.maxMemory + } + else if (to_compare.compareTo(process.maxMemory as nextflow.util.MemoryUnit) == 1) { + return process.maxMemory as nextflow.util.MemoryUnit + } + else { + return to_compare + } + } catch (all) { + println "Error processing memory resources. Please check that process.maxMemory '${process.maxMemory}' and process.maxRetries '${process.maxRetries}' are valid!" + System.exit(1) + } +} + +// set tracing file +trace { + enabled = true + overwrite = true + file = "${params.publish_dir}/trace.txt" +} + +aws.batch.maxSpotAttempts = 5 +google.batch.maxSpotAttempts = 5 diff --git a/src/data_processors/process_dataset/script.py b/src/data_processors/process_dataset/script.py index b311e70db..aeb2b63eb 100644 --- a/src/data_processors/process_dataset/script.py +++ b/src/data_processors/process_dataset/script.py @@ -172,6 +172,61 @@ def crop_shapes_by_global_xy(shapes, x0, x1, y0, y1): set_transformation(new, trans, set_all=True) return new +def rasterize_boundaries_to_labels(sdata, shapes_key="cell_boundaries", image_key="image"): + """Rasterize polygon cell boundaries into a ``cell_labels`` labels element. + + Some loaders (e.g. allen_brain_cell_atlas_merfish) provide the vendor + segmentation only as polygon shapes and never rasterize it into a label + image. The ``custom_segmentation`` method hard-requires ``labels["cell_labels"]``, + so synthesize it here from the boundaries when it is absent — done post-crop + (see call site) so only the retained region is rasterized. + + Mirrors the vizgen_merscope loader: ``sd.rasterize`` labels regions 1..N + positionally and caps a single pass at 65535 regions, so rasterize in chunks + and offset each chunk's labels past that; then promote to a multiscale pyramid + so downstream ``["scale0"]`` indexing on the copied segmentation works. + """ + import dask.array as da + from spatialdata.models import Labels2DModel + + # sd.rasterize(return_regions_as_labels=True) numbers regions 1..n positionally + # into a uint16 array, so a single pass encodes at most 65535 distinct cells. + UINT16_MAX = 65535 + img_extent = sd.get_extent(sdata[image_key]) + n_cells = len(sdata[shapes_key]) + n_iter = n_cells // UINT16_MAX + bool(n_cells % UINT16_MAX) + + rasterize_args = { + "min_coordinate": [int(img_extent["x"][0]), int(img_extent["y"][0])], + "max_coordinate": [int(img_extent["x"][1]), int(img_extent["y"][1])], + "target_coordinate_system": "global", + "target_unit_to_pixels": 1, + "return_regions_as_labels": True, + } + + if n_iter <= 1: + labels_image = sd.rasterize(sdata[shapes_key], ["x", "y"], **rasterize_args) + else: + combined = None + template = None + for i in range(n_iter): + start = i * UINT16_MAX + end = min((i + 1) * UINT16_MAX, n_cells) + chunk = sd.rasterize(sdata[shapes_key].iloc[start:end], ["x", "y"], **rasterize_args) + chunk_np = np.asarray(chunk.data) + if combined is None: + combined = chunk_np.astype("uint32") + template = chunk + else: + mask = chunk_np > 0 + combined[mask] = chunk_np[mask].astype("uint32") + start + labels_image = template.copy(data=da.from_array(combined, chunks=template.data.chunksize)) + + # rasterize tags the labels with a shape->category map that the Labels model + # does not expect; drop it before parsing (matches the vizgen loader). + labels_image.attrs.pop("label_index_to_category", None) + return Labels2DModel.parse(labels_image, scale_factors=[2, 2, 2, 2]) + def rechunk_sdata(sdata, CHUNK_SIZE=1024): """Rechunk the sdata to the given chunk size @@ -383,6 +438,15 @@ def subsample_adata_group_balanced(adata, group_key, n_samples, seed=0): else: sdata_output = sdata +# Synthesize cell_labels from polygon boundaries for loaders that only provide +# shapes (e.g. allen_brain_cell_atlas_merfish). Done here — post-crop — so we +# rasterize only the retained region (rasterizing whole-brain labels would OOM) +# and avoid re-running the expensive stitching loader. custom_segmentation +# requires labels["cell_labels"]; loaders that already provide it are untouched. +if "cell_labels" not in sdata_output.labels and "cell_boundaries" in sdata_output.shapes: + print("No cell_labels found; rasterizing cell_boundaries -> cell_labels", flush=True) + sdata_output["cell_labels"] = rasterize_boundaries_to_labels(sdata_output) + # Rechunk to uniform chunks before writing (NOTE: rechunking currently needed, # https://github.com/scverse/spatialdata/issues/929). Run unconditionally so # that uncropped datasets (e.g. 10x Atera, whose store has rectilinear chunk diff --git a/src/methods_segmentation/stardist/NOTES.md b/src/methods_segmentation/stardist/NOTES.md index 724ecdab7..cde39c4d8 100644 --- a/src/methods_segmentation/stardist/NOTES.md +++ b/src/methods_segmentation/stardist/NOTES.md @@ -48,31 +48,41 @@ detector**: the script feeds it `image[0]` only (see Tier 0 below). (prob=0.479071, nms=0.3 for that model) as the fallback thresholds. 4. **Percentile normalizer** (`:64-77`) — a csbdeep `Normalizer` subclass that min-max scales the image to its **1st / 99.8th percentiles** (`normalize_mi_ma`), the - recommended StarDist preprocessing but with fixed percentile bounds. `block_size` - and `context` are derived from the image width so a large panel is processed in - tiles; **`min_overlap` is derived from object size, not block size** (see step 6 and - the min_overlap gotcha below). + recommended StarDist preprocessing but with fixed percentile bounds. The **image + size then selects the segmentation path** (single-pass vs tiled — see step 6 and the + min_overlap gotcha below). 5. **Build eval-params** (`:85-89`) — collects the newly exposed tunables (`prob_thresh, nms_thresh, scale`) from `par`, **dropping any that are `None`**. A dropped key ⇒ `predict_instances` uses the model's own optimized value (so the no-args call is byte-for-byte the pre-tuning behaviour). Mirrors the cellposev4 eval-params pattern. -6. **Segment** (`:92-131`) — `model.predict_instances_big(image[0], axes='YX', - block_size=…, min_overlap=…, context=…, normalizer=…, **eval_params)`. - `predict_instances_big` splits the image into `block_size` blocks, calls - `predict_instances` on each (forwarding `**eval_params` unchanged — it only - overrides `axes/overlap_label/return_labels/return_predict`), and reassembles the - labels into global coordinates. `image[0]` = first channel → a single 2D plane. - **The stitching invariant is that every predicted object is smaller than - `min_overlap`** (an object bigger than the overlap can span a block seam and can't be - uniquely assigned → `RuntimeError: ...violates the assumption of being smaller than - 'min_overlap'`). So `min_overlap` is set from an **object-size** bound - (`max_object_diameter`, default **192 px**), *not* from `block_size` (the old - `block_size // 5.5` shrank it to 64 px on small panels while real blobs reached - ~110 px → crash). `block_size` is then grown if needed to satisfy - `min_overlap + 2*context < block_size`, and the call is wrapped in a **retry that - doubles `min_overlap` on that specific error** so a rare oversized blob self-heals - instead of failing the run. +6. **Segment** (`:91-141`) — **two paths chosen by image size** (`image[0]` = first + channel → a single 2D plane): + - **Fits (largest side ≤ `BIG_PX`=4096) → `model.predict_instances(...)`.** One pass, + **no block stitching**, so there is *no* `min_overlap`/block-geometry constraint and + objects of any size are fine. `n_tiles` only sub-tiles the **forward pass** to bound + GPU memory (that tiling has its own automatic context and no overlap requirement). + This is the path all current benchmark panels take. + - **Large whole-slide (largest side > 4096) → `model.predict_instances_big(...)`.** + Tiles into `block_size` blocks and stitches, under **two** constraints: (1) a + *stitching* invariant that every object be smaller than `min_overlap` (else + `RuntimeError: ...violates the assumption of being smaller than 'min_overlap'`), and + (2) a *block-geometry* one — per-axis stride is `size − (min_overlap + 2·context)` + and stardist's `Block.cover` **asserts** consecutive write-regions overlap by ≥ + `min_overlap`, which requires `block_size` **comfortably** larger than + `min_overlap + 2·context`. Here `block_size = BIG_PX = 4096 ≫ min_overlap(192) + + 2·context`, so the geometry holds. `min_overlap` is **object-size-based** + (`max_object_diameter`, default 192 px, measured in original px — `predict_instances` + undoes `scale`), and the call is wrapped in a **retry that doubles `min_overlap`** + on that `RuntimeError` so a rare oversized blob self-heals. + + Why the split: the old code *always* used `predict_instances_big` with + `block_size = image.shape[1] // 3` (~336 px on a ~1000 px panel), which (a) forced + tiling even on tiny images and (b) left only a ~16 px margin over `min_overlap + + 2·context` → `Block.cover` `AssertionError`. A too-small `min_overlap` (`block_size // + 5.5` = 64 px) had earlier caused the *stitching* `RuntimeError` on a 110 px blob. Both + bug classes only exist on the tiled path; single-pass `predict_instances` sidesteps + them entirely. 7. **Post-process** (`:100-104`) — `convert_to_lower_dtype` downcasts the label array to the smallest uint that holds `max label`; wrap as an `xarray.DataArray`, `Labels2DModel.parse` with the copied transform, store as @@ -129,10 +139,11 @@ wrong for a different `--model` whose optimized thresholds differ. **They need `viash ns build` + a container rebuild to take effect** (see `check-component`). `max_object_diameter` is a **geometry/robustness knob, not a quality knob** — it only -sizes `predict_instances_big`'s `min_overlap`; it does not change which pixels get -segmented. It is **not part of the quality sweep**; leave it at the 192 px default -unless you hit the min_overlap `RuntimeError` (the script also auto-doubles it), or your -nuclei are unusually large. +sizes `min_overlap` on the **tiled (>4096 px) path**; on the single-pass path it is +unused, and it never changes which pixels get segmented. It is **not part of the quality +sweep**; leave it at the 192 px default unless a large whole-slide image hits the +min_overlap `RuntimeError` (the script also auto-doubles it) or its nuclei are unusually +large. Not exposed: - `--n_tiles` — a pure GPU-memory tiling knob. `predict_instances_big` **already** tiles @@ -219,16 +230,23 @@ That is exactly the sweep encoded in `scripts/run_benchmark/stardist_params.yaml `viash ns build` + a container rebuild; a stale image silently ignores them. The sweep has **not yet been run end-to-end** with the new args — validated only by `viash config view` + a script `ast.parse`. -- **`min_overlap` must exceed the largest object.** `predict_instances_big`'s block - stitching asserts every object is smaller than `min_overlap`; a bigger object throws - `RuntimeError: ...violates the assumption of being smaller than 'min_overlap'`. The old - code tied it to `block_size // 5.5`, so on small/narrow panels it fell to 64 px while - real blobs reached ~110 px → crash. Now `min_overlap` is **object-size-based** - (`max_object_diameter`, default 192 px), `block_size` is grown to keep - `min_overlap + 2*context < block_size`, and a **retry doubles `min_overlap`** on that - error. `scale` is forwarded per-block via `**kwargs` but objects are measured in - original pixels (predict_instances undoes `scale`), so `min_overlap` is in original px - regardless of `scale`; still sanity-check masks at block seams for extreme `scale`. +- **The `predict_instances_big` tiling has TWO independent failure modes — which is why + small images now bypass it entirely** (single-pass `predict_instances`, see step 6): + 1. *Stitching* — `RuntimeError: ...violates the assumption of being smaller than + 'min_overlap'` when an object is bigger than `min_overlap`. The old + `min_overlap = block_size // 5.5` fell to 64 px on small panels while blobs reached + ~110 px. + 2. *Block geometry* — `AssertionError` in `Block.cover` (per-axis + `stride = size − (min_overlap + 2·context)`; consecutive write-regions must overlap + by ≥ `min_overlap`). Fires when `block_size` is only *marginally* above + `min_overlap + 2·context`. Over-correcting fix #1 to `min_overlap=192` with + `block_size = image//3 ≈ 336` left a 16 px margin → this assertion tripped. + + On the surviving tiled path both are avoided by construction: `block_size = 4096 ≫ + min_overlap(192) + 2·context`, and the retry doubles `min_overlap` for a rare huge + object. `scale` is forwarded per-block but objects are measured in original pixels + (`predict_instances` undoes `scale`), so `min_overlap` is scale-independent; still + sanity-check masks at block seams for extreme `scale`. - **Only channel 0 is segmented** (`image[0]`). Fine for single-channel iST morphology; StarDist2D has no multi-channel mode anyway. - **Whole image loaded into RAM** (`:53`) — full-res plane; big panels are why the label diff --git a/src/methods_segmentation/stardist/script.py b/src/methods_segmentation/stardist/script.py index 4014f2772..d3fd12dd4 100644 --- a/src/methods_segmentation/stardist/script.py +++ b/src/methods_segmentation/stardist/script.py @@ -77,7 +77,7 @@ def do_after(self): mi, ma = np.percentile(image, [1,99.8]) normalizer = MyNormalizer(mi, ma) -# Tunable knobs forwarded through predict_instances_big -> predict_instances. +# Tunable knobs forwarded to predict_instances / predict_instances_big. # A value left as None (i.e. omitted from par) means "use the model's own optimized # value": thresholds.json for prob_thresh/nms_thresh, no rescaling for scale. This # keeps the default (no-args) call identical to the pre-tuning behaviour. @@ -86,46 +86,67 @@ def do_after(self): for k in ("prob_thresh", "nms_thresh", "scale") if par.get(k) is not None } -print(f"predict_instances_big overrides: {eval_params}", flush=True) - -# predict_instances_big tiles the image and stitches the per-block predictions. Its -# stitching invariant is that EVERY predicted object is smaller than `min_overlap`; -# an object spanning a block seam that is larger than the overlap can't be assigned to -# a single block, which raises "Found object of shape (...), which violates the -# assumption of being smaller than 'min_overlap'". So `min_overlap` must be -# OBJECT-SIZE-based, not block-size-based — the old `block_size // 5.5` shrank the -# overlap below real nuclei/blobs on small panels (min_overlap fell to 64 px while -# objects reached ~110 px). Objects are measured in ORIGINAL image pixels -# (predict_instances undoes `scale` internally), so this bound is in original px and is -# independent of `scale`. `context` is only the receptive-field margin discarded around -# each block, so deriving it from the image size is fine. -block_size = min(image.shape[1] // 3, 4096) -context = int(min(block_size // 5.5, 128)) -min_overlap = int(par.get("max_object_diameter") or 192) # px; must exceed largest object -# predict_instances_big asserts: min_overlap + 2*context < block_size. -block_size = max(block_size, min_overlap + 2 * context + 1) - -# Self-heal: if a rare oversized blob (merged nuclei / debris) still exceeds -# `min_overlap`, double it (and grow block_size to keep the geometry constraint) and -# retry, rather than failing the whole segmentation. -while True: - try: - labels, _ = model.predict_instances_big( - image[0, :, :], axes='YX', block_size=block_size, - min_overlap=min_overlap, context=context, - normalizer=normalizer, **eval_params, # n_tiles left to block_size - ) - break - except RuntimeError as e: - if "min_overlap" not in str(e) or min_overlap >= 2048: - raise - min_overlap *= 2 - block_size = max(block_size, min_overlap + 2 * context + 1) - print( - "predict_instances_big: an object exceeded min_overlap; retrying with " - f"min_overlap={min_overlap}, block_size={block_size}", - flush=True, - ) +print(f"stardist overrides: {eval_params}", flush=True) + +# Segmentation strategy — two paths, chosen by image size: +# +# * predict_instances_big TILES the image and stitches the per-block predictions +# under TWO strict constraints. (1) a *stitching* invariant that every object be +# smaller than `min_overlap`; and (2) a *block-geometry* one, since the per-axis +# stride is `size - (min_overlap + 2*context)` and stardist's `Block.cover` +# asserts consecutive blocks' write-regions overlap by >= min_overlap — which +# needs `block_size` to be *comfortably* larger than `min_overlap + 2*context`, +# not just larger (the old `block_size = image.shape[1] // 3` made ~336 px blocks +# on a ~1000 px panel, leaving a 16 px margin, so Block.cover's assertion failed). +# * predict_instances processes the whole image in ONE pass: no block stitching, so +# NEITHER constraint exists and objects of any size are fine. `n_tiles` only +# sub-tiles the forward pass to bound GPU memory (that tiling has its own automatic +# context and no min_overlap requirement). +# +# So: if the image fits (largest side <= BIG_PX) use the constraint-free single-pass +# `predict_instances`; only genuinely large whole-slide images take the tiled path, +# where block_size=BIG_PX >> min_overlap+2*context keeps Block.cover's geometry valid. +# `min_overlap` (the tiled path only) is OBJECT-SIZE-based: it must exceed the largest +# object, measured in ORIGINAL image pixels (predict_instances undoes `scale`), so it is +# independent of `scale`. +def _n_tiles(px): + # Forward-pass tiling to bound GPU memory (~2048 px/tile); no stitching constraint. + n = max(1, int(px) // 2048) + return (n, n) + +BIG_PX = 4096 +max_dim = int(max(image.shape[1], image.shape[2])) + +if max_dim <= BIG_PX: + labels, _ = model.predict_instances( + image[0, :, :], axes='YX', normalizer=normalizer, + n_tiles=_n_tiles(max_dim), **eval_params, + ) +else: + # Large image -> tile + stitch. Self-heal: if a rare oversized blob (merged nuclei / + # debris) exceeds `min_overlap`, double it (and grow block_size to keep the geometry + # constraint) and retry, rather than failing the whole segmentation. + block_size = BIG_PX + min_overlap = int(par.get("max_object_diameter") or 192) # px; must exceed largest object + context = int(min(min_overlap, 128)) + while True: + try: + labels, _ = model.predict_instances_big( + image[0, :, :], axes='YX', block_size=block_size, + min_overlap=min_overlap, context=context, n_tiles=_n_tiles(block_size), + normalizer=normalizer, **eval_params, + ) + break + except RuntimeError as e: + if "min_overlap" not in str(e) or min_overlap >= 2048: + raise + min_overlap *= 2 + block_size = max(block_size, min_overlap + 2 * context + 256) + print( + "predict_instances_big: an object exceeded min_overlap; retrying with " + f"min_overlap={min_overlap}, block_size={block_size}", + flush=True, + ) diff --git a/src/methods_transcript_assignment/pciseq/script.py b/src/methods_transcript_assignment/pciseq/script.py index 264c61718..b3f67b96c 100644 --- a/src/methods_transcript_assignment/pciseq/script.py +++ b/src/methods_transcript_assignment/pciseq/script.py @@ -124,7 +124,14 @@ def eta_update_no_assert(self): #TODO this will immediately break when the name of the gene isn't feature_name # Materialize the full transcripts once, in the same row order as the transformed # x/y coordinates above, so every downstream filter stays positionally aligned. -transcripts_full = sdata[par['transcripts_key']].compute() +# .compute() collapses the multi-partition transcripts (each merscope partition is 0-indexed, +# so large datasets like the kuppe merscope have a globally non-unique index) into one frame. +# Every sibling transcript-assignment method resets the index here; pciSeq previously reset only +# the transform copy (transcripts_reset above) and built its OUTPUT from this frame, so the +# duplicate index reached PointsModel.parse and broke the final write with +# "cannot reindex on an axis with duplicate labels". reset_index preserves row order, so the +# x_coords/y_coords computed in the same .compute() order stay positionally aligned. +transcripts_full = sdata[par['transcripts_key']].compute().reset_index(drop=True) transcripts_dataframe = transcripts_full[['feature_name']].copy() transcripts_dataframe['x'] = x_coords transcripts_dataframe['y'] = y_coords