From f203841141764da12c30a9a3e1da05c06b41f59e Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 25 Aug 2026 06:50:37 +0000 Subject: [PATCH 1/3] Avoid dense n x n allocations in upper_gather() upper_gather() dominated the cost of every pairwise function that calls it (distance, cosine similarity, KL distance, pair counting). Measured on a 4000 x 5 input, stats::dist() itself took 0.071s while the upper_gather() step took 8.3s and peaked at 1.1GB. The dist vector path built a dense n x n triangular matrix, the row() and col() index matrices for it, a logical mask, a transpose, and one more n x n copy inside mat_to_df(): six n x n allocations, about 60 * n^2 bytes, which is 6GB at n = 10000. It now builds the long form directly from the dist vector. 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 keep their original order and no n x n object is allocated at all. The matrix path had a second copy of the same problem. With zero.rm = FALSE it evaluated is.na(tmat) | !is.na(tmat), an all TRUE matrix of the same shape, purely to hand every cell to which(). For a sparse input that forces the sparse class to store all n^2 cells, which is heavier than a plain dense logical. The selected indices are known in closed form there, so only the surviving half is generated. The output contract is unchanged: same row order, row count, column names, column types, and NA positions. bench/micro/upper_gather_equivalence.R compares this implementation against the previous one with identical() over the six argument combinations the callers use plus degenerate, non square, multibyte and duplicated name cases: 47 of 49 identical. The two remaining cases are duplicated dimension names on the vector path, where the previous implementation raised an error from tibble inside mat_to_df() and this one returns the natural result. mat_to_df() itself is deliberately left untouched; the vector path simply no longer calls it. --- R/util.R | 128 +++++++++++---- bench/micro/upper_gather_equivalence.R | 212 +++++++++++++++++++++++++ 2 files changed, 305 insertions(+), 35 deletions(-) create mode 100644 bench/micro/upper_gather_equivalence.R 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_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") From 437429619b500e4127aace89e286a2a65e24243c Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 25 Aug 2026 06:52:53 +0000 Subject: [PATCH 2/3] Pin the upper_gather() output contract and add a micro benchmark The five callers of upper_gather() depend on details that were not covered by any test: the full n x n grid returned when na.rm and zero.rm are FALSE, the row major order, which half of the grid carries values, the column types, and the fact that dimension names are copied into the output verbatim. Adds two test_that blocks to test_util.R covering the six argument combinations the callers use, plus the degenerate cases (n = 1, n = 2, all NA) and dimension names that are multibyte, punctuation heavy, near colliding or duplicated. Adds bench/micro/upper_gather_bench.R next to the equivalence harness. It compares both paths against the previous implementation over a range of n and reports elapsed time, gc peak and row count. --- bench/micro/upper_gather_bench.R | 101 ++++++++++++++++++++++++++++++ tests/testthat/test_util.R | 104 +++++++++++++++++++++++++++++++ 2 files changed, 205 insertions(+) create mode 100644 bench/micro/upper_gather_bench.R diff --git a/bench/micro/upper_gather_bench.R b/bench/micro/upper_gather_bench.R new file mode 100644 index 00000000..fb404feb --- /dev/null +++ b/bench/micro/upper_gather_bench.R @@ -0,0 +1,101 @@ +# 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) + # columns 6 and 5 of gc() are max used Mb for Vcells and Ncells + peak <- sum(g[, 6]) + 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/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), From 0505c19efc9b81b802450f3d12bff4275955f5ae Mon Sep 17 00:00:00 2001 From: Kei Saito Date: Thu, 27 Aug 2026 22:34:33 -0700 Subject: [PATCH 3/3] Read gc()'s peak memory from the last column in the benchmark gc() inserts a "limit (Mb)" column when a memory limit is set, which shifts column 6 from "max used (Mb)" to "max used" in cells. The benchmark indexed column 6 directly, so on a machine with a limit set it reported cell counts instead of megabytes (about 9.5e6 where the value is 178 Mb). "max used (Mb)" is the last column in both shapes, so index from the end. --- bench/micro/upper_gather_bench.R | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/bench/micro/upper_gather_bench.R b/bench/micro/upper_gather_bench.R index fb404feb..e00a6486 100644 --- a/bench/micro/upper_gather_bench.R +++ b/bench/micro/upper_gather_bench.R @@ -24,8 +24,11 @@ measure <- function(fn, ...) { invisible(gc(reset = TRUE, full = TRUE)) t <- system.time(res <- fn(...))[["elapsed"]] g <- gc(full = TRUE) - # columns 6 and 5 of gc() are max used Mb for Vcells and Ncells - peak <- sum(g[, 6]) + # "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)) }