---
title: "Density-sensitive paths with Fermat distances"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Density-sensitive paths with Fermat distances}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5)
library(dgraphs)
```

## The estimator and its graph restriction

The sample Fermat distance between observed points minimizes the sum of
powered Euclidean step lengths over **all sequences of observed points**:

\[
D_p(i,j)=\min_{\pi:i\leadsto j}\sum_{(a,b)\in\pi}\|x_a-x_b\|^p.
\]

`fermat.distances(points = X, p = 2)` computes this complete-graph reference.
It does not first remove long edges. For `p > 1`, several short steps can cost
less than a single long step. This favors densely sampled routes without
explicitly fitting a density estimate. For `p = 1`, the complete-graph answer
is Euclidean distance, not distance around a circle.

The default result is the unrooted sum. `rooted = TRUE` takes its `1/p` power:
this changes reported distances but not minimizing paths. Neither option
applies sample-size normalization or estimates a limiting constant.

```{r small-example}
X <- matrix(c(0, 1, 2), ncol = 1)
fermat.distances(points = X, p = 2)
```

The distance from 0 to 2 is 2: the two unit steps cost less than the direct
step, which costs 4.

Two explicit restrictions are also available:

```{r restrictions}
fermat.distances(points = X, graph.type = "sknn", k = 1, p = 2)
graph <- create.sknn.graph(X, k = 1, neighbor.method = "exact",
                          prune.edges = FALSE, connect.components = FALSE)
