# 使用R模拟网络扩散

1. 生成一个网络:g(V, E)。
2. 随机选择一个或几个节点作为种子（seeds）。
3. 每个感染者以概率p（可视作该节点的传染能力,通常表示为$\beta$）影响与其相连的节点。

require(igraph)
# generate a social graph
size = 50

# 规则网
g = graph.tree(size, children = 2); plot(g)
g = graph.star(size); plot(g)
g = graph.full(size); plot(g)
g = graph.ring(size); plot(g)
g = connect.neighborhood(graph.ring(size), 2); plot(g) # 最近邻耦合网络

# 随机网络
g = erdos.renyi.game(size, 0.1)

# 小世界网络
g = rewire.edges(erdos.renyi.game(size, 0.1), prob = 0.8 )
# 无标度网络
g = barabasi.game(size) ; plot(g)


#  initiate the diffusers
seeds_num = 1
set.seed(2014); diffusers = sample(V(g),seeds_num) ; diffusers
infected =list()
infected[[1]]= diffusers


# for example, set percolation probability = 0.5
coins = c(0,1)
n = length(coins)
sample(coins, 1, replace=TRUE, prob=rep(1/n, n))


p = 0.128
coins = c(rep(1, p*1000), rep(0,(1-p)*1000))
n = length(coins)
sample(coins, 1, replace=TRUE, prob=rep(1/n, n))


# function for updating the diffusers
update_diffusers = function(diffusers){
nearest_neighbors = neighborhood(g, 1, diffusers)
nearest_neighbors = data.frame(table(unlist(nearest_neighbors)))
nearest_neighbors = subset(nearest_neighbors, !(nearest_neighbors[,1]%in%diffusers))
# toss the coins
toss = function(freq) {
tossing = NULL
for (i in 1:freq ) tossing[i] = sample(coins, 1, replace=TRUE, prob=rep(1/n, times=n))
tossing = sum(tossing)
return (tossing)
}
keep = unlist(lapply(nearest_neighbors[,2], toss))
new_infected = as.numeric(as.character(nearest_neighbors[,1][keep >= 1]))
diffusers = unique(c(diffusers, new_infected))
return(diffusers)
}


## Start the contagion!
i = 1
while(length(infected[[i]]) < size){
infected[[i+1]] = sort(update_diffusers(infected[[i]]))
cat(length(infected[[i+1]]), "\n")
i = i + 1
}


# "growth_curve"
num_cum = unlist(lapply(1:i, function(x) length(infected［x］) ))
p_cum = num_cum/max(num_cum)
time = 1:i

png(file = "./temporal_growth_curve.png",
width=5, height=5,
units="in", res=300)
plot(p_cum~time, type = "b")
dev.off()


# generate a palette
E(g)$color = "blueviolet" V(g)$color = "white"
set.seed(2014); layout.old = layout.fruchterman.reingold(g)
V(g)$color[V(g)%in%diffusers] = "red" plot(g, layout =layout.old)  使用谢益辉开发的animation的R包可视化。 library(animation) saveGIF({ ani.options(interval = 0.5, convert = shQuote("C:/Program Files/ImageMagick-6.8.8-Q16/convert.exe")) # start the plot m = 1 while(m <= length(infected)){ V(g)$color = "white"
V(g)$color[V(g)%in%infected[[m]]] = "red" plot(g, layout =layout.old) m = m + 1} })  如同在Netlogo里一样，我们可以把网络扩散与增长曲线同时展示出来： saveGIF({ ani.options(interval = 0.5, convert = shQuote("C:/Program Files/ImageMagick-6.8.8-Q16/convert.exe")) # start the plot m = 1 while(m <= length(infected)){ # start the plot layout(matrix(c(1, 2, 1, 3), 2,2, byrow = TRUE), widths=c(3,1), heights=c(1, 1)) V(g)$color = "white"
V(g)\$color[V(g)%in%infected[[m]]] = "red"
num_cum = unlist(lapply(1:m, function(x) length(infected[[x]]) ))
p_cum = num_cum/size
p = diff(c(0, p_cum))
time = 1:m
plot(g, layout =layout.old, edge.arrow.size=0.2)
title(paste("Scale-free Network \n Day", m))
plot(p_cum~time, type = "b", ylab = "CDF", xlab = "Time",
xlim = c(0,i), ylim =c(0,1))
plot(p~time, type = "h", ylab = "PDF", xlab = "Time",
xlim = c(0,i), ylim =c(0,1), frame.plot = FALSE)
m = m + 1}
}, ani.width = 800, ani.height = 500)


Tags:

Categories:

Updated: