Internals

sweep!(perm, lg_table, model_dof, log_two_outlier_halfwidth, dg, scores, n; ...) -> (S*, k*)

Low-level sort-and-sweep kernel for the scale-free marginal score (Algorithm 1).

The per-hypothesis path — deliberately not the KernelAbstractions kernel

The batched marginal GRID is one KernelAbstractions kernel (ig_sweep_kernel! under grid_evaluate(::MarginalScoring, …), CPU and GPU alike). This function is kept beside it on purpose, as the ONE per-hypothesis path, and is not a twin of it:

• It is latency-bound. A RANSAC loop scores one hypothesis at a time, and routing one through the kernel on KernelAbstractions.CPU() (squares sorted into a preallocated flat buffer, then a 1×1×1 launch) measured slower (2026-09-21, BenchmarkTools, Float64): n = 200: 1.88 → 2.57 µs (1.37×); n = 715: 6.16 → 10.61 µs (1.72×); n = 2500: 24.5 → 31.7 µs (1.29×).
• It does more. The kernel scores pre-sorted, penalty-free Jeffreys/IG columns; this also sorts (the radix below), carries the predictive per-point penalties, the profiled-γ prior_table, and the exact parametric sweep for heterogeneous gains.

Sorts points (ascending by scores, or sort_key if provided), then sweeps prefixes k = m+1,...,n evaluating the partition score via _sweep_score(scale_prior, ...). Returns the best score S* and optimal inlier count k*. Mutates perm.

Complexity — depends on which sort path is taken (see Sizing)

• Θ(n) on the exact-size fast path: _radix_sortperm! is an LSD radix over monotone unsigned keys — one fused histogram scan plus at most sizeof(K) stable counting iterations (4 for the default UInt32), with trivial digits skipped. No comparisons, so no log factor. The sweep itself is a single Θ(n) scan, so the whole kernel is linear.
• Θ(n log n) on the partial-buffer path, which falls back to the comparison-based Base.sortperm!.

This is not a constant-factor distinction: an oversized workspace silently changes the complexity class. Measured on PT-E-SIFT-4K (2026-07-26), the same scorer cost 51.5 µs/hypothesis through the partial path and 21.7 through the radix path, because a caller had sized one workspace to the pool-wide maximum instead of per-problem. RobustVisualGeometry's _set_n! resizes every sortable buffer to n for exactly this reason.

Sizing

Two invariants are handled:

• Exact-size fast path (length(perm) == n and length(sort_target) == n): sorts via _radix_sortperm! when radix buffers of capacity ≥ n are present (the emitted permutation is bitwise identical to sortperm!), else via sortperm! directly. Hot path for RANSAC scoring, where ws.scores, sws.perm, sws.scratch are all sized to data_size(problem).
• Partial-buffer path (length(perm) > n or length(sort_target) > n): sorts views of the first-n prefix through sortperm!. The radix is NOT reachable here — it requires exact-size concrete Vectors, and a @view is neither. Used by BlobBoards blob adaptation (one sws sized to max_n, per-blob n varies) and the CUDA-ext host-side seed loop.

scratch (if provided) is forwarded to sortperm! only when length(scratch) == n, because Base.sortperm! rejects SubArray scratch.

sweep!(sws::ScoringWorkspace, sq_residuals, n;
       model_dof, log_two_outlier_halfwidth, codimension=1) -> (S*, k*)

The sort-and-sweep on a bare ScoringWorkspace, for callers without a MarginalScoring. Student-t / posterior IRLS callers have a ScoringWorkspace and raw spec arguments (model_dof, log_two_outlier_halfwidth, codimension).

The scale prior is not a parameter here: it is a property of sws, committed at allocation time (see ScoringWorkspace). Callers that want IG choose the prior when building the sws via default_scoring_workspace(n; dg, scale_prior=InverseGammaPrior(...)).

Used by BlobBoards blob adaptation (posterior_adapt.jl) and the CUDA extension's host-side seeding loop. Both of these share one sws across many calls with varying n, so the partial-buffer path in the low-level kernel is the common case.

Forwarder — no sweep loop of its own.

default_scoring_workspace(n; dg=1, scale_prior=JeffreysPrior()) -> ScoringWorkspace

Vanilla scoring workspace, flat γ prior. Length-n buffers for a problem of n points with codimension dg. The scale_prior kwarg determines the lg_table's α₀; the returned sws carries that prior and cannot be repurposed for a different one.