fermat.distances(graph, p = 2)
```

Coordinate kNN construction uses union symmetrization: either endpoint may
select an edge. Supplied graphs use their stored lengths as base lengths before
powering. With the same retained weights, restricting edges cannot decrease
the distance. Missing routes return `Inf`; no connectivity repair is implicit.
Use `vertices = c(i, j)` for a selected submatrix; paths still use other vertices.

## Noisy circles with different sampling densities

We generate 100 observations around a unit circle with radial Gaussian noise
of standard deviation 0.04. One example has uniform angular sampling; the
other places half the observations near angle zero. Both retain a uniform
angular component. These are fixed illustrations, not a simulation study or a
proof of geometric recovery. Radial noise changes the sampled geometry from
an exact circle; circle arc length is therefore not a universal truth target.

```{r circle-data}
set.seed(4101)
n <- 100
angles <- list(
  Uniform = sort(runif(n, 0, 2*pi)),
  Concentrated = sort(c(runif(n/2, 0, 2*pi),
                        rnorm(n/2, 0, 0.4) %% (2*pi)))
)
radial.noise <- rnorm(n, 0, 0.04)
circles <- lapply(angles, function(theta) {
  cbind((1 + radial.noise)*cos(theta), (1 + radial.noise)*sin(theta))
})
powers <- c(1, 2, 3)
reference <- lapply(circles, function(X) {
  lapply(powers, function(p) fermat.distances(points = X, p = p))
})
```

To visualize the minimizing routes, the next chunk uses igraph on the same
complete graph and powered Euclidean weights. It checks each plotted path's
cost against `fermat.distances()`. Black points mark the selected endpoints;
red segments show a shortest path. Tied minimizers may have different shapes.
Notice how the preferred route changes with the power and sampling pattern;
a dense route need not be the shortest circular arc.

```{r paths, fig.height=5.5, fig.cap="Complete-graph shortest paths at powers 1, 2 and 3. Rows show uniform and concentrated angular sampling; axes are the observed Cartesian coordinates.", fig.alt="Six circle panels compare complete-graph paths for three powers and two angular sampling distributions. Endpoints are black, observations gray, and minimizing paths red."}
old <- par(mfrow = c(2, 3), mar = c(3, 3, 2, 0.5))
for (s in seq_along(circles)) {
  X <- circles[[s]]
  ends <- c(which.min(rowSums(sweep(X, 2, c(1,0), "-")^2)),
            which.min(rowSums(sweep(X, 2, c(-1,0), "-")^2)))
  g <- igraph::make_full_graph(n, directed = FALSE)
  e <- igraph::as_edgelist(g, names = FALSE)
  eu <- as.matrix(dist(X))
  for (j in seq_along(powers)) {
    p <- powers[j]
    weights <- eu[e]^p
    path <- as.integer(igraph::shortest_paths(g, from = ends[1], to = ends[2],
                                             weights = weights)$vpath[[1]])
    cost <- sum(eu[cbind(head(path, -1), tail(path, -1))]^p)
    stopifnot(isTRUE(all.equal(cost, reference[[s]][[j]][ends[1],ends[2]])))
    plot(X, asp = 1, pch = 16, cex = 0.5, col = "gray60",
         xlab = "x", ylab = "y", main = paste(names(circles)[s], "p =", p))
    lines(X[path, , drop = FALSE], col = "firebrick", lwd = 2)
    points(X[ends, ], pch = 16, cex = 1.1)
  }
}
par(old)
```

## Does a sparse graph give the same distances?

For each circle and power, compare symmetric kNN restrictions at `k = 4, 8, 16`
against **all unordered pairs** of the complete reference. The table reports
unreachable pairs separately. Relative error is measured only for finite
restricted distances; it does not hide disconnections. Numerical agreement
uses tolerance `1e-10 * max(1, reference distance)`.

```{r comparison}
rows <- list()
for (s in seq_along(circles)) for (j in seq_along(powers)) {
  full <- reference[[s]][[j]]
  use <- upper.tri(full)
  for (k in c(4, 8, 16)) {
    restricted <- fermat.distances(points = circles[[s]], p = powers[j],
                                   graph.type = "sknn", k = k)
    finite <- use & is.finite(restricted)
    stopifnot(all(restricted[use] >= full[use] - 1e-10))
    agree <- is.finite(restricted[use]) &
      abs(restricted[use] - full[use]) <= 1e-10*pmax(1, full[use])
    rows[[length(rows)+1L]] <- data.frame(
      Sampling = names(circles)[s], p = powers[j], k = k,
      Unreachable = sum(use & !is.finite(restricted)),
      Agreement_percent = 100*mean(agree),
      Max_relative_error = if (any(finite))
        max(abs(restricted[finite] - full[finite])/full[finite]) else NA_real_
    )
  }
}
knitr::kable(do.call(rbind, rows), digits = 4)
```

Agreement in this table establishes equivalence only for these finite samples,
powers and neighborhood sizes. In particular, the `p = 1` rows demonstrate
that circle-following kNN distances can differ from the Euclidean reference.
Increasing `k` to `n - 1` always restores the complete graph, at quadratic
storage cost. Merely observing no change between two smaller `k` values does
not certify equivalence to the reference.

## Two noisy arms with diffuse background points

The circle illustrates a change in distance. Two arms ask a different question:
when does density make a path follow populated structures instead of crossing
between them? We use two perpendicular segments of length one, with 80 samples
per arm. Radial coordinates are uniform on `[0,1]`; transverse displacements
are independent absolute Gaussian draws with scale 0.005. There is no supplied
junction vertex. Background points are uniform in the unit square, not placed
on translated rows and not assigned to an arm.

We compare 20 and 120 background points. The smaller set is contained in the
larger, and all arm observations and endpoint identities stay fixed. The seed,
width, background levels and endpoint-selection rule were fixed before inspecting
paths. This is the first seed of a separate 12-seed exploratory comparison at
three widths (0.005, 0.02, 0.05) and three background counts (20, 60, 120), not
a seed selected for its appearance.

The installed teaching helper below contains the generator, deterministic pair
selection, geometric scoring and plotting code. It is a script, not an exported
package API. Its source can be inspected at the path returned by `system.file()`.

```{r arm-data}
source(system.file("examples", "fermat-two-arms.R", package = "dgraphs"))
arm.examples <- lapply(c(20L, 120L), function(count) {
  fermat.two.arms(seed = 6101L, width = 0.005, background = count)
})
arm.routes <- lapply(arm.examples, function(sample) {
  lapply(1:3, function(p) fermat.arm.paths(sample, p))
})
```

We evaluate ten arm-endpoint pairs away from the junction and up to twelve
background pairs from the first twenty background points. Background endpoints
must lie more than 0.08 from either axis and have different nearest axes. Pairs
are ranked by closeness to targets `(0.65,0.18)` and `(0.18,0.65)`, before any
path calculation. Only the first pair of each kind is plotted; the distance
checks below use every evaluated pair. The background endpoints' distances to
the nearest arm are shown here, in units where an arm has length one.

```{r arm-endpoints}
sample <- arm.examples[[1]]
ends <- sample$background.pairs[1, ]
knitr::kable(data.frame(Endpoint = c("Start", "End"),
                       x = sample$X[ends,1], y = sample$X[ends,2],
                       Distance_to_arm = apply(sample$X[ends,,drop=FALSE],1,min)),
             digits = 3)
