R/BiGraph.R

In keepsimpler/StabEco: Stability of Ecological Systems

library(R6)
library(igraph)
# library(bipartite) # use sortweb function

###############################################################################
#' @title Swap links of bipartite network, that will keep the node degree distribution.
#'
#' @param B incidence matrix of bipartite network, rows and cols represent two groups of nodes/species
#' @param HowManyToTry the times to try for swapping links
#' @return the incidence matrix whose links being randomly swapped.
#' @details .  
swaplinks <- function(B, HowManyToTry = 5000) {
  count1 <- 0
  NumP <- dim(B)[1]
  NumA <- dim(B)[2]
  flag = FALSE
  while (count1 < HowManyToTry){
    count1 <- count1 + 1
    ## pick two rows and two columns
    row1 <- sample(1:NumP, 1)
    row2 <- sample(1:NumP, 1)
    col1 <- sample(1:NumA, 1)
    col2 <- sample(1:NumA, 1)
    ## check swappable
    if (B[row1, col1] == 0.0 && B[row1, col2] != 0.0 && B[row2, col1] != 0.0 && B[row2, col2] == 0.0){
      ## swap
      B[row1, col1] <- B[row1, col2]
      B[row1, col2] <- 0.0
      B[row2, col2] <- B[row2, col1]
      B[row2, col1] <- 0.0
      flag = TRUE
    }
    else{
      if (B[row1, col1] != 0.0 && B[row1, col2] == 0.0 && B[row2, col1] == 0.0 && B[row2, col2] != 0.0){
        ## swap
        B[row1, col2] <- B[row1, col1]
        B[row1, col1] <- 0.0
        B[row2, col1] <- B[row2, col2]
        B[row2, col2] <- 0.0
        flag = TRUE
      }
    }
  }
  if (flag == FALSE) warning('swap links fail!')
  return(B)
}



#' @title bipartite graph
#' @field type type of bipartite graph
#' \describe{
#' \item{bipartite_regular}{a regular bipartite graph, input params: n1, k, restriction: n1 = n2}
#' \item{bipartite_gnm}{a bipartite graph generated by ER model, input params: n1, n2, edges}
#' \item{bipartite_sf}{a bipartite graph generated by fitness model, input params: n1, k, beta, restriction: n1 = n2}
#' \item{bipartite_real}{a bipartite graph gotten from the incidency matrix of an empirical network}
#' \item{bipartite_empty}{an empty bipartite graph}
#' }
#' @examples 
#' BiGraph$new(type = 'bipartite_regular', n1 = 50, k = 5)
BiGraph <- R6Class('BiGraph', 
  inherit = Graph,
  public = list(
    n1 = NULL, # node number in group 1
    n2 = NULL, # node number in group 2
    k = NULL, # average degree of nodes
    m = NULL, # number of (directed) edges
    beta = NULL, # fitness of node, gamma = 1 / beta + 1, gamma is scale-free exponent
    # set graph object from an incidency matrix
    set_graph_inc = function(graph_inc) {
      self$n1 = dim(graph_inc)[1]
      self$n2 = dim(graph_inc)[2]
      self$m = sum(graph_inc)
      self$k = 2 * self$m / (self$n1 + self$n2)
      self$n = self$n1 + self$n2
      private$G = graph_from_incidence_matrix(graph_inc)
      self$type = 'bipartite_real'
    },
    initialize = function(type = c('bipartite_regular', 'bipartite_gnm', 'bipartite_sf', 'bipartite_empty'), 
                          n1 = NULL, n2 = NULL, k = NULL, m = NULL, beta = NULL, ...) {
      type <- match.arg(type)
      switch (type,
              bipartite_regular  = {
                G = sample_k_regular(n1, k, directed = TRUE)
                graph <- as.matrix(as_adjacency_matrix(G))
                private$G = graph_from_incidence_matrix(graph) # tricky!
                self$n1 <- n1
                self$n2 <- n1
                self$type <- type
                self$k <- k
                self$m <- self$k * (self$n1 + self$n2) / 2
                self$n <- self$n1 + self$n2
                # shuffle the rows and cols, to simulate a bipartite regular graph,
                # rather than a unipartite regular graph which have 0 value diagonal elements
                # this permutation does not change the eigenvalues
                # graph <- graph[sample.int(s[1]), sample.int(s[1])]
              },
              bipartite_gnm = {
                private$G = sample_bipartite(n1 = n1, n2 = n2, type = 'gnm', m = m) # m = k * (n1 + n2) / 2
                self$n1 <- n1
                self$n2 <- n2
                self$type <- type
                self$m <- m
                self$k <- round(2 * m / (n1 + n2))
                self$n <- self$n1 + self$n2
              },
              bipartite_sf = {
                stopifnot(! is.null(beta))
                G = sample_fitness(n1 * k, (1:n1)^-beta, (1:n1)^-beta)
                graph <- as.matrix(as_adjacency_matrix(G))
                private$G = igraph::graph_from_incidence_matrix(graph) #, add.names = NA
                self$n1 <- n1
                self$n2 <- n1
                self$type <- type
                self$k <- k
                self$m <- self$k * (self$n1 + self$n2) / 2
                self$n <- self$n1 + self$n2
                self$beta <- beta
              },
              bipartite_empty = {
                # do nothing
                private$G = graph.empty()
              }
      )
    },
    get_graph = function(is_adj = TRUE) {
      #cat('BiGraph')
      if (is_adj == TRUE) {
        graph <- as.matrix(as_adjacency_matrix(private$G))
        return(graph)
      }
      else {
        graph <- as.matrix(as_incidence_matrix(private$G))
        return(graph)
      }
    },
    plot = function() {
      plot(private$G, layout = layout_as_bipartite,
           vertex.color = c("green", "cyan")[V(private$G)$type+1])
    },
    swaplinks = function(HowManyToTry = 5000) {
      graph_inc <- as.matrix(as_incidence_matrix(private$G))
      graph_inc = swaplinks(B = graph_inc, HowManyToTry = HowManyToTry)
      private$G = graph_from_incidence_matrix(graph_inc)
    }
  ))