enable & expose ground distance power

NEWS.md

@@ -1,5 +1,16 @@
 # phutil (development version)
 
+## New features
+
+- Activation of the `internal_p` parameter in Hera's Wasserstein distance
+  function for alternative Minkowski ground distances with a binding to the new
+  argument `q` for `wasserstein_distance()` and
+  `wasserstein_pairwise_distances()`. The implementation is efficient for
+  `p < 6` (roughly); higher values of `p` are accepted with a warning and may
+  cause R to stall. The bottleneck distance is deployed only when
+  `q == p == Inf`; Hera's implementation does not provide support for
+  `internal_p`.
+
 # phutil 0.0.2
 
 ## New features

R/cpp11.R

@@ -8,10 +8,10 @@ bottleneckPairwiseDistances <- function(x, delta, ncores) {
   .Call(`_phutil_bottleneckPairwiseDistances`, x, delta, ncores)
 }
 
-wassersteinDistance <- function(x, y, delta, wasserstein_power) {
-  .Call(`_phutil_wassersteinDistance`, x, y, delta, wasserstein_power)
+wassersteinDistance <- function(x, y, delta, wasserstein_power, internal_p) {
+  .Call(`_phutil_wassersteinDistance`, x, y, delta, wasserstein_power, internal_p)
 }
 
-wassersteinPairwiseDistances <- function(x, delta, wasserstein_power, ncores) {
-  .Call(`_phutil_wassersteinPairwiseDistances`, x, delta, wasserstein_power, ncores)
+wassersteinPairwiseDistances <- function(x, delta, wasserstein_power, internal_p, ncores) {
+  .Call(`_phutil_wassersteinPairwiseDistances`, x, delta, wasserstein_power, internal_p, ncores)
 }

R/distances.R

@@ -8,26 +8,37 @@
 #'
 #' A matching \eqn{\varphi : D_1 \to D_2} between persistence diagrams is a
 #' bijection of multisets, where both diagrams are assumed to have all points on
-#' the diagonal with infinite multiplicity. The _\eqn{p}-Wasserstein distance_
+#' the diagonal with infinite multiplicity. The _\eqn{pq}-Wasserstein distance_
 #' between \eqn{D_1} and \eqn{D_2} is defined as the infimum over all matchings
 #' of the expression
 #'
-#' \deqn{W_p(D_1,D_2) = \inf_{\varphi: D_1 \to D_2}
-#' \left( \sum_{x \in D_1}{\lVert x - \varphi(x) \rVert^p}
+#' \deqn{W^q_p(D_1,D_2) = \inf_{\varphi: D_1 \to D_2}
+#' \left( \sum_{x \in D_1}{{\lVert x - \varphi(x) \rVert_q}^p}
 #' \right)^{\frac{1}{p}}}
 #'
 #' that can be thought of as the Minkowski distance between the diagrams viewed
 #' as vectors on the shared coordinates defined by the matching \eqn{\varphi}.
-#' The norm \eqn{\lVert \cdot \rVert} can be arbitrary; as implemented here, it
-#' is the infinity norm \eqn{\lVert (x_1,x_2) \rVert_\infty = \max(x_1,x_2)}. In
-#' the limit \eqn{p \to \infty}, the Wasserstein distance becomes the
-#' _bottleneck distance_:
+#' The norm \eqn{\lVert \cdot \rVert_q} is the Minkowski metric with exponent
+#' \eqn{q}, used to measure distances in the plane. In the limit \eqn{p \to
+#' \infty}, the \eqn{pq}-Wasserstein distance becomes the \eqn{q}-_bottleneck
+#' distance_:
 #'
 #' \deqn{B(D_1,D_2) = \inf_{\varphi: D_1 \to D_2}
 #' \sup_{x \in D_1}{\lVert x - \varphi(x) \rVert}.}
 #'
 #' The Wasserstein metric is also called the Kantorovich metric in recognition
-#' of the originator of the metric.
+#' of its originator.
+#'
+#' phutil provides support for the \eqn{\infty}-bottleneck distance and for
+#' \eqn{pq}-Wasserstein distances satisfying \eqn{p,q < \infty}, through a
+#' slight modification of and bindings to the Hera C++ library (Kerber, Morozov,
+#' & Nigmetov, 2017).
+#'
+
+#' @references Kerber, M., Morozov, D., & Nigmetov, A. (2017). Geometry Helps to
+#'   Compare Persistence Diagrams. ACM J. Exp. Algorithmics, 22, 1.4:1-1.4:20.
+#'   https://doi.org/10.1145/3064175
+
 #'
 #' @param x Either a matrix of shape \eqn{n \times 2} or an object of class
 #'   [`persistence`] specifying the first persistence diagram.
@@ -40,6 +51,8 @@
 #'   Wasserstein distance, it must be strictly positive.
 #' @param p A numeric value specifying the power for the Wasserstein distance.
 #'   Defaults to `1.0`.
+#' @param q A numeric value specifying the power of the internal Minkowski
+#'   metric. Defaults to `Inf`.
 #' @param validate A boolean value specifying whether to validate the input
 #'   persistence diagrams. Defaults to `TRUE`. If `FALSE`, the function will not
 #'   check if the input persistence diagrams are valid. This can be useful for
@@ -112,6 +125,7 @@ wasserstein_distance <- function(
   y,
   tol = sqrt(.Machine$double.eps),
   p = 1.0,
+  q = Inf,
   validate = TRUE,
   dimension = 0L
 ) {
@@ -127,7 +141,21 @@ wasserstein_distance <- function(
     y <- y[y[, 1] < y[, 2], , drop = FALSE]
   }
 
-  if (p > 20) {
+  if (p == Inf && q != Inf) cli::cli_abort(
+    "q-bottleneck distances (`q < Inf`) are not yet supported."
+  )
+  # TODO: Should the p threshold (currently 6) depend on the tolerance?
+  if (p >= 6 && p < Inf) {
+    cli::cli_alert_warning(
+      paste(
+        "Values `p ≥ 6` can crash or stall the Wasserstein calculation;",
+        "for `p = Inf`, use the bottleneck distance (where also `q = Inf`)."
+      ),
+      wrap = TRUE
+    )
+  }
+
+  if (p == Inf && q == Inf) {
     return(bottleneck_distance(
       x = x,
       y = y,
@@ -141,7 +169,8 @@ wasserstein_distance <- function(
     x = x,
     y = y,
     delta = tol,
-    wasserstein_power = p
+    wasserstein_power = p,
+    internal_p = q
   )
 }
 
@@ -152,6 +181,7 @@ kantorovich_distance <- function(
   y,
   tol = sqrt(.Machine$double.eps),
   p = 1.0,
+  q = Inf,
   validate = TRUE,
   dimension = 0L
 ) {
@@ -160,6 +190,7 @@ kantorovich_distance <- function(
     y = y,
     tol = tol,
     p = p,
+    q = q,
     validate = validate,
     dimension = dimension
   )
@@ -237,10 +268,12 @@ wasserstein_pairwise_distances <- function(
   x,
   tol = sqrt(.Machine$double.eps),
   p = 1.0,
+  q = Inf,
   validate = TRUE,
   dimension = 0L,
   ncores = 1L
 ) {
+
   indices <- seq_along(x)
   if (validate) {
     for (i in indices) {
@@ -250,7 +283,17 @@ wasserstein_pairwise_distances <- function(
     }
   }
 
-  if (p > 20) {
+  if (p == Inf && q != Inf) cli::cli_abort(
+    "q-bottleneck distances are not yet supported."
+  )
+  # TODO: Should the p threshold (currently 6) depend on the tolerance?
+  if (p >= 6 && p < Inf) {
+    cli::cli_alert_warning(
+      "Values `p ≥ 6` can crash or stall the Wasserstein calculation."
+    )
+  }
+
+  if (p == Inf && q == Inf) {
     return(bottleneck_pairwise_distances(
       x = x,
       tol = tol,
@@ -264,6 +307,7 @@ wasserstein_pairwise_distances <- function(
     x = x,
     delta = tol,
     wasserstein_power = p,
+    internal_p = q,
     ncores = ncores
   )
   attr(distance_matrix, "Size") <- length(x)
@@ -281,6 +325,7 @@ kantorovich_pairwise_distances <- function(
   x,
   tol = sqrt(.Machine$double.eps),
   p = 1.0,
+  q = Inf,
   validate = TRUE,
   dimension = 0L,
   ncores = 1L
@@ -289,6 +334,7 @@ kantorovich_pairwise_distances <- function(
     x = x,
     tol = tol,
     p = p,
+    q = q,
     validate = validate,
     dimension = dimension,
     ncores = ncores

inst/tinytest/test-distances.R

@@ -38,7 +38,10 @@ expect_error(
   'Wasserstein_degree was "0.000000", must be a number >= 1.0. Cannot proceed.'
 )
 expect_error(wasserstein_distance(x, y, tol = 0.0, p = 1))
-expect_equal(wasserstein_distance(x, y, p = 21), 1)
+expect_message(
+  wasserstein_distance(x, y, p = 7),
+  'crash or stall the Wasserstein calculation'
+)
 expect_equal(wasserstein_distance(x, y, p = 1), 2)
 expect_equal(round(wasserstein_distance(x, y, p = 2), digits = 6L), 1.414214)
 
@@ -79,3 +82,91 @@ expect_equal(
   kantorovich_pairwise_distances(mod_sample),
   wasserstein_pairwise_distances(mod_sample)
 )
+
+# alternative internal norms
+X <- rbind(
+  c(1, 3),
+  c(3, 5)
+)
+Y <- rbind(
+  c(3, 4)
+)
+Z <- matrix(NA_real_, nrow = 0L, ncol = 2L)
+# Manhattan
+expect_equal(wasserstein_distance(X, Y, p = 1, q = 1), 3, tol = 1e-6)
+expect_equal(wasserstein_distance(X, Y, p = 2, q = 1), sqrt(5), tol = 1e-6)
+# expect_equal(wasserstein_distance(X, Y, p = Inf, q = 1), 2, tol = 1e-6)
+expect_error(
+  wasserstein_distance(X, Y, p = Inf, q = 1),
+  "q-bottleneck distances (`q < Inf`) are not yet supported."
+)
+# Pythagorean
+expect_equal(wasserstein_distance(X, Y, p = 1, q = 2), 1+sqrt(2), tol = 1e-6)
+expect_equal(wasserstein_distance(X, Y, p = 2, q = 2), sqrt(3), tol = 1e-6)
+# expect_equal(wasserstein_distance(X, Y, p = Inf, q = 2), sqrt(2), tol = 1e-4)
+expect_error(
+  wasserstein_distance(X, Y, p = Inf, q = 2),
+  "q-bottleneck distances (`q < Inf`) are not yet supported."
+)
+# supremum (default)
+expect_equal(wasserstein_distance(X, Y, p = 1, q = Inf), 2, tol = 1e-6)
+expect_equal(wasserstein_distance(X, Y, p = 2, q = Inf), sqrt(2), tol = 1e-6)
+expect_equal(wasserstein_distance(X, Y, p = Inf, q = Inf), 1, tol = 1e-6)
+# pairwise
+XYZ <- list(X, Y, Z)
+# Manhattan
+expect_equal(
+  as.vector(wasserstein_pairwise_distances(XYZ, p = 1, q = 1)),
+  c(3, 4, 1),
+  tol = 1e-6
+)
+expect_equal(
+  as.vector(wasserstein_pairwise_distances(XYZ, p = 2, q = 1)),
+  c(sqrt(5), 2*sqrt(2), 1),
+  tol = 1e-6
+)
+# expect_equal(
+#   as.vector(wasserstein_pairwise_distances(XYZ, p = Inf, q = 1)),
+#   c(2, 2, 1),
+#   tol = 1e-1
+# )
+expect_error(
+  as.vector(wasserstein_pairwise_distances(XYZ, p = Inf, q = 1)),
+  "q-bottleneck distances (`q < Inf`) are not yet supported."
+)
+# Pythagorean
+expect_equal(
+  as.vector(wasserstein_pairwise_distances(XYZ, p = 1, q = 2)),
+  c(1+sqrt(2), 2*sqrt(2), 1/sqrt(2)),
+  tol = 1e-6
+)
+expect_equal(
+  as.vector(wasserstein_pairwise_distances(XYZ, p = 2, q = 2)),
+  c(sqrt(3), 2, 1/sqrt(2)),
+  tol = 1e-6
+)
+# expect_equal(
+#   as.vector(wasserstein_pairwise_distances(XYZ, p = Inf, q = 2)),
+#   c(sqrt(2), sqrt(2), 1/sqrt(2)),
+#   tol = 1e-1
+# )
+expect_error(
+  as.vector(wasserstein_pairwise_distances(XYZ, p = Inf, q = 2)),
+  "q-bottleneck distances (`q < Inf`) are not yet supported."
+)
+# supremum
+expect_equal(
+  as.vector(wasserstein_pairwise_distances(XYZ, p = 1, q = Inf)),
+  c(2, 2, 1/2),
+  tol = 1e-6
+)
+expect_equal(
+  as.vector(wasserstein_pairwise_distances(XYZ, p = 2, q = Inf)),
+  c(sqrt(2), sqrt(2), 1/2),
+  tol = 1e-6
+)
+expect_equal(
+  as.vector(wasserstein_pairwise_distances(XYZ, p = Inf, q = Inf)),
+  c(1, 1, 1/2),
+  tol = 1e-6
+)