```

For each polygonal path `gamma`, the junction indicator is
\(J=\mathbf{1}\{\min_{z\in\gamma}\|z\|\leq0.1\}\).
This checks the whole path, not only its vertices. Let
\(T=([0,1]\times[0,0.05])\cup([0,0.05]\times[0,1])\).
Arm coverage is
\(A=\operatorname{length}(\gamma\cap T)/\operatorname{length}(\gamma)\),
using Euclidean segment length, not the powered cost. The score is between zero
and one. Higher coverage means more travel near the arms; neither score is a
universal measure of correctness. Segment intersections with the two rectangles
are computed exactly up to floating-point arithmetic.

Blue points are arm observations, gray crosses are diffuse background, black
circles are endpoints, and red lines are minimizing paths. Shading marks `T`;
the dashed quarter-circle marks the junction disk. All panels use observed
coordinates and complete-graph powered Euclidean costs.

```{r arm-paths, fig.height=5.5, fig.cap="Paths between arm observations. Rows contain 20 and 120 background points; columns show powers 1, 2 and 3. The arm sample and endpoints are fixed across panels.", fig.alt="Six two-arm panels show complete-graph paths between fixed arm endpoints, comparing two nested background levels and three powers."}
old <- par(mfrow = c(2,3), mar = c(3,3,2,0.5))
for (s in 1:2) for (p in 1:3) {
  fermat.plot.arm.path(arm.examples[[s]], arm.routes[[s]][[p]]$paths[[1]],
                      paste(arm.examples[[s]]$background, "background; p =", p))
}
par(old)
```

```{r background-paths, fig.height=5.5, fig.cap="Paths between background observations near different arms, with the same samples and visual encodings. Endpoints are fixed and are not assigned generating-arm labels.", fig.alt="Six panels show paths between two fixed diffuse background observations near different arms, with 20 or 120 background points and powers 1, 2 and 3."}
old <- par(mfrow = c(2,3), mar = c(3,3,2,0.5))
for (s in 1:2) for (p in 1:3) {
  index <- nrow(arm.examples[[s]]$arm.pairs) + 1L
  fermat.plot.arm.path(arm.examples[[s]], arm.routes[[s]][[p]]$paths[[index]],
                      paste(arm.examples[[s]]$background, "background; p =", p))
}
par(old)
```

The next table describes the twelve displayed paths. Junction entry is a yes/no
indicator, and arm coverage is a percentage of the entire path's Euclidean
length, including its endpoint attachments. Thus background paths can follow
the arms between attachments without reaching 100% coverage.

```{r arm-scores}
shown <- do.call(rbind, lapply(1:2, function(s) do.call(rbind,lapply(1:3,function(p) {
  scores <- arm.routes[[s]][[p]]$scores[c(1,11),]
  data.frame(Background=arm.examples[[s]]$background, Power=p,
             Endpoints=scores$type, Junction=ifelse(scores$junction==1,"Yes","No"),
             Arm_coverage_percent=100*scores$arm.fraction)
}))))
knitr::kable(shown, digits=1)
```

These plots separate density preference from obligatory attachment. Adding
background creates potential travel routes between the arms. Such a route is
not automatically an estimator error: the background is part of the sampled
support. In particular, a finite density-sensitive distance need not take
nearby background points to an arm and back.

We also check the same endpoint pairs with symmetric kNN restrictions. Agreement
means absolute error at most `1e-10 * max(1, complete distance)`. The table gives
percent agreement over the selected pairs, separately for `k=8` and `k=16`;
nonfinite distances count as disagreements. The unreachable count sums across
both restrictions. Each igraph path cost is independently checked against the
package's complete-graph calculation.

```{r arm-restriction-check}
arm.agreement <- list()
for (s in 1:2) for (p in 1:3) {
  result <- arm.routes[[s]][[p]]
  vertices <- sort(unique(as.vector(result$pairs)))
  index <- cbind(match(result$pairs[,1],vertices),match(result$pairs[,2],vertices))
  full <- fermat.distances(points=arm.examples[[s]]$X,p=p,vertices=vertices)[index]
  stopifnot(max(abs(full-result$scores$cost)) < 1e-10)
  restricted <- lapply(c(8,16),function(k) {
    fermat.distances(points=arm.examples[[s]]$X,p=p,vertices=vertices,
                     graph.type="sknn",k=k)[index]
  })
  agree <- vapply(restricted,function(d) {
    stopifnot(all(d >= full-1e-10))
    100*mean(is.finite(d) & abs(d-full)<=1e-10*pmax(1,full))
  },0)
  arm.agreement[[length(arm.agreement)+1L]] <- data.frame(
    Background=arm.examples[[s]]$background, Power=p, Pairs=length(full),
    Agreement_k8_percent=agree[1], Agreement_k16_percent=agree[2],
    Unreachable=sum(!is.finite(unlist(restricted))))
}
knitr::kable(do.call(rbind,arm.agreement),digits=1)
```

Any equivalence shown here applies to these selected pairs only. The circle
comparison above checks all pairs. A connected graph or an attractive path
plot does not certify that a restriction preserves Fermat distances.

### How much does the first seed conceal?

The displayed pair follows the arms at powers 2 and 3 under both background
levels. To show its limits, the next calculation repeats power 2 over the
predeclared twelve seeds, at the narrowest and broadest arm widths and both
background levels. Each seed contributes its mean over ten arm pairs or twelve
background pairs; endpoint pairs are not treated as independent replicates.
The table gives means and 95% percentile bootstrap intervals from 2,000 seed
resamples. Width is the pre-folding Gaussian scale; both path scores are
percentages. Degenerate intervals do not imply certainty for future samples.

```{r arm-seed-summary}
seed.scores <- list()
for (seed in 6101:6112) for (width in c(0.005,0.05)) for (count in c(20L,120L)) {
  sample <- fermat.two.arms(seed, width, count)
  scores <- fermat.arm.paths(sample, p=2)$scores
  means <- aggregate(scores[c("junction","arm.fraction")],list(Endpoints=scores$type),mean)
  means$Seed <- seed; means$Width <- width; means$Background <- count
  seed.scores[[length(seed.scores)+1L]] <- means
}
seed.scores <- do.call(rbind,seed.scores)
set.seed(8201)
mean.interval <- function(x) {
  boot <- replicate(2000,mean(sample(x,replace=TRUE)))
  sprintf("%.1f (%.1f–%.1f)",100*mean(x),
          100*quantile(boot,.025),100*quantile(boot,.975))
}
seed.summary <- do.call(rbind,lapply(split(seed.scores,
  interaction(seed.scores[c("Width","Background","Endpoints")],drop=TRUE)),function(d) {
    data.frame(d[1,c("Width","Background","Endpoints")],
               Junction_percent=mean.interval(d$junction),
               Arm_coverage_percent=mean.interval(d$arm.fraction))
}))
rownames(seed.summary) <- NULL
knitr::kable(seed.summary)
```

For narrow arms, every tested arm-to-arm path at power 2 visits the junction,
whereas background-to-background visits average only 45.1% with 20 background
points and 22.2% with 120. Wider arms weaken junction routing further. The
fixed width-0.05 coverage rectangles contain less of the broad clouds, so
falling coverage alone is not evidence of additional shortcuts.

The first-seed distance-agreement table likewise establishes only what was
checked in that realization. It is not a guarantee for subsequent datasets.

## Complete-graph distances without storing the graph

`backend = "implicit"` retains the complete-graph estimator but computes edge
costs during dense Dijkstra searches. It avoids quadratic adjacency storage.
Selecting source vertices (pivots) also avoids an all-pairs output matrix;
every observation can still occur inside a minimizing path.

```{r implicit-pivots}
X <- circles[[1]]
pivots <- c(1L, 14L, 27L, 40L, 53L, 66L, 79L, 92L)
pivot.distances <- fermat.distances(points = X, p = 2,
    backend = "implicit", sources = pivots)
