R 中带有 igraph 的线性阈值模型

Linear Threshold Model in R with igraph

我正在尝试在 R 中实现一个函数,它接受一个图形对象,并且在同一个图形上有两个 "opinions",从随机点开始,我应该研究扩散模式。 [0, 1] 之间有一个顽固性 (alpha) 参数,对于已经感染的节点,意见的变化设置为 (感染邻居的总和)/[总权重之和 * (1+alpha)] > 阈值。

我明白了这个想法,但由于我不是真正的 R 用户,所以我对如何实现它有点迷茫。我在互联网上找到了这个 Cascade and Linear threshold Models,线性阈值模型(未加权,下面的代码)看起来很完美,但我对如何修改代码以满足要求有点迷茫。

library('igraph')

# Independent Cascade model
#   - graph: igraph object
#   - activated: list of initially activated nodes
IC <- function (g, activated) {
  # Defining the graph layout to preserve it
  # the same throughout the visualization
  l <- layout.fruchterman.reingold(g)
  # Setting the activated nodes
  V(g)$activated <- F
  for (v in activated) {
    V(g)[v]$activated <- T
  }
  # Marking all activations (edges) as "not yet tried"
  E(g)$tried <- F
  possible.activations = ecount(g)
  # The process goes on until there are possible activations
  while(possible.activations > 0) {
    # Network visualization (at each simulation step)
    V(g)$color <- ifelse(V(g)$activated, "red", "lightblue")
    plot(g, layout=l, edge.width=E(g)$weight*5)
    # Iterating through activated nodes
    for(v in V(g)) {
      if(V(g)[v]$activated) {
        # Finding activations for each note that have not been tried yet
        for(w in neighbors(g, V(g)[v]$name, mode="out")) {
          e <- E(g)[from(V(g)[v]$name) & to(V(g)[w]$name)]
          if (! e$tried) {
            # Activation attempt
            if (runif(1, 0, 1) <= e$weight) {
              V(g)[w]$activated <- T
            }
            e$tried <- T
            possible.activations <- possible.activations - 1
          }
        }
      }
    }
  }
  # Network visualization after the process has terminated
  V(g)$color <- ifelse(V(g)$activated, "red", "lightblue")
  plot(g, layout=l, edge.width=E(g)$weight*5)
}

以一种相当残酷的方式,但正如我所说,我不是 R 用户也不是粉丝,代码如下:

# Linear Threshold model (unweighted version)
#   - graph: igraph object
#   - activated_1: list of initially activated nodes opinion 1
#   - activated_2: list of initially activated nodes opinion 2
#   - activated_2: alpha value, between 0 and 1
LT <- function (g, activated_1, activated_2, stubborness) {
  # Defining the graph layout to preserve it
  # the same throughout the visualization
  l <- layout.fruchterman.reingold(g)
  # calculating the threshold randomly setting it between 0 and 1
  g <- set.vertex.attribute(g, "threshold", value = stubborness)
  # Setting the activated nodes for each opinion
  g <- set.vertex.attribute(g, "activated", value = 0)
  for (v in activated_1) {
    V(g)[v]$activated <- 1
  }
  for (v in activated_2) {
    V(g)[v]$activated <- 2
  }
  # Indicator of whether simulation should stop
  any.changes <- TRUE
  while(any.changes) {
    # Network visualization (at each simulation step)
    for(v in V(g)) {
      if (V(g)[v]$activated == 0) {
        V(g)[v]$color <- "lightblue"
      }
      if (V(g)[v]$activated == 1) {
        V(g)[v]$color <- "red"
      }
      if (V(g)[v]$activated == 2) {
        V(g)[v]$color <- "green"
      }
    }
    plot(g, layout=l)
    any.changes <- FALSE
    # Iterating through non-activated nodes
    for(v in V(g)) {
        # Calculating the fraction of activated neighbors
        neighborhood <- neighbors(g, V(g)[v]$name)
        activated.neighbors_1 <- 0
        activated.neighbors_2 <- 0
        total.neighbors <- length(neighborhood)
        for(w in neighborhood) {
          if(V(g)[w]$activated == 1) {
            activated.neighbors_1 <- activated.neighbors_1 + 1
          }
          if(V(g)[w]$activated == 2) {
            activated.neighbors_2 <- activated.neighbors_2 + 1
          }
        }
        if (V(g)[v]$activated == 0){
          activation.power_1 <- activated.neighbors_1 / total.neighbors
          activation.power_2 <- activated.neighbors_2 / total.neighbors
        }
        else {
          activation.power_1 <- activated.neighbors_1 / (total.neighbors * (1+V(g)[v]$threshold))
          activation.power_2 <- activated.neighbors_2 / (total.neighbors * (1+V(g)[v]$threshold))
        }
        # Checking if the opinions spread in neutral nodes or opinionated nodes
        if(V(g)[v]$activated == 0){
          if ((activated.neighbors_1 + activated.neighbors_2) != 0){
          if (activation.power_1 > activation.power_2){
            if (activation.power_1 > V(g)[v]$threshold){
              V(g)[v]$activated <- 1
              any.changes <- TRUE
          }}
          else {
            if (activation.power_2 > V(g)[v]$threshold){
              V(g)[v]$activated <- 2
              any.changes <- TRUE
            }
          }
        }
        }
        if(V(g)[v]$activated == 1){
          if ((activated.neighbors_1 + activated.neighbors_2) != 0 &&
              activated.neighbors_1 != total.neighbors){
          if (activation.power_2 > activation.power_1){
            if (activation.power_2 > V(g)[v]$threshold){
              V(g)[v]$activated <- 2
              any.changes <- TRUE
          }
        }
        }}
        if(V(g)[v]$activated == 2){
          if ((activated.neighbors_1 + activated.neighbors_2) != 0 &&
              activated.neighbors_2 != total.neighbors){
          if (activation.power_1 > activation.power_2){
            if (activation.power_1 > V(g)[v]$threshold){
              V(g)[v]$activated <- 1
              any.changes <- TRUE
          }
        }
        }
      }
  # Network visualization after the process has terminated
    }
  }
  for(v in V(g)) {
    if (V(g)[v]$activated == 0) {
      V(g)[v]$color <- "lightblue"
    }
    if (V(g)[v]$activated == 1) {
      V(g)[v]$color <- "red"
    }
    if (V(g)[v]$activated == 2) {
      V(g)[v]$color <- "green"
    }
  }
  plot(g, layout=l)
  # Recap of initial and final number for each opinion
  count_1 <- 0
  count_2 <- 0
  for(v in V(g)) {

    if (V(g)[v]$activated == 1) {
      count_1 <- count_1 + 1
    }
    if (V(g)[v]$activated == 2) {
      count_2 <- count_2 + 1
    }
  }
  cat("Initial number of red was ", length(activated_1), "final number is ", count_1, "\n")
  cat("Initial number of green was ", length(activated_2), "final number is ", count_2, "\n")
}