run_bagol(smld::SMLD; kwargs...) -> (BasicSMLD, BaGoLDiagnostics)

Run BaGoL analysis on an SMLD. Returns grouped emitters as BasicSMLD and diagnostics.

Uses precision-weighted DBSCAN to partition localizations, then runs parallel collapsed Gibbs MCMC on each partition with global hierarchical updates.

Count Model

n_j ~ Gamma(shape, μ/shape) where:

• μ = mean locs per emitter
• shape = 1: exponential (dSTORM), shape > 1: peaked (DNA-PAINT)

Partitioning Arguments

• partition_sigma=3.0: DBSCAN threshold in sigma units (Inf = no partitioning)
• min_partition_size=0: Minimum locs per partition (smaller clusters dropped as noise)
• max_partition_size=1000: Split partitions larger than this
• skip_partition_size=typemax(Int): Skip partitions larger than this

Uncertainty Correction (standalone — leave 0 in the integrated pipeline)

• se_adjust=0.0: Independent extra position error τ (μm) added in quadrature to the per-loc CRLB σ (σ²eff = σ² + τ²). Accepts a scalar (both axes), a (τx, τy) tuple (per-axis), a length-N vector (per-loc), or `(τxvec, τyvec). Corrects BaGoL under-grouping when CRLB underestimates the true localization error. If the input SMLD is already σ-corrected (metadata["sigmacorrected"]==true), se_adjust is skipped with a warning unlessforceseadjust=true`.
• force_se_adjust=false: Override the already-corrected guard (intentional double-count).

MCMC Arguments

• sync_interval=100: Iterations between global μ/shape updates
• n_iterations=4000: Total MCMC iterations
• burn_in=2000: Burn-in iterations before recording
• shape=2.0: Initial Gamma shape (1=exponential, higher=more peaked)
• learn_distribution=true: Control count distribution learning. true=learn both μ and shape, false=fix both, :mu=learn μ only (fix shape), :shape=learn shape only (fix μ)
• learn_rho=true: Learn Poisson-K emitter density ρ independently of μ/shape. For fully fixed hyperparameters, set both learn_distribution=false and learn_rho=false.
• rho=nothing: Initial ρ in emitters per unit area. If provided with learn_rho=false, hold this fixed for Fazel-style fixed-ρ-over-bounding-box flat Poisson-K runs.
• allocation_model=:dm: :dm (Dirichlet-Multinomial), :decoupled (no partition prior), or :categorical (labeled K^(-N) — Fazel-equivalent when paired with spatial_model=:flat + k_prior=:none)
• spatial_model=:locmix: :locmix (localization mixture) or :flat
• k_prior=:auto: K-prior gating. :auto uses spatial-model default (Poisson(ρA) under :flat, none under :locmix); :poisson always include (only valid with :flat); :none always disable
• verbose=true: Print progress

Posterior Image

• posterior_pixel_size=0.001: Rao-Blackwellized posterior image pixel size in μm (0.0 to disable)
• posterior_xlim=nothing: Override x bounds for posterior image
• posterior_ylim=nothing: Override y bounds for posterior image

Archive

• archive_path=nothing: Enable mmap chain archive for post-hoc analysis

Returns

• BasicSMLD: Grouped emitter positions with uncertainties
• BaGoLDiagnostics: nemitters, posteriork, acceptance_rates, final parameters

run_bagol(smld, cfg::BaGoLConfig; posterior_xlim=nothing, posterior_ylim=nothing)

Run BaGoL from a config struct. Runtime-only kwargs (posterior bounds) are separate.

run_bagol(locs::Vector{<:AbstractEmitter}; camera, kwargs...) -> (BasicSMLD, BaGoLDiagnostics)

Run BaGoL on a vector of localizations. Requires camera argument.

See run_bagol(smld::SMLD; ...) for full documentation.

BaGoLConfig <: SMLMData.AbstractSMLMConfig

Configuration for BaGoL analysis. All fields correspond 1:1 to run_bagol kwargs.

Count Model

• μ=10.0: Mean localizations per emitter
• shape=2.0: Gamma shape (1=exponential/dSTORM, >1=peaked/DNA-PAINT)
• learn_distribution=true: Learn count params. true=both, false=fix both, :mu=learn μ only, :shape=learn shape only
• gamma=nothing: DM concentration. nothing=use shape (default), Float64=fixed
• learn_rho=true: Learn Poisson-K emitter density ρ independently of μ/shape. For fully fixed hyperparameters, set both learn_distribution=false and learn_rho=false.
• rho=nothing: Initial ρ in emitters per unit area. If provided with learn_rho=false, hold this fixed for Fazel-style fixed-ρ-over-bounding-box flat Poisson-K runs.

MCMC

