diff --git a/R/util.R b/R/util.R index a2e43ad0..b24bd354 100644 --- a/R/util.R +++ b/R/util.R @@ -219,18 +219,57 @@ upper_gather <- function(mat, names=NULL, diag=NULL, cnames = c("Var1", "Var2", } } - # create a triangler matrix to melt - # use NA_real_ for performance - trimat <- matrix(data=NA_real_, nrow=length(names), ncol=length(names)) - # fill only lower half of the matrix (transpose later to keep the order) - trimat[row(trimat)>col(trimat)] <- as.numeric(mat) - colnames(trimat) <- names - rownames(trimat) <- names + # Build the long form directly, without ever allocating an n x n matrix. + # The previous implementation created a dense n x n triangular matrix, its + # row()/col() index matrices, a logical mask, a transpose, and a data frame + # copy inside mat_to_df. That is 6 n x n allocations (about 60 * n^2 bytes), + # which is what made this the dominant cost of the pairwise functions. + # + # Output contract that is reproduced here exactly: + # - one row per cell of the full n x n grid, in row major order, so the + # first name varies slowly and the second name varies fast + # - the value is filled only on the strict upper half (i < j); the rest + # stays NA unless diag is given + # - stats::dist stores the lower half in column major order, which is the + # same sequence as the upper half in row major order, so the values of + # mat are placed in their original order + n <- length(names) + names_chr <- as.character(names) + + value <- rep(NA_real_, n * n) + if (n >= 2L) { + # i is the slow index, j the fast one, over the strict upper half + i <- rep.int(seq_len(n - 1L), times = (n - 1L):1L) + j <- sequence(nvec = (n - 1L):1L, from = 2:n) + value[(i - 1L) * n + j] <- as.numeric(mat) + } if(!is.null(diag)){ # fill diagonal elements - trimat[row(trimat)==col(trimat)] = rep(diag, length(names)) + value[seq.int(from = 1L, by = n + 1L, length.out = n)] <- rep(diag, length.out = n) + } + + var1 <- rep(names_chr, each = n) + var2 <- rep.int(names_chr, times = n) + + # mat_to_df() applies na.rm first and zero.rm after it; both keep the row + # order, so a single combined mask gives the same result. + keep <- NULL + if(na.rm){ + keep <- !is.na(value) + } + if(zero.rm){ + not_zero <- is.na(value) | value != 0 + keep <- if(is.null(keep)) not_zero else keep & not_zero } - mat_to_df(t(trimat), na.rm=na.rm, cnames=cnames, zero.rm = zero.rm) + if(!is.null(keep)){ + var1 <- var1[keep] + var2 <- var2[keep] + value <- value[keep] + } + + df <- tibble::tibble(Var1 = var1, Var2 = var2, value = value) + colnames(df) <- cnames + df } else { # diag can be NULL or FALSE if(is.null(diag)){ @@ -247,34 +286,53 @@ upper_gather <- function(mat, names=NULL, diag=NULL, cnames = c("Var1", "Var2", r_names <- seq(nrow(tmat)) } - # remove 0 if zero.rm is TRUE - ind_mat <- if(zero.rm){ - tmat != 0 - } else { - is.na(tmat) | !is.na(tmat) # Just return matrix of same shape with TRUE for all the values. - } - # preserve NA if na.rm is FALSE - if(!na.rm){ - ind_mat <- is.na(ind_mat) | ind_mat - } - # get indice of matrix - ind <- Matrix::which(ind_mat, arr.ind = TRUE) + if(zero.rm){ + # remove 0. tmat != 0 keeps a sparse matrix sparse, so this stays cheap. + ind_mat <- tmat != 0 + # preserve NA if na.rm is FALSE + if(!na.rm){ + ind_mat <- is.na(ind_mat) | ind_mat + } + # get indice of matrix + ind <- Matrix::which(ind_mat, arr.ind = TRUE) - # remove duplicated pairs - # by comparing indice - filtered <- if(diag) { - ind[ind[,2] <= ind[,1], ] - } else { - ind[ind[,2] < ind[,1], ] - } + # remove duplicated pairs + # by comparing indice + filtered <- if(diag) { + ind[ind[,2] <= ind[,1], ] + } else { + ind[ind[,2] < ind[,1], ] + } - # when there is only one index pairs, - # filtered becomes a vector, not matrix - # but matrix is expected later - # so should be converted to matrix with - # one row - if(is.vector(filtered)){ - filtered <- t(as.matrix(filtered)) + # when there is only one index pairs, + # filtered becomes a vector, not matrix + # but matrix is expected later + # so should be converted to matrix with + # one row + if(is.vector(filtered)){ + filtered <- t(as.matrix(filtered)) + } + } else { + # zero.rm FALSE means "keep the zeros too", so every cell is selected and + # na.rm cannot remove anything either. The previous implementation spelled + # that out as is.na(tmat) | !is.na(tmat), an all TRUE matrix of the same + # shape, and then called which() on it, which materialises every one of + # the n^2 cells (for a sparse input it even forces the sparse class to + # store all of them) and builds an n^2 x 2 index matrix on top. + # The selected indices are known in closed form, so build only the half + # that survives the row/column comparison below. + # Matrix::which(, arr.ind = TRUE) walks in column major order, so the + # pairs are ordered by column and then by row, which is reproduced here. + n_row <- nrow(tmat) + n_col <- ncol(tmat) + col_seq <- seq_len(n_col) + first_row <- if(diag) col_seq else col_seq + 1L + counts <- pmax(n_row - first_row + 1L, 0L) + has_any <- counts > 0L + filtered <- cbind( + row = sequence(nvec = counts[has_any], from = first_row[has_any]), + col = rep.int(col_seq[has_any], times = counts[has_any]) + ) } # this creates pairs of row and column indices diff --git a/bench/micro/upper_gather_bench.R b/bench/micro/upper_gather_bench.R new file mode 100644 index 00000000..e00a6486 --- /dev/null +++ b/bench/micro/upper_gather_bench.R @@ -0,0 +1,104 @@ +# Micro benchmark for the upper_gather() dist vector path. +# +# Compares the implementation on this branch against the verbatim previous one +# on the argument combination the pairwise distance functions actually use +# (distinct = TRUE, so na.rm = FALSE and zero.rm = FALSE), for a range of n. +# +# Usage: +# Rscript bench/micro/upper_gather_bench.R [path/to/old_impl.R] [n1,n2,...] +# +# See upper_gather_equivalence.R for how to produce the old implementation file. +# Peak memory is gc()'s "max used" across the call, in MB, which counts R heap +# only. Timings are elapsed seconds from system.time(). + +args <- commandArgs(trailingOnly = TRUE) +OLD_IMPL <- if (length(args) >= 1) args[1] else "/tmp/upper_gather_old_impl.R" +NS <- if (length(args) >= 2) as.integer(strsplit(args[2], ",")[[1]]) else c(500L, 1000L, 2000L, 4000L) + +suppressMessages(library(Matrix)) +suppressMessages(library(magrittr)) +source(OLD_IMPL) +source("R/util.R", local = FALSE) + +measure <- function(fn, ...) { + invisible(gc(reset = TRUE, full = TRUE)) + t <- system.time(res <- fn(...))[["elapsed"]] + g <- gc(full = TRUE) + # "max used (Mb)" is always the LAST column of gc(), but its INDEX is not + # fixed: gc() adds a "limit (Mb)" column when a memory limit is set (see + # mem.maxVSize / R_MAX_VSIZE), which shifts column 6 from "max used (Mb)" to + # "max used" in cells. Index from the end so both shapes report Mb. + peak <- sum(g[, ncol(g)]) + list(sec = t, peak_mb = peak, rows = nrow(res)) +} + +cat(sprintf("%6s %6s %10s %10s %10s %10s %12s %12s %7s\n", + "n", "dist_s", "old_s", "new_s", "old_peakMB", "new_peakMB", + "old_rows", "new_rows", "speedup")) + +for (n in NS) { + set.seed(1) + m <- matrix(stats::runif(n * 5), nrow = n) + rownames(m) <- paste0("r", seq_len(n)) + + invisible(gc(reset = TRUE, full = TRUE)) + t_dist <- system.time(d <- as.vector(stats::dist(m)))[["elapsed"]] + + old <- try(measure(upper_gather_old, d, rownames(m), diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE), + silent = TRUE) + if (inherits(old, "try-error")) { + cat(sprintf("%6d %6.3f OLD FAILED: %s\n", n, t_dist, + conditionMessage(attr(old, "condition")))) + old <- list(sec = NA_real_, peak_mb = NA_real_, rows = NA_integer_) + } + + new <- try(measure(upper_gather, d, rownames(m), diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE), + silent = TRUE) + if (inherits(new, "try-error")) { + cat(sprintf("%6d %6.3f NEW FAILED: %s\n", n, t_dist, + conditionMessage(attr(new, "condition")))) + next + } + + cat(sprintf("%6d %6.3f %10.3f %10.3f %10.1f %10.1f %12s %12s %6.1fx\n", + n, t_dist, old$sec, new$sec, old$peak_mb, new$peak_mb, + format(old$rows, big.mark = ","), format(new$rows, big.mark = ","), + old$sec / new$sec)) +} + +# The matrix path with zero.rm = FALSE, which is what do_cosine_sim.kv uses. +# A sparse input is the interesting case: the previous code forced the sparse +# class to store all n^2 cells before calling which() on them. +cat("\nmatrix path, sparse input, na.rm = FALSE / zero.rm = FALSE (do_cosine_sim.kv)\n") +cat(sprintf("%6s %10s %10s %10s %10s %12s %12s %7s\n", + "n", "old_s", "new_s", "old_peakMB", "new_peakMB", "old_rows", "new_rows", "speedup")) +for (n in NS) { + set.seed(2) + nnz <- max(1L, as.integer(n * n * 0.01)) + sm <- Matrix::sparseMatrix( + i = sample.int(n, nnz, replace = TRUE), + j = sample.int(n, nnz, replace = TRUE), + x = stats::runif(nnz), dims = c(n, n), + dimnames = list(paste0("r", seq_len(n)), paste0("r", seq_len(n)))) + + old <- try(measure(upper_gather_old, sm, diag = FALSE, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE), + silent = TRUE) + if (inherits(old, "try-error")) { + cat(sprintf("%6d OLD FAILED: %s\n", n, conditionMessage(attr(old, "condition")))) + old <- list(sec = NA_real_, peak_mb = NA_real_, rows = NA_integer_) + } + new <- try(measure(upper_gather, sm, diag = FALSE, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE), + silent = TRUE) + if (inherits(new, "try-error")) { + cat(sprintf("%6d NEW FAILED: %s\n", n, conditionMessage(attr(new, "condition")))) + next + } + cat(sprintf("%6d %10.3f %10.3f %10.1f %10.1f %12s %12s %6.1fx\n", + n, old$sec, new$sec, old$peak_mb, new$peak_mb, + format(old$rows, big.mark = ","), format(new$rows, big.mark = ","), + old$sec / new$sec)) +} diff --git a/bench/micro/upper_gather_equivalence.R b/bench/micro/upper_gather_equivalence.R new file mode 100644 index 00000000..bdf11f3f --- /dev/null +++ b/bench/micro/upper_gather_equivalence.R @@ -0,0 +1,212 @@ +# Equivalence harness for upper_gather(). +# +# Compares the implementation on this branch against the verbatim previous +# implementation (extracted from the base revision, see OLD_IMPL below) on the +# six argument combinations used by the callers, plus degenerate and +# multibyte-name cases. +# +# Usage: +# Rscript bench/micro/upper_gather_equivalence.R [path/to/old_impl.R] +# +# The old implementation file must define upper_gather_old() and mat_to_df(). +# Produce it with: +# git show :R/util.R > /tmp/util_old.R +# and copy mat_to_df() and upper_gather() (renamed to upper_gather_old) out of +# it verbatim. Do not retype either function from memory. + +args <- commandArgs(trailingOnly = TRUE) +OLD_IMPL <- if (length(args) >= 1) args[1] else "/tmp/upper_gather_old_impl.R" + +suppressMessages(library(Matrix)) +suppressMessages(library(magrittr)) + +source(OLD_IMPL) # upper_gather_old(), mat_to_df() +source("R/util.R", local = FALSE) # upper_gather() under test (and mat_to_df, unchanged) + +results <- list() + +# Cases where the two implementations are known and intended to differ. +# The previous vector path routed through mat_to_df(), whose +# tibble::rownames_to_column() step rejects duplicated dimension names with +# "The `.data` argument of `add_column()` must have unique names". The new +# vector path never builds that intermediate frame, so duplicated names now +# produce the natural result instead of an error. Not returning an error is the +# point of the change, so these two are listed here rather than counted as +# regressions. Every other case must be identical(). +EXPECTED_DIVERGENCE <- c("D22 vec duplicated names na=F zero=F", + "D23 vec duplicated names defaults") + +check <- function(label, ...) { + old <- try(upper_gather_old(...), silent = TRUE) + new <- try(upper_gather(...), silent = TRUE) + old_err <- inherits(old, "try-error") + new_err <- inherits(new, "try-error") + if (old_err || new_err) { + same <- old_err && new_err && + identical(conditionMessage(attr(old, "condition")), + conditionMessage(attr(new, "condition"))) + verdict <- if (same) "SAME-ERROR" else "ERROR-MISMATCH" + cat(sprintf("%-14s %-58s old_err=%s new_err=%s\n", verdict, label, old_err, new_err)) + if (!same) { + if (old_err) cat(" old: ", conditionMessage(attr(old, "condition")), "\n") + if (new_err) cat(" new: ", conditionMessage(attr(new, "condition")), "\n") + } + results[[length(results) + 1L]] <<- list(label = label, ok = same) + return(invisible(NULL)) + } + ok <- identical(old, new) + detail <- "" + if (!ok) { + ae <- all.equal(old, new) + detail <- paste(utils::head(as.character(ae), 5), collapse = " ; ") + } + cat(sprintf("%-14s %-58s rows=%s\n", if (ok) "identical" else "DIFFERENT", + label, nrow(new))) + if (!ok) cat(" all.equal: ", detail, "\n") + results[[length(results) + 1L]] <<- list(label = label, ok = ok) + invisible(NULL) +} + +# ---------------------------------------------------------------- vector path +set.seed(1) +m4 <- matrix(c(1, 2, 3, 4, 2, 1, 5, 1, 9, 3, 2, 4, 0, 1, 1, 1), nrow = 4, byrow = TRUE) +rownames(m4) <- c("a", "b", "c", "d") +d4 <- as.vector(stats::dist(m4)) + +m5 <- matrix(stats::runif(25), nrow = 5) +rownames(m5) <- c("n1", "n2", "n3", "n4", "n5") +d5 <- as.vector(stats::dist(m5)) + +# case 1: do_dist.kv_ / do_kl_dist.kv_ / do_dist.cols, distinct=TRUE diag=FALSE +check("C1 vec names diag=NULL na=F zero=F", d4, rownames(m4), diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +# case 2: same callers, diag=TRUE -> diag=0 +check("C2 vec names diag=0 na=F zero=F", d4, rownames(m4), diag = 0, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +# case 3: the same shape with the caller specific cnames and names=NULL +check("C3 vec names=NULL diag=NULL na=F zero=F", d4, NULL, diag = NULL, + cnames = c("p.x", "p.y", "value"), na.rm = FALSE, zero.rm = FALSE) +check("C3b vec names=NULL diag=0 na=F zero=F", d4, NULL, diag = 0, + cnames = c("p.x", "p.y", "value"), na.rm = FALSE, zero.rm = FALSE) +check("C3c vec defaults na=T zero=T", d4, rownames(m4), diag = 0) +check("C3d vec defaults diag=NULL", d4, rownames(m4)) +check("C3e vec n=5 na=F zero=F", d5, rownames(m5), diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("C3f vec n=5 na=T zero=F", d5, rownames(m5), diag = 0, + cnames = c("Var1", "Var2", "value"), na.rm = TRUE, zero.rm = FALSE) +check("C3g vec n=5 na=F zero=T", d5, rownames(m5), diag = 0, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = TRUE) + +# ---------------------------------------------------------------- matrix path +sim <- matrix(c(1, 0.5, 0, 0.2, 0.5, 1, 0.3, 0, 0, 0.3, 1, 0.9, 0.2, 0, 0.9, 1), nrow = 4) +dimnames(sim) <- list(c("a", "b", "c", "d"), c("a", "b", "c", "d")) + +# case 4: do_cosine_sim.kv (the zero.rm=FALSE branch) +check("C4 mat diag=FALSE na=F zero=F", sim, rownames(sim), diag = FALSE, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("C4b mat diag=NULL na=F zero=F", sim, rownames(sim), diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("C4c mat diag=TRUE na=F zero=F", sim, rownames(sim), diag = TRUE, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) + +sim_na <- sim +sim_na[3, 1] <- NA; sim_na[1, 3] <- NA; sim_na[2, 4] <- 0; sim_na[4, 2] <- 0 +check("C4d mat with NA and 0, na=F zero=F", sim_na, rownames(sim_na), diag = FALSE, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("C4e mat with NA and 0, na=T zero=F", sim_na, rownames(sim_na), diag = FALSE, + cnames = c("Var1", "Var2", "value"), na.rm = TRUE, zero.rm = FALSE) + +# case 5 / 6: pair_count_ (sparse co-occurrence, defaults) +cm <- Matrix::Matrix(c(3, 1, 0, 2, 1, 4, 0, 0, 0, 0, 2, 1, 2, 0, 1, 5), nrow = 4, sparse = TRUE) +dimnames(cm) <- list(c("a", "b", "c", "d"), c("a", "b", "c", "d")) +check("C5 dgCMatrix diag=FALSE defaults", cm, diag = FALSE, cnames = c("v.x", "v.y", "value")) +check("C6 dgCMatrix diag=TRUE defaults", cm, diag = TRUE, cnames = c("v.x", "v.y", "value")) +check("C6b dgCMatrix diag=FALSE na=F zero=F", cm, diag = FALSE, + cnames = c("v.x", "v.y", "value"), na.rm = FALSE, zero.rm = FALSE) +check("C6c dgCMatrix diag=TRUE na=F zero=F", cm, diag = TRUE, + cnames = c("v.x", "v.y", "value"), na.rm = FALSE, zero.rm = FALSE) +check("C6d dgCMatrix no dimnames na=F zero=F", Matrix::Matrix(matrix(1:16, 4), sparse = TRUE), + diag = FALSE, na.rm = FALSE, zero.rm = FALSE) + +# ------------------------------------------------------- degenerate and names +# n < 2 +check("D1 vec length 0 (n=1) na=F zero=F", numeric(0), NULL, diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("D2 vec length 0 (n=1) diag=0", numeric(0), NULL, diag = 0, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("D3 vec length 0 defaults", numeric(0), NULL) +check("D4 vec length 1 (n=2) na=F zero=F", 2.5, c("x", "y"), diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("D5 vec length 1 (n=2) defaults", 2.5, c("x", "y")) +# all values NA +check("D6 vec all NA na=F zero=F", rep(NA_real_, 6), letters[1:4], diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("D7 vec all NA defaults", rep(NA_real_, 6), letters[1:4]) +# all values 0 +check("D8 vec all zero na=F zero=T", rep(0, 6), letters[1:4], diag = 0, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = TRUE) +# name mismatch must still stop() with the same message +check("D9 vec wrong name count", d4, c("a", "b"), diag = NULL) +# 1x1 matrix path +one <- matrix(1, 1, 1, dimnames = list("a", "a")) +check("D10 mat 1x1 diag=FALSE defaults", one, diag = FALSE) +check("D11 mat 1x1 diag=TRUE defaults", one, diag = TRUE) +check("D12 mat 1x1 diag=FALSE na=F zero=F", one, diag = FALSE, na.rm = FALSE, zero.rm = FALSE) +check("D13 mat 1x1 diag=TRUE na=F zero=F", one, diag = TRUE, na.rm = FALSE, zero.rm = FALSE) +# single surviving pair (the is.vector(filtered) path) +two <- matrix(c(0, 7, 7, 0), 2, dimnames = list(c("a", "b"), c("a", "b"))) +check("D14 mat 2x2 single pair defaults", two, diag = FALSE) +check("D15 mat 2x2 single pair na=F zero=F", two, diag = FALSE, na.rm = FALSE, zero.rm = FALSE) +# all-zero matrix, zero.rm=TRUE -> empty result +zmat <- matrix(0, 3, 3, dimnames = list(letters[1:3], letters[1:3])) +check("D16 mat all zero defaults (empty)", zmat, diag = FALSE) + +# multibyte and near-colliding names (project convention stress name). +# Non ASCII here is test data only; the package sources stay ASCII. +stress <- c("航空 会社 !\"#$%&'()*+, -./:;<=>?@[]^_`{|}~ 表", + "col", "col ", "col\n") +check("D17 vec multibyte names na=F zero=F", d4, stress, diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("D18 vec multibyte names diag=0", d4, stress, diag = 0, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("D19 vec multibyte names defaults", d4, stress, diag = 0) +sim_stress <- sim +dimnames(sim_stress) <- list(stress, stress) +check("D20 mat multibyte names na=F zero=F", sim_stress, stress, diag = FALSE, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("D21 mat multibyte names defaults", sim_stress, stress, diag = FALSE) +# duplicated names +dup <- c("a", "a", "b", "c") +check("D22 vec duplicated names na=F zero=F", d4, dup, diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +check("D23 vec duplicated names defaults", d4, dup, diag = 0) +sim_dup <- sim +dimnames(sim_dup) <- list(dup, dup) +check("D24 mat duplicated names na=F zero=F", sim_dup, dup, diag = FALSE, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +# numeric-looking names +check("D25 vec numeric names na=F zero=F", d4, c(10, 20, 30, 40), diag = NULL, + cnames = c("Var1", "Var2", "value"), na.rm = FALSE, zero.rm = FALSE) +# non square matrix on the matrix path +rect <- matrix(1:12, nrow = 4, dimnames = list(letters[1:4], c("p", "q", "r"))) +check("D26 mat non square defaults", rect, diag = FALSE) +check("D27 mat non square na=F zero=F", rect, diag = FALSE, na.rm = FALSE, zero.rm = FALSE) +check("D28 mat non square diag=TRUE na=F zero=F", rect, diag = TRUE, na.rm = FALSE, zero.rm = FALSE) +rect2 <- matrix(1:12, nrow = 3, dimnames = list(letters[1:3], c("p", "q", "r", "s"))) +check("D29 mat wide na=F zero=F", rect2, diag = FALSE, na.rm = FALSE, zero.rm = FALSE) +check("D30 mat wide diag=TRUE na=F zero=F", rect2, diag = TRUE, na.rm = FALSE, zero.rm = FALSE) + +ok_n <- sum(vapply(results, function(r) isTRUE(r$ok), logical(1))) +cat(sprintf("\n%d / %d cases match\n", ok_n, length(results))) +diverged <- Filter(function(r) !isTRUE(r$ok), results) +expected <- Filter(function(r) r$label %in% EXPECTED_DIVERGENCE, diverged) +unexpected <- Filter(function(r) !(r$label %in% EXPECTED_DIVERGENCE), diverged) +if (length(expected)) { + cat("expected divergence (old implementation errored, new one does not): ", + paste(vapply(expected, function(r) r$label, character(1)), collapse = ", "), "\n") +} +if (length(unexpected)) { + cat("FAILED: ", paste(vapply(unexpected, function(r) r$label, character(1)), collapse = ", "), "\n") + quit(status = 1) +} +cat("no unexpected divergence\n") diff --git a/tests/testthat/test_util.R b/tests/testthat/test_util.R index 8eb9fa11..0996cf49 100644 --- a/tests/testthat/test_util.R +++ b/tests/testthat/test_util.R @@ -349,6 +349,110 @@ test_that("test upper_gather with vector diag true", { expect_equal(nrow(result), 10) }) +test_that("upper_gather output contract for the pairwise callers", { + # do_dist.kv_, do_kl_dist.kv_ and do_dist.cols all call upper_gather with + # na.rm = FALSE and zero.rm = FALSE, which means nothing is dropped and the + # result is the full n x n grid in row major order, with values only on the + # strict upper half. Downstream display and do_cmdscale_ depend on that + # shape, so pin it here. + names <- c("a", "b", "c", "d") + vec <- seq(6) + + res <- upper_gather(vec, names, diag = NULL, + cnames = c("Var1", "Var2", "value"), + na.rm = FALSE, zero.rm = FALSE) + expect_equal(nrow(res), 16) + expect_equal(colnames(res), c("Var1", "Var2", "value")) + expect_true(is.character(res[[1]])) + expect_true(is.character(res[[2]])) + expect_true(is.numeric(res[[3]])) + # row major: the first name varies slowly, the second one fast + expect_equal(res[[1]], rep(names, each = 4)) + expect_equal(res[[2]], rep(names, times = 4)) + # values sit on the strict upper half, in the order stats::dist stores them + expect_equal(res[[3]][c(2, 3, 4, 7, 8, 12)], as.numeric(vec)) + expect_equal(which(is.na(res[[3]])), c(1L, 5L, 6L, 9L, 10L, 11L, 13L, 14L, 15L, 16L)) + + # diag = 0 is what those callers pass when the diagonal is requested + res_diag <- upper_gather(vec, names, diag = 0, + cnames = c("Var1", "Var2", "value"), + na.rm = FALSE, zero.rm = FALSE) + expect_equal(nrow(res_diag), 16) + expect_equal(res_diag[[3]][c(1, 6, 11, 16)], rep(0, 4)) + + # do_cosine_sim.kv uses the matrix path with the same na.rm / zero.rm + sim <- matrix(c(1, 0.5, 0, 0.2, 0.5, 1, 0.3, 0, 0, 0.3, 1, 0.9, 0.2, 0, 0.9, 1), nrow = 4) + dimnames(sim) <- list(names, names) + res_mat <- upper_gather(sim, names, diag = FALSE, + cnames = c("Var1", "Var2", "value"), + na.rm = FALSE, zero.rm = FALSE) + # zeros are kept because zero.rm is FALSE + expect_equal(nrow(res_mat), 6) + expect_equal(res_mat[[1]], c("a", "a", "a", "b", "b", "c")) + expect_equal(res_mat[[2]], c("b", "c", "d", "c", "d", "d")) + expect_equal(res_mat[[3]], c(0.5, 0, 0.2, 0.3, 0, 0.9)) + + # pair_count_ uses the matrix path with the defaults, where zeros and NAs go + cm <- Matrix::Matrix(c(3, 1, 0, 2, 1, 4, 0, 0, 0, 0, 2, 1, 2, 0, 1, 5), + nrow = 4, sparse = TRUE) + dimnames(cm) <- list(names, names) + res_pc <- upper_gather(cm, diag = FALSE, cnames = c("v.x", "v.y", "value")) + expect_equal(nrow(res_pc), 3) + expect_equal(res_pc[[1]], c("a", "a", "c")) + expect_equal(res_pc[[2]], c("b", "d", "d")) + expect_equal(res_pc[[3]], c(1, 2, 1)) + res_pc_diag <- upper_gather(cm, diag = TRUE, cnames = c("v.x", "v.y", "value")) + expect_equal(nrow(res_pc_diag), 7) +}) + +test_that("upper_gather degenerate and awkward dimension names", { + names <- c("a", "b", "c", "d") + vec <- seq(6) + + # n = 1 (empty dist vector) must not fall over + res1 <- upper_gather(numeric(0), NULL, diag = NULL, + cnames = c("Var1", "Var2", "value"), + na.rm = FALSE, zero.rm = FALSE) + expect_equal(nrow(res1), 1) + expect_true(is.na(res1[[3]][1])) + expect_equal(nrow(upper_gather(numeric(0), NULL)), 0) + + # n = 2 + res2 <- upper_gather(2.5, c("x", "y"), diag = NULL, + cnames = c("Var1", "Var2", "value"), + na.rm = FALSE, zero.rm = FALSE) + expect_equal(nrow(res2), 4) + expect_equal(res2[[3]][2], 2.5) + + # every value NA + expect_equal(nrow(upper_gather(rep(NA_real_, 6), names, diag = NULL, + cnames = c("Var1", "Var2", "value"), + na.rm = FALSE, zero.rm = FALSE)), 16) + expect_equal(nrow(upper_gather(rep(NA_real_, 6), names)), 0) + + # a name count that does not match the vector length is still an error + expect_error(upper_gather(vec, c("a", "b"))) + + # dimension names are used as output values, so they must survive verbatim, + # including multibyte characters, punctuation and near collisions + stress <- c("航空 会社 !\"#$%&'()*+, -./:;<=>?@[]^_`{|}~ 表", + "col", "col ", "col\n") + res_s <- upper_gather(vec, stress, diag = NULL, + cnames = c("Var1", "Var2", "value"), + na.rm = FALSE, zero.rm = FALSE) + expect_equal(unique(res_s[[1]]), stress) + expect_equal(res_s[[2]][1:4], stress) + + # duplicated names must not be mangled or rejected + dup <- c("a", "a", "b", "c") + res_d <- upper_gather(vec, dup, diag = NULL, + cnames = c("Var1", "Var2", "value"), + na.rm = FALSE, zero.rm = FALSE) + expect_equal(nrow(res_d), 16) + expect_equal(res_d[[1]], rep(dup, each = 4)) + expect_equal(res_d[[2]], rep(dup, times = 4)) +}) + test_that("sparse_cast", { test_df <- data.frame( row = rep(paste("row", 6-seq(5)), each=4),