change output in video (for AUDPC_rel)

R/output.R

@@ -450,7 +450,7 @@ epid_output <- function(types = "all", time_param, Npatho, area, rotation, cropt
           ylab_output <- expression("Equivalent nb of diseased hosts / time step / m"^2) ## expression('title'^2)
           round_ylab_output <- 2
         }, "audpc_rel" = {
-          main_output <- expression("Relative AUDPC")
+          main_output <- "Relative AUDPC"
           ylab_output <- "Proportion of diseased hosts: (I+R)/N"
           round_ylab_output <- 2
         }, "gla" = {
@@ -873,7 +873,6 @@ evol_output <- function(types = "all", time_param, Npoly, cultivars_param, genes
 #' as if cultivated in pure crops,
 #' \item cultivars_genes_list = a list containing, for each host genotype, the indices of carried resistance genes.
 #' }
-#' @param audpc100S AUDPC of a 100% susceptible landscape, used as a reference.
 #' @param keyDates a vector of times (in time steps) where to draw vertical lines in the AUDPC graphic. Usually
 #' used to delimit durabilities of the resistance genes. No line is drawn if keyDates=NULL (default).
 #' @param nMapPY an integer specifying the number of epidemic maps per year to generate.
@@ -895,66 +894,85 @@ evol_output <- function(types = "all", time_param, Npoly, cultivars_param, genes
 #' @importFrom parallel detectCores
 #' @importFrom grDevices colorRampPalette dev.off graphics.off gray png tiff
 #' @export
-video <- function(audpc, time_param, Npatho, landscape, area, rotation, croptypes, croptype_names = c(), cultivars_param, audpc100S, keyDates = NULL,
-                  nMapPY = 5, path = getwd()) {
-
+video <- function(audpc, time_param, Npatho, landscape, area, rotation, croptypes, croptype_names = c(), cultivars_param
+                  # , audpc100S
+                  , keyDates = NULL,nMapPY = 5, path = getwd()) {
+  
   if (system("ffmpeg -version", ignore.stdout = TRUE) == 127) {
     stop("You need to install ffmpeg before generating videos. Go to https://ffmpeg.org/download.html")
   }
-
+  
   ## Time & graphic parameters
   Nyears <- time_param$Nyears
   nTSpY <- time_param$nTSpY
   nTS <- Nyears * nTSpY ## Total number of time-steps
-
+  
   ## Landscape
   rotation <- data.frame(rotation)
   Npoly <- length(area)
-
+  
   ## croptype parameters (for legend)
   if (length(croptype_names) == 0) {
     croptype_names <- paste("Croptype", unique(croptypes$croptypeID))
   }
   Ncroptypes <- length(croptype_names)
-
+  
   ## Host parameters
   cultivar_names <- cultivars_param$name
-  cultivars_genes_list <- cultivars_param$cultivars_genes_list
-  max_density <- cultivars_param$max_density
-  Nhost <- length(max_density)
-
+  # cultivars_genes_list <- cultivars_param$cultivars_genes_list
+  # max_density <- cultivars_param$max_density
+  Nhost <- length(cultivar_names)
+  
   ## Calculation of the carrying capacity
-  K <- array(dim = c(Npoly, Nhost, Nyears))
-  for (y in 1:Nyears) {
-    if (ncol(rotation) == 1) {
-      index_year <- 1
-    } else {
-      index_year <- y
-    }
-    ts_year <- ((y - 1) * nTSpY + 1):(y * nTSpY)
-
-    for (poly in 1:Npoly) {
-      indices_croptype <- grep(rotation[poly, index_year], croptypes$croptypeID)
-      for (i in indices_croptype) {
-        index_host <- croptypes[i, "cultivarID"] + 1 ## +1 because C indices start at 0
-        prop <- croptypes[i, "proportion"]
-        K[poly, index_host, y] <- floor(area[poly] * max_density[index_host] * prop)
-      }
-    } ## for poly
-  } ## for y
-
+  # K <- array(dim = c(Npoly, Nhost, Nyears))
+  # for (y in 1:Nyears) {
+  #   if (ncol(rotation) == 1) {
+  #     index_year <- 1
+  #   } else {
+  #     index_year <- y
+  #   }
+  #   ts_year <- ((y - 1) * nTSpY + 1):(y * nTSpY)
+  #   
+  #   for (poly in 1:Npoly) {
+  #     indices_croptype <- grep(rotation[poly, index_year], croptypes$croptypeID)
+  #     for (i in indices_croptype) {
+  #       index_host <- croptypes[i, "cultivarID"] + 1 ## +1 because C indices start at 0
+  #       prop <- croptypes[i, "proportion"]
+  #       K[poly, index_host, y] <- floor(area[poly] * max_density[index_host] * prop)
+  #     }
+  #   } ## for poly
+  # } ## for y
+  
   #### IMPORTATION OF THE SIMULATION OUTPUT
+  H <- as.list(1:nTS)
+  L <- as.list(1:nTS)
   I <- as.list(1:nTS)
   R <- as.list(1:nTS)
   index <- 0
-
+  
+  # print("Reading binary files to compute epidemiological model outputs")
   for (year in 1:Nyears) {
+    binfileH <- file(paste(path, sprintf("/H-%02d", year), ".bin", sep = ""), "rb")
+    H.tmp <- readBin(con = binfileH, what = "int", n = Npoly * Nhost * nTSpY, size = 4, signed = T, endian = "little")
+    close(binfileH)
+    
+    binfileL <- file(paste(path, sprintf("/L-%02d", year), ".bin", sep = ""), "rb")
+    L.tmp <- readBin(con = binfileL, what = "int", n = Npoly * Npatho * Nhost * nTSpY, size = 4, signed = T, endian = "little")
+    close(binfileL)
     binfileI <- file(paste(path, sprintf("/I-%02d", year), ".bin", sep = ""), "rb")
     I.tmp <- readBin(con = binfileI, what = "int", n = Npoly * Npatho * Nhost * nTSpY, size = 4, signed = T, endian = "little")
+    close(binfileI)
+    
     binfileR <- file(paste(path, sprintf("/R-%02d", year), ".bin", sep = ""), "rb")
     R.tmp <- readBin(con = binfileR, what = "int", n = Npoly * Npatho * Nhost * nTSpY, size = 4, signed = T, endian = "little")
-
+    close(binfileR)
+    
     for (t in 1:nTSpY) {
+      H[[t + index]] <- matrix(H.tmp[((Nhost * Npoly) * (t - 1) + 1):(t * (Nhost * Npoly))], ncol = Nhost, byrow = T)
+      L[[t + index]] <- array(
+        data = L.tmp[((Npatho * Npoly * Nhost) * (t - 1) + 1):(t * (Npatho * Npoly * Nhost))],
+        dim = c(Nhost, Npatho, Npoly)
+      )
       I[[t + index]] <- array(
         data = I.tmp[((Npatho * Npoly * Nhost) * (t - 1) + 1):(t * (Npatho * Npoly * Nhost))],
         dim = c(Nhost, Npatho, Npoly)
@@ -963,12 +981,10 @@ video <- function(audpc, time_param, Npatho, landscape, area, rotation, croptype
         data = R.tmp[((Npatho * Npoly * Nhost) * (t - 1) + 1):(t * (Npatho * Npoly * Nhost))],
         dim = c(Nhost, Npatho, Npoly)
       )
-    }
-
+    } ## for t
+    
     index <- index + nTSpY
-    close(binfileI)
-    close(binfileR)
-  }
+  } ## for year
 
 
   print("Generate video (and create directory maps)")
@@ -1008,17 +1024,21 @@ video <- function(audpc, time_param, Npatho, landscape, area, rotation, croptype
       index_year <- y
     }
 
-    K_poly <- apply(as.data.frame(K[, , index_year]), 1, sum, na.rm = TRUE) ## as.data.frame in case Nhost==1 & Npatho==1
+    # K_poly <- apply(as.data.frame(K[, , index_year]), 1, sum, na.rm = TRUE) ## as.data.frame in case Nhost==1 & Npatho==1
     for (d in round(seq(1, nTSpY, length.out = nMapPY))) {
       subtitle <- paste("Year =", y, "   Day =", d)
       ## Proportion of each type of host relative to carrying capacity
       ts <- (y - 1) * nTSpY + d
 
-      propI <- apply(I[[ts]], 3, sum) / K_poly
-      propR <- apply(R[[ts]], 3, sum) / K_poly
-      propIR <- (propI + propR)
-      intvlsIR <- findInterval(propIR, intvls)
+      # propI <- apply(I[[ts]], 3, sum) / K_poly
+      # propR <- apply(R[[ts]], 3, sum) / K_poly
+      # propIR <- (propI + propR)
+      # intvlsIR <- findInterval(propIR, intvls)
 
+      IR_poly <- apply(I[[ts]], 3, sum) + apply(R[[ts]], 3, sum) 
+      N_poly  <- apply(H[[ts]], 1, sum) + apply(L[[ts]], 3, sum) + apply(I[[ts]], 3, sum) + apply(R[[ts]], 3, sum)
+      intvlsIR <- findInterval(IR_poly / N_poly, intvls)
+      
       png(
         filename = paste(path, "/maps/HLIR_", sprintf("%02d", y), "-", sprintf("%03d", d), ".png", sep = ""),
         width = 700, height = 350, units = "mm", res = 72, type = "cairo", bg = "white", antialias = "default"
@@ -1027,17 +1047,20 @@ video <- function(audpc, time_param, Npatho, landscape, area, rotation, croptype
 
       ## moving AUDPC
       plot(0, 0,
-        type = "n", bty = "n", xlim = c(1, Nyears), ylim = c(0, 1), xaxt = "n", yaxt = "n", xlab = "", ylab = "",
-        main = "Disease severity relative to a fully susceptible landscape"
+        type = "n", bty = "n", xlim = c(1, Nyears), ylim = c(0, 1), xaxt = "n", yaxt = "n"
+        , xlab = "", ylab = "", main = "Average disease severity"
       )
       axis(1, at = seq(1, Nyears, by = (Nyears - 1) %/% 10 + 1))
       ticksMarks <- round(seq(0, 1, length.out = 5), 2)
       axis(2, at = ticksMarks, labels = paste(100 * ticksMarks, "%"), las = 1)
       title(xlab = "Years", mgp = c(2, 1, 0))
+      title(ylab = "Proportion of diseased hosts: (I+R)/N", mgp = c(3.5, 1, 0))
       for (host in 1:Nhost) {
-        lines(1:y, audpc[1:y, host] / audpc100S, type = "o", lwd = 2, lty = host, pch = PCH[host], col = GRAY[host])
+        # lines(1:y, audpc[1:y, host] / audpc100S, type = "o", lwd = 2, lty = host, pch = PCH[host], col = GRAY[host])
+        lines(1:y, audpc[1:y, host], type = "o", lwd = 2, lty = host, pch = PCH[host], col = GRAY[host])
       }
-      lines(1:y, audpc[1:y, "total"] / audpc100S, type = "o", lwd = 2, lty = LTY.tot, pch = PCH.tot, col = COL.tot)
+      # lines(1:y, audpc[1:y, "total"] / audpc100S, type = "o", lwd = 2, lty = LTY.tot, pch = PCH.tot, col = COL.tot)
+      lines(1:y, audpc[1:y, "total"], type = "o", lwd = 2, lty = LTY.tot, pch = PCH.tot, col = COL.tot)
       ## Add lines for durability
       if (!is.null(keyDates)) {
         for (k in 1:length(keyDates)) {

R/simul_landsepi.R

@@ -348,8 +348,10 @@ simul_landsepi <- function(seed = 12345, time_param = list(Nyears = 20, nTSpY =
   ## Video
   if (videoMP4 & !is.null(epid_res[["audpc"]])) {
     video(
-      epid_res[["audpc"]], time_param, Npatho, landscape, area, rotation, croptypes_cultivars_prop,
-      croptype_names, cultivars_param, audpc100S, evol_res$durability
+      epid_res[["audpc_rel"]], time_param, Npatho, landscape, area, rotation, croptypes_cultivars_prop,
+      croptype_names, cultivars_param
+      # , audpc100S
+      , evol_res$durability
     )
   }