Reorganized the code and updated the docs

.vscode/.env

@@ -13,7 +13,7 @@ R_PDFVIEWER=/usr/bin/xdg-open
 R_PLATFORM=x86_64-redhat-linux-gnu
 R_PRINTCMD=
 R_RD4PDF=times,inconsolata,hyper
-R_SESSION_TMPDIR=/tmp/RtmpMNylqA
+R_SESSION_TMPDIR=/tmp/RtmpHERMsP
 R_SHARE_DIR=/usr/share/R
 R_STRIP_SHARED_LIB=strip --strip-unneeded
 R_STRIP_STATIC_LIB=strip --strip-debug

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: diffusr
 Type: Package
 Title: Network Diffusion Algorithms
-Version: 0.2.0
+Version: 0.2.1
 Date: 2018-04-20
 Authors@R: person("Simon", "Dirmeier",
                   email = "[email protected]",

NEWS.md

@@ -1,3 +1,8 @@
+# diffusr v0.2.1
+
+* Allow sparse matrix as input
+* Lint the code
+
 # diffusr v0.2.0
 
 * Reorganized the code

R/heat_diffusion.R

@@ -75,7 +75,6 @@ heat.diffusion <- function(h0, graph, t = 0.5, ...) {
   diag(graph) <- 0
   if (is.dgCMatrix(graph)) {
     assert_dgCMatrix(graph)
-    # TODO: sparse matrix
   } else {
     assert(
       test_matrix(graph, mode = "numeric", nrows = n_elements,
@@ -84,7 +83,7 @@ heat.diffusion <- function(h0, graph, t = 0.5, ...) {
       any(graph >= 0),
       combine = "and"
     )
-    heat <- heat_diffusion_(h0, laplacian_(graph), t)
   }
+  heat <- heat_diffusion_(h0, laplacian_(graph), t)
   return(heat)
 }

R/mat_util.R

@@ -22,9 +22,12 @@
 #'
 #' @export
 #'
-#' @param obj  matrix/vector that is stochstically normalized
+#' @param obj  \code{\link[base]{matrix}} (or
+#'   \code{\link[Matrix:dgCMatrix-class]{dgCMatrix}}, \link[base]{vector}) that
+#'   is stochstically normalized
 #' @param ...  additional params
-#' @return  returns the normalized matrix/vector
+#' @return  returns the normalized \code{\link[base]{matrix}} (or
+#'   \code{\link[Matrix:dgCMatrix-class]{dgCMatrix}}, \link[base]{vector})
 #'
 #' @importFrom checkmate assert check_matrix test_numeric test_atomic_vector
 #'
@@ -74,7 +77,9 @@ normalize.stochastic <- function(obj, ...) {
 #'
 #' @export
 #'
-#' @param obj  matrix for which the Laplacian is calculated
+#' @param obj  \code{\link[base]{matrix}} (or
+#'   \code{\link[Matrix:dgCMatrix-class]{dgCMatrix}}) for which the Laplacian is
+#'   calculated
 #' @param ...  additional params
 #' @return  returns the Laplacian
 #'
@@ -104,9 +109,12 @@ normalize.laplacian <- function(obj, ...) {
 #'
 #' @export
 #'
-#' @param obj  matrix for which hubs are corrected
+#' @param obj  \code{\link[base]{matrix}} (or
+#'   \code{\link[Matrix:dgCMatrix-class]{dgCMatrix}}) for which hubs
+#'   are corrected
 #'
-#' @return  returns the matrix with hub correction
+#' @return  returns the \code{\link[base]{matrix}} (or
+#'   \code{\link[Matrix:dgCMatrix-class]{dgCMatrix}}) with hub correction
 #'
 #' @useDynLib diffusr
 #'

R/mrw.R

@@ -31,7 +31,8 @@
 #'  (does not need to sum to one).
 #'
 #' @param graph  an (\code{n x n})-dimensional numeric non-negative adjacence
-#' matrix representing the graph
+#'   \code{\link[base]{matrix}} (or
+#'   \code{\link[Matrix:dgCMatrix-class]{dgCMatrix}}) representing the graph
 #'
 #' @param r  a scalar between (0, 1). restart probability if a Markov random
 #' walk with restart is desired

R/neighbors.R

@@ -32,7 +32,9 @@
 #'   for which the algorithm is applied iteratively
 #'
 #' @param graph an (\eqn{n \times n})-dimensional numeric non-negative adjacence
-#'   matrix representing the graph
+#'   \code{\link[base]{matrix}} (or
+#'   \code{\link[Matrix:dgCMatrix-class]{dgCMatrix}}, \link[base]{vector})
+#'   representing the graph
 #'
 #' @param k the depth of the nearest neighbor search, e.g. the depth of the
 #'   graph traversal

inst/include/diffusr.h

@@ -1,9 +1,3 @@
-#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64)
-#define WINDOWS 1
-#else
-#define WINDOWS 0
-#endif
-
 #ifdef INTEL_MKL_VERSION 
 #define EIGEN_USE_MKL_ALL
 #endif

src/heat_diffusion.cpp

@@ -34,17 +34,14 @@
 // [[Rcpp::export]]
 MatrixXd heat_diffusion_(const MatrixXd& v0,
                          const MatrixXd& W,
-                         const double    t)
-{
+                         const double    t) {
     SelfAdjointEigenSolver<MatrixXd> es(W);
     MatrixXd                         V  = es.eigenvectors();
     VectorXd                         D  = es.eigenvalues();
     MatrixXd                         co = V.transpose() * v0;
     // solution to the heat equation at time t
-    for (int i = 0; i < co.rows(); ++i)
-    {
-        for (int j = 0; j < co.cols(); ++j)
-        {
+    for (int i = 0; i < co.rows(); ++i) {
+        for (int j = 0; j < co.cols(); ++j) {
             if (j % 25 == 0)
                 Rcpp::checkUserInterrupt();
             co(i, j) *= std::exp(-D(i) * t);
@@ -53,4 +50,4 @@ MatrixXd heat_diffusion_(const MatrixXd& v0,
     co = V * co;
 
     return co;
-}
+}

src/mat_util.cpp

@@ -29,8 +29,7 @@
 //' @return  returns the normalized matrix
 // [[Rcpp::interfaces(r, cpp)]]
 // [[Rcpp::export]]
-MatrixXd stoch_col_norm_(const MatrixXd& W)
-{
+MatrixXd stoch_col_norm_(const MatrixXd& W) {
     MatrixXd res(W.rows(), W.cols());
     VectorXd colsums      = W.colwise().sum();
     const double    empt_col_val = 1.0 / W.cols();
@@ -51,6 +50,8 @@ template <typename T> MatrixXd compute_laplacian(const T &W, const ArrayXd rowsu
 
 template <> MatrixXd compute_laplacian(const MatrixXd &W, const ArrayXd rowsums, const size_t P) {
     ArrayXd rowsums_m = rowsums.replicate(1, P).array().sqrt().inverse();
+    // remove 1/0 = Inf
+    rowsums_m = (rowsums_m.isInf()).select(0, rowsums_m);
     MatrixXd res = - W.array() * (rowsums_m.transpose() * rowsums_m) * W.array().cast<bool>().cast<double>();
     return res;
 }
@@ -61,23 +62,18 @@ template <> MatrixXd compute_laplacian(const SpMat &W, const ArrayXd rowsums, co
     return res;
 }
 
-template <typename T> MatrixXd laplacian_t(const T& W)
-{
+template <typename T> MatrixXd laplacian_t(const T& W) {
     const size_t       P = W.rows();
     //MatrixXd res(P, P);
     // Equivalent with: VectorXd rowsums = W.rowwise().sum();
     ArrayXd rowsums = (W * VectorXd::Ones(P)).array();
 
     // for diagnoals
-    ArrayXd W_diags = 1 - W.diagonal().array() / rowsums;
+    ArrayXd W_diags = 1 - W.diagonal().array() / (rowsums == 0).select(1, rowsums);
     W_diags *= (rowsums != 0.0).cast<double>();
 
     // for others
-    //ArrayXd rowsums_m = rowsums.replicate(1, P).array().sqrt().inverse();
-    //MatrixXd res = - W.array() * (rowsums_m.transpose() * rowsums_m) * W.array().cast<bool>().cast<double>();
     MatrixXd res = compute_laplacian(W, rowsums, P);
-    //res.noalias() = rowsums_m.transpose().cwiseProduct(rowsums_m).cwiseProduct(-W).cwiseProduct(W.cast<bool>().cast<double>());
-    // Rcout << - W * rowsums_m.transpose().cwiseProduct(rowsums_m) << endl << endl;
     res.diagonal().array() = W_diags;
     return res;
 }

src/mrw.cpp

@@ -61,8 +61,7 @@ VectorXd mrwr_t(const MatrixXd& p0,
 
         int iter = 0;
         pt = p0;
-        do
-        {
+        do {
             if (iter % 25 == 0)
                 Rcpp::checkUserInterrupt();
             pold = pt;

src/neighbors.cpp

@@ -22,25 +22,21 @@
 
 #include "../inst/include/diffusr.h"
 
-struct distance_comparator
-{
-    bool operator()(pair<int, double>& lhs, pair<int, double>& rhs)
-    {
+struct distance_comparator {
+    bool operator()(pair<int, double>& lhs, pair<int, double>& rhs) {
         return lhs.second < rhs.second;
     }
 };
 
-bool equals(const double val, const double cmp, const double delta)
-{
+bool equals(const double val, const double cmp, const double delta) {
     return val <= cmp + delta && val >= cmp - delta;
 }
 
 vector<pair<int, double>> current_neighbors(
   priority_queue<pair<int, double>,
                  vector<pair<int, double>>,
                  distance_comparator>& queue,
-  vector<uint8_t>& visited)
-{
+  vector<uint8_t>& visited) {
 
     // get the nearest neighbor
     pair<int, double> cn = queue.top();
@@ -49,14 +45,12 @@ vector<pair<int, double>> current_neighbors(
 
     vector<pair<int, double>> curr_nei;
     // add the node itself to the list, since we iterate over the list later
-    if (!visited[cn.first])
-    {
+    if (!visited[cn.first]) {
         curr_nei.push_back(cn);
     }
 
     // add neighbors that are as close as the first neighbor `cn`
-    while (queue.size() && equals(queue.top().second, cn.second, .001))
-    {
+    while (queue.size() && equals(queue.top().second, cn.second, .001)) {
         pair<int, double> nn = queue.top();
         if (!visited[nn.first])
             curr_nei.push_back(nn);
@@ -71,8 +65,7 @@ void add_neighbor_to_queue(
                  vector<pair<int, double>>,
                  distance_comparator>& queue,
   const NumericMatrix& W,
-  const pair<int, double>& cn)
-{
+  const pair<int, double>& cn) {
     for (int i = 0; i < W.cols(); ++i)
     {
         if (i != cn.first && W(cn.first, i) > 0)
@@ -145,8 +138,7 @@ List neighbors_t(const vector<int> &node_idxs,
     // parallelize node search
     // Fucking R bug: you can never have R objects within the parallel sections
     // #pragma omp parallel for
-    for (uint32_t i = 0; i < len; ++i)
-    {
+    for (uint32_t i = 0; i < len; ++i) {
         // substract one, cause R was one-based
         int node_idx = node_idxs[i] - 1;
         // neighbors of current node