stopifnot(isTRUE(all.equal(pivot.distances,
                          reference[[1]][[2]][pivots, , drop = FALSE])))
dim(pivot.distances)
# Rectangular queries preserve the supplied row and column order.
fermat.distances(points = X, backend = "implicit", sources = pivots[1:2],
                targets = c(25L, 50L, 75L))
```

With h sources, n observations and d coordinates, worst-case work is
O(h n² d), plus an O(n² d) numeric-range check. Memory is O(n d + h n)
for source-to-all queries. Requesting all sources still produces an n-by-n
matrix. `max.workspace.bytes` checks an estimate of input, output and native
scratch storage; it is not a cap on total R process memory. This backend trades
repeated edge evaluation for less storage; it does not promise subquadratic
search time or approximate the missing distances.

### Sparse metric MDS with consistent Fermat targets

The teaching helper below constructs local symmetric kNN pairs and pivot
interactions. Both receive **complete-graph Fermat distances**. In particular,
a local pair's target can be smaller than its powered direct-edge length.
The helper uses nearest-pivot regions and half-distance region counts for
asymmetric long-range weights, with unit endpoint multiplicities for local
pairs, following the weighting construction used by grip's sparse mode.
Pivots here are fixed observation indices, not graph-distance random samples.

```{r fermat-constraints}
source(system.file("examples", "fermat-sparse-mds.R", package = "dgraphs"))
constraints <- fermat.mds.constraints(X, pivots, k = 4, p = 2)
stopifnot(max(abs(constraints$targets -
                 reference[[1]][[2]][constraints$pairs])) < 1e-10)
