Migrating analyses to the 0.2.0 development interface

Numerical audit repairs in 0.2.0.9002

The exact network-simplex backend now normalizes active transport costs before optimization. Returned objectives and potentials retain their input units; ot_inspect() exposes the normalized dual checks and their tolerance. This repairs failures at small coordinate scales, but does not establish relative accuracy for an arbitrarily tiny optimum within a large cost matrix.

Histogram grids must have exactly equal numerical break values. A reused histogram estimate supplies its display grid, while inspection describes the new computation. Scalar covariance arrays work through gaussmedpd(), and PW retains successful initializations when another attempted start fails. The attempts and their individual statuses remain available in initialization metadata.

Sinkhorn reports cost + lambda * sum(plan * (log(plan) - 1)); omitting the minus-one term shifts a unit-mass objective by lambda. Overflow of this full objective is a numerical failure even when marginal_converged is true. Density quadrature scales spacings before addition to avoid unnecessary overflow.

Quantile summaries preserve selected input measures and evaluate final diagnostics on the returned probability representation. For free-support summaries, coincident candidate atoms can leave a zero relocation residual without a local minimum. Such unresolved cases now report unresolved_degeneracy instead of success. Inspect the returned fit and retry with distinct initial locations when appropriate; this status change does not add a global optimality guarantee.

The frozen 0.2.0.9001 manuscript archive and saved results have not been replaced. Reproducing that analysis and testing the repaired development version are separate workflows.

Additive toolbox interfaces in 0.2.0.9001

ot_inspect(result) provides a common view of existing numerical metadata; it leaves the native result unchanged. It accepts individual comparisons, barycenters, medians and interpolation results. For several interpolation times use lapply(path, ot_inspect). Constructors, display helpers and wassboot() are outside this inspection interface. Older saved objects that lack identifying metadata may need to be regenerated with their recorded inputs and controls.

swdist() and pwdist() additionally accept ot_measure() objects. Their stored masses replace raw-input mass arguments; providing both is an error. pwbary() accepts a collection or a list of measures. Explicit outer weights override collection weights; weights = NULL retains them.

ot_gaussian_summary() accepts Gaussian parameter objects, with controls maxiter, atol, and rtol. Multidimensional medians additionally accept inner_maxiter, inner_atol, and inner_rtol. These names map to the existing algorithms; legacy gauss* functions retain their original control names and defaults. Analytic branches report controls they did not use.

Recheck the target and numerical outcome

The development release keeps historical entry-point names but makes several correctness repairs and changes some implementations. Existing scripts may produce different answers or reject controls that used to be ignored. Preserve an old analysis and its package version when comparing results across releases. A matching function name does not establish an identical numerical method.

Previous workflow	Explicit current workflow
wasserstein(x, y, p = 2)	ot_distance(ot_measure(x), ot_measure(y), p = 2, return_plan = TRUE)
Cross-distance matrix	ot_cost(D, type = "distance", metric = TRUE) passed to ot_distance()
Already powered cost matrix	ot_cost(C, type = "cost"); do not power again
ECDF/histogram summary	Convert with as_ot_measure() and use ot_summary() when the full finite estimate is needed
Free-support summary	Set feasible = "free_support_fixed_mass"; separate initial locations from fixed candidate masses
Fixed-grid summary	Set feasible = "fixed_support" and supply locations; use control$init_mass for initialization
Gaussian parameters	Use ot_gaussian()/ot_gaussian_distance() and the Gaussian-family summary functions

Within-measure mass and collection weights are separate. Nonnegative values are normalized, but negative, missing, nonfinite, or all-zero masses are errors. Unknown controls, malformed dimensions, invalid iteration budgets, and unsupported combinations are rejected. Common interfaces take a named control list; legacy interfaces generally retain their documented ... controls, including names such as abstol instead of atol.

x <- c(0, 2); y <- c(0, 3)
old_entry <- wasserstein(x, y, p = 2)
new_entry <- ot_distance(ot_measure(x), ot_measure(y), p = 2, return_plan = TRUE)
stopifnot(old_entry$status == "success", new_entry$status == "success")
c(legacy = old_entry$distance, common = new_entry$wasserstein_distance)
#>    legacy    common 
#> 0.7071068 0.7071068

Inspect status and stop_reason rather than inferring convergence from a returned number. An exact solver can return a feasible limited plan without proving optimality. Entropic transport exposes its transport component and regularized objective separately. A rooted transport component from an entropic plan is not presented as an exact Wasserstein metric.

