Adds hub correction

Author: dirmeier

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: diffusr
 Type: Package
 Title: Network Diffusion Algorithms
-Version: 0.1.2
+Version: 0.1.3
 Date: 2017-10-29
 Authors@R: person("Simon", "Dirmeier",
                   email = "simon.dirmeier@gmx.de",

NAMESPACE

@@ -1,8 +1,10 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(hub.correction,numeric)
 S3method(normalize.laplacian,numeric)
 S3method(normalize.stochastic,numeric)
 export(heat.diffusion)
+export(hub.correction)
 export(nearest.neighbors)
 export(normalize.laplacian)
 export(normalize.stochastic)

R/RcppExports.R

@@ -13,6 +13,14 @@ laplacian_ <- function(W) {
     .Call('_diffusr_laplacian_', PACKAGE = 'diffusr', W)
 }
 
+node_degrees_ <- function(W) {
+    .Call('_diffusr_node_degrees_', PACKAGE = 'diffusr', W)
+}
+
+hub_normalize_ <- function(W) {
+    .Call('_diffusr_hub_normalize_', PACKAGE = 'diffusr', W)
+}
+
 mrwr_ <- function(p0, W, r, thresh, niter, do_analytical) {
     .Call('_diffusr_mrwr_', PACKAGE = 'diffusr', p0, W, r, thresh, niter, do_analytical)
 }

R/mat_util.R

@@ -26,13 +26,14 @@
 #' @param ...  additional params
 #' @return  returns the normalized matrix/vector
 #' @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, ...)
 {
   UseMethod("normalize.stochastic")
 }
 
+
 #' @export
 #' @method normalize.stochastic numeric
 normalize.stochastic.numeric <- function(obj, ...)
@@ -58,6 +59,7 @@ normalize.stochastic.numeric <- function(obj, ...)
   return(obj)
 }
 
+
 #' Calculate the Laplacian of a matrix
 #'
 #' @export
@@ -66,13 +68,14 @@ normalize.stochastic.numeric <- function(obj, ...)
 #' @param ...  additional params
 #' @return  returns the Laplacian
 #' @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, ...)
 {
   UseMethod("normalize.laplacian")
 }
 
+
 #' @export
 #' @method normalize.laplacian numeric
 normalize.laplacian.numeric <- function(obj, ...)
@@ -82,16 +85,48 @@ normalize.laplacian.numeric <- function(obj, ...)
   if (any(obj < 0.0))
     stop('please provide a matrix with only non-negative alues!')
   lapl <- laplacian_(obj)
-  return(lapl)
+
+  lapl
+}
+
+
+
+#' Correct for hubs in an adjacency matrix
+#'
+#' @export
+#'
+#' @param obj  matrix for which hubs are corrected
+#' @return  returns the matrix with hub correction
+#' @examples
+#'  W <- matrix(abs(rnorm(10000)), 100, 100)
+#'  cor.hub <- hub.correction(W)
+hub.correction <- function(obj)
+{
+  UseMethod("hub.correction")
+}
+
+
+#' @export
+#' @method hub.correction numeric
+hub.correction.numeric <- function(obj)
+{
+  if (!is.matrix(obj)) stop('please provide a matrix object!')
+  if (nrow(obj) != ncol(obj)) stop('please provide a square matrix!')
+  if (any(obj < 0.0))
+    stop('please provide a matrix with only non-negative alues!')
+  message("Correcting for hub degrees.")
+  hub.mat <- hub_normalize_(obj)
+
+  hub.mat
 }
 
+
 #' @noRd
 #' @importFrom igraph graph_from_adjacency_matrix components
 .is.ergodic <- function(obj)
 {
-  adj   <- igraph::graph_from_adjacency_matrix(obj,
-                                               mode="directed",
-                                               weighted=TRUE)
+  adj   <- igraph::graph_from_adjacency_matrix(
+    obj, mode="directed", weighted=TRUE)
   comps <- igraph::components(adj)
   ifelse(length(comps$csize) == 1, TRUE, FALSE)
 }

R/mrw.R

@@ -46,8 +46,14 @@
 #' stop criterion.
 #' @param do.analytical  boolean if the stationary distribution shall be
 #'  computed solving the analytical solution or rather iteratively
-#' @param ...  additional parameters
-#' @return  returns the stationary distribution as numeric vector
+#' @param correct.for.hubs if \code{TRUE} multiplies a correction factor to the
+#'  nodes, such that the random walk gets not biased to nodes with high
+#'  degree.
+#' @return  returns a list with the following elements
+#'  \itemize{
+#'   \item p.inf  the stationary distribution as numeric vector
+#'   \item transition.matrix the column normalized transition matrix used for the random walk
+#'  }
 #'
 #' @references
 #' Tong, H., Faloutsos, C., & Pan, J. Y. (2006),
@@ -57,32 +63,38 @@
 #' \emph{The American Journal of Human Genetics}\cr \cr
 #'
 #' @examples
-#' # count of nodes
-#' n <- 5
-#' # starting distribution (has to sum to one)
-#' p0    <- as.vector(rmultinom(1, 1, prob=rep(.2, n)))
-#' # adjacency matrix (either normalized or not)
-#' graph <- matrix(abs(rnorm(n*n)), n, n)
-#' # computation of stationary distribution
-#' pt    <- random.walk(p0, graph)
+#'  # count of nodes
+#'  n <- 5
+#'  # starting distribution (has to sum to one)
+#'  p0    <- as.vector(rmultinom(1, 1, prob=rep(.2, n)))
+#'  # adjacency matrix (either normalized or not)
+#'  graph <- matrix(abs(rnorm(n*n)), n, n)
+#'  # computation of stationary distribution
+#'  pt    <- random.walk(p0, graph)
+#'
 setGeneric(
   "random.walk",
-  function(p0, graph, r=.5, niter=1e4, thresh=1e-4, do.analytical=FALSE, ...)
+  function(p0, graph, r=.5,
+           niter=1e4, thresh=1e-4, do.analytical=FALSE,
+           correct.for.hubs=FALSE)
   {
     standardGeneric("random.walk")
   },
   package="diffusr"
 )
 
