---
title: "Distribution representations and transport objectives"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Distribution representations and transport objectives}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", cache = FALSE)
library(T4transport)
```

## One observation is a probability distribution

Use `ot_measure()` for finite weighted atoms. A numeric vector gives scalar
locations; a matrix gives one multivariate location per row. Probability masses
within an observation and weights across observations answer different questions.

```{r representations}
mu <- ot_measure(c(0, 2), mass = c(1, 3))
nu <- as_ot_measure(ecdf(c(0, 0, 4)))
data.frame(location = mu$support[, 1], probability = mu$mass)
observations <- ot_collection(list(first = mu, second = nu), weights = c(2, 1))
observations
```

Here `mu` places probabilities 1/4 and 3/4 at its two locations. The collection
gives the *whole distribution* `mu` weight 2/3. Duplicating every observation in
an empirical sample leaves that empirical probability distribution unchanged;
it does not automatically increase its outer weight. Input masses and outer
weights must be finite, nonnegative, and have positive total. Zero masses and
duplicate support rows retain their original indices in a requested plan.

## A cost and a distance have different units

For order `p`, the Euclidean transport cost is distance raised to `p`, and the
Wasserstein distance is the `p`th root of the optimal cost. With `p = 2`, support
in meters gives a cost in squared meters and a distance in meters.

```{r exact-comparison}
comparison <- ot_distance(mu, nu, p = 2, return_plan = TRUE)
comparison
stopifnot(comparison$status == "success")
c(row_error = max(abs(rowSums(comparison$plan) - mu$mass)),
  column_error = max(abs(colSums(comparison$plan) - nu$mass)))
```

Automatic method selection uses monotone finite transport for scalar supports
and network simplex otherwise. A numerical answer is an exact *optimization
target* evaluated in floating-point arithmetic. Inspect status and residuals
before treating a computed plan as optimal.

A supplied cross-matrix must declare whether it contains distances or already
powered costs. Its rows and columns refer to the original source and target
indices. A metric declaration is supplied by the caller; the package cannot
verify a triangle inequality from one cross-matrix alone.

```{r supplied-cost}
D <- abs(outer(mu$support[, 1], nu$support[, 1], "-"))
from_distance <- ot_distance(mu, nu, p = 2,
                            cost = ot_cost(D, type = "distance", metric = TRUE))
from_cost <- ot_distance(mu, nu, p = 2,
                        cost = ot_cost(D^2, type = "cost", metric = TRUE))
stopifnot(isTRUE(all.equal(from_distance$transport_cost, from_cost$transport_cost)))
```

For arbitrary nonmetric costs, keep `metric = FALSE`. The result then reports
a transport cost without claiming that it defines a Wasserstein metric.

## Entropy regularization changes the objective

With plan entries `P`, the Sinkhorn target is

\[
\langle C,P\rangle + \varepsilon\sum_{ij} P_{ij}(\log P_{ij}-1),
\]

where a zero entry contributes zero. `epsilon` has the same units as the
powered cost. Scaling coordinates by a positive factor `s` requires scaling
`epsilon` by `s^p` to preserve the corresponding regularized problem.

```{r regularized-comparison}
regularized <- ot_distance(mu, nu, method = "sinkhorn", epsilon = 0.5)
regularized
stopifnot(regularized$status == "success")
regularized[c("transport_cost", "transport_cost_root", "regularized_objective")]
```

The transport component is evaluated at the regularized plan. Its root is not
claimed to be an exact Wasserstein distance or a metric. The full regularized
objective can be negative, so taking its square root is not a median construction.
`return_plan = FALSE` reduces the returned object, but internal dense allocations
may still be needed.

## Scalar barycenters and medians

The barycenter minimizes `sum(weights * W2^2)`. The median minimizes
`sum(weights * W2)`. Both use inner Wasserstein order two. The second problem is
a geometric median of complete quantile functions, not separate pointwise
quantile medians.

```{r scalar-summaries}
barycenter <- ot_summary(observations, target = "barycenter")
median <- ot_summary(observations, target = "median")
barycenter
median
stopifnot(barycenter$status == "success", median$status == "success")
data.frame(location = barycenter$estimate$support[, 1],
           mass = barycenter$estimate$mass)
```

The finite engine combines all cumulative-mass breakpoints and integrates the
resulting quantile intervals exactly, up to floating-point arithmetic. It needs
neither a uniform probability grid nor simulated repeated samples. The median
uses modified Weiszfeld updates and verifies a collision before accepting an
input as a solution. An unrepresentable positive mass interval produces an
explicit error rather than silently disappearing.

## Choose a conversion deliberately

| Input | Conversion | Meaning |
|---|---|---|
| Samples or weighted atoms | `ot_measure()` | Finite atomic probability measure |
| Base ECDF or fitted weighted CDF | `as_ot_measure()` | Exact finite jumps, preserving ties |
| Histogram | `as_ot_measure()` | Normalized bin counts at bin midpoints |
| Density values on a grid | `ot_density()` | Trapezoidal node masses, a finite approximation |
| Image | `ot_image()` | Normalized pixel masses at explicit pixel centers |
| Gaussian parameters | `ot_gaussian()` | A separate parametric Gaussian object |

A histogram density height is not a probability mass for unequal-width bins.
The input converter uses counts; displayed fitted densities equal mass divided
by bin width. Histogram summary displays can differ from the separately
retained unprojected finite estimate.

```{r density-and-image}
grid <- seq(-3, 3, length.out = 31)
density_measure <- ot_density(grid, dnorm(grid))
image_measure <- ot_image(matrix(c(1, 0, 2, 1, 0, 0), 2, 3, byrow = TRUE))
density_measure$metadata$approximation
image_measure$metadata[c("x", "y", "order")]
```

The density calculation normalizes the discretized curve on the supplied grid;
truncation and grid refinement remain modeling choices. Image intensities are
pixel masses, even with nonuniform coordinates. Normalizing discards total
intensity. Default image coordinates increase rightwards and upwards, with
row-major indexing. Legacy image functions use other documented coordinate
conventions; see the migration tutorial.

## Work directly with Gaussian parameters

```{r gaussian}
g1 <- ot_gaussian(mean = 0, covariance = 1)
g2 <- ot_gaussian(mean = 4, covariance = 9)
ot_gaussian_distance(g1, g2)
gaussbary1d(c(0, 4), c(1, 9), weights = c(1, 3))[c("mean", "var", "objective", "status")]
gaussmed1d(c(0, 1, 2), c(4, 16, 9))[c("mean", "var", "objective", "status")]
```

Univariate Gaussian W2 geometry is Euclidean geometry of `(mean, standard
deviation)`. Consequently the barycenter averages these two parameters and
squares the averaged standard deviation to recover variance. The median is
their joint geometric median. Zero univariate variance represents a point mass.

For multivariate summaries, `gaussbarypd()` and `gaussmedpd()` take an `n` by `p`
mean matrix and a `p` by `p` by `n` covariance array. Every covariance must be
symmetric and strictly positive definite, including when using these `pd`
functions with `p = 1`. No ridge or PSD projection is applied. General covariance
medians use reweighted Gaussian barycenter subproblems and report numerical
stopping criteria without a global-optimality claim. These parameter objects
are not automatically converted into random finite samples.