c(retained = nrow(constraints$pairs), all.pairs = choose(nrow(X), 2))
```

With a grip version supporting explicit `constraints`, fit these targets as
follows. Sparse fitting approximates the embedding objective; the supplied
Fermat targets themselves are exact sample distances, up to arithmetic.

```{r sparse-fermat-fit, eval=FALSE}
fit <- grip::metric.mds(n = nrow(X), approximation = "sparse",
                       constraints = constraints, seed = 1, max.iter = 100)
plot(fit$coords, asp = 1, pch = 16)
```

The helper is a readable teaching implementation, not a fast neighbor-search
API. Its local targets require additional source searches with early stopping;
worst-case preparation can still be cubic in n. It never allocates an all-pairs
distance matrix. Inverse-squared MDS weights require positive targets, so
consolidate coincident observations first; the Fermat distance function itself
continues to support them. Evaluate embedding fidelity on independent pairs,
not only the constraints used for fitting.

## What the theory does and does not establish

Groisman, Jonckheere and Sapienza define the unrooted complete-graph estimator
and establish a normalized continuum limit under their sampling assumptions.
They also define directed nearest-neighbor restrictions and prove a
high-probability equality with the unrestricted estimator when `k` exceeds
a bound of the form `c1 * log(n / epsilon) + c2`, under the stated hypotheses
(Definitions 2.1 and 2.11, Proposition 2.12). This does not select a certified
finite-sample `k` for these noisy examples.

Union-symmetric kNN contains the directed edges, so equality for the directed
restriction implies equality for the union restriction by inclusion of paths.
Arbitrary pruning, mutual-neighbor intersections and other supplied graphs
need separate justification. A rooted output is a transformation of the
unrooted estimator and must be interpreted with its corresponding scaling.

This example tests graph restriction directly rather than assuming it is
harmless. Further assessment should use multiple sample sizes and seeds, and
should distinguish recovering ordinary geometry from deliberately changing
distances through sampling density. Neither a visually appealing path nor a
connected graph establishes the continuum theorem's assumptions.

The explicit complete-graph backend and all-pairs output require quadratic
storage. The implicit backend avoids adjacency storage. Both `fermat.distances()` coordinate graph
constructions still evaluate all Euclidean pairs; kNN restriction reduces the
stored graph, not that neighbor-search work. Start with small reference examples. Positive lengths that overflow or underflow
to zero after powering are rejected; rescale the coordinates if necessary.

## Exact nearest neighbors without all-pairs distances

When the desired output is a neighborhood graph, `fermat.knn()` finds each
observation's nearest neighbors in the complete-graph Fermat distance. It caches
exact Euclidean neighbor lists and uses truncated Dijkstra searches. This is a
separate operation from asking `fermat.distances(graph.type = "sknn")` for a
matrix of arbitrary pairwise distances.

McKenzie and Damelin's Lemma 5.6 justifies pruning for these local queries:
the neighbor distances are those of the complete graph. It does not assert that
the resulting sparse graph preserves every pairwise distance. Our `k` counts
other observations, while the paper includes the source in its count.

```{r exact-fermat-neighbors}
set.seed(193)
neighbor.points <- matrix(runif(120), 60, 2)
neighbors <- fermat.knn(neighbor.points, k = 6, p = 2)
neighbor.graph <- create.fermat.sknn.graph(neighbor.points, k = 6, p = 2)
reference <- fermat.distances(points = neighbor.points, p = 2)
selected <- cbind(rep(seq_len(nrow(neighbor.points)), 6),
                  as.vector(neighbors$index))
