注1：記事として公開していますが、シミュレーションの出来はかなり微妙です。アニメーションがこういう風にできるというヒントになれば幸いです。

注2：かなり処理が重いので、実行の際はご注意ください。

ワシントンポストのシミュレーションに感動したので、Rでマネできないかを試みるその2。

一つ前の記事では、ボールをはねさせるという動画の作成ができたので、後は、

• ボールをたくさんランダムに射出
• あたり判定(あたったら色が変わる）
• 跳ね返る挙動
• 一定時間経過後に色が変わる
• ソーシャルディスタンシング中(動かない)ボールを設定する

あたりの処理ができればそれなりのアニメーションはできるはず。

ボールの位置や動いている方向、ボールの状態等は一つのTibbleで管理してみる方針でやります。

ボールの状態を生成する

canvas_x <- 1200 #ボールが動ける範囲を指定
canvas_y <- 900  #ボールが動ける範囲を指定

#ボールの状態を作成する関数
generate_ball_status <- function(ball_numbers = 10, percent_social_distance = 0, xsize, ysize){

  tibble(
    id = 1:ball_numbers,
    xpos = runif(ball_numbers, 0, xsize),
    ypos = runif(ball_numbers, 0, ysize),
    direction = runif(ball_numbers, 0, 2*pi), #ラジアンで現在進んでいる方向を表す。
    velocity  = sample(x = c(0,1), 
                       size = ball_numbers, 
                       replace = TRUE, 
                       prob = c(percent_social_distance, 1 - percent_social_distance)),
    status = c("I",rep("N", ball_numbers - 1)) 
  ) %>% 

    mutate(status = factor(status, levels = c("I","N","R"))) %>% 

    mutate(time_after_infected = 0) %>% 

    mutate(
      xvelocity = velocity * cos(direction),
      yvelocity = velocity * sin(direction)

      ) %>% 
    select(-direction, -velocity)
}

dat <- generate_ball_status(10, percent_social_distance = 0, xsize = canvas_x, ysize = canvas_y)

状態更新を行う関数の作成

(尚、ワシントンポストで表現されていた、都市閉鎖・隔離のための中の壁は今回作りません。)

update_ball_status <- function(dat, timespan = 10, xsize, ysize){

  #timespan時間経過後の、各々のボールの位置を計算。
  dat2 <- dat %>% 
    mutate(
     xpos =  xpos + (xvelocity * timespan),
     ypos =  ypos + (yvelocity * timespan)
    )

  #範囲外に出ているボールがあれば、位置をはみ出した分だけ
  #内側に戻して、速度を逆方向にする。
  dat3 <- dat2 %>% 
    mutate(
      new_xpos = case_when(
        xpos < 0     ~ -1 * xpos,
        xpos > xsize ~ xpos - 2*(xpos - xsize),
        TRUE         ~ xpos
      )
    ) %>% 
    mutate(
      xvelocity = case_when(
        xpos < 0 | xpos > xsize ~ xvelocity * -1,
        TRUE                      ~ xvelocity
      )
    ) %>% 
    mutate(
      new_ypos = case_when(
        ypos < 0     ~ -1 * ypos,
        ypos > ysize ~ ysize - 2 * (ypos-ysize),
        TRUE         ~ ypos
      )
    ) %>% 
    mutate(
      yvelocity = case_when(
        ypos < 0 | ypos > ysize ~ yvelocity * -1,
        TRUE                      ~ yvelocity
      )
    ) %>% 
    mutate(
      ypos = new_ypos,
      xpos = new_xpos
    ) %>% 
    select(-new_ypos, -new_xpos)

  #互いに一定範囲内に入っていれば反発する。
  #かつ、どちらかが感染者であればステータスを更新する。

  position <- dat3 %>% 
    select(id, xpos, ypos, status, xvelocity, yvelocity)

  dat3 <- dat3 %>% 
    mutate(position = list(position)) %>% 
    mutate(new_state = pmap(.l = list(id, xpos, ypos, status, position), .f = ~{

      adat <- ..5

      #ここであたり判定を単純な距離でもとめて、
      #それが一定以下の場合は、相手のxy方向のベクトルの値に応じて、新しいベクトルを設定する。
      adat <- adat %>% 
        filter(id != ..1) %>% 
        mutate(is_collide =  sqrt((..2 - xpos)^2 + (..3 - ypos)^2) <= 15 ) %>% 
        filter(is_collide)

      if(nrow(adat)==0){
        return(
          tibble(
            xvel_mult = 1, 
            yvel_mult = 1,
            any_infected = FALSE
        ))
      }else{

        return(tibble(
          xvel_mult = if_else(sum(adat$xvelocity) > 0, -1, 1),
          yvel_mult = if_else(sum(adat$yvelocity) > 0, -1, 1),
          any_infected = any(adat$status=="I")
        ))
      }

    }))

  #最終的な反発後の速度を修正。
  dat4 <- dat3 %>% 
    select(-position) %>% 
    unnest(new_state) %>% 
    mutate(
      xvelocity = xvelocity * xvel_mult,
      yvelocity = yvelocity * yvel_mult
    ) %>% 
    mutate(
      status = case_when(
        status == "R"                ~ "R",
        status == "N" & any_infected ~ "I",
        TRUE                         ~ as.character(status)
      )
    ) %>% 
    select(-xvel_mult, -yvel_mult, -any_infected)

  #I(Infection)状態のボールの感染時間を足して、一定時間経過後(timespanの100倍後) に、回復させる。
  dat5 <- dat4 %>% 
    mutate(time_after_infected = if_else(status == "I",
                                         time_after_infected + timespan,
                                         time_after_infected),
           status = if_else(time_after_infected >= timespan*100,
                            "R",
                            status),
           status = factor(status, levels = c("I","N","R")) )

  return(dat5)
}

#ここまでの関数でシミュレーションを行って結果を保存する関数を作成する。
simulate_covid19 <- function(number_of_ball = 200, 
                             percent_social_distance = 0,
                             output_file_dir,
                             output_mp4_name, 
                             max_repeat = 200){

  dat <- generate_ball_status(number_of_ball, 
                            percent_social_distance = percent_social_distance, 
                            xsize = canvas_x*2, 
                            ysize = canvas_y*2)

  dir.create(output_file_dir)
  file.remove(list.files(output_file_dir, full.names = TRUE))

  baseg <- ggplot() + 
    geom_point(aes(x = xpos, y = ypos, color = status)) +
    scale_x_continuous(breaks = NULL, labels = NULL) +
    scale_y_continuous(breaks = NULL, labels = NULL) +
    labs(x = NULL, y = NULL) +
    scale_color_discrete(drop = FALSE) +
    coord_cartesian(xlim = c(0, canvas_x), ylim = c(0,canvas_y))

  res <- list()
  stop_loop <- TRUE
  i <- 1
  while(stop_loop){
    current_infection <- sum(dat$status=="I")
    current_recovery  <- sum(dat$status=="R")
    current_normal    <- nrow(dat)- current_infection - current_recovery
    print(str_glue("{i}:I/R/N->{current_infection}/{current_recovery}/{current_normal}"))

    dat <- dat %>% update_ball_status(.,20,canvas_x, canvas_y)
    res[[i]] <- dat

    if(i == max_repeat){stop_loop <- FALSE}
    if(sum(dat$status == "N")/nrow(dat) == 0){stop_loop <- FALSE}
    if(sum(dat$status == "I") == 0){stop_loop <- FALSE}
    i <- i + 1
  }

  restib <- tibble(res) %>% 
    mutate(graph = map(res, ~{baseg %+% .}))

  map(1:(i-1), ~{
    print(str_glue("{.}/{i}"))
    ggsave(restib$graph[[.]], 
           filename = sprintf("%s/sim1%04d.png",output_file_dir , .), 
           width = 120.1, height = 90, units = "mm")
  })

  input_files <- list.files(output_file_dir, full.names = TRUE, pattern = "png$")
  av::av_encode_video(input = input_files, 
                  output = str_c(output_file_dir,"/",output_mp4_name,".mp4"))

  write_rds(x = restib, path = str_c(output_file_dir,"/simulated.rds"), compress = "gz")
}

# restib
# i <- nrow(restib)
# output_file_dir <- "simulated_data/sim200_00"
# output_mp4_name <- "sim200_00"

ということで,simulation!(以下の関数は実行すると結構時間がかかるのでコメントアウトしてあります。あと、シミュレーション結果の保存場所等は、適宜環境に応じて変更する必要があるかもしれません。)

simulate_covid19(200, 0  , "simulated_data/sim200_00", "sim200_00", 500)

simulate_covid19(200, 0.2, "simulated_data/sim200_20", "sim200_20", 500)

simulate_covid19(200, 0.6, "simulated_data/sim200_70", "sim200_70", 500)

simulate_covid19(200, 0.7, "simulated_data/sim200_90", "sim200_90", 500)

simulate_covid19(200, 0.8, "simulated_data/sim200_90", "sim200_90", 500)

simulate_covid19(200, 0.9, "simulated_data/sim200_90", "sim200_90", 500)

静止ボール0%の場合の動画：

うーん、左下のコーナーにボールが寄ってしまっていたり、
静止しているボールとの衝突が衝突として表現できていなかったり、
改善するべきところはたくさんありますが、
やりたいことがおおむね実行はできたので満足です。

最後に、シミュレーションの結果として、時間経過における、ボールの状態の個数をそれぞれの静止ボールの割合毎に描画をして終わりにしましょう。

obtain_summarised_data <- function(sd = "00"){
  getwd()
  dat <- read_rds(str_c(getwd(),"/simulated_data/sim200_",sd,"/simulated.rds"))

  return(dat %>% 
    select(res) %>% 
    mutate(timelapse = 1:n()) %>% 
    mutate(num_i = map_int(res,~{sum(.$status=="I")})) %>% 
    mutate(num_n = map_int(res,~{sum(.$status=="N")})) %>% 
    mutate(num_r = map_int(res,~{sum(.$status=="R")})) %>% 
    select(-res) %>% 
    mutate(social_distance = sd)
  )
}

gdat <- tibble(
  tgt = c("00","20","60","70","80","90")
) %>% 
  mutate(simulated = map(tgt, ~{obtain_summarised_data(.)}))

gdat <- gdat %>% 
  unnest(simulated)

gdat %>% 
  pivot_longer(cols = c(num_i, num_n, num_r)) %>% 
ggplot() +
  geom_col(aes(x = timelapse, y = value, fill = name )) + 
  facet_wrap(~social_distance)

＊WordPressへの移行に伴い表示していません。

以上です。

感想：

ggplotでもアニメーションを作成できました。
とはいっても、基本的には「紙芝居」をまとめて動画ファイルにするだけなのですが、もしアニメーションをRで作成するという珍しい状況がありましたら、何かの参考になりますと幸いです。