ECDF summaries preserve weighted finite distributions

ecdfbary() and ecdfmed() still return callable CDFs. They now have class t4transport_cdf and inherit stepfun, rather than claiming to be base ecdf objects built from invented repeated samples. Cumulative jumps retain ties and unequal sample sizes. Their finite quantiles use the generalized-inverse type = 1 convention; interpolation types are rejected.

F <- ecdf(c(0, 0, 2))
G <- ecdf(c(1, 3))
fit_cdf <- ecdfbary(list(F, G))
class(fit_cdf)
#> [1] "t4transport_cdf" "stepfun"         "function"
fit_cdf(c(0, 1, 2, 3))
#> [1] 0.0000000 0.5000000 0.6666667 1.0000000
quantile(fit_cdf, probs = c(0, .5, 1), type = 1)
#>   0%  50% 100% 
#>  0.5  0.5  2.5
as_ot_measure(fit_cdf)
#> Finite probability measure: 3 atoms in 1 dimension(s)
#> Representation: finite_cdf | total mass: 1
attr(fit_cdf, "t4transport_fit")[c("objective", "status", "stop_reason")]
#> $objective
#> [1] 0.5833333
#> 
#> $status
#> [1] "success"
#> 
#> $stop_reason
#> [1] "quantile_average"

Use as_ot_measure() instead of inspecting a CDF closure’s private environment. The fitted CDF can be passed into subsequent ecdfbary() or ecdfmed() calls. Median collisions are checked using the appropriate combined weight and subgradient; being close to one input is not automatically a stopping rule.

Histogram displays and estimates are separate

Histograms represent discrete probability masses at bin midpoints. For ordinary input histograms, normalized counts determine masses, including with unequal bin widths. Fitted histograms store probabilities in both counts and mass; density is mass divided by bin width. effective_sample_size is NA because an estimated distribution has no implied number of observations. Do not round these probabilities into artificial counts.

breaks <- c(0, 1, 3)
h1 <- hist(c(.2, .4, 2), breaks = breaks, plot = FALSE)
h2 <- hist(c(.5, 2, 2.5), breaks = breaks, plot = FALSE)
hfit <- histbary(list(h1, h2))
data.frame(midpoint = hfit$mids, mass = hfit$mass, density = hfit$density)
#>   midpoint      mass   density
#> 1      0.5 0.3333333 0.3333333
#> 2      2.0 0.6666667 0.3333333
stopifnot(isTRUE(all.equal(hfit$density * diff(hfit$breaks), hfit$mass)))
as_ot_measure(hfit)
#> Finite probability measure: 2 atoms in 1 dimension(s)
#> Representation: histogram_midpoints | total mass: 1
as_ot_measure(hfit$unprojected_measure)
#> Finite probability measure: 3 atoms in 1 dimension(s)
#> Representation: empirical | total mass: 1

histbary(), histmed(), and histinterp() retain the unprojected finite result in $unprojected_measure. The histogram display bins this result onto the original breaks and represents the resulting masses at bin midpoints. It can differ from the unrestricted summary or displacement geodesic. Their objective and history refer to the unprojected fit. Later histogram operations use the displayed midpoint measure consistently; use the separate unprojected measure explicitly if that is the intended input. The compatibility argument L is validated but no longer controls a probability grid for these finite operations.

Historical regularized names now share one checked implementation

fbary14C(), fbary15B(), their dist variants, and corresponding histogram and image wrappers use objective-checked mirror descent with log-domain Sinkhorn gradients. The historical names no longer identify the original accelerated or Bregman-projection iteration. Numerical trajectories, iteration counts, and results at a fixed small budget may therefore change.

The target has fixed candidate locations and variable candidate masses:

\[ \sum_i \alpha_i \min_{P_i} \{\langle D_i^p,P_i\rangle+ \lambda\sum_{jk}(P_i)_{jk}(\log(P_i)_{jk}-1)\}. \]

Input distances are powered exactly once. The log(P)-1 convention differs from log(P) by a constant for normalized plans, giving the same minimizer but a different reported objective. It is not interchangeable with a KL penalty relative to candidate mass times input mass when candidate masses vary.

