refactor constrained layouts

NAMESPACE

@@ -25,6 +25,7 @@ export(layout_with_centrality_group)
 export(layout_with_constrained_stress)
 export(layout_with_constrained_stress3D)
 export(layout_with_eigen)
+export(layout_with_fixed_coords)
 export(layout_with_focus)
 export(layout_with_focus_group)
 export(layout_with_pmds)

NEWS.md

@@ -4,6 +4,7 @@
   work for disconnected graphs
 * internal code refactoring
 * added `layout_as_metromap()`
+* added `layout_with_fixed_coords()`
 
 # graphlayouts 1.0.2
 

R/RcppExports.R

@@ -9,6 +9,10 @@ constrained_stress_major <- function(y, dim, W, D, iter, tol) {
     .Call(`_graphlayouts_constrained_stress_major`, y, dim, W, D, iter, tol)
 }
 
+fixed_stress_major <- function(y, fixedCoords, W, D, iter, tol) {
+    .Call(`_graphlayouts_fixed_stress_major`, y, fixedCoords, W, D, iter, tol)
+}
+
 constrained_stress3D <- function(x, W, D) {
     .Call(`_graphlayouts_constrained_stress3D`, x, W, D)
 }

R/layout_stress.R

@@ -438,6 +438,77 @@ layout_with_constrained_stress3D <- function(g, coord, fixdim = "x", weights = N
         ))
     }
 }
+
+#------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------#
+#' Layout with fixed coordinates
+#'
+#' @name layout_fixed_coords
+#' @description force-directed graph layout based on stress majorization with
+#' fixed coordinates for some nodes
+#' @param g igraph object
+#' @param coords numeric n x 2 matrix, where n is the number of nodes. values
+#' are either NA or fixed coordinates. coordinates are only calculated for the
+#' NA values.
+#' @param weights possibly a numeric vector with edge weights. If this is NULL and the graph has a weight edge attribute, then the attribute is used. If this is NA then no weights are used (even if the graph has a weight attribute). By default, weights are ignored. See details for more.
+#' @param iter number of iterations during stress optimization
+#' @param tol stopping criterion for stress optimization
+#' @param mds should an MDS layout be used as initial layout (default: TRUE)
+#' @param bbox constrain dimension of output. Only relevant to determine the placement of disconnected graphs
+#' @details Be careful when using weights. In most cases, the inverse of the edge weights should be used to ensure that the endpoints of an edges with higher weights are closer together (weights=1/E(g)$weight).
+#'
+#' The layout_igraph_* function should not be used directly. It is only used as an argument for plotting with 'igraph'.
+#' 'ggraph' natively supports the layout.
+#' @seealso [layout_constrained_stress]
+#' @return matrix of xy coordinates
+#' @examples
+#' library(igraph)
+#' set.seed(12)
+#' g <- sample_bipartite(10, 5, "gnp", 0.5)
+#' fxy <- cbind(c(rep(0, 10), rep(1, 5)), NA)
+#' xy <- layout_with_fixed_coords(g, fxy)
+#' @export
+layout_with_fixed_coords <- function(g, coords, weights = NA,
+                                     iter = 500, tol = 0.0001, mds = TRUE, bbox = 30) {
+    ensure_igraph(g)
+
+    oldseed <- get_seed()
+    set.seed(42) # stress is deterministic and produces same result up to translation. This keeps the layout fixed
+    on.exit(restore_seed(oldseed))
+
+    if (missing(coords)) {
+        stop('"coords" is missing with no default.')
+    }
+    if (nrow(coords) != igraph::vcount(g) && ncol(coords) != 2) {
+        stop("coords has the wrong dimensions")
+    }
+    if (all(!is.na(coords))) {
+        warning("all coordinates fixed")
+        return(coords)
+    }
+    comps <- igraph::components(g, "weak")
+    if (comps$no == 1) {
+        if (igraph::vcount(g) == 1) {
+            return(matrix(c(0, 0), 1, 2))
+        } else {
+            D <- igraph::distances(g, weights = weights)
+            W <- 1 / D^2
+            diag(W) <- 0
+            n <- igraph::vcount(g)
+            xinit <- .init_layout(g, D, mds, n, dim = 2)
+            not_na_idx <- which(!is.na(coords), arr.ind = TRUE)
+            xinit[not_na_idx] <- coords[not_na_idx]
+            return(fixed_stress_major(xinit, coords, W, D, iter, tol))
+        }
+    } else {
+        # return(.component_layouter(
+        #     g = g, weights = weights, comps = comps, dim = 2, mds = mds,
+        #     bbox = bbox, iter = iter, tol = tol, FUN = constrained_stress_major, fixdim = fixdim, coord = coord
+        # ))
+        stop("only connected graphs are supported")
+    }
+}
+
 #' radial focus group layout
 #'
 #' @description arrange nodes in concentric circles around a focal node according to their distance from the focus and keep predefined groups in the same angle range.