+
 #' @rdname random-walk-methods
 #' @aliases random.walk,numeric,matrix-method
 setMethod(
   "random.walk",
   signature = signature(p0="numeric", graph="matrix"),
-  function(p0, graph, r=.5, niter=1e4, thresh=1e-4, do.analytical=FALSE, ...)
+  function(p0, graph, r=.5,
+           niter=1e4, thresh=1e-4, do.analytical=FALSE,
+           correct.for.hubs=FALSE)
   {
     p0 <- as.matrix(p0, ncol=1)
-    random.walk(p0, graph, r, niter, thresh, do.analytical, ...)
+    random.walk(p0, graph, r, niter, thresh, do.analytical, correct.for.hubs)
   }
 )
 
@@ -91,9 +103,12 @@ setMethod(
 setMethod(
   "random.walk",
   signature = signature(p0="matrix", graph="matrix"),
-  function(p0, graph, r=.5, niter=1e4, thresh=1e-4, do.analytical=FALSE, ...)
+  function(p0, graph, r=.5,
+           niter=1e4, thresh=1e-4, do.analytical=FALSE,
+           correct.for.hubs=FALSE)
   {
     stopifnot(length(r) == 1)
+    stopifnot(is.logical(correct.for.hubs))
     .check.restart(r)
     .check.starting.matrix(p0)
     .check.graph(graph, p0)
@@ -103,14 +118,17 @@ setMethod(
       message("setting diag of graph to zero")
       diag(graph) <- 0
     }
-
+    if (correct.for.hubs) graph <- hub.correction(graph)
     stoch.graph <- normalize.stochastic(graph)
     if(!.is.ergodic(stoch.graph))
       stop(paste("the provided graph has more than one component.",
                  "It is likely not ergodic."))
 
-    invisible(
-      mrwr_(normalize.stochastic(p0),
-            stoch.graph, r, thresh, niter, do.analytical))
+    l <- list(
+       p.inf=mrwr_(normalize.stochastic(p0),
+                   stoch.graph, r, thresh, niter, do.analytical),
+       transition.matrix=stoch.graph)
+
+    l
   }
 )

diffusr.Rproj

@@ -0,0 +1,20 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: knitr
+LaTeX: pdfLaTeX
+
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source

inst/include/diffusr_RcppExports.h

@@ -40,7 +40,7 @@ namespace diffusr {
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();
         if (rcpp_result_gen.inherits("try-error"))
-            throw Rcpp::exception(as<std::string>(rcpp_result_gen).c_str());
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
         return Rcpp::as<Eigen::MatrixXd >(rcpp_result_gen);
     }
 
@@ -59,7 +59,7 @@ namespace diffusr {
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();
         if (rcpp_result_gen.inherits("try-error"))
-            throw Rcpp::exception(as<std::string>(rcpp_result_gen).c_str());
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
         return Rcpp::as<Eigen::MatrixXd >(rcpp_result_gen);
     }
 
@@ -78,7 +78,45 @@ namespace diffusr {
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();
         if (rcpp_result_gen.inherits("try-error"))
-            throw Rcpp::exception(as<std::string>(rcpp_result_gen).c_str());
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<Eigen::MatrixXd >(rcpp_result_gen);
+    }
+
+    inline std::vector<double> node_degrees_(const Eigen::MatrixXd& W) {
+        typedef SEXP(*Ptr_node_degrees_)(SEXP);
+        static Ptr_node_degrees_ p_node_degrees_ = NULL;
+        if (p_node_degrees_ == NULL) {
+            validateSignature("std::vector<double>(*node_degrees_)(const Eigen::MatrixXd&)");
+            p_node_degrees_ = (Ptr_node_degrees_)R_GetCCallable("diffusr", "_diffusr_node_degrees_");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_node_degrees_(Shield<SEXP>(Rcpp::wrap(W)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<std::vector<double> >(rcpp_result_gen);
+    }
+
+    inline Eigen::MatrixXd hub_normalize_(const Eigen::MatrixXd& W) {
+        typedef SEXP(*Ptr_hub_normalize_)(SEXP);
+        static Ptr_hub_normalize_ p_hub_normalize_ = NULL;
+        if (p_hub_normalize_ == NULL) {
+            validateSignature("Eigen::MatrixXd(*hub_normalize_)(const Eigen::MatrixXd&)");
+            p_hub_normalize_ = (Ptr_hub_normalize_)R_GetCCallable("diffusr", "_diffusr_hub_normalize_");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_hub_normalize_(Shield<SEXP>(Rcpp::wrap(W)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
         return Rcpp::as<Eigen::MatrixXd >(rcpp_result_gen);
     }
 
@@ -97,7 +135,7 @@ namespace diffusr {
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();
         if (rcpp_result_gen.inherits("try-error"))
-            throw Rcpp::exception(as<std::string>(rcpp_result_gen).c_str());
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
         return Rcpp::as<Eigen::MatrixXd >(rcpp_result_gen);
     }
 
@@ -116,7 +154,7 @@ namespace diffusr {
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();
         if (rcpp_result_gen.inherits("try-error"))
-            throw Rcpp::exception(as<std::string>(rcpp_result_gen).c_str());
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
         return Rcpp::as<Rcpp::List >(rcpp_result_gen);
     }