• n_iterations=4000: Total MCMC iterations
• burn_in=2000: Burn-in before accumulating
• sync_interval=100: Iterations between global hierarchical updates
• allocation_model=:dm: :dm (Dirichlet-Multinomial), :decoupled (no partition prior), or :categorical (labeled K^(-N) — Fazel-equivalent when paired with spatial_model=:flat + k_prior=:none)
• spatial_model=:locmix: :locmix (localization mixture) or :flat
• k_prior=:auto: K-prior gating. :auto uses spatial-model default (Poisson(ρA) under :flat, none under :locmix); :poisson always include (only valid with :flat); :none always disable
• n_restricted_scans=5: Jain-Neal restricted Gibbs scans for split/merge
• n_bd_substeps=5: Birth/death substeps per BD selection

Partitioning

• partition_sigma=3.0: DBSCAN threshold in sigma units (Inf=no partitioning)
• min_partition_size=0: Minimum locs per partition
• max_partition_size=1000: METIS-split partitions larger than this
• skip_partition_size=typemax(Int): Skip partitions larger than this
• overlap=:auto: Overlap width for bisected partitions (:auto or Float64 in μm)
• bridge_ratio=0.0: Optional bridge refinement strength for DBSCAN clusters
• min_split_size=3: Minimum core component size for bridge refinement

Uncertainty Correction (standalone; 0 in the integrated pipeline)

• se_adjust=:auto: :auto runs the finder (estimate_se_adjust) and uses τ̂ (skipped if the SMLD is already σ-corrected); 0.0 = no correction; a numeric value = manual τ (μm) added in quadrature to per-loc σ (σ²+τ²) — scalar, (τx,τy) tuple, length-N vector, or (τx_vec,τy_vec).
• force_se_adjust=false: Apply even if the SMLD is already σ-corrected.

Output

• posterior_pixel_size=0.001: Rao-Blackwellized posterior image pixel size in μm (1 nm; 0.0=disable)
• archive_path=nothing: Mmap chain archive path (nothing=disable)
• progress_file=nothing: Write progress to file (nothing=disable)
• verbose=true: Print progress to stdout

Diagnostics from BaGoL analysis for QC and visualization.

BaGoLResult

Complete result from collapsed BaGoL analysis.

estimate_se_adjust(smld::SMLMData.SMLD; kwargs...) -> NamedTuple

Estimate the per-blink systematic localization error τ (the se_adjust value) directly from data, by distribution comparison to Rayleigh(1) via an adaptive over-merge descent. The "find" companion to apply_se_adjust / run_bagol(; se_adjust).

smld must hold raw-σ localizations (e.g. frame-connected blinks). τ is the correction to find, so an already-σ-corrected input (metadata["sigma_corrected"]==true) is rejected.

Returns

A NamedTuple with τ values in μm:

• tau_hat_um — point estimate (feeds run_bagol(; se_adjust = tau_hat_um))
• ci_lo_um, ci_hi_um — spatial-block-bootstrap 95% CI
• ks_at_hat — KS-to-Rayleigh(1) at the estimate (lower = better fit)
• n_bagol — number of BaGoL E-steps run (the descent length)
• path_um — the descent trajectory as (g, m) pairs (grouping τ, M-step τ)
• grouping — the grouping backend used
• diagnostics — nothing, or (with return_diagnostics=true) the arrays the BaGoLMakieExt plots consume, so the extension never re-runs BaGoL.

Keyword arguments

