weighted jaccard functions

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",
@@ -18,14 +18,15 @@ 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: 
+Imports:
+    Matrix,
     Rcpp (>= 1.0.6)
-Suggests: 
+Suggests:
     spelling,
     testthat (>= 3.0.0)
 LinkingTo: 
     Rcpp
 Config/testthat/edition: 3
 Encoding: UTF-8
 Language: en-GB
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3

NAMESPACE

@@ -10,5 +10,7 @@ export(c_sub2ind)
 export(c_topntail)
 export(c_topntail_list)
 export(c_total_cable)
+export(c_weighted_jaccard_dense)
+export(c_weighted_jaccard_sparse)
 import(Rcpp)
 useDynLib(natcpp, .registration = TRUE)

R/RcppExports.R

@@ -135,3 +135,25 @@ c_EdgeListFromSegList <- function(L) {
     .Call(`_natcpp_c_EdgeListFromSegList`, L)
 }
 
+#' @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++
+#'
+#' 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
+c_weighted_jaccard_dense <- function(x, transpose = FALSE) {
+    .Call(`_natcpp_c_weighted_jaccard_dense`, x, transpose)
+}
+

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
+}

src/RcppExports.cpp

@@ -135,6 +135,31 @@ BEGIN_RCPP
     return rcpp_result_gen;
 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) {
+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));
+    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) {
+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));
+    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", (DL_FUNC) &_natcpp_weighted_jaccard_sparse_fill, 3},
+    {"_natcpp_c_weighted_jaccard_dense", (DL_FUNC) &_natcpp_c_weighted_jaccard_dense, 2},
     {NULL, NULL, 0}
 };