Used by standalone callers (BlobBoards blob adaptation, the CUDA extension's host-side seed loop, raw low-level RansacWorkspace construction in tests) where no MarginalScoring is in hand.

alloc_scoring_workspace(problem, scoring::MarginalScoring) -> ScoringWorkspace

Allocate scratch buffers for Algorithm 1 sweep. Depends on:

• data_size(problem) — buffer lengths
• codimension(problem) — lg_table step size
• scoring.scale_prior — lg_table α₀ offset (Jeffreys → 0.0)
• scoring.profile_prior — whether the profiled-γ correction table is built

Called by RansacWorkspace(problem, scoring) so buffers live on RansacWorkspace.scoring_ws.

alloc_scoring_workspace(n::Int, dg::Int, scoring::MarginalScoring) -> ScoringWorkspace

Size-only overload: build a ScoringWorkspace from raw n / dg instead of a live problem instance. Useful for benchmarks and pre-allocation contexts where the sweep buffers must be sized for the maximum dataset a thread will see, before any particular problem is in hand.

The resulting sws carries scoring.scale_prior, so subsequent sweep!(sws, ...) calls use that prior without any kwarg plumbing.

init_score(scoring::AbstractScoring)

Return the initial "best" score value for the scoring strategy. MarginalScoring: (-Inf, -Inf) (global, local score tuple to maximize)

_build_lg_table(n, d_g=1; alpha0=0.0) -> Vector{Float64}

Precompute log Γ(k·dg/2 + α₀) for k = 1..n via SpecialFunctions.loggamma — single uniform path for Jeffreys (α₀=0) and Inverse-Gamma (α₀>0) priors.

Cost: ~5000 host-side loggamma calls per probe (~0.5 ms total) — table is read-only on the device for the entire sweep. An earlier version used a log-recurrence specialised on Jeffreys; we standardised on loggamma so there's a single code path with uniform Float64-precision values, no recurrence drift, and no separate code branch to maintain.

outlier_halfwidth(scoring::MarginalScoring) -> Float64

The half-width of the uniform outlier support per axis. The sweep charges log(2·outlier_halfwidth) per point and axis.

scale_prior(sws::ScoringWorkspace) -> IGPrior
scale_prior(scoring::MarginalScoring) -> IGPrior

The scale prior the workspace is committed to (and that sweep!(sws, …) and posterior_weights! score under), or the one a MarginalScoring carries.

sweep_order(sws::ScoringWorkspace) -> Vector{Int}

The order the last sweep! visited the points in: sweep_order(sws)[1:k] is the support of the partition it returned k for (the SeedPartition/mask contract).

penalties(sws::ScoringWorkspace) -> Vector{Float64}

Per-point log-determinant penalties ℓᵢ = log|Vᵢ| of predictive scoring (MarginalScoring{Predictive}), filled by the scoring glue before a sweep.

sort_key(sws::ScoringWorkspace) -> Vector{Float64}

The per-point sweep ordering key of predictive scoring (see ScoringWorkspace), filled by the scoring glue before a sweep.

lg_table(sws::ScoringWorkspace) -> Vector{Float64}

The log Γ(k·dg/2 + α₀) table the workspace was built with; its length is the largest n the workspace can sweep (see resize!).

_build_score_weight_matrix(scoring, loss_widths, N_bins, d2_min, d2_max)
    -> Matrix{Float32}  size (N_bins, T)

Precompute W[bin, t] = residual_quality(scoring, d²_bin, loss_width_t) on the host at the bin centres of a log-spaced grid. The centres are d²_bin = d²_min · (d²_max/d²_min)^((bin-0.5)/N_bins), and W holds the SAME scalar quality the Exact kernel evaluates per residual. Returned column-major so the kernel can index W[bin + (t-1)·N_bins].

_predictive_score_penalty(rᵢ::T, ∂θgᵢ_w::SVector{n_θ,T}, Σ̃_unit) -> (qᵢ, ℓᵢ)

Model-uncertain score and penalty for a scalar constraint (dg=1). It uses the MULTIPLICATIVE predictive-covariance convention of the generalized sweep lemma (uncertain_evidence 35ae25c): the full predictive variance is Ṽᵢ = s²·Vᵢ with the UNIT-SCALE shape

Vᵢ = 1 + (∂θgᵢw)ᵀ Σ̃unit (∂θgᵢ_w),

where Σ̃unit is the propagation of UNIT measurement noise through the solver/fit (the physical Σθ = s²·Σ̃_unit — the hypothesis derives from the same noise process the measurements carry). Returned qᵢ = rᵢ²/Vᵢ and ℓᵢ = log Vᵢ are therefore s-FREE; the s-dependence factors into the common scale marginal exactly as the lemma's separation requires. (The former s2 + argument slot read as an ADDITIVE convention with the model covariance frozen at a reference scale — a latent modeling inconsistency while s is profiled; numerically identical at s2 = 1, semantically corrected here.)

_predictive_score_penalty(rᵢ::StaticVector{dg,T}, ∂θgᵢ_w::StaticMatrix{dg,n_θ,T}, Σ̃_unit)
    -> (qᵢ, ℓᵢ)

Vector-constraint (dg≥2) form of the MULTIPLICATIVE convention (see the scalar method). Here Vᵢ = I + (∂θgᵢw) Σ̃unit (∂θgᵢ_w)ᵀ (unit-scale dg×dg shape), qᵢ = rᵢᵀ Vᵢ⁻¹ rᵢ and ℓᵢ = log|Vᵢ| — both s-free; Ṽᵢ = s²·Vᵢ.

_score_improved(scoring, new, old) -> Bool

Check whether new score is strictly better than old. All strategies use higher = better convention.

_sq_norm(r) -> T

Squared norm: rᵢ² for scalar, rᵢᵀrᵢ for vector. Computes ‖rᵢ‖² from the shape-whitened residual (Table 2).

_wrap_solver_cov(scoring, problem, model, sample_indices)

Wrap model with estimation covariance Σ̃_θ from the minimal solver.

• MarginalScoring{Nothing}: returns plain model (model-certain, no wrapping)
• MarginalScoring{Predictive}: returns Uncertain(model, CovMat(Σ̃_θ)). A problem with no solver_covariance throws an ArgumentError. A provider's nothing (a DEGENERATE sample) returns nothing: no hypothesis, which the caller skips like a failed solve. It returned the plain model until 2026-09-22, scoring a degenerate sample's model certain.

RadixPermBuffers{K}(n)
RadixPermBuffers(n)          # K = UInt32 (Float32 key precision) — THE default

Scratch buffers for _radix_sortperm! — the allocation-free stable LSD radix sortperm used by the sweep! single-sort fast path. Owned by ScoringWorkspace (sws.radix); sized to the data count and trimmed alongside the other sweep buffers (RVG set_n!).

Key width K — the sort-precision knob (2026-07-26)

The radix does sizeof(K) byte iterations, each scattering a K key plus an Int perm entry, so the key width sets both the iteration count and the scatter bandwidth. K = UInt32 (sort keys at Float32 precision) is the default: measured 1.65× faster than UInt64 on real PhotoTourism residual columns (10.12 → 6.15 µs/hypothesis; ≈1.16× on the whole marginal per-pair scorer, then score_marginal), and it loses nothing, because the scored residuals are Float32 to begin with — a pool's residual matrix is a Matrix{Float32}, so the Float64 sweep buffer holds an exactly-representable widening and the round trip back to Float32 is lossless. Verified on that pool: the emitted permutation is BYTE-IDENTICAL to the UInt64 path in 150/150 hypotheses (worst residual delta exactly 0.0; experiments/f32_radix_probe.jl).

K = UInt64 remains available and is the right choice for a caller whose residuals carry genuine Float64 information (they are then ordered at full precision). Both widths run the SAME generic code below — the width is a type parameter, not a second implementation.

Ordering semantics at K = UInt32: the key is monotone but not injective on Float64 (values within one Float32 ulp collapse to equal keys). Stability then orders such a tie group by index, exactly as a stable comparison sort under an isless that cannot separate them. Every sweep prefix boundary inside such a group is between residuals equal to Float32 precision.

_radix_sortperm!(perm, target, rb::RadixPermBuffers, n) -> nothing

Fill perm[1:n] with the stable ascending-isless permutation of target[1:n] — bitwise identical to sortperm!(perm, target) for every input (see the stability argument in the header comment). O(n): one fused histogram scan, then ≤ 8 stable counting iterations over byte digits, skipping digits on which all keys agree (a stable no-op). Allocation-free; buffers may have capacity ≥ n.

_radix_key(x::Float64, ::Type{K}) -> K

Monotone key: _radix_key(a, K) < _radix_key(b, K) ⟹ isless(a, b). For K = UInt64 the converse holds too — _radix_key(a) == _radix_key(b) ⟺ neither isless(a, b) nor isless(b, a) (every NaN maps to the single maximum key). Total-order equivalent of isless on all Float64 values, including ±0.0 (distinct, -0.0 < +0.0 — exactly isless's ordering) and ±Inf.

For K = UInt32 the value is first rounded to Float32 (the precision the scored residuals actually carry); the map stays monotone non-decreasing, so the emitted order is a valid ascending order, with ties broken by index through LSD stability.

_build_prior_table(n)

Precompute profiled-γ correction for k = 1..N:

prior_table[k] = k·log(k/N) + (N−k)·log((N−k)/N)

This is the log of the profiled label prior γ̂^{nᵢ}(1−γ̂)^{nₒ} with γ̂ = k/N (MLE of the inlier ratio for a given partition size k). Equivalent to −N·H(k/N) where H is binary entropy.

Convention: 0·log(0) = 0, so prior_table[N] = 0.

The Student-t mixture of a hard sweep! partition: its posterior predictive, partition_student_t — the SAME (ν, ŝ²) rule a PosteriorWeightState seeds from.

The RSS is summed in Float64 REGARDLESS of the residuals' storage eltype (Float32 storage upstream is fine, Float32 accumulation is not). prior adds its pseudo-observations — 2β₀ to the RSS, 2α₀ to the dof — so the mixture is the posterior predictive (see posterior_weights!); at Jeffreys both are + 0.0. The inlier prior is ½, the fixed prior of the sweep.

Until 2026-09-30 this restated the rule with two differences: it floored the RSS (max(RSS + 2β₀, eps)/ν) where fit_scale floors the quotient, and it did not clamp ν ≥ 2 — so at a thin partition (ν_data + 2α₀ < 2) the mixture RobustVisualGeometry's _finalize weighs with was a different Student-t from the one its LO seeded with.

resize!(sws::ScoringWorkspace, n) -> sws

Size every per-point buffer of sws for a problem of n points. They are the sort permutation and its scratch, the predictive penalties (zeroed, as a fresh workspace has them) and sort_key, and the radix perm-sort scratch — keeping the exact-size radix fast path of sweep! reachable. A profiled-γ prior_table is REBUILT for n, since k·log(k/N) + (N−k)·log((N−k)/N) depends on N: truncating a table built for a larger N would keep that N's values. At n = 0 the table is left as it is, because its non-emptiness is the workspace's only record that the scoring is profiled.

lg_table is not resized: its entries depend on k alone, so a table of capacity ≥ n serves any n. It cannot grow — the workspace does not record the codimension it was built for — so n > length(lg_table(sws)) throws; build a larger workspace instead. (Growing the other buffers past it used to leave sweep! reading beyond the table under @inbounds.)

Replaces RobustVisualGeometry's _set_n!(::ScoringWorkspace, n), which reached into seven of these fields by name and truncated the prior_table.

PrecomputedPair

Sorted squared residuals per model of one pool pair — the staging unit of stage_ig_shared_from_pairs: sorting once here makes every device sweep an O(n) linear scan over prefix sums.

ig_sweep_kernel!(scores_out, sorted_sq_flat, pair_offset, K_per_pair, M_per_pair,
                 lg_table, log_two_outlier_halfwidth_grid, model_dof, dg, prior, p_offset)

One work-item per (pair, model, grid-point): ONE sequential sweep over the model's K sorted squared residuals, writing the best partition score max_k S(k) to scores_out[p_local, m, t] (Float64). S(k) is sweep.jl's _partition_score — the one per-partition formula every sweep shares, host and device — at the grid point's log(2·outlier_halfwidth). RSS accumulates in Float64. Work-items beyond a pair's model count write nothing (the chunked driver pre-fills −Inf).

ScaleMarginal{K}(alpha0, beta0; outlier_halfwidth, t0 = nothing)

The scale-marginal spec of loss shape K at the IG prior (alpha0, beta0) on the squared loss width — e.g. one fit by type-II MLE (empirical Bayes), where beta0 is not tied to a prior centre. outlier_halfwidth bounds the density normalizer Z̃ (REQUIRED: it is the declared outlier support, POOL_FRAME_R_MAX in the pool frame's unit); t0 the centre of the log-scale domain, by default the prior mode log(beta0/(alpha0 + 1)) (or log((outlier_halfwidth/4)²) for an improper prior). THE one constructor of every family's spec: the keyword forms below only name the prior centre.

SmoothTLSMarginal(; alpha0=1.0, prior_loss_width, outlier_halfwidth)

Marginal smooth-TLS scoring spec. prior_loss_width is the IG-prior centre on the smooth-TLS loss width (beta0 = alpha0·prior_loss_width² ⇒ 2·alpha0 pseudo-observations at pseudo-squared-width prior_loss_width²); outlier_halfwidth bounds the density normalizer Z̃ (set to POOL_FRAME_R_MAX, the declared outlier support in the pool frame's unit — required, no default). alpha0/beta0 match IGPrior, as for WelschMarginal. The evidence is the trapezoid scale quadrature of loss_spec(::ScaleMarginal).

family_tag(family::ScalePriorFamily) -> Symbol

Compact symbol identifying the prior family. :jeffreys (parameterless) or :ig (inverse-gamma). Used in display + sidecar-key construction.

fitting_tag(fitting::PriorFitting) -> Symbol

Compact symbol identifying the fitting algorithm. :none for parameterless families, :robust for median-based, :mle for closed-form MLE.

tabulate!(tab::TabulatedLift, m::StudentTMixture; warm = false) -> w_max

Tabulate VisualGeometryCore's half-quadratic conjugate κ² of the mixture's IRLS weight — which has no closed form — from the weight and its slope, posterior_weight and posterior_weight_derivative (the pair a Newton inversion of the weight reads). warm starts from the previous table.