• g_start_um = 0.008 — over-merged starting τ for the descent (μm)
• stop_tol_um = 2e-4 — self-consistency tolerance |m−g| to stop (μm). Tight enough to avoid a premature stop when g_start sits just above the fixed point.
• max_steps = 12 — descent step cap (:dahl descends slower than the old proxy, so it needs more steps to reach the fixed point from an over-merged start)
• ks_noise = 0.01 — KS tie band for the right-biased M-step
• grid_um = 0.0:0.0001:0.014 — fine M-step τ grid (μm)
• n_iterations = 4000, burn_in = 2000 — per-E-step BaGoL chain length
• n_boot = 200 — bootstrap resamples for the CI
• block_um = 1.0 — spatial-block tile size (μm)
• seed = 1 — bootstrap RNG seed (a local RNG; BaGoL's own E-step stochasticity is only controlled as far as the sampler permits)
• grouping = :dahl — :dahl (BaGoL consensus assignment, default) or :mapn_proxy (legacy nearest-emitter proxy)
• return_diagnostics = false — attach the diagnostic arrays
• bagol_kwargs... — forwarded to run_bagol (e.g. partition_sigma, shape, allocation_model). Reserved (managed by the finder, do not pass): se_adjust, force_se_adjust, posterior_pixel_size, verbose.

apply_se_adjust(smld, se_adjust; force_se_adjust=false) -> BasicSMLD

Return a copy of smld with per-localization σ inflated in quadrature by se_adjust (σ² → σ² + τ², per axis; see run_bagol's se_adjust). No-op (returns smld unchanged) if se_adjust is 0 or the SMLD is already σ-corrected (unless force_se_adjust). Used by run_bagol and by renders to display the localization uncertainty BaGoL actually used.

run_collapsed_chain(locs; kwargs...) -> CollapsedChainResult

Run the collapsed Gibbs sampler on a set of localizations.

Arguments

• locs: Vector of localizations (Emitter2DFit or similar)
• n_iterations=10000: Total MCMC iterations
• burn_in=2000: Burn-in iterations before accumulating
• initial_assignments=nothing: Optional initial assignment vector. If provided, the chain starts from this state instead of all-in-one. Used for diagnostics (transition matrix test, kernel verification).
• shape=2.0: Initial Gamma shape for count distribution
• learn_distribution=true: Control count distribution learning. true=learn both μ and shape, false=fix both, :mu=learn μ only (fix shape), :shape=learn shape only (fix μ)
• learn_rho=true: Learn Poisson-K emitter density ρ independently of μ/shape. For fully fixed hyperparameters, set both learn_distribution=false and learn_rho=false.
• rho=nothing: Initial ρ in emitters per unit area. If provided with learn_rho=false, hold this fixed for Fazel-style fixed-ρ-over-bounding-box flat Poisson-K runs.
• hierarchical_interval=100: Iterations between μ/shape MH updates
• accumulators=AbstractAccumulator[]: List of accumulators to update after burn-in
• verbose=false: Print progress
• callback: Optional callback (iter, state, μ, shape) -> nothing
• callback_interval=1: How often to call callback

Move distribution

• 50%: Allocation Gibbs sweep (full sweep per iteration)
• 25%: Split/merge (K ± 1, restricted Gibbs scan)
• 25%: Birth/death (K ± 1, singleton detach/absorb)
• n_restricted_scans=5: Jain-Neal restricted Gibbs scans per split/merge. 0 = sequential allocation only (Round 3 behavior).

  0 = launch + (n-1) intermediate sweeps + 1 final sweep with density.

• n_bd_substeps=5: Number of birth/death attempts per BD selection. Multiple substeps increase K-transition throughput without changing the target distribution. Each substep is independent MH with proper acceptance.

initialize_from_assignments(assignments, locs) -> CollapsedState

Initialize a CollapsedState from an explicit assignment vector. Used for diagnostics (transition matrix test, kernel verification).

assignments[i] is the 1-based cluster label for localization i. All workspace buffers and precomputed data (LocPrecision) are built fresh.

CollapsedState

MCMC state for the collapsed Gibbs sampler. Stores only assignments (which locs belong to which cluster) — emitter positions are derived from ClusterStats sufficient statistics.

Cluster slots are pre-allocated and reused via the active bitvector. Workspace buffers (_perm, _active_slots, _log_probs, _rollback_*) are pre-allocated for zero-allocation hot-path operation.

CollapsedChainResult

Result from a single collapsed chain (one partition). Lightweight — no sample storage, just final state + accumulator results.

ClusterStats{D,L}

Immutable sufficient statistics for a cluster of localizations in a D-dimensional Gaussian feature (L == D*D).

Stores the posterior precision matrix Λ (D×D symmetric), the natural parameter vector η = Σ_i Λ_i μ_i, the quadratic form Σ_i μ_iᵀ Λ_i μ_i, the localization count, and the sum of per-loc log|Σ_i|. These are all that is needed to compute marginal likelihood, predictive probability, and posterior position/covariance without storing individual localization data.

All operations (add_loc, remove_loc) are O(1) and exact inverses.

estimate_dahl(samples, locs, psm) -> (emitters, posterior_k, stability, assignments)

Dahl's method: select the MCMC sample whose association matrix is closest (squared Frobenius) to the PSM. Equivalent to minimizing posterior expected Binder loss restricted to visited partitions.

Returns emitters with ClusterStats posteriors, the K histogram, per-cluster stability scores, and the chosen sample's assignment vector.

estimate_mapn_collapsed(samples, locs; n_refine=10) -> (Vector{Emitter2DFit}, Vector{Int})

Estimate MAP-N emitters from collapsed chain assignment samples.

Uses iterative Hungarian matching with median-based reference positions to handle label switching. Positions are derived from ClusterStats posterior means (deterministic given assignments).

Algorithm

1. Count K per assignment sample
2. Find MAP-N (most common K)
3. Filter to samples with K = MAP-N
4. For each sample, build ClusterStats per cluster → posterior mean positions
5. Iterative Hungarian matching → update reference to median
6. Final positions: median (robust to label switching)
7. Final σ: ClusterStats posterior covariance (analytic)

Arguments

• samples: Vector of assignment vectors from PartitionSamples
• locs: Original localizations
• n_refine: Number of iterative refinement steps (default 10)

Returns

• emitters: Vector of Emitter2DFit with positions and uncertainties
• posterior_k: Histogram of K values (index k+1 = count of K=k)

estimate_mapn_psm(samples, locs, psm; threshold=0.5, n_refine=10)
    -> (Vector{Emitter2DFit}, Vector{Int})

PSM-corrected MAP-N: determines K from the PSM block structure (threshold + connected components), then uses the standard Hungarian matching pipeline from estimate_mapn_collapsed on samples with that K.

This fixes the linear K bias of histogram-mode MAP-N at large K, while keeping the well-tested Hungarian + median position estimation.

estimate_vi_greedy(samples, locs, psm; n_restarts=3, k_up=0) -> (emitters, posterior_k, stability)

Greedy search under Variation of Information loss (Rastelli & Friel 2018).

Uses cached contingency tables for O(1) per-sample delta computation. Total cost per sweep: O(T × N × K_up) where T = samples, N = locs.

Algorithm

1. Pre-compute T contingency tables between candidate partition and each sample
2. For each loc i (random order), try all clusters + new singleton
3. Compute ΔVI in O(T) by updating 2 contingency entries per sample
4. Accept the move with minimum expected VI
5. Repeat sweeps until no improvement; multiple restarts

Returns emitters with ClusterStats posteriors plus per-cluster stability scores.

estimate_mapn_overlap(samples, locs, dahl_assignments; min_overlap_frac=0.5)
    -> (Vector{Emitter2DFit}, Vector{Int})

Fast MAP-N extraction using the Dahl partition as reference template. Matches each MCMC sample to the Dahl partition via overlap (contingency matrix) Hungarian matching, requiring only ONE O(K³) call per sample instead of 11 iterative position-based rounds.

Algorithm

1. K from Dahl assignments
2. Filter samples to K = K_dahl
3. For each matching sample: overlap Hungarian → per-cluster overlap filter
4. Only include a sample's posterior mean for cluster i if its per-cluster overlap fraction >= minoverlapfrac (rejects label-switching mismatches)
5. Position: mean of well-matched posterior means
6. Uncertainty: Term 1 (Dahl analytic covariance) + Term 2 (allocation variance from well-matched samples only). This is the law of total variance applied to the collapsed posterior, equivalent to the marginal posterior variance an RJMCMC sampler would produce.
7. Fallback: if no samples match, use Dahl assignments directly

Returns

• emitters: Vector of Emitter2DFit with positions and uncertainties
• posterior_k: Histogram of K values (index k+1 = count of K=k)

AbstractAccumulator

Base type for MCMC accumulators. Subtypes must implement:

• update!(acc, state, locs, μ, shape, iter) — add one iteration
• result(acc) — extract final result
• merge!(acc_target, acc_source) — combine results from two partitions

EmitterCountHist

Histogram of emitter count K across iterations. Result: Vector{Int} where index k+1 = count of iterations with K=k.

PosteriorImage

Rao-Blackwellized posterior image accumulator. Each iteration adds a Gaussian blob N(μpostj, Σpostj) per active cluster, rendering only within 4σ of each blob center. This integrates over position uncertainty analytically, producing much smoother images than point-delta histogramming.

NNDistHist

Histogram of nearest-neighbor distances between emitter posterior means, computed per iteration. Useful for detecting clustering patterns. Pre-allocates position buffers for zero-allocation updates.

PartitionSamples

Stores thinned assignment vectors from the MCMC chain. Used for MAP-N estimation: find the posterior mode of K, filter to matching samples, and Hungarian-match to get consistent emitter positions.

PSMAccumulator

Accumulates the Posterior Similarity Matrix (PSM) during the MCMC chain. C[i,j] = fraction of post-burn-in iterations where locs i and j are co-assigned to the same cluster.

Uses cluster-based iteration for efficiency: O(Σ n_k²) per iteration rather than O(N²), since only locs in the same cluster contribute.

Result: (psm=Matrix{Float64}, n_samples=Int)

Partition{E<:SMLMData.AbstractEmitter}

A partition of localizations for parallel processing.

Fields

• id: Unique partition identifier
• locs: Vector of localizations in this partition (typed for dispatch)
• original_indices: Indices mapping back to original input
• is_boundary: Flags for localizations near partition edge
• parent_id: 0 for original DBSCAN clusters, >0 if split from oversized
• is_overlap: Flags for localizations included from sibling sub-partitions for spatial context. Emitters primarily composed of overlap locs are discarded after MAP-N extraction.

partition_locs(locs; partition_sigma, min_size, max_size, skip_size, boundary_margin, overlap)

Partition localizations using precision-weighted DBSCAN.

Two localizations are neighbors if ||p_i - p_j|| / (σ_i + σ_j) < partition_sigma.

Arguments

• locs: Vector of localizations
• partition_sigma=3.0: DBSCAN threshold in sigma units (Inf = no partitioning)
• min_size=0: Minimum locs per partition (clusters below this are noise)
• max_size=1000: Split partitions larger than this
• skip_size=typemax(Int): Skip partitions larger than this (Inf = never skip)
• boundary_margin=0.0: Distance from edge to flag as boundary (0 = auto: partition_sigma×median(σ))
• overlap=:auto: Overlap width for bisected oversized partitions. :auto uses partition_sigma × median(σ) (μm). Numeric value in μm. 0.0 disables overlap.
• bridge_ratio=0.0: Split DBSCAN clusters joined only by low-density bridges. Disabled by default. Values in (0, 1] prune points whose local DBSCAN-neighbor count is below this fraction of the cluster's median neighbor count.
• min_split_size=3: Minimum core component size retained by bridge refinement.

Returns

• partitions: Vector of Partition for valid clusters
• skipped: Vector of Partition for skipped oversized clusters

Base type for spatial prior models (FlatSpatial, LocmixSpatial).

FlatSpatial <: AbstractSpatialModel

Flat uniform spatial prior over a region of area exp(log_area). Each cluster's marginal likelihood includes a -log(A) Occam factor.

LocmixSpatial <: AbstractSpatialModel

Localization mixture spatial prior via grid-based bilinear interpolation.

Cluster marginal likelihood under the spatial model.

Predictive probability for adding a loc to a cluster.

Base type for allocation models (DMAllocation, DecoupledAllocation).

DMAllocation(gamma=nothing)

DM/Polya partition prior: P(z|K,N) ∝ ∏k Γ(nk + γ).

• gamma=nothing: use count-distribution shape as γ (default, preserves NegBin identity)
• gamma=<Float64>: fixed DM concentration, independent of count shape

Decoupled allocation: no partition prior. Spatial ML determines z.

CategoricalAllocation

Labeled categorical allocation with equal-weight categories: each loc independently picks a label uniformly from {1..K}. Yields P(z|K) = K^(-N) regardless of resulting cluster sizes — matches Fazel et al. supplement (after analytically integrating θ_k under uniform prior, the labeled allocation factor is exactly K^(-N)).

At fixed K, behavior matches DecoupledAllocation: predictive-only Gibbs, no DM correction. At K-changing MH (split/merge, birth/death), this allocation contributes Δalloc = -N·log(Knew/K_old), which previous :decoupled was missing (it contributed 0 across K).

Uniform spatial prior over rectangular region.

count_model_map_k(n_locs, mu, shape) -> Int

MAP estimate of K from count model alone (no spatial information). Maximizes P(K | N, μ, shape) ∝ NegBin(N; K×shape, shape/(shape+μ)).

This is the Q-PAINT baseline: best you can do without resolving emitters.

save_posterior_png(filename, post; percentile=0.99, colormap=:inferno)

Save posterior image as a color PNG with percentile scaling.

Scales non-zero pixels so the percentile value maps to 255. Values above the percentile are clipped to 255.

Arguments

• filename: Output path (should end in .png)
• post: NamedTuple from posterior_image()
• percentile: Percentile of non-zero pixels for max scaling (default 0.99)
• colormap: :inferno (default), :hot, or :gray

BaGoLArchive

Writer/reader for binary assignment chain archives.

Write mode: created with BaGoLArchive(path, n_partitions, partition_sizes) Read mode: created with BaGoLArchive(path)

compute_report(result_smld, diagnostics; true_positions=nothing, locs_smld=nothing)

Compute standard analysis report from BaGoL results.

Category 1 (no GT): always computed

Returns metrics for convergence assessment, cluster statistics, and emitter output.

Category 2 (with GT): when true_positions is provided

Adds matching metrics (Jaccard, RMSE, F1), calibration, and oracle comparison.

Returns

NamedTuple with all computed metrics. Pass to write_report() for disk output, plot_report() (CairoMakie extension) for figures, or render_report() (SMLMRender extension) for spatial visualizations.

write_report(report; output_dir="output")

Write standard report files to disk. No plotting dependencies required.

Creates:

• summary.txt — human-readable summary
• emitters.csv — emitter positions and uncertainties
• metrics.json — machine-readable metrics (if GT available)
• posterior_image.png — posterior density image (if available)

match_positions(estimated, true_positions; threshold=0.020) -> NamedTuple

Match estimated emitters to true positions using Hungarian algorithm. Threshold in μm (default 20 nm).

Returns (assignments, matched_distances, cost_matrix) where assignments[i] = matched true index for estimated emitter i (0 if unmatched).

compute_fov(smld; margin=nothing) -> (x_min, x_max, y_min, y_max)

Compute standard field of view from localizations. Bounds cover all localization 1σ circles plus margin (μm) on each side.

If margin is not specified, defaults to 5× median σ.

Use this to set consistent bounds for run_bagol(; posterior_xlim, posterior_ylim) and render_report(; fov=...).

plot_report(report; output_dir="output")

Plot standard report figures. Requires using CairoMakie.

render_report(locs_smld, bagol_smld; output_dir="output", true_positions=...)

Render SMLMRender visualization suite. Requires using SMLMRender.

SimulationResult

Result of SMLM simulation containing localizations, SMLD, and ground truth.

Fields:

• smld: Ready-to-use BasicSMLD with auto-sized camera
• true_positions: Ground truth emitter positions (μm)
• true_counts: Localizations generated per emitter (all pass min_photons when photon sampling is used)
• n_emitters: Number of true emitters

nmer_positions(n, diameter; center=(0.0, 0.0))

Place n emitters on a circle of given diameter (μm). n=1 returns a single point at center.

nmer_grid_positions(; n_per_cluster, cluster_diameter, grid_nx, grid_ny, grid_spacing)

Place n-mer clusters on an grid_nx × grid_ny grid with grid_spacing (μm) between centers. Grid is centered so that the center of mass is at the field center.

Returns (all_positions, cluster_centers).

nmer_random_positions(; n_per_cluster, cluster_diameter, density, field_size)

Place n-mer clusters randomly at given density (clusters/μm²) within a square field. Number of clusters = round(Int, density * field_size²).

Returns (all_positions, cluster_centers).

simulate_localizations(true_positions; kwargs...) -> SimulationResult

Generate synthetic SMLM localizations from known emitter positions.

Count models (count_model)

• :poisson — NegBin with shape=1000 (≈ Poisson)
• :exp — NegBin with shape=1 (geometric/exponential blinking)
• <number> — NegBin with that shape (e.g., 5.0)
• :fixed — exactly round(Int, mean_count) per emitter

Photophysics modes

• Default: photons from Exponential(mean_photons), redrawn until min_photons passes, precision σ = psf_sigma / √photons
• fixed_sigma mode: constant σ and photons, no photon sampling

Returns

SimulationResult with count_params = (μ=empirical_mean, shape=effective_shape). Use these directly in run_bagol(; μ=sim.count_params.μ, shape=sim.count_params.shape).

Keywords

• mean_count=10.0: mean localizations per emitter
• count_model=:poisson: :poisson, :exp, :fixed, or numeric NegBin shape
• psf_sigma=0.130: PSF standard deviation (μm)
• mean_photons=500.0: mean photons per localization
• min_photons=100.0: redraw photon samples below this threshold
• fixed_sigma=nothing: if set, use this constant σ (skips photophysics)
• fixed_photons=1000.0: photon count used with fixed_sigma
• background=10.0: background photons per pixel
• pixel_size=0.100: camera pixel size (μm)
• field_size=nothing: field of view (μm); auto-computed if nothing

simulate_nmer(; n, diameter, center=(0.0, 0.0), kwargs...) -> SimulationResult

Simulate a single n-mer cluster. All other kwargs passed to simulate_localizations.

simulate_nmer_grid(; n_per_cluster, cluster_diameter, grid_nx, grid_ny,
                     grid_spacing, kwargs...) -> SimulationResult

Simulate a grid of identical n-mer clusters. All other kwargs passed to simulate_localizations.

make_camera(field_size, pixel_size) -> IdealCamera

Create a square camera with pixels of pixel_size (μm) covering field_size (μm).

print_simulation_summary(result::SimulationResult)

Print summary statistics of simulation result.

run_optimality_sweep(; kwargs...) -> NamedTuple

Run systematic d/σ sweep across (n, μ) conditions with parallel trials.

Conditions

• n_values=[2, 8]: cluster sizes (dimer, octamer)
• mu_values=[5.0, 10.0]: mean blinks per emitter
• d_sigma_range: separation / localization precision grid

Per condition, per d/σ

Runs n_trials independent simulations with count-only, nohier, and hier methods. Trials run in parallel via Threads.@threads.

Returns

NamedTuple with:

• trials: Vector of per-trial results
• curves: Dict{(n, μ) => aggregated curves}
• scorecard: Vector of pass/fail checks
• config: sweep parameters

run_speed_test(; kwargs...) -> Vector{NamedTuple}

Benchmark BaGoL throughput at different dataset sizes.

Generates grids of 8-mers at varying density to control localization count. Measures wall time for full run_bagol pipeline.

Returns

Vector of (n_locs_target, n_locs_actual, n_partitions, wall_time_s, locs_per_s, times).

write_sweep(sweep; output_dir="output")

Write sweep results to disk: sweep_data.csv and scorecard.txt.

write_speed(speed; output_dir="output")

Write speed test results to speed_data.csv.

plot_sweep(sweep; output_dir="output")

Plot optimality sweep results. Requires using CairoMakie.

plot_speed(speed; output_dir="output")

Plot speed test results. Requires using CairoMakie.

plot_se_adjust(result; output_dir="output", n_boot=300)

Plot the estimate_se_adjust τ-finder diagnostics (4 panels: frozen KS-vs-τ landscape + descent path + CI, empirical CDF − Rayleigh(1) at τ̂ with a spatial-block-bootstrap band, Q-Q vs Rayleigh(1) with an upper-tail inset, and σ-stratified ⟨z²⟩/2). Uses ONLY the finder's frozen-grouping output (no extra BaGoL). Requires estimate_se_adjust(...; return_diagnostics=true) and using CairoMakie. Ported from SMLMClustering diagnose_tau.jl.

MultiClusterStats{T} <: AbstractClusterStats

Block-diagonal sufficient statistics: a tuple of per-feature ClusterStats blocks plus the shared membership count n.

FeatureSet{T} <: AbstractSpatialModel

The per-feature priors (one AbstractSpatialModel per block). Occupies the sampler's spatial slot so the engine's spatial_ml/spatial_pred calls dispatch to the block-diagonal sums below.

run_multicue_chain(positions, values, σ_values, feature_set; kwargs...) -> CollapsedState

Run the collapsed Gibbs sampler grouping localizations JOINTLY on 2D position + a scalar cue (e.g. spectral wavelength). Returns the final CollapsedState; read out per-emitter joint estimates with extract_multicue.

build_multicue_precisions(positions, values, σ_values) -> Vector{MultiLocPrecision}

Build per-localization (position-2D + scalar-cue) contributions from a vector of 2D localizations and a parallel scalar cue (value ± σ) per localization.

extract_multicue(state) -> Vector{NamedTuple}

Per-emitter joint estimate from a multi-cue CollapsedState: (x, y, σ_x, σ_y, value, σ_value, n) — position from block 1, scalar cue from block 2.

gaussian_contribution(val, σ) -> LocPrecision{1,1}

Precompute a 1-D Gaussian feature contribution from an observed value and its standard deviation. Used for scalar cues (spectral wavelength, lifetime, ...).

AbstractTargetDensity

Base type for target distributions over the partition space. Implement log_target(td, z, loc_precs, log_area, μ, shape, ρ) for new variants.

DecoupledTarget <: AbstractTargetDensity

Count model × Poisson K prior × collapsed spatial likelihood.

P(z, K | data) ∝ P_count(N | K, μ, α) × P_K(K | ρ, A) × ∏_k ML_flat_k(z)

The Poisson(ρA) K prior combined with flat spatial -log(A) per cluster makes the target area-independent: the A^K from the K prior cancels the A^{-K} from K clusters' uniform position priors.

DecoupledLocmixTarget <: AbstractTargetDensity

Count model × locmix spatial likelihood, with no explicit K prior and no DM/Polya allocation term.

P(z, K | data) ∝ P_count(N | K, μ, α) × ∏_k ML_locmix_k(z)

This matches allocation_model=:decoupled, spatial_model=:locmix.

DMFlatTarget <: AbstractTargetDensity

Collapsed target implied by:

• a Poisson emitter-count prior K ~ Poisson(ρ A),
• iid per-emitter localization counts n_k ~ NegBin(shape, p),
• a flat spatial prior on emitter positions over area A,
• and positions integrated out analytically.

For a labeled assignment vector z with active cluster sizes n_1, ..., n_K,

P(z, K | data) ∝ P_K(K | ρ, A) × P_count(N | K, μ, shape) ×
                 P_DM(z | N, K, shape) × ∏_k ML_flat_k(z)

This is the collapsed form used to test whether the DM partition factor is the correct consequence of the NegBin generative model.

DMLocmixTarget <: AbstractTargetDensity

Default collapsed locmix target:

P(z, K | data) ∝ P_count(N | K, μ, shape) ×
                 P_DM(z | N, K, shape) × ∏_k ML_locmix_k(z)

No separate Poisson K prior is included.

DirectNegBinFlatTarget <: AbstractTargetDensity

Direct form of the same flat/NegBin target before rewriting the count model into P_count(N | K) × P_DM(z | N, K).

For a labeled assignment vector z with active cluster sizes n_1, ..., n_K,

P(z, K | data) ∝ P_K(K | ρ, A) × (∏_k NegBin(n_k; shape, p)) ×
                 (∏_k n_k! / N!) × ∏_k ML_flat_k(z)

The factor ∏_k n_k! / N! is the probability of a specific labeled assignment vector given the count vector. This target is algebraically identical to DMFlatTarget and is useful as an independent implementation check.

DirectNegBinLocmixTarget <: AbstractTargetDensity

Direct NegBin-count form of DMLocmixTarget, useful as an independent algebraic check of the DM decomposition under locmix spatial likelihoods.

FazelFlatTarget <: AbstractTargetDensity

Fazel-equivalent no-drift collapsed target after analytically integrating θ_k ~ Uniform(R):

π(K, Z | Y, Σ) ∝ P_count(N | K, μ, shape) × P(Z | K) × ∏_k ML_flat(D_k)

with labeled equal-weight allocation P(Z|K) = K^(-N). Spatial = flat (uniform prior). NO K prior (ρA) and no DM/Polya term. Matches:

spatial_model=:flat, allocation_model=:categorical, k_prior=:none.

Count distribution: NegBin(N; K·shape, shape/(shape+μ)). (Fazel's original supplement uses Poisson(Kμ); here we keep NegBin and note that Poisson is recovered in the limit shape→∞ for fixed μ.)

log_target(td, z, loc_precs, log_area, μ, shape, ρ) -> Float64

Evaluate the unnormalized log target density for assignment vector z.

evaluate_target(td, z, locs; μ, shape, ρ) -> Float64

Convenience wrapper that builds LocPrecision automatically.

enumerate_canonical_partitions(N, K_max) -> Vector{Vector{Int}}

Enumerate all canonical partitions of N items into at most Kmax groups. Labels assigned in order of first appearance. Count = Σ S(N,k) for k=1..Kmax.

enumerate_labeled_partitions(N, K_max) -> Vector{Vector{Int}}

All labeled (dense) partitions: z ∈ {1,...,K}^N where every label 1..K appears at least once, for K = 1..K_max. Label permutations ARE distinct.

Count = Σ{K=1}^{Kmax} S(N,K) × K! (ordered Bell numbers truncated at K_max).

canonicalize(z) -> Vector{Int}

Relabel so labels appear in order of first occurrence.

exact_posterior(partitions, locs, td; μ, shape, ρ,
                label_multiplicity=:factorial) -> (probs, log_targets)

Normalized exact posterior over canonical partitions.

label_multiplicity controls how canonical partitions are weighted:

• :factorial multiplies each K-cluster partition by K!
• :none uses a single canonical representative with no multiplicity factor

run_kernel_invariance_test(locs, td; μ, shape, K_max=3,
                            n_steps_per_state=1500,
                            label_multiplicity=:factorial,
                            seed=nothing, verbose=false) -> NamedTuple

Build empirical one-step transition matrix P and verify:

1. Invariance: π^T P ≈ π^T (Wald test, not arbitrary threshold)
2. Irreducibility: support graph is strongly connected (empirical evidence)
3. MCMC vs exact: K-marginal z-tests with ESS correction

Also reports overflow fraction and within-K conditional TV.

Uses initial_assignments + learn_distribution=false to hold (μ, shape) fixed. The transition matrix is for the z-marginal kernel only.

run_labeled_invariance_test(locs, td; μ, shape, ρ, K_max,
                             n_steps_per_state, seed, verbose) -> NamedTuple

Kernel invariance test on the FULL labeled state space (no canonicalization). Each of the K! relabelings of a canonical partition is a distinct state. The target density is evaluated directly — no K! multiplicity correction.

Tests:

1. Invariance: πL^T P ≈ πL^T (Wald test)
2. Irreducibility on labeled space
3. Label symmetry: π^T P gives same value for all relabelings of each partition

DetailedBalanceResult

Result of checking detailed balance for a single (zfrom, zto) transition.

check_detailed_balance(z_from, z_to, locs, td;
                       μ, shape, ρ, log_q_fwd, log_q_rev,
                       log_ratio_code=nothing, tol=1e-8) -> DetailedBalanceResult

Verify detailed balance for a specific transition pair.

Computes: Δ = [log π(z') + log q(z'→z)] - [log π(z) + log q(z→z')]

If log_ratio_code is provided, checks |Δ - Δ_code| < tol.

check_detailed_balance(z_from, z_to, loc_precs, log_area, td; ...) -> DetailedBalanceResult

Precomputed-arguments version for checking multiple transitions efficiently.

ChainDiagnosticAccumulator <: AbstractAccumulator

Records K trace during the MCMC chain for post-hoc mixing diagnostics.

effective_sample_size(trace; method=:initial_sequence) -> Float64

Geyer (1992) initial positive sequence ESS estimator.

autocorrelation(trace; max_lag=nothing) -> Vector{Float64}

Normalized ACF. Denominator n·σ² (not (n-lag)·σ²) for PSD guarantee.

split_gelman_rubin(chains) -> Float64

Split each chain in half, compute R-hat on 2m half-chains.

indicator_ess(trace; k_values=nothing) -> Dict{Int, Float64}

ESS on 1[K=k] for MAP-K ±2 neighbors. Catches P(K=k) estimation issues that K-trace ESS misses.

run_mixing_test(locs; n_chains=4, n_iterations=10000, burn_in=2000,
                shape=2.0, seed=nothing, kwargs...) -> NamedTuple

Run multiple independent chains and assess mixing quality. Returns ESS, split R-hat, indicator ESS, chain means/vars, acceptance rates.