Pools, grids and tuning

best_pool_model(scoring::LossWidthScoring, pair::PoolPair, R, loss_width::Real) -> Int
best_pool_model(scoring::MarginalScoring, pair::PoolPair, R; sws, q) -> Int
best_pool_model([problem_type::Type,] scorer, pair::PoolPair, R; prior) -> Int

Index of the best-scoring stored pool model of pair, scored on its residual matrix R. R is num_points(pair) × num_models(pair), Float32, in NormalizedImageFrame units — every width here is in that unit too; 0 means no model (see the sentinel contract above). A matrix of another shape is a DimensionMismatch.

• A loss-width family scores each model by Σᵢ residual_quality(scoring, dᵢ², loss_width) at the loss_width (for MAGSAC++ its cutoff), exactly the Exact grid's (grid_evaluate) per-model score on this one pair.
• The marginal takes the max over partitions (the sort-and-sweep). sws::ScoringWorkspace and q::Vector{Float64} are scratch keywords, allocated for pair when omitted: supply per-thread buffers sized for the largest pair to reuse them.
• The typed pool scorers carry their hyperparameter: MarginalPoolScorer(outlier_halfwidth), MsacPoolScorer(loss_width), RansacPoolScorer(loss_width). The marginal one needs the problem_type (its minimal sample size and codimension) and takes the scale prior; the loss-width ones accept both and ignore them, so a caller can supply them uniformly.

pool_model_score(scoring::LossWidthScoring, pair::PoolPair, R, j, loss_width[, q0]) -> Float64

