Functions to generate geometric networks: random geometric graphs, connected random geometric graphs, random geometric trees.

networks.R

# Functions to generate networks ("graphs") for spatial ecology.
#
# Stores the graph in a rather inefficient (but easy-to-use) adjacency matrix
# of boolean values.
#
# Key functions:
#  - geograph returns a random geometric graph.
#  - cgeograph returns a connected random geometric graph.
#  - geotree returns a random geometric tree.
#
# Warning: not idiomatic R.
#
# By: Philippe Desjardins-Proulx ([email protected])
# License: MIT

geograph = function(n = 64, r = 0.32) {
  # Generates a random geometric graph.
  #
  # Args:
  #   n: Number of vertices.
  #   r: Threshold distance to connect vertices.
  #
  # Returns:
  #   A list with the xy coordinates and the adjacency matrix.
  xy <- cbind(runif(n), runif(n))
  distMat <- as.matrix(dist(xy, method = 'euclidean', upper = T, diag = T))
  adjMat <- matrix(F, nr = n, nc = n)
  adjMat[distMat < r] <- T
  diag(adjMat) <- F
  return(list(xy, adjMat))
}

cgeograph = function(n = 64, r = 0.32) {
  # Generates a connected random geometric graph.
  #
  # Args:
  #   n: Number of vertices.
  #   r: Threshold distance to connect vertices.
  #
  # Returns:
  #   A list with the xy coordinates, the adjacency matrix, and
  # the number of attempts that was required to find a connected
  # graph.
  attempts <- 0
  repeat {
    attempts <- attempts + 1
    g <- geograph(n, r)
    if (isConnected(g[[2]])) {
      return(append(g, attempts))
    }
  }
}

geotree = function(n = 64, r = 0.32) {
  # Generates a random geometric tree.
  #
  # Args:
  #   n: Number of vertices.
  #   r: Threshold distance to connect vertices.
  #
  # Returns:
  #   A list with the xy coordinates, the adjacency matrix, the number of attempts
  # to build the geometric graph, and the ID of the source.
  source <- sample(1:n, 1)
  g <- cgeograph(n, r)
  return(list(g[[1]], spt(g[[2]], source), g[[3]], source))
}

plotGeoGraph = function(graph) {
  # Simple plot for geometric graph. Adapted from code by Hedvig Nenzen.
  #
  # Args:
  #   g: A list with g[[1]] for the coordinate and g[[2]] for the adjacency
  #      matrix.
  xy = g[[1]]
  adjMat = g[[2]]
  plot(xy[,1], xy[,2], ylab = '', axes = F, cex.lab = 1.5,
       xlim = c(0,1), xlab = '', main = '', family = 'sans', col = 1)
  m = stack(as.data.frame(adjMat))[,1]
  xx <- expand.grid(xy[,1], xy[,1])
  yy <- expand.grid(xy[,2], xy[,2])
  xx <- subset(xx, m == T)
  yy <- subset(yy, m == T)
  arrows(x0 = xx[,1], x1 = xx[,2], y0 = yy[,1], y1 = yy[,2], length = 0, lwd = 0.5)
  points(xy[,1], xy[,2], pch = 21, cex = 2, bg = 1)
}

isConnected = function(adjMat) {
  # Tests if the graph is fully connected (i.e.: there is a path between all nodes).
  #
  # Args:
  #   adjMat: An adjacency matrix (made of boolean values).
  #
  # Returns:
  #   TRUE if the graph is connected, FALSE otherwise.
  diag(adjMat) <- F
  n = nrow(adjMat)
  for (v in 1:n) {
    inPath <- vector('logical', n)
    inPath[v] <- T
    inPath <- testCon(adjMat, inPath, v)
    if (sum(inPath == T) < n) {
      return(F)
    }
  }
  return(T)
}

testCon = function(adjMat, inPath, v) {
  # Helper recursive function for 'isConnected'.
  for (i in 1:nrow(adjMat)) {
    if (adjMat[v, i] && inPath[i] == F) {
      inPath[i] <- T
      inPath <- testCon(adjMat, inPath, i)
    }
  }
  return(inPath)
}

spt = function(adjMat0, source = 1) {
  # Generates the graph of the shortest path tree using the Dijkstra
  # algorithm for a given graph and source vertex.
  #
  # Args:
  #   adjMat0: The adjacency matrix for the graph.
  #   source: The starting vertex for the algorithm.
  #
  # Returns:
  #   The shortest path tree as an adjacency matrix (a matrix of bool).
  n <- nrow(adjMat0)
  distance <- rep(Inf, n)
  distance[source] <- 0
  previous <- rep(0, n)
  q <- 1:n

  while (length(q) != 0) {
    u <- q[1]
    for (i in q) {
      if (distance[i] < distance[u]) {
        u <- i
      }
    }
    q <- setdiff(q, u)
    for (i in q) {
      if (adjMat0[u, i] && distance[u] + 1 < distance[i]) {
        distance[i] <- distance[u] + 1
        previous[i] <- u
      }
    }
  }
  # Build the adjacency matrix:
  adjMat <- matrix(F, nr = n, nc = n)
  for (i in 1:n) {
    adjMat[i, previous[i]] <- adjMat[previous[i], i] <- T
  }
  return(adjMat)
}