probability <- fbary14C(support = 0:1, atoms = list(0), lambda = 0.5)
probability
#> [1] 0.8807971 0.1192029
#> Status: SUCCESS ( single_target_columnwise_entropy_optimum )
attr(probability, "diagnostics")[c("method", "status", "objective", "entropy_convention")]
#> $method
#> [1] "objective_checked_mirror_descent_log_sinkhorn"
#> 
#> $status
#> [1] "success"
#> 
#> $objective
#> [1] -0.563464
#> 
#> $entropy_convention
#> [1] "sum(plan * (log(plan) - 1))"

Initial masses supplied as init.vec, or as the corresponding image initial mass matrix, must be strictly positive. Input measures may contain zero mass. Controls for the shared solver include rtol, step0, max_backtrack, inner_maxiter, and inner_tol; return_result = TRUE is available where specified in the function help. Histogram/image wrappers exposing nthread accept only nthread = 1. The convex bound and actual accepted objectives, rather than iteration count alone, support interpretation of a completed fit. Choose lambda explicitly for reproducibility: legacy defaults can depend on the support scale.

Image methods have explicit grids and projection error

Image intensities are normalized probability masses. Three coordinate conventions are intentionally distinguished:

Converter or method	Coordinates	Ordering
ot_image()	Default x = column - 1/2; y = number of rows - row + 1/2	Row-major; x right, y up
imagedist() and image summary/interpolation wrappers	Columns span [0,1]; rows go from y = 1 to y = 0	Row-major
img2measure()	x = row - 1/2; y = column - 1/2, in pixel units	Legacy base image() orientation; nonzero pixels

A singleton dimension uses the corresponding first endpoint in the legacy endpoint grid. Transport values from different coordinate conventions need not agree, particularly on rectangular images. Normalize units and orientation deliberately before comparing or pooling representations. img2measure() reports mass removed by optional thresholding; threshold = FALSE preserves positive masses.

imagemed() now optimizes the sum of W2 distances over masses on the prescribed grid, using the shared objective-checked selected-subgradient solver. This mass-space problem is generally nonconvex. The old denominator-floor IRLS and its delta/bary.* controls have been removed; these arguments are errors. Inspect the returned matrix’s diagnostics attribute, or request the full result with return_result = TRUE. An arbitrary supplied cost matrix changes the interpretation to a sum of rooted transport costs without a metric robustness guarantee.

imageinterp() now defaults to eps = 0, preserving input zero masses. A positive floor perturbs the normalized inputs and its L1 change is recorded. The unprojected measure follows displacement interpolation of an exact optimal plan. Nearest-grid deposition is a separate approximation.

left <- matrix(c(1, 0), 1, 2)
right <- matrix(c(0, 1), 1, 2)
midpoint <- imageinterp(left, right, t = .5)
attr(midpoint, "unprojected_measure")
#> Finite probability measure: 1 atoms in 2 dimension(s)
#> Representation: displacement_interpolation | total mass: 1
attr(midpoint, "diagnostics")[c("projection_w2_upper_bound",
  "input_mass_perturbation_l1", "displayed_image_is_exact_geodesic")]
#> $projection_w2_upper_bound
#> [1] 0.5
#> 
#> $input_mass_perturbation_l1
#> [1] 0 0
#> 
#> $displayed_image_is_exact_geodesic
#> [1] FALSE

The displayed two-pixel image cannot represent a point halfway between pixels; the retained unprojected measure can. The reported bound is the W2 cost of the specified deposition coupling, not a claim that the deposition is the optimal projection onto grid-supported distributions.

Gaussian means, variances, and covariance domains

The Gaussian summaries preserve list fields mean and var and add actual objective history, status, controls, and stopping diagnostics. Univariate barycenters average means and standard deviations, then square the latter. Univariate medians solve a joint geometric-median problem in those coordinates; separate coordinate medians are generally different.

Zero variance is supported in gaussbary1d(), gaussmed1d(), and univariate ot_gaussian(). The pd summary functions require a finite n by p mean matrix and an exactly matching p by p by n array of symmetric strictly positive-definite covariances. The multivariate constructor has the same strict SPD domain. No automatic symmetrization of user input, ridge, or PSD projection is performed.

General SPD median updates now use arithmetic reweighted means, which are the mean minimizers of their squared-distance subproblems. A collision without a valid certificate is reported as unresolved, and inner failures propagate. Numerical iterate/objective stopping does not certify a global Gaussian-family median. The dedicated ot_gaussian_distance() compares parameters directly and does not create simulated empirical measures.