Sparse matrix

Author: randef1ned

DESCRIPTION

@@ -25,7 +25,8 @@ Imports:
     igraph,
     methods,
     checkmate,
-    matrixStats
+    matrixStats,
+    sparseMatrixStats
 Suggests:
     knitr,
     rmarkdown,

NAMESPACE

@@ -14,6 +14,7 @@ importFrom(checkmate,assert_int)
 importFrom(checkmate,assert_integer)
 importFrom(checkmate,assert_logical)
 importFrom(checkmate,assert_number)
+importFrom(checkmate,assert_vector)
 importFrom(checkmate,check_matrix)
 importFrom(checkmate,check_numeric)
 importFrom(checkmate,test_atomic_vector)
@@ -23,4 +24,5 @@ importFrom(igraph,components)
 importFrom(igraph,graph_from_adjacency_matrix)
 importFrom(matrixStats,colSums2)
 importFrom(methods,is)
+importFrom(sparseMatrixStats,colAnyNAs)
 useDynLib(diffusr)

NEWS.md

@@ -42,6 +42,6 @@
 	* Codecov
 * Initial submission to CRAN
 
-# Author
+## Author
 
 * Simon Dirmeier <a href="mailto:simon.dirmeier@gmx.de">simon.dirmeier@gmx.de</a>

R/RcppExports.R

@@ -13,22 +13,42 @@ laplacian_ <- function(W) {
     .Call('_diffusr_laplacian_', PACKAGE = 'diffusr', W)
 }
 
+laplacian_s <- function(W) {
+    .Call('_diffusr_laplacian_s', PACKAGE = 'diffusr', W)
+}
+
 node_degrees_ <- function(W) {
     .Call('_diffusr_node_degrees_', PACKAGE = 'diffusr', W)
 }
 
+node_degrees_s <- function(W) {
+    .Call('_diffusr_node_degrees_s', PACKAGE = 'diffusr', W)
+}
+
 hub_normalize_ <- function(W) {
     .Call('_diffusr_hub_normalize_', PACKAGE = 'diffusr', W)
 }
 
+hub_normalize_s <- function(W) {
+    .Call('_diffusr_hub_normalize_s', PACKAGE = 'diffusr', W)
+}
+
 mrwr_ <- function(p0, W, r, thresh, niter, do_analytical) {
     .Call('_diffusr_mrwr_', PACKAGE = 'diffusr', p0, W, r, thresh, niter, do_analytical)
 }
 
+mrwr_s <- function(p0, W, r, thresh, niter, do_analytical) {
+    .Call('_diffusr_mrwr_s', PACKAGE = 'diffusr', p0, W, r, thresh, niter, do_analytical)
+}
+
 neighbors_ <- function(node_idxs, W, k) {
     .Call('_diffusr_neighbors_', PACKAGE = 'diffusr', node_idxs, W, k)
 }
 
