Fix Spearman correlation of CSR and symmetric input

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: dynutils
 Type: Package
 Title: Common Functionality for the 'dynverse' Packages
-Version: 1.0.12
+Version: 1.0.13
 Authors@R: c(
     person(
       "Robrecht", 

NAMESPACE

@@ -75,6 +75,7 @@ importFrom(magrittr,set_rownames)
 importFrom(methods,as)
 importFrom(methods,formalArgs)
 importFrom(methods,is)
+importFrom(methods,new)
 importFrom(proxyC,dist)
 importFrom(proxyC,simil)
 importFrom(purrr,"%>%")

NEWS.md

@@ -1,3 +1,7 @@
+# dynutils 1.0.13
+
+* BUG FIX `calculate_similarity()`, `calculate_distance()`: Fix Spearman correlations of RsparseMatrix, TsparseMatrix and symmetric inputs (PR #37).
+
 # dynutils 1.0.12
 
 * DOCUMENTATION: Fix roxygen docs (PR #36).

R/calculate_distance.R

@@ -107,37 +107,38 @@ calculate_similarity <- function(
 #' @export
 list_similarity_methods <- function() eval(formals(calculate_similarity)$method)
 
-#' @importFrom Matrix t
-#' @importFrom methods as
 .process_input_matrices <- function(x, y, margin) {
-  if (is.data.frame(x)) x <- as.matrix(x)
-  if (is.matrix(x)) x <- as(x, "CsparseMatrix")
-
   if (!is.null(y)) {
-    if (is.data.frame(y)) y <- as.matrix(y)
-    if (is.matrix(y)) y <- as(y, "CsparseMatrix")
-  }
-
-  assert_that(
-    is_sparse(x),
-    is.null(y) || is_sparse(y)
-  )
-
-  if (margin == 1) {
-    if (!is.null(y)) {
+    if (margin == 1) {
       assert_that(ncol(x) == ncol(y))
-      y <- Matrix::t(y)
-    }
-    x <- Matrix::t(x)
-  } else {
-    if (!is.null(y)) {
+    } else {
       assert_that(nrow(x) == nrow(y))
     }
   }
 
+  x <- .as_column_vectors(x, margin = margin)
+  if (!is.null(y)) {
+    y <- .as_column_vectors(y, margin = margin)
+  }
+
   list(x = x, y = y)
 }
 
+# convert to a dgCMatrix whose columns are the vectors to compare
+#' @importFrom Matrix t
+#' @importFrom methods as is new
+.as_column_vectors <- function(x, margin) {
+  if (is.data.frame(x)) x <- as.matrix(x)
+
+  # a dgRMatrix holds the same slots as the dgCMatrix of its transpose
+  if (margin == 1 && is(x, "dgRMatrix")) {
+    return(new("dgCMatrix", i = x@j, p = x@p, x = x@x, Dim = rev(x@Dim), Dimnames = rev(x@Dimnames)))
+  }
+
+  x <- as(as(as(x, "CsparseMatrix"), "generalMatrix"), "dMatrix")
+  if (margin == 1) Matrix::t(x) else x
+}
+
 .process_pairwise_matrix <- function(mat, x, y) {
   assert_that(
     ncol(x) == nrow(mat),

R/spearman_rank_sparse.R

@@ -1,4 +1,5 @@
 spearman_rank_sparse <- function(x) {
+  x <- .as_column_vectors(x, margin = 2)
   x@x <- spearman_rank_sparse_rcpp(x@x, x@p, nrow(x))
   x
 }

tests/testthat/test-calculate_distance.R

@@ -119,3 +119,76 @@ test_that("calculate_distance works on matrices, data frames and sparse matrices
     expect_true(sum(abs(na.omit(out3-out))) < 1e-10)
   }
 })
+
+
+# the same data in every storage format that calculate_distance() and calculate_similarity() accept
+as_formats <- function(x) {
+  list(
+    dense = x,
+    data_frame = as.data.frame(x),
+    csc = as(x, "CsparseMatrix"),
+    csr = as(x, "RsparseMatrix"),
+    triplet = as(x, "TsparseMatrix"),
+    symmetric = Matrix::forceSymmetric(as(x, "CsparseMatrix"))
+  )
+}
+
+set.seed(1)
+x_sym <- as.matrix(Matrix::forceSymmetric(Matrix::rsparsematrix(8, 8, .5)))
+# forceSymmetric() symmetrises dimnames, so use the same row and column names
+dimnames(x_sym) <- list(paste0("c", 1:8), paste0("c", 1:8))
+y_any <- as.matrix(Matrix::rsparsematrix(8, 8, .5))
+dimnames(y_any) <- list(paste0("s", 1:8), paste0("d", 1:8))
+
+test_that("calculate_similarity gives the same result for every storage format and margin", {
+  for (margin in c(1, 2)) {
+    for (method in list_similarity_methods()) {
+      expected <- calculate_similarity(x_sym, y_any, margin = margin, method = method)
+
+      for (fmt in names(as_formats(x_sym))) {
+        x <- as_formats(x_sym)[[fmt]]
+        y <- as_formats(y_any)[[if (fmt == "symmetric") "csr" else fmt]]
+        observed <- calculate_similarity(x, y, margin = margin, method = method)
+        expect_equal(observed, expected, tolerance = 1e-10, info = paste(method, fmt, margin))
+      }
+    }
+  }
+})
+
+test_that("calculate_similarity computes the correct spearman correlation for every storage format and margin", {
+  for (margin in c(1, 2)) {
+    xd <- if (margin == 1) t(x_sym) else x_sym
+    yd <- if (margin == 1) t(y_any) else y_any
+    gold <- suppressWarnings(cor(xd, yd, method = "spearman"))
+    gold[is.na(gold)] <- 0
+
+    for (fmt in names(as_formats(x_sym))) {
+      x <- as_formats(x_sym)[[fmt]]
+      y <- as_formats(y_any)[[if (fmt == "symmetric") "csr" else fmt]]
+      observed <- calculate_similarity(x, y, margin = margin, method = "spearman")
+      expect_equal(observed, gold, tolerance = 1e-10, info = paste(fmt, margin))
+    }
+  }
+})
+
+test_that("calculate_distance gives the same result for every storage format and margin", {
+  for (margin in c(1, 2)) {
+    for (method in methods) {
+      expected <- calculate_distance(x_sym, y_any, margin = margin, method = method)
+
+      for (fmt in names(as_formats(x_sym))) {
+        x <- as_formats(x_sym)[[fmt]]
+        y <- as_formats(y_any)[[if (fmt == "symmetric") "csr" else fmt]]
+        observed <- calculate_distance(x, y, margin = margin, method = method)
+        expect_equal(observed, expected, tolerance = 1e-10, info = paste(method, fmt, margin))
+      }
+    }
+  }
+})
+
+test_that(".as_column_vectors reuses the slots of a CSR matrix for margin 1", {
+  x <- as(y_any, "RsparseMatrix")
+  expect_is(x, "dgRMatrix")
+  expect_identical(.as_column_vectors(x, margin = 1), Matrix::t(as(x, "CsparseMatrix")))
+  expect_identical(.as_column_vectors(x, margin = 2), as(x, "CsparseMatrix"))
+})

tests/testthat/test-spearman_rank_sparse.R

@@ -112,3 +112,19 @@ test_that("spearman_rank_sparse works with only positive integers", {
 
 
 
+test_that("spearman_rank_sparse ranks columns regardless of the storage format", {
+  x <- as.matrix(Matrix::forceSymmetric(Matrix::rsparsematrix(50, 50, .3)))
+  gold <- cor(x, method = "spearman") %>% set_nas_to_two()
+
+  formats <- list(
+    csc = as(x, "CsparseMatrix"),
+    csr = as(x, "RsparseMatrix"),
+    triplet = as(x, "TsparseMatrix"),
+    symmetric = Matrix::forceSymmetric(as(x, "CsparseMatrix"))
+  )
+  for (fmt in names(formats)) {
+    ranked <- spearman_rank_sparse(formats[[fmt]])
+    rcpp_imp <- cor(as.matrix(ranked), method = "pearson") %>% set_nas_to_two()
+    expect_true(sum(abs(gold - rcpp_imp)) < 1e-8, info = fmt)
+  }
+})