From 185cb911b95a901d79034b1dc763caefe3cd7fc6 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 31 Mar 2026 17:57:43 +0100 Subject: [PATCH 01/10] Rename cosine-specific functions to generic connectivity names - prepare_cosine_matrix -> prepare_similarity_matrix (old name kept as alias) - cosine_heatmap -> connectivity_heatmap (old name kept as alias) - Add distfun parameter to connectivity_heatmap for pluggable distance --- NAMESPACE | 4 ++++ R/cosine.R | 18 +++++++++++++----- R/heatmap.R | 17 ++++++++++------- ..._matrix.Rd => prepare_similarity_matrix.Rd} | 16 ++++++++++------ 4 files changed, 37 insertions(+), 18 deletions(-) rename man/{prepare_cosine_matrix.Rd => prepare_similarity_matrix.Rd} (52%) diff --git a/NAMESPACE b/NAMESPACE index 988e48c..43da9ab 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,15 +2,19 @@ export(add_cluster_info) export(colScaleM) +export(connectivity_similarity) export(cosine_sim) export(dataset_names) export(geomScaleM) export(id2char) +export(jaccard_sim) export(partner_summary2adjacency_matrix) export(prepare_cosine_matrix) +export(prepare_similarity_matrix) export(register_dataset) export(rowScaleM) importFrom(grDevices,hcl.colors) +importFrom(methods,as) importFrom(stats,as.dist) importFrom(stats,hclust) importFrom(stats,heatmap) diff --git a/R/cosine.R b/R/cosine.R index dadd65d..3ec1e56 100644 --- a/R/cosine.R +++ b/R/cosine.R @@ -34,21 +34,22 @@ cosine_sim <- function(x, sparse=FALSE, transpose=FALSE) { } -#' Cosine matrix utility functions +#' Prepare a similarity matrix from input/output connectivity #' #' @description These functions are intended for use by package authors rather -#' than end users. +#' than end users. \code{prepare_cosine_matrix} is an alias for +#' \code{prepare_similarity_matrix} retained for backwards compatibility. #' #' @param x A matrix or a named list of input/output matrices #' @param partners Whether to select input or output matrices when both are #' available -#' @param action Whether to zero out or drop any NA values in the cosine matrix -#' (these may be present when some columns have no entries) +#' @param action Whether to zero out or drop any NA values in the similarity +#' matrix (these may be present when some columns have no entries) #' #' @return A matrix. When both inputs and outputs are used these will be #' weighted by the total number of input and output synapses. #' @export -prepare_cosine_matrix <- function(x, partners=c("inputs", "outputs"), action=c("zero", 'drop')) { +prepare_similarity_matrix <- function(x, partners=c("inputs", "outputs"), action=c("zero", 'drop')) { x <- fix_nas(x, action=action) if(is.list(x)) { x <- if(length(partners)==2) { @@ -63,6 +64,13 @@ prepare_cosine_matrix <- function(x, partners=c("inputs", "outputs"), action=c(" x } +#' @rdname prepare_similarity_matrix +#' @usage prepare_cosine_matrix(x, partners = c("inputs", "outputs"), action = c("zero", "drop")) +#' @export +prepare_cosine_matrix <- function(x, partners=c("inputs", "outputs"), action=c("zero", 'drop')) { + prepare_similarity_matrix(x, partners=partners, action=action) +} + fix_nas <- function(x, action=c("zero", 'drop')) { action=match.arg(action) if(is.list(x)) { diff --git a/R/heatmap.R b/R/heatmap.R index 017b3f7..8307136 100644 --- a/R/heatmap.R +++ b/R/heatmap.R @@ -32,13 +32,14 @@ custom_interactive_heatmap <- function(hm) { shiny::shinyApp(ui, server) } -# private function to draw a cosine heatmap using either the basic stats::heatmap -# or InteractiveComplexHeatmap +# private function to draw a connectivity heatmap using either the basic +# stats::heatmap or InteractiveComplexHeatmap # #' @importFrom stats heatmap as.dist hclust #' @importFrom grDevices hcl.colors -cosine_heatmap <- function(x, labRow=rownames(x), interactive=FALSE, +connectivity_heatmap <- function(x, labRow=rownames(x), interactive=FALSE, heatmap=TRUE, col=hcl.colors(12, "YlOrRd", rev = TRUE), + distfun=function(x) as.dist(1-x), method=c("ward.D", "single", "complete", "average", "mcquitty", "median", "centroid", "ward.D2"), ...) { @@ -57,17 +58,19 @@ cosine_heatmap <- function(x, labRow=rownames(x), interactive=FALSE, x, row_labels=labRow, col=col, - cluster_rows=function(x,...) hclust(as.dist(1-x), method=method,...), - cluster_columns=function(x,...) hclust(as.dist(1-x), method=method,...), + cluster_rows=function(x,...) hclust(distfun(x), method=method,...), + cluster_columns=function(x,...) hclust(distfun(x), method=method,...), ... ) custom_interactive_heatmap(hm) } else if(isTRUE(heatmap)) { FUN(x, - distfun = function(x) as.dist(1-x), + distfun = distfun, hclustfun = function(...) hclust(..., method=method), symm = T, keep.dendro = T, labRow=labRow, col=col, ...) } else { - hclust(as.dist(1-x), method = method, ...) + hclust(distfun(x), method = method, ...) } } + +cosine_heatmap <- connectivity_heatmap diff --git a/man/prepare_cosine_matrix.Rd b/man/prepare_similarity_matrix.Rd similarity index 52% rename from man/prepare_cosine_matrix.Rd rename to man/prepare_similarity_matrix.Rd index dc8710e..e67c19f 100644 --- a/man/prepare_cosine_matrix.Rd +++ b/man/prepare_similarity_matrix.Rd @@ -1,14 +1,17 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/cosine.R -\name{prepare_cosine_matrix} +\name{prepare_similarity_matrix} +\alias{prepare_similarity_matrix} \alias{prepare_cosine_matrix} -\title{Cosine matrix utility functions} +\title{Prepare a similarity matrix from input/output connectivity} \usage{ -prepare_cosine_matrix( +prepare_similarity_matrix( x, partners = c("inputs", "outputs"), action = c("zero", "drop") ) + +prepare_cosine_matrix(x, partners = c("inputs", "outputs"), action = c("zero", "drop")) } \arguments{ \item{x}{A matrix or a named list of input/output matrices} @@ -16,8 +19,8 @@ prepare_cosine_matrix( \item{partners}{Whether to select input or output matrices when both are available} -\item{action}{Whether to zero out or drop any NA values in the cosine matrix -(these may be present when some columns have no entries)} +\item{action}{Whether to zero out or drop any NA values in the similarity +matrix (these may be present when some columns have no entries)} } \value{ A matrix. When both inputs and outputs are used these will be @@ -25,5 +28,6 @@ A matrix. When both inputs and outputs are used these will be } \description{ These functions are intended for use by package authors rather - than end users. + than end users. \code{prepare_cosine_matrix} is an alias for + \code{prepare_similarity_matrix} retained for backwards compatibility. } From 2aa51da5f21028dd8c4cb9749fcb4f1da7f34ffd Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 31 Mar 2026 17:58:27 +0100 Subject: [PATCH 02/10] Add methods dependency, update wordlist and build config - Add methods to Imports (needed for as() in similarity code) - Add samples/ to .gitignore and .Rbuildignore - Add Jaccard, binarised etc to spelling wordlist --- .Rbuildignore | 1 + .gitignore | 1 + DESCRIPTION | 5 +++-- inst/WORDLIST | 6 ++++-- 4 files changed, 9 insertions(+), 4 deletions(-) diff --git a/.Rbuildignore b/.Rbuildignore index 22badcd..afdc1d0 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -7,3 +7,4 @@ ^\.github$ ^codecov\.yml$ ^README\.Rmd$ +^samples$ diff --git a/.gitignore b/.gitignore index 0d7f03b..ade91e5 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ .Ruserdata docs inst/doc +samples/ diff --git a/DESCRIPTION b/DESCRIPTION index 8cb6052..9ae0c1a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,12 +13,13 @@ Description: Provides data source agnostic utility functions to support License: GPL (>= 3) Encoding: UTF-8 LazyData: true -RoxygenNote: 7.3.2 -Imports: +RoxygenNote: 7.3.3 +Imports: bit64, checkmate, glue, Matrix (>= 1.5-3), + methods, grDevices, stats Suggests: diff --git a/inst/WORDLIST b/inst/WORDLIST index eebe043..985f70b 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,16 +1,18 @@ CMD Codecov Connectomics +Jaccard Lifecycle Natverse ORCID PNs VNC +banc +binarised +coconatfly connectome connectomics etc -fafbseg flywire natverse -neuprintr org From a09a46673258c24d6d58f8d14da862a2cac0bbf5 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 31 Mar 2026 17:58:40 +0100 Subject: [PATCH 03/10] Add jaccard_sim and connectivity_similarity dispatcher - jaccard_sim() with weighted parameter for binary and weighted Jaccard - connectivity_similarity() as generic dispatcher (cosine/jaccard/weighted_jaccard) - Binary Jaccard uses sparse crossprod with in-place slot transformation - Weighted Jaccard uses threshold decomposition (sum of sparse crossprods) - Comprehensive tests including hand-computed values and naive reference --- R/similarity.R | 124 +++++++++++++++ man/connectivity_similarity.Rd | 45 ++++++ man/jaccard_sim.Rd | 47 ++++++ tests/testthat/_snaps/similarity.md | 24 +++ tests/testthat/test-similarity.R | 226 ++++++++++++++++++++++++++++ 5 files changed, 466 insertions(+) create mode 100644 R/similarity.R create mode 100644 man/connectivity_similarity.Rd create mode 100644 man/jaccard_sim.Rd create mode 100644 tests/testthat/_snaps/similarity.md create mode 100644 tests/testthat/test-similarity.R diff --git a/R/similarity.R b/R/similarity.R new file mode 100644 index 0000000..17e9ce4 --- /dev/null +++ b/R/similarity.R @@ -0,0 +1,124 @@ +#' Generic connectivity similarity between columns (or rows) of a matrix +#' +#' @description \code{connectivity_similarity} computes pairwise similarity +#' between columns (or rows) of a matrix using a pluggable metric. This is the +#' main generic function for connectivity-based clustering; individual metric +#' functions like \code{\link{cosine_sim}} and \code{\link{jaccard_sim}} can +#' also be used directly. +#' +#' @param x A (sparse) matrix, typically an adjacency matrix from +#' \code{\link{partner_summary2adjacency_matrix}} +#' @param metric Character specifying the similarity metric. One of +#' \code{"cosine"}, \code{"jaccard"}, or \code{"weighted_jaccard"}. +#' @param sparse Whether to return a sparse matrix (default \code{FALSE}) +#' @param transpose When \code{FALSE} (the default) calculates similarity +#' between columns. When \code{TRUE} calculates similarity between rows. +#' +#' @return A square similarity matrix with values in \code{[0,1]}. +#' @export +#' @seealso \code{\link{cosine_sim}}, \code{\link{jaccard_sim}} +#' @examples +#' da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) +#' am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') +#' connectivity_similarity(am, metric="cosine") +#' connectivity_similarity(am, metric="jaccard") +#' connectivity_similarity(am, metric="weighted_jaccard") +connectivity_similarity <- function(x, metric = c("cosine", "jaccard", "weighted_jaccard"), + sparse = FALSE, transpose = FALSE) { + metric <- match.arg(metric) + switch(metric, + cosine = cosine_sim(x, sparse = sparse, transpose = transpose), + jaccard = jaccard_sim(x, weighted = FALSE, sparse = sparse, transpose = transpose), + weighted_jaccard = jaccard_sim(x, weighted = TRUE, sparse = sparse, transpose = transpose) + ) +} + + +#' Jaccard similarity for sparse or dense matrices +#' +#' @description Computes pairwise Jaccard similarity between columns (or rows) +#' of a matrix. When \code{weighted=FALSE}, uses binary Jaccard (presence/ +#' absence). When \code{weighted=TRUE}, uses the generalised (weighted) Jaccard +#' index: \code{sum(min(a,b)) / sum(max(a,b))}. +#' +#' @details Both variants are optimised for sparse matrices. The binary variant +#' uses \code{Matrix::crossprod} on the binarised matrix for efficient +#' intersection/union computation. The weighted variant uses the identity +#' \code{min(a,b) = (a + b - |a - b|) / 2} to leverage sparse matrix +#' arithmetic. +#' +#' @param x A data matrix suitable for clustering (non-negative values expected) +#' @param weighted If \code{FALSE} (the default), compute binary Jaccard +#' similarity. If \code{TRUE}, compute weighted (generalised) Jaccard +#' similarity. +#' @param sparse Whether to return a sparse matrix (default \code{FALSE}) +#' @param transpose When \code{FALSE} (the default) calculates similarity +#' between columns. When \code{TRUE} calculates similarity between rows. +#' +#' @return A square similarity matrix with values in \code{[0,1]}. +#' @importFrom methods as +#' @export +#' @seealso \code{\link{cosine_sim}}, \code{\link{connectivity_similarity}} +#' @examples +#' da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) +#' am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') +#' # Binary Jaccard +#' jaccard_sim(am) +#' # Weighted Jaccard +#' jaccard_sim(am, weighted=TRUE) +jaccard_sim <- function(x, weighted = FALSE, sparse = FALSE, transpose = FALSE) { + cx <- class(x) + if (!is.matrix(x) && !isTRUE(attr(cx, "package") == "Matrix")) + stop("I don't recognise that as a matrix!") + if (!inherits(x, "dgCMatrix")) + x <- as(x, "dgCMatrix") + crossfun <- if(transpose) Matrix::tcrossprod else Matrix::crossprod + n <- if(transpose) nrow(x) else ncol(x) + nms <- if(transpose) rownames(x) else colnames(x) + + if (!weighted) { + # Binary Jaccard: J = I / (s_i + s_j - I) + # Binarise, crossprod for intersection, transform values in place + b <- x + b@x <- rep(1, length(b@x)) + A <- as(crossfun(b), "generalMatrix") + sizes <- if(transpose) Matrix::rowSums(b) else diff(b@p) + col_idx <- rep(seq_along(diff(A@p)), diff(A@p)) + row_idx <- A@i + 1L + isect <- A@x + A@x <- isect / (sizes[row_idx] + sizes[col_idx] - isect) + Matrix::diag(A) <- 1 + sim <- A + } else { + # Weighted Jaccard: sum(min(a,b)) / sum(max(a,b)) + # For non-negative values: min(a,b) = sum_{t=1}^{max} I(a>=t)*I(b>=t) + # Each threshold gives a sparse crossprod (same fast op as cosine). + max_val <- if (length(x@x)) max(x@x) else 0 + if (max_val == 0) { + sim <- Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, + dims = c(n, n), dimnames = list(nms, nms)) + } else { + cs <- if(transpose) Matrix::rowSums(x) else Matrix::colSums(x) + vals <- x@x + p <- x@p + ri <- x@i + col_idx <- rep(seq_len(ncol(x)), diff(p)) + min_sums <- matrix(0, n, n) + for (t in seq_len(max_val)) { + keep <- vals >= t + bt <- Matrix::sparseMatrix( + i = ri[keep] + 1L, j = col_idx[keep], + x = 1, dims = x@Dim, dimnames = x@Dimnames) + min_sums <- min_sums + as.matrix(crossfun(bt)) + } + max_sums <- outer(cs, cs, "+") - min_sums + sim <- min_sums / max_sums + sim[is.nan(sim)] <- 0 + diag(sim) <- 1 + dimnames(sim) <- list(nms, nms) + sim <- Matrix::Matrix(sim, sparse = TRUE) + } + } + + if (sparse) sim else as.matrix(sim) +} diff --git a/man/connectivity_similarity.Rd b/man/connectivity_similarity.Rd new file mode 100644 index 0000000..6a838a3 --- /dev/null +++ b/man/connectivity_similarity.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/similarity.R +\name{connectivity_similarity} +\alias{connectivity_similarity} +\title{Generic connectivity similarity between columns (or rows) of a matrix} +\usage{ +connectivity_similarity( + x, + metric = c("cosine", "jaccard", "weighted_jaccard"), + sparse = FALSE, + transpose = FALSE +) +} +\arguments{ +\item{x}{A (sparse) matrix, typically an adjacency matrix from +\code{\link{partner_summary2adjacency_matrix}}} + +\item{metric}{Character specifying the similarity metric. One of +\code{"cosine"}, \code{"jaccard"}, or \code{"weighted_jaccard"}.} + +\item{sparse}{Whether to return a sparse matrix (default \code{FALSE})} + +\item{transpose}{When \code{FALSE} (the default) calculates similarity +between columns. When \code{TRUE} calculates similarity between rows.} +} +\value{ +A square similarity matrix with values in \code{[0,1]}. +} +\description{ +\code{connectivity_similarity} computes pairwise similarity + between columns (or rows) of a matrix using a pluggable metric. This is the + main generic function for connectivity-based clustering; individual metric + functions like \code{\link{cosine_sim}} and \code{\link{jaccard_sim}} can + also be used directly. +} +\examples{ +da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) +am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') +connectivity_similarity(am, metric="cosine") +connectivity_similarity(am, metric="jaccard") +connectivity_similarity(am, metric="weighted_jaccard") +} +\seealso{ +\code{\link{cosine_sim}}, \code{\link{jaccard_sim}} +} diff --git a/man/jaccard_sim.Rd b/man/jaccard_sim.Rd new file mode 100644 index 0000000..891889d --- /dev/null +++ b/man/jaccard_sim.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/similarity.R +\name{jaccard_sim} +\alias{jaccard_sim} +\title{Jaccard similarity for sparse or dense matrices} +\usage{ +jaccard_sim(x, weighted = FALSE, sparse = FALSE, transpose = FALSE) +} +\arguments{ +\item{x}{A data matrix suitable for clustering (non-negative values expected)} + +\item{weighted}{If \code{FALSE} (the default), compute binary Jaccard +similarity. If \code{TRUE}, compute weighted (generalised) Jaccard +similarity.} + +\item{sparse}{Whether to return a sparse matrix (default \code{FALSE})} + +\item{transpose}{When \code{FALSE} (the default) calculates similarity +between columns. When \code{TRUE} calculates similarity between rows.} +} +\value{ +A square similarity matrix with values in \code{[0,1]}. +} +\description{ +Computes pairwise Jaccard similarity between columns (or rows) + of a matrix. When \code{weighted=FALSE}, uses binary Jaccard (presence/ + absence). When \code{weighted=TRUE}, uses the generalised (weighted) Jaccard + index: \code{sum(min(a,b)) / sum(max(a,b))}. +} +\details{ +Both variants are optimised for sparse matrices. The binary variant + uses \code{Matrix::crossprod} on the binarised matrix for efficient + intersection/union computation. The weighted variant uses the identity + \code{min(a,b) = (a + b - |a - b|) / 2} to leverage sparse matrix + arithmetic. +} +\examples{ +da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) +am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') +# Binary Jaccard +jaccard_sim(am) +# Weighted Jaccard +jaccard_sim(am, weighted=TRUE) +} +\seealso{ +\code{\link{cosine_sim}}, \code{\link{connectivity_similarity}} +} diff --git a/tests/testthat/_snaps/similarity.md b/tests/testthat/_snaps/similarity.md new file mode 100644 index 0000000..bf1ce83 --- /dev/null +++ b/tests/testthat/_snaps/similarity.md @@ -0,0 +1,24 @@ +# jaccard_sim snapshot on sample data + + Code + jaccard_sim(am) + Output + 1796817841 1797505019 1796818119 1827516355 818983130 + 1796817841 1.00000000 0.06451613 0.00000000 0.00000000 0.02500000 + 1797505019 0.06451613 1.00000000 0.00000000 0.00000000 0.06666667 + 1796818119 0.00000000 0.00000000 1.00000000 0.00000000 0.04545455 + 1827516355 0.00000000 0.00000000 0.00000000 1.00000000 0.06896552 + 818983130 0.02500000 0.06666667 0.04545455 0.06896552 1.00000000 + +--- + + Code + jaccard_sim(am, weighted = TRUE) + Output + 1796817841 1797505019 1796818119 1827516355 818983130 + 1796817841 1.00000000 0.05385996 0.00000000 0.00000000 0.02013423 + 1797505019 0.05385996 1.00000000 0.00000000 0.00000000 0.05575540 + 1796818119 0.00000000 0.00000000 1.00000000 0.00000000 0.03699284 + 1827516355 0.00000000 0.00000000 0.00000000 1.00000000 0.08255159 + 818983130 0.02013423 0.05575540 0.03699284 0.08255159 1.00000000 + diff --git a/tests/testthat/test-similarity.R b/tests/testthat/test-similarity.R new file mode 100644 index 0000000..4f66621 --- /dev/null +++ b/tests/testthat/test-similarity.R @@ -0,0 +1,226 @@ +# Naive reference implementation for verification +naive_jaccard <- function(x, weighted = FALSE) { + n <- ncol(x) + sim <- matrix(0, n, n, dimnames = list(colnames(x), colnames(x))) + for (i in seq_len(n)) { + for (j in seq(i, n)) { + a <- x[, i] + b <- x[, j] + if (!weighted) { + ab <- (a != 0) & (b != 0) + aub <- (a != 0) | (b != 0) + sim[i, j] <- sim[j, i] <- if (sum(aub) == 0) 0 else sum(ab) / sum(aub) + } else { + min_sum <- sum(pmin(a, b)) + max_sum <- sum(pmax(a, b)) + sim[i, j] <- sim[j, i] <- if (max_sum == 0) 0 else min_sum / max_sum + } + } + } + diag(sim) <- 1 + sim +} + +test_that("jaccard_sim matches naive implementation on sample data", { + da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) + am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') + + # Binary + jb <- jaccard_sim(am) + jb_naive <- naive_jaccard(as.matrix(am), weighted = FALSE) + expect_equal(jb, jb_naive) + + # Weighted + jw <- jaccard_sim(am, weighted = TRUE) + jw_naive <- naive_jaccard(as.matrix(am), weighted = TRUE) + expect_equal(jw, jw_naive) +}) + +test_that("jaccard_sim hand-computed values are correct", { + # 3 rows (partners), 3 columns (neurons) with known structure + # n1 n2 n3 + # p1 [ 4 2 0 ] + # p2 [ 0 3 3 ] + # p3 [ 1 1 0 ] + m <- matrix(c(4,0,1, 2,3,1, 0,3,0), nrow = 3, ncol = 3) + colnames(m) <- c("n1", "n2", "n3") + + # Binary Jaccard between columns (nonzero patterns): + # n1: {p1, p3} n2: {p1, p2, p3} n3: {p2} + # J(n1,n2) = |{p1,p3}| / |{p1,p2,p3}| = 2/3 + # J(n1,n3) = |{}| / |{p1,p2,p3}| = 0/3 = 0 + # J(n2,n3) = |{p2}| / |{p1,p2,p3}| = 1/3 + jb <- jaccard_sim(m) + expect_equal(jb["n1","n2"], 2/3) + expect_equal(jb["n1","n3"], 0) + expect_equal(jb["n2","n3"], 1/3) + + # Weighted Jaccard: + # J_w(n1,n2) = (min(4,2)+min(0,3)+min(1,1)) / (max(4,2)+max(0,3)+max(1,1)) + # = (2+0+1) / (4+3+1) = 3/8 + # J_w(n1,n3) = (min(4,0)+min(0,3)+min(1,0)) / (max(4,0)+max(0,3)+max(1,0)) + # = (0+0+0) / (4+3+1) = 0/8 = 0 + # J_w(n2,n3) = (min(2,0)+min(3,3)+min(1,0)) / (max(2,0)+max(3,3)+max(1,0)) + # = (0+3+0) / (2+3+1) = 3/6 = 1/2 + jw <- jaccard_sim(m, weighted = TRUE) + expect_equal(jw["n1","n2"], 3/8) + expect_equal(jw["n1","n3"], 0) + expect_equal(jw["n2","n3"], 1/2) +}) + +test_that("jaccard_sim: identical columns give 1, disjoint columns give 0", { + # n1 n2 n3 + # p1 [ 5 5 0 ] + # p2 [ 3 3 0 ] + # p3 [ 0 0 7 ] + m <- Matrix::Matrix(c(5,3,0, 5,3,0, 0,0,7), nrow = 3, ncol = 3, sparse = TRUE) + colnames(m) <- c("n1", "n2", "n3") + + # n1 and n2 are identical → both Jaccard variants should be 1 + # n1/n2 and n3 are completely disjoint → both should be 0 + jb <- jaccard_sim(m) + expect_equal(jb["n1","n2"], 1) + expect_equal(jb["n1","n3"], 0) + expect_equal(jb["n2","n3"], 0) + + jw <- jaccard_sim(m, weighted = TRUE) + expect_equal(jw["n1","n2"], 1) + expect_equal(jw["n1","n3"], 0) + expect_equal(jw["n2","n3"], 0) +}) + +test_that("jaccard_sim: proportional columns differ from cosine", { + # a = 2*b: cosine similarity = 1, but weighted Jaccard = 0.5 + m <- matrix(c(2,4,6, 1,2,3), nrow = 3, ncol = 2) + colnames(m) <- c("a", "b") + + # Cosine should be 1 (proportional vectors have cosine sim = 1) + cs <- cosine_sim(m) + expect_equal(cs["a","b"], 1) + + # Binary Jaccard should be 1 (same nonzero pattern) + jb <- jaccard_sim(m) + expect_equal(jb["a","b"], 1) + + # Weighted Jaccard: sum(min)/sum(max) = (1+2+3)/(2+4+6) = 6/12 = 0.5 + jw <- jaccard_sim(m, weighted = TRUE) + expect_equal(jw["a","b"], 0.5) +}) + +test_that("jaccard_sim: binary matrix gives same result for binary and weighted", { + m <- Matrix::Matrix(c(1,0,1,0, 0,1,1,0, 1,1,0,1), nrow = 4, ncol = 3, sparse = TRUE) + colnames(m) <- c("a", "b", "c") + + jb <- jaccard_sim(m) + jw <- jaccard_sim(m, weighted = TRUE) + expect_equal(jb, jw) +}) + +test_that("jaccard_sim works with dense matrix input", { + m_sparse <- Matrix::Matrix(c(4,0,1, 2,3,1, 0,3,0), nrow = 3, ncol = 3, sparse = TRUE) + colnames(m_sparse) <- c("n1", "n2", "n3") + m_dense <- as.matrix(m_sparse) + + expect_equal(jaccard_sim(m_dense), jaccard_sim(m_sparse)) + expect_equal(jaccard_sim(m_dense, weighted = TRUE), + jaccard_sim(m_sparse, weighted = TRUE)) +}) + +test_that("jaccard_sim sparse output parameter works", { + m <- Matrix::Matrix(c(4,0,1, 2,3,1, 0,3,0), nrow = 3, ncol = 3, sparse = TRUE) + colnames(m) <- c("n1", "n2", "n3") + + js <- jaccard_sim(m, sparse = TRUE) + expect_true(inherits(js, "Matrix")) + expect_equal(as.matrix(js), jaccard_sim(m, sparse = FALSE)) + + jws <- jaccard_sim(m, weighted = TRUE, sparse = TRUE) + expect_true(inherits(jws, "Matrix")) + expect_equal(as.matrix(jws), jaccard_sim(m, weighted = TRUE, sparse = FALSE)) +}) + +test_that("jaccard_sim transpose parameter works", { + m <- Matrix::Matrix(c(4,0,1, 2,3,1, 0,3,0), nrow = 3, ncol = 3, sparse = TRUE) + colnames(m) <- c("n1", "n2", "n3") + rownames(m) <- c("p1", "p2", "p3") + + # Default: similarity between columns (3x3) + jc <- jaccard_sim(m) + expect_equal(dim(jc), c(3L, 3L)) + expect_equal(rownames(jc), c("n1", "n2", "n3")) + + # Transpose: similarity between rows (3x3 here too, but different values) + jt <- jaccard_sim(m, transpose = TRUE) + expect_equal(dim(jt), c(3L, 3L)) + expect_equal(rownames(jt), c("p1", "p2", "p3")) + + # Transpose should give same result as manually transposing + expect_equal(jt, jaccard_sim(Matrix::t(m))) +}) + +test_that("jaccard_sim snapshot on sample data", { + da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) + am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') + + expect_snapshot(jaccard_sim(am)) + expect_snapshot(jaccard_sim(am, weighted = TRUE)) +}) + +test_that("jaccard_sim handles edge cases", { + # Single column matrix + m1 <- Matrix::Matrix(c(1, 0, 3), ncol = 1, sparse = TRUE) + colnames(m1) <- "a" + j1 <- jaccard_sim(m1) + expect_equal(j1, matrix(1, 1, 1, dimnames = list("a", "a"))) + + j1w <- jaccard_sim(m1, weighted = TRUE) + expect_equal(j1w, matrix(1, 1, 1, dimnames = list("a", "a"))) + + # All-zero columns: diagonal should still be 1, off-diagonal 0 + m2 <- Matrix::Matrix(0, nrow = 3, ncol = 2, sparse = TRUE) + colnames(m2) <- c("a", "b") + j2 <- jaccard_sim(m2) + expect_true(all(diag(j2) == 1)) + expect_equal(j2["a", "b"], 0) + + j2w <- jaccard_sim(m2, weighted = TRUE) + expect_true(all(diag(j2w) == 1)) + expect_equal(j2w["a", "b"], 0) + + # Rejects non-matrix input + expect_error(jaccard_sim(1:10), "matrix") + expect_error(jaccard_sim(data.frame(a = 1:3)), "matrix") +}) + +test_that("connectivity_similarity dispatches correctly", { + da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) + am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') + + expect_equal(connectivity_similarity(am, metric = "cosine"), cosine_sim(am)) + expect_equal(connectivity_similarity(am, metric = "jaccard"), jaccard_sim(am)) + expect_equal(connectivity_similarity(am, metric = "weighted_jaccard"), + jaccard_sim(am, weighted = TRUE)) +}) + +test_that("prepare_similarity_matrix is an alias for prepare_cosine_matrix", { + da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) + am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') + cm <- cosine_sim(am, transpose = TRUE) + expect_equal(prepare_similarity_matrix(cm), prepare_cosine_matrix(cm)) +}) + +test_that("connectivity_heatmap works with distfun", { + da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) + am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') + jw <- jaccard_sim(am, weighted = TRUE, transpose = TRUE) + pm <- prepare_similarity_matrix(jw) + + # Default distfun (1-x) + expect_silent(cl <- connectivity_heatmap(pm, heatmap = FALSE)) + expect_true(inherits(cl, "hclust")) + + # Custom distfun + expect_silent(cl2 <- connectivity_heatmap(pm, heatmap = FALSE, + distfun = function(x) as.dist(1 - x))) + expect_equal(cl, cl2) +}) From 62b6dfe197f574d3856e543ecf5b6db34ad33d9f Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 31 Mar 2026 18:53:39 +0100 Subject: [PATCH 04/10] Add a Tanimoto similarity function A form of weighted jaccard but potentially faster to compute --- NAMESPACE | 1 + R/similarity.R | 63 ++++++++++++++++++++++++-- man/connectivity_similarity.Rd | 9 ++-- man/tanimoto_sim.Rd | 35 +++++++++++++++ tests/testthat/test-similarity.R | 76 ++++++++++++++++++++++++++++++++ 5 files changed, 177 insertions(+), 7 deletions(-) create mode 100644 man/tanimoto_sim.Rd diff --git a/NAMESPACE b/NAMESPACE index 43da9ab..c95d85f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,6 +13,7 @@ export(prepare_cosine_matrix) export(prepare_similarity_matrix) export(register_dataset) export(rowScaleM) +export(tanimoto_sim) importFrom(grDevices,hcl.colors) importFrom(methods,as) importFrom(stats,as.dist) diff --git a/R/similarity.R b/R/similarity.R index 17e9ce4..01dd78d 100644 --- a/R/similarity.R +++ b/R/similarity.R @@ -9,27 +9,31 @@ #' @param x A (sparse) matrix, typically an adjacency matrix from #' \code{\link{partner_summary2adjacency_matrix}} #' @param metric Character specifying the similarity metric. One of -#' \code{"cosine"}, \code{"jaccard"}, or \code{"weighted_jaccard"}. +#' \code{"cosine"}, \code{"jaccard"}, \code{"weighted_jaccard"}, or +#' \code{"tanimoto"}. #' @param sparse Whether to return a sparse matrix (default \code{FALSE}) #' @param transpose When \code{FALSE} (the default) calculates similarity #' between columns. When \code{TRUE} calculates similarity between rows. #' #' @return A square similarity matrix with values in \code{[0,1]}. #' @export -#' @seealso \code{\link{cosine_sim}}, \code{\link{jaccard_sim}} +#' @seealso \code{\link{cosine_sim}}, \code{\link{jaccard_sim}}, +#' \code{\link{tanimoto_sim}} #' @examples #' da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) #' am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') #' connectivity_similarity(am, metric="cosine") #' connectivity_similarity(am, metric="jaccard") #' connectivity_similarity(am, metric="weighted_jaccard") -connectivity_similarity <- function(x, metric = c("cosine", "jaccard", "weighted_jaccard"), +#' connectivity_similarity(am, metric="tanimoto") +connectivity_similarity <- function(x, metric = c("cosine", "jaccard", "weighted_jaccard", "tanimoto"), sparse = FALSE, transpose = FALSE) { metric <- match.arg(metric) switch(metric, cosine = cosine_sim(x, sparse = sparse, transpose = transpose), jaccard = jaccard_sim(x, weighted = FALSE, sparse = sparse, transpose = transpose), - weighted_jaccard = jaccard_sim(x, weighted = TRUE, sparse = sparse, transpose = transpose) + weighted_jaccard = jaccard_sim(x, weighted = TRUE, sparse = sparse, transpose = transpose), + tanimoto = tanimoto_sim(x, sparse = sparse, transpose = transpose) ) } @@ -122,3 +126,54 @@ jaccard_sim <- function(x, weighted = FALSE, sparse = FALSE, transpose = FALSE) if (sparse) sim else as.matrix(sim) } + + +#' Tanimoto (extended Jaccard) similarity for sparse or dense matrices +#' +#' @description Computes pairwise Tanimoto similarity between columns (or rows) +#' of a matrix. Also known as the extended Jaccard coefficient, defined as: +#' \code{dot(a,b) / (||a||^2 + ||b||^2 - dot(a,b))}. On binary data this +#' reduces to the standard Jaccard index. Computed via a single +#' \code{crossprod} call, so performance is comparable to \code{cosine_sim}. +#' +#' @param x A data matrix suitable for clustering (non-negative values expected) +#' @param sparse Whether to return a sparse matrix (default \code{FALSE}) +#' @param transpose When \code{FALSE} (the default) calculates similarity +#' between columns. When \code{TRUE} calculates similarity between rows. +#' +#' @return A square similarity matrix with values in \code{[0,1]}. +#' @importFrom methods as +#' @export +#' @seealso \code{\link{jaccard_sim}}, \code{\link{cosine_sim}}, +#' \code{\link{connectivity_similarity}} +#' @examples +#' da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) +#' am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') +#' tanimoto_sim(am) +tanimoto_sim <- function(x, sparse = FALSE, transpose = FALSE) { + cx <- class(x) + if (!is.matrix(x) && !isTRUE(attr(cx, "package") == "Matrix")) + stop("I don't recognise that as a matrix!") + if (!inherits(x, "dgCMatrix")) + x <- as(x, "dgCMatrix") + crossfun <- if (transpose) Matrix::tcrossprod else Matrix::crossprod + + # dot(a,b) for all pairs via crossprod (returns symmetric dsCMatrix) + A <- crossfun(x) + norms2 <- Matrix::diag(A) + + # T(a,b) = dot(a,b) / (||a||^2 + ||b||^2 - dot(a,b)) + # Work directly on the symmetric matrix (upper triangle stored) + if (length(A@x)) { + col_idx <- rep(seq_along(diff(A@p)), diff(A@p)) + row_idx <- A@i + 1L + denom <- norms2[row_idx] + norms2[col_idx] - A@x + nonzero <- denom != 0 + A@x[nonzero] <- A@x[nonzero] / denom[nonzero] + A@x[!nonzero] <- 0 + } + Matrix::diag(A) <- 1 + sim <- A + + if (sparse) sim else as.matrix(sim) +} diff --git a/man/connectivity_similarity.Rd b/man/connectivity_similarity.Rd index 6a838a3..c433e21 100644 --- a/man/connectivity_similarity.Rd +++ b/man/connectivity_similarity.Rd @@ -6,7 +6,7 @@ \usage{ connectivity_similarity( x, - metric = c("cosine", "jaccard", "weighted_jaccard"), + metric = c("cosine", "jaccard", "weighted_jaccard", "tanimoto"), sparse = FALSE, transpose = FALSE ) @@ -16,7 +16,8 @@ connectivity_similarity( \code{\link{partner_summary2adjacency_matrix}}} \item{metric}{Character specifying the similarity metric. One of -\code{"cosine"}, \code{"jaccard"}, or \code{"weighted_jaccard"}.} +\code{"cosine"}, \code{"jaccard"}, \code{"weighted_jaccard"}, or +\code{"tanimoto"}.} \item{sparse}{Whether to return a sparse matrix (default \code{FALSE})} @@ -39,7 +40,9 @@ am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = ' connectivity_similarity(am, metric="cosine") connectivity_similarity(am, metric="jaccard") connectivity_similarity(am, metric="weighted_jaccard") +connectivity_similarity(am, metric="tanimoto") } \seealso{ -\code{\link{cosine_sim}}, \code{\link{jaccard_sim}} +\code{\link{cosine_sim}}, \code{\link{jaccard_sim}}, + \code{\link{tanimoto_sim}} } diff --git a/man/tanimoto_sim.Rd b/man/tanimoto_sim.Rd new file mode 100644 index 0000000..2d65085 --- /dev/null +++ b/man/tanimoto_sim.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/similarity.R +\name{tanimoto_sim} +\alias{tanimoto_sim} +\title{Tanimoto (extended Jaccard) similarity for sparse or dense matrices} +\usage{ +tanimoto_sim(x, sparse = FALSE, transpose = FALSE) +} +\arguments{ +\item{x}{A data matrix suitable for clustering (non-negative values expected)} + +\item{sparse}{Whether to return a sparse matrix (default \code{FALSE})} + +\item{transpose}{When \code{FALSE} (the default) calculates similarity +between columns. When \code{TRUE} calculates similarity between rows.} +} +\value{ +A square similarity matrix with values in \code{[0,1]}. +} +\description{ +Computes pairwise Tanimoto similarity between columns (or rows) + of a matrix. Also known as the extended Jaccard coefficient, defined as: + \code{dot(a,b) / (||a||^2 + ||b||^2 - dot(a,b))}. On binary data this + reduces to the standard Jaccard index. Computed via a single + \code{crossprod} call, so performance is comparable to \code{cosine_sim}. +} +\examples{ +da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) +am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') +tanimoto_sim(am) +} +\seealso{ +\code{\link{jaccard_sim}}, \code{\link{cosine_sim}}, + \code{\link{connectivity_similarity}} +} diff --git a/tests/testthat/test-similarity.R b/tests/testthat/test-similarity.R index 4f66621..30f3763 100644 --- a/tests/testthat/test-similarity.R +++ b/tests/testthat/test-similarity.R @@ -192,6 +192,81 @@ test_that("jaccard_sim handles edge cases", { expect_error(jaccard_sim(data.frame(a = 1:3)), "matrix") }) +test_that("tanimoto_sim matches naive implementation", { + naive_tanimoto <- function(x) { + n <- ncol(x) + sim <- matrix(0, n, n, dimnames = list(colnames(x), colnames(x))) + for (i in seq_len(n)) { + for (j in seq(i, n)) { + a <- x[, i]; b <- x[, j] + d <- sum(a * b) + denom <- sum(a^2) + sum(b^2) - d + sim[i, j] <- sim[j, i] <- if (denom == 0) 0 else d / denom + } + } + diag(sim) <- 1 + sim + } + + da2ds15 <- readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) + am <- partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') + expect_equal(tanimoto_sim(am), naive_tanimoto(as.matrix(am))) +}) + +test_that("tanimoto_sim hand-computed values are correct", { + m <- matrix(c(4,0,1, 2,3,1, 0,3,0), nrow = 3, ncol = 3) + colnames(m) <- c("n1", "n2", "n3") + + # T(n1,n2) = dot(n1,n2) / (||n1||^2 + ||n2||^2 - dot(n1,n2)) + # dot = 4*2 + 0*3 + 1*1 = 9 + # ||n1||^2 = 16+0+1 = 17, ||n2||^2 = 4+9+1 = 14 + # T = 9 / (17 + 14 - 9) = 9/22 + ts <- tanimoto_sim(m) + expect_equal(ts["n1","n2"], 9/22) + expect_equal(ts["n1","n3"], 0) # disjoint + # dot(n2,n3) = 0+9+0 = 9, ||n3||^2 = 9 + # T = 9 / (14 + 9 - 9) = 9/14 + expect_equal(ts["n2","n3"], 9/14) +}) + +test_that("tanimoto_sim equals binary jaccard on 0/1 data", { + m <- Matrix::Matrix(c(1,0,1,0, 0,1,1,0, 1,1,0,1), nrow = 4, ncol = 3, sparse = TRUE) + colnames(m) <- c("a", "b", "c") + expect_equal(tanimoto_sim(m), jaccard_sim(m)) +}) + +test_that("tanimoto_sim handles edge cases", { + # Single column + m1 <- Matrix::Matrix(c(1, 0, 3), ncol = 1, sparse = TRUE) + colnames(m1) <- "a" + expect_equal(tanimoto_sim(m1), matrix(1, 1, 1, dimnames = list("a", "a"))) + + # All-zero columns + m2 <- Matrix::Matrix(0, nrow = 3, ncol = 2, sparse = TRUE) + colnames(m2) <- c("a", "b") + t2 <- tanimoto_sim(m2) + expect_true(all(diag(t2) == 1)) + expect_equal(t2["a", "b"], 0) + + # Rejects non-matrix input + expect_error(tanimoto_sim(1:10), "matrix") +}) + +test_that("tanimoto_sim sparse and transpose parameters work", { + m <- Matrix::Matrix(c(4,0,1, 2,3,1, 0,3,0), nrow = 3, ncol = 3, sparse = TRUE) + colnames(m) <- c("n1", "n2", "n3") + rownames(m) <- c("p1", "p2", "p3") + + ts <- tanimoto_sim(m, sparse = TRUE) + expect_true(inherits(ts, "Matrix")) + expect_equal(as.matrix(ts), tanimoto_sim(m)) + + tt <- tanimoto_sim(m, transpose = TRUE) + expect_equal(dim(tt), c(3L, 3L)) + expect_equal(rownames(tt), c("p1", "p2", "p3")) + expect_equal(tt, tanimoto_sim(Matrix::t(m))) +}) + test_that("connectivity_similarity dispatches correctly", { da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') @@ -200,6 +275,7 @@ test_that("connectivity_similarity dispatches correctly", { expect_equal(connectivity_similarity(am, metric = "jaccard"), jaccard_sim(am)) expect_equal(connectivity_similarity(am, metric = "weighted_jaccard"), jaccard_sim(am, weighted = TRUE)) + expect_equal(connectivity_similarity(am, metric = "tanimoto"), tanimoto_sim(am)) }) test_that("prepare_similarity_matrix is an alias for prepare_cosine_matrix", { From c02fb5b6be590e00e7c392c962dfcb8fd7ff0eba Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Wed, 1 Apr 2026 13:03:34 +0100 Subject: [PATCH 05/10] clean up / don't generate pdfs insides tests --- .gitignore | 1 + tests/testthat/Rplots.pdf | Bin 56005 -> 0 bytes tests/testthat/setup.R | 3 +++ 3 files changed, 4 insertions(+) delete mode 100644 tests/testthat/Rplots.pdf create mode 100644 tests/testthat/setup.R diff --git a/.gitignore b/.gitignore index ade91e5..4de88c5 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,4 @@ docs inst/doc samples/ +tests/testthat/Rplots.pdf diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf deleted file mode 100644 index 5619e06475fc0d7acf45caaa5d092cff62ebb161..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 56005 zcmd?S30MbE_7Hzs)L_HaMbpO!jU z&#|_dV>f5PoCR~|&6_`Sl+{|y&OyKP+#MbWa;-OA5xRZfuH8GrOx+{G_eF$*Z0kZX zJ0o`Zg+kxXv$V0C3!cGsVW3AdM>#u>@(bPxwZdB68dEFJn0rZx6%-m|Csd z9^e;dY6B9iH`QdY`zQmH4oV&B7aXokZnXpx91i{rGhF~@w_5KPhKUH>0m_>LnOe`cx3)L6TE06l9K3?n@<7l#OZ|3WcKSgL4+{sK7BnjQ!ljml$rf8l znL?>-g2k<2glN*G?*c#nB5={u&nHYw+Bx&v^5K&*@1hwS!j9j2kY}tHMP)6LJii(s zzxHsCFK&8VVKPbHQot-w^m4KlJT6|wFJQKP5G;^7#7`bevP|UDIs_;MlOz?pljIV~ z9~?dhCyA)MK(b7eCVR=Q;YE6iUbI3?Rdgpai~g*#k_kUV;JtA&3SQ1tynO$vJ0kua zTmX}XUt69@%JMp%T);fkF6)VK-d>LNNuFRrS?0K{Jl-j5Vm3N)rM5kJ{3p-MvKM90VbRAy&a(fB>y z8?)20yM*Zt9h@!owYl^Myl_bw^Xr9KX~Wv<-b7d8CX3>;ve5y~*_5EZM{ItmqKvtg z9qoj3;@t@=Z><~~&9rY#9iK$VVs>^%v=tXN>j%f%rb_-gSiqIZBVKVN3viizz5EiK zPha3HVWLs?3D(#UHb^%uD;n#b*gzXs!0p&4abUVDZZ%KKkG6HAIJGAQFmE;YW`jR+ z`EY`8g+rUbDKi;uAnmYVS^T)ssM#;m5E3T#}yzk-jUZcD|WSfhq z`cn2NxrdpeG(q86gU^4Zer#th24Ci`hVrPt+|UimRzKD*3aPSxQM8$g6NakR{701| z)gtBRd&swTz1Hi|ep}@JU;NHC@`Xvnxqa^Ar0kcn+^^(e@8up}D(oYn&&&kHB%)$M z4gP8i5E$v1p~BC<8lRQkjrK2tWD6cEzv%+fMDGy-nOWVIfDe#dT9XV2d`&0x0g}BO z2`R*SqRxEktUB`nrO@QFkvEJHscTSIg4BSdR40bo+yy?`jYDckb6}A|R+uSHtc2cH z`wHH9`E=;bmoHUmqlcVddqy$OhK-``AE3Ma)aMln)>_CSv~=*#KQw00Uw8Xkhb{w9 z7M0n}IOn2jjmB_3p>l2Ob6}DhN-Po-6Yk-YTJT3-p~zjpw6k{zq0iJlYN5paluBxu zJhZ>Rl$}0P;g_Qn=6}@H=7Wj|=z+F8|4u|qRNYA}$s;=0Pnc>w7@%p|HC*+O8%5a; zEk_27I+-|rC>k|uY}nF%_Keu}zu2mQBTF;=-}B0wj%(%yU=pC^F_F1Jf2TNL%7DxY zBC7FeA^Z`l;YrJ#sXf|MmHz=wueP?2YU3fNxauCBkdb46TT%~e8@knF(@(B3L%RS}P%VLIFYjg)Pu zd-kPkf~ArN#&giQRR0mI|3t3ecm%wzjjA;tP5u6^a%9X?Dno!YYO>wWAC$wB)*MAW zO_`~lRM^NVwZ|2XHY`$Y`2RsygBk}Bk?t9UP7wu2a}5X8n7`LF?FkeRA9MLixzQ;p zn2~8;HItRYDl^K3hw0t^_4usMm2+;D9Q$>p?y30lWrNKtfU+Xu)iP1=2HI0`WL6KR zQgaWgT5}JY*=yzabfbQRr__xhDK-zKcdo`0LY90nBpsRgkaJ`vb|Z$QBO|Zv=fV(o z_Ahp1)%axIu;D|{s0q@CA0*lOLuh&eW@NaZhMFT2PRgAO%}X>>FnF5?+9^V&;;)qx zrgz5`);3INYtYI@IWPW{sg!ya2BsgXdH8^Q1E%6SCsxXx1w+N|-&DIAT6oGkQT*LH2Aeomx*<(>_iyX||tG=FIBl#fBu zehD;>+dAgB4HK!h!I8~y=n8gG?WT(0>!J*&}Nr<@vB;98R#8$1PoaO#OPuOF~!8 z3T512m8f3bEy;1D&I&3;2YMSLG%GaJci9X!hIHjuhFd*R|EEt^{#& zN=KTYs6OlB^MSmd{pwejy;7u)%dta zhz1gK9QGecOu0ei@L$>-e?$2AsntN$_5BX{+`w7k zpvQu`_J$A_1T>_2%wT?1DaZaB6a$(FA3M^zNc@*(L(!-$ExaG%r=>vCH0|qA^=9Ie zkDW2H-7{cMCzB6ZG;}%y=2mLwt{S`;SG4y1z|pX-acywxR95)k**?67-(fSfhUvXR zMz<&SFRlVP#2r$eOhOYK`LwNufxDCg_BG#KfG1?|ZyI#?DeXN|@XG^THh}!#;cv#9 z|IU$alZ@cl;>!bzZmmO*!OnC)TK)7tCr5!-fyS)#bY6DAT^%bg?K534{D8A?bpfhn zn)X4v?x{t6A8`HZ)tv6Nmg@Hc$Z)Ez1fIVrXFBG6+_LWu46WIyzEtQvZqbla@aD$k zK3Tt?E5N^1fPSukoCmEObpA9zEy4T&>_;HyfzU5Ql)9w)*$S9JsHfU%G6yjN_y;b^ zIoSWr6Cfb$J_n{yF>@#l*F;5X?q$xAPug4plV_FSbc@3c>^dor1a-}@`)w`_vM5u| z=Li3x(|_w^_a2@gs5v*Z+}UlhUY3yo_)+K9-}PM%CHWjWdT}v1fAO+)zwhvZK2xte zF1)txl;^po4R1%8yWTMOu8H)v zZ6ok`ra$I<(>L!%jab_`B4YT#oJc5VSqfC3j|)^FROZFHXZdTr9{ca1pxxE*B z8uV|aH{I?S)#Taw4N>O1SSP<8?h~@SY0{}%SHJ3a_53m@X=M|1Jq+c{OR@FI zz15LEGiZ;$TSnIB@Bp|RHcH)S68cQcge_uQ(Np(76VtQY?|szT}7#C&@L)pR@}+zX#15B zbY*3}=a*hdT37P~ixtzopzBgd$Qx?TLaCq08}C@e-^m-C|`9%c=t2_2p_;rZJfHbAMtV>|Qn;F`M{D!Koc11dy{ zQq`bc6fL-n_4TU+hS!y^a>Cx{czoGlAK5T7p<@ye`b@YNp41Y4bZ6PKV9*wyJ$B#& zs413GSy_UOz9wFl0cuT`J)1kN3zJhn2mIA%`G^LavDfI6; zSa@pTX!@V+$~GJ*vE|EP(I8)O&UDOMqR|7DC=qEPVi}}tWihb$Ra${mjHp!qqSTtE ztwUSZM)_XpzQ{cr`95dObx7LktNJGF%Z7?Q^MRCGu5Hx11F{8=p=ZdQGJ)3P`g2|L zJs^Dmd+Xr^TU2+T%e`*otjV0)vNz{^*M=F8TEI^_bK;)w(Yk|o@*z$7`^VtL25;Yg z+U$gmV9t2G^x&n|nN!bcU){#J+D%P|uY(zw+U}6IYy^kft0ncz{Q_f2|4-`hX?aSI^74*u1p(9iYk3Xyq1wZ zGz2KW%J;}=7y}Jh8lR{cuyl?|s^Lvp3_QpU;sX1ReIHxMUb~irgoRjc;k0x0g&Dwl z)f%oD*izupNi~pY4+#ASHqAhiY3<`*rW`<2dI$z!XdndM+`3_Hn!^qUFp8^38FeCf z!2Em50nXz;$)|fj9FK2nFW>fR6=Ui9l}|ms>^PxyoterO(`QPm_(IV04X(&?Os zU99{0=YN{-X&0rg0iNk0Q@XO=sBhB}3#Va?il&u3EhMfV=iTwr=jvDa;qRdjnBF-= zdYPxd)Af62ugWIheyz7Cc^inMJg$vGwW3xY zy6ST%{~~mj{45_#5y7eM?-p0)+`a}TieRmwTyAulhYyF2i@{vdLzxp@bM)#*fsB)cv-^bxy;Yu5RJh? zCg+Uw`;zZ|081C~B#DC+TS9jq1ubDVm>n!jfR;OXJHaaFin$?3k+&LJ)_}(wYQXwu z-3Dv0^6@-72BeP2w0?WeMe|kv9C;~FK?8COKZrcmD*xRmu9q_{H{$7=!P6nSr#k7`S&1FcAT+03G zY1n(nUD!uL{$Uak`b=<%P4bC7I=$kVC8)@EPa61SanAHuH(Unwt?7h5p!my?P{r$s zJ>I#Oz!7tYU(On8#Om1e4&|5X?}5L$>I2#V%2O8!+G%xr4fuf8rZp+cMG95t1C7_t zw;=^;$0z4gWkp^=nNR-`!7-3&x6&d7ddq2A*Q2Ht<5tB&*Pqp02KA)F3@7S-*kzik zdtKSR^{c>7{Gam2taOiXi_Q9cK!0OA&wxsiURMhvZ#qan)yR1)nkx{9Csw?<2@ zx{3llhCGIq11YLBvSap)zQYWQ8$(XkQCm%3uk0sfMB~Wv2N0`vcpybfV3K%^S8=wYdWIRog1$2Kt*T z>^Ts#Lu_1~Puq8Ze1p3Oo%z&NYS-s^vege8JGI}e2blrsnESwn{O&@J<1nT39AGw( z$-4sEHg$S#-w`uH=}narC)2+%Nu91%mq~(|*fNdR)@}NTGuP=qRIWcR-YNell`10& z4PBjbb!!N7L;n|t5hI}FHeN_FRbFk6nl%VM* z)a9*L`nTr)aNN2UdKZJ)s5N3C!xc1QnL3dzotE_o*g|GNi?zf+18S@v{ zVn^;70M4%59wdDJ`*ub6JL;ZQc!LcXCS&%Lg8 z^AnqG%Qh7KVX+Bb06>%Aq@&O^ubt8HihndAJX19KudFGmeb85{z4l!;FMX=-eRVao z4W5@`;dAe7{;Sm+LN8mCtUv1w_Q*ZWpwvFsEZVr0yGB6U&&a;6&YAl^Z};xSqeDs? zd_K@K^zWTbYB_uK)oOCr>Sg;b?}&mvQ{f5!Dev?qcz-@K#LvIX!Y9`|J-rLG)o)J} z_<*Ka-zI`aLx~>c&%&c-+1TXQ*4CnpdTY;E*nG6L4_*t}+81$2Z82tEIrSw^ZC*Wk zHb{XwPX)j-qvhVD_;=vSMw$Jz8MH4RmwDt2bHP7^S=B>7l7@l61;JRT|PQIzfhLo!o z4?WNBd^ZYuoR{M24q$rx2P@{MGm+3xR;+%1@vDxvswszt(8;=ont^-e-fWV>C zrc6+R|CsaXhLDP=524iOP+-mWPls1vA@`r;Q=fZ*Zqwug%fKBK3s1YOm=3+F_ErDh z#c3%n(De`Hs*(oO+5geS^XWjgca_V`0X@owMLuY%$@=h0RZBU3gBGLz;nH>Z!K9-q zO&qvN)m4kT#sV;=Q>#M@wa3B#pF@fRzJXFBcc3-c=X)MM zGdx)7|Kk|(`SdxEDxq)x0L~+OHeA!X0}AjdkAXZjrEBh#O&nPLU@cd^5O4$Z*tF>? zi%0ApzG~ECsmj@Sj1bhGahb2nxgk3bdG&^n5e%A2X}%7?jk>}qpO5*6ZHxZ)3~zFP z4pr*aC+>iN4pTM{-sV8|GdRax(nfYJ#1eC-O?*~pr0opU{=i&)qHF_;*U853F*t2T zY*lTp$U!qO7RXbir_Y7l+xc|m(X-kLmDUUT8>rMg8=Db=G0>qjVw`IG1a1u+O6#lu zKFNlPGdKKC@(nlySJnrbK@o3T-OukMBf-Dj`_DXR8O5IVL-CV+EsfF3IL828v6Vns<+v6s^tLx3UA|i z&3SuemwJy+r+4eG;RD8I#7qLllr^8M@0lz}+P01yyl&Z?D?98!7XaPYsQ*gU4I%t*n66YxalHct*!V%N&eKM- zNBar#Lp+qlMCGgr*}&E93AJZ717@Vip03RZALHiWp)2+R&aerF+IEiI4*_eJF4F6h zT|7Djt$ok@w-`8bFvvytrUTLxDt` zt?-jZNf{QmEjgq$z9!G)9R+Ab1FtvUnVUvtNl^-MuA(=~1I6B?PX5Dz_ zBMkAe{2(j)s)J2ja3bI;;TW&PIZ4PA#7acY^~}b1K_+B6ud;bZGWiI07Ii6^a4odd zmrLbi-voMBdr;|?mcoV*`>W#vm=0M7S|dZsoh_3aI(U4a=IZfJxIWEe6D8YYzBpep zRmh9)!8MeaWF5PmM9t z>^lqtCaz*ujuQ$mzggJ+Ha3nF((IO8)geE(;CCUqr45tD88oG1qw+hb1P(@QLmRa$1SMJz!FEBag$In+HS&^XpI zm&L=_Uqdk-BosD;gk__w+=GPVT1#UzjV;eD^3FUG@KV5GWYZI~D7-@XPF&?NmKiyQ z??`IiWkqFOm&Y4WoLH4HS!kRwfpwrq%%@lTVTX4&zh{P4ws38PFZGFO!Z>YLOR-Hq4_{90u_@)P}a$Mc7zlY!*3uLYG-3ByX@7Z8Oi}bP7j`6YBE>Ms~T^F;r&%x)qnEk~G z))8E!fsmJziZNmz3QiOq<4wc_i(CR?L}kti^mCja0=S4rO`NzE!&66^ydb0w&sn@JRu1^GL|!}%GUW2kKeie*tJ#+^eW6uyR2ro~!9(L@xoCwa`EkQ zSLdup2{JS9v4zr|tp>Z+SXN0A>YUMuo93Ng5_m^?TfY3hxItXGjgUx=KGI#e2+Gtt zK09N023GETf4ARj|9dY}_CEPOP00Zex2*i?D~wR zU%-tCd{v$-}L^`l4R2D&n8I3sr6f5K)E9~qg?5}w6o`$H;kXMB0!_kq14hx>xUM_lZI zJLUNggu{tW&B-M07VcL07L&J0$(^9LfiV0cc&bh%G<20d-0D2*<>dD6cu{CEt13Zu zDE>MvGHOb4vis=X#F8#!|8iDCLZ|VIW8r;~Q}nXGf3FuJ#U6}d0I3A)Q3Xtz?2seQ z`=--HU%V|)B>Rewh1BM7G^1N$T|PS+--B~HXJX7ryZ{Pd z2}bcH39^QJOo4TD@v^v*AcJx)UpU2*gP9Wa{qfn^Kv8d&zb_36?rrVfeoy`_hg**^ zfS$ka{jiBad5L~VArH@3GAYfZZM(!YXLk83`l0gc7_Tm?*PT@fn8XEL(QlwqdT;t( zs>0%@=m~*xCq_H=y4N*91upibp)WMgnh141=%6?k7F_2CY44@pX6N7I(LnqA7KK4? zatj!hMEo@c>A~J_h3qbWF4D^)PI<+Ky`yKBVDX=o-XQ%_47x&E8HGpRY^mxfDP@-5 zDet|xb)Q1MoW5H&`*>YKA&_OAaXqC@UM32O8*TN{AI(MbQ^bZ$`Qj{C891&m^4gH}-CBQOL<6c;e9ck#PFTvC3dj!3;?wl1@ zk$1^(@*YcWS>fIK7wFtgS)Z{cHKR ziGmx?5rMTBzoEYF`j&E9(mQMt+`1`jDo8;VWY5AVmG#5OfY7MJtQy&!Qcesdakyl( zJ^n`_kfWk97->8)5w+knafx>}G3^$t(~Ik~Yh+XMiPQs#h09+-GJ1*igjDf;?;QM( zQBYGcbfX)PHzu8=wHVDutRYumN3k1So$5aNPF&wXo6?2%zewryC`R6h9}tQ^Op&iz zt9}`{?O|~$lP<*7RLd7T!*VK@-e279I`kr%kjQ7#vxqw+1Spv#m7F9qn?P){s1>-ljYOi;a8z{pzd!MR%ap-lnDz8Dj?m7 z`wII_mE?;cNWxR7oPq~kmalnYsJb*+oi9>dq2v{HYg2-NW=|loJJ772Y63&tY{^$h zTZ+mE3N*>2OkRA1rQZ_K$c#n@$d=t^oy}KpLoVOR@Lf|)isuLIBc+{gmEjna$cDW& ze&e8l}&5SMko`j(|1t_dj44K8_yIXc@!ksWQVyp6HD;q%NbI zaWilwP|?9wG+uN0=Ay=)Vl$2}m!TJ~D1S{q@Y@FM>Vm6+Gt4s9qPBwK-_j4{GvkKY zsaqrS&0=IdvIFJqd@H3NaKjhxB&13F;#a(sz%AU|D{js%6RyS95QaHbM|~HZ3coh5 zUM62eu42q;>|}+ABTqL4iW^fT(eovpi;}5pDkr~{EaC4EzLlLeV?GBd(1@$D{DQBV)Z4{xvu&@oyz{aJN{u-lKya_;O8e zTnvtQI?(pH89i+#g_2qoeOWFq@*{Q#Q}p6CGq!Fu6yHhf%wep_C*h~c&UPHS8}xidpLWkTY>&@gKjn5-% zBsnIQ;M-etN=$qz0qJ`n8Su`wm7hl;uOV%m-_*$pZAP8N*vSGRt-y>dOXBicX2F7F zHgGT_9prT>p!X%Y-JPJs(j-iN%jCE5eHmuNidk0r!NiPbP+PCZVdGJcnS7vhWgMSs zAtHLm%7_^JMcLU_tHw^hn5#EQtAcDnJ!Leya9E48l_>HeSFeuH@~6-}^}>sS{VmGg zGt9;YB5F$)7p#MHH&OtLpVs{N?cP8?)cG26b<*-CLIO>p_I?@!X#lG3g-4r50(8!5vn7RD^GZpkKdD-HHc+pp@+#D&9C;88PI)4Z-L;6#1S-vC zl_oxo1bW;20xu}`fqSVp)S;8{10xtXhpN62Qb>hw?j)sA2n%ccJgAy8g ze~*FmmgRpMAJP+KHTsANu_DkGynN{_?n01c#KY();2Sd5V%9jnXT?8{JHj;(VjJiw zl2I=@Qwn$cgc{>5eu31Do+6|aF+h*F;=w8xq;KslvfseDx-egn5~k^qu|%Xt(!T8- z=^(2Eg3nKr86AXV+%4G~%P5Ih02(XvHD=(B18?E|2TE}h^_Y)Iecu}@k?)J9JN+Wx zR#QadGwBEF5hJ&i|4{D_8-#r)+B^_wlgziWqR*b|dq=vfCP>;wXgSbDfP1w22R~tC z57bqYs?q31UOFw00Afj@zH@2ehx#vs};oa=7e!6VKLh2H`xtIcE|Mzbs!lwtbn+ zV#}1kzBkABPgt^^>HpK27l%jZq)x*S{;{gzi5~-p@!q}^cEGi_--x6Nt&iNHn=~R;MinkIi4+vb zStuzTFXNyT6jTx=X)LLSQALx5h`?5*h@{w2e2_xpU4@HCMl0s?BIU8ZbDb-b42}vL zXrfiCX5U zM`p4(J+h^I6G_tJ{Ls z*dzjlc8FRX>Q%yP+|)H%55QFXuL0kiOaCPqJ8+9S{DbM6-s!a4_-b+rS)Tgoc6IR^ihgy)@g4zLe#RMX|u(s<{ zpNOJm4sG@hCfR5M$I*bqf~zXt2Z$|tE)Z>0ZkdS~#ABRHvx4d(GA9q1l^(!Ut1EGDu;WISU zm+7ilWdTT-IuOY#y#or!2I9#5C!yqk6x9*Y2ILzMkV6+OMp@r~DNf8g5Hf7?*-s0s z;0Dwhgd3~{!GAy?IOLyre0r?Jeo+hi#jYC=Tu2uR35j|3-wD*nIj(SH=704`q3l5X z3OL+P4v_f3a41Mr4J0n0TAl+!@@R*$`}=4dFv3S?D4+jU&_ZC<&5c9LkpZK2>~b55 zMxDB_cAQKVYV~hpk^Wn(8aOhag#Sg3Z0`R&P@|#3r@GGQlcE6uC-QF!HC^)gX%&U+ zZ~til2-yY)+)M`n8Sg4ld;Mo$e-{9i^z)|%d>ni*1YY;Q5Pdb1$B1ZNdhe5>eAV|s zah;%`$<_@633G!7Y*Gh{>KEnf6C%y^AJ&w^A&3s<|4>JWyuqV`4VVZiWAOrPq0y92 z3s9MfMA1zLv6DyWgeCrM0MveQr29>3)Zsf3P5HkbeRZZ{ww;Ys9c7C*B>3j?BbvXgM-Ng;TZOgrR8E-P5YL_)qzlkN=wK zNjIbK_aW!V!Y)J(NnX|cl%X|beP95Zi4Q~pM*kwvzH0uQ*2mp&4GQ<{ds__dK-7W` z-`nqIgP-$j)~>m4xdJQ&j8`fEhUB+D5HAFj)u1fiTff z+2F1Kbu1v~+J_C-polx@+j*Kh5S&Nzhvx31$Eel}mqP_xpcp=Vc^kCuK*XQHA45s| zB|!dw9_Lmq_~z4~KOZ432kdFr!<&mRM_LEuq91S*%zDA%LKSi56L#2>)>lV3wGF+BOzJn(-DylJ@NMBNY@rd#k755dt>fqh22vVaq=Iw7s<98qk`~jr;#C za7m!xp#i9NfIHD5?;U-GB6mR_1ZXr>q$UcfF)g}(oZ=w|vBZ$`KwXTY4O$>S@B@e} zU5q|$*3kdW6ClqI3?}C|Ly@}@`n|>ptt;q><;dCJQjL)NFDk(TTIi?r5^+BF;ndHENfIje%pe|-6_#Kb>=xu#~l&xz~h8*hY z)+K0tMtb&Nf(Q>+w7i(HV`#ZE+B;Nqus8ToC+5(`FZCTF_01DTO(u>)kS`kQ%X(E9 z=DdAo^bpaimU=Q5Wp%UtsU?Tn7aM;j-Ew7g55-j(Elv|K4vy~@zRauj1wYnAcyYr~ zF8d0xn;+AeFO}MG^S>nd=O{@lYS;=0cn2Y0B7+hpfM`z? z{nXv?b&oRKbLLYB;&L_JJ}7pClE2a)%K`&j(g#na4C(zGF5v%hWN_-8wA{q;IKuB! zOF>w0d#x*o4?gLXH*_5bg${(m7WPu_ct`^1`5zi1&9AIG>j|hpAVzmyiU#%q)11K$ z<5^$x^!&I7HX1Ms5}f|TQDuPrvN)$}UJSzgmYClk9jx?f zD$KAfJP%L^>Z|&uew-`BX#uDLfNIg)eS~j85@!1dvlcKWC;kkW7%!mP$i?W}eK-P_X1u#brbte(b)5icIm1QGzr~byPEJnq6 z(*CkPdE;+*L6Hr!_3qCQ%TX)lWbUv7bQ?7gL(RnSnBJe*q#+~?puT7*D+sws19^h8 zATTD{%m%-l3`hY0{YH&oQL`;R0ki?$=0C{?R6(0Uqioc@@U*jIhLofeIzTj+8Rxa> zMaccXnPD_~R98<0Icb{%IT%>4TEo>!p1V@RBvUID!UGzN&Vq1?5R92RixNHLkwbB{ z^Lcl?!*ST4E?k>Wq{5l~9WBO9dgvHG;lK>+_`+$cKzw(*shTSY z^lo$tLOK|YgMbPUT>wf>tG$mddFblZMy3zn^7Gm2fENOVb57BMy#T0-B*_$*dh| zh>dHgFDo^eCJ340w#v#zH-rg>Hl?5P^)Ew6PFm0>7p^BfI_5qJu=@;y0ZlObY5;`^BEZ;~!E`4;KQSCg z)Eb=7(w~6xw}ciWebr11X$_2j(uOl3G(O$rE|u%&@&qM+s0e-t#}B@02v+?^3SBfR zq?H~J0VW52SzY{_zag~zGouE9S%DGj&;G5^!j7p(FrPARSng!-aZWo?y2&@CX};ItT>`P5Zv(_+k} zf+TU;=c&ZkUz6N7pMrkpZKQS8oYowM{)sg7&)?e$qQa|UMnJUbl`lSw(Y-}!o_A9w zJ6-Fkk8>{Q^J>d;29fWVQTOX_r~%|cHSKYM`4$NIaulC?kZLjNX~Y|=ui}2w!nfQz zoA-JE>_y4K?gx7bzmfkdN%DuG8lII(vs_h|E}%}?WqJs0;wu_JJH9YFRD%(nzDV)_A+i#OwHudbNhO2DFn z;7F7|0G$}l*@ge{sjh_TSnWMkF~bky$`xxLdD>{rD%wH-LS_L$`qMig)Cf3O8hEw@ z*e-zZ0Z=bcK=*kzEg+6RIBL?IHX&FBu>JxF2nh88sR#&wZvjyr)W8NL(!K5qRK+2I z^FKKs0Dd3@3OLfevOds^WI>nf5ac_2eb6^GLzrH9wYkt}Uf~abZJ@rYZ=Q`)6I2|w z!vGhE-~;yNfxzBg+o)#1Y=o2m0>Zq5&@SUgVC!6AE&{gBQ~3>)8c72>9-_H;1i*-w zw|o8nf>a46zXU(_x}kLkI0v8d7}^O?m#(=}HW8`xqsJ2%^3Y9!QBFhthVwl4ZRiq* zhNx!?5t7-xi%Bg?2!w|Eq9OUBPt5)aMmgbBXWkP022l335lptgSw?D~8MS?+^rwewb#1hp{43F7Id@^*! zp0BrVr?#CV_d~$iyN5Zw5~{d`L(AHqUsT9xC#0w40fF)7jO*R$#lM`9BrHi0yJ~7<#d#s zhW{b4uRNLA-E#|_Ft%I3Du%evfZSKPFS-kIE&Iy8&-Ni7H>j6CYL>rG|A=g zyZP7)ZsL|hG(Fa>=I8kYx`3U#gGFt-Q!$BMC5wrXN=Fi{--hB8S0`oD=&aa2kzQ}O zA$50#M9L43&B`C2+vAW$$Dxl)1V!h1rBYE<_0c~$owDkXg-&C0B|F-GZ6({Z$VbMN zT*wD}nH`SI3P!CJTawV4?sI%brO;;<^$_L8-u8g*NmM*$1~aBK%QE*!Vz>UfG~q%z zOJ344mP)6Zl-HAafhDMHHlOzapk1;Yu%zP2sTpFbu(0{(d}3l2-LG^4ij;jJO;W=b zcamc()?RqLY%`lL?-DOC=A!T~7gz~T$Jrav$oCC>o<+QL?3ty9OJvQL`?Z(N5|;&L zvD0X|Y0T|=OaPswFy?LPg;X*LrnQW&b|cqe<>ijITTE<-(eXanlI%yJ5kF>arxoDw z!6vbmWU~0QNCb!~lh}29Y3QF7j%meWzRxV-8OfU%oG7L_b>veH6=Rep(0>?KoJ_&?s$EWd5CN6KgPZm3h>srlU*rBNoM&zo_!t!~W zjB5gffT9vC>Z)vhp2ZeVl0V`i*LONN=DNQXlrDMO zeSID=na<>d7rx6LJuh{JxY4JXoUx@c=z>AnLp)r^CK%Yb0?QF!MMu=l}V;A;bTu4v&m34>% zHp`+y%J7+C*2Yu{|1sm|zFWswxKI+w*Tw2$7B0z`@u!n7gW*$3ZMA+u*i$fq$}DYl zWEVUPMa5J)eQ54VM!h75$rgz^1vQ}RA|=OdN#{e?l2j0d3zU=e!96Y4tx52xnrx3B}q=kiv-y?nUJx!8ejf8 zbfVQBAqdgD95Ystq1eJ{<#C9;IWrqirSelYe=#Yiwa-H%)9bj^oae zB(xYCGtC}Bk@PvVSoiqg?*V&6@@?+fQvi3cCVv?qDICNdXagLQS)+rm-LV_pDg1aX+R zS9vA7@BrUmtXKakSV0!dekpXc-Ph{oo>ksiZz2@fOK0yL_NHYSG2bMJYuR)H@KIb3lp`Jc(WSIY^c!7=IWp%e5fw^-a>Qk$}>Y zJmWf2JnDs42K`v%LfpfQOhy-Hs{>{;oZC+}Cn)ZTBSJ7E=a6+FHdOo)N&H zo#u(Fh1-h(K4%MK`&QH)dF3`@Vgs)nXUzv7p6b&)UeV^8jcm@hG_pBNqA2^GVer0p zA4*=i7Z48iVlOUqz615P_%iyg}@O}hbep(r4Avh6#CLz4mm0@3?S+= z0#LAe+D1R*2*4Mp`#j+{Zy>-2(aDjPj0Vt_zL0yp-?Am*CIOHvf`kyHWknivM~Weg z21UGgzp6@r?g(uKkQ5QbXQ;O73WhlqprvTW>K31e4J>9t=M=YOn7RP5KltLUpF_Ii z^Fe-?7yH?hf;I>U#9J5dOc*Xnyhstdh#_q$1NsQR9YN3_g}4#IX|+N?ASs{Lj-_xZ zAMU)8J&_2u;i9jbghB<-ENZ!n92hZL4kbp7hE~@k?1>g?7(Up_i zOPJ{pW1_Jt_S9G*AmKElAf5~FC2Ia&n8Fd`gAak}ArMo9(lz`nl)9ygPe0uhCi%4+ z3$&j(tWL%G=;F+nrR>2a&Dq3^3WTV^F5vfWyfORMt6tgXukXl@)A2!om+ZP@K>LtVeGj(~?6x4>4VYb$CdjT=gLb?xk;Cv)OsEiXXusYmg^ihT8@(JAA_E7(f+@I?N}GEadn8 z+ywwXhstkUWcWN}f&yhzdaZW5)6-4|fmWqCEvSHyH-IW*@Sqh&83_d*6eWpdN>8#B zK&}u0`r_|JK?D^@fE>Ah0Srg^;Eg}B4Lofd&TyXRbi?E_qThkN`|Kx0nS{i^WSy8$ zQgE=6wbB}CwHp<0J0`KD+&QuZg8INb7Ri!HF9;hEluBgWITQf2%yRFd>+LItfjk#b zr$GuDzNlQKtO$!|hLHXypRL8>DVIyiR2UVA<1#h`+vWQ-Rs(`5>2e;~t>(2raFm~6 zRaBySKvRP-x0Ga}D9p2H*cvXpV9a7bh#Mu8rt90$EAP11Y0+aWB{Sv+sa^(dJGB!6 zsX!)87MmGKMvM9&*T4~Tqg22cz-Va$_LIoY1qL@sh?)fzKuBbFLEx4x03^kuAAA`i z`AG1$RJhrRzkF@DLwOnqqqh4 zp<@B*dm`~9nFCel-D?+zR7dVQh{Eh1{Q}KymS^SS&o&q_aGWq=t%4s#dDkk(%W!YS z)E`s5!fAfVOl+pW;CZT7nVdvE+Y}r-+fdP6N4O{Rb&C~BuvYZ71!if>_Lu4hMU@-x zD!*1_|JWZeBRW|NTK{aC)d+D+LNAkinz?{tva39R0H6*Xi4cjUQ)t~<)vPA505D7@ z+0MSg3eEyBFpRPDS_3qs7Vv$tP->@5fuSza#C*_%o%|IO7S}% zut#A2LsuGa@^*w%v$K&ChJV*;Cop4-?U~3}RRV|~;I*q43MIoP>%mYMatD3w9*9cQ z#eDuAz+hCEkz^P%a@-4~^yz|AvbV?68WD_2%<-s&y(7yICJDbU=Prc9KzUoY${J1H zhR8C^Bo%Ls+HDaWP4q~N23G+99-e7|=$u(rTgNibiY1g(2<+ng-Vcij@_S47y#w$e zrE>X+>z#lMaw2LCpo|120CvhGavT-G-K@{=AsmxH;F=S{_I`YxL|*3ketaJ49*tyV zzSRL<+K3%CZ*waGjscPpM_#1&jc7(*L&@t2QFvCi00Iv}fpMrLD=9F?-3v&Q-}c!- z+?uZJM}Rz1^4g)YG&1TUzYW4_6qS3HbOx^nu#K%nT>_p=iO$K7GGri>7JlCgmh3>U zS1buYd03bKlK*%cA^l(g=|=*QekxJ(7?H`XsLi2UNMA%Md(SVvw-8a=U=C?9q$3c= zs5F^zwedSPOxa<)MxFu!M<@DXFhnl$i?b9T2xRds!cFVl9(D0g=QK2_N)q2ZX%_@EVYUehE(_x}lNN`> z7gp0Bv;>MPc9zacrt~gq;-Q%o%nQ4hjIp;`dg5wgBX?U-)*KRCm13d;#nT!9?YXmNpOsL6UKu3puG~$0U0zXAknxOPB0omyR5DfA5K{$Mj5nEZf(v>i z!;D^WAuVloa_57VVDXcbLK;Olth1npg&7;pP)vjcP%eiELs@_@1du=&_wjVl`=JvV zUPVriB{`L#OYFVQI-xA1%8iNH$6 zGL1e!L?sA+V#(;5o&jh{Z+`D9;1RxdqbI(B^fonPt^C0gLdC4S1@Ix4|@zC zKg|fdL`LMuG6sO`1Y1$|o@t>!fxBuV9_t72G%$}x4DooLr3L!_Ae~=X7$>7F=UVC8 z0t?=Z|K`n?ahoTW;p%%U(;)DRZ5j(p<~+c@8=tJ~7KOvYKPMb^nRr;v-TmUq^H#+# zd)IHeaAen8?uT<%&2HUUS`^aJ_~YZdtooy3va^)txB9?lwmD$3%QwE6<~GYJog128YW0L{`ja^h;2 zrGdWFNAmtKoMN!bHa=|mqN7u;hS*Gi&@3Cg)|(DL)$`ta?#jxMzxp&YzrAd_UNFIE zBuldY%ZC4|{EQvEw7mi;YgqOPHzEn%z6Y zOt+fa!B`gv2D9GO$JA;GCL%c86d+sH`t97keK985bSngYnLlS9^q-HZv-2n=KW3h~ z;89lVLNPl${K8GQTCH2U9N=X9qQiYgSuMwaq8uHqmRp<7wT3GIxviE(hp+Gm-yZG< zzJiJaRR&^0J@#$i0Y9+{Tfbtl>3mgkO?}WKh%U3$YSoe@i?@gQ?KHK8Qu(MF1pS5J zXH2bJ0>k}6gP}^Smik5R-r)xoYUO5XV`}B<7aS12%hY>Z+YL3>oDy+`0^gi4eCpyMo#hz{S2O`|PCj;PAg( zH+{W-%0iXZJRX=nzR#edY^At>gELNji^vG^tEcebA-cNJ?&zW=XJu@HvznOWS31y%y&7jpH zcHMLLzqszv0_`$4;+Zj$n&J^{+WWXkxO$i1Tb@0}%5i%ydqUnuh~3z229gj9tI9ld zJEGK6sSrTR6xND_`yHZ&qu*-j@TrG)dF~zEZLe`4Jxt{j@P(03RP-|;ar;;*$n=@? zqiB~KQXhDDO@+@&-`=NSYJ%7!ZZw-|UnpH1Gn0j#IzOdM{k3-2ep6qB8YZ-u9=&OSzzMDz1cP#qcd2xe2VAJ0{xg(%yWl()hzK*R6?;V>a#z3c=H z;e$^tOp*fR`GsXmZiwXW!5O{2aU@r!CLWPg@a()8z}y^i_VjImP8sp&iR?RB&ROuB zVmrRGX@w^ToCcIH0@s9XEG}iU`nsofQC0OM(-PDzkeLs94^LH039ssK!~v06=tKzP zQ0QrCMe&L~6_OQFY=&;MWP2hizue_Fy&>t@`2B6jQut+ekK7$^4bb{*l zSm%PTb;JqlrR z45gNJv)4k$y4pIxSYBi;?o6v+c35?c%c&D)Cv2q^8WijhSb4ZgxV$Vuc3vq>?u}eM z=+<)!Z$a{jbEn5lh>97q<4TQona5m??N$CH=6rlWnR?94E+9V+spYup==R&)GMB6X zv&)F3A`y@0h-;;-cCLAjIaJ5#nEpzr)3AlEqpoRS1vz1v_YX@k)wZ|^(|309%3I1l z%4o-^Vq1(u%>Zq{V&>e8&VEzPeBNGrTwW!~H?mSiJU%|ns?sXm>U!s8ah1baQhhvA z)2mPpibLb*JE=L@IyKBF4#EeI7ncowEDS7M8nhaCJ!mm7R-`#VA7BqKi~Nhm*+#g# zxTs4XF3$Dolq^vX5Exz33#m7CV<$(4qF^koqgqJc1-5i=>QMjCGP|?CxL+@4&GEtS z^I2)d^`=`<>$wMd@79gm8MX;RvTHX8PpzR7Ze>d zAGEQ`&r$s<|Mf7Zp5p+}1IPk+0i^)u{-qez7O@uhn6EKq*AuUMCBhHUBn_a>ng=wx zwTZBXOD6Ek+UKzSo@AC7x)W-n`GycQ{0qE7+Zs{=k%Gpu4!_|2SjUQ0>pghpph*HC zLF%*!`_{b{Jci-guNJ8$pjM}vaZb{u?L47khRN|IkkW{6eCZuAUZTDZzHIbLN$EYc zmyg@mU%M^4jn|Y+Ug@C|eTW^*5wt%${i2LlV{^mHKO-Y^BA;ckQ2l*PkiPPvB|t(u z5KX|Q`J3nm>Xu&7`?xwZ+Z{KZ68YBmtrNO(0^3(MSALzmUNciecAIt2BmTlP zW7ZSx{oVZ;!I4`$f^@-j%Uc~$1w}r^lg2$h9zz~-We3V$;(gYw)>TNYBRzGk!M(8) zqJKEdQ^t`!+QZC2|IVH<2j)Is!UwcOcNEUdVQjcy^p|i2@0Z z>S&GMiRP18I`24g@67`f)K!OEQX4imMcNcn6)Z*iu0gJ;U0c-6N(FSU_J`$FSr=5} z7t|L-)A>kzq??PK;7MEg8pUKK1FI`GW9cXEO{Ci=EjuMWkghzF4TEjg9VXbD5 zJfQdA-*LFLqNuI6fk-E(Ymhu^hL>%+l$Pl$jUM5y&F+2YGtWPK>`Bd=?BqE07Uz`& zA!p0fT@4!ANv!ays4Cp0P~}Q>Rq)jsb*#H>l5^kv?DDmc+4BkF@#byX@M?VFNZf1X zS31?Wh5TY08S4NA3F}hp<5uIOkHm(}WEL{{0{!*PO3cF-S@<637tZg7W=BSuhn11) zJz$xUBaiASd}zz(j1>KpxhH2)NIPc2L+eM@u99xU?wGSxRWo#&)f(bD(k8eRJD>mh zWnExh+1%1AwsV0K(shGa=DElUYUFeSa$@4?cj>xYN8lFP>!B=8ZOy5iqRQ!(kar=| z;c+~47K6o7yRU{R^()DxEqgv~?ivpWWxr|1myimWDM$9i%EulFeXu#d)?yqR{J)2c zO1jUcJAtn%Cs%aFZA^DCxnE(PF?;0XuK|_Jj=X~R!Ci|RGHe59irXIn!~9RLWs7E< z`IwPD{c|OVSvk zTDJaq-RCK1PU$3*LhKnmN0Fd2Kcx59$8X)=kP6ZITJAmrGB`P7p!dnMIZ_`zTCQY> z`L-JJWJ>lksx7M9Y~|>PWrreN32f-G8Sv%Q(5-s~ZHC9cxoy4pFznwL@CB+Nt`Gau zZS@Pg?7@Y?##i6^7ch&|SF2x`g~}0 zX`vv|HY0>*BWh*q{|q+-^4*65!H$7`TnlPWM3cO6L>~}%+qeEZL*j_WL<;U70SGq3 zVX+u8hKR;&TXQWKn1CbFC?H)p_%e+|#o#emDiETj1@@)bAKjNm^g)qn1U!mH1%XLEBqHV#w;F0EGiBSly957q>N{rtA1D5I z2X;^Z>_GG4eiFl7#kL6H`f(J2JK^HNPGA0xL16G#FjO2G#U+J*+P2eD`~R><0QYqt zPy+T(h4hd3x)Q~T$9+)Zu&So4QC7z_s2`B$ToQ8+w? z4CL|@*Z~)W;r0atK1IZmK-;&XJqT=1B5gC>_VEtfpvcsn4sz4L;Xnlib6X4GUmF=J A(*OVf diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R new file mode 100644 index 0000000..ba25b12 --- /dev/null +++ b/tests/testthat/setup.R @@ -0,0 +1,3 @@ +# Route implicit graphics output to an in-memory/null PDF device during tests. +# Without this, heatmap tests can create tests/testthat/Rplots.pdf. +options(device = function(...) grDevices::pdf(file = NULL)) From b3c801a6aeb640c3abddca2ef325c6bd8effd38e Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Fri, 3 Apr 2026 23:03:52 +0100 Subject: [PATCH 06/10] Rcpp and faster R implementations --- .gitignore | 4 + DESCRIPTION | 3 + R/RcppExports.R | 11 ++ R/similarity.R | 256 +++++++++++++++++++++++++---- man/jaccard_sim.Rd | 22 ++- src/RcppExports.cpp | 48 ++++++ src/weighted_jaccard.cpp | 272 +++++++++++++++++++++++++++++++ tests/testthat/test-similarity.R | 38 +++++ 8 files changed, 620 insertions(+), 34 deletions(-) create mode 100644 R/RcppExports.R create mode 100644 src/RcppExports.cpp create mode 100644 src/weighted_jaccard.cpp diff --git a/.gitignore b/.gitignore index 4de88c5..c198f33 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,7 @@ docs inst/doc samples/ tests/testthat/Rplots.pdf +src/*.o +src/*.so +src/*.dll +src/*.dylib diff --git a/DESCRIPTION b/DESCRIPTION index 9ae0c1a..8e3f6e9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -20,8 +20,11 @@ Imports: glue, Matrix (>= 1.5-3), methods, + Rcpp, grDevices, stats +LinkingTo: + Rcpp Suggests: covr, spelling, diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..c20f3fe --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,11 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +weighted_jaccard_sparse_fill_cpp <- function(x, pattern, transpose = FALSE) { + .Call(`_coconat_weighted_jaccard_sparse_fill_cpp`, x, pattern, transpose) +} + +weighted_jaccard_dense_cpp <- function(x, transpose = FALSE) { + .Call(`_coconat_weighted_jaccard_dense_cpp`, x, transpose) +} + diff --git a/R/similarity.R b/R/similarity.R index 01dd78d..c5f21ef 100644 --- a/R/similarity.R +++ b/R/similarity.R @@ -47,14 +47,21 @@ connectivity_similarity <- function(x, metric = c("cosine", "jaccard", "weighted #' #' @details Both variants are optimised for sparse matrices. The binary variant #' uses \code{Matrix::crossprod} on the binarised matrix for efficient -#' intersection/union computation. The weighted variant uses the identity -#' \code{min(a,b) = (a + b - |a - b|) / 2} to leverage sparse matrix -#' arithmetic. +#' intersection/union computation. The weighted variant accumulates +#' \code{sum(min(a,b))} across shared features for each pair; compiled C++ +#' backends are available for both dense and sparse output. #' #' @param x A data matrix suitable for clustering (non-negative values expected) #' @param weighted If \code{FALSE} (the default), compute binary Jaccard #' similarity. If \code{TRUE}, compute weighted (generalised) Jaccard #' similarity. +#' @param weighted_method Algorithm to use when \code{weighted=TRUE}. +#' \code{"auto"} (the default) selects the best available backend: compiled +#' C++ dense for moderate output sizes, compiled C++ sparse when the output +#' exceeds \code{10000 x 10000}, or pure R dense as a fallback when compiled +#' code is unavailable. Explicit choices: \code{"cpp_dense"} and +#' \code{"cpp_sparse"} for the compiled backends, \code{"dense"} and +#' \code{"sparse"} for the pure R implementations. #' @param sparse Whether to return a sparse matrix (default \code{FALSE}) #' @param transpose When \code{FALSE} (the default) calculates similarity #' between columns. When \code{TRUE} calculates similarity between rows. @@ -70,12 +77,15 @@ connectivity_similarity <- function(x, metric = c("cosine", "jaccard", "weighted #' jaccard_sim(am) #' # Weighted Jaccard #' jaccard_sim(am, weighted=TRUE) -jaccard_sim <- function(x, weighted = FALSE, sparse = FALSE, transpose = FALSE) { +jaccard_sim <- function(x, weighted = FALSE, sparse = FALSE, transpose = FALSE, + weighted_method = c("auto", "cpp_dense", "cpp_sparse", + "dense", "sparse")) { cx <- class(x) if (!is.matrix(x) && !isTRUE(attr(cx, "package") == "Matrix")) stop("I don't recognise that as a matrix!") if (!inherits(x, "dgCMatrix")) x <- as(x, "dgCMatrix") + weighted_method <- match.arg(weighted_method) crossfun <- if(transpose) Matrix::tcrossprod else Matrix::crossprod n <- if(transpose) nrow(x) else ncol(x) nms <- if(transpose) rownames(x) else colnames(x) @@ -94,39 +104,225 @@ jaccard_sim <- function(x, weighted = FALSE, sparse = FALSE, transpose = FALSE) Matrix::diag(A) <- 1 sim <- A } else { - # Weighted Jaccard: sum(min(a,b)) / sum(max(a,b)) - # For non-negative values: min(a,b) = sum_{t=1}^{max} I(a>=t)*I(b>=t) - # Each threshold gives a sparse crossprod (same fast op as cosine). - max_val <- if (length(x@x)) max(x@x) else 0 - if (max_val == 0) { - sim <- Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, - dims = c(n, n), dimnames = list(nms, nms)) - } else { - cs <- if(transpose) Matrix::rowSums(x) else Matrix::colSums(x) - vals <- x@x - p <- x@p - ri <- x@i - col_idx <- rep(seq_len(ncol(x)), diff(p)) - min_sums <- matrix(0, n, n) - for (t in seq_len(max_val)) { - keep <- vals >= t - bt <- Matrix::sparseMatrix( - i = ri[keep] + 1L, j = col_idx[keep], - x = 1, dims = x@Dim, dimnames = x@Dimnames) - min_sums <- min_sums + as.matrix(crossfun(bt)) + if (weighted_method == "auto") { + has_cpp <- is.function(tryCatch(weighted_jaccard_dense_cpp, + error = function(e) NULL)) + ncomp <- if (transpose) nrow(x) else ncol(x) + if (has_cpp) { + weighted_method <- if (ncomp > 10000L) "cpp_sparse" else "cpp_dense" + } else { + weighted_method <- "dense" } - max_sums <- outer(cs, cs, "+") - min_sums - sim <- min_sums / max_sums - sim[is.nan(sim)] <- 0 - diag(sim) <- 1 - dimnames(sim) <- list(nms, nms) - sim <- Matrix::Matrix(sim, sparse = TRUE) } + sim <- switch(weighted_method, + cpp_dense = jaccard_sim_weighted_cpp_dense(x, sparse = sparse, transpose = transpose), + cpp_sparse = jaccard_sim_weighted_cpp_sparse(x, sparse = sparse, transpose = transpose), + dense = jaccard_sim_weighted_dense_r(x, sparse = sparse, transpose = transpose), + sparse = jaccard_sim_weighted_sparse_r(x, sparse = TRUE, transpose = transpose) + ) } if (sparse) sim else as.matrix(sim) } +jaccard_weighted_feature_view <- function(x, transpose = FALSE) { + if (transpose) { + feat_idx <- rep.int(seq_len(ncol(x)), diff(x@p)) + items <- x@i + 1L + nfeat <- ncol(x) + totals <- Matrix::rowSums(x) + } else { + feat_idx <- x@i + 1L + items <- rep.int(seq_len(ncol(x)), diff(x@p)) + nfeat <- nrow(x) + totals <- Matrix::colSums(x) + } + + vals <- x@x + counts <- tabulate(feat_idx, nbins = nfeat) + offsets <- cumsum(c(0L, counts)) + fill <- offsets[-length(offsets)] + 1L + feat_items <- integer(length(vals)) + feat_vals <- numeric(length(vals)) + + for (k in seq_along(vals)) { + f <- feat_idx[k] + dest <- fill[f] + feat_items[dest] <- items[k] + feat_vals[dest] <- vals[k] + fill[f] <- dest + 1L + } + + list( + ncomp = if (transpose) nrow(x) else ncol(x), + nfeat = nfeat, + totals = totals, + offsets = offsets, + items = feat_items, + vals = feat_vals + ) +} + +jaccard_sim_weighted_cpp_sparse <- function(x, sparse = TRUE, transpose = FALSE) { + crossfun <- if (transpose) Matrix::tcrossprod else Matrix::crossprod + n <- if (transpose) nrow(x) else ncol(x) + nms <- if (transpose) rownames(x) else colnames(x) + if (length(x@x) == 0L) { + sim <- Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, + dims = c(n, n), dimnames = list(nms, nms)) + return(if (sparse) sim else as.matrix(sim)) + } + + # Get sparsity pattern via crossprod of binarised matrix + b <- x + b@x <- rep(1, length(b@x)) + A <- as(crossfun(b), "generalMatrix") + + # Column sums for denominator + totals <- if (transpose) Matrix::rowSums(x) else Matrix::colSums(x) + + # Fill in min_sums using C++ + A@x <- weighted_jaccard_sparse_fill_cpp(x, A, transpose = transpose) + + # Convert min_sums to similarity: sim = ms / (s_i + s_j - ms) + col_idx <- rep(seq_along(diff(A@p)), diff(A@p)) + row_idx <- A@i + 1L + denom <- totals[row_idx] + totals[col_idx] - A@x + nonzero <- denom != 0 + A@x[nonzero] <- A@x[nonzero] / denom[nonzero] + A@x[!nonzero] <- 0 + Matrix::diag(A) <- 1 + + if (sparse) A else as.matrix(A) +} + +jaccard_sim_weighted_cpp_dense <- function(x, sparse = FALSE, transpose = FALSE) { + n <- if (transpose) nrow(x) else ncol(x) + nms <- if (transpose) rownames(x) else colnames(x) + if (length(x@x) == 0L) { + sim <- diag(n) + dimnames(sim) <- list(nms, nms) + if (sparse) return(Matrix::Matrix(sim, sparse = TRUE)) + return(sim) + } + sim <- weighted_jaccard_dense_cpp(x, transpose = transpose) + dimnames(sim) <- list(nms, nms) + if (sparse) Matrix::Matrix(sim, sparse = TRUE) else sim +} + +jaccard_sim_weighted_sparse_r <- function(x, sparse = TRUE, transpose = FALSE) { + warning("Using slow pure R sparse weighted Jaccard. ", + "Install the natcpp package for much faster compiled code.", + call. = FALSE) + n <- if (transpose) nrow(x) else ncol(x) + nms <- if (transpose) rownames(x) else colnames(x) + if (length(x@x) == 0L) { + sim <- Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, + dims = c(n, n), dimnames = list(nms, nms)) + return(if (sparse) sim else as.matrix(sim)) + } + + fv <- jaccard_weighted_feature_view(x, transpose = transpose) + counts <- diff(fv$offsets) + npairs <- sum((counts * pmax.int(counts - 1L, 0L)) %/% 2L) + + ii <- integer(npairs) + jj <- integer(npairs) + mins <- numeric(npairs) + pos <- 1L + + for (f in seq_len(fv$nfeat)) { + start <- fv$offsets[f] + 1L + end <- fv$offsets[f + 1L] + k <- end - start + 1L + if (k <= 1L) next + + if (k == 2L) { + va <- fv$vals[start] + vb <- fv$vals[end] + ii[pos] <- fv$items[start] + jj[pos] <- fv$items[end] + mins[pos] <- if (va < vb) va else vb + pos <- pos + 1L + next + } + + for (a in start:(end - 1L)) { + ca <- fv$items[a] + va <- fv$vals[a] + for (b in (a + 1L):end) { + vb <- fv$vals[b] + ii[pos] <- ca + jj[pos] <- fv$items[b] + mins[pos] <- if (va < vb) va else vb + pos <- pos + 1L + } + } + } + + if (npairs == 0L) { + sim <- Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, + dims = c(n, n), dimnames = list(nms, nms)) + return(if (sparse) sim else as.matrix(sim)) + } + + offdiag <- ii != jj + min_sums <- Matrix::sparseMatrix( + i = c(ii, jj[offdiag]), + j = c(jj, ii[offdiag]), + x = c(mins, mins[offdiag]), + dims = c(n, n), + dimnames = list(nms, nms) + ) + + cs <- fv$totals + col_idx <- rep(seq_len(ncol(min_sums)), diff(min_sums@p)) + row_idx <- min_sums@i + 1L + denom <- cs[row_idx] + cs[col_idx] - min_sums@x + nz <- denom != 0 + min_sums@x[nz] <- min_sums@x[nz] / denom[nz] + min_sums@x[!nz] <- 0 + Matrix::diag(min_sums) <- 1 + + if (sparse) min_sums else as.matrix(min_sums) +} + +jaccard_sim_weighted_dense_r <- function(x, sparse = FALSE, transpose = FALSE) { + n <- if (transpose) nrow(x) else ncol(x) + nms <- if (transpose) rownames(x) else colnames(x) + if (length(x@x) == 0L) { + sim <- diag(1, n) + dimnames(sim) <- list(nms, nms) + return(if (sparse) Matrix::Matrix(sim, sparse = TRUE) else sim) + } + + fv <- jaccard_weighted_feature_view(x, transpose = transpose) + ncomp <- fv$ncomp + out <- matrix(0, ncomp, ncomp) + + for (f in seq_len(fv$nfeat)) { + start <- fv$offsets[f] + 1L + end <- fv$offsets[f + 1L] + k <- end - start + 1L + if (k <= 1L) next + items_f <- fv$items[start:end] + vals_f <- fv$vals[start:end] + out[items_f, items_f] <- out[items_f, items_f] + outer(vals_f, vals_f, pmin) + } + + # Convert min_sums to similarity column by column to avoid n*n temporary + totals <- fv$totals + for (j in seq_len(ncomp)) { + denom <- totals + totals[j] - out[, j] + nz <- denom != 0 + out[nz, j] <- out[nz, j] / denom[nz] + out[!nz, j] <- 0 + out[j, j] <- 1 + } + dimnames(out) <- list(nms, nms) + + if (sparse) Matrix::Matrix(out, sparse = TRUE) else out +} #' Tanimoto (extended Jaccard) similarity for sparse or dense matrices #' diff --git a/man/jaccard_sim.Rd b/man/jaccard_sim.Rd index 891889d..c73af94 100644 --- a/man/jaccard_sim.Rd +++ b/man/jaccard_sim.Rd @@ -4,7 +4,13 @@ \alias{jaccard_sim} \title{Jaccard similarity for sparse or dense matrices} \usage{ -jaccard_sim(x, weighted = FALSE, sparse = FALSE, transpose = FALSE) +jaccard_sim( + x, + weighted = FALSE, + sparse = FALSE, + transpose = FALSE, + weighted_method = c("auto", "cpp_dense", "cpp_sparse", "dense", "sparse") +) } \arguments{ \item{x}{A data matrix suitable for clustering (non-negative values expected)} @@ -17,6 +23,14 @@ similarity.} \item{transpose}{When \code{FALSE} (the default) calculates similarity between columns. When \code{TRUE} calculates similarity between rows.} + +\item{weighted_method}{Algorithm to use when \code{weighted=TRUE}. +\code{"auto"} (the default) selects the best available backend: compiled +C++ dense for moderate output sizes, compiled C++ sparse when the output +exceeds \code{10000 x 10000}, or pure R dense as a fallback when compiled +code is unavailable. Explicit choices: \code{"cpp_dense"} and +\code{"cpp_sparse"} for the compiled backends, \code{"dense"} and +\code{"sparse"} for the pure R implementations.} } \value{ A square similarity matrix with values in \code{[0,1]}. @@ -30,9 +44,9 @@ Computes pairwise Jaccard similarity between columns (or rows) \details{ Both variants are optimised for sparse matrices. The binary variant uses \code{Matrix::crossprod} on the binarised matrix for efficient - intersection/union computation. The weighted variant uses the identity - \code{min(a,b) = (a + b - |a - b|) / 2} to leverage sparse matrix - arithmetic. + intersection/union computation. The weighted variant accumulates + \code{sum(min(a,b))} across shared features for each pair; compiled C++ + backends are available for both dense and sparse output. } \examples{ da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..fb0ccbd --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,48 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include + +using namespace Rcpp; + +#ifdef RCPP_USE_GLOBAL_ROSTREAM +Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); +Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); +#endif + +// weighted_jaccard_sparse_fill_cpp +NumericVector weighted_jaccard_sparse_fill_cpp(const S4& x, const S4& pattern, bool transpose); +RcppExport SEXP _coconat_weighted_jaccard_sparse_fill_cpp(SEXP xSEXP, SEXP patternSEXP, SEXP transposeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const S4& >::type x(xSEXP); + Rcpp::traits::input_parameter< const S4& >::type pattern(patternSEXP); + Rcpp::traits::input_parameter< bool >::type transpose(transposeSEXP); + rcpp_result_gen = Rcpp::wrap(weighted_jaccard_sparse_fill_cpp(x, pattern, transpose)); + return rcpp_result_gen; +END_RCPP +} +// weighted_jaccard_dense_cpp +NumericMatrix weighted_jaccard_dense_cpp(const S4& x, bool transpose); +RcppExport SEXP _coconat_weighted_jaccard_dense_cpp(SEXP xSEXP, SEXP transposeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const S4& >::type x(xSEXP); + Rcpp::traits::input_parameter< bool >::type transpose(transposeSEXP); + rcpp_result_gen = Rcpp::wrap(weighted_jaccard_dense_cpp(x, transpose)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_coconat_weighted_jaccard_sparse_fill_cpp", (DL_FUNC) &_coconat_weighted_jaccard_sparse_fill_cpp, 3}, + {"_coconat_weighted_jaccard_dense_cpp", (DL_FUNC) &_coconat_weighted_jaccard_dense_cpp, 2}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_coconat(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp new file mode 100644 index 0000000..853621f --- /dev/null +++ b/src/weighted_jaccard.cpp @@ -0,0 +1,272 @@ +#include +#include +#include + +using namespace Rcpp; + +// Sparse accumulation: given the sparsity pattern from crossprod(binarise(x)) +// as a dgCMatrix, fill in min_sums by iterating one output column at a time. +// For each output column cb, populate a row->position lookup, then iterate +// over all features where cb is nonzero, accumulating min(x[f,ca], x[f,cb]) +// for each co-occurring item ca. Writes to out[] are contiguous per column. +// [[Rcpp::export]] +NumericVector weighted_jaccard_sparse_fill_cpp( + const S4& x, const S4& pattern, bool transpose = false) { + + // Original matrix slots (dgCMatrix: nr x nc) + IntegerVector dims = x.slot("Dim"); + IntegerVector xi = x.slot("i"); + IntegerVector xp = x.slot("p"); + NumericVector xx = x.slot("x"); + + const int nr = dims[0]; + const int nc = dims[1]; + const int nfeat = transpose ? nc : nr; + + // Pattern matrix slots (dgCMatrix: ncomp x ncomp) + IntegerVector Ai = pattern.slot("i"); + IntegerVector Ap = pattern.slot("p"); + IntegerVector Adims = pattern.slot("Dim"); + const int ncomp = Adims[0]; + const int annz = Ai.size(); + + // Output: min_sums values aligned with pattern's (i, p) structure + NumericVector out(annz, 0.0); + + // Build CSR for feature dimension: for each feature f, the items (columns + // or rows being compared) that are nonzero in that feature, with their values. + std::vector feat_count(nfeat, 0); + if (!transpose) { + for (int col = 0; col < nc; ++col) + for (int idx = xp[col]; idx < xp[col + 1]; ++idx) + feat_count[xi[idx]]++; + } else { + for (int col = 0; col < nc; ++col) + feat_count[col] = xp[col + 1] - xp[col]; + } + + std::vector feat_offsets(nfeat + 1, 0); + for (int f = 0; f < nfeat; ++f) + feat_offsets[f + 1] = feat_offsets[f] + feat_count[f]; + + const int total_nnz = feat_offsets[nfeat]; + std::vector feat_items(total_nnz); // which item (col/row being compared) + std::vector feat_vals(total_nnz); // value of x at that position + std::vector feat_fill(feat_offsets.begin(), feat_offsets.end()); + + if (!transpose) { + for (int col = 0; col < nc; ++col) { + for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { + const int row = xi[idx]; + const int dest = feat_fill[row]++; + feat_items[dest] = col; + feat_vals[dest] = xx[idx]; + } + } + } else { + for (int col = 0; col < nc; ++col) { + for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { + const int dest = feat_fill[col]++; + feat_items[dest] = xi[idx]; + feat_vals[dest] = xx[idx]; + } + } + } + + // Also need item->features CSR: for each compared item, which features + // it participates in and with what value. + std::vector item_count(ncomp, 0); + for (int f = 0; f < nfeat; ++f) { + for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) + item_count[feat_items[idx]]++; + } + + std::vector item_offsets(ncomp + 1, 0); + for (int c = 0; c < ncomp; ++c) + item_offsets[c + 1] = item_offsets[c] + item_count[c]; + + std::vector item_feats(total_nnz); // which feature + std::vector item_vals(total_nnz); // value + std::vector item_fill(item_offsets.begin(), item_offsets.end()); + + for (int f = 0; f < nfeat; ++f) { + for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) { + const int item = feat_items[idx]; + const int dest = item_fill[item]++; + item_feats[dest] = f; + item_vals[dest] = feat_vals[idx]; + } + } + + // row_to_pos[row] = position in out[] for entry (row, current_col). + // Populated once per output column, then cleared. + std::vector row_to_pos(ncomp, -1); + + // Iterate over output columns + for (int cb = 0; cb < ncomp; ++cb) { + // Populate row_to_pos for pattern column cb + for (int idx = Ap[cb]; idx < Ap[cb + 1]; ++idx) + row_to_pos[Ai[idx]] = idx; + + // For each feature f where cb is nonzero + for (int fi = item_offsets[cb]; fi < item_offsets[cb + 1]; ++fi) { + const int f = item_feats[fi]; + const double vb = item_vals[fi]; + + // For each other item ca also in this feature, accumulate min into (ca, cb) + for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) { + const int ca = feat_items[idx]; + const int pos = row_to_pos[ca]; + if (pos >= 0) + out[pos] += std::min(feat_vals[idx], vb); + } + } + + // Clear row_to_pos + for (int idx = Ap[cb]; idx < Ap[cb + 1]; ++idx) + row_to_pos[Ai[idx]] = -1; + } + + return out; +} + +// Dense accumulation with adaptive strategy: +// - Small ncomp (output fits in cache): feature-oriented loop, no extra setup +// - Large ncomp: column-oriented loop with item->features CSR for stride-1 writes +// [[Rcpp::export]] +NumericMatrix weighted_jaccard_dense_cpp(const S4& x, bool transpose = false) { + IntegerVector dims = x.slot("Dim"); + IntegerVector xi = x.slot("i"); + IntegerVector xp = x.slot("p"); + NumericVector xx = x.slot("x"); + + const int nr = dims[0]; + const int nc = dims[1]; + const int ncomp = transpose ? nr : nc; + const int nfeat = transpose ? nc : nr; + + // Totals per compared item + std::vector totals(ncomp, 0.0); + + // Build feature CSR: for each feature f, which items are nonzero + std::vector feat_count(nfeat, 0); + if (!transpose) { + for (int col = 0; col < nc; ++col) { + for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { + totals[col] += xx[idx]; + feat_count[xi[idx]]++; + } + } + } else { + for (int col = 0; col < nc; ++col) { + feat_count[col] = xp[col + 1] - xp[col]; + for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { + totals[xi[idx]] += xx[idx]; + } + } + } + + std::vector feat_offsets(nfeat + 1, 0); + for (int f = 0; f < nfeat; ++f) + feat_offsets[f + 1] = feat_offsets[f] + feat_count[f]; + + const int total_nnz = feat_offsets[nfeat]; + std::vector feat_items(total_nnz); + std::vector feat_vals(total_nnz); + std::vector feat_fill(feat_offsets.begin(), feat_offsets.end()); + + if (!transpose) { + for (int col = 0; col < nc; ++col) { + for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { + const int row = xi[idx]; + const int dest = feat_fill[row]++; + feat_items[dest] = col; + feat_vals[dest] = xx[idx]; + } + } + } else { + for (int col = 0; col < nc; ++col) { + for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { + const int dest = feat_fill[col]++; + feat_items[dest] = xi[idx]; + feat_vals[dest] = xx[idx]; + } + } + } + + NumericMatrix out(ncomp, ncomp); + + // Threshold: ncomp^2 * 8 bytes > ~12MB (typical L3 share per core) + const bool use_colwise = (static_cast(ncomp) * ncomp * 8 > 12 * 1024 * 1024); + + if (use_colwise) { + // Column-oriented: build item->features CSR, iterate output columns + std::vector item_count(ncomp, 0); + for (int f = 0; f < nfeat; ++f) + for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) + item_count[feat_items[idx]]++; + + std::vector item_offsets(ncomp + 1, 0); + for (int c = 0; c < ncomp; ++c) + item_offsets[c + 1] = item_offsets[c] + item_count[c]; + + std::vector item_feats(total_nnz); + std::vector item_vals(total_nnz); + std::vector item_fill(item_offsets.begin(), item_offsets.end()); + + for (int f = 0; f < nfeat; ++f) { + for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) { + const int item = feat_items[idx]; + const int dest = item_fill[item]++; + item_feats[dest] = f; + item_vals[dest] = feat_vals[idx]; + } + } + + for (int cb = 0; cb < ncomp; ++cb) { + double* col = &out[static_cast(cb) * ncomp]; + for (int fi = item_offsets[cb]; fi < item_offsets[cb + 1]; ++fi) { + const int f = item_feats[fi]; + const double vb = item_vals[fi]; + for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) { + col[feat_items[idx]] += std::min(feat_vals[idx], vb); + } + } + } + } else { + // Feature-oriented: iterate features, scatter into output + for (int f = 0; f < nfeat; ++f) { + const int start = feat_offsets[f]; + const int end = feat_offsets[f + 1]; + const int k = end - start; + + for (int a = 0; a < k; ++a) { + const int ca = feat_items[start + a]; + const double va = feat_vals[start + a]; + out[static_cast(ca) * ncomp + ca] += va; + for (int b = a + 1; b < k; ++b) { + const int cb = feat_items[start + b]; + const double mv = std::min(va, feat_vals[start + b]); + out[static_cast(ca) * ncomp + cb] += mv; + out[static_cast(cb) * ncomp + ca] += mv; + } + } + } + } + + // Convert min_sums to similarity + for (int cb = 0; cb < ncomp; ++cb) { + for (int ca = 0; ca < ncomp; ++ca) { + if (ca == cb) { + out[static_cast(cb) * ncomp + ca] = 1.0; + } else { + const std::size_t pos = static_cast(cb) * ncomp + ca; + const double ms = out[pos]; + const double denom = totals[ca] + totals[cb] - ms; + out[pos] = (denom > 0.0) ? ms / denom : 0.0; + } + } + } + + return out; +} diff --git a/tests/testthat/test-similarity.R b/tests/testthat/test-similarity.R index 30f3763..8c9ec99 100644 --- a/tests/testthat/test-similarity.R +++ b/tests/testthat/test-similarity.R @@ -139,6 +139,44 @@ test_that("jaccard_sim sparse output parameter works", { expect_equal(as.matrix(jws), jaccard_sim(m, weighted = TRUE, sparse = FALSE)) }) +test_that("all weighted jaccard methods agree", { + m_sparse <- Matrix::Matrix(c(4,0,1, 2,3,1, 0,3,0), nrow = 3, ncol = 3, sparse = TRUE) + colnames(m_sparse) <- c("n1", "n2", "n3") + rownames(m_sparse) <- c("p1", "p2", "p3") + m_dense <- as.matrix(m_sparse) + + ref <- jaccard_sim(m_sparse, weighted = TRUE, weighted_method = "dense") + ref_t <- jaccard_sim(m_sparse, weighted = TRUE, transpose = TRUE, + weighted_method = "dense") + + for (method in c("sparse", "dense", "cpp_dense", "cpp_sparse")) { + suppressWarnings({ + # Default (dense return, transpose=FALSE) + expect_equal( + jaccard_sim(m_sparse, weighted = TRUE, weighted_method = method), + ref, info = paste(method, "sparse input") + ) + # Dense matrix input + expect_equal( + jaccard_sim(m_dense, weighted = TRUE, weighted_method = method), + ref, info = paste(method, "dense input") + ) + # transpose=TRUE + expect_equal( + jaccard_sim(m_sparse, weighted = TRUE, transpose = TRUE, + weighted_method = method), + ref_t, info = paste(method, "transpose") + ) + # sparse return + expect_equal( + as.matrix(jaccard_sim(m_sparse, weighted = TRUE, sparse = TRUE, + weighted_method = method)), + ref, info = paste(method, "sparse return") + ) + }) + } +}) + test_that("jaccard_sim transpose parameter works", { m <- Matrix::Matrix(c(4,0,1, 2,3,1, 0,3,0), nrow = 3, ncol = 3, sparse = TRUE) colnames(m) <- c("n1", "n2", "n3") From c7e2c6733600b6a08ab4ca7551441335cb8f1107 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Sat, 4 Apr 2026 19:12:26 +0100 Subject: [PATCH 07/10] move cpp methods to natcpp --- DESCRIPTION | 6 +- NAMESPACE | 1 + R/RcppExports.R | 11 -- R/similarity.R | 37 +---- src/RcppExports.cpp | 48 ------ src/weighted_jaccard.cpp | 272 ------------------------------- tests/testthat/test-similarity.R | 6 +- 7 files changed, 14 insertions(+), 367 deletions(-) delete mode 100644 R/RcppExports.R delete mode 100644 src/RcppExports.cpp delete mode 100644 src/weighted_jaccard.cpp diff --git a/DESCRIPTION b/DESCRIPTION index 8e3f6e9..3672d37 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -20,13 +20,11 @@ Imports: glue, Matrix (>= 1.5-3), methods, - Rcpp, grDevices, stats -LinkingTo: - Rcpp -Suggests: +Suggests: covr, + natcpp (>= 0.3.0), spelling, testthat (>= 3.0.0), ComplexHeatmap, diff --git a/NAMESPACE b/NAMESPACE index c95d85f..0d9a8f8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,6 +13,7 @@ export(prepare_cosine_matrix) export(prepare_similarity_matrix) export(register_dataset) export(rowScaleM) +export(sim2dist) export(tanimoto_sim) importFrom(grDevices,hcl.colors) importFrom(methods,as) diff --git a/R/RcppExports.R b/R/RcppExports.R deleted file mode 100644 index c20f3fe..0000000 --- a/R/RcppExports.R +++ /dev/null @@ -1,11 +0,0 @@ -# Generated by using Rcpp::compileAttributes() -> do not edit by hand -# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 - -weighted_jaccard_sparse_fill_cpp <- function(x, pattern, transpose = FALSE) { - .Call(`_coconat_weighted_jaccard_sparse_fill_cpp`, x, pattern, transpose) -} - -weighted_jaccard_dense_cpp <- function(x, transpose = FALSE) { - .Call(`_coconat_weighted_jaccard_dense_cpp`, x, transpose) -} - diff --git a/R/similarity.R b/R/similarity.R index c5f21ef..17d091e 100644 --- a/R/similarity.R +++ b/R/similarity.R @@ -105,8 +105,8 @@ jaccard_sim <- function(x, weighted = FALSE, sparse = FALSE, transpose = FALSE, sim <- A } else { if (weighted_method == "auto") { - has_cpp <- is.function(tryCatch(weighted_jaccard_dense_cpp, - error = function(e) NULL)) + has_cpp <- requireNamespace("natcpp", quietly = TRUE) && + utils::packageVersion("natcpp") >= "0.3.0" ncomp <- if (transpose) nrow(x) else ncol(x) if (has_cpp) { weighted_method <- if (ncomp > 10000L) "cpp_sparse" else "cpp_dense" @@ -164,36 +164,11 @@ jaccard_weighted_feature_view <- function(x, transpose = FALSE) { } jaccard_sim_weighted_cpp_sparse <- function(x, sparse = TRUE, transpose = FALSE) { - crossfun <- if (transpose) Matrix::tcrossprod else Matrix::crossprod n <- if (transpose) nrow(x) else ncol(x) nms <- if (transpose) rownames(x) else colnames(x) - if (length(x@x) == 0L) { - sim <- Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, - dims = c(n, n), dimnames = list(nms, nms)) - return(if (sparse) sim else as.matrix(sim)) - } - - # Get sparsity pattern via crossprod of binarised matrix - b <- x - b@x <- rep(1, length(b@x)) - A <- as(crossfun(b), "generalMatrix") - - # Column sums for denominator - totals <- if (transpose) Matrix::rowSums(x) else Matrix::colSums(x) - - # Fill in min_sums using C++ - A@x <- weighted_jaccard_sparse_fill_cpp(x, A, transpose = transpose) - - # Convert min_sums to similarity: sim = ms / (s_i + s_j - ms) - col_idx <- rep(seq_along(diff(A@p)), diff(A@p)) - row_idx <- A@i + 1L - denom <- totals[row_idx] + totals[col_idx] - A@x - nonzero <- denom != 0 - A@x[nonzero] <- A@x[nonzero] / denom[nonzero] - A@x[!nonzero] <- 0 - Matrix::diag(A) <- 1 - - if (sparse) A else as.matrix(A) + sim <- natcpp::c_weighted_jaccard_sparse(x, transpose = transpose) + dimnames(sim) <- list(nms, nms) + if (sparse) sim else as.matrix(sim) } jaccard_sim_weighted_cpp_dense <- function(x, sparse = FALSE, transpose = FALSE) { @@ -205,7 +180,7 @@ jaccard_sim_weighted_cpp_dense <- function(x, sparse = FALSE, transpose = FALSE) if (sparse) return(Matrix::Matrix(sim, sparse = TRUE)) return(sim) } - sim <- weighted_jaccard_dense_cpp(x, transpose = transpose) + sim <- natcpp::c_weighted_jaccard_dense(x, transpose = transpose) dimnames(sim) <- list(nms, nms) if (sparse) Matrix::Matrix(sim, sparse = TRUE) else sim } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp deleted file mode 100644 index fb0ccbd..0000000 --- a/src/RcppExports.cpp +++ /dev/null @@ -1,48 +0,0 @@ -// Generated by using Rcpp::compileAttributes() -> do not edit by hand -// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 - -#include - -using namespace Rcpp; - -#ifdef RCPP_USE_GLOBAL_ROSTREAM -Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); -Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); -#endif - -// weighted_jaccard_sparse_fill_cpp -NumericVector weighted_jaccard_sparse_fill_cpp(const S4& x, const S4& pattern, bool transpose); -RcppExport SEXP _coconat_weighted_jaccard_sparse_fill_cpp(SEXP xSEXP, SEXP patternSEXP, SEXP transposeSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const S4& >::type x(xSEXP); - Rcpp::traits::input_parameter< const S4& >::type pattern(patternSEXP); - Rcpp::traits::input_parameter< bool >::type transpose(transposeSEXP); - rcpp_result_gen = Rcpp::wrap(weighted_jaccard_sparse_fill_cpp(x, pattern, transpose)); - return rcpp_result_gen; -END_RCPP -} -// weighted_jaccard_dense_cpp -NumericMatrix weighted_jaccard_dense_cpp(const S4& x, bool transpose); -RcppExport SEXP _coconat_weighted_jaccard_dense_cpp(SEXP xSEXP, SEXP transposeSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const S4& >::type x(xSEXP); - Rcpp::traits::input_parameter< bool >::type transpose(transposeSEXP); - rcpp_result_gen = Rcpp::wrap(weighted_jaccard_dense_cpp(x, transpose)); - return rcpp_result_gen; -END_RCPP -} - -static const R_CallMethodDef CallEntries[] = { - {"_coconat_weighted_jaccard_sparse_fill_cpp", (DL_FUNC) &_coconat_weighted_jaccard_sparse_fill_cpp, 3}, - {"_coconat_weighted_jaccard_dense_cpp", (DL_FUNC) &_coconat_weighted_jaccard_dense_cpp, 2}, - {NULL, NULL, 0} -}; - -RcppExport void R_init_coconat(DllInfo *dll) { - R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); - R_useDynamicSymbols(dll, FALSE); -} diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp deleted file mode 100644 index 853621f..0000000 --- a/src/weighted_jaccard.cpp +++ /dev/null @@ -1,272 +0,0 @@ -#include -#include -#include - -using namespace Rcpp; - -// Sparse accumulation: given the sparsity pattern from crossprod(binarise(x)) -// as a dgCMatrix, fill in min_sums by iterating one output column at a time. -// For each output column cb, populate a row->position lookup, then iterate -// over all features where cb is nonzero, accumulating min(x[f,ca], x[f,cb]) -// for each co-occurring item ca. Writes to out[] are contiguous per column. -// [[Rcpp::export]] -NumericVector weighted_jaccard_sparse_fill_cpp( - const S4& x, const S4& pattern, bool transpose = false) { - - // Original matrix slots (dgCMatrix: nr x nc) - IntegerVector dims = x.slot("Dim"); - IntegerVector xi = x.slot("i"); - IntegerVector xp = x.slot("p"); - NumericVector xx = x.slot("x"); - - const int nr = dims[0]; - const int nc = dims[1]; - const int nfeat = transpose ? nc : nr; - - // Pattern matrix slots (dgCMatrix: ncomp x ncomp) - IntegerVector Ai = pattern.slot("i"); - IntegerVector Ap = pattern.slot("p"); - IntegerVector Adims = pattern.slot("Dim"); - const int ncomp = Adims[0]; - const int annz = Ai.size(); - - // Output: min_sums values aligned with pattern's (i, p) structure - NumericVector out(annz, 0.0); - - // Build CSR for feature dimension: for each feature f, the items (columns - // or rows being compared) that are nonzero in that feature, with their values. - std::vector feat_count(nfeat, 0); - if (!transpose) { - for (int col = 0; col < nc; ++col) - for (int idx = xp[col]; idx < xp[col + 1]; ++idx) - feat_count[xi[idx]]++; - } else { - for (int col = 0; col < nc; ++col) - feat_count[col] = xp[col + 1] - xp[col]; - } - - std::vector feat_offsets(nfeat + 1, 0); - for (int f = 0; f < nfeat; ++f) - feat_offsets[f + 1] = feat_offsets[f] + feat_count[f]; - - const int total_nnz = feat_offsets[nfeat]; - std::vector feat_items(total_nnz); // which item (col/row being compared) - std::vector feat_vals(total_nnz); // value of x at that position - std::vector feat_fill(feat_offsets.begin(), feat_offsets.end()); - - if (!transpose) { - for (int col = 0; col < nc; ++col) { - for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { - const int row = xi[idx]; - const int dest = feat_fill[row]++; - feat_items[dest] = col; - feat_vals[dest] = xx[idx]; - } - } - } else { - for (int col = 0; col < nc; ++col) { - for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { - const int dest = feat_fill[col]++; - feat_items[dest] = xi[idx]; - feat_vals[dest] = xx[idx]; - } - } - } - - // Also need item->features CSR: for each compared item, which features - // it participates in and with what value. - std::vector item_count(ncomp, 0); - for (int f = 0; f < nfeat; ++f) { - for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) - item_count[feat_items[idx]]++; - } - - std::vector item_offsets(ncomp + 1, 0); - for (int c = 0; c < ncomp; ++c) - item_offsets[c + 1] = item_offsets[c] + item_count[c]; - - std::vector item_feats(total_nnz); // which feature - std::vector item_vals(total_nnz); // value - std::vector item_fill(item_offsets.begin(), item_offsets.end()); - - for (int f = 0; f < nfeat; ++f) { - for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) { - const int item = feat_items[idx]; - const int dest = item_fill[item]++; - item_feats[dest] = f; - item_vals[dest] = feat_vals[idx]; - } - } - - // row_to_pos[row] = position in out[] for entry (row, current_col). - // Populated once per output column, then cleared. - std::vector row_to_pos(ncomp, -1); - - // Iterate over output columns - for (int cb = 0; cb < ncomp; ++cb) { - // Populate row_to_pos for pattern column cb - for (int idx = Ap[cb]; idx < Ap[cb + 1]; ++idx) - row_to_pos[Ai[idx]] = idx; - - // For each feature f where cb is nonzero - for (int fi = item_offsets[cb]; fi < item_offsets[cb + 1]; ++fi) { - const int f = item_feats[fi]; - const double vb = item_vals[fi]; - - // For each other item ca also in this feature, accumulate min into (ca, cb) - for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) { - const int ca = feat_items[idx]; - const int pos = row_to_pos[ca]; - if (pos >= 0) - out[pos] += std::min(feat_vals[idx], vb); - } - } - - // Clear row_to_pos - for (int idx = Ap[cb]; idx < Ap[cb + 1]; ++idx) - row_to_pos[Ai[idx]] = -1; - } - - return out; -} - -// Dense accumulation with adaptive strategy: -// - Small ncomp (output fits in cache): feature-oriented loop, no extra setup -// - Large ncomp: column-oriented loop with item->features CSR for stride-1 writes -// [[Rcpp::export]] -NumericMatrix weighted_jaccard_dense_cpp(const S4& x, bool transpose = false) { - IntegerVector dims = x.slot("Dim"); - IntegerVector xi = x.slot("i"); - IntegerVector xp = x.slot("p"); - NumericVector xx = x.slot("x"); - - const int nr = dims[0]; - const int nc = dims[1]; - const int ncomp = transpose ? nr : nc; - const int nfeat = transpose ? nc : nr; - - // Totals per compared item - std::vector totals(ncomp, 0.0); - - // Build feature CSR: for each feature f, which items are nonzero - std::vector feat_count(nfeat, 0); - if (!transpose) { - for (int col = 0; col < nc; ++col) { - for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { - totals[col] += xx[idx]; - feat_count[xi[idx]]++; - } - } - } else { - for (int col = 0; col < nc; ++col) { - feat_count[col] = xp[col + 1] - xp[col]; - for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { - totals[xi[idx]] += xx[idx]; - } - } - } - - std::vector feat_offsets(nfeat + 1, 0); - for (int f = 0; f < nfeat; ++f) - feat_offsets[f + 1] = feat_offsets[f] + feat_count[f]; - - const int total_nnz = feat_offsets[nfeat]; - std::vector feat_items(total_nnz); - std::vector feat_vals(total_nnz); - std::vector feat_fill(feat_offsets.begin(), feat_offsets.end()); - - if (!transpose) { - for (int col = 0; col < nc; ++col) { - for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { - const int row = xi[idx]; - const int dest = feat_fill[row]++; - feat_items[dest] = col; - feat_vals[dest] = xx[idx]; - } - } - } else { - for (int col = 0; col < nc; ++col) { - for (int idx = xp[col]; idx < xp[col + 1]; ++idx) { - const int dest = feat_fill[col]++; - feat_items[dest] = xi[idx]; - feat_vals[dest] = xx[idx]; - } - } - } - - NumericMatrix out(ncomp, ncomp); - - // Threshold: ncomp^2 * 8 bytes > ~12MB (typical L3 share per core) - const bool use_colwise = (static_cast(ncomp) * ncomp * 8 > 12 * 1024 * 1024); - - if (use_colwise) { - // Column-oriented: build item->features CSR, iterate output columns - std::vector item_count(ncomp, 0); - for (int f = 0; f < nfeat; ++f) - for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) - item_count[feat_items[idx]]++; - - std::vector item_offsets(ncomp + 1, 0); - for (int c = 0; c < ncomp; ++c) - item_offsets[c + 1] = item_offsets[c] + item_count[c]; - - std::vector item_feats(total_nnz); - std::vector item_vals(total_nnz); - std::vector item_fill(item_offsets.begin(), item_offsets.end()); - - for (int f = 0; f < nfeat; ++f) { - for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) { - const int item = feat_items[idx]; - const int dest = item_fill[item]++; - item_feats[dest] = f; - item_vals[dest] = feat_vals[idx]; - } - } - - for (int cb = 0; cb < ncomp; ++cb) { - double* col = &out[static_cast(cb) * ncomp]; - for (int fi = item_offsets[cb]; fi < item_offsets[cb + 1]; ++fi) { - const int f = item_feats[fi]; - const double vb = item_vals[fi]; - for (int idx = feat_offsets[f]; idx < feat_offsets[f + 1]; ++idx) { - col[feat_items[idx]] += std::min(feat_vals[idx], vb); - } - } - } - } else { - // Feature-oriented: iterate features, scatter into output - for (int f = 0; f < nfeat; ++f) { - const int start = feat_offsets[f]; - const int end = feat_offsets[f + 1]; - const int k = end - start; - - for (int a = 0; a < k; ++a) { - const int ca = feat_items[start + a]; - const double va = feat_vals[start + a]; - out[static_cast(ca) * ncomp + ca] += va; - for (int b = a + 1; b < k; ++b) { - const int cb = feat_items[start + b]; - const double mv = std::min(va, feat_vals[start + b]); - out[static_cast(ca) * ncomp + cb] += mv; - out[static_cast(cb) * ncomp + ca] += mv; - } - } - } - } - - // Convert min_sums to similarity - for (int cb = 0; cb < ncomp; ++cb) { - for (int ca = 0; ca < ncomp; ++ca) { - if (ca == cb) { - out[static_cast(cb) * ncomp + ca] = 1.0; - } else { - const std::size_t pos = static_cast(cb) * ncomp + ca; - const double ms = out[pos]; - const double denom = totals[ca] + totals[cb] - ms; - out[pos] = (denom > 0.0) ? ms / denom : 0.0; - } - } - } - - return out; -} diff --git a/tests/testthat/test-similarity.R b/tests/testthat/test-similarity.R index 8c9ec99..d6c9592 100644 --- a/tests/testthat/test-similarity.R +++ b/tests/testthat/test-similarity.R @@ -149,7 +149,11 @@ test_that("all weighted jaccard methods agree", { ref_t <- jaccard_sim(m_sparse, weighted = TRUE, transpose = TRUE, weighted_method = "dense") - for (method in c("sparse", "dense", "cpp_dense", "cpp_sparse")) { + r_methods <- c("sparse", "dense") + cpp_methods <- if (requireNamespace("natcpp", quietly = TRUE)) + c("cpp_dense", "cpp_sparse") else character(0) + + for (method in c(r_methods, cpp_methods)) { suppressWarnings({ # Default (dense return, transpose=FALSE) expect_equal( From 1adb342081d7a4066fcb21ebcb390387d018ba13 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Fri, 3 Jul 2026 08:27:36 +0100 Subject: [PATCH 08/10] Add triangle=/distance= options to similarity functions cosine_sim(), jaccard_sim(), tanimoto_sim(), and connectivity_similarity() gain: - triangle=TRUE: return a dist object (lower triangle only) instead of a square matrix, matching the c_weighted_jaccard_* API in natcpp 0.3.0 - distance=TRUE: return 1 - similarity in place of similarity New helpers sim_to_output() and dsCMatrix_to_dist() centralise the conversion between full matrix / sparse triangular / dist layouts, avoiding a densification step when converting a symmetric sparse similarity to a dist. A new exported sim2dist() exposes this sparse-to-dist conversion directly. The pure-R weighted Jaccard fallback previously warned on every call; now warn_natcpp() uses memoise to rate-limit the message to once per hour per distinct text. Adds memoise to Imports. --- DESCRIPTION | 1 + R/cosine.R | 15 +- R/similarity.R | 260 +++++++++++++++++++++++++++++---- man/connectivity_similarity.Rd | 13 +- man/cosine_sim.Rd | 18 ++- man/jaccard_sim.Rd | 13 +- man/sim2dist.Rd | 21 +++ man/tanimoto_sim.Rd | 17 ++- 8 files changed, 312 insertions(+), 46 deletions(-) create mode 100644 man/sim2dist.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 3672d37..4d06285 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,6 +19,7 @@ Imports: checkmate, glue, Matrix (>= 1.5-3), + memoise, methods, grDevices, stats diff --git a/R/cosine.R b/R/cosine.R index 3ec1e56..90ef272 100644 --- a/R/cosine.R +++ b/R/cosine.R @@ -7,9 +7,13 @@ #' @param sparse Whether to return a sparse or dense matrix (default dense) #' @param transpose When \code{F} (the default) calculates the cosine distance #' between columns. When \code{T} calculates the distance between rows. +#' @param triangle If \code{TRUE}, return a \code{\link{dist}} object (lower +#' triangle only, half memory). Default \code{FALSE}. +#' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) +#' instead of similarity. Default \code{FALSE}. #' -#' @return A square matrix, dense unless \code{sparse=TRUE} and \code{x} is -#' sparse. +#' @return A square matrix, or a \code{\link{dist}} object when +#' \code{triangle = TRUE}. #' @export #' #' @examples @@ -24,13 +28,14 @@ #' kckc.cos=cosine_sim(fam_pnkc2) #' pnpn.cos=cosine_sim(fam_pnkc2, transpose=T) #' } -cosine_sim <- function(x, sparse=FALSE, transpose=FALSE) { +cosine_sim <- function(x, sparse=FALSE, transpose=FALSE, + triangle=FALSE, distance=FALSE) { cx=class(x) if(!is.matrix(x) && !isTRUE(attr(cx, "package") == "Matrix")) stop("I don't recognise that as a matrix!") cpx <- if(transpose) Matrix::tcrossprod(x) else Matrix::crossprod(x) - cosx=Matrix::cov2cor(cpx) - if(sparse) cosx else as.matrix(cosx) + sim=Matrix::cov2cor(cpx) + sim_to_output(sim, sparse=sparse, triangle=triangle, distance=distance) } diff --git a/R/similarity.R b/R/similarity.R index 17d091e..24da0e1 100644 --- a/R/similarity.R +++ b/R/similarity.R @@ -14,8 +14,13 @@ #' @param sparse Whether to return a sparse matrix (default \code{FALSE}) #' @param transpose When \code{FALSE} (the default) calculates similarity #' between columns. When \code{TRUE} calculates similarity between rows. +#' @param triangle If \code{TRUE}, return a \code{\link{dist}} object (lower +#' triangle only, half memory). Default \code{FALSE}. +#' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) +#' instead of similarity. Default \code{FALSE}. #' -#' @return A square similarity matrix with values in \code{[0,1]}. +#' @return A square similarity matrix, or a \code{\link{dist}} object when +#' \code{triangle = TRUE}. #' @export #' @seealso \code{\link{cosine_sim}}, \code{\link{jaccard_sim}}, #' \code{\link{tanimoto_sim}} @@ -27,17 +32,46 @@ #' connectivity_similarity(am, metric="weighted_jaccard") #' connectivity_similarity(am, metric="tanimoto") connectivity_similarity <- function(x, metric = c("cosine", "jaccard", "weighted_jaccard", "tanimoto"), - sparse = FALSE, transpose = FALSE) { + sparse = FALSE, transpose = FALSE, + triangle = FALSE, distance = FALSE) { metric <- match.arg(metric) switch(metric, - cosine = cosine_sim(x, sparse = sparse, transpose = transpose), - jaccard = jaccard_sim(x, weighted = FALSE, sparse = sparse, transpose = transpose), - weighted_jaccard = jaccard_sim(x, weighted = TRUE, sparse = sparse, transpose = transpose), - tanimoto = tanimoto_sim(x, sparse = sparse, transpose = transpose) + cosine = cosine_sim(x, sparse = sparse, transpose = transpose, + triangle = triangle, distance = distance), + jaccard = jaccard_sim(x, weighted = FALSE, sparse = sparse, transpose = transpose, + triangle = triangle, distance = distance), + weighted_jaccard = jaccard_sim(x, weighted = TRUE, sparse = sparse, transpose = transpose, + triangle = triangle, distance = distance), + tanimoto = tanimoto_sim(x, sparse = sparse, transpose = transpose, + triangle = triangle, distance = distance) ) } +warn_hourly_check <- memoise::memoise(function(msg) TRUE, + ~ memoise::timeout(3600)) + +warn_hourly <- function(msg) { + cached <- memoise::has_cache(warn_hourly_check)(msg) + warn_hourly_check(msg) # populate cache before warning + if (!cached) + warning(msg, call. = FALSE, immediate. = TRUE) +} + +warn_natcpp <- function() { + has <- requireNamespace("natcpp", quietly = TRUE) + if (has) { + msg <- paste0("natcpp (>= 0.3.0) is required for fast weighted Jaccard. ", + "You have ", utils::packageVersion("natcpp"), ". ", + "Please update with: ", + "install.packages('natcpp', repos='https://natverse.r-universe.dev')") + } else { + msg <- paste0("Install the natcpp package for much faster weighted Jaccard: ", + "install.packages('natcpp', repos='https://natverse.r-universe.dev')") + } + warn_hourly(msg) +} + #' Jaccard similarity for sparse or dense matrices #' #' @description Computes pairwise Jaccard similarity between columns (or rows) @@ -65,8 +99,13 @@ connectivity_similarity <- function(x, metric = c("cosine", "jaccard", "weighted #' @param sparse Whether to return a sparse matrix (default \code{FALSE}) #' @param transpose When \code{FALSE} (the default) calculates similarity #' between columns. When \code{TRUE} calculates similarity between rows. +#' @param triangle If \code{TRUE}, return a \code{\link{dist}} object (lower +#' triangle only, half memory). Default \code{FALSE}. +#' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) +#' instead of similarity. Default \code{FALSE}. #' -#' @return A square similarity matrix with values in \code{[0,1]}. +#' @return A square similarity matrix, or a \code{\link{dist}} object when +#' \code{triangle = TRUE}. #' @importFrom methods as #' @export #' @seealso \code{\link{cosine_sim}}, \code{\link{connectivity_similarity}} @@ -79,7 +118,8 @@ connectivity_similarity <- function(x, metric = c("cosine", "jaccard", "weighted #' jaccard_sim(am, weighted=TRUE) jaccard_sim <- function(x, weighted = FALSE, sparse = FALSE, transpose = FALSE, weighted_method = c("auto", "cpp_dense", "cpp_sparse", - "dense", "sparse")) { + "dense", "sparse"), + triangle = FALSE, distance = FALSE) { cx <- class(x) if (!is.matrix(x) && !isTRUE(attr(cx, "package") == "Matrix")) stop("I don't recognise that as a matrix!") @@ -102,7 +142,9 @@ jaccard_sim <- function(x, weighted = FALSE, sparse = FALSE, transpose = FALSE, isect <- A@x A@x <- isect / (sizes[row_idx] + sizes[col_idx] - isect) Matrix::diag(A) <- 1 - sim <- A + dimnames(A) <- list(nms, nms) + return(sim_to_output(A, sparse = sparse, triangle = triangle, + distance = distance)) } else { if (weighted_method == "auto") { has_cpp <- requireNamespace("natcpp", quietly = TRUE) && @@ -111,18 +153,24 @@ jaccard_sim <- function(x, weighted = FALSE, sparse = FALSE, transpose = FALSE, if (has_cpp) { weighted_method <- if (ncomp > 10000L) "cpp_sparse" else "cpp_dense" } else { + warn_natcpp() weighted_method <- "dense" } } sim <- switch(weighted_method, - cpp_dense = jaccard_sim_weighted_cpp_dense(x, sparse = sparse, transpose = transpose), - cpp_sparse = jaccard_sim_weighted_cpp_sparse(x, sparse = sparse, transpose = transpose), - dense = jaccard_sim_weighted_dense_r(x, sparse = sparse, transpose = transpose), - sparse = jaccard_sim_weighted_sparse_r(x, sparse = TRUE, transpose = transpose) + cpp_dense = jaccard_sim_weighted_cpp_dense(x, sparse = sparse, transpose = transpose, + triangle = triangle, distance = distance), + cpp_sparse = jaccard_sim_weighted_cpp_sparse(x, sparse = sparse, transpose = transpose, + triangle = triangle, distance = distance), + dense = jaccard_sim_weighted_dense_r(x, sparse = sparse, transpose = transpose, + triangle = triangle, distance = distance), + sparse = jaccard_sim_weighted_sparse_r(x, sparse = TRUE, transpose = transpose, + triangle = triangle, distance = distance) ) } - if (sparse) sim else as.matrix(sim) + if (inherits(sim, "dist")) sim + else if (sparse) sim else as.matrix(sim) } jaccard_weighted_feature_view <- function(x, transpose = FALSE) { @@ -163,36 +211,58 @@ jaccard_weighted_feature_view <- function(x, transpose = FALSE) { ) } -jaccard_sim_weighted_cpp_sparse <- function(x, sparse = TRUE, transpose = FALSE) { +jaccard_sim_weighted_cpp_sparse <- function(x, sparse = TRUE, transpose = FALSE, + triangle = FALSE, distance = FALSE) { n <- if (transpose) nrow(x) else ncol(x) nms <- if (transpose) rownames(x) else colnames(x) - sim <- natcpp::c_weighted_jaccard_sparse(x, transpose = transpose) + sim <- natcpp::c_weighted_jaccard_sparse(x, transpose = transpose, + triangle = triangle, distance = distance) dimnames(sim) <- list(nms, nms) if (sparse) sim else as.matrix(sim) } -jaccard_sim_weighted_cpp_dense <- function(x, sparse = FALSE, transpose = FALSE) { +jaccard_sim_weighted_cpp_dense <- function(x, sparse = FALSE, transpose = FALSE, + triangle = FALSE, distance = FALSE) { n <- if (transpose) nrow(x) else ncol(x) nms <- if (transpose) rownames(x) else colnames(x) if (length(x@x) == 0L) { - sim <- diag(n) + if (triangle) { + d <- rep(if (distance) 1 else 0, n * (n - 1L) / 2L) + return(structure(d, Size = n, Labels = nms, Diag = FALSE, Upper = FALSE, + class = "dist")) + } + sim <- if (distance) matrix(1, n, n) - diag(n) else diag(n) dimnames(sim) <- list(nms, nms) if (sparse) return(Matrix::Matrix(sim, sparse = TRUE)) return(sim) } - sim <- natcpp::c_weighted_jaccard_dense(x, transpose = transpose) - dimnames(sim) <- list(nms, nms) - if (sparse) Matrix::Matrix(sim, sparse = TRUE) else sim + res <- natcpp::c_weighted_jaccard_dense(x, transpose = transpose, + triangle = triangle, distance = distance) + if (triangle) { + attr(res, "Labels") <- nms + return(res) + } + dimnames(res) <- list(nms, nms) + if (sparse) Matrix::Matrix(res, sparse = TRUE) else res } -jaccard_sim_weighted_sparse_r <- function(x, sparse = TRUE, transpose = FALSE) { - warning("Using slow pure R sparse weighted Jaccard. ", - "Install the natcpp package for much faster compiled code.", - call. = FALSE) +jaccard_sim_weighted_sparse_r <- function(x, sparse = TRUE, transpose = FALSE, + triangle = FALSE, distance = FALSE) { + warn_natcpp() + if (distance) + warning("distance=TRUE with sparse output produces a mostly-dense matrix; ", + "consider using the dense backend with triangle=TRUE instead.", + call. = FALSE) n <- if (transpose) nrow(x) else ncol(x) nms <- if (transpose) rownames(x) else colnames(x) if (length(x@x) == 0L) { - sim <- Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, + if (triangle) { + d <- rep(if (distance) 1 else 0, n * (n - 1L) / 2L) + return(structure(d, Size = n, Labels = nms, Diag = FALSE, Upper = FALSE, + class = "dist")) + } + sim_val <- if (distance) 0 else 1 + sim <- Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = sim_val, dims = c(n, n), dimnames = list(nms, nms)) return(if (sparse) sim else as.matrix(sim)) } @@ -259,14 +329,31 @@ jaccard_sim_weighted_sparse_r <- function(x, sparse = TRUE, transpose = FALSE) { min_sums@x[!nz] <- 0 Matrix::diag(min_sums) <- 1 + if (triangle) { + dm <- as.matrix(min_sums) + d <- dm[lower.tri(dm)] + if (distance) d <- 1 - d + return(structure(d, Size = n, Labels = nms, Diag = FALSE, Upper = FALSE, + class = "dist")) + } + + if (distance) + min_sums@x <- 1 - min_sums@x + if (sparse) min_sums else as.matrix(min_sums) } -jaccard_sim_weighted_dense_r <- function(x, sparse = FALSE, transpose = FALSE) { +jaccard_sim_weighted_dense_r <- function(x, sparse = FALSE, transpose = FALSE, + triangle = FALSE, distance = FALSE) { n <- if (transpose) nrow(x) else ncol(x) nms <- if (transpose) rownames(x) else colnames(x) if (length(x@x) == 0L) { - sim <- diag(1, n) + if (triangle) { + d <- rep(if (distance) 1 else 0, n * (n - 1L) / 2L) + return(structure(d, Size = n, Labels = nms, Diag = FALSE, Upper = FALSE, + class = "dist")) + } + sim <- if (distance) matrix(1, n, n) - diag(n) else diag(1, n) dimnames(sim) <- list(nms, nms) return(if (sparse) Matrix::Matrix(sim, sparse = TRUE) else sim) } @@ -294,8 +381,17 @@ jaccard_sim_weighted_dense_r <- function(x, sparse = FALSE, transpose = FALSE) { out[!nz, j] <- 0 out[j, j] <- 1 } - dimnames(out) <- list(nms, nms) + if (triangle) { + d <- out[lower.tri(out)] + if (distance) d <- 1 - d + return(structure(d, Size = ncomp, Labels = nms, Diag = FALSE, Upper = FALSE, + class = "dist")) + } + + if (distance) out <- 1 - out + + dimnames(out) <- list(nms, nms) if (sparse) Matrix::Matrix(out, sparse = TRUE) else out } @@ -311,8 +407,13 @@ jaccard_sim_weighted_dense_r <- function(x, sparse = FALSE, transpose = FALSE) { #' @param sparse Whether to return a sparse matrix (default \code{FALSE}) #' @param transpose When \code{FALSE} (the default) calculates similarity #' between columns. When \code{TRUE} calculates similarity between rows. +#' @param triangle If \code{TRUE}, return a \code{\link{dist}} object (lower +#' triangle only, half memory). Default \code{FALSE}. +#' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) +#' instead of similarity. Default \code{FALSE}. #' -#' @return A square similarity matrix with values in \code{[0,1]}. +#' @return A square similarity matrix, or a \code{\link{dist}} object when +#' \code{triangle = TRUE}. #' @importFrom methods as #' @export #' @seealso \code{\link{jaccard_sim}}, \code{\link{cosine_sim}}, @@ -321,13 +422,15 @@ jaccard_sim_weighted_dense_r <- function(x, sparse = FALSE, transpose = FALSE) { #' da2ds15=readRDS(system.file('sampledata/da2ds15.rds', package = 'coconat')) #' am=partner_summary2adjacency_matrix(da2ds15, inputcol = 'partner', outputcol = 'bodyid') #' tanimoto_sim(am) -tanimoto_sim <- function(x, sparse = FALSE, transpose = FALSE) { +tanimoto_sim <- function(x, sparse = FALSE, transpose = FALSE, + triangle = FALSE, distance = FALSE) { cx <- class(x) if (!is.matrix(x) && !isTRUE(attr(cx, "package") == "Matrix")) stop("I don't recognise that as a matrix!") if (!inherits(x, "dgCMatrix")) x <- as(x, "dgCMatrix") crossfun <- if (transpose) Matrix::tcrossprod else Matrix::crossprod + nms <- if (transpose) rownames(x) else colnames(x) # dot(a,b) for all pairs via crossprod (returns symmetric dsCMatrix) A <- crossfun(x) @@ -344,7 +447,100 @@ tanimoto_sim <- function(x, sparse = FALSE, transpose = FALSE) { A@x[!nonzero] <- 0 } Matrix::diag(A) <- 1 - sim <- A + dimnames(A) <- list(nms, nms) + + sim_to_output(A, sparse = sparse, triangle = triangle, distance = distance) +} + +#' @noRd +sim_to_output <- function(sim, sparse = FALSE, triangle = FALSE, distance = FALSE) { + if (triangle) { + if (inherits(sim, "dsCMatrix")) { + d <- dsCMatrix_to_dist(sim, distance = distance) + } else { + dm <- as.matrix(sim) + d <- dm[lower.tri(dm)] + if (distance) d <- 1 - d + nms <- rownames(dm) + d <- structure(d, Size = nrow(dm), Labels = nms, Diag = FALSE, + Upper = FALSE, class = "dist") + } + return(d) + } + if (distance) { + sim <- if (is.matrix(sim)) 1 - sim else 1 - as.matrix(sim) + return(sim) + } if (sparse) sim else as.matrix(sim) } + +#' Extract lower triangle from a dsCMatrix as a dist object without densifying +#' @noRd +dsCMatrix_to_dist <- function(x, distance = FALSE) { + n <- nrow(x) + fill_val <- if (distance) 1 else 0 + d <- rep(fill_val, n * (n - 1L) / 2L) + + col0 <- rep(seq_along(diff(x@p)), diff(x@p)) - 1L + row0 <- x@i + vals <- x@x + + # Keep only off-diagonal entries (upper triangle: row0 < col0) + offdiag <- row0 != col0 + row0 <- row0[offdiag] + col0 <- col0[offdiag] + vals <- vals[offdiag] + + if (length(vals)) { + # dsCMatrix stores upper triangle (row < col). + # dist layout is lower triangle column-major: for pair (i,j), i > j: + # pos = (j-1)*n - j*(j-1)/2 + (i-j) [1-based] + # Here row0 < col0 (0-based), so i=col0+1, j=row0+1: + r <- row0 + 1L + cc <- col0 + 1L + pos <- (r - 1L) * n - r * (r - 1L) / 2L + (cc - r) + d[pos] <- if (distance) 1 - vals else vals + } + + structure(d, Size = n, Labels = rownames(x), + Diag = FALSE, Upper = FALSE, class = "dist") +} + +#' Convert a symmetric sparse similarity matrix to a dist object +#' +#' Efficiently converts a symmetric sparse similarity matrix (\code{dsCMatrix}) +#' to a \code{\link{dist}} object without materialising the full dense matrix. +#' Similarity values are converted to distances as \code{1 - similarity}. +#' +#' @param x A symmetric sparse matrix of class \code{dsCMatrix} with values in +#' \code{[0,1]}, as returned by e.g. \code{\link{jaccard_sim}} with +#' \code{sparse=TRUE}. +#' @return A \code{\link{dist}} object +#' @export +sim2dist <- function(x) { + if (!inherits(x, "dsCMatrix")) + stop("x must be a dsCMatrix (symmetric sparse column-compressed matrix)") + n <- nrow(x) + # Pre-fill with distance=1 (similarity=0 for absent pairs) + d <- rep(1, n * (n - 1L) / 2L) + + # Extract upper triangle entries (excluding diagonal) + col0 <- rep(seq_along(diff(x@p)), diff(x@p)) - 1L + row0 <- x@i + vals <- x@x + + offdiag <- row0 != col0 + row0 <- row0[offdiag] + col0 <- col0[offdiag] + vals <- vals[offdiag] + + # Map upper triangle (r,c) r Date: Fri, 3 Jul 2026 08:27:38 +0100 Subject: [PATCH 09/10] Ignore .claude/ in package build --- .Rbuildignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.Rbuildignore b/.Rbuildignore index afdc1d0..58b14db 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -8,3 +8,4 @@ ^codecov\.yml$ ^README\.Rmd$ ^samples$ +^\.claude$ From 4b1289a451d8417c6fac704cf8809805f21700f6 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Fri, 3 Jul 2026 09:22:40 +0100 Subject: [PATCH 10/10] Move natcpp to Enhances and pin Remotes to feature branch MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Suggests(natcpp) triggered pak dependency resolution failures on CI because natcpp 0.3.0 is not yet on CRAN or on natverse/natcpp master (currently 0.2). Following the nat convention, natcpp is more accurately an Enhances relationship — coconat works fine without it, but gets much faster weighted Jaccard when present. Remotes points at the feature/weighted-jaccard branch for now; drop the @feature/... suffix once that branch is merged. --- DESCRIPTION | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 4d06285..f8bd41d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -25,7 +25,6 @@ Imports: stats Suggests: covr, - natcpp (>= 0.3.0), spelling, testthat (>= 3.0.0), ComplexHeatmap, @@ -36,6 +35,8 @@ Suggests: rmarkdown, dplyr, dendroextras +Enhances: natcpp (>= 0.3.0) +Remotes: natverse/natcpp@feature/weighted-jaccard Language: en-GB Config/testthat/edition: 3 URL: https://github.com/natverse/coconat,