+neighbors_s <- function(node_idxs, W, k) {
+    .Call('_diffusr_neighbors_s', PACKAGE = 'diffusr', node_idxs, W, k)
+}
+
 # Register entry points for exported C++ functions
 methods::setLoadAction(function(ns) {
     .Call('_diffusr_RcppExport_registerCCallable', PACKAGE = 'diffusr')

R/diffusr-package.R

@@ -54,4 +54,10 @@
 #' @useDynLib diffusr
 #'
 #' @importFrom Rcpp sourceCpp
+#' @importFrom checkmate assert_vector
+#' @importFrom sparseMatrixStats colAnyNAs
 NULL
+
+.onAttach <- function(libname, pkgname) {
+  assert_vector(colAnyNAs(matrix(0)))
+}

R/mat_util.R

@@ -27,11 +27,10 @@
 #' @return  returns the normalized matrix/vector
 #'
 #' @importFrom checkmate assert check_matrix test_numeric test_atomic_vector
-#' @importFrom matrixStats colSums2
 #'
 #' @examples
-#'  W <- matrix(abs(rnorm(10000)), 100, 100)
-#'  stoch.W <- normalize.stochastic(W)
+#' W <- matrix(abs(rnorm(10000)), 100, 100)
+#' stoch.W <- normalize.stochastic(W)
 normalize.stochastic <- function(obj, ...) {
   is_matrix <- FALSE
   if (test_numeric(obj, lower = 0, finite = TRUE, any.missing = FALSE,
@@ -53,13 +52,17 @@ normalize.stochastic <- function(obj, ...) {
     is_matrix <- TRUE
   }
   if (is_matrix) {
-    sums <- colSums2(obj)
+    sums <- colSums3(obj, !is_matrix)
     if (!all(.equals.double(sums, 1, .001))) {
       message("normalizing column vectors!")
       empt_col_val <- 1.0 / ncol(obj)
 
       obj <- obj / sums[col(obj)]
-      obj[, which(sums < empt_col_val)] <- 0.00001
+      # check if need wipe zeros
+      zeros <- which(sums < empt_col_val)
+      if (length(zeros)) {
+        obj[, zeros] <- 0.00001
+      }
     }
   } else {
     obj <- obj / sum(obj)
@@ -79,12 +82,12 @@ normalize.stochastic <- function(obj, ...) {
 #' @importFrom Rcpp sourceCpp
 #'
 #' @examples
-#'  W <- matrix(abs(rnorm(10000)), 100, 100)
-#'  lapl.W <- normalize.laplacian(W)
+#' W <- matrix(abs(rnorm(10000)), 100, 100)
+#' lapl.W <- normalize.laplacian(W)
 normalize.laplacian <- function(obj, ...) {
   if (is.dgCMatrix(obj)) {
     assert_dgCMatrix(obj)
-    # TODO: sparse matrix
+    return(laplacian_s(obj))
   } else {
     assert(
       check_matrix(obj, mode = "numeric", nrows = ncol(obj), ncols = nrow(obj),
@@ -117,7 +120,7 @@ hub.correction <- function(obj) {
   n_elements <- nrow(obj)
   if (is.dgCMatrix(obj)) {
     assert_dgCMatrix(obj)
-    # TODO: sparse matrix
+    return(hub_normalize_s(obj))
   } else {
     assert(
       check_matrix(obj, mode = "numeric", min.rows = 3, nrows = n_elements,
@@ -137,3 +140,16 @@ hub.correction <- function(obj) {
   comps <- components(adj)
   ifelse(length(comps$csize) == 1, TRUE, FALSE)
 }
+
+#' @noRd
+#' @importFrom matrixStats colSums2
+colSums3 <- function(mat, is.sparse = NULL) {
+  if (is.null(is.sparse)) {
+    is.sparse <- is.dgCMatrix(mat)
+  }
+  if (is.sparse) {
+    sums <- sparseMatrixStats::colSums2(mat)
+  } else {
+    sums <- colSums2(mat)
+  }
+}

R/mrw.R

@@ -55,7 +55,7 @@
 #'  \deqn{ P(j | i) =  1 /degree(i) *  min(1, degree(j)/degree(j))}
 #'  \emph{Note that this will not consider edge weights.}
 #'
-#' @param ergodic.tolerance Tolerance ergodic of the graph.
+#' @param allow.ergodic Allow multiple components in a graph.
 #'
 #' @param return.pt.only Return pt only.
 #'
@@ -94,7 +94,7 @@
 #'
 random.walk <- function(p0, graph, r = 0.5, niter = 1e4, thresh = 1e-4,
                         do.analytical = FALSE, correct.for.hubs = FALSE,
-                        ergodic.tolerance = FALSE, return.pt.only = FALSE) {
+                        allow.ergodic = FALSE, return.pt.only = FALSE) {
   ## Check the fucking inputs
   assert_number(r, lower = 0, upper = 1, na.ok = FALSE, finite = TRUE,
                 null.ok = FALSE)
@@ -111,7 +111,6 @@ random.walk <- function(p0, graph, r = 0.5, niter = 1e4, thresh = 1e-4,
   if (is.dgCMatrix(graph)) {
     assert_dgCMatrix(graph)
     sparse <- TRUE
-    # TODO: sparse matrix
   } else {
     assert(
       test_matrix(graph, mode = "numeric", min.rows = 3, nrows = n_elements,
@@ -142,13 +141,20 @@ random.walk <- function(p0, graph, r = 0.5, niter = 1e4, thresh = 1e-4,
     graph <- hub.correction(graph)
   }
   stoch.graph <- normalize.stochastic(graph)
-  if ((!ergodic.tolerance) && (!.is.ergodic(stoch.graph))) {
+  if ((!allow.ergodic) && (!.is.ergodic(stoch.graph))) {
     stop(paste("the provided graph has more than one component.",
                "It is likely not ergodic."))
   }
 
-  l <- mrwr_(normalize.stochastic(p0),
-             stoch.graph, r, thresh, niter, do.analytical)
+  if (sparse) {
+    # sparse matrix
+    l <- mrwr_s(normalize.stochastic(p0),
+                stoch.graph, r, thresh, niter, do.analytical)
+  } else {
+    # dense matrix
+    l <- mrwr_(normalize.stochastic(p0),
+               stoch.graph, r, thresh, niter, do.analytical)
+  }
   if (!return.pt.only) {
     l <- list(p.inf = l, transition.matrix = stoch.graph)
   }

R/neighbors.R

@@ -68,7 +68,7 @@ nearest.neighbors <- function(nodes, graph, k = 1L, ...) {
   diag(graph) <- 0
   if (is.dgCMatrix(graph)) {
     assert_dgCMatrix(graph)
-    # TODO: sparse matrix
+    neighbors <- neighbors_s(nodes, graph, k)
   } else {
     assert(
       test_matrix(graph, mode = "numeric", min.rows = 3, nrows = n_elements,

inst/include/diffusr.h

@@ -46,4 +46,5 @@ using Eigen::SparseMatrix;
 typedef Eigen::MappedSparseMatrix<double> MSpMat;
 typedef Eigen::Map<MatrixXd> MMatrixXd;
 typedef Eigen::SparseMatrix<double> SpMat;
+typedef Eigen::SparseMatrix<int> SpMati;
 typedef Eigen::Map<VectorXd> MVectorXd;

inst/include/diffusr_RcppExports.h

@@ -88,11 +88,32 @@ namespace diffusr {
         return Rcpp::as<MatrixXd >(rcpp_result_gen);
     }
 
-    inline vector<double> node_degrees_(const MatrixXd& W) {
+    inline MatrixXd laplacian_s(const SpMat& W) {
+        typedef SEXP(*Ptr_laplacian_s)(SEXP);
+        static Ptr_laplacian_s p_laplacian_s = NULL;
+        if (p_laplacian_s == NULL) {
+            validateSignature("MatrixXd(*laplacian_s)(const SpMat&)");
+            p_laplacian_s = (Ptr_laplacian_s)R_GetCCallable("diffusr", "_diffusr_laplacian_s");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_laplacian_s(Shield<SEXP>(Rcpp::wrap(W)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<MatrixXd >(rcpp_result_gen);
+    }
+
+    inline VectorXd node_degrees_(const MatrixXd& W) {
         typedef SEXP(*Ptr_node_degrees_)(SEXP);
         static Ptr_node_degrees_ p_node_degrees_ = NULL;
         if (p_node_degrees_ == NULL) {
-            validateSignature("vector<double>(*node_degrees_)(const MatrixXd&)");
+            validateSignature("VectorXd(*node_degrees_)(const MatrixXd&)");
             p_node_degrees_ = (Ptr_node_degrees_)R_GetCCallable("diffusr", "_diffusr_node_degrees_");
         }
         RObject rcpp_result_gen;
@@ -106,7 +127,28 @@ namespace diffusr {
             throw Rcpp::LongjumpException(rcpp_result_gen);
         if (rcpp_result_gen.inherits("try-error"))
             throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
-        return Rcpp::as<vector<double> >(rcpp_result_gen);
+        return Rcpp::as<VectorXd >(rcpp_result_gen);
+    }
+
+    inline VectorXd node_degrees_s(const MSpMat& W) {
+        typedef SEXP(*Ptr_node_degrees_s)(SEXP);
+        static Ptr_node_degrees_s p_node_degrees_s = NULL;
+        if (p_node_degrees_s == NULL) {
+            validateSignature("VectorXd(*node_degrees_s)(const MSpMat&)");
+            p_node_degrees_s = (Ptr_node_degrees_s)R_GetCCallable("diffusr", "_diffusr_node_degrees_s");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_node_degrees_s(Shield<SEXP>(Rcpp::wrap(W)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<VectorXd >(rcpp_result_gen);
     }
 
     inline MatrixXd hub_normalize_(const MatrixXd& W) {
@@ -130,6 +172,27 @@ namespace diffusr {
         return Rcpp::as<MatrixXd >(rcpp_result_gen);
     }
 
+    inline MatrixXd hub_normalize_s(const SpMat& W) {
+        typedef SEXP(*Ptr_hub_normalize_s)(SEXP);
+        static Ptr_hub_normalize_s p_hub_normalize_s = NULL;
+        if (p_hub_normalize_s == NULL) {
+            validateSignature("MatrixXd(*hub_normalize_s)(const SpMat&)");
+            p_hub_normalize_s = (Ptr_hub_normalize_s)R_GetCCallable("diffusr", "_diffusr_hub_normalize_s");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_hub_normalize_s(Shield<SEXP>(Rcpp::wrap(W)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<MatrixXd >(rcpp_result_gen);
+    }
+
     inline VectorXd mrwr_(const MatrixXd& p0, const MatrixXd& W, const double r, const double thresh, const int niter, const bool do_analytical) {
         typedef SEXP(*Ptr_mrwr_)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_mrwr_ p_mrwr_ = NULL;
@@ -151,11 +214,32 @@ namespace diffusr {
         return Rcpp::as<VectorXd >(rcpp_result_gen);
     }
 
-    inline List neighbors_(const vector<int>& node_idxs, const NumericMatrix& W, const int k) {
+    inline VectorXd mrwr_s(const MatrixXd& p0, const SpMat& W, const double r, const double thresh, const int niter, const bool do_analytical) {
+        typedef SEXP(*Ptr_mrwr_s)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
+        static Ptr_mrwr_s p_mrwr_s = NULL;
+        if (p_mrwr_s == NULL) {
+            validateSignature("VectorXd(*mrwr_s)(const MatrixXd&,const SpMat&,const double,const double,const int,const bool)");
+            p_mrwr_s = (Ptr_mrwr_s)R_GetCCallable("diffusr", "_diffusr_mrwr_s");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_mrwr_s(Shield<SEXP>(Rcpp::wrap(p0)), Shield<SEXP>(Rcpp::wrap(W)), Shield<SEXP>(Rcpp::wrap(r)), Shield<SEXP>(Rcpp::wrap(thresh)), Shield<SEXP>(Rcpp::wrap(niter)), Shield<SEXP>(Rcpp::wrap(do_analytical)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<VectorXd >(rcpp_result_gen);
+    }
+
+    inline List neighbors_(const vector<int>& node_idxs, const MatrixXd& W, const int& k) {
         typedef SEXP(*Ptr_neighbors_)(SEXP,SEXP,SEXP);
         static Ptr_neighbors_ p_neighbors_ = NULL;
         if (p_neighbors_ == NULL) {
-            validateSignature("List(*neighbors_)(const vector<int>&,const NumericMatrix&,const int)");
+            validateSignature("List(*neighbors_)(const vector<int>&,const MatrixXd&,const int&)");
             p_neighbors_ = (Ptr_neighbors_)R_GetCCallable("diffusr", "_diffusr_neighbors_");
         }
         RObject rcpp_result_gen;
@@ -172,6 +256,27 @@ namespace diffusr {
         return Rcpp::as<List >(rcpp_result_gen);
     }
 
+    inline List neighbors_s(const vector<int>& node_idxs, const MSpMat& W, const int& k) {
+        typedef SEXP(*Ptr_neighbors_s)(SEXP,SEXP,SEXP);
+        static Ptr_neighbors_s p_neighbors_s = NULL;
+        if (p_neighbors_s == NULL) {
+            validateSignature("List(*neighbors_s)(const vector<int>&,const MSpMat&,const int&)");
+            p_neighbors_s = (Ptr_neighbors_s)R_GetCCallable("diffusr", "_diffusr_neighbors_s");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_neighbors_s(Shield<SEXP>(Rcpp::wrap(node_idxs)), Shield<SEXP>(Rcpp::wrap(W)), Shield<SEXP>(Rcpp::wrap(k)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<List >(rcpp_result_gen);
+    }
+
 }
 
 #endif // RCPP_diffusr_RCPPEXPORTS_H_GEN_