diff --git a/NEWS.md b/NEWS.md index 8e58f5e..589a0de 100644 --- a/NEWS.md +++ b/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 diff --git a/R/cpp11.R b/R/cpp11.R index 8884e7e..7104c58 100644 --- a/R/cpp11.R +++ b/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) } diff --git a/R/distances.R b/R/distances.R index 06c6f75..b67c97a 100644 --- a/R/distances.R +++ b/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 diff --git a/inst/tinytest/test-distances.R b/inst/tinytest/test-distances.R index aca386b..46c4d99 100644 --- a/inst/tinytest/test-distances.R +++ b/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 +) diff --git a/man/distances.Rd b/man/distances.Rd index bfe5761..bf96d66 100644 --- a/man/distances.Rd +++ b/man/distances.Rd @@ -20,6 +20,7 @@ wasserstein_distance( y, tol = sqrt(.Machine$double.eps), p = 1, + q = Inf, validate = TRUE, dimension = 0L ) @@ -29,6 +30,7 @@ kantorovich_distance( y, tol = sqrt(.Machine$double.eps), p = 1, + q = Inf, validate = TRUE, dimension = 0L ) @@ -57,6 +59,9 @@ are objects of class \code{\link{persistence}}.} \item{p}{A numeric value specifying the power for the Wasserstein distance. Defaults to \code{1.0}.} + +\item{q}{A numeric value specifying the power of the internal Minkowski +metric. Defaults to \code{Inf}.} } \value{ A numeric value storing either the Bottleneck or the Wasserstein @@ -72,26 +77,31 @@ contain the birth values and the death values, respectively. \details{ 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 \emph{\eqn{p}-Wasserstein distance} +the diagonal with infinite multiplicity. The \emph{\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 -\emph{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}-\emph{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). } \examples{ bottleneck_distance( @@ -114,6 +124,11 @@ wasserstein_distance( persistence_sample[[2]] ) +} +\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 } \seealso{ \href{https://github.com/anigmetov/hera}{the Hera C++ library} diff --git a/man/pairwise-distances.Rd b/man/pairwise-distances.Rd index 15b8e5f..e2f6fb6 100644 --- a/man/pairwise-distances.Rd +++ b/man/pairwise-distances.Rd @@ -19,6 +19,7 @@ wasserstein_pairwise_distances( x, tol = sqrt(.Machine$double.eps), p = 1, + q = Inf, validate = TRUE, dimension = 0L, ncores = 1L @@ -28,6 +29,7 @@ kantorovich_pairwise_distances( x, tol = sqrt(.Machine$double.eps), p = 1, + q = Inf, validate = TRUE, dimension = 0L, ncores = 1L @@ -57,6 +59,9 @@ parallel computation. Defaults to \code{1L}.} \item{p}{A numeric value specifying the power for the Wasserstein distance. Defaults to \code{1.0}.} + +\item{q}{A numeric value specifying the power of the internal Minkowski +metric. Defaults to \code{Inf}.} } \value{ An object of class 'dist' containing the pairwise distance matrix diff --git a/src/cpp11.cpp b/src/cpp11.cpp index 828187a..ff50b12 100644 --- a/src/cpp11.cpp +++ b/src/cpp11.cpp @@ -20,17 +20,17 @@ extern "C" SEXP _phutil_bottleneckPairwiseDistances(SEXP x, SEXP delta, SEXP nco END_CPP11 } // wasserstein.cpp -double wassersteinDistance(const cpp11::doubles_matrix<>& x, const cpp11::doubles_matrix<>& y, const double delta, const double wasserstein_power); -extern "C" SEXP _phutil_wassersteinDistance(SEXP x, SEXP y, SEXP delta, SEXP wasserstein_power) { +double wassersteinDistance(const cpp11::doubles_matrix<>& x, const cpp11::doubles_matrix<>& y, const double delta, const double wasserstein_power, const double internal_p); +extern "C" SEXP _phutil_wassersteinDistance(SEXP x, SEXP y, SEXP delta, SEXP wasserstein_power, SEXP internal_p) { BEGIN_CPP11 - return cpp11::as_sexp(wassersteinDistance(cpp11::as_cpp&>>(x), cpp11::as_cpp&>>(y), cpp11::as_cpp>(delta), cpp11::as_cpp>(wasserstein_power))); + return cpp11::as_sexp(wassersteinDistance(cpp11::as_cpp&>>(x), cpp11::as_cpp&>>(y), cpp11::as_cpp>(delta), cpp11::as_cpp>(wasserstein_power), cpp11::as_cpp>(internal_p))); END_CPP11 } // wasserstein.cpp -cpp11::doubles wassersteinPairwiseDistances(const cpp11::list& x, const double delta, const double wasserstein_power, const unsigned int ncores); -extern "C" SEXP _phutil_wassersteinPairwiseDistances(SEXP x, SEXP delta, SEXP wasserstein_power, SEXP ncores) { +cpp11::doubles wassersteinPairwiseDistances(const cpp11::list& x, const double delta, const double wasserstein_power, const double internal_p, const unsigned int ncores); +extern "C" SEXP _phutil_wassersteinPairwiseDistances(SEXP x, SEXP delta, SEXP wasserstein_power, SEXP internal_p, SEXP ncores) { BEGIN_CPP11 - return cpp11::as_sexp(wassersteinPairwiseDistances(cpp11::as_cpp>(x), cpp11::as_cpp>(delta), cpp11::as_cpp>(wasserstein_power), cpp11::as_cpp>(ncores))); + return cpp11::as_sexp(wassersteinPairwiseDistances(cpp11::as_cpp>(x), cpp11::as_cpp>(delta), cpp11::as_cpp>(wasserstein_power), cpp11::as_cpp>(internal_p), cpp11::as_cpp>(ncores))); END_CPP11 } @@ -38,8 +38,8 @@ extern "C" { static const R_CallMethodDef CallEntries[] = { {"_phutil_bottleneckDistance", (DL_FUNC) &_phutil_bottleneckDistance, 3}, {"_phutil_bottleneckPairwiseDistances", (DL_FUNC) &_phutil_bottleneckPairwiseDistances, 3}, - {"_phutil_wassersteinDistance", (DL_FUNC) &_phutil_wassersteinDistance, 4}, - {"_phutil_wassersteinPairwiseDistances", (DL_FUNC) &_phutil_wassersteinPairwiseDistances, 4}, + {"_phutil_wassersteinDistance", (DL_FUNC) &_phutil_wassersteinDistance, 5}, + {"_phutil_wassersteinPairwiseDistances", (DL_FUNC) &_phutil_wassersteinPairwiseDistances, 5}, {NULL, NULL, 0} }; } diff --git a/src/wasserstein.cpp b/src/wasserstein.cpp index 11738c1..220cfa6 100644 --- a/src/wasserstein.cpp +++ b/src/wasserstein.cpp @@ -63,24 +63,28 @@ double wassersteinDist(PairVector& diagramA, double wassersteinDistance(const cpp11::doubles_matrix<>& x, const cpp11::doubles_matrix<>& y, const double delta = 0.01, - const double wasserstein_power = 1.0) + const double wasserstein_power = 1.0, + const double internal_p = std::numeric_limits::infinity()) { PairVector diagramA, diagramB; parseMatrix(x, diagramA); parseMatrix(y, diagramB); - return wassersteinDist(diagramA, diagramB, wasserstein_power, delta); + double hera_internal_p = std::isinf(internal_p) ? hera::get_infinity() : internal_p; + return wassersteinDist(diagramA, diagramB, wasserstein_power, delta, hera_internal_p); } [[cpp11::register]] cpp11::doubles wassersteinPairwiseDistances(const cpp11::list& x, const double delta = 0.01, const double wasserstein_power = 1.0, + const double internal_p = std::numeric_limits::infinity(), const unsigned int ncores = 1) { unsigned int N = x.size(); unsigned int K = N * (N - 1) / 2; cpp11::writable::doubles result(K); std::vector pairs(N); + double hera_internal_p = std::isinf(internal_p) ? hera::get_infinity() : internal_p; for (int n = 0;n < N;++n) { @@ -95,7 +99,7 @@ cpp11::doubles wassersteinPairwiseDistances(const cpp11::list& x, { unsigned int i = N - 2 - std::floor(std::sqrt(-8 * k + 4 * N * (N - 1) - 7) / 2.0 - 0.5); unsigned int j = k + i + 1 - N * (N - 1) / 2 + (N - i) * ((N - i) - 1) / 2; - result[k] = wassersteinDist(pairs[i], pairs[j], wasserstein_power, delta); + result[k] = wassersteinDist(pairs[i], pairs[j], wasserstein_power, delta, hera_internal_p); } return result; diff --git a/vignettes/validation-benchmark.qmd b/vignettes/validation-benchmark.qmd index 5d61357..4051d07 100644 --- a/vignettes/validation-benchmark.qmd +++ b/vignettes/validation-benchmark.qmd @@ -170,17 +170,19 @@ $X$ and $Y$ for several choices of half-plane and "matched space" metrics. The results make intuitive sense; for example, the values change monotonically along each row and column. -Let us now validate the bottom row---using the $L^\infty$ distance on the -half-plane, giving the popular _bottleneck distance_---using both Hera, as -exposed through {phutil}, and Dionysus, as exposed through {TDA}: +Let us now validate them using Hera, as exposed through {phutil}---at least, those that are currently supported: ```{r validate small PDs with Hera} -wasserstein_distance(X, Y, p = 1) -wasserstein_distance(X, Y, p = 2) -bottleneck_distance(X, Y) +ps <- qs <- c(1, 2, Inf) +sapply(ps, \(p) sapply(qs, \(q) { + as.numeric(try(wasserstein_distance(X, Y, p = p, q = q))) +})) ``` -In order to compute distances with {TDA}, we must restructure the PDs to include +We also validate these calculations using Dionysus, as exposed through {TDA}; +though only the bottom row---using the $L^\infty$ distance on the +half-plane---can be checked for lack of user control of the internal metric. +In order to compute distances, we must restructure the PDs to include a `"dimension"` column. Note also that `TDA::wasserstein()` does not take the $1/p$th power after computing the sum of $p$th powers; we do this manually to get comparable