--- title: "Migrating analyses to the 0.2.0 development interface" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Migrating analyses to the 0.2.0 development interface} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", cache = FALSE) library(T4transport) ``` ## 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`. ```{r comparison-migration} 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) ``` 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. ```{r cdf-migration} F <- ecdf(c(0, 0, 2)) G <- ecdf(c(1, 3)) fit_cdf <- ecdfbary(list(F, G)) class(fit_cdf) fit_cdf(c(0, 1, 2, 3)) quantile(fit_cdf, probs = c(0, .5, 1), type = 1) as_ot_measure(fit_cdf) attr(fit_cdf, "t4transport_fit")[c("objective", "status", "stop_reason")] ``` 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. ```{r histogram-migration} 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) stopifnot(isTRUE(all.equal(hfit$density * diff(hfit$breaks), hfit$mass))) as_ot_measure(hfit) as_ot_measure(hfit$unprojected_measure) ``` `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. ```{r legacy-entropy} probability <- fbary14C(support = 0:1, atoms = list(0), lambda = 0.5) probability attr(probability, "diagnostics")[c("method", "status", "objective", "entropy_convention")] ``` 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. ```{r image-interpolation} left <- matrix(c(1, 0), 1, 2) right <- matrix(c(0, 1), 1, 2) midpoint <- imageinterp(left, right, t = .5) attr(midpoint, "unprojected_measure") attr(midpoint, "diagnostics")[c("projection_w2_upper_bound", "input_mass_perturbation_l1", "displayed_image_is_exact_geodesic")] ``` 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.