---
title: "Ordinary and robust summaries with explicit numerical status"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ordinary and robust summaries with explicit numerical status}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Compare targets on the same observations

A barycenter minimizes a weighted sum of squared W2 distances. A median
minimizes a weighted sum of W2 distances, so its numerical objective is in
different units. Compare their fitted distributions and per-observation
residuals; a smaller raw objective across these different targets is not a
model-selection criterion.

```{r point-mass-example}
observations <- ot_collection(lapply(c(0, 0, 0, 2, 8), ot_measure))
ordinary <- ot_summary(observations, target = "barycenter")
robust <- ot_summary(observations, target = "median")
c(barycenter_location = ordinary$estimate$support[1, 1],
  median_location = robust$estimate$support[1, 1])
summary(robust)$per_measure
```

In this deliberately simple example, the majority of input distributions are
point masses at zero. The barycenter lies at two; the median stays at zero.
This analytical illustration is not evidence of robustness for every data
representation, contamination pattern, or local solver.

## State what is allowed to change

| Feasible set | Optimized quantities | Main limitation |
|---|---|---|
| `unrestricted_1d` | Full scalar quantile function on its common refinement | Scalar measures only |
| `free_support_fixed_mass` | A fixed number of candidate locations | Candidate masses stay fixed; optimization is generally local |
| `fixed_support` | Probability masses on supplied locations | Estimate cannot leave the supplied grid |

Here is a small two-dimensional collection of translated four-atom measures.
Three of the five observations coincide. We provide the same support and
candidate masses to both fits so the comparison has the same feasible set.

```{r free-support}
shape <- rbind(c(0, 0), c(1, 0), c(0, 1), c(1, 1))
shifts <- c(0, 0, 0, 2, 8)
clouds <- lapply(shifts, function(s) ot_measure(sweep(shape, 2, c(s, 0), "+")))
ordinary_cloud <- ot_summary(clouds, target = "barycenter",
  feasible = "free_support_fixed_mass", support = shape, mass = rep(1/4, 4))
robust_cloud <- ot_summary(clouds, target = "median",
  feasible = "free_support_fixed_mass", support = shape, mass = rep(1/4, 4))
stopifnot(ordinary_cloud$status == "success", robust_cloud$status == "success")
rbind(barycenter_center = colSums(ordinary_cloud$estimate$support * ordinary_cloud$estimate$mass),
      median_center = colSums(robust_cloud$estimate$support * robust_cloud$estimate$mass))
robust_cloud$diagnostics$optimality_certificate
```

The barycenter is the shape translated by the average shift; the median is the
majority shape. The reported median certificate concerns this feasible
majority input. More general fits need not have such a certificate.

```{r compare-plots, fig.cap="The ordinary summary follows the mean translation; the median retains the majority shape."}
op <- par(mfrow = c(1, 2))
plot(ordinary_cloud$estimate, main = "Barycenter", xlim = c(-.5, 3.5), ylim = c(-.5, 1.5))
plot(robust_cloud$estimate, main = "Median", xlim = c(-.5, 3.5), ylim = c(-.5, 1.5))
par(op)
```

## Initialization and support-size sensitivity

Without supplied locations, the first free-support start uses deterministic
systematic selection from pooled input atom indices weighted by both outer and
within-measure masses. Later starts sample pooled atoms with those probabilities.
This initialization depends on atom order. Supply `support` when the initial
geometry must be controlled.

```{r multiple-starts}
multiple <- ot_summary(clouds, target = "barycenter",
  feasible = "free_support_fixed_mass", support = shape, mass = rep(1/4, 4),
  nstart = 3, seed = 23)
data.frame(start = seq_len(multiple$multistart$nstart),
           objective = multiple$multistart$objectives,
           status = multiple$multistart$statuses,
           reason = multiple$multistart$stop_reasons)
multiple$multistart$chosen
multiple$initialization
```

The selected run has the smallest objective among successful runs. If every
run fails, the best finite failed result is retained with a non-success status.
The local seed restores the caller's prior RNG state. Saving only a seed is
insufficient for an audit: retain the initial supports, masses, controls,
statuses, and fitted objective as well.

Increasing the candidate support size changes the feasible approximation.
The table below reports both the objective and numerical outcome; a failed run
must not be silently dropped from a sensitivity analysis.

```{r support-size}
sizes <- c(1L, 2L, 4L)
size_fits <- lapply(sizes, function(m) ot_summary(clouds, target = "barycenter",
  feasible = "free_support_fixed_mass", control = list(num_support = m),
  nstart = 3, seed = 23))
data.frame(atoms = sizes,
           objective = vapply(size_fits, function(x) x$objective, numeric(1)),
           status = vapply(size_fits, function(x) x$status, character(1)))
```

## Read failures and stopping criteria

`history` contains accepted objective values, starting with initialization.
The plot and printed result show status. Relocation residuals use the currently
selected optimal plans; a small residual does not generally certify global
optimality over moving supports.

```{r deliberately-limited}
limited <- ot_summary(clouds, target = "barycenter",
  feasible = "free_support_fixed_mass", support = shape, mass = rep(1/4, 4),
  control = list(maxiter = 1, alpha = 0.1))
limited[c("status", "stop_reason", "objective")]
plot(limited, type = "history")
```

The warning above is intentional: a single damped update does not meet the
relocation criterion. An `inner_solver_failure` instead means a required
transport subproblem did not succeed. A `step_rejected` outcome means no
acceptable trial update was established. An unresolved collision is a failure
condition, not proof that the current observation is a median. Increasing an
iteration budget may address a limit; it is not a universal repair for other
failure reasons.

A positive median `control$smoothing` changes the target to
`sum(weights * sqrt(W2^2 + smoothing^2))`, reported as `smoothed_median`.
This is different from entropy regularization of an inner transport plan.
Compare runs only when target, smoothing, candidate masses, support size, and
regularization are held fixed or explicitly part of the sensitivity study.

## Fixed-grid summaries

On a prescribed grid, the candidate masses vary instead of the locations.
Initial masses use `control$init_mass` and must be strictly positive for the
mirror-descent implementation. Input measures themselves may contain zeros.

```{r fixed-grid}
fixed <- ot_summary(list(ot_measure(0), ot_measure(3)), weights = c(.75, .25),
  target = "barycenter", feasible = "fixed_support", support = 0:3)
fixed
fixed$diagnostics[c("convex_lower_bound", "objective_gap")]
regularized <- ot_summary(list(ot_measure(0)), target = "barycenter",
  feasible = "fixed_support", support = 0:1, epsilon = 0.5)
regularized
```

The unregularized barycenter objective is convex in the candidate mass vector;
its reported valid dual lower bound can support an objective-gap stopping rule.
The fixed-grid median objective is generally nonconvex and has no corresponding
global convex-gap claim. Regularized median requests are rejected. The
regularized barycenter above optimizes the complete entropy objective, including
the candidate-mass entropy contribution; it is not the unregularized solution.

These tutorials run offline and use only small deterministic examples. Full
application studies, data provenance records, and discretization sensitivity
results are maintained in the manuscript's separate replication bundle.