R語言模擬疫情傳播圖RVirusBroadcast展示疫情資料

  ###name:RVirusBroadcast 
  ###author:hxj7(hxj5hxj5@126.com)  
  ###version:202002010  
  ###note:本程式是"VirusBroadcast (in Java)"的R版本  
  ###      VirusBroadcast (in Java) 專案連結：
  ###      https://github.com/KikiLetGo/VirusBroadcast/tree/master/src  
  library(tibble)  
  library(dplyr) 
  ########## 模擬引數 ########## 
  ORIGINAL_COUNT <- 50     # 初始感染數量 
  BROAD_RATE <- 0.8        # 傳播率 
  SHADOW_TIME <- 140       # 潛伏時間，14天為140 
  HOSPITAL_RECEIVE_TIME <- 10   # 醫院收治響應時間 
  BED_COUNT <- 1000        # 醫院床位 
  MOVE_WISH_MU <- -0.99   # 流動意向平均值，建議調整範圍：[-0.99,0.99]; 
                       #   -0.99 人群流動最慢速率，甚至完全控制疫情傳播; 
                       #   0.99為人群流動最快速率, 可導致全城感染 
  CITY_PERSON_SIZE <- 5000    # 城市總人口數量 
FATALITY_RATE <- 0.02       # 病死率，根據2月6日資料估算（病死數/確診數）為0.02 
  SHADOW_TIME_SIGMA <- 25     # 潛伏時間方差 
  CURED_TIME <- 50            # 治癒時間均值，從入院開始計時 
  CURED_SIGMA <- 10           # 治癒時間標準差 
  DIE_TIME <- 300             # 死亡時間均值，30天，從發病（確診）時開始計時 
  DIE_SIGMA <- 50             # 死亡時間標準差 
  CITY_WIDTH <- 700           # 城市大小即視窗邊界，限制不允許出城 
  CITY_HEIGHT <- 800 
  MAX_TRY <- 300             # 最大模擬次數，300代表30天 
  ########## 生成人群點，用不同顏色代表不同健康狀態。 ########## 
  # 用正態分佈刻畫人群點的分佈 
  CITY_CENTERX <- 400         # x軸的mu值 
  CITY_CENTERY <- 400 
  PERSON_DIST_X_SIGMA <- 100  # x軸的sigma值 
  PERSON_DIST_Y_SIGMA <- 100 
  # 市民狀態應該需要細分，雖然有的狀態暫未納入模擬，但是細分狀態應該保留 
  STATE_NORMAL <- 0            # 正常人，未感染的健康人 
  STATE_SUSPECTED <- STATE_NORMAL + 1   # 有暴露感染風險 
  STATE_SHADOW <- STATE_SUSPECTED + 1   # 潛伏期 
  STATE_CONFIRMED <- STATE_SHADOW + 1   # 發病且已確診為感染病人 
  STATE_FREEZE <- STATE_CONFIRMED + 1   # 隔離治療，禁止位移 
  STATE_DEATH <- STATE_FREEZE + 1    # 病死者 
  STATE_CURED <- STATE_DEATH + 1   # 治癒數量用於計算治癒出院後歸還床位數量，該狀態是否存續待定 
  worldtime <- 0 
  NTRY_PER_DAY <- 10   # 一天模擬幾次 
  getday <- function(t) (t - 1) %/% NTRY_PER_DAY + 1 
  # 生成人群資料 
  format_coord <- function(coord, boundary) { 
    if (coord < 0) return(runif(1, 0, 10)) 
    else if  (coord > boundary) return(runif(1, boundary - 10, boundary)) 
    else return(coord) 
  } 
  set.seed(123) 
  people <- tibble( 
    id = 1:CITY_PERSON_SIZE, 
    x = sapply(rnorm(CITY_PERSON_SIZE, CITY_CENTERX, PERSON_DIST_X_SIGMA),  
             format_coord, boundary = CITY_WIDTH),    # (x, y) 為人群點座標 
    y = sapply(rnorm(CITY_PERSON_SIZE, CITY_CENTERY, PERSON_DIST_Y_SIGMA),  
             format_coord, boundary = CITY_HEIGHT), 
    state = STATE_NORMAL,    # 健康狀態 
    infected_time = 0,     # 感染時刻 
    confirmed_time = 0,    # 確診時刻 
    freeze_time = 0,       # 隔離時刻 
    cured_moment = 0,      # 痊癒時刻，為0代表不確定 
    die_moment = 0         # 死亡時刻，為0代表未確定，-1代表不會病死 
  ) %>% 
    mutate(tx = rnorm(CITY_PERSON_SIZE, x, PERSON_DIST_X_SIGMA),  # target x 
           ty = rnorm(CITY_PERSON_SIZE, y, PERSON_DIST_Y_SIGMA), 
           has_target = T, is_arrived = F) 
 
