R/RR.R

#' Randomness Recycler 
#' 
#' Randomness recycler method for sampling from the Random Cluster model. Introduced by
#' Fill and Huber (2000).
#' 
#' @param n.nodes The number of nodes in the network
#' @param p The wiring probability
#' @param q The clustering coefficient of the Random Cluster model
#' 
#' @return A adjacency matrix of a graph sampled from the Random Cluster model.
#' 
#' @export
#' 


RandomnessRecycler <- function(n.nodes, p, q) {
  
  i <- 0
  
  # First create the empty set
  x.mat <- matrix(0, n.nodes, n.nodes)
  
  # Than create Et
  Et <- matrix(NA,nrow = 1, ncol = 2)
  
  # Create E
  E <- which(upper.tri(x.mat) == TRUE, arr.ind=T)


  # Start while loop
  #!all(apply(E, 1, paste, collapse = "") %in% 
   #      apply(Et, 1, paste, collapse = ""))
  while (!all(apply(E, 1, paste, collapse = "") %in% 
                    apply(Et, 1, paste, collapse = ""))) {
    i <- i + 1
    # Sample an edge
   
    e <- sample(which(!apply(E, 1, paste, collapse = "") %in% 
                        apply(Et, 1, paste, collapse = "")),1)
    
    # Make graph object
    graph <- graph_from_adjacency_matrix(x.mat, mode = "undirected")
    
    # Sample state
    state <- sample(c(0,1),1, prob = c((1-p),p))
    # Set state
    x.mat[E[e,1], E[e,2]] <- x.mat[E[e,2], E[e,1]] <- state
  
    # Get in which component the nodes are
    comp <- components(graph)$membership
    
    if (state == 0) {
      # Accept state
      Et <- rbind(Et, E[e, ])
    }
    
    if (state == 1 ) {
      if (comp[E[e,1]] == comp[E[e,2]]) {
        # If connected, accept the proposed edge
        Et <- rbind(Et, E[e, ])
      } else if (((1/q) >= runif(1))) {
        # Unconnected with prob 1/q Accept the edge
        Et <- rbind(Et, E[e, ])
      } else {
        # Reject; but we have information about the  clusters
        
        # Sample w
        w <- sample(c(E[e,1], E[e,2]),1 )
        # Get all vertices in same component W
        comp.w <- which(comp == comp[w], arr.ind = T)
        
        # Delete from Et all edges in or adjacent to W
        # Colour them zero in x
        # Set edge to zero
        x.mat[Et[apply(Et, 1, function(r) any(r %in% as.vector(comp.w))),, drop = FALSE]] <- 0
        x.mat <- SymMat(x.mat)
        
        # Create a spanning tree of the nodes that are deleted
        
        # Make a edgelist
        x.edge <- as_edgelist(graph)
        a <- x.edge[apply(x.edge, 1, function(r) any(r %in% comp.w)),, drop = FALSE]
        # Add edge (v,w) to the spanning tree
        a <- rbind(a, E[e,])
        # Create graph
        a <- graph_from_edgelist(a, directed = FALSE)
        # use mst to create uinique spanning tree of the edges
        a <- mst(a)
        # Change back into edge list
        a <- as_edgelist(a)
     
        # Delete all edges adjacent to or in C
        Et <- Et[!apply(Et, 1, function(r) any(r %in% w)),, drop = FALSE]
        Et <- Et[!apply(Et, 1, function(r) any(r %in% as.vector(comp.w))),, drop = FALSE]
        
        # Sample edges from the spanning tree
        rho <- p/q
        t <- sample(c(1,0), size = nrow(a), replace = TRUE, 
                 prob = c(rho/(1-p+rho),(1-p)/(1-p+rho)))
        
        # Set edge states
        x.mat[a] <- t
        x.mat <- SymMat(x.mat)

        
        # Add the spanning tree back into Et
        Et <- rbind(Et, a)
        
      }
    }
    
    if(!isSymmetric(x.mat)) {
      print("not symmetric at")
      print(i)
    }
  }
  
  return(x.mat)
}


  
TimoMatzen/RandomCluster documentation built on May 15, 2019, 11:48 a.m.