From 5f16e4c1f90a5cb6de74f596dca73ed8ac531f18 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Fri, 3 Apr 2026 23:32:00 +0100 Subject: [PATCH 1/8] weighted jaccard functions for coconat --- DESCRIPTION | 4 +- NAMESPACE | 2 + R/RcppExports.R | 35 +++ man/weighted_jaccard_dense_cpp.Rd | 24 ++ man/weighted_jaccard_sparse_fill_cpp.Rd | 26 +++ src/RcppExports.cpp | 27 +++ src/weighted_jaccard.cpp | 291 ++++++++++++++++++++++++ tests/testthat/test-weighted-jaccard.R | 82 +++++++ 8 files changed, 489 insertions(+), 2 deletions(-) create mode 100644 man/weighted_jaccard_dense_cpp.Rd create mode 100644 man/weighted_jaccard_sparse_fill_cpp.Rd create mode 100644 src/weighted_jaccard.cpp create mode 100644 tests/testthat/test-weighted-jaccard.R diff --git a/DESCRIPTION b/DESCRIPTION index 5dff4ca..71d25a7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: natcpp Title: Fast C++ Primitives for the 'NeuroAnatomy Toolbox' -Version: 0.2 +Version: 0.3.0 Authors@R: person(given = "Gregory", family = "Jefferis", @@ -28,4 +28,4 @@ LinkingTo: Config/testthat/edition: 3 Encoding: UTF-8 Language: en-GB -RoxygenNote: 7.3.2 +RoxygenNote: 7.3.3 diff --git a/NAMESPACE b/NAMESPACE index 504bfff..f614ccb 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,5 +10,7 @@ export(c_sub2ind) export(c_topntail) export(c_topntail_list) export(c_total_cable) +export(weighted_jaccard_dense_cpp) +export(weighted_jaccard_sparse_fill_cpp) import(Rcpp) useDynLib(natcpp, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index c90d42d..90fcc19 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -135,3 +135,38 @@ c_EdgeListFromSegList <- function(L) { .Call(`_natcpp_c_EdgeListFromSegList`, L) } +#' Fill min_sums into a sparse pattern for weighted Jaccard similarity +#' +#' Given a dgCMatrix \code{x} and a sparsity pattern \code{pattern} (from +#' \code{crossprod(binarise(x))}), compute the sum of element-wise minima +#' for each pair of columns (or rows) that co-occur in at least one feature. +#' Intended for internal use by \code{coconat::jaccard_sim}. +#' +#' @param x A dgCMatrix (sparse column-compressed matrix) +#' @param pattern A dgCMatrix giving the output sparsity pattern +#' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare +#' rows +#' @return A numeric vector of min_sums aligned with the pattern's slot +#' structure +#' @export +weighted_jaccard_sparse_fill_cpp <- function(x, pattern, transpose = FALSE) { + .Call(`_natcpp_weighted_jaccard_sparse_fill_cpp`, x, pattern, transpose) +} + +#' Dense weighted Jaccard similarity via C++ +#' +#' Compute the full weighted Jaccard similarity matrix for a dgCMatrix using +#' an adaptive dense accumulation strategy. For small output matrices, uses a +#' feature-oriented loop; for large outputs, switches to a column-oriented +#' loop for better cache performance. Intended for internal use by +#' \code{coconat::jaccard_sim}. +#' +#' @param x A dgCMatrix (sparse column-compressed matrix) +#' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare +#' rows +#' @return A dense numeric similarity matrix +#' @export +weighted_jaccard_dense_cpp <- function(x, transpose = FALSE) { + .Call(`_natcpp_weighted_jaccard_dense_cpp`, x, transpose) +} + diff --git a/man/weighted_jaccard_dense_cpp.Rd b/man/weighted_jaccard_dense_cpp.Rd new file mode 100644 index 0000000..5269d34 --- /dev/null +++ b/man/weighted_jaccard_dense_cpp.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{weighted_jaccard_dense_cpp} +\alias{weighted_jaccard_dense_cpp} +\title{Dense weighted Jaccard similarity via C++} +\usage{ +weighted_jaccard_dense_cpp(x, transpose = FALSE) +} +\arguments{ +\item{x}{A dgCMatrix (sparse column-compressed matrix)} + +\item{transpose}{If \code{FALSE}, compare columns; if \code{TRUE}, compare +rows} +} +\value{ +A dense numeric similarity matrix +} +\description{ +Compute the full weighted Jaccard similarity matrix for a dgCMatrix using +an adaptive dense accumulation strategy. For small output matrices, uses a +feature-oriented loop; for large outputs, switches to a column-oriented +loop for better cache performance. Intended for internal use by +\code{coconat::jaccard_sim}. +} diff --git a/man/weighted_jaccard_sparse_fill_cpp.Rd b/man/weighted_jaccard_sparse_fill_cpp.Rd new file mode 100644 index 0000000..00b597e --- /dev/null +++ b/man/weighted_jaccard_sparse_fill_cpp.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{weighted_jaccard_sparse_fill_cpp} +\alias{weighted_jaccard_sparse_fill_cpp} +\title{Fill min_sums into a sparse pattern for weighted Jaccard similarity} +\usage{ +weighted_jaccard_sparse_fill_cpp(x, pattern, transpose = FALSE) +} +\arguments{ +\item{x}{A dgCMatrix (sparse column-compressed matrix)} + +\item{pattern}{A dgCMatrix giving the output sparsity pattern} + +\item{transpose}{If \code{FALSE}, compare columns; if \code{TRUE}, compare +rows} +} +\value{ +A numeric vector of min_sums aligned with the pattern's slot + structure +} +\description{ +Given a dgCMatrix \code{x} and a sparsity pattern \code{pattern} (from +\code{crossprod(binarise(x))}), compute the sum of element-wise minima +for each pair of columns (or rows) that co-occur in at least one feature. +Intended for internal use by \code{coconat::jaccard_sim}. +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 3415f85..6db9ee8 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -135,6 +135,31 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// weighted_jaccard_sparse_fill_cpp +NumericVector weighted_jaccard_sparse_fill_cpp(const S4& x, const S4& pattern, bool transpose); +RcppExport SEXP _natcpp_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 _natcpp_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[] = { {"_natcpp_c_seglengths", (DL_FUNC) &_natcpp_c_seglengths, 4}, @@ -147,6 +172,8 @@ static const R_CallMethodDef CallEntries[] = { {"_natcpp_c_topntail", (DL_FUNC) &_natcpp_c_topntail, 1}, {"_natcpp_c_topntail_list", (DL_FUNC) &_natcpp_c_topntail_list, 1}, {"_natcpp_c_EdgeListFromSegList", (DL_FUNC) &_natcpp_c_EdgeListFromSegList, 1}, + {"_natcpp_weighted_jaccard_sparse_fill_cpp", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill_cpp, 3}, + {"_natcpp_weighted_jaccard_dense_cpp", (DL_FUNC) &_natcpp_weighted_jaccard_dense_cpp, 2}, {NULL, NULL, 0} }; diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp new file mode 100644 index 0000000..58f4fae --- /dev/null +++ b/src/weighted_jaccard.cpp @@ -0,0 +1,291 @@ +#include +#include +#include + +using namespace Rcpp; + +//' Fill min_sums into a sparse pattern for weighted Jaccard similarity +//' +//' Given a dgCMatrix \code{x} and a sparsity pattern \code{pattern} (from +//' \code{crossprod(binarise(x))}), compute the sum of element-wise minima +//' for each pair of columns (or rows) that co-occur in at least one feature. +//' Intended for internal use by \code{coconat::jaccard_sim}. +//' +//' @param x A dgCMatrix (sparse column-compressed matrix) +//' @param pattern A dgCMatrix giving the output sparsity pattern +//' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare +//' rows +//' @return A numeric vector of min_sums aligned with the pattern's slot +//' structure +//' @export +// [[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 weighted Jaccard similarity via C++ +//' +//' Compute the full weighted Jaccard similarity matrix for a dgCMatrix using +//' an adaptive dense accumulation strategy. For small output matrices, uses a +//' feature-oriented loop; for large outputs, switches to a column-oriented +//' loop for better cache performance. Intended for internal use by +//' \code{coconat::jaccard_sim}. +//' +//' @param x A dgCMatrix (sparse column-compressed matrix) +//' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare +//' rows +//' @return A dense numeric similarity matrix +//' @export +// [[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-weighted-jaccard.R b/tests/testthat/test-weighted-jaccard.R new file mode 100644 index 0000000..207ab32 --- /dev/null +++ b/tests/testthat/test-weighted-jaccard.R @@ -0,0 +1,82 @@ +test_that("weighted_jaccard_dense_cpp computes correct similarity", { + skip_if_not_installed("Matrix") + m <- Matrix::sparseMatrix( + i = c(1L, 2L, 1L, 2L, 3L, 3L), + j = c(1L, 1L, 2L, 2L, 2L, 3L), + x = c(4, 2, 1, 3, 3, 1), + dims = c(3L, 3L) + ) + + # Naive reference: sim[a,b] = sum(pmin(col_a, col_b)) / sum(pmax(col_a, col_b)) + dm <- as.matrix(m) + n <- ncol(dm) + ref <- matrix(0, n, n) + for (a in seq_len(n)) { + for (b in seq_len(n)) { + mins <- sum(pmin(dm[, a], dm[, b])) + maxs <- sum(pmax(dm[, a], dm[, b])) + ref[a, b] <- if (maxs > 0) mins / maxs else 0 + } + } + diag(ref) <- 1 + + sim <- weighted_jaccard_dense_cpp(m, transpose = FALSE) + expect_equal(sim, ref, tolerance = 1e-12) +}) + +test_that("weighted_jaccard_dense_cpp transpose works", { + skip_if_not_installed("Matrix") + m <- Matrix::sparseMatrix( + i = c(1L, 2L, 1L, 2L, 3L, 3L), + j = c(1L, 1L, 2L, 2L, 2L, 3L), + x = c(4, 2, 1, 3, 3, 1), + dims = c(3L, 3L) + ) + + dm <- as.matrix(m) + nr <- nrow(dm) + ref_t <- matrix(0, nr, nr) + for (a in seq_len(nr)) { + for (b in seq_len(nr)) { + mins <- sum(pmin(dm[a, ], dm[b, ])) + maxs <- sum(pmax(dm[a, ], dm[b, ])) + ref_t[a, b] <- if (maxs > 0) mins / maxs else 0 + } + } + diag(ref_t) <- 1 + + sim_t <- weighted_jaccard_dense_cpp(m, transpose = TRUE) + expect_equal(sim_t, ref_t, tolerance = 1e-12) +}) + +test_that("weighted_jaccard_sparse_fill_cpp matches dense", { + skip_if_not_installed("Matrix") + m <- Matrix::sparseMatrix( + i = c(1L, 2L, 1L, 2L, 3L, 3L), + j = c(1L, 1L, 2L, 2L, 2L, 3L), + x = c(4, 2, 1, 3, 3, 1), + dims = c(3L, 3L) + ) + + # Dense reference + sim_dense <- weighted_jaccard_dense_cpp(m, transpose = FALSE) + + # Sparse fill: build pattern from binarised crossprod + b <- m + b@x <- rep(1, length(b@x)) + A <- methods::as(Matrix::crossprod(b), "generalMatrix") + totals <- Matrix::colSums(m) + + A@x <- weighted_jaccard_sparse_fill_cpp(m, A, transpose = FALSE) + + # Convert min_sums to similarity + 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 + + expect_equal(as.matrix(A), sim_dense, tolerance = 1e-12) +}) From cc7e20f21c75dcb1a9dd2025d5a7b867c6851d21 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Sat, 4 Apr 2026 09:25:36 +0100 Subject: [PATCH 2/8] suggest Matrix package --- DESCRIPTION | 1 + 1 file changed, 1 insertion(+) diff --git a/DESCRIPTION b/DESCRIPTION index 71d25a7..3d39b65 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -21,6 +21,7 @@ BugReports: https://github.com/natverse/natcpp/issues Imports: Rcpp (>= 1.0.6) Suggests: + Matrix, spelling, testthat (>= 3.0.0) LinkingTo: From 1895e525871b987a35f89a4290e438ba88dcd7db Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Sat, 4 Apr 2026 16:03:30 +0100 Subject: [PATCH 3/8] make c_weighted_jaccard_sparse a full implementation * this mean that it can be tested for equivalence * and rename both functions --- DESCRIPTION | 6 +-- NAMESPACE | 4 +- R/RcppExports.R | 25 +++------ R/weighted_jaccard.R | 51 +++++++++++++++++++ ...nse_cpp.Rd => c_weighted_jaccard_dense.Rd} | 8 +-- man/c_weighted_jaccard_sparse.Rd | 34 +++++++++++++ man/weighted_jaccard_sparse_fill_cpp.Rd | 26 ---------- src/RcppExports.cpp | 20 ++++---- src/weighted_jaccard.cpp | 21 ++------ tests/testthat/test-weighted-jaccard.R | 39 ++++---------- 10 files changed, 125 insertions(+), 109 deletions(-) create mode 100644 R/weighted_jaccard.R rename man/{weighted_jaccard_dense_cpp.Rd => c_weighted_jaccard_dense.Rd} (77%) create mode 100644 man/c_weighted_jaccard_sparse.Rd delete mode 100644 man/weighted_jaccard_sparse_fill_cpp.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 3d39b65..10934b6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -18,10 +18,10 @@ Description: Fast functions implemented in C++ via 'Rcpp' to support the License: GPL (>= 3) URL: https://github.com/natverse/natcpp, https://natverse.org/natcpp/ BugReports: https://github.com/natverse/natcpp/issues -Imports: - Rcpp (>= 1.0.6) -Suggests: +Imports: Matrix, + Rcpp (>= 1.0.6) +Suggests: spelling, testthat (>= 3.0.0) LinkingTo: diff --git a/NAMESPACE b/NAMESPACE index f614ccb..355a815 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,7 +10,7 @@ export(c_sub2ind) export(c_topntail) export(c_topntail_list) export(c_total_cable) -export(weighted_jaccard_dense_cpp) -export(weighted_jaccard_sparse_fill_cpp) +export(c_weighted_jaccard_dense) +export(c_weighted_jaccard_sparse) import(Rcpp) useDynLib(natcpp, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index 90fcc19..6474e92 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -135,22 +135,9 @@ c_EdgeListFromSegList <- function(L) { .Call(`_natcpp_c_EdgeListFromSegList`, L) } -#' Fill min_sums into a sparse pattern for weighted Jaccard similarity -#' -#' Given a dgCMatrix \code{x} and a sparsity pattern \code{pattern} (from -#' \code{crossprod(binarise(x))}), compute the sum of element-wise minima -#' for each pair of columns (or rows) that co-occur in at least one feature. -#' Intended for internal use by \code{coconat::jaccard_sim}. -#' -#' @param x A dgCMatrix (sparse column-compressed matrix) -#' @param pattern A dgCMatrix giving the output sparsity pattern -#' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare -#' rows -#' @return A numeric vector of min_sums aligned with the pattern's slot -#' structure -#' @export -weighted_jaccard_sparse_fill_cpp <- function(x, pattern, transpose = FALSE) { - .Call(`_natcpp_weighted_jaccard_sparse_fill_cpp`, x, pattern, transpose) +#' @noRd +weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE) { + .Call(`_natcpp_weighted_jaccard_sparse_fill`, x, pattern, transpose) } #' Dense weighted Jaccard similarity via C++ @@ -159,14 +146,14 @@ weighted_jaccard_sparse_fill_cpp <- function(x, pattern, transpose = FALSE) { #' an adaptive dense accumulation strategy. For small output matrices, uses a #' feature-oriented loop; for large outputs, switches to a column-oriented #' loop for better cache performance. Intended for internal use by -#' \code{coconat::jaccard_sim}. +#' \code{coconat::jaccard_sim}. See also \code{\link{c_weighted_jaccard_sparse}}. #' #' @param x A dgCMatrix (sparse column-compressed matrix) #' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare #' rows #' @return A dense numeric similarity matrix #' @export -weighted_jaccard_dense_cpp <- function(x, transpose = FALSE) { - .Call(`_natcpp_weighted_jaccard_dense_cpp`, x, transpose) +c_weighted_jaccard_dense <- function(x, transpose = FALSE) { + .Call(`_natcpp_c_weighted_jaccard_dense`, x, transpose) } diff --git a/R/weighted_jaccard.R b/R/weighted_jaccard.R new file mode 100644 index 0000000..eb4eb35 --- /dev/null +++ b/R/weighted_jaccard.R @@ -0,0 +1,51 @@ +#' Sparse weighted Jaccard similarity via C++ +#' +#' Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a +#' sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute +#' min-sums only for column (or row) pairs that share at least one non-zero +#' feature, then normalises to similarity. +#' +#' @param x A dgCMatrix (sparse column-compressed matrix) +#' @param transpose If \code{FALSE} (default), compare columns; if +#' \code{TRUE}, compare rows. +#' @return A sparse dgCMatrix similarity matrix +#' @export +#' @seealso \code{\link{c_weighted_jaccard_dense}} for the dense equivalent +#' @examples +#' \dontrun{ +#' library(Matrix) +#' m <- sparseMatrix(i = c(1,2,1,2,3,3), j = c(1,1,2,2,2,3), +#' x = c(4,2,1,3,3,1), dims = c(3,3)) +#' c_weighted_jaccard_sparse(m) +#' } +c_weighted_jaccard_sparse <- function(x, transpose = FALSE) { + crossfun <- if (transpose) Matrix::tcrossprod else Matrix::crossprod + n <- if (transpose) nrow(x) else ncol(x) + + if (length(x@x) == 0L) { + return(Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, + dims = c(n, n))) + } + + # Sparsity pattern from binarised crossprod + b <- x + b@x <- rep(1, length(b@x)) + A <- as(crossfun(b), "generalMatrix") + + # Column/row totals for denominator + totals <- if (transpose) Matrix::rowSums(x) else Matrix::colSums(x) + + # Fill in min_sums using C++ + A@x <- weighted_jaccard_sparse_fill(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 + + A +} diff --git a/man/weighted_jaccard_dense_cpp.Rd b/man/c_weighted_jaccard_dense.Rd similarity index 77% rename from man/weighted_jaccard_dense_cpp.Rd rename to man/c_weighted_jaccard_dense.Rd index 5269d34..d63a502 100644 --- a/man/weighted_jaccard_dense_cpp.Rd +++ b/man/c_weighted_jaccard_dense.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/RcppExports.R -\name{weighted_jaccard_dense_cpp} -\alias{weighted_jaccard_dense_cpp} +\name{c_weighted_jaccard_dense} +\alias{c_weighted_jaccard_dense} \title{Dense weighted Jaccard similarity via C++} \usage{ -weighted_jaccard_dense_cpp(x, transpose = FALSE) +c_weighted_jaccard_dense(x, transpose = FALSE) } \arguments{ \item{x}{A dgCMatrix (sparse column-compressed matrix)} @@ -20,5 +20,5 @@ Compute the full weighted Jaccard similarity matrix for a dgCMatrix using an adaptive dense accumulation strategy. For small output matrices, uses a feature-oriented loop; for large outputs, switches to a column-oriented loop for better cache performance. Intended for internal use by -\code{coconat::jaccard_sim}. +\code{coconat::jaccard_sim}. See also \code{\link{c_weighted_jaccard_sparse}}. } diff --git a/man/c_weighted_jaccard_sparse.Rd b/man/c_weighted_jaccard_sparse.Rd new file mode 100644 index 0000000..7f41e52 --- /dev/null +++ b/man/c_weighted_jaccard_sparse.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/weighted_jaccard.R +\name{c_weighted_jaccard_sparse} +\alias{c_weighted_jaccard_sparse} +\title{Sparse weighted Jaccard similarity via C++} +\usage{ +c_weighted_jaccard_sparse(x, transpose = FALSE) +} +\arguments{ +\item{x}{A dgCMatrix (sparse column-compressed matrix)} + +\item{transpose}{If \code{FALSE} (default), compare columns; if +\code{TRUE}, compare rows.} +} +\value{ +A sparse dgCMatrix similarity matrix +} +\description{ +Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a +sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute +min-sums only for column (or row) pairs that share at least one non-zero +feature, then normalises to similarity. +} +\examples{ +\dontrun{ +library(Matrix) +m <- sparseMatrix(i = c(1,2,1,2,3,3), j = c(1,1,2,2,2,3), + x = c(4,2,1,3,3,1), dims = c(3,3)) +c_weighted_jaccard_sparse(m) +} +} +\seealso{ +\code{\link{c_weighted_jaccard_dense}} for the dense equivalent +} diff --git a/man/weighted_jaccard_sparse_fill_cpp.Rd b/man/weighted_jaccard_sparse_fill_cpp.Rd deleted file mode 100644 index 00b597e..0000000 --- a/man/weighted_jaccard_sparse_fill_cpp.Rd +++ /dev/null @@ -1,26 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{weighted_jaccard_sparse_fill_cpp} -\alias{weighted_jaccard_sparse_fill_cpp} -\title{Fill min_sums into a sparse pattern for weighted Jaccard similarity} -\usage{ -weighted_jaccard_sparse_fill_cpp(x, pattern, transpose = FALSE) -} -\arguments{ -\item{x}{A dgCMatrix (sparse column-compressed matrix)} - -\item{pattern}{A dgCMatrix giving the output sparsity pattern} - -\item{transpose}{If \code{FALSE}, compare columns; if \code{TRUE}, compare -rows} -} -\value{ -A numeric vector of min_sums aligned with the pattern's slot - structure -} -\description{ -Given a dgCMatrix \code{x} and a sparsity pattern \code{pattern} (from -\code{crossprod(binarise(x))}), compute the sum of element-wise minima -for each pair of columns (or rows) that co-occur in at least one feature. -Intended for internal use by \code{coconat::jaccard_sim}. -} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 6db9ee8..a0e9c71 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -135,28 +135,28 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// weighted_jaccard_sparse_fill_cpp -NumericVector weighted_jaccard_sparse_fill_cpp(const S4& x, const S4& pattern, bool transpose); -RcppExport SEXP _natcpp_weighted_jaccard_sparse_fill_cpp(SEXP xSEXP, SEXP patternSEXP, SEXP transposeSEXP) { +// weighted_jaccard_sparse_fill +NumericVector weighted_jaccard_sparse_fill(const S4& x, const S4& pattern, bool transpose); +RcppExport SEXP _natcpp_weighted_jaccard_sparse_fill(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)); + rcpp_result_gen = Rcpp::wrap(weighted_jaccard_sparse_fill(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 _natcpp_weighted_jaccard_dense_cpp(SEXP xSEXP, SEXP transposeSEXP) { +// c_weighted_jaccard_dense +NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose); +RcppExport SEXP _natcpp_c_weighted_jaccard_dense(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)); + rcpp_result_gen = Rcpp::wrap(c_weighted_jaccard_dense(x, transpose)); return rcpp_result_gen; END_RCPP } @@ -172,8 +172,8 @@ static const R_CallMethodDef CallEntries[] = { {"_natcpp_c_topntail", (DL_FUNC) &_natcpp_c_topntail, 1}, {"_natcpp_c_topntail_list", (DL_FUNC) &_natcpp_c_topntail_list, 1}, {"_natcpp_c_EdgeListFromSegList", (DL_FUNC) &_natcpp_c_EdgeListFromSegList, 1}, - {"_natcpp_weighted_jaccard_sparse_fill_cpp", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill_cpp, 3}, - {"_natcpp_weighted_jaccard_dense_cpp", (DL_FUNC) &_natcpp_weighted_jaccard_dense_cpp, 2}, + {"_natcpp_weighted_jaccard_sparse_fill", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill, 3}, + {"_natcpp_c_weighted_jaccard_dense", (DL_FUNC) &_natcpp_c_weighted_jaccard_dense, 2}, {NULL, NULL, 0} }; diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp index 58f4fae..297bd32 100644 --- a/src/weighted_jaccard.cpp +++ b/src/weighted_jaccard.cpp @@ -4,22 +4,9 @@ using namespace Rcpp; -//' Fill min_sums into a sparse pattern for weighted Jaccard similarity -//' -//' Given a dgCMatrix \code{x} and a sparsity pattern \code{pattern} (from -//' \code{crossprod(binarise(x))}), compute the sum of element-wise minima -//' for each pair of columns (or rows) that co-occur in at least one feature. -//' Intended for internal use by \code{coconat::jaccard_sim}. -//' -//' @param x A dgCMatrix (sparse column-compressed matrix) -//' @param pattern A dgCMatrix giving the output sparsity pattern -//' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare -//' rows -//' @return A numeric vector of min_sums aligned with the pattern's slot -//' structure -//' @export +//' @noRd // [[Rcpp::export]] -NumericVector weighted_jaccard_sparse_fill_cpp( +NumericVector weighted_jaccard_sparse_fill( const S4& x, const S4& pattern, bool transpose = false) { // Original matrix slots (dgCMatrix: nr x nc) @@ -145,7 +132,7 @@ NumericVector weighted_jaccard_sparse_fill_cpp( //' an adaptive dense accumulation strategy. For small output matrices, uses a //' feature-oriented loop; for large outputs, switches to a column-oriented //' loop for better cache performance. Intended for internal use by -//' \code{coconat::jaccard_sim}. +//' \code{coconat::jaccard_sim}. See also \code{\link{c_weighted_jaccard_sparse}}. //' //' @param x A dgCMatrix (sparse column-compressed matrix) //' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare @@ -153,7 +140,7 @@ NumericVector weighted_jaccard_sparse_fill_cpp( //' @return A dense numeric similarity matrix //' @export // [[Rcpp::export]] -NumericMatrix weighted_jaccard_dense_cpp(const S4& x, bool transpose = false) { +NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose = false) { IntegerVector dims = x.slot("Dim"); IntegerVector xi = x.slot("i"); IntegerVector xp = x.slot("p"); diff --git a/tests/testthat/test-weighted-jaccard.R b/tests/testthat/test-weighted-jaccard.R index 207ab32..e1c1f11 100644 --- a/tests/testthat/test-weighted-jaccard.R +++ b/tests/testthat/test-weighted-jaccard.R @@ -1,5 +1,4 @@ -test_that("weighted_jaccard_dense_cpp computes correct similarity", { - skip_if_not_installed("Matrix") +test_that("c_weighted_jaccard_dense computes correct similarity", { m <- Matrix::sparseMatrix( i = c(1L, 2L, 1L, 2L, 3L, 3L), j = c(1L, 1L, 2L, 2L, 2L, 3L), @@ -20,12 +19,11 @@ test_that("weighted_jaccard_dense_cpp computes correct similarity", { } diag(ref) <- 1 - sim <- weighted_jaccard_dense_cpp(m, transpose = FALSE) + sim <- c_weighted_jaccard_dense(m, transpose = FALSE) expect_equal(sim, ref, tolerance = 1e-12) }) -test_that("weighted_jaccard_dense_cpp transpose works", { - skip_if_not_installed("Matrix") +test_that("c_weighted_jaccard_dense transpose works", { m <- Matrix::sparseMatrix( i = c(1L, 2L, 1L, 2L, 3L, 3L), j = c(1L, 1L, 2L, 2L, 2L, 3L), @@ -45,12 +43,11 @@ test_that("weighted_jaccard_dense_cpp transpose works", { } diag(ref_t) <- 1 - sim_t <- weighted_jaccard_dense_cpp(m, transpose = TRUE) + sim_t <- c_weighted_jaccard_dense(m, transpose = TRUE) expect_equal(sim_t, ref_t, tolerance = 1e-12) }) -test_that("weighted_jaccard_sparse_fill_cpp matches dense", { - skip_if_not_installed("Matrix") +test_that("c_weighted_jaccard_sparse matches dense", { m <- Matrix::sparseMatrix( i = c(1L, 2L, 1L, 2L, 3L, 3L), j = c(1L, 1L, 2L, 2L, 2L, 3L), @@ -58,25 +55,11 @@ test_that("weighted_jaccard_sparse_fill_cpp matches dense", { dims = c(3L, 3L) ) - # Dense reference - sim_dense <- weighted_jaccard_dense_cpp(m, transpose = FALSE) + sim_dense <- c_weighted_jaccard_dense(m, transpose = FALSE) + sim_sparse <- c_weighted_jaccard_sparse(m, transpose = FALSE) + expect_equal(as.matrix(sim_sparse), sim_dense, tolerance = 1e-12) - # Sparse fill: build pattern from binarised crossprod - b <- m - b@x <- rep(1, length(b@x)) - A <- methods::as(Matrix::crossprod(b), "generalMatrix") - totals <- Matrix::colSums(m) - - A@x <- weighted_jaccard_sparse_fill_cpp(m, A, transpose = FALSE) - - # Convert min_sums to similarity - 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 - - expect_equal(as.matrix(A), sim_dense, tolerance = 1e-12) + sim_dense_t <- c_weighted_jaccard_dense(m, transpose = TRUE) + sim_sparse_t <- c_weighted_jaccard_sparse(m, transpose = TRUE) + expect_equal(as.matrix(sim_sparse_t), sim_dense_t, tolerance = 1e-12) }) From 8eae444be743d9d9b6f3e16a0f7614e2f676f491 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Sat, 4 Apr 2026 22:56:29 +0100 Subject: [PATCH 4/8] switch to symmetric sparse output * this should halve memory usage --- R/weighted_jaccard.R | 22 +++++++++++++--------- man/c_weighted_jaccard_sparse.Rd | 7 ++++--- 2 files changed, 17 insertions(+), 12 deletions(-) diff --git a/R/weighted_jaccard.R b/R/weighted_jaccard.R index eb4eb35..7d72808 100644 --- a/R/weighted_jaccard.R +++ b/R/weighted_jaccard.R @@ -1,14 +1,15 @@ #' Sparse weighted Jaccard similarity via C++ #' #' Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a -#' sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute +#' symmetric sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute #' min-sums only for column (or row) pairs that share at least one non-zero -#' feature, then normalises to similarity. +#' feature, then normalises to similarity. Only the upper triangle is computed, +#' taking advantage of the symmetry of the Jaccard index. #' -#' @param x A dgCMatrix (sparse column-compressed matrix) +#' @param x A \code[Matrix]{Matrix} object of class dgCMatrix (sparse column-compressed matrix) #' @param transpose If \code{FALSE} (default), compare columns; if #' \code{TRUE}, compare rows. -#' @return A sparse dgCMatrix similarity matrix +#' @return A symmetric sparse dsCMatrix similarity matrix #' @export #' @seealso \code{\link{c_weighted_jaccard_dense}} for the dense equivalent #' @examples @@ -24,18 +25,18 @@ c_weighted_jaccard_sparse <- function(x, transpose = FALSE) { if (length(x@x) == 0L) { return(Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, - dims = c(n, n))) + dims = c(n, n), symmetric = TRUE)) } - # Sparsity pattern from binarised crossprod + # Sparsity pattern from binarised crossprod — dsCMatrix (upper triangle) b <- x b@x <- rep(1, length(b@x)) - A <- as(crossfun(b), "generalMatrix") + A <- crossfun(b) # Column/row totals for denominator totals <- if (transpose) Matrix::rowSums(x) else Matrix::colSums(x) - # Fill in min_sums using C++ + # Fill in min_sums using C++ (only upper triangle entries) A@x <- weighted_jaccard_sparse_fill(x, A, transpose = transpose) # Convert min_sums to similarity: sim = ms / (s_i + s_j - ms) @@ -45,7 +46,10 @@ c_weighted_jaccard_sparse <- function(x, transpose = FALSE) { nonzero <- denom != 0 A@x[nonzero] <- A@x[nonzero] / denom[nonzero] A@x[!nonzero] <- 0 - Matrix::diag(A) <- 1 + + # Set diagonal to 1 — diagonal entries are in the upper triangle pattern + diag_pos <- which(A@i + 1L == col_idx) + A@x[diag_pos] <- 1 A } diff --git a/man/c_weighted_jaccard_sparse.Rd b/man/c_weighted_jaccard_sparse.Rd index 7f41e52..efa9d3a 100644 --- a/man/c_weighted_jaccard_sparse.Rd +++ b/man/c_weighted_jaccard_sparse.Rd @@ -13,13 +13,14 @@ c_weighted_jaccard_sparse(x, transpose = FALSE) \code{TRUE}, compare rows.} } \value{ -A sparse dgCMatrix similarity matrix +A symmetric sparse dsCMatrix similarity matrix } \description{ Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a -sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute +symmetric sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute min-sums only for column (or row) pairs that share at least one non-zero -feature, then normalises to similarity. +feature, then normalises to similarity. Only the upper triangle is computed, +taking advantage of the symmetry of the Jaccard index. } \examples{ \dontrun{ From ff967b895152c184988443bc83ffc5206876fb57 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Sat, 4 Apr 2026 23:07:36 +0100 Subject: [PATCH 5/8] add progress to c_weighted_jaccard_sparse --- DESCRIPTION | 5 +++-- R/RcppExports.R | 4 ++-- R/weighted_jaccard.R | 9 ++++++--- man/c_weighted_jaccard_sparse.Rd | 5 ++++- src/RcppExports.cpp | 9 +++++---- src/weighted_jaccard.cpp | 8 +++++++- tests/testthat/test-weighted-jaccard.R | 4 ++-- 7 files changed, 29 insertions(+), 15 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 10934b6..9a5066a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,8 +24,9 @@ Imports: Suggests: spelling, testthat (>= 3.0.0) -LinkingTo: - Rcpp +LinkingTo: + Rcpp, + RcppProgress Config/testthat/edition: 3 Encoding: UTF-8 Language: en-GB diff --git a/R/RcppExports.R b/R/RcppExports.R index 6474e92..ae83b03 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -136,8 +136,8 @@ c_EdgeListFromSegList <- function(L) { } #' @noRd -weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE) { - .Call(`_natcpp_weighted_jaccard_sparse_fill`, x, pattern, transpose) +weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE, display_progress = TRUE) { + .Call(`_natcpp_weighted_jaccard_sparse_fill`, x, pattern, transpose, display_progress) } #' Dense weighted Jaccard similarity via C++ diff --git a/R/weighted_jaccard.R b/R/weighted_jaccard.R index 7d72808..f31ea2f 100644 --- a/R/weighted_jaccard.R +++ b/R/weighted_jaccard.R @@ -6,9 +6,11 @@ #' feature, then normalises to similarity. Only the upper triangle is computed, #' taking advantage of the symmetry of the Jaccard index. #' -#' @param x A \code[Matrix]{Matrix} object of class dgCMatrix (sparse column-compressed matrix) +#' @param x A dgCMatrix (sparse column-compressed matrix) #' @param transpose If \code{FALSE} (default), compare columns; if #' \code{TRUE}, compare rows. +#' @param display_progress Whether to show a text progress bar (default +#' \code{TRUE}). #' @return A symmetric sparse dsCMatrix similarity matrix #' @export #' @seealso \code{\link{c_weighted_jaccard_dense}} for the dense equivalent @@ -19,7 +21,7 @@ #' x = c(4,2,1,3,3,1), dims = c(3,3)) #' c_weighted_jaccard_sparse(m) #' } -c_weighted_jaccard_sparse <- function(x, transpose = FALSE) { +c_weighted_jaccard_sparse <- function(x, transpose = FALSE, display_progress = TRUE) { crossfun <- if (transpose) Matrix::tcrossprod else Matrix::crossprod n <- if (transpose) nrow(x) else ncol(x) @@ -37,7 +39,8 @@ c_weighted_jaccard_sparse <- function(x, transpose = FALSE) { totals <- if (transpose) Matrix::rowSums(x) else Matrix::colSums(x) # Fill in min_sums using C++ (only upper triangle entries) - A@x <- weighted_jaccard_sparse_fill(x, A, transpose = transpose) + A@x <- weighted_jaccard_sparse_fill(x, A, transpose = transpose, + display_progress = display_progress) # Convert min_sums to similarity: sim = ms / (s_i + s_j - ms) col_idx <- rep(seq_along(diff(A@p)), diff(A@p)) diff --git a/man/c_weighted_jaccard_sparse.Rd b/man/c_weighted_jaccard_sparse.Rd index efa9d3a..e744c12 100644 --- a/man/c_weighted_jaccard_sparse.Rd +++ b/man/c_weighted_jaccard_sparse.Rd @@ -4,13 +4,16 @@ \alias{c_weighted_jaccard_sparse} \title{Sparse weighted Jaccard similarity via C++} \usage{ -c_weighted_jaccard_sparse(x, transpose = FALSE) +c_weighted_jaccard_sparse(x, transpose = FALSE, display_progress = TRUE) } \arguments{ \item{x}{A dgCMatrix (sparse column-compressed matrix)} \item{transpose}{If \code{FALSE} (default), compare columns; if \code{TRUE}, compare rows.} + +\item{display_progress}{Whether to show a text progress bar (default +\code{TRUE}).} } \value{ A symmetric sparse dsCMatrix similarity matrix diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index a0e9c71..48bf24e 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -136,15 +136,16 @@ BEGIN_RCPP END_RCPP } // weighted_jaccard_sparse_fill -NumericVector weighted_jaccard_sparse_fill(const S4& x, const S4& pattern, bool transpose); -RcppExport SEXP _natcpp_weighted_jaccard_sparse_fill(SEXP xSEXP, SEXP patternSEXP, SEXP transposeSEXP) { +NumericVector weighted_jaccard_sparse_fill(const S4& x, const S4& pattern, bool transpose, bool display_progress); +RcppExport SEXP _natcpp_weighted_jaccard_sparse_fill(SEXP xSEXP, SEXP patternSEXP, SEXP transposeSEXP, SEXP display_progressSEXP) { 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(x, pattern, transpose)); + Rcpp::traits::input_parameter< bool >::type display_progress(display_progressSEXP); + rcpp_result_gen = Rcpp::wrap(weighted_jaccard_sparse_fill(x, pattern, transpose, display_progress)); return rcpp_result_gen; END_RCPP } @@ -172,7 +173,7 @@ static const R_CallMethodDef CallEntries[] = { {"_natcpp_c_topntail", (DL_FUNC) &_natcpp_c_topntail, 1}, {"_natcpp_c_topntail_list", (DL_FUNC) &_natcpp_c_topntail_list, 1}, {"_natcpp_c_EdgeListFromSegList", (DL_FUNC) &_natcpp_c_EdgeListFromSegList, 1}, - {"_natcpp_weighted_jaccard_sparse_fill", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill, 3}, + {"_natcpp_weighted_jaccard_sparse_fill", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill, 4}, {"_natcpp_c_weighted_jaccard_dense", (DL_FUNC) &_natcpp_c_weighted_jaccard_dense, 2}, {NULL, NULL, 0} }; diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp index 297bd32..affcc0e 100644 --- a/src/weighted_jaccard.cpp +++ b/src/weighted_jaccard.cpp @@ -1,13 +1,15 @@ #include #include #include +#include using namespace Rcpp; //' @noRd // [[Rcpp::export]] NumericVector weighted_jaccard_sparse_fill( - const S4& x, const S4& pattern, bool transpose = false) { + const S4& x, const S4& pattern, bool transpose = false, + bool display_progress = true) { // Original matrix slots (dgCMatrix: nr x nc) IntegerVector dims = x.slot("Dim"); @@ -98,8 +100,12 @@ NumericVector weighted_jaccard_sparse_fill( // Populated once per output column, then cleared. std::vector row_to_pos(ncomp, -1); + Progress prog(ncomp, display_progress); + // Iterate over output columns for (int cb = 0; cb < ncomp; ++cb) { + if (Progress::check_abort()) return NumericVector(annz); + prog.increment(); // Populate row_to_pos for pattern column cb for (int idx = Ap[cb]; idx < Ap[cb + 1]; ++idx) row_to_pos[Ai[idx]] = idx; diff --git a/tests/testthat/test-weighted-jaccard.R b/tests/testthat/test-weighted-jaccard.R index e1c1f11..d595252 100644 --- a/tests/testthat/test-weighted-jaccard.R +++ b/tests/testthat/test-weighted-jaccard.R @@ -56,10 +56,10 @@ test_that("c_weighted_jaccard_sparse matches dense", { ) sim_dense <- c_weighted_jaccard_dense(m, transpose = FALSE) - sim_sparse <- c_weighted_jaccard_sparse(m, transpose = FALSE) + sim_sparse <- c_weighted_jaccard_sparse(m, transpose = FALSE, display_progress = FALSE) expect_equal(as.matrix(sim_sparse), sim_dense, tolerance = 1e-12) sim_dense_t <- c_weighted_jaccard_dense(m, transpose = TRUE) - sim_sparse_t <- c_weighted_jaccard_sparse(m, transpose = TRUE) + sim_sparse_t <- c_weighted_jaccard_sparse(m, transpose = TRUE, display_progress = FALSE) expect_equal(as.matrix(sim_sparse_t), sim_dense_t, tolerance = 1e-12) }) From 4abf9ef380096c252fd2189a629f05fe80d4332f Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Sun, 5 Apr 2026 23:10:47 +0100 Subject: [PATCH 6/8] threaded versions of weighted_jaccard with progress * 4 threads are a good compromise in tests --- DESCRIPTION | 2 +- R/RcppExports.R | 8 +- R/weighted_jaccard.R | 8 +- man/c_weighted_jaccard_dense.Rd | 2 +- man/c_weighted_jaccard_sparse.Rd | 10 +- src/Makevars | 1 + src/RcppExports.cpp | 19 ++-- src/weighted_jaccard.cpp | 167 +++++++++++++++++++++++-------- 8 files changed, 159 insertions(+), 58 deletions(-) create mode 100644 src/Makevars diff --git a/DESCRIPTION b/DESCRIPTION index 9a5066a..7378e8c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -26,7 +26,7 @@ Suggests: testthat (>= 3.0.0) LinkingTo: Rcpp, - RcppProgress + RcppThread Config/testthat/edition: 3 Encoding: UTF-8 Language: en-GB diff --git a/R/RcppExports.R b/R/RcppExports.R index ae83b03..3f25b30 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -136,8 +136,8 @@ c_EdgeListFromSegList <- function(L) { } #' @noRd -weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE, display_progress = TRUE) { - .Call(`_natcpp_weighted_jaccard_sparse_fill`, x, pattern, transpose, display_progress) +weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE, display_progress = TRUE, threads = 4L) { + .Call(`_natcpp_weighted_jaccard_sparse_fill`, x, pattern, transpose, display_progress, threads) } #' Dense weighted Jaccard similarity via C++ @@ -153,7 +153,7 @@ weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE, display_ #' rows #' @return A dense numeric similarity matrix #' @export -c_weighted_jaccard_dense <- function(x, transpose = FALSE) { - .Call(`_natcpp_c_weighted_jaccard_dense`, x, transpose) +c_weighted_jaccard_dense <- function(x, transpose = FALSE, threads = 4L) { + .Call(`_natcpp_c_weighted_jaccard_dense`, x, transpose, threads) } diff --git a/R/weighted_jaccard.R b/R/weighted_jaccard.R index f31ea2f..e9c3bee 100644 --- a/R/weighted_jaccard.R +++ b/R/weighted_jaccard.R @@ -11,6 +11,8 @@ #' \code{TRUE}, compare rows. #' @param display_progress Whether to show a text progress bar (default #' \code{TRUE}). +#' @param threads Number of threads for parallel computation (default 4). +#' Set to 0 to use all available cores. #' @return A symmetric sparse dsCMatrix similarity matrix #' @export #' @seealso \code{\link{c_weighted_jaccard_dense}} for the dense equivalent @@ -21,7 +23,8 @@ #' x = c(4,2,1,3,3,1), dims = c(3,3)) #' c_weighted_jaccard_sparse(m) #' } -c_weighted_jaccard_sparse <- function(x, transpose = FALSE, display_progress = TRUE) { +c_weighted_jaccard_sparse <- function(x, transpose = FALSE, display_progress = TRUE, + threads = 4L) { crossfun <- if (transpose) Matrix::tcrossprod else Matrix::crossprod n <- if (transpose) nrow(x) else ncol(x) @@ -40,7 +43,8 @@ c_weighted_jaccard_sparse <- function(x, transpose = FALSE, display_progress = T # Fill in min_sums using C++ (only upper triangle entries) A@x <- weighted_jaccard_sparse_fill(x, A, transpose = transpose, - display_progress = display_progress) + display_progress = display_progress, + threads = threads) # Convert min_sums to similarity: sim = ms / (s_i + s_j - ms) col_idx <- rep(seq_along(diff(A@p)), diff(A@p)) diff --git a/man/c_weighted_jaccard_dense.Rd b/man/c_weighted_jaccard_dense.Rd index d63a502..e9aaa0a 100644 --- a/man/c_weighted_jaccard_dense.Rd +++ b/man/c_weighted_jaccard_dense.Rd @@ -4,7 +4,7 @@ \alias{c_weighted_jaccard_dense} \title{Dense weighted Jaccard similarity via C++} \usage{ -c_weighted_jaccard_dense(x, transpose = FALSE) +c_weighted_jaccard_dense(x, transpose = FALSE, threads = 4L) } \arguments{ \item{x}{A dgCMatrix (sparse column-compressed matrix)} diff --git a/man/c_weighted_jaccard_sparse.Rd b/man/c_weighted_jaccard_sparse.Rd index e744c12..11b3d52 100644 --- a/man/c_weighted_jaccard_sparse.Rd +++ b/man/c_weighted_jaccard_sparse.Rd @@ -4,7 +4,12 @@ \alias{c_weighted_jaccard_sparse} \title{Sparse weighted Jaccard similarity via C++} \usage{ -c_weighted_jaccard_sparse(x, transpose = FALSE, display_progress = TRUE) +c_weighted_jaccard_sparse( + x, + transpose = FALSE, + display_progress = TRUE, + threads = 4L +) } \arguments{ \item{x}{A dgCMatrix (sparse column-compressed matrix)} @@ -14,6 +19,9 @@ c_weighted_jaccard_sparse(x, transpose = FALSE, display_progress = TRUE) \item{display_progress}{Whether to show a text progress bar (default \code{TRUE}).} + +\item{threads}{Number of threads for parallel computation (default 4). +Set to 0 to use all available cores.} } \value{ A symmetric sparse dsCMatrix similarity matrix diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..054912b --- /dev/null +++ b/src/Makevars @@ -0,0 +1 @@ +PKG_LIBS = $(shell "${R_HOME}/bin/Rscript" -e "RcppThread::LdFlags()") diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 48bf24e..2aa76ba 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1,6 +1,7 @@ // Generated by using Rcpp::compileAttributes() -> do not edit by hand // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +#include #include using namespace Rcpp; @@ -136,8 +137,8 @@ BEGIN_RCPP END_RCPP } // weighted_jaccard_sparse_fill -NumericVector weighted_jaccard_sparse_fill(const S4& x, const S4& pattern, bool transpose, bool display_progress); -RcppExport SEXP _natcpp_weighted_jaccard_sparse_fill(SEXP xSEXP, SEXP patternSEXP, SEXP transposeSEXP, SEXP display_progressSEXP) { +NumericVector weighted_jaccard_sparse_fill(const S4& x, const S4& pattern, bool transpose, bool display_progress, int threads); +RcppExport SEXP _natcpp_weighted_jaccard_sparse_fill(SEXP xSEXP, SEXP patternSEXP, SEXP transposeSEXP, SEXP display_progressSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -145,19 +146,21 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const S4& >::type pattern(patternSEXP); Rcpp::traits::input_parameter< bool >::type transpose(transposeSEXP); Rcpp::traits::input_parameter< bool >::type display_progress(display_progressSEXP); - rcpp_result_gen = Rcpp::wrap(weighted_jaccard_sparse_fill(x, pattern, transpose, display_progress)); + Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); + rcpp_result_gen = Rcpp::wrap(weighted_jaccard_sparse_fill(x, pattern, transpose, display_progress, threads)); return rcpp_result_gen; END_RCPP } // c_weighted_jaccard_dense -NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose); -RcppExport SEXP _natcpp_c_weighted_jaccard_dense(SEXP xSEXP, SEXP transposeSEXP) { +NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose, int threads); +RcppExport SEXP _natcpp_c_weighted_jaccard_dense(SEXP xSEXP, SEXP transposeSEXP, SEXP threadsSEXP) { 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(c_weighted_jaccard_dense(x, transpose)); + Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); + rcpp_result_gen = Rcpp::wrap(c_weighted_jaccard_dense(x, transpose, threads)); return rcpp_result_gen; END_RCPP } @@ -173,8 +176,8 @@ static const R_CallMethodDef CallEntries[] = { {"_natcpp_c_topntail", (DL_FUNC) &_natcpp_c_topntail, 1}, {"_natcpp_c_topntail_list", (DL_FUNC) &_natcpp_c_topntail_list, 1}, {"_natcpp_c_EdgeListFromSegList", (DL_FUNC) &_natcpp_c_EdgeListFromSegList, 1}, - {"_natcpp_weighted_jaccard_sparse_fill", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill, 4}, - {"_natcpp_c_weighted_jaccard_dense", (DL_FUNC) &_natcpp_c_weighted_jaccard_dense, 2}, + {"_natcpp_weighted_jaccard_sparse_fill", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill, 5}, + {"_natcpp_c_weighted_jaccard_dense", (DL_FUNC) &_natcpp_c_weighted_jaccard_dense, 3}, {NULL, NULL, 0} }; diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp index affcc0e..ab8bf0b 100644 --- a/src/weighted_jaccard.cpp +++ b/src/weighted_jaccard.cpp @@ -1,15 +1,57 @@ #include #include +#include +#include #include -#include +#include using namespace Rcpp; +// Progress bar that suppresses output for the first `delaySecs` seconds, +// then shows a bar with time estimate (using RcppThread's Rcout). +class DelayedProgressBar : public RcppThread::ProgressPrinter { +public: + DelayedProgressBar(size_t numIt, size_t printEvery, size_t delaySecs, + bool display = true) + : ProgressPrinter(numIt, printEvery), delaySecs_(delaySecs), + display_(display) {} + +private: + size_t delaySecs_; + bool display_; + bool started_{false}; + + void printProgress() override { + if (!display_ || isDone_) return; + if (it_ == numIt_) isDone_ = true; + + using namespace std::chrono; + auto elapsed = duration(steady_clock::now() - startTime_).count(); + if (!started_ && elapsed < delaySecs_) return; + started_ = true; + + double pct = std::round(it_ * 100.0 / numIt_); + std::ostringstream msg; + msg << "\rComputing: " << makeBar(pct) << progressString(); + RcppThread::Rcout << msg.str(); + } + + std::string makeBar(size_t pct, size_t numBars = 40) { + std::ostringstream msg; + msg << "["; + size_t i = 0; + for (; i < pct / 100.0 * numBars; i++) msg << "="; + for (; i < numBars; i++) msg << " "; + msg << "] "; + return msg.str(); + } +}; + //' @noRd // [[Rcpp::export]] NumericVector weighted_jaccard_sparse_fill( const S4& x, const S4& pattern, bool transpose = false, - bool display_progress = true) { + bool display_progress = true, int threads = 4) { // Original matrix slots (dgCMatrix: nr x nc) IntegerVector dims = x.slot("Dim"); @@ -96,38 +138,64 @@ NumericVector weighted_jaccard_sparse_fill( } } - // 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); - - Progress prog(ncomp, display_progress); + // Raw pointer for thread-safe access (Rcpp NumericVector not thread-safe) + double* out_ptr = out.begin(); + const int* Ai_ptr = &Ai[0]; + const int* Ap_ptr = &Ap[0]; + const int* item_feats_ptr = item_feats.data(); + const double* item_vals_ptr = item_vals.data(); + const int* feat_items_ptr = feat_items.data(); + const double* feat_vals_ptr = feat_vals.data(); + const int* feat_offsets_ptr = feat_offsets.data(); + const int* item_offsets_ptr = item_offsets.data(); + + DelayedProgressBar bar(ncomp, 1, 2, display_progress); + + const size_t nThreads = (threads > 0) ? static_cast(threads) + : std::thread::hardware_concurrency(); + + // Pre-allocate per-thread row_to_pos vectors (freed when function exits) + std::vector> per_thread_r2p(nThreads, std::vector(ncomp, -1)); + + // Thread-local index into per_thread_r2p, assigned once per thread via atomic + static thread_local int tl_tid = -1; + static thread_local int tl_gen = -1; + static std::atomic generation{0}; + const int cur_gen = ++generation; + std::atomic next_tid{0}; + + // Parallel over output columns — each column writes to its own slice of out[] + RcppThread::parallelFor(0, ncomp, [&](int cb) { + if (tl_gen != cur_gen) { + tl_tid = next_tid++; + tl_gen = cur_gen; + } + auto& row_to_pos = per_thread_r2p[tl_tid]; - // Iterate over output columns - for (int cb = 0; cb < ncomp; ++cb) { - if (Progress::check_abort()) return NumericVector(annz); - prog.increment(); // 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 (int idx = Ap_ptr[cb]; idx < Ap_ptr[cb + 1]; ++idx) + row_to_pos[Ai_ptr[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 (int fi = item_offsets_ptr[cb]; fi < item_offsets_ptr[cb + 1]; ++fi) { + const int f = item_feats_ptr[fi]; + const double vb = item_vals_ptr[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]; + // For each other item ca also in this feature, accumulate min + for (int idx = feat_offsets_ptr[f]; idx < feat_offsets_ptr[f + 1]; ++idx) { + const int ca = feat_items_ptr[idx]; const int pos = row_to_pos[ca]; if (pos >= 0) - out[pos] += std::min(feat_vals[idx], vb); + out_ptr[pos] += std::min(feat_vals_ptr[idx], vb); } } // Clear row_to_pos - for (int idx = Ap[cb]; idx < Ap[cb + 1]; ++idx) - row_to_pos[Ai[idx]] = -1; - } + for (int idx = Ap_ptr[cb]; idx < Ap_ptr[cb + 1]; ++idx) + row_to_pos[Ai_ptr[idx]] = -1; + + ++bar; + }, nThreads); return out; } @@ -146,7 +214,8 @@ NumericVector weighted_jaccard_sparse_fill( //' @return A dense numeric similarity matrix //' @export // [[Rcpp::export]] -NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose = false) { +NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose = false, + int threads = 4) { IntegerVector dims = x.slot("Dim"); IntegerVector xi = x.slot("i"); IntegerVector xp = x.slot("p"); @@ -235,16 +304,26 @@ NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose = false) { } } - 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); + const size_t nT = (threads > 0) ? static_cast(threads) + : std::thread::hardware_concurrency(); + double* out_ptr = &out[0]; + const int* feat_items_ptr = feat_items.data(); + const double* feat_vals_ptr = feat_vals.data(); + const int* feat_offsets_ptr = feat_offsets.data(); + const int* item_feats_ptr = item_feats.data(); + const double* item_vals_ptr = item_vals.data(); + const int* item_offsets_ptr = item_offsets.data(); + + RcppThread::parallelFor(0, ncomp, [&](int cb) { + double* col = out_ptr + static_cast(cb) * ncomp; + for (int fi = item_offsets_ptr[cb]; fi < item_offsets_ptr[cb + 1]; ++fi) { + const int f = item_feats_ptr[fi]; + const double vb = item_vals_ptr[fi]; + for (int idx = feat_offsets_ptr[f]; idx < feat_offsets_ptr[f + 1]; ++idx) { + col[feat_items_ptr[idx]] += std::min(feat_vals_ptr[idx], vb); } } - } + }, nT); } else { // Feature-oriented: iterate features, scatter into output for (int f = 0; f < nfeat; ++f) { @@ -267,17 +346,23 @@ NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose = false) { } // 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 size_t nT = (threads > 0) ? static_cast(threads) + : std::thread::hardware_concurrency(); + double* out_ptr = &out[0]; + const double* totals_ptr = totals.data(); + RcppThread::parallelFor(0, ncomp, [=](int cb) { + for (int ca = 0; ca < ncomp; ++ca) { 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; + if (ca == cb) { + out_ptr[pos] = 1.0; + } else { + const double ms = out_ptr[pos]; + const double denom = totals_ptr[ca] + totals_ptr[cb] - ms; + out_ptr[pos] = (denom > 0.0) ? ms / denom : 0.0; + } } - } + }, nT); } return out; From 4f75a012afae316f00e58fd48cb1c889fc1f8b1d Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Wed, 24 Jun 2026 18:05:43 +0100 Subject: [PATCH 7/8] Add triangle/distance options and portable Makevars MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit c_weighted_jaccard_sparse() gains: - triangle=FALSE: return dgCMatrix (general) by default, dsCMatrix on TRUE - distance=TRUE: return 1 - similarity (with warning about densification) Companion changes: - importFrom(methods, as); methods added to Imports - expanded test coverage (+70 lines) - src/Makevars: replace $(shell Rscript ... RcppThread::LdFlags()) with $(SHLIB_OPENMP_CXXFLAGS) — drops GNU make extension flagged by R CMD check WARNING, while still pulling in pthread/OpenMP per platform --- DESCRIPTION | 1 + NAMESPACE | 1 + R/RcppExports.R | 16 +- R/weighted_jaccard.R | 42 +++- man/c_weighted_jaccard_dense.Rd | 23 +- man/c_weighted_jaccard_sparse.Rd | 23 +- src/Makevars | 2 +- src/RcppExports.cpp | 10 +- src/weighted_jaccard.cpp | 307 +++++++++++++++++-------- tests/testthat/test-weighted-jaccard.R | 70 ++++++ 10 files changed, 368 insertions(+), 127 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7378e8c..5bca7b7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -20,6 +20,7 @@ URL: https://github.com/natverse/natcpp, https://natverse.org/natcpp/ BugReports: https://github.com/natverse/natcpp/issues Imports: Matrix, + methods, Rcpp (>= 1.0.6) Suggests: spelling, diff --git a/NAMESPACE b/NAMESPACE index 355a815..b0708b5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,4 +13,5 @@ export(c_total_cable) export(c_weighted_jaccard_dense) export(c_weighted_jaccard_sparse) import(Rcpp) +importFrom(methods,as) useDynLib(natcpp, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index 3f25b30..8ea311b 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -145,15 +145,21 @@ weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE, display_ #' Compute the full weighted Jaccard similarity matrix for a dgCMatrix using #' an adaptive dense accumulation strategy. For small output matrices, uses a #' feature-oriented loop; for large outputs, switches to a column-oriented -#' loop for better cache performance. Intended for internal use by -#' \code{coconat::jaccard_sim}. See also \code{\link{c_weighted_jaccard_sparse}}. +#' loop for better cache performance. #' #' @param x A dgCMatrix (sparse column-compressed matrix) #' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare #' rows -#' @return A dense numeric similarity matrix +#' @param threads Number of threads (default 4). Set to 0 for all cores. +#' @param triangle If \code{TRUE}, return only the lower triangle as a flat +#' numeric vector in \code{\link{dist}} layout. If \code{FALSE} (default), +#' return a full square matrix. +#' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) +#' instead of similarity. Default \code{FALSE}. +#' @return A dense numeric similarity matrix, or a numeric vector in dist +#' layout when \code{triangle = TRUE}. #' @export -c_weighted_jaccard_dense <- function(x, transpose = FALSE, threads = 4L) { - .Call(`_natcpp_c_weighted_jaccard_dense`, x, transpose, threads) +c_weighted_jaccard_dense <- function(x, transpose = FALSE, threads = 4L, triangle = FALSE, distance = FALSE) { + .Call(`_natcpp_c_weighted_jaccard_dense`, x, transpose, threads, triangle, distance) } diff --git a/R/weighted_jaccard.R b/R/weighted_jaccard.R index e9c3bee..6b73b15 100644 --- a/R/weighted_jaccard.R +++ b/R/weighted_jaccard.R @@ -1,10 +1,10 @@ #' Sparse weighted Jaccard similarity via C++ #' #' Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a -#' symmetric sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute -#' min-sums only for column (or row) pairs that share at least one non-zero -#' feature, then normalises to similarity. Only the upper triangle is computed, -#' taking advantage of the symmetry of the Jaccard index. +#' sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute min-sums +#' only for column (or row) pairs that share at least one non-zero feature, then +#' normalises to similarity. Only the upper triangle is computed, taking +#' advantage of the symmetry of the Jaccard index. #' #' @param x A dgCMatrix (sparse column-compressed matrix) #' @param transpose If \code{FALSE} (default), compare columns; if @@ -13,7 +13,15 @@ #' \code{TRUE}). #' @param threads Number of threads for parallel computation (default 4). #' Set to 0 to use all available cores. -#' @return A symmetric sparse dsCMatrix similarity matrix +#' @param triangle If \code{TRUE}, return a symmetric \code{dsCMatrix} +#' (upper triangle only). If \code{FALSE} (default), return a general +#' \code{dgCMatrix}. +#' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) +#' instead of similarity. Default \code{FALSE}. A warning is issued since +#' sparse distance matrices are typically dense. +#' @return A sparse similarity (or distance) matrix: \code{dsCMatrix} when +#' \code{triangle = TRUE}, \code{dgCMatrix} otherwise. +#' @importFrom methods as #' @export #' @seealso \code{\link{c_weighted_jaccard_dense}} for the dense equivalent #' @examples @@ -24,13 +32,22 @@ #' c_weighted_jaccard_sparse(m) #' } c_weighted_jaccard_sparse <- function(x, transpose = FALSE, display_progress = TRUE, - threads = 4L) { + threads = 4L, triangle = FALSE, + distance = FALSE) { + if (distance) + warning("distance=TRUE with sparse output produces a mostly-dense matrix; ", + "consider using c_weighted_jaccard_dense() with triangle=TRUE instead.") + crossfun <- if (transpose) Matrix::tcrossprod else Matrix::crossprod n <- if (transpose) nrow(x) else ncol(x) if (length(x@x) == 0L) { - return(Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, - dims = c(n, n), symmetric = TRUE)) + sim_val <- if (distance) 0 else 1 + m <- Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = sim_val, + dims = c(n, n), symmetric = triangle) + if (!triangle) + m <- as(m, "generalMatrix") + return(m) } # Sparsity pattern from binarised crossprod — dsCMatrix (upper triangle) @@ -54,9 +71,16 @@ c_weighted_jaccard_sparse <- function(x, transpose = FALSE, display_progress = T A@x[nonzero] <- A@x[nonzero] / denom[nonzero] A@x[!nonzero] <- 0 - # Set diagonal to 1 — diagonal entries are in the upper triangle pattern + # Set diagonal to 1 for similarity diag_pos <- which(A@i + 1L == col_idx) A@x[diag_pos] <- 1 + if (distance) { + A@x <- 1 - A@x # diagonal 1→0, off-diag sim→1-sim + } + + if (!triangle) + A <- as(A, "generalMatrix") + A } diff --git a/man/c_weighted_jaccard_dense.Rd b/man/c_weighted_jaccard_dense.Rd index e9aaa0a..1d074e0 100644 --- a/man/c_weighted_jaccard_dense.Rd +++ b/man/c_weighted_jaccard_dense.Rd @@ -4,21 +4,36 @@ \alias{c_weighted_jaccard_dense} \title{Dense weighted Jaccard similarity via C++} \usage{ -c_weighted_jaccard_dense(x, transpose = FALSE, threads = 4L) +c_weighted_jaccard_dense( + x, + transpose = FALSE, + threads = 4L, + triangle = FALSE, + distance = FALSE +) } \arguments{ \item{x}{A dgCMatrix (sparse column-compressed matrix)} \item{transpose}{If \code{FALSE}, compare columns; if \code{TRUE}, compare rows} + +\item{threads}{Number of threads (default 4). Set to 0 for all cores.} + +\item{triangle}{If \code{TRUE}, return only the lower triangle as a flat +numeric vector in \code{\link{dist}} layout. If \code{FALSE} (default), +return a full square matrix.} + +\item{distance}{If \code{TRUE}, return distance (\code{1 - similarity}) +instead of similarity. Default \code{FALSE}.} } \value{ -A dense numeric similarity matrix +A dense numeric similarity matrix, or a numeric vector in dist + layout when \code{triangle = TRUE}. } \description{ Compute the full weighted Jaccard similarity matrix for a dgCMatrix using an adaptive dense accumulation strategy. For small output matrices, uses a feature-oriented loop; for large outputs, switches to a column-oriented -loop for better cache performance. Intended for internal use by -\code{coconat::jaccard_sim}. See also \code{\link{c_weighted_jaccard_sparse}}. +loop for better cache performance. } diff --git a/man/c_weighted_jaccard_sparse.Rd b/man/c_weighted_jaccard_sparse.Rd index 11b3d52..ffcefe3 100644 --- a/man/c_weighted_jaccard_sparse.Rd +++ b/man/c_weighted_jaccard_sparse.Rd @@ -8,7 +8,9 @@ c_weighted_jaccard_sparse( x, transpose = FALSE, display_progress = TRUE, - threads = 4L + threads = 4L, + triangle = FALSE, + distance = FALSE ) } \arguments{ @@ -22,16 +24,25 @@ c_weighted_jaccard_sparse( \item{threads}{Number of threads for parallel computation (default 4). Set to 0 to use all available cores.} + +\item{triangle}{If \code{TRUE}, return a symmetric \code{dsCMatrix} +(upper triangle only). If \code{FALSE} (default), return a general +\code{dgCMatrix}.} + +\item{distance}{If \code{TRUE}, return distance (\code{1 - similarity}) +instead of similarity. Default \code{FALSE}. A warning is issued since +sparse distance matrices are typically dense.} } \value{ -A symmetric sparse dsCMatrix similarity matrix +A sparse similarity (or distance) matrix: \code{dsCMatrix} when + \code{triangle = TRUE}, \code{dgCMatrix} otherwise. } \description{ Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a -symmetric sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute -min-sums only for column (or row) pairs that share at least one non-zero -feature, then normalises to similarity. Only the upper triangle is computed, -taking advantage of the symmetry of the Jaccard index. +sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute min-sums +only for column (or row) pairs that share at least one non-zero feature, then +normalises to similarity. Only the upper triangle is computed, taking +advantage of the symmetry of the Jaccard index. } \examples{ \dontrun{ diff --git a/src/Makevars b/src/Makevars index 054912b..4c6d109 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1 +1 @@ -PKG_LIBS = $(shell "${R_HOME}/bin/Rscript" -e "RcppThread::LdFlags()") +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 2aa76ba..e6ff2f2 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -152,15 +152,17 @@ BEGIN_RCPP END_RCPP } // c_weighted_jaccard_dense -NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose, int threads); -RcppExport SEXP _natcpp_c_weighted_jaccard_dense(SEXP xSEXP, SEXP transposeSEXP, SEXP threadsSEXP) { +SEXP c_weighted_jaccard_dense(const S4& x, bool transpose, int threads, bool triangle, bool distance); +RcppExport SEXP _natcpp_c_weighted_jaccard_dense(SEXP xSEXP, SEXP transposeSEXP, SEXP threadsSEXP, SEXP triangleSEXP, SEXP distanceSEXP) { 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::traits::input_parameter< int >::type threads(threadsSEXP); - rcpp_result_gen = Rcpp::wrap(c_weighted_jaccard_dense(x, transpose, threads)); + Rcpp::traits::input_parameter< bool >::type triangle(triangleSEXP); + Rcpp::traits::input_parameter< bool >::type distance(distanceSEXP); + rcpp_result_gen = Rcpp::wrap(c_weighted_jaccard_dense(x, transpose, threads, triangle, distance)); return rcpp_result_gen; END_RCPP } @@ -177,7 +179,7 @@ static const R_CallMethodDef CallEntries[] = { {"_natcpp_c_topntail_list", (DL_FUNC) &_natcpp_c_topntail_list, 1}, {"_natcpp_c_EdgeListFromSegList", (DL_FUNC) &_natcpp_c_EdgeListFromSegList, 1}, {"_natcpp_weighted_jaccard_sparse_fill", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill, 5}, - {"_natcpp_c_weighted_jaccard_dense", (DL_FUNC) &_natcpp_c_weighted_jaccard_dense, 3}, + {"_natcpp_c_weighted_jaccard_dense", (DL_FUNC) &_natcpp_c_weighted_jaccard_dense, 5}, {NULL, NULL, 0} }; diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp index ab8bf0b..4ae611e 100644 --- a/src/weighted_jaccard.cpp +++ b/src/weighted_jaccard.cpp @@ -200,41 +200,31 @@ NumericVector weighted_jaccard_sparse_fill( return out; } -//' Dense weighted Jaccard similarity via C++ -//' -//' Compute the full weighted Jaccard similarity matrix for a dgCMatrix using -//' an adaptive dense accumulation strategy. For small output matrices, uses a -//' feature-oriented loop; for large outputs, switches to a column-oriented -//' loop for better cache performance. Intended for internal use by -//' \code{coconat::jaccard_sim}. See also \code{\link{c_weighted_jaccard_sparse}}. -//' -//' @param x A dgCMatrix (sparse column-compressed matrix) -//' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare -//' rows -//' @return A dense numeric similarity matrix -//' @export -// [[Rcpp::export]] -NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose = false, - int threads = 4) { - 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; +// Helper: build feature CSR and totals from dgCMatrix slots. +// Shared by both dense matrix and dense triangle functions. +struct FeatureCSR { + std::vector totals; + std::vector feat_offsets; + std::vector feat_items; + std::vector feat_vals; + int ncomp; + int nfeat; + int total_nnz; +}; - // Totals per compared item - std::vector totals(ncomp, 0.0); +static FeatureCSR build_feature_csr(const IntegerVector& xi, const IntegerVector& xp, + const NumericVector& xx, int nr, int nc, + bool transpose) { + FeatureCSR csr; + csr.ncomp = transpose ? nr : nc; + csr.nfeat = transpose ? nc : nr; + csr.totals.resize(csr.ncomp, 0.0); - // Build feature CSR: for each feature f, which items are nonzero - std::vector feat_count(nfeat, 0); + std::vector feat_count(csr.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]; + csr.totals[col] += xx[idx]; feat_count[xi[idx]]++; } } @@ -242,128 +232,249 @@ NumericMatrix c_weighted_jaccard_dense(const S4& x, bool transpose = false, 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]; + csr.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]; + csr.feat_offsets.resize(csr.nfeat + 1, 0); + for (int f = 0; f < csr.nfeat; ++f) + csr.feat_offsets[f + 1] = csr.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()); + csr.total_nnz = csr.feat_offsets[csr.nfeat]; + csr.feat_items.resize(csr.total_nnz); + csr.feat_vals.resize(csr.total_nnz); + std::vector feat_fill(csr.feat_offsets.begin(), csr.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]; + csr.feat_items[dest] = col; + csr.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]; + csr.feat_items[dest] = xi[idx]; + csr.feat_vals[dest] = xx[idx]; } } } - NumericMatrix out(ncomp, ncomp); + return csr; +} + +//' Dense weighted Jaccard similarity via C++ +//' +//' Compute the full weighted Jaccard similarity matrix for a dgCMatrix using +//' an adaptive dense accumulation strategy. For small output matrices, uses a +//' feature-oriented loop; for large outputs, switches to a column-oriented +//' loop for better cache performance. +//' +//' @param x A dgCMatrix (sparse column-compressed matrix) +//' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare +//' rows +//' @param threads Number of threads (default 4). Set to 0 for all cores. +//' @param triangle If \code{TRUE}, return only the lower triangle as a flat +//' numeric vector in \code{\link{dist}} layout. If \code{FALSE} (default), +//' return a full square matrix. +//' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) +//' instead of similarity. Default \code{FALSE}. +//' @return A dense numeric similarity matrix, or a numeric vector in dist +//' layout when \code{triangle = TRUE}. +//' @export +// [[Rcpp::export]] +SEXP c_weighted_jaccard_dense(const S4& x, bool transpose = false, + int threads = 4, bool triangle = false, + bool distance = false) { + IntegerVector dims = x.slot("Dim"); + IntegerVector xi = x.slot("i"); + IntegerVector xp = x.slot("p"); + NumericVector xx = x.slot("x"); - // Threshold: ncomp^2 * 8 bytes > ~12MB (typical L3 share per core) - const bool use_colwise = (static_cast(ncomp) * ncomp * 8 > 12 * 1024 * 1024); + const int nr = dims[0]; + const int nc = dims[1]; - if (use_colwise) { - // Column-oriented: build item->features CSR, iterate output columns + FeatureCSR csr = build_feature_csr(xi, xp, xx, nr, nc, transpose); + const int ncomp = csr.ncomp; + const int nfeat = csr.nfeat; + const size_t nT = (threads > 0) ? static_cast(threads) + : std::thread::hardware_concurrency(); + + if (triangle) { + // --- Triangle path: output is a flat vector in dist layout --- + // dist position for pair (r,c) with r>c (1-based): + // pos = (c-1)*n - c*(c-1)/2 + (r-c) [1-based] + // We use 0-based: for (ca, cb) with ca > cb: + // pos = cb*n - cb*(cb+1)/2 + (ca-cb) - 1 [0-based] + // which simplifies to: cb*(ncomp-1) - cb*(cb-1)/2 + (ca-1) ... let's use a + // direct formula. + + const std::size_t dist_len = static_cast(ncomp) * (ncomp - 1) / 2; + NumericVector out(dist_len, 0.0); + double* out_ptr = out.begin(); + const double* totals_ptr = csr.totals.data(); + + // Build item->features CSR (always needed for colwise triangle) 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]]++; + for (int idx = csr.feat_offsets[f]; idx < csr.feat_offsets[f + 1]; ++idx) + item_count[csr.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_feats(csr.total_nnz); + std::vector item_vals(csr.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]; + for (int idx = csr.feat_offsets[f]; idx < csr.feat_offsets[f + 1]; ++idx) { + const int item = csr.feat_items[idx]; const int dest = item_fill[item]++; item_feats[dest] = f; - item_vals[dest] = feat_vals[idx]; + item_vals[dest] = csr.feat_vals[idx]; } } - const size_t nT = (threads > 0) ? static_cast(threads) - : std::thread::hardware_concurrency(); - double* out_ptr = &out[0]; - const int* feat_items_ptr = feat_items.data(); - const double* feat_vals_ptr = feat_vals.data(); - const int* feat_offsets_ptr = feat_offsets.data(); + const int* feat_items_ptr = csr.feat_items.data(); + const double* feat_vals_ptr = csr.feat_vals.data(); + const int* feat_offsets_ptr = csr.feat_offsets.data(); const int* item_feats_ptr = item_feats.data(); const double* item_vals_ptr = item_vals.data(); const int* item_offsets_ptr = item_offsets.data(); - RcppThread::parallelFor(0, ncomp, [&](int cb) { - double* col = out_ptr + static_cast(cb) * ncomp; + // For column cb, accumulate min-sums for rows ca > cb only, + // writing directly into dist layout. + // dist index for (ca, cb) with ca > cb (0-based items, 0-based output): + // pos = cb * ncomp - cb*(cb+1)/2 + (ca - cb - 1) + RcppThread::parallelFor(0, ncomp - 1, [=](int cb) { + const std::size_t base = static_cast(cb) * ncomp + - static_cast(cb) * (cb + 1) / 2 + - static_cast(cb) - 1; + // Thread-local accumulator for this column's min-sums (only rows > cb) + std::vector col_ms(ncomp, 0.0); + for (int fi = item_offsets_ptr[cb]; fi < item_offsets_ptr[cb + 1]; ++fi) { const int f = item_feats_ptr[fi]; const double vb = item_vals_ptr[fi]; for (int idx = feat_offsets_ptr[f]; idx < feat_offsets_ptr[f + 1]; ++idx) { - col[feat_items_ptr[idx]] += std::min(feat_vals_ptr[idx], vb); + const int ca = feat_items_ptr[idx]; + if (ca > cb) + col_ms[ca] += std::min(feat_vals_ptr[idx], vb); } } + + // Convert to similarity or distance and store + const double tb = totals_ptr[cb]; + for (int ca = cb + 1; ca < ncomp; ++ca) { + const double ms = col_ms[ca]; + const double denom = totals_ptr[ca] + tb - ms; + double sim = (denom > 0.0) ? ms / denom : 0.0; + out_ptr[base + ca] = distance ? (1.0 - sim) : sim; + } }, nT); + + // Return as a dist object + out.attr("class") = "dist"; + out.attr("Size") = ncomp; + out.attr("Diag") = false; + out.attr("Upper") = false; + return out; + } 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; + // --- Full matrix path (existing) --- + 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) { + std::vector item_count(ncomp, 0); + for (int f = 0; f < nfeat; ++f) + for (int idx = csr.feat_offsets[f]; idx < csr.feat_offsets[f + 1]; ++idx) + item_count[csr.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(csr.total_nnz); + std::vector item_vals(csr.total_nnz); + std::vector item_fill(item_offsets.begin(), item_offsets.end()); + + for (int f = 0; f < nfeat; ++f) { + for (int idx = csr.feat_offsets[f]; idx < csr.feat_offsets[f + 1]; ++idx) { + const int item = csr.feat_items[idx]; + const int dest = item_fill[item]++; + item_feats[dest] = f; + item_vals[dest] = csr.feat_vals[idx]; } } - } - } - // Convert min_sums to similarity - { - const size_t nT = (threads > 0) ? static_cast(threads) - : std::thread::hardware_concurrency(); - double* out_ptr = &out[0]; - const double* totals_ptr = totals.data(); - RcppThread::parallelFor(0, ncomp, [=](int cb) { - for (int ca = 0; ca < ncomp; ++ca) { - const std::size_t pos = static_cast(cb) * ncomp + ca; - if (ca == cb) { - out_ptr[pos] = 1.0; - } else { - const double ms = out_ptr[pos]; - const double denom = totals_ptr[ca] + totals_ptr[cb] - ms; - out_ptr[pos] = (denom > 0.0) ? ms / denom : 0.0; + double* out_ptr = &out[0]; + const int* feat_items_ptr = csr.feat_items.data(); + const double* feat_vals_ptr = csr.feat_vals.data(); + const int* feat_offsets_ptr = csr.feat_offsets.data(); + const int* item_feats_ptr = item_feats.data(); + const double* item_vals_ptr = item_vals.data(); + const int* item_offsets_ptr = item_offsets.data(); + + RcppThread::parallelFor(0, ncomp, [&](int cb) { + double* col = out_ptr + static_cast(cb) * ncomp; + for (int fi = item_offsets_ptr[cb]; fi < item_offsets_ptr[cb + 1]; ++fi) { + const int f = item_feats_ptr[fi]; + const double vb = item_vals_ptr[fi]; + for (int idx = feat_offsets_ptr[f]; idx < feat_offsets_ptr[f + 1]; ++idx) { + col[feat_items_ptr[idx]] += std::min(feat_vals_ptr[idx], vb); + } + } + }, nT); + } else { + // Feature-oriented: iterate features, scatter into output + for (int f = 0; f < nfeat; ++f) { + const int start = csr.feat_offsets[f]; + const int end = csr.feat_offsets[f + 1]; + const int k = end - start; + + for (int a = 0; a < k; ++a) { + const int ca = csr.feat_items[start + a]; + const double va = csr.feat_vals[start + a]; + out[static_cast(ca) * ncomp + ca] += va; + for (int b = a + 1; b < k; ++b) { + const int cb = csr.feat_items[start + b]; + const double mv = std::min(va, csr.feat_vals[start + b]); + out[static_cast(ca) * ncomp + cb] += mv; + out[static_cast(cb) * ncomp + ca] += mv; + } } } - }, nT); - } + } - return out; + // Convert min_sums to similarity (or distance) + { + double* out_ptr = &out[0]; + const double* totals_ptr = csr.totals.data(); + RcppThread::parallelFor(0, ncomp, [=](int cb) { + for (int ca = 0; ca < ncomp; ++ca) { + const std::size_t pos = static_cast(cb) * ncomp + ca; + if (ca == cb) { + out_ptr[pos] = distance ? 0.0 : 1.0; + } else { + const double ms = out_ptr[pos]; + const double denom = totals_ptr[ca] + totals_ptr[cb] - ms; + double sim = (denom > 0.0) ? ms / denom : 0.0; + out_ptr[pos] = distance ? (1.0 - sim) : sim; + } + } + }, nT); + } + + return out; + } } diff --git a/tests/testthat/test-weighted-jaccard.R b/tests/testthat/test-weighted-jaccard.R index d595252..06875ca 100644 --- a/tests/testthat/test-weighted-jaccard.R +++ b/tests/testthat/test-weighted-jaccard.R @@ -63,3 +63,73 @@ test_that("c_weighted_jaccard_sparse matches dense", { sim_sparse_t <- c_weighted_jaccard_sparse(m, transpose = TRUE, display_progress = FALSE) expect_equal(as.matrix(sim_sparse_t), sim_dense_t, tolerance = 1e-12) }) + +test_that("c_weighted_jaccard_dense triangle returns dist object", { + m <- Matrix::sparseMatrix( + i = c(1L, 2L, 1L, 2L, 3L, 3L), + j = c(1L, 1L, 2L, 2L, 2L, 3L), + x = c(4, 2, 1, 3, 3, 1), + dims = c(3L, 3L) + ) + + full <- c_weighted_jaccard_dense(m) + + # triangle=TRUE returns a dist object + tri <- c_weighted_jaccard_dense(m, triangle = TRUE) + expect_s3_class(tri, "dist") + expect_equal(attr(tri, "Size"), ncol(m)) + + # triangle similarity matches as.dist of full similarity + # (as.dist extracts lower triangle; dist stores 1-sim by convention, + # but here we just compare the raw values) + ref_lower <- full[lower.tri(full)] + expect_equal(as.numeric(tri), ref_lower, tolerance = 1e-12) + + # triangle + distance matches as.dist(1 - full) + tri_d <- c_weighted_jaccard_dense(m, triangle = TRUE, distance = TRUE) + expect_s3_class(tri_d, "dist") + expect_equal(as.numeric(tri_d), 1 - ref_lower, tolerance = 1e-12) +}) + +test_that("c_weighted_jaccard_dense distance works for full matrix", { + m <- Matrix::sparseMatrix( + i = c(1L, 2L, 1L, 2L, 3L, 3L), + j = c(1L, 1L, 2L, 2L, 2L, 3L), + x = c(4, 2, 1, 3, 3, 1), + dims = c(3L, 3L) + ) + + full <- c_weighted_jaccard_dense(m) + full_d <- c_weighted_jaccard_dense(m, distance = TRUE) + expect_equal(full_d, 1 - full, tolerance = 1e-12) +}) + +test_that("c_weighted_jaccard_sparse triangle and distance args", { + m <- Matrix::sparseMatrix( + i = c(1L, 2L, 1L, 2L, 3L, 3L), + j = c(1L, 1L, 2L, 2L, 2L, 3L), + x = c(4, 2, 1, 3, 3, 1), + dims = c(3L, 3L) + ) + + # Default: dgCMatrix + sp <- c_weighted_jaccard_sparse(m, display_progress = FALSE) + expect_true(inherits(sp, "dgCMatrix")) + + # triangle=TRUE: dsCMatrix + sp_tri <- c_weighted_jaccard_sparse(m, display_progress = FALSE, triangle = TRUE) + expect_true(inherits(sp_tri, "dsCMatrix")) + expect_equal(as.matrix(sp_tri), as.matrix(sp), tolerance = 1e-12) + + # distance=TRUE warns + expect_warning( + sp_d <- c_weighted_jaccard_sparse(m, display_progress = FALSE, distance = TRUE), + "distance=TRUE" + ) + # Stored entries are correct (1-sim); structural zeros remain 0 instead of 1 + # — this is the inherent limitation the warning describes + full <- c_weighted_jaccard_dense(m) + ref_d <- 1 - full + ref_d[full == 0 & row(full) != col(full)] <- 0 # zero out unstored positions + expect_equal(as.matrix(sp_d), ref_d, tolerance = 1e-12) +}) From 055fa9ccf909bf6e021c4ae283cb889ab0b6a030 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Wed, 24 Jun 2026 18:18:30 +0100 Subject: [PATCH 8/8] Cross-link dgCMatrix mentions to Matrix:dgCMatrix-class Use \link[Matrix:dgCMatrix-class]{dgCMatrix} in roxygen for the hand-written R/weighted_jaccard.R and the cpp roxygen blocks in src/weighted_jaccard.cpp, regenerating RcppExports and the .Rd files. WORDLIST picks up the new spelling terms surfaced by spelling::spell_check_package(). --- R/RcppExports.R | 9 +++++---- R/weighted_jaccard.R | 9 +++++---- inst/WORDLIST | 2 ++ man/c_weighted_jaccard_dense.Rd | 9 +++++---- man/c_weighted_jaccard_sparse.Rd | 9 +++++---- src/weighted_jaccard.cpp | 9 +++++---- 6 files changed, 27 insertions(+), 20 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 8ea311b..c113355 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -142,12 +142,13 @@ weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE, display_ #' Dense weighted Jaccard similarity via C++ #' -#' Compute the full weighted Jaccard similarity matrix for a dgCMatrix using +#' Compute the full weighted Jaccard similarity matrix for a +#' \link[Matrix:dgCMatrix-class]{dgCMatrix} using #' an adaptive dense accumulation strategy. For small output matrices, uses a #' feature-oriented loop; for large outputs, switches to a column-oriented #' loop for better cache performance. #' -#' @param x A dgCMatrix (sparse column-compressed matrix) +#' @param x A \link[Matrix:dgCMatrix-class]{dgCMatrix} (sparse column-compressed matrix) #' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare #' rows #' @param threads Number of threads (default 4). Set to 0 for all cores. @@ -156,8 +157,8 @@ weighted_jaccard_sparse_fill <- function(x, pattern, transpose = FALSE, display_ #' return a full square matrix. #' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) #' instead of similarity. Default \code{FALSE}. -#' @return A dense numeric similarity matrix, or a numeric vector in dist -#' layout when \code{triangle = TRUE}. +#' @return A dense numeric similarity matrix, or a numeric vector in +#' \code{\link{dist}} layout when \code{triangle = TRUE}. #' @export c_weighted_jaccard_dense <- function(x, transpose = FALSE, threads = 4L, triangle = FALSE, distance = FALSE) { .Call(`_natcpp_c_weighted_jaccard_dense`, x, transpose, threads, triangle, distance) diff --git a/R/weighted_jaccard.R b/R/weighted_jaccard.R index 6b73b15..870a61e 100644 --- a/R/weighted_jaccard.R +++ b/R/weighted_jaccard.R @@ -1,12 +1,13 @@ #' Sparse weighted Jaccard similarity via C++ #' -#' Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a +#' Compute the weighted Jaccard similarity matrix for a +#' \link[Matrix:dgCMatrix-class]{dgCMatrix}, returning a #' sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute min-sums #' only for column (or row) pairs that share at least one non-zero feature, then #' normalises to similarity. Only the upper triangle is computed, taking #' advantage of the symmetry of the Jaccard index. #' -#' @param x A dgCMatrix (sparse column-compressed matrix) +#' @param x A \link[Matrix:dgCMatrix-class]{dgCMatrix} (sparse column-compressed matrix) #' @param transpose If \code{FALSE} (default), compare columns; if #' \code{TRUE}, compare rows. #' @param display_progress Whether to show a text progress bar (default @@ -15,12 +16,12 @@ #' Set to 0 to use all available cores. #' @param triangle If \code{TRUE}, return a symmetric \code{dsCMatrix} #' (upper triangle only). If \code{FALSE} (default), return a general -#' \code{dgCMatrix}. +#' \link[Matrix:dgCMatrix-class]{dgCMatrix}. #' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) #' instead of similarity. Default \code{FALSE}. A warning is issued since #' sparse distance matrices are typically dense. #' @return A sparse similarity (or distance) matrix: \code{dsCMatrix} when -#' \code{triangle = TRUE}, \code{dgCMatrix} otherwise. +#' \code{triangle = TRUE}, \link[Matrix:dgCMatrix-class]{dgCMatrix} otherwise. #' @importFrom methods as #' @export #' @seealso \code{\link{c_weighted_jaccard_dense}} for the dense equivalent diff --git a/inst/WORDLIST b/inst/WORDLIST index 93c16a1..29595b8 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,4 +1,6 @@ CMD +Codecov +Jaccard Lifecycle NeuroAnatomy Nx diff --git a/man/c_weighted_jaccard_dense.Rd b/man/c_weighted_jaccard_dense.Rd index 1d074e0..e1b1f9c 100644 --- a/man/c_weighted_jaccard_dense.Rd +++ b/man/c_weighted_jaccard_dense.Rd @@ -13,7 +13,7 @@ c_weighted_jaccard_dense( ) } \arguments{ -\item{x}{A dgCMatrix (sparse column-compressed matrix)} +\item{x}{A \link[Matrix:dgCMatrix-class]{dgCMatrix} (sparse column-compressed matrix)} \item{transpose}{If \code{FALSE}, compare columns; if \code{TRUE}, compare rows} @@ -28,11 +28,12 @@ return a full square matrix.} instead of similarity. Default \code{FALSE}.} } \value{ -A dense numeric similarity matrix, or a numeric vector in dist - layout when \code{triangle = TRUE}. +A dense numeric similarity matrix, or a numeric vector in + \code{\link{dist}} layout when \code{triangle = TRUE}. } \description{ -Compute the full weighted Jaccard similarity matrix for a dgCMatrix using +Compute the full weighted Jaccard similarity matrix for a +\link[Matrix:dgCMatrix-class]{dgCMatrix} using an adaptive dense accumulation strategy. For small output matrices, uses a feature-oriented loop; for large outputs, switches to a column-oriented loop for better cache performance. diff --git a/man/c_weighted_jaccard_sparse.Rd b/man/c_weighted_jaccard_sparse.Rd index ffcefe3..cadadda 100644 --- a/man/c_weighted_jaccard_sparse.Rd +++ b/man/c_weighted_jaccard_sparse.Rd @@ -14,7 +14,7 @@ c_weighted_jaccard_sparse( ) } \arguments{ -\item{x}{A dgCMatrix (sparse column-compressed matrix)} +\item{x}{A \link[Matrix:dgCMatrix-class]{dgCMatrix} (sparse column-compressed matrix)} \item{transpose}{If \code{FALSE} (default), compare columns; if \code{TRUE}, compare rows.} @@ -27,7 +27,7 @@ Set to 0 to use all available cores.} \item{triangle}{If \code{TRUE}, return a symmetric \code{dsCMatrix} (upper triangle only). If \code{FALSE} (default), return a general -\code{dgCMatrix}.} +\link[Matrix:dgCMatrix-class]{dgCMatrix}.} \item{distance}{If \code{TRUE}, return distance (\code{1 - similarity}) instead of similarity. Default \code{FALSE}. A warning is issued since @@ -35,10 +35,11 @@ sparse distance matrices are typically dense.} } \value{ A sparse similarity (or distance) matrix: \code{dsCMatrix} when - \code{triangle = TRUE}, \code{dgCMatrix} otherwise. + \code{triangle = TRUE}, \link[Matrix:dgCMatrix-class]{dgCMatrix} otherwise. } \description{ -Compute the weighted Jaccard similarity matrix for a dgCMatrix, returning a +Compute the weighted Jaccard similarity matrix for a +\link[Matrix:dgCMatrix-class]{dgCMatrix}, returning a sparse result. Uses \code{weighted_jaccard_sparse_fill} to compute min-sums only for column (or row) pairs that share at least one non-zero feature, then normalises to similarity. Only the upper triangle is computed, taking diff --git a/src/weighted_jaccard.cpp b/src/weighted_jaccard.cpp index 4ae611e..396b542 100644 --- a/src/weighted_jaccard.cpp +++ b/src/weighted_jaccard.cpp @@ -270,12 +270,13 @@ static FeatureCSR build_feature_csr(const IntegerVector& xi, const IntegerVector //' Dense weighted Jaccard similarity via C++ //' -//' Compute the full weighted Jaccard similarity matrix for a dgCMatrix using +//' Compute the full weighted Jaccard similarity matrix for a +//' \link[Matrix:dgCMatrix-class]{dgCMatrix} using //' an adaptive dense accumulation strategy. For small output matrices, uses a //' feature-oriented loop; for large outputs, switches to a column-oriented //' loop for better cache performance. //' -//' @param x A dgCMatrix (sparse column-compressed matrix) +//' @param x A \link[Matrix:dgCMatrix-class]{dgCMatrix} (sparse column-compressed matrix) //' @param transpose If \code{FALSE}, compare columns; if \code{TRUE}, compare //' rows //' @param threads Number of threads (default 4). Set to 0 for all cores. @@ -284,8 +285,8 @@ static FeatureCSR build_feature_csr(const IntegerVector& xi, const IntegerVector //' return a full square matrix. //' @param distance If \code{TRUE}, return distance (\code{1 - similarity}) //' instead of similarity. Default \code{FALSE}. -//' @return A dense numeric similarity matrix, or a numeric vector in dist -//' layout when \code{triangle = TRUE}. +//' @return A dense numeric similarity matrix, or a numeric vector in +//' \code{\link{dist}} layout when \code{triangle = TRUE}. //' @export // [[Rcpp::export]] SEXP c_weighted_jaccard_dense(const S4& x, bool transpose = false,