Score of the stored pool model j of pair under a loss-width family (MsacScoring, RansacScoring, MagsacScoring, the loss scorings) at loss_width. It is Σᵢ residual_quality(scoring, dᵢ², loss_width) over the model's column of the pair's residual matrix R (num_points(pair) × num_models(pair), Float32, in NormalizedImageFrame units — the pool frame's, as is loss_width). The ONE per-model loss-width score: best_pool_model is its argmax, the Exact grid (grid_evaluate with Exact()) computes the same value, and a running-best selector over a pool calls it per (pair, j) rather than re-summing the quality.

The stack's minimal-sample rule, the same in every scorer: the model's minimal-sample points (pair.sample_indices[:, j]) are fit exactly by construction and carry no information, so they enter at ZERO residual whatever was stored for them — here each contributes the constant q0 = residual_quality(scoring, 0, loss_width), identical for every model of the pool, so the ranking is that of the non-sample points alone. (The marginal scorers zero the same entries and subtract the model's m dof in the partition evidence, ν = dg·(k − m).) A non-finite stored d is an outlier (quality 0); finiteness is tested before squaring, so a finite d whose square overflows also scores 0. An invalid sample index throws. q0 is an optional precomputed cache.

Stu-tU/IG marginal pool scorer at an outlier half-width.

MSAC pool scorer at a loss width: Σ max(0, 1 − d²/loss_width²).

RANSAC pool scorer at a loss width: the inlier count Σ 1(d² < loss_width²).

precompute_pair(pair::PoolPair, R) -> PrecomputedPair

Per-pair staging precompute of pair and its residual matrix R (num_points × num_models, Float32, NormalizedImageFrame units). Each model column is squared in Float32, its minimal-sample points zeroed (exact inliers; an out-of-range index THROWS, as in every other scorer) and sorted; the selection error pool_imc is kept per model.

IGSharedData{Backend, F32A, F64A, IntA}

Backend-agnostic flat-layout view of every PrecomputedPair in a pool. Built once per pool by stage_ig_shared_from_pairs and reused by every grid evaluator.

Memory layout — all distinct pairs flattened:

• sorted_sq_flat::F32A (a vector of the pool's residual eltype on backend) contiguous concatenation of all (Ki × Mi) sortedsq matrices in column-major order: pair 1 model 1 (length K1), pair 1 model 2, …, pair 1 model M_1, pair 2 model 1, …
• pair_offset::IntA (length P+1) pair_offset[p] = absolute starting index of pair p in sorted_sq_flat. pair_offset[P+1] = total length (sentinel).
• K_per_pair::IntA (length P) — # residuals for each pair
• M_per_pair::IntA (length P) — # models for each pair
• M_max::Int — max # of models across all pairs (kernel grid sizing)
• errors_imc_flat::F64A — concatenated per-model selection imc = pool_imc (H: Δrot; E/F: max(Δrot,Δtrans)), indexed by errors_offset[p] + j
• errors_offset::IntA — errors_offset[p] = start of pair p's imc segment

P is the number of pairs in the pool.

stage_ig_shared_from_pairs(pairs, residuals; backend) -> IGSharedData

Stage a pool once for every grid evaluator. The inputs are the light records pairs (AbstractVector{<:PoolPair}) with their residual matrices residuals (residuals[p] is pair p's num_points × num_models Float32 matrix, in NormalizedImageFrame units — the caller loads them; staging opens no file), each pair precomputed by precompute_pair into the flat IGSharedData layout on backend. No prior fit, no width grid baked in.

grid_evaluate(scoring::MarginalScoring, staged::IGSharedData, outlier_halfwidths;
              pair_chunk_size = 256, return_idxs = false) -> Matrix{Float64}

The Student-t marginal grid (the IG/Jeffreys sort-and-sweep) on the staged pool's backend. It returns a (P, T) matrix of each pair's error at the model scoring selects at each outlier half-width in outlier_halfwidths (NormalizedImageFrame units). The grid REPLACES the scoring's own outlier half-width; its minimal sample size, codimension and scale prior are the scoring's. With return_idxs, also the (P, T) Int32 selected model indices.

Runs ig_sweep_kernel! on the shared chunked driver (_chunked_argmax_eval): pair-chunks of pair_chunk_size bound the score block, and the model argmax is the shared argmax_models_kernel!, which takes the FIRST strictly best model and SKIPS a NaN score. The device sweep scores the model-certain marginal with the flat label prior only: a PredictiveMarginalScoring or a scoring with profile_prior = true is an ArgumentError, never scored as something else.

grid_evaluate([S = Float32,] scoring::LossWidthScoring, staged::IGSharedData, loss_widths;
              aggregation = LogBinned(), return_idxs = false, pair_chunk_size = 256)
    -> Matrix{Float64}

Per-pair selection errors of a loss-width family over a loss-width grid. They form a (P, T) matrix of the error (imc) of the model scoring selects at each loss width in loss_widths (NormalizedImageFrame units; with return_idxs, also the selected pool indices). Runs on the staged data's backend — KernelAbstractions CPU() or a GPU.

aggregation chooses the algorithm: LogBinned (default, 200 bins per decade; the fast cross-validation path) or Exact (every residual; the reference). Both evaluate the family's one per-residual quality. S is the LogBinned kernel's bin-index precision (Float32, the default, is the production setting); Exact evaluates in Float64. Scores accumulate in Float64 either way.

scoring is a LossWidthScoring member — build it with loss_width_scoring, supplying the problem's codimension.

grid_evaluate([S = Float32,] m::ScaleMarginal, staged::IGSharedData, prior_loss_widths;
              NG = _MARG_NG, return_idxs = false, mem_budget = 6_000_000_000)
    -> Matrix{Float64}

Full-data scale-marginal scoring (WelschMarginal, SmoothTLSMarginal) at precision S (Float32 by default; Float64 for full precision), on the staged data's backend — KernelAbstractions CPU() or a GPU. The scale is integrated out by the package's one scale quadrature, the trapezoid logsumexp over NG log-scale nodes spanning t0 ± 7.0 around each grid point's prior centre (loss_spec(::ScaleMarginal)). Returns a (P, T) per-pair error matrix at the model each grid-point's marginal selects.

prior_loss_widths is THIS FAMILY'S tuned parameter — the IG-prior centre on the loss width (prior_loss_width, NormalizedImageFrame units), not a decision threshold. Grid point t is m's family with m's prior shape alpha0, Z̃ domain m.outlier_halfwidth (the declared outlier support, POOL_FRAME_R_MAX on pool data), centred at prior_loss_widths[t]: the grid REPLACES m's own centre. return_idxs=true also returns the selected pool indices. The pair-chunk is sized so the NG-per-work-item scratch and the score block stay under mem_budget bytes.

grid_evaluate([S = Float32,] spec::LossSpec, staged::IGSharedData, outlier_halfwidths;
              model_dof, dg, return_idxs = false, mem_budget = 6_000_000_000)
    -> Matrix{Float64}

Track-B mixture-marginal grid evaluation, on the staged data's backend (KernelAbstractions CPU() or a GPU), the ONLY evaluation of this evidence. For the depth-Δ quadrature selector spec (built by loss_spec — Z̃ over the FIXED [0, r_max] domain), score every staged (pair, model) at every outlier half-width in outlier_halfwidths (NormalizedImageFrame units) and return the (P, T) per-pair error matrix at each grid-point's selected model (return_idxs=true also returns the selected pool indices). Three kernels per pair-chunk: marg_mix_evidence_kernel! (E over k; the NG-node trapezoid quadrature), marg_mix_profile_kernel! (the outlier-half-width profile re-argmax over k), and the shared argmax_models_kernel!. Precision S drives the evidence arithmetic — Float64 for full precision; Float32 is the fast production setting. The pair-chunk is sized so the E block + the NG-per-work-item Dscratch stay under mem_budget bytes. model_dof (the minimal sample size, whose points are the exact inliers at the front of each column) and dg (the codimension) are the problem's, and REQUIRED: a LossSpec carries neither.

Aggregation

How a loss-width family's per-residual quality is summed over a model's residuals: Exact or LogBinned. Both evaluate the SAME scalar quality; only the aggregation differs.

Exact()

Every residual through the family's scalar quality — the reference. The quality is evaluated in Float64 whatever the kernel precision S of grid_evaluate.

LogBinned(bins_per_decade = 200)

The fast cross-validation aggregation: residuals histogrammed on a log-spaced d² grid. The grid has bins_per_decade bins per decade, spanning 0.01·minlosswidth² to 100·maxlosswidth², and the family's quality is evaluated at the bin centres. Out-of-range and non-finite residuals score as outliers.

native_axis(::Type{<:AbstractScoring}) -> WidthAxis
native_axis(scoring) -> WidthAxis

The width axis a scoring's pool grid is native on — what its grid's columns are, and therefore which to_native method maps a scale multiple onto it:

native_map(axis::Union{LossWidth, OutlierHalfwidth}; d_g) -> (scale_multiple, scale) -> Float64

Curried to_native on axis at codimension d_g — THE per-axis map closure every window and coverage grid is cut with. d_g is required for both axes (the loss-width axis ignores it), so a grid and the window cut inside it are always cut at the same codimension.

to_native(::OutlierHalfwidth, scale_multiple, scale; d_g) -> Float64

A scale multiple as the Student-t marginal's native width. That width is the outlier half-width whose MAP inlier/outlier decision flips at the residual scale_multiple · scale — decision_width of the Gaussian-uniform mixture GauInlier(scale²; d_g), the package-wide definition of the crossover. scale is the estimated inlier scale (linear, in the units of the residuals; the pool frame's on pool data), and the result is in those units.

d_g is the residual codimension (1 for E and F, 2 for H) and has no default: it is the EXPONENT of the outlier density, so it moves the outlier half-width by decades (at scale = 1, a scale multiple of 6 maps to 8.2293e7 at d_g = 1 and to 1.0156e4 at d_g = 2). VisualGeometryCore owns the generic and the LossWidth method.

outlier_halfwidth_to_loss_width(outlier_halfwidth; operating_scale_multiple, d_g) -> Float64

Convert a marginal's outlier half-width to the truncation loss width at the SAME effective decision boundary. It inverts the Gaussian-uniform crossover at the operating scale multiple: ŝ_implied = outlier_halfwidth / _unit_halfwidth(operating_scale_multiple, d_g), then loss_width = operating_scale_multiple·ŝ_implied.

d_g, the codimension, is REQUIRED, as it is for GauInlier: the crossover moves with it, so no value is a safe default (1 for F/E, 2 for H).

operating_scale_multiple (σ-units) is caller policy — LocalOptimizationDoneRight pins it at OPERATING_SCALE_MULTIPLE (its refit layer's consumed_native). It is built on decision_width, the ONE definition of the crossover. This is NOT a wrapper over it — it answers a different question: not "where is the boundary" but "what is this outlier half-width worth in the other family's units, holding the operating scale multiple fixed". That is a PROTOCOL choice (compare methods, not parameterisations), which is why it takes operating_scale_multiple and the boundary functions do not.

loss_width_to_outlier_halfwidth(loss_width; operating_scale_multiple, d_g) -> Float64

Inverse of outlier_halfwidth_to_loss_width: the marginal's outlier half-width whose decision boundary matches a truncation at loss_width, at the operating scale multiple. Written as the literal reciprocal of outlier_halfwidth_to_loss_width around the same _unit_halfwidth(operating_scale_multiple, d_g), so the two are inverses by construction rather than by algebra. d_g is required.

WelschMarginal(; alpha0=1.0, prior_loss_width, outlier_halfwidth)

Marginal-Welsch scoring spec. prior_loss_width is the IG-prior centre on the Welsch 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̃ — REQUIRED, as for every scale-marginal spec. alpha0/beta0 are the same two IG hyperparameters IGPrior and InverseGammaMarginal name. For an (alpha0, beta0) prior not tied to a centre, ScaleMarginal{WelschShape}(alpha0, beta0; outlier_halfwidth).

The evidence is the trapezoid scale quadrature of loss_spec(::ScaleMarginal).

LossSpec{K<:AbstractLossShape}

The depth-Δ quadrature-marginal spec of the redescending loss shape K: the shape, its well depth Δ, the quadrature grid with Z̃ baked at that depth over [0, r_max], and r_max.

Build it with loss_spec, which keeps Δ and the grid in step, and evaluate a staged pool with it through grid_evaluate.

loss_spec(shape, prior; Δ, r_max=POOL_FRAME_R_MAX, t_lo, t_hi, NG=1537) -> LossSpec

Build a depth-Δ quadrature-marginal spec for the loss shape under IG scale prior. Z̃ is normalised over the FIXED domain [0, r_max] (the declared outlier support POOL_FRAME_R_MAX, in the pool frame's half-diagonal unit), so the outlier half-width is decoupled from the quadrature. [t_lo, t_hi] are the log-scale grid bounds (use loss_scale_bounds(sq_residuals, n)).

loss_spec(m::ScaleMarginal{K}; NG=_MARG_NG) -> LossSpec

The quadrature of the full-data scale marginal m: shape K at its canonical depth, m's IG prior, Z̃ over [0, m.outlier_halfwidth] and NG log-scale nodes spanning m.t0 ± _MARG_L. The nodes lie on the ABSOLUTE lattice t = i·h, h = 2·_MARG_L/(NG − 1) (the domain snapped by at most h/2), so the specs of every grid point share their nodes and grid_evaluate(::ScaleMarginal, …) sums each model's data term once for the whole grid.

loss_scale_bounds(sq_residuals, n; pad=4.0) -> (t_lo, t_hi)

Data-driven grid bounds on t = log sq_loss_width for a redescending-loss sweep over squared residuals sq_residuals[1:n]. The scale-posterior peak for partition k lands near sq_loss_width ≈ RSS_k/(k−m), which ranges from ≈ the tightest residual scale (small k) to ≈ the all-point scale (k=n); the returned [t_lo, t_hi] brackets that range with pad decades (natural-log) of margin on each side so the trapezoid grid covers every k's peak. Supply to loss_spec(...; t_lo, t_hi).

PriorSpec(family::ScalePriorFamily, fitting::PriorFitting)

Compound spec: the prior family + the algorithm that turns σ̂² data into its parameters. The named specs are the constructor functions JeffreysPriorSpec and IGRobustSpec; the set a sweep iterates is ALL_PRIOR_SPECS.

fit_prior(::InverseGamma, ::MLE, σ²s) -> InverseGammaPrior

MLE of the inverse-gamma shape/scale (α, β). Distributions.jl does not implement fit_mle(InverseGamma, ·) directly, so we use the closed-form transformation

X ~ InverseGamma(α, β)  ⇔  1/X ~ Gamma(shape = α, scale = 1/β)

and call fit_mle(Gamma, 1 ./ x). Returns nothing for n < 2.

fit_prior(::InverseGamma, fit::RobustFit, σ²s) -> InverseGammaPrior

Robust empirical prior: fixed shape α₀ = fit.alpha0, scale β₀ = α₀ · median(σ̂²). The IG prior mean is β₀/(α₀−1); with α₀=1 the mean is undefined but the mode β₀/(α₀+1) = α₀ · median(σ̂²) / 2 anchors near the data median.

Used for HEB-like datasets where per-pair σ̂² estimates are heavy-tailed (low-n_in pairs and non-planar GT-inlier sets contaminate MLE Gamma fits). The median is unaffected by individual outlier σ̂² values, so the prior location stays stable even with single-digit val-pool sizes.

Improper Jeffreys prior — equivalent to PriorSpec(Jeffreys(), NoFit()).

Robust IG prior with shape α₀ and scale fit from median(σ̂²).

ALL_PRIOR_SPECS :: Tuple of canonical (family, fitting) specs probed + evaluated.

Single source of truth for the prior specs a sweep iterates.

ScalePriorFamily

Abstract supertype for the scale-prior distribution family. Concrete subtypes identify the kind of prior (its parameter shape and admissible values); the values come from a PriorFitting applied to σ̂² data.

Jeffreys()

Improper Jeffreys prior — parameterless, equivalent to IGPrior(0, 0).

InverseGamma()

Proper inverse-gamma prior IG(α₀, β₀); hyperparameters come from a fitting.

PriorFitting

Abstract supertype for "how do we derive (α₀, β₀) from the per-pair σ̂² data?" Concrete subtypes implement fit_prior for one or more ScalePriorFamily values.

NoFit()

Parameterless families take NoFit — there's nothing to fit.

RobustFit(alpha0 = 1.0)

Robust median-based fitting of an inverse-gamma prior: β₀ = α₀ · median(σ²s) with fixed shape α₀ (default 1.0). Median is robust to heavy-tailed σ̂² distributions (HEB).

MLE()

Closed-form MLE on the inverse-gamma family via the Gamma reparameterisation 1/X ~ Gamma(α, 1/β). Requires length(σ²s) ≥ 2.

SimFitPP(; alpha = 0.99, max_iter = 4, ftol = 0.01, train_fraction = 0.5,
           min_loss_width = 0.25, max_loss_width = 8.0, dg = 1)

SIMFIT++ configuration (the reference's constructor defaults; the reference calls the loss width its threshold). alpha is the inlier coverage the loss width targets (losswidth² = ppf{χ²dg}(alpha)·σ²); `minlosswidthandmaxloss_widthbound an acceptable loss width in the residual unit of the residuals handed to the estimator (pixels in the reference);dg` is the residual codimension (the χ² degrees of freedom).

simfitpp_loss_width(fit, residuals_sq, n, initial_loss_width, cfg = SimFitPP();
                    num_previous_est = 0, rng = Random.default_rng())
    -> (loss_width, success::Bool)

The reference's SIMFITPP.__call__ over n correspondences. fit(loss_width, train) estimates the geometry from the correspondences indexed by train at loss_width and returns a model handle (or nothing on failure); residuals_sq(model, val) returns the squared residuals of that model on the held-out indices val (or nothing). initial_loss_width seeds the iteration and is returned unchanged when the estimate fails (success == false); with num_previous_est = p > 0 earlier estimates behind initial_loss_width, the successful estimate is filtered as exp(log(loss_width)/(p+1) + p·log(initial_loss_width)/(p+1)).

simfitpp_select(residuals::AbstractMatrix, initial_loss_width, cfg = SimFitPP();
                num_previous_est = 0, rng = default_rng())
    -> (best_idx, loss_width, success)
simfitpp_select(pair::PoolPair, initial_loss_width_px, cfg = SimFitPP(); ...)

SIMFIT++ on a fixed hypothesis pool: residuals is the pool's n×M residual matrix (points × models) in the unit initial_loss_width and the bounds of cfg are expressed in. The geometry "estimator" is the MSAC-best pool hypothesis on the train split at the current loss width, scored by its validation-split residuals; the returned index is the MSAC argmax over ALL points at the final loss width.

The pair method is PIXELS in, pixels out: it runs on pixel_residuals(pair), the pool models' residuals in pixels, so initial_loss_width_px, the bounds of cfg (the reference's pixel defaults, 0.25–8 px) and the returned loss width are all pixels. To work in the pool frame's unit instead, call the matrix method on the pair's stored residual matrix with bounds in that unit (with_residual_unit).

unbiased_median_variance(residuals_sq, sq_loss_width; dg = 1) -> σ²

The reference's UnbiasedMed: the noise VARIANCE from the squared residuals below sq_loss_width, corrected for the truncation. With med the median of the kept residuals and ratio = sq_loss_width / med, the quantile q that the truncated median corresponds to under χ²dg solves q = F(ratio · F⁻¹(q)) / 2, iterated five times from q = 1/2; the estimate is med / F⁻¹(q). For dg == 1 a ratio below 5 has no useful fixed point and the naive med / F⁻¹(1/2) is returned (the reference's `minratio). ReturnsNaN` when no residual lies below the loss width.

It returns σ², not σ — which is why the name no longer says "scale". It was unbiased_median_scale, sitting in the same package as posterior_scale, which returns a LINEAR σ̂. Two exported *_scale names, one squared and one not, is the ambiguity VisualGeometryCore's s2 convention exists to remove.

pixel_residuals(pair::PoolPair) -> Matrix{Float32}

The image-plane Sampson residuals of the pair's stored models in the PIXEL frame, num_points × num_models. Pixels are the frame the SIMFIT++ reference works in, so its pixel bounds apply to them without any scalar approximation (the pool frame scales the two images independently, cf. residual_unit_scale). One n×M scan of the light record: the correspondences come back to pixels through the stored maps (pixel_csponds), each model by composition with the pool normalization's two sides, inv(D) ∘ M ∘ S (a pool F → T₂ᵀ F T₁, a pool H → T₂⁻¹ H T₁). A model type with neither law (a bare matrix) has no sampson_distance and is a MethodError.

It never reads the pool's stored residual matrix, so it does not depend on that matrix's residual kind. The stored residuals of an essential pool are TANGENT Sampson distances (the pool generator's residual_kind = :tangent_sampson), those of an F or H pool image-plane Sampson (:sampson); this exit is image-plane Sampson for every pool, the quantity SIMFIT++'s pixel bounds are stated in. (It used to return the stored matrix itself when both maps were the identity — the image-plane quantity for an F/H pool, but the tangent one for an E pool: two routes to two quantities under one name.)

residual_unit_scale(pair::PoolPair) -> Float64

Pool-frame residual units per pixel for the pair: the mean isotropic scale of the two pool-frame maps (pixel → NormalizedImageFrame). A CONVENTION, not an exact conversion — the pool frame scales the two images independently and the Sampson denominator mixes both images' gradients, so no single scalar relates the two frames per point. Use it only for soft quantities (an initial loss-width guess from a pool-frame-unit grid); anything the bounds gate must run on pixel_residuals.

with_residual_unit(cfg::SimFitPP, unit_scale) -> SimFitPP

The same configuration with the loss-width bounds expressed in another residual unit. min_loss_width/max_loss_width are multiplied by unit_scale (e.g. residual_unit_scale(pair) to carry the reference's pixel bounds into a pool stored in the pool frame). The argument was a bare s, the letter this package otherwise spends on a noise scale.