stopifnot(max(abs(as.vector(neighbors$distance) - reference[selected])) < 1e-10)
```

The graph takes the union of directed selections: an edge remains if either
endpoint selected the other. Its stored lengths are the computed Fermat
distances. Do not power those lengths again when calculating graph geodesics.
Union symmetrization can give a vertex more than `k` neighbors; connectivity is
not forced. Equal-distance ties at the cutoff may produce different valid
neighbor selections. Coincident observations remain separate vertices.

Candidate storage grows as `n * k`. Once candidates are available, binary-heap
path searches cost approximately `n * k^2 * log(k + 1)`. Exact Euclidean search
can still be expensive in high dimensions. The search offers exact k-d tree
and brute-force options; neither uses approximate neighbors. Pivot-to-all
queries still use the separate implicit backend.

## References

Groisman, P., Jonckheere, M. and Sapienza, F. (2022).
*Nonhomogeneous Euclidean first-passage percolation and distance learning.*
Bernoulli 28(1), 255--276. [Published paper](https://doi.org/10.3150/21-BEJ1341)
and [preprint with definitions and propositions](https://arxiv.org/html/1810.09398v2).

For ordinary graph distances and graph construction, see
[Constructing and Diagnosing Data-Derived Graphs](data-derived-graph-workflow.html).

McKenzie, D. and Damelin, S. (2019). *Power Weighted Shortest Paths for
Clustering Euclidean Data.* [Preprint, section 5](https://arxiv.org/abs/1905.13345).