changes inline to make web-version possible

This commit is contained in:
Swapnil Mishra 2020-04-10 14:06:54 +01:00
Родитель e70be3b78b
Коммит a1e97e54e5
9 изменённых файлов: 686 добавлений и 31 удалений

2
.gitignore поставляемый
Просмотреть файл

@ -46,6 +46,8 @@ results/*.png
figures/*.png
figures/*.pdf
Rplots.pdf
web/*
results/*.csv

9
base.r
Просмотреть файл

@ -223,11 +223,16 @@ system(paste0("Rscript covariate-size-effects.r ", filename,'-stanfit.Rdata'))
mu = (as.matrix(out$mu))
colnames(mu) = countries
g = (mcmc_intervals(mu,prob = .9))
ggsave(sprintf("results/%s-mu.pdf",filename),g,width=4,height=6)
ggsave(sprintf("results/%s-mu.png",filename),g,width=4,height=6)
tmp = lapply(1:length(countries), function(i) (out$Rt_adj[,stan_data$N[i],i]))
Rt_adj = do.call(cbind,tmp)
colnames(Rt_adj) = countries
g = (mcmc_intervals(Rt_adj,prob = .9))
ggsave(sprintf("results/%s-final-rt.pdf",filename),g,width=4,height=6)
ggsave(sprintf("results/%s-final-rt.png",filename),g,width=4,height=6)
system(paste0("Rscript plot-3-panel.r ", filename,'-stanfit.Rdata'))
system(paste0("Rscript plot-forecast.r ",filename,'-stanfit.Rdata'))
system(paste0("Rscript make-table.r results/",filename,'-stanfit.Rdata'))
verify_result <- system(paste0("Rscript web-verify-output.r ", filename,'.Rdata'),intern=FALSE)
if(verify_result != 0){
stop("Verification of web output failed!")
}

Просмотреть файл

@ -3,6 +3,8 @@ library(scales)
library(ggpubr)
library(bayesplot)
library(matrixStats)
library(cowplot)
library(svglite)
args <- commandArgs(trailingOnly = TRUE)
filename <- args[1]
@ -72,5 +74,11 @@ p = ggplot(data) +theme_pubr() + geom_point(aes(x=m,y=parameter,colour=type),po
theme(plot.margin = margin(0, 2, 0, .5, "cm"))
#+ guides(fill=guide_legend(nrow=2))
p
ggsave(filename = "results/covars-alpha-reduction.pdf",
p,height=4,width=8)
ggsave(filename = "results/covars-alpha-reduction.png",
p,height=4,width=8)
dir.create("web/figures/desktop/", showWarnings = FALSE, recursive = TRUE)
dir.create("web/figures/mobile/", showWarnings = FALSE, recursive = TRUE)
save_plot(filename = paste0("web/figures/desktop/", "covars-alpha-reduction.svg"),
p, base_height = 4, base_asp = 1.618 * 2 * 8/12)
save_plot(filename = paste0("web/figures/mobile/", "covars-alpha-reduction.svg"),
p, base_height = 4, base_asp = 1.1)

Разница между файлами не показана из-за своего большого размера Загрузить разницу

Двоичные данные
data/COVID-19-up-to-date.rds

Двоичный файл не отображается.

133
make-table.r Normal file
Просмотреть файл

@ -0,0 +1,133 @@
library(matrixStats)
args <- commandArgs(trailingOnly = TRUE)
filename <- args[1]
print(sprintf("Running %s",filename))
load(filename)
df_pop= read.csv("data/popt_ifr.csv", stringsAsFactors = FALSE)
df_pop$country[df_pop$country == "United Kingdom"] = "United_Kingdom"
dates_italy <- dates[[which(countries == "Italy")]]
len_dates <- length(dates_italy)
date_till_percentage = date_iso <- as.character(Sys.Date())
cases <- vector("list", length = length(countries))
total_cases <- vector("list", length = length(countries))
total_cases_ui <- vector("list", length = length(countries))
total_cases_li <- vector("list", length = length(countries))
deaths <- vector("list", length = length(countries))
total_deaths <- vector("list", length = length(countries))
rt <- vector("list", length = length(countries))
fraction_infected <- vector("list", length = length(countries))
fraction_infected_li <- vector("list", length = length(countries))
fraction_infected_ui <- vector("list", length = length(countries))
fraction_obs_infected <- vector("list", length = length(countries))
fraction_total_obs_infected <- vector("list", length = length(countries))
y <- vector("list", length = length(countries))
for(i in 1:length(countries)) {
Country = countries[i]
x = dates[[i]]
N = length(x)
forecast = 7
x = c(x,x[length(x)]+1:forecast)
padding <- len_dates - length(dates[[i]])
y[[i]] = c(rep(0, padding),reported_cases[[i]], rep(NA, forecast))
cases[[i]] = c(rep(0, padding), round(colMeans(prediction[,1:length(x),i])))
total_cases[[i]] = c( round(cumsum(colMeans(prediction[,1:length(x),i]))))
# chk = c(round((colMeans(rowCumsums(prediction[,1:length(x),i])))))
total_cases_li[[i]] = c(
round((colQuantiles(rowCumsums(prediction[,1:length(x),i]),probs=.025))))
total_cases_ui[[i]] = c(
round((colQuantiles(rowCumsums(prediction[,1:length(x),i]),probs=.975))))
deaths[[i]] = c(rep(0, padding), round(colMeans(estimated.deaths[,1:length(x),i])))
total_deaths[[i]] = c(rep(0, padding), round(cumsum(colMeans(estimated.deaths[,1:length(x),i]))))
rt[[i]] = c(rep(NA, padding), colMeans(out$Rt[,1:length(x),i]))
fraction_infected[[i]] = c(rep(0, padding), total_cases[[i]]/ df_pop[df_pop$country==Country,]$popt)
fraction_infected_li[[i]] = c(rep(0, padding),
total_cases_li[[i]]/ df_pop[df_pop$country==Country,]$popt)
fraction_infected_ui[[i]] = c(rep(0, padding),
total_cases_ui[[i]]/ df_pop[df_pop$country==Country,]$popt)
fraction_obs_infected[[i]] = c(rep(0, padding), y[[i]] / cases[[i]])
fraction_total_obs_infected[[i]] = c(rep(0, padding), cumsum(y[[i]]) / cases[[i]])
total_cases[[i]] = c(rep(0, padding),total_cases[[i]])
}
dates_italy = c(dates_italy,dates_italy[length(dates_italy)]+1:forecast)
cases <- do.call(rbind, cases)
cases_df <- as.data.frame(cases)
names(cases_df) <- dates_italy
cases_df$countries <- countries
# write.csv(cases_df, "figures/cases.csv")
total_cases <- do.call(rbind, total_cases)
total_cases_df <- as.data.frame(total_cases)
names(total_cases_df) <- dates_italy
total_cases_df$countries <- countries
# write.csv(total_cases_df, "figures/total_cases.csv")
deaths <- do.call(rbind, deaths)
deaths_df <- as.data.frame(deaths)
names(deaths_df) <- dates_italy
deaths_df$countries <- countries
# write.csv(deaths_df, "figures/deaths.csv")
total_deaths <- do.call(rbind, total_deaths)
total_deaths_df <- as.data.frame(total_deaths)
names(total_deaths_df) <- dates_italy
total_deaths_df$countries <- countries
# write.csv(total_deaths_df, "figures/total_deaths.csv")
rt <- do.call(rbind, rt)
rt_df <- as.data.frame(rt)
names(rt_df) <- dates_italy
rt_df$countries <- countries
# write.csv(rt_df, "figures/rt.csv")
fraction_infected <- do.call(rbind, fraction_infected)
fraction_infected_df <- as.data.frame(fraction_infected)
names(fraction_infected_df) <- dates_italy
fraction_infected_df$countries <- countries
# write.csv(fraction_infected_df, "figures/fraction_infected.csv")
fraction_infected_li <- do.call(rbind, fraction_infected_li)
fraction_infected_li_df <- as.data.frame(fraction_infected_li)
names(fraction_infected_li_df) <- dates_italy
fraction_infected_li_df$countries <- countries
# write.csv(fraction_infected_li_df, "figures/fraction_infected_li.csv")
fraction_infected_ui <- do.call(rbind, fraction_infected_ui)
fraction_infected_ui_df <- as.data.frame(fraction_infected_ui)
names(fraction_infected_ui_df) <- dates_italy
fraction_infected_ui_df$countries <- countries
# write.csv(fraction_infected_ui_df, "figures/fraction_infected_ui.csv")
total_infected = data.frame(countries=countries,mean=fraction_infected[,dates_italy == date_till_percentage],
li=fraction_infected_li[,dates_italy == date_till_percentage],ui=fraction_infected_ui[,dates_italy == date_till_percentage])
total_infected$value = sprintf("%.02f%% [%.02f%%-%.02f%%]",
total_infected$mean*100,total_infected$li*100,total_infected$ui*100)
total_infected[order(total_infected$countries),c("countries","value")]
total_infected <- total_infected[,c("countries","value")]
write.csv(total_infected,paste0("results/total_infected_",date_till_percentage,".csv"),row.names=F)
fraction_obs_infected <- do.call(rbind, fraction_obs_infected)
fraction_obs_infected_df <- as.data.frame(fraction_obs_infected)
names(fraction_obs_infected_df) <- dates_italy
fraction_obs_infected_df$countries <- countries
# write.csv(fraction_obs_infected_df, "figures/fraction_obs_infected.csv")
fraction_total_obs_infected <- do.call(rbind, fraction_total_obs_infected)
fraction_total_obs_infected_df <- as.data.frame(fraction_total_obs_infected)
names(fraction_total_obs_infected_df) <- dates_italy
fraction_total_obs_infected_df$countries <- countries
# write.csv(fraction_total_obs_infected_df, "figures/fraction_total_obs_infected.csv")

Просмотреть файл

@ -12,6 +12,7 @@ library(gridExtra)
library(ggpubr)
library(bayesplot)
library(cowplot)
library(svglite)
source("utils/geom-stepribbon.r")
#---------------------------------------------------------------------------
@ -19,11 +20,19 @@ make_three_pannel_plot <- function(){
args <- commandArgs(trailingOnly = TRUE)
filename2 <- args[1]
if (length(args)==1){
filename2 = args[1]
percent_pop = FALSE
} else {
filename2 = args[1]
percent_pop = args[2]
}
load(paste0("results/", filename2))
print(sprintf("loading: %s",paste0("results/",filename2)))
covariates = read.csv('data/interventions.csv', stringsAsFactors = FALSE)
names_covariates = c('Schools + Universities','Self-isolating if ill', 'Public events', 'Lockdown', 'Social distancing encouraged')
names_covariates = c('Schools + Universities','Self-isolating if ill', 'Public events',
'Lockdown', 'Social distancing encouraged')
covariates <- covariates %>%
filter((Type %in% names_covariates))
covariates <- covariates[,c(1,2,4)]
@ -113,15 +122,18 @@ make_three_pannel_plot <- function(){
make_plots(data_country = data_country,
covariates_country_long = covariates_country_long,
filename2 = filename2,
country = country)
country = country,
percent_pop = percent_pop)
}
}
#---------------------------------------------------------------------------
make_plots <- function(data_country, covariates_country_long,
filename2, country){
filename2, country, percent_pop){
if (country == 'United_Kingdom')
country = 'United Kingdom'
data_cases_95 <- data.frame(data_country$time, data_country$predicted_min,
data_country$predicted_max)
names(data_cases_95) <- c("time", "cases_min", "cases_max")
@ -139,14 +151,15 @@ make_plots <- function(data_country, covariates_country_long,
geom_ribbon(data = data_cases,
aes(x = time, ymin = cases_min, ymax = cases_max, fill = key)) +
xlab("") +
ylab("Daily number of infections") +
ylab("Daily number of infections\n") +
scale_x_date(date_breaks = "weeks", labels = date_format("%e %b")) +
scale_y_continuous(expand = c(0, 0), labels = comma) +
scale_fill_manual(name = "", labels = c("50%", "95%"),
values = c(alpha("deepskyblue4", 0.55),
alpha("deepskyblue4", 0.45))) +
theme_pubr() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "None") +
legend.position = "None") + ggtitle(country) +
guides(fill=guide_legend(ncol=1))
data_deaths_95 <- data.frame(data_country$time, data_country$death_min,
@ -168,9 +181,12 @@ make_plots <- function(data_country, covariates_country_long,
data = data_deaths,
aes(ymin = death_min, ymax = death_max, fill = key)) +
scale_x_date(date_breaks = "weeks", labels = date_format("%e %b")) +
scale_y_continuous(expand = c(0, 0), labels = comma) +
scale_fill_manual(name = "", labels = c("50%", "95%"),
values = c(alpha("deepskyblue4", 0.55),
alpha("deepskyblue4", 0.45))) +
ylab("Daily number of deaths\n") +
xlab("") +
theme_pubr() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "None") +
@ -218,17 +234,38 @@ make_plots <- function(data_country, covariates_country_long,
scale_x_date(date_breaks = "weeks", labels = date_format("%e %b"),
limits = c(data_country$time[1],
data_country$time[length(data_country$time)])) +
scale_y_continuous(expand = expansion(mult=c(0,0.1))) +
theme_pubr() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(legend.position="right")
if (country == 'United Kingdom')
country = 'United_Kingdom'
# Special plot settings for mobile
p3_mobile <- p3 +
theme(legend.position="below")
# Plots for Web, Desktop version
dir.create("web/figures/desktop/", showWarnings = FALSE, recursive = TRUE)
save_plot(filename = paste0("web/figures/desktop/", country, "_infections", ".svg"),
p1, base_height = 4, base_asp = 1.618)
save_plot(filename = paste0("web/figures/desktop/", country, "_deaths", ".svg"),
p2, base_height = 4, base_asp = 1.618)
save_plot(filename = paste0("web/figures/desktop/", country, "_rt", ".svg"),
p3, base_height = 4, base_asp = 1.618 * 2)
# Plots for Web, Mobile version
dir.create("web/figures/mobile/", showWarnings = FALSE, recursive = TRUE)
save_plot(filename = paste0("web/figures/mobile/", country, "_infections", ".svg"),
p1, base_height = 4, base_asp = 1.1)
save_plot(filename = paste0("web/figures/mobile/", country, "_deaths", ".svg"),
p2, base_height = 4, base_asp = 1.1)
save_plot(filename = paste0("web/figures/mobile/", country, "_rt", ".svg"),
p3_mobile, base_height = 4, base_asp = 1.1)
p <- plot_grid(p1, p2, p3, ncol = 3, rel_widths = c(1, 1, 2))
save_plot(filename = paste0("figures/", country, "_three_pannel_", filename2, ".pdf"),
save_plot(filename = paste0("figures/", country, "_three_pannel_", filename2, ".png"),
p, base_width = 14)
}
#-----------------------------------------------------------------------------------------------
#filename <- "base-joint-1236305.pbs.Rdata"
# make_three_pannel_plot(filename)
make_three_pannel_plot()
make_three_pannel_plot()

Просмотреть файл

@ -132,8 +132,16 @@ make_single_plot <- function(data_country, data_country_forecast, filename, coun
color="black")
print(p)
ggsave(file= paste0("figures/", country, "_forecast_", filename, ".pdf"),
ggsave(file= paste0("figures/", country, "_forecast_", filename, ".png"),
p, width = 10)
# Produce plots for Website
dir.create("web/figures/desktop/", showWarnings = FALSE, recursive = TRUE)
save_plot(filename = paste0("web/figures/desktop/", country, "_forecast", ".svg"),
p, base_height = 4, base_asp = 1.618 * 2 * 8/12)
dir.create("web/figures/mobile/", showWarnings = FALSE, recursive = TRUE)
save_plot(filename = paste0("web/figures/mobile/", country, "_forecast", ".svg"),
p, base_height = 4, base_asp = 1.1)
}
#-----------------------------------------------------------------------------------------------
make_forecast_plot()

53
web-verify-output.r Normal file
Просмотреть файл

@ -0,0 +1,53 @@
library(jsonlite)
d=readRDS('data/COVID-19-up-to-date.rds')
countries <- c(
"Denmark",
"Italy",
"Germany",
"Spain",
"United_Kingdom",
"France",
"Norway",
"Belgium",
"Austria",
"Sweden",
"Switzerland",
"Greece",
"Portugal",
"Netherlands"
)
verify_web_output <- function(){
plot_names <- c("deaths", "forecast", "infections", "rt")
plot_versions <- c("mobile", "desktop")
args <- commandArgs(trailingOnly = TRUE)
filename2 <- args[1]
load(paste0("results/", filename2))
print(sprintf("loading: %s",paste0("results/",filename2)))
date_results <- list()
for(country in countries) {
for (plot_version in plot_versions) {
for (plot_name in plot_names) {
path = sprintf("web/figures/%s/%s_%s.svg", plot_version, country, plot_name)
if (! file.exists(path)) {
stop(sprintf("Missing web output during verification: %s", path))
}
}
}
d1=d[d$Countries.and.territories==country,]
d1$date = as.Date(d1$DateRep,format='%d/%m/%Y')
latest_date = max(d1$date)
date_results[[country]] = latest_date
}
dir.create("web/data/", showWarnings = FALSE, recursive = TRUE)
write_json(date_results, "web/data/latest-updates.json", auto_unbox=TRUE)
}
verify_web_output()