src/RcppExports.cpp

@@ -40,6 +40,22 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// fixed_stress_major
+NumericMatrix fixed_stress_major(NumericMatrix y, NumericMatrix fixedCoords, NumericMatrix W, NumericMatrix D, int iter, double tol);
+RcppExport SEXP _graphlayouts_fixed_stress_major(SEXP ySEXP, SEXP fixedCoordsSEXP, SEXP WSEXP, SEXP DSEXP, SEXP iterSEXP, SEXP tolSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type y(ySEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type fixedCoords(fixedCoordsSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type W(WSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type D(DSEXP);
+    Rcpp::traits::input_parameter< int >::type iter(iterSEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    rcpp_result_gen = Rcpp::wrap(fixed_stress_major(y, fixedCoords, W, D, iter, tol));
+    return rcpp_result_gen;
+END_RCPP
+}
 // constrained_stress3D
 double constrained_stress3D(NumericMatrix x, NumericMatrix W, NumericMatrix D);
 RcppExport SEXP _graphlayouts_constrained_stress3D(SEXP xSEXP, SEXP WSEXP, SEXP DSEXP) {
@@ -267,6 +283,7 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_graphlayouts_constrained_stress", (DL_FUNC) &_graphlayouts_constrained_stress, 3},
     {"_graphlayouts_constrained_stress_major", (DL_FUNC) &_graphlayouts_constrained_stress_major, 6},
+    {"_graphlayouts_fixed_stress_major", (DL_FUNC) &_graphlayouts_fixed_stress_major, 6},
     {"_graphlayouts_constrained_stress3D", (DL_FUNC) &_graphlayouts_constrained_stress3D, 3},
     {"_graphlayouts_constrained_stress_major3D", (DL_FUNC) &_graphlayouts_constrained_stress_major3D, 6},
     {"_graphlayouts_criterion_angular_resolution", (DL_FUNC) &_graphlayouts_criterion_angular_resolution, 2},

src/constrained_stress.cpp

@@ -2,79 +2,127 @@
 using namespace Rcpp;
 
 // [[Rcpp::export]]
-double constrained_stress(NumericMatrix x, NumericMatrix W, NumericMatrix D){
-  double fct=0;
-  int n=x.nrow();
-  for(int i=0;i<(n-1);++i){
-    for(int j=(i+1);j<n;++j){
-      double denom = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+(x(i,1)-x(j,1))*(x(i,1)-x(j,1)));
-      fct+=W(i,j)*(denom-D(i,j))*(denom-D(i,j));
+double constrained_stress(NumericMatrix x, NumericMatrix W, NumericMatrix D) {
+  double fct = 0;
+  int n = x.nrow();
+  for (int i = 0; i < (n - 1); ++i) {
+    for (int j = (i + 1); j < n; ++j) {
+      double denom = sqrt((x(i, 0) - x(j, 0)) * (x(i, 0) - x(j, 0)) +
+                          (x(i, 1) - x(j, 1)) * (x(i, 1) - x(j, 1)));
+      fct += W(i, j) * (denom - D(i, j)) * (denom - D(i, j));
     }
   }
   return fct;
 }
 
 // [[Rcpp::export]]
-NumericMatrix constrained_stress_major(NumericMatrix y,
-                           int dim,
-                           NumericMatrix W,
-                           NumericMatrix D,
-                           int iter,
-                           double tol) {
+NumericMatrix constrained_stress_major(NumericMatrix y, int dim,
+                                       NumericMatrix W, NumericMatrix D,
+                                       int iter, double tol) {
   int n = y.nrow();
+  NumericMatrix x(clone(y));
 
-  NumericMatrix x(n,2);
-  for(int i=0;i<n;++i){
-    for(int d=0;d<2;++d){
-      x(i,d)=y(i,d);
+  NumericVector wsum = rowSums(W);
+
+  double stress_old = constrained_stress(x, W, D);
+
+  for (int k = 0; k < iter; ++k) {
+    NumericMatrix xnew(n, 2);
+    xnew(_, dim - 1) = y(_, dim - 1);
+
+    for (int i = 0; i < n; ++i) {
+      for (int j = 0; j < n; ++j) {
+        if (i != j) {
+          double denom = sqrt(sum(pow(x(i, _) - x(j, _), 2)));
+          if (denom > 0.00001) {
+            int updateDim = 1 - (dim - 1);
+            xnew(i, updateDim) +=
+                W(i, j) *
+                (x(j, updateDim) +
+                 D(i, j) * (x(i, updateDim) - x(j, updateDim)) / denom);
+          }
+        }
+      }
+      int updateDim = 1 - (dim - 1);
+      xnew(i, updateDim) /= wsum[i];
     }
-  }
 
-  NumericVector wsum(n);
-  for(int i=0;i<n;++i){
-    for(int j=0;j<n;++j){
-      wsum[i]+=W(i,j);
+    double stress_new = constrained_stress(xnew, W, D);
+    double eps = (stress_old - stress_new) / stress_old;
+
+    if (eps <= tol) {
+      break;
     }
+
+    stress_old = stress_new;
+    x = clone(xnew);
   }
+  return x;
+}
 
-  double stress_old = constrained_stress(x,W,D);
+NumericMatrix replaceNA(NumericMatrix x, NumericMatrix y) {
+  int nrow = x.nrow(), ncol = x.ncol();
 
-  for(int k=0; k<iter; ++k){
-    NumericMatrix xnew(n,2); //out or in?
-    for(int i=0;i<n;++i){
-      if(dim==1){
-        xnew(i,0) = y(i,0);
-      } else{
-        xnew(i,1) = y(i,1);
+  for (int i = 0; i < nrow; ++i) {
+    for (int j = 0; j < ncol; ++j) {
+      if (!NumericMatrix::is_na(y(i, j))) {
+        x(i, j) = y(i, j);
       }
-      for(int j=0; j<n;++j){
-        if(i!=j){
-          double denom = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+(x(i,1)-x(j,1))*(x(i,1)-x(j,1)));
-          if(denom>0.00001){
-            if(dim==2){
-              xnew(i,0) += W(i,j)*(x(j,0)+D(i,j)*(x(i,0)-x(j,0))/denom);
-            } else{
-              xnew(i,1) += W(i,j)*(x(j,1)+D(i,j)*(x(i,1)-x(j,1))/denom);
+    }
+  }
+
+  return x;
+}
+
+// [[Rcpp::export]]
+NumericMatrix fixed_stress_major(NumericMatrix y, NumericMatrix fixedCoords,
+                                 NumericMatrix W, NumericMatrix D, int iter,
+                                 double tol) {
+  int n = y.nrow();
+  NumericMatrix x(clone(y));
+
+  NumericVector wsum = rowSums(W);
+
+  double stress_old = constrained_stress(x, W, D);
+
+  for (int k = 0; k < iter; ++k) {
+    NumericMatrix xnew(n, 2);
+    std::fill(xnew.begin(), xnew.end(), 0);
+    xnew = replaceNA(xnew, fixedCoords);
+
+    for (int i = 0; i < n; ++i) {
+      for (int j = 0; j < n; ++j) {
+        if (i != j) {
+          double denom = sqrt(sum(pow(x(i, _) - x(j, _), 2)));
+          if ((denom > 0.00001)) {
+            if (!NumericMatrix::is_na(fixedCoords(i, 0))) {
+              xnew(i, 0) +=
+                  W(i, j) * (x(j, 0) + D(i, j) * (x(i, 0) - x(j, 0)) / denom);
+            }
+            if (!NumericMatrix::is_na(fixedCoords(i, 1))) {
+              xnew(i, 1) +=
+                  W(i, j) * (x(j, 1) + D(i, j) * (x(i, 1) - x(j, 1)) / denom);
             }
           }
         }
       }
-      if(dim==2){
-        xnew(i,0) = xnew(i,0)/wsum[i];
-      } else{
-        xnew(i,1) = xnew(i,1)/wsum[i];
+      if (!NumericMatrix::is_na(fixedCoords(i, 0))) {
+        xnew(i, 0) = xnew(i, 0) / wsum[i];
+      }
+      if (!NumericMatrix::is_na(fixedCoords(i, 1))) {
+        xnew(i, 1) = xnew(i, 1) / wsum[i];
       }
     }
-    double stress_new=constrained_stress(xnew,W,D);
-    double eps=(stress_old-stress_new)/stress_old;
-    if(eps<= tol){
+
+    double stress_new = constrained_stress(xnew, W, D);
+    double eps = (stress_old - stress_new) / stress_old;
+
+    if (eps <= tol) {
       break;
     }
-    stress_old=stress_new;
-    for(int i=0;i<n;++i){
-      x(i,0) = xnew(i,0);
-      x(i,1) = xnew(i,1);
-    }
+
+    stress_old = stress_new;
+    x = clone(xnew);
   }
   return x;
-}
+}