Add trout distance

- Add trout distance function
- Add trout distance as option
  to weight schemes

Author: michaelweylandt

R/weights.R

@@ -181,7 +181,7 @@ dense_rbf_kernel_weights <- function(phi = "auto",
 
   tryCatch(dist.method <- match.arg(dist.method),
            error = function(e){
-             crv_error("Unsupported choice of ", sQuote("weight.dist;"),
+             crv_error("Unsupported choice of ", sQuote("dist.method;"),
                        " see the ", sQuote("method"), " argument of ",
                        sQuote("stats::dist"), " for supported distances.")
            })
@@ -192,6 +192,12 @@ dense_rbf_kernel_weights <- function(phi = "auto",
               " argument of ", sQuote("stats::dist"), " for details.")
   }
 
+  if(dist.method == "trout"){
+    dist_f <- trout_dist
+  } else {
+    dist_f <- function(X) dist(X, method = dist.method, p = p)
+  }
+
   function(X){
     user_phi <- (phi != "auto")
 
@@ -200,7 +206,7 @@ dense_rbf_kernel_weights <- function(phi = "auto",
       ##         necessary...
       phi_range <- 10^(seq(-10, 10, length.out = 21))
       weight_vars <- vapply(phi_range,
-                            function(phi) var(exp((-1) * phi * (dist(X, method = dist.method, p = p)[TRUE])^2)),
+                            function(phi) var(exp((-1) * phi * (dist_f(X)[TRUE])^2)),
                             numeric(1))
 
       phi <- phi_range[which.max(weight_vars)]
@@ -214,7 +220,7 @@ dense_rbf_kernel_weights <- function(phi = "auto",
       crv_error(sQuote("phi"), " must be positive.")
     }
 
-    dist_mat <- as.matrix(dist(X, method = dist.method, p = p))
+    dist_mat <- as.matrix(dist_f(X))
     dist_mat <- exp(-1 * phi * dist_mat^2)
 
     check_weight_matrix(dist_mat)

src/utils.cpp

@@ -299,3 +299,26 @@ Eigen::MatrixXcd align_phase(const Eigen::MatrixXcd& U,
 
   return V;
 }
+
+// [[Rcpp::export(rng = false)]]
+Rcpp::NumericVector trout_dist(const Eigen::MatrixXcd& X){
+  Eigen::Index n = X.rows();
+  Rcpp::NumericVector distances(n * (n - 1) / 2);
+
+  Eigen::Index ix = 0;
+  for(Eigen::Index i = 0; i < n; i++){
+    for(Eigen::Index j = 0; j < i; j++){
+
+      distances(ix) = std::sqrt((X.row(i) - align_phase_v(X.row(j), X.row(i))).squaredNorm());
+      ix++; // FIXME - Would be better to calculate this as a function of (i, j) explicitly
+    }
+  }
+
+  distances.attr("Size") = n;
+  distances.attr("Diag") = false;
+  distances.attr("Upper") = false;
+  distances.attr("method") = "trout";
+  distances.attr("call") = R_NilValue;
+  distances.attr("class") = "dist";
+
+  return distances;

tests/testthat/test_trout_distance.R

@@ -0,0 +1,30 @@
+context("Test trout_distance")
+
+test_that("Trout distance doesn't depend on signs", {
+  trout_dist <- clustRviz:::trout_dist
+
+  X <- matrix(c(1, 1, -1, -1), ncol = 2)
+  expect_equal(as.vector(trout_dist(X)), 0)
+
+  X <- matrix(c(1, -1, 1, -1), ncol = 2)
+  expect_equal(as.vector(trout_dist(X)), 0)
+
+  X <- matrix(c(1 + 1i, -1 -1i, 1 + 1i, -1 -1i), ncol = 2)
+  expect_equal(as.vector(trout_dist(X)), 0)
+})
+
+test_that("Trout distance minimizes distance", {
+  trout_dist <- clustRviz:::trout_dist
+  set.seed(125)
+
+  X1 <- rnorm(25) + (0 + 1i) * rnorm(25)
+  X2 <- rnorm(25) + (0 + 1i) * rnorm(25)
+  X <- rbind(X1, X2)
+
+  theta_grid <- seq(0, 2 * pi, length.out = 501)
+  d <- Vectorize(function(theta) sum(Mod(X1 - exp((0 + 1i) * theta) * X2)^2))
+  min_d <- sqrt(min(d(theta_grid)))
+
+  td <- as.vector(trout_dist(X))
+  expect_lte(td, min_d)
+})