R/sample_kcores.R

In schochastics/netUtils: A Collection of Tools for Network Analysis

#' Generate random graphs with a given coreness sequence
#' @description Similar to \link[igraph]{sample_degseq} just with \link[igraph]{coreness}
#' @param cores coreness sequence
#' @details The code is an adaption of the python code from https://github.com/ktvank/Random-Graphs-with-Prescribed-K-Core-Sequences/
#' @return igraph object of graph with the same coreness sequence as the input
#' @references
#' Van Koevering, Katherine, Austin R. Benson, and Jon Kleinberg. 2021. ‘Random Graphs with Prescribed K-Core Sequences: A New Null Model for Network Analysis’. ArXiv:2102.12604. https://doi.org/10.1145/3442381.3450001.
#' @author David Schoch
#' @examples
#' library(igraph)
#' g1 <- make_graph("Zachary")
#' kcores1 <- coreness(g1)
#' g2 <- sample_coreseq(kcores1)
#' kcores2 <- coreness(g2)
#'
#' # the sorted arrays are the same
#' all(sort(kcores1) == sort(kcores2))
#'
#' @export
sample_coreseq <- function(cores) {
    cores <- sort(cores, decreasing = TRUE)

    if (!is_kcoreseq(cores)) {
        stop("sequence is not a valid kcore sequence")
    }
    K <- cores[1]

    indices <- lapply(0:K, function(x) which(cores == x))

    if (K == 0) {
        n0 <- length(indices[[1]])
        g <- igraph::make_empty_graph(n = n0, directed = FALSE)
        return(g)
    } else if (K == 1) {
        n0 <- length(indices[[1]])
        n1 <- length(indices[[2]])
        g <- igraph::make_empty_graph(n = n1, directed = FALSE)
        g <- igraph::add_edges(g, edges = c(t(cbind(1:(n1 - 1), 2:n1))))
        g <- igraph::add_vertices(g, nv = n0)
        return(g)
    } else if (K == 2) {
        n0 <- length(indices[[1]])
        n1 <- length(indices[[2]])
        n2 <- length(indices[[3]])
        g <- igraph::graph.empty(n = n2, directed = FALSE)
        g <- igraph::add_edges(g, edges = c(t(cbind(1:(n2 - 1), 2:n2))))
        g <- igraph::add_edges(g, edges = c(1, n2))
        g <- add_lower_nodes(cores, indices, g, 1)
        return(g)
    } else {
        g <- igraph::graph.empty(directed = FALSE)
        g <- generate_k_graph(K, sum(cores == K), g)
        g <- add_lower_nodes(cores, indices, g, K - 1)
        return(g)
    }
}

is_kcoreseq <- function(cores) {
    cores <- sort(cores, decreasing = TRUE)
    num_max_core_val <- sum(cores == max(cores))
    return(num_max_core_val >= cores[1])
}

add_lower_nodes <- function(cores, indices, g, k) {
    K <- cores[1]
    higher_nodes <- sort(unname(unlist(indices[(k + 2):(K + 1)])))
    g <- igraph::add_vertices(g, nv = length(indices[[k + 1]]))
    if (k == 0) {
        return(g)
    }
    for (v in indices[[k + 1]]) {
        randv <- sample(higher_nodes, k, replace = FALSE)
        g <- igraph::add_edges(g, c(t(cbind(randv, v))))
    }
    return(add_lower_nodes(cores, indices, g, k - 1))
}

generate_k_graph <- function(C, N, g) {
    if (!((C %% 2) == (N %% 2))) {
        g <- igraph::add_vertices(g, nv = N)
        g <- igraph::add_edges(g, edges = c(t(cbind(1:(N - 1), 2:N))))
        g <- igraph::add_edges(g, edges = c(1, N))
        z <- ceiling((N - C + 1) / 2)
        for (i in 0:(N - 1)) {
            start <- (i + z + 1) %% (N) # BUG POTENTIAL!!!
            stop <- (i - z) %% (N) # BUG POTENTIAL!!!

            if (stop >= start) {
                listv <- start:(stop - 1) # BUG POTENTIAL!!!
                edges <- c(t(cbind(listv, i)))
            } else {
                listv <- c((0:(N - 1))[0:(stop)], (0:(N - 1))[(start + 1):N])
                edges <- c(t(cbind(listv, i)))
            }
            g <- igraph::add_edges(g, edges + 1)
            g <- igraph::simplify(g)
        }
    } else {
        g <- generate_k_graph(C, N - 1, g)
        g <- igraph::add_vertices(g, 1)
        randv <- sample(1:(N - 1), C)
        g <- igraph::add_edges(g, c(t(cbind(randv, N))))
    }
    return(g)
}