Fixes for overhaul to 0.1.2

Author: dirmeier

NAMESPACE

@@ -2,8 +2,7 @@
 
 S3method(normalize.laplacian,numeric)
 S3method(normalize.stochastic,numeric)
-export(insulated.heat.diffusion)
-export(laplacian.heat.diffusion)
+export(heat.diffusion)
 export(nearest.neighbors)
 export(normalize.laplacian)
 export(normalize.stochastic)

R/heat_diffusion.R

@@ -53,7 +53,7 @@ setGeneric(
   "heat.diffusion",
   function(h0, graph, t=.5, ...)
   {
-    standardGeneric("laplacian.heat.diffusion")
+    standardGeneric("heat.diffusion")
   },
   package="diffusr"
 )

R/mrw.R

@@ -36,6 +36,16 @@
 #' matrix representing the graph
 #' @param r  a scalar between (0, 1). restart probability if a Markov random
 #' walk with restart is desired
+#' @param thresh  threshold for breaking the iterative computation of the
+#'  stationary distribution. If the absolute difference of the distribution at
+#'  time point $t-1$ and $t$ is less than \code{thresh}, then the algorithm stops.
+#'  If \code{thresh} is not reached before \code{niter}, then the algorithm stops
+#'  as well.
+#' @param niter  maximal number of iterations for computation of the
+#'  Markov chain. If \code{thresh} is not reached, then \code{niter} is used as
+#' 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
 #'
@@ -57,7 +67,7 @@
 #' pt    <- random.walk(p0, graph)
 setGeneric(
   "random.walk",
-  function(p0, graph, r=.5, ...)
+  function(p0, graph, r=.5, niter=1e4, thresh=1e-4, do.analytical=FALSE, ...)
   {
     standardGeneric("random.walk")
   },
@@ -69,10 +79,10 @@ setGeneric(
 setMethod(
   "random.walk",
   signature = signature(p0="numeric", graph="matrix"),
-  function(p0, graph, r=.5, ...)
+  function(p0, graph, r=.5, niter=1e4, thresh=1e-4, do.analytical=FALSE, ...)
   {
     p0 <- as.matrix(p0, ncol=1)
-    random.walk(p0, graph, r, ...)
+    random.walk(p0, graph, r, niter, thresh, do.analytical, ...)
   }
 )
 
@@ -81,23 +91,32 @@ setMethod(
 setMethod(
   "random.walk",
   signature = signature(p0="matrix", graph="matrix"),
-  function(p0, graph, r=.5, ...)
+  function(p0, graph, r=.5, niter=1e4, thresh=1e-4, do.analytical=FALSE, ...)
   {
     stopifnot(length(r) == 1)
     .check.restart(r)
     .check.starting.matrix(p0)
     .check.graph(graph, p0)
+
     if (any(diag(graph) != 0))
     {
       message("setting diag of graph to zero")
       diag(graph) <- 0
     }
+
     stoch.graph <- normalize.stochastic(graph)
     if(!.is.ergodic(stoch.graph))
       stop("the provided graph has more than one component. It is likely not ergodic.")
-    invisible(mrwr_(normalize.stochastic(p0),
-                    stoch.graph,
-                    r)
-              )
+
+    print(p0)
+    print(stoch.graph)
+    print(r)
+    print(thresh)
+    print(niter)
+    print(do.analytical)
+
+    invisible(
+      mrwr_(normalize.stochastic(p0),
+            stoch.graph, r, thresh, niter, do.analytical))
   }
 )

R/util.R

@@ -50,7 +50,7 @@
     stop(paste0("'", name.graph, "' has to be non-negative"))
   if (dim(m)[1] != dim(m)[2])
     stop(paste(name.graph, "has to be of dimension (n x n)!"))
-  if (!is.null(v) && dim(m)[1] != length(v))
+  if (!is.null(v) && dim(m)[1] != dim(v)[1])
     stop(paste(name.vector, "has to have same dimension as", name.graph))
 }
 

man/laplacian-diffusion-methods.Rd man/heat-diffusion-methods.Rdman/laplacian-diffusion-methods.Rd renamed to man/heat-diffusion-methods.Rd

@@ -47,7 +47,8 @@ Eigen::MatrixXd heat_diffusion_(const Eigen::MatrixXd& v0,
   Eigen::VectorXd co =  V.transpose() * v0;
   for (int i = 0; i < co.rows(); ++i)
   {
-    for (int j = 0; j < co.cols(); ++j) {
+    for (int j = 0; j < co.cols(); ++j)
+    {
       if (j % 25 == 0) Rcpp::checkUserInterrupt();
       co(i, j) *= std::exp(-D(i) * t);
     }

man/insulated-diffusion-methods.Rd

@@ -20,8 +20,10 @@
  * along with diffusr. If not, see <http://www.gnu.org/licenses/>.
  */
 
+
 // [[Rcpp::depends(RcppEigen)]]
 #include <RcppEigen.h>
+// [[Rcpp::plugins(cpp11)]]
 #include <cmath>
 
 //' Column normalize a matrix, so that it is stochastic.
@@ -42,7 +44,7 @@ Eigen::MatrixXd stoch_col_norm_(const Eigen::MatrixXd& W)
     if ((W.col(i)).sum() <= zero_col) res.col(i).fill(empt_col_val);
     else res.col(i) = W.col(i) / colsums(i);
   }
-  
+
   return res;
 }
 

man/random-walk-methods.Rd

@@ -23,6 +23,7 @@
 // [[Rcpp::depends(RcppEigen)]]
 #include <RcppEigen.h>
 #include <vector>
+#include <iostream>
 
 //' Do a Markon random walk (with restart) on an column-normalised adjacency
 //' matrix.
@@ -49,7 +50,7 @@ Eigen::MatrixXd mrwr_(const Eigen::MatrixXd& p0,
     Eigen::MatrixXd pt;
     if (!do_analytical)
     {
-      Eigen::MatrixXd pt = p0;
+      pt = p0;
       Eigen::MatrixXd pold;
       int iter   = 0;
       do
@@ -62,8 +63,9 @@ Eigen::MatrixXd mrwr_(const Eigen::MatrixXd& p0,
     }
     else
     {
+      Eigen::MatrixXd I = Eigen::MatrixXd::Identity(W.rows(), W.cols());
       Eigen::MatrixXd T  = r * (I - (1 - r) * W).inverse();
-      pt = T %*% p0;
+      pt = T * p0;
     }
 
     return pt;