# 隨機選擇初始感染者 
  peop_id <- sample(people$id, ORIGINAL_COUNT) 
  people$state[peop_id] <- STATE_SHADOW 
  people$infected_time[peop_id] <- worldtime 
  people$confirmed_time[peop_id] <- worldtime +  
    max(rnorm(length(peop_id), SHADOW_TIME / 2, SHADOW_TIME_SIGMA), 0) 
 
  ########## 生成床位點 ########## 
  HOSPITAL_X <- 720   # 第一張床位的x座標 
  HOSPITAL_Y <- 80    # 第一張床位的y座標 
  NBED_PER_COLUMN <- 100   # 醫院每一列有多少張床位 
  BED_ROW_SPACE <- 6       # 一行中床位的間距 
  BED_COLUMN_SPACE <- 6    # 一列中床位的間距 
  bed_ncolumn <- ceiling(BED_COUNT / NBED_PER_COLUMN) 
  hosp_beds <- tibble(id = 1, x = 0, y = 0, is_empty = T, state = STATE_NORMAL) %>%  
    slice(-1) 
  if (BED_COUNT > 0) { 
    hosp_beds <- tibble( 
      id = 1:BED_COUNT, 
      x = HOSPITAL_X + rep(((1:bed_ncolumn) - 1) * BED_ROW_SPACE, 
                         each = NBED_PER_COLUMN)[1:BED_COUNT],
      y = HOSPITAL_Y + 10 - BED_COLUMN_SPACE + 
        rep((1:NBED_PER_COLUMN) * BED_COLUMN_SPACE, bed_ncolumn)[1:BED_COUNT],
      is_empty = T,
      person_id = 0       # 佔用床位的患者的序號，床位為空時為0
    )
  }

  ########## 準備畫圖的資料 ##########
  npeople_total <- CITY_PERSON_SIZE
  npeople_shadow <- ORIGINAL_COUNT
  npeople_confirmed <- npeople_freeze <- npeople_cured <- npeople_death <- 0
  nbed_need <- 0

  ########## 畫出初始資料 ##########
  # 設定畫圖引數
  person_color <- data.frame(   # 不同健康狀態的顏色不同
    label = c("健康", "潛伏", "確診", "隔離", "治癒", "死亡"),
    state = c(STATE_NORMAL, STATE_SHADOW, STATE_CONFIRMED, STATE_FREEZE, 
              STATE_CURED, STATE_DEATH),
    color = c(
      "lightgreen",   # 健康
      "#EEEE00",      # 潛伏期
      "red",          # 確診
      "#FFC0CB",      # 隔離
      "green",        # 治癒
      "black"         # 死亡
    ), stringsAsFactors = F
  )
  bed_color <- data.frame(  
    is_empty = c(T, F), color = c("#F8F8FF", "#FFC0CB"), stringsAsFactors = F  
  ) 
  x11(width = 5, height = 7, xpos = 0, ypos = 0, title = "人群變化模擬")
  window_hist <- dev.cur()
  x11(width = 7, height = 7, xpos = 460, ypos = 0, title = "疫情傳播模擬")
  window_scatter <- dev.cur()
  max_plot_x <- ifelse(BED_COUNT > 0, max(hosp_beds$x), CITY_WIDTH) + 10

  # 疫情傳播模擬散點圖
  dev.set(window_scatter)
  plot(x = people$x, y = people$y, cex = .8, pch = 20, xlab = NA, ylab = NA,
       xlim = c(5, max_plot_x), xaxt = "n", yaxt = "n", bty = "n", main = "疫情傳播模擬", 
       sub = paste0("世界時間第 ", getday(worldtime), " 天"),
       col = (people %>% left_join(person_color, by = "state") %>%
              select(color))$color)
  points(x = hosp_beds$x, y = hosp_beds$y, cex = .8, pch = 20,
         col = (hosp_beds %>% left_join(bed_color, by = "is_empty") %>%
              select(color))$color)
  rect(HOSPITAL_X - BED_ROW_SPACE / 2, HOSPITAL_Y + 10 - BED_COLUMN_SPACE, 
       max(hosp_beds$x) + BED_ROW_SPACE / 2, max(hosp_beds$y + BED_COLUMN_SPACE))
  legend(x = 150, y = -30, legend = person_color$label, col = person_color$color,
         pch = 20, horiz = T, bty = "n", xpd = T)
  
  # 人群變化模擬條形圖
  dev.set(window_hist)
  bp_data <- c(npeople_death, npeople_cured, nbed_need, npeople_freeze, 
               npeople_confirmed, npeople_shadow)
  bp_color <- c("black", "green", "#FFE4E1", "#FFC0CB", "red", "#EEEE00")
  bp_labels <- c("死亡", "治癒", "不足n床位", "隔離", "累計n確診", "潛伏")
  bp <- barplot(bp_data, horiz = T, border = NA, col = bp_color, 
                xlim = c(0, CITY_PERSON_SIZE * 1), main = "人群變化模擬", 
                sub = paste0("世界時間第 ", getday(worldtime), " 天"))
  abline(v = BED_COUNT, col = "gray", lty = 3)
  abline(v = CITY_PERSON_SIZE, col = "gray", lty = 1)
  text(x = -350, y = bp, labels = bp_labels, xpd = T)
  text(x = bp_data + CITY_PERSON_SIZE / 15, y = bp, xpd = T,
     labels = ifelse(bp_data > 0, bp_data, ""))
  legend(x = 300, y = -.6, legend = c("總床位數", "城市總人口"), col = "gray",
       lty = c(3, 1), bty = "n", horiz = T, xpd = T)
  Sys.sleep(5)  # 手動調整視窗大小
  
  ########## 更新人群資料 ##########
  # 市民流動意願以及移動位置引數174  MOVE_WISH_SIGMA <- 1
  MOVE_DIST_SIGMA <- 50
  SAFE_DIST <- 2   # 安全距離
  worldtime <- worldtime + 1
  get_min_dist <- function(person, peop) {  # 一個人和一群人之間的最小距離
    min(sqrt((person["x"] - peop$x) ^ 2 + (person["y"] - peop$y) ^ 2))
  }
  for (i in 1:MAX_TRY) {
    # 如果已經隔離或者死亡了，就不需要處理了
    #
    # 處理已經確診的感染者（即患者）
    peop_id <- people$id[people$state == STATE_CONFIRMED & 
                                 people$die_moment == 0]
    if ((npeop <- length(peop_id)) > 0) {
      people$die_moment[peop_id] <- ifelse(
        runif(npeop, 0, 1) < FATALITY_RATE,     # 用均勻分佈模擬確診患者是否會死亡
        people$confirmed_time + max(rnorm(npeop, DIE_TIME, DIE_SIGMA), 0),  # 發病後確定死亡時刻
          -1                                      # 逃過了死神的魔爪
        )
    }
    # 如果患者已經確診，且（世界時刻-確診時刻）大於醫院響應時間，
    # 即醫院準備好病床了，可以擡走了
      peop_id <- people$id[people$state == STATE_CONFIRMED & 
                    worldtime - people$confirmed_time >= HOSPITAL_RECEIVE_TIME]
    if ((npeop <- length(peop_id)) > 0) {
      if ((nbed_empty <- sum(hosp_beds$is_empty)) > 0) {  # 有空餘床位
        nbed_use <- min(npeop, nbed_empty)
        bed_id <- hosp_beds$id[hosp_beds$is_empty][1:nbed_use]
       # 更新患者資訊
        peop_id2 <- sample(peop_id, nbed_use)   # 這裡是隨機選擇，理論上應該按症狀輕重
          people$x[peop_id2] <- hosp_beds$x[bed_id]
        people$y[peop_id2] <- hosp_beds$y[bed_id]
        people$state[peop_id2] <- STATE_FREEZE
        people$freeze_time[peop_id2] <- worldtime
       # 更新床位資訊
        hosp_beds$is_empty[bed_id] <- F
        hosp_beds$person_id[bed_id] <- peop_id2
      } 
    }
    # TODO 需要確定一個變數用於治癒時長。
    # 為了說明問題，暫時用一個正態分佈模擬治癒時長並且假定治癒的人不會再被感染
    peop_id <- people$id[people$state == STATE_FREEZE & 
                           people$cured_moment == 0]
    if ((npeop <- length(peop_id)) > 0) { # 正態分佈模擬治癒時間
      people$cured_moment[peop_id] <- people$freeze_time[peop_id] + 
        max(rnorm(npeop, CURED_TIME, CURED_SIGMA), 0)
    }
    peop_id <- people$id[people$state == STATE_FREEZE & people$cured_moment > 0 &
                           worldtime >= people$cured_moment]
    if ((npeop <- length(peop_id)) > 0) {  # 歸還床位
      people$state[peop_id] <- STATE_CURED
      hosp_beds$is_empty[! hosp_beds$is_empty & hosp_beds$person_id %in% peop_id] <- T
      people$x[peop_id] <- sapply(rnorm(npeop, CITY_CENTERX, PERSON_DIST_X_SIGMA), 
               format_coord, boundary = CITY_WIDTH)    # (x, y) 為人群點座標
      people$y[peop_id] <- sapply(rnorm(npeop, CITY_CENTERY, PERSON_DIST_Y_SIGMA), 
               format_coord, boundary = CITY_HEIGHT)
      people$tx[peop_id] <- rnorm(npeop, people$x[peop_id], PERSON_DIST_X_SIGMA)
      people$ty[peop_id] <- rnorm(npeop, people$y[peop_id], PERSON_DIST_Y_SIGMA)
      people$has_target[peop_id] <- T
      people$is_arrived[peop_id] <- F
    }
    # 處理病死者
    peop_id <- people$id[people$state %in% c(STATE_CONFIRMED, STATE_FREEZE) & 
        worldtime >= people$die_moment & people$die_moment > 0]
    if (length(peop_id) > 0) {  # 歸還床位
      people$state[peop_id] <- STATE_DEATH
      hosp_beds$is_empty[! hosp_beds$is_empty & hosp_beds$person_id %in% peop_id] <- T
    }
    # 處理髮病的潛伏期感染者
    peop_id <- people$id[people$state == STATE_SHADOW &
                        worldtime >= people$confirmed_time]
    if ((npeop <- length(peop_id)) > 0) {
      people$state[peop_id] <- STATE_CONFIRMED   # 潛伏者發病
    }
    # 處理未隔離者的移動問題
    peop_id <- people$id[
      ! people$state %in% c(STATE_FREEZE, STATE_DEATH) & 
      rnorm(CITY_PERSON_SIZE, MOVE_WISH_MU, MOVE_WISH_SIGMA) > 0] # 流動意願
    if ((npeop <- length(peop_id)) > 0) {  # 正態分佈模擬要移動到的目標點
      pp_id <- peop_id[! people$has_target[peop_id] | people$is_arrived[peop_id]]
      if ((npp <- length(pp_id)) > 0) {
        people$tx[pp_id] <- rnorm(npp, people$tx[pp_id], PERSON_DIST_X_SIGMA)
        people$ty[pp_id] <- rnorm(npp, people$ty[pp_id], PERSON_DIST_Y_SIGMA)
        people$has_target[pp_id] <- T
        people$is_arrived[pp_id] <- F
      }
      # 計算運動位移262      dx <- people$tx[peop_id] - people$x[peop_id]
      dy <- people$ty[peop_id] - people$y[peop_id]
      move_dist <- sqrt(dx ^ 2 + dy ^ 2)
      people$is_arrived[peop_id][move_dist < 1] <- T  # 判斷是否到達目標點266    pp_id <- peop_id[move_dist >= 1]
      if ((npp <- length(pp_id)) > 0) {
        udx <- sign(dx[move_dist >= 1])  # x軸運動方向269        udy <- sign(dy[move_dist >= 1])
        # 是否到了邊界
        pid_x <- (1:npp)[people$x[pp_id] + udx < 0 | people$x[pp_id] + udx > CITY_WIDTH]
        pid_y <- (1:npp)[people$y[pp_id] + udy < 0 | people$y[pp_id] + udy > CITY_HEIGHT]
        # 更新到了邊界的點的資訊
        people$x[pp_id[pid_x]] <- people$x[pp_id[pid_x]] - udx[pid_x]
        people$y[pp_id[pid_y]] <- people$y[pp_id[pid_y]] - udy[pid_y]
        people$has_target[unique(c(pp_id[pid_x], pp_id[pid_y]))] <- F
        # 更新沒有到邊界的點的資訊278        people$x[pp_id[! pp_id %in% pid_x]] <- people$x[pp_id[! pp_id %in% pid_x]] + 
          udx[! pp_id %in% pid_x]
        people$y[pp_id[! pp_id %in% pid_y]] <- people$y[pp_id[! pp_id %in% pid_y]] + 
          udy[! pp_id %in% pid_y]
      }
    }
    # 處理健康人被感染的問題
    # 通過一個隨機幸運值和安全距離決定感染其他人286    normal_peop_id <- people$id[people$state == STATE_NORMAL]
    other_peop_id <- people$id[! people$state %in% c(STATE_NORMAL, STATE_CURED)]
    if (length(normal_peop_id) > 0) {
      normal_other_dist <- apply(people[normal_peop_id, ], 1, get_min_dist,
                               peop = people[other_peop_id, ])
      normal2other_id <- normal_peop_id[normal_other_dist < SAFE_DIST &
                          runif(length(normal_peop_id), 0, 1) < BROAD_RATE]
      if ((n2other <- length(normal2other_id)) > 0) {
        people$state[normal2other_id] <- STATE_SHADOW
        people$infected_time[normal2other_id] <- worldtime
        people$confirmed_time[normal2other_id] <- worldtime + 
          max(rnorm(n2other, SHADOW_TIME / 2, SHADOW_TIME_SIGMA), 0)
      }
    }
    # 畫出更新後的資料
      npeople_confirmed <- sum(people$state >= STATE_CONFIRMED)
    npeople_death <- sum(people$state == STATE_DEATH)
    npeople_freeze <- sum(people$state == STATE_FREEZE)
    npeople_shadow <- sum(people$state == STATE_SHADOW)
   npeople_cured <- sum(people$state == STATE_CURED)
    nbed_need <- npeople_confirmed - npeople_cured - npeople_death - BED_COUNT
    nbed_need <- ifelse(nbed_need > 0, nbed_need, 0)  # 不足病床數
    # 疫情傳播模擬散點圖
    dev.set(window_scatter)
    plot(x = people$x, y = people$y, cex = .8, pch = 20, xlab = NA, ylab = NA,
         xlim = c(5, max_plot_x), xaxt = "n", yaxt = "n", bty = "n", main = "疫情傳播模擬", 
         sub = paste0("世界時間第 ", getday(worldtime), " 天"),
         col = (people %>% left_join(person_color, by = "state") %>%
                select(color))$color)
    points(x = hosp_beds$x, y = hosp_beds$y, cex = .8, pch = 20,
           col = (hosp_beds %>% left_join(bed_color, by = "is_empty") %>%
                  select(color))$color)
    rect(HOSPITAL_X - BED_ROW_SPACE / 2, HOSPITAL_Y + 10 - BED_COLUMN_SPACE, 
         max(hosp_beds$x) + BED_ROW_SPACE / 2, max(hosp_beds$y + BED_COLUMN_SPACE))
    legend(x = 150, y = -30, legend = person_color$label, col = person_color$color,
           pch = 20, horiz = T, bty = "n", xpd = T)
    # 人群變化模擬條形圖
    dev.set(window_hist)
    bp_data <- c(npeople_death, npeople_cured, nbed_need, npeople_freeze, 
                 npeople_confirmed, npeople_shadow)
    bp <- barplot(bp_data, horiz = T, border = NA, col = bp_color, 
                  xlim = c(0, CITY_PERSON_SIZE * 1), main = "人群變化模擬", 
                  sub = paste0("世界時間第 ", getday(worldtime), " 天"))
    abline(v = BED_COUNT, col = "gray", lty = 3)
    abline(v = CITY_PERSON_SIZE, col = "gray", lty = 1)
    text(x = -350, y = bp, labels = bp_labels, xpd = T)
    text(x = bp_data + CITY_PERSON_SIZE / 15, y = bp, xpd = T,
         labels = ifelse(bp_data > 0, bp_data, ""))
   legend(x = 300, y = -.6, legend = c("總床位數", "城市總人口"), col = "gray",
           lty = c(3, 1), bty = "n", horiz = T, xpd = T)
  # 更新世界時間
    worldtime <- worldtime + 1
  }