#### Loading Package #### ## Installation from CRAN # install.packages('hmer') ## For newest version, install from github # remotes::install_github("andy-iskauskas/hmer") library("hmer") ### Load in the required and suggested packages # deSolve is only required to run ode_results library("lhs") library("deSolve") library("ggplot2") set.seed(12) ################# HELPER FUNCTIONS ################### # Simple code for generating the output of the SIR deterministic model ode_results <- function(parms, end_time = 25) { des = function(time, state, parms) { with(as.list(c(state, parms)), { dS <- aSR*R-aSI*I*S/(S+I+R) dI <- aSI*I*S/(S+I+R)-aIR*I dR <- aIR*I-aSR*R return(list(c(dS, dI, dR))) }) } yini = c(S = 950, I = 50, R = 0) times = seq(0, end_time, by = 1) out = ode(yini, times, des, parms) return(out) } # Collates the results to give the outputs for emulation get_res <- function(inputs) { ode_out <- data.frame(t(apply(inputs, 1, function(x) ode_results(x)[11, c('S', 'I', 'R')]))) return(cbind(inputs, ode_out) |> setNames(c(names(inputs), c('nS', 'nI', 'nR')))) } # Very simple wrapper for lhs::randomLHS that scales the inputs to our ranges get_lhs_design <- function(npoints, ranges) { initial_design <- lhs::maximinLHS(npoints, length(ranges)) return(setNames(data.frame( t(apply(initial_design, 1, function(x) x*sapply(ranges, diff)+sapply(ranges, "[[", 1)))), names(ranges))) } ## Usage of the above functions ## # initial_points <- get_lhs_design(90, ranges) # initial_results <- get_res(initial_points) ########### FIGURE 1 ############## func <- function(x) { 2 * x + 3 * x * sin(5 * pi * (x - 0.1)/0.4) } data1d <- data.frame(x = seq(0.05, 0.5, by = 0.05), f = func(seq(0.05, 0.5, by = 0.05))) ranges1d <- list(x = c(0, 0.6)) em1d <- emulator_from_data(data1d, c('f'), ranges1d) test1d <- data.frame(x = seq(0, 0.6, by = 0.001)) em1d_exp <- em1d$f$get_exp(test1d) em1d_var <- em1d$f$get_cov(test1d) plot1d <- data.frame( x = test1d$x, f = func(test1d$x), E = em1d_exp, max = em1d_exp + 3*sqrt(abs(em1d_var)), min = em1d_exp - 3*sqrt(abs(em1d_var)) ) plot(data = plot1d, f ~ x, ylim = c(min(plot1d[,-1]), max(plot1d[-1])), type = 'l', main = "Emulation of a Simple 1d Function", xlab = "x", ylab = "f(x)") lines(data = plot1d, E ~ x, col = 'blue') lines(data = plot1d, max ~ x, col = 'red', lty = 2) lines(data = plot1d, min ~ x, col = 'red', lty = 2) points(data = data1d, f ~ x, pch = 16, cex = 1) legend('bottomleft', inset = c(0.05, 0.05), legend = c("Function value", "Emulated value", "Uncertainty Bounds"), col = c('black', 'blue', 'red'), lty = c(1,1,2)) ########### FIGURE 2 ############## # Illustrative schematic for using point generation using Proto_emulator # The exact figure in the manuscript can be replicated with v109i10-cardioid.R c_ranges <- list(x = c(-1, 1), y = c(0, 1.5)) c_membership <- function(x) { a <- 1/2 r2 <- apply(x, 1, function(x) x[[1]]^2+(x[[2]]-1)^2) theta <- apply(x, 1, function(x) atan2(x[[2]]-1, x[[1]])) card_membership <- r2 <= a*(1-sin(theta)) c_rad <- 0.1 r2 <- apply(x, 1, function(x) (x[[1]]+0.1)^2+(x[[2]]-0.6)^2) circ_membership <- r2 <= c_rad return(card_membership & !circ_membership) } c_em <- Proto_emulator$new( ranges = c_ranges, output_name = "Z", predict_func = function(x) return(x), variance_func = function(x) return(x) ) cust_imp <- function(em, x, z, cutoff, ...) { return(c_membership(x)) } op <- par(mfrow = c(2, 2)) lhs_points <- generate_new_design(c_em, 1000, NULL, opts = list(accept_measure = cust_imp)) plot(lhs_points, pch = 16, cex = 0.5) line_points <- generate_new_design(c_em, 1000, NULL, method = 'line', plausible_set = lhs_points, opts = list(accept_measure = cust_imp)) plot(lhs_points, pch = 16, cex = 0.3, col = 'grey') points(line_points[!row.names(line_points) %in% row.names(lhs_points),], pch = 16, cex = 0.5) importance_points <- generate_new_design(c_em, 1000, NULL, method = "importance", plausible_set = rbind(lhs_points, line_points), opts = list(accept_measure = cust_imp)) plot(importance_points, pch = 16, cex = 0.5) plot(rbind.data.frame(lhs_points, line_points), pch = 16, cex = 0.3, col = 'grey') points(importance_points[!row.names(importance_points) %in% row.names(rbind.data.frame(lhs_points, line_points)),], pch = 16, cex = 0.5) par(op) ########### Section 3 ########### example_points <- data.frame(rate2 = c(1, -1), rate1 = c(1, 2)) points <- matrix(c(1, 1, 2, -1), nrow = 2, byrow = TRUE) ########### EMULATION SET UP ############## # Define the training and validation runs training <- SIRSample$training valid <- SIRSample$validation # Set up the ranges ranges <- list(aSI = c(0.1, 0.8), aIR = c(0, 0.5), aSR = c(0, 0.05)) targets <- list(nS = c(580, 651), nI = list(val = 169, sigma = 8.45), nR = c(199, 221)) ######### FIGURE 3 ######### plot(do.call('rbind.data.frame', SIRSample)[,names(ranges)], col = rep(c('black', 'blue'), times = c(nrow(SIRSample$training), nrow(SIRSample$validation))), pch = c(rep(16, nrow(training)), rep(4, nrow(valid)))) # Train the emulators using the emulator_from_data defaults ems_wave1 <- emulator_from_data(training, names(targets), ranges) # Examine the emulator structure ems_wave1$nI # Demonstrating the behaviour of emulators at training points ems_wave1$nS$get_exp(training[c(7, 19), ]) SIRSample$training[c(7, 19), 'nS'] ems_wave1$nS$get_cov(training[c(7, 19), ]) # Diagnostics on the emulators ###### FIGURE 4 ####### validation_diagnostics(ems_wave1, targets, SIRSample$validation) ###### FIGURE 5a ####### validation_diagnostics(ems_wave1, validation = valid) # Zooming in on the lower part of nR ###### FIGURE 5b ####### validation_diagnostics(ems_wave1$nR, validation = valid[SIRSample$valid$nR < 150,], which_diag = "cd", row = 1) ########### VISUALISATION ########### # Plot of emulator expectations ####### FIGURE 6 ####### emulator_plot(ems_wave1, include_legend = TRUE) ems_wave1$nR$plot(plot_type = 'sd', params = c('aSI', 'aIR'), fixed_vals = c(aSR = 0.045)) + geom_point(data = training, aes(x = aSI, y = aIR)) # Demonstration of the difference between an untrained and trained emulator ###### FIGURE 7a ######## plot(ems_wave1$nR) ###### FIGURE 7b ######## plot(ems_wave1$nR$o_em) # nth implausibility plot ###### FIGURE 8a ####### emulator_plot(ems_wave1$nS, plot_type = 'imp', targets = targets, ppd = 40) ###### FIGURE 8b ####### emulator_plot(ems_wave1, plot_type = "nimp", targets = targets, ppd = 40, cb = TRUE) # Plot lattice: minimum implausibility and optical depth ####### FIGURE 9 ###### plot_lattice(ems_wave1, targets, ppd = 35) # Diagnostic for 'optimal' cutoff ###### FIGURE 10 ####### space_removed(ems_wave1, targets, ppd = 20) + geom_vline(xintercept = 3, lty = 2) # Same idea, but w.r.t. emulator variance ### NOT USED IN PAPER # space_removed(ems_wave1, targets, modified = 'var', ppd = 20, # u_mod = seq(0.5, 1.5, by = 0.1)) ############# PROPOSING POINTS #################### # Simple generate_new_runs call proposal <- generate_new_design(ems_wave1, 500, targets) # Plotting the proposal ######## FIGURE 11 ######## plot_wrap(proposal, ranges) ############# MULTIWAVE PLOTTING ################ data("SIRMultiWaveData", package = "hmer") # Plot the outputs across waves ####### FIGURE 12 ######### wave_values(SIRMultiWaveData, targets, l_wid = 0.8) # Examining the emulator uncertainties sapply(SIRMultiWaveEmulators[[3]], function(x) x$u_sigma) # Simulator plot: checking that we can hit all targets with some model runs ######## FIGURE 13 ######## simulator_plot(SIRMultiWaveData, targets, barcol = 'black') ############# STOCHASTIC EXAMPLE #################### # Setting up the parameters for the Gillespie algorithm gillespied <- function (N, T = 400, dt = 1, ...) { tt <- 0 n <- T%/%dt x <- N$M S <- t(N$Post - N$Pre) u <- nrow(S) v <- ncol(S) xmat <- matrix(0, ncol = u, nrow = n) i <- 1 target <- 0 repeat { h <- N$h(x, tt, ...) h0 <- sum(h) if (h0 < 1e-10) tt <- 1e99 else tt <- tt + rexp(1, h0) while (tt >= target) { xmat[i, ] <- x i <- i + 1 target <- target + dt if (i > n) return(xmat) } j <- sample(v, 1, prob = h) x <- x + S[, j] } } get_results <- function(params, nreps = 100, outs, times, raw = FALSE) { tseq <- 0:max(times) arra <- array(0, dim = c(max(tseq) + 1, 3, nreps)) for(i in 1:nreps) arra[,,i] <- gillespied(N, T = max(times) + 1 + 0.001, dt = 1, th = params) if(raw) return(arra) collected <- list() for (i in 1:nreps) { relev <- c(arra[times+1, which(c("S","I","R") %in% outs), i]) names <- c(outer(outs, times, function(x, y) paste0(x, y))) relev <- setNames(relev, names) collected[[i]] <- relev } input_dat <- setNames(data.frame(matrix(rep(params, nreps), ncol = length(params), byrow = TRUE)), names(params)) return(cbind(input_dat, do.call('rbind', collected))) } Num <- 1000 N <- list() N$M <- c(950, 50, 0) N$Pre <- matrix(c(1,0,0,0,1,0,0,0,1), nrow = 3, byrow = TRUE) N$Post <- matrix(c(0,1,0,0,0,1,1,0,0), nrow = 3, byrow = TRUE) N$h <- function(x, t, th = rep(1, 3)) { return( c(th[1] * x[1] * x[2]/Num, th[2] * x[2], th[3] * x[3]) ) } ## Create training and validation sets using Gillespie algorithm set.seed(123) training_stoch <- do.call('rbind.data.frame', apply(training[,names(ranges)], 1, function(x) get_results(unlist(x), outs = c("S", "I", "R"), times = c(10)))) |> setNames(c(names(ranges), "nS", "nI", "nR")) valid_stoch <- do.call('rbind.data.frame', apply(valid[,names(ranges)], 1, function(x) get_results(unlist(x), outs = c("S", "I", "R"), times = c(10)))) |> setNames(c(names(ranges), "nS", "nI", "nR")) # Create some emulators for the data stoch_emulators <- emulator_from_data(training_stoch, names(targets), ranges, emulator_type = "variance", order = 3) # Plot the outputs... ####### FIGURE 14 ####### emulator_plot(subset_emulators(stoch_emulators, c("nS", "nR"))$variance) #### NOT USED IN PAPER ##### # emulator_plot(subset_emulators(stoch_emulators, c("nS", "nR")), 'var') # ..note that, even at training points, the uncertainty is non-zero. train_sample <- unique(training_stoch[,names(ranges)])[1:10,] stoch_emulators$expectation$nI$get_cov(train_sample) # Validation diagnostics call is identical to deterministic case ###### FIGURE 15a ####### validation_diagnostics(stoch_emulators, targets, valid_stoch) # As is point generation stoch_points1 <- generate_new_design(stoch_emulators, 90, targets) ######### FIGURE 15b ######### # Long (maybe excessive) runtime! ## Multiple waves of history matching with stochastic model set.seed(123) training_stoch2 <- do.call('rbind.data.frame', apply(stoch_points1[1:30,], 1, function(x) get_results(unlist(x), outs = c("S", "I", "R"), times = c(10)))) |> setNames(c(names(ranges), "nS", "nI", "nR")) valid_stoch2 <- do.call('rbind.data.frame', apply(stoch_points1[31:90,], 1, function(x) get_results(unlist(x), outs = c("S", "I", "R"), times = c(10)))) |> setNames(c(names(ranges), "nS", "nI", "nR")) stoch_emulators2 <- emulator_from_data(training_stoch2, names(targets), ranges, check.ranges = T, emulator_type = "variance", order = 3) validation_diagnostics(stoch_emulators2, targets, valid_stoch2) stoch_points2 <- generate_new_design(list(stoch_emulators2, stoch_emulators), 90, targets) wave_points(list(stoch_points2, SIRMultiWaveData[[3]][,names(ranges)]), names(ranges), p_size = 0.8) + viridis::scale_fill_viridis(direction = -1, discrete = TRUE, alpha = 0.5, name = "Method", labels = c("Stochastic", "Deterministic"))