output.R 48.7 KB
Newer Older
Jean-Francois Rey's avatar
Jean-Francois Rey committed
1
# Part of the landsepi R package.
2
3
4
# Copyright (C) 2017 Loup Rimbaud <loup.rimbaud@inrae.fr>
#                    Julien Papaix <julien.papaix@inrae.fr>
#                    Jean-François Rey <jean-francois.rey@inrae.fr>
Jean-Francois Rey's avatar
Jean-Francois Rey committed
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software Foundation, Inc.,i
# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#
Loup Rimbaud's avatar
Loup Rimbaud committed
20
#' @title Generation of epidemiological and economic model outputs
Jean-Francois Rey's avatar
Jean-Francois Rey committed
21
#' @name epid_output
Loup Rimbaud's avatar
Loup Rimbaud committed
22
#' @description Generates epidemiological and economic outputs from model simulations.
Jean-Francois Rey's avatar
Jean-Francois Rey committed
23
#' @param types a character string (or a vector of character strings if several outputs are to be computed)
Loup Rimbaud's avatar
Loup Rimbaud committed
24
#' specifying the type of outputs to generate (see details):\itemize{
25
#' \item "audpc": Area Under Disease Progress Curve
26
27
#' \item "audpc_rel": Relative Area Under Disease Progress Curve
#' \item "gla": Green Leaf Area
28
#' \item "gla_rel": Relative Green Leaf Area
29
30
31
32
#' \item "eco_yield": Total crop yield
#' \item "eco_cost": Operational crop costs
#' \item "eco_product": Crop products
#' \item "eco_margin": Margin (products - operational costs)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
33
#' \item "HLIR_dynamics", "H_dynamics", "L_dynamics", "IR_dynamics", "HLI_dynamics", etc.: Epidemic dynamics
Loup Rimbaud's avatar
Loup Rimbaud committed
34
35
#' related to the specified sanitary status (H, L, I or R and all their combinations). 
#' Graphics only, works only if graphic=TRUE.
Loup Rimbaud's avatar
Loup Rimbaud committed
36
37
#' \item "all": compute all these outputs (default).
#' }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
38
#' @param time_param list of simulation parameters:\itemize{
Loup Rimbaud's avatar
Loup Rimbaud committed
39
40
41
42
#' \item Nyears = number cropping seasons,
#' \item nTSpY = number of time-steps per cropping season.
#' }
#' @param Npatho number of pathogen genotypes.
Loup Rimbaud's avatar
Loup Rimbaud committed
43
#' @param area a vector containing polygon areas (must be in square meters).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
44
45
46
#' @param rotation a dataframe containing for each field (rows) and year (columns, named "year_1", "year_2", etc.),
#' the index of the cultivated croptype. Importantly, the matrix must contain 1 more column than the real number
#' of simulated years.
Loup Rimbaud's avatar
Loup Rimbaud committed
47
48
#' @param croptypes a dataframe with three columns named 'croptypeID' for croptype index,
#' 'cultivarID' for cultivar index and 'proportion' for the proportion of the cultivar within the croptype.
Loup Rimbaud's avatar
Loup Rimbaud committed
49
50
#' @param cultivars_param list of parameters associated with each host genotype (i.e. cultivars):  
#' \itemize{
Jean-Francois Rey's avatar
Jean-Francois Rey committed
51
#' \item name = vector of cultivar names,
52
53
54
55
#' \item initial_density = vector of host densities (per square meter) at the beginning of the cropping season 
#' as if cultivated in pure crop,
#' \item max_density = vector of maximum host densities (per square meter) at the end of the cropping season 
#' as if cultivated in pure crop,
Loup Rimbaud's avatar
Loup Rimbaud committed
56
57
58
#' \item cultivars_genes_list = a list containing, for each host genotype, the indices of carried resistance genes.
#' }
#' @param eco_param a list of economic parameters for each host genotype as if cultivated in pure crop:\itemize{
59
60
61
62
63
#' \item yield_perHa = a dataframe of 4 columns for the theoretical yield associated with hosts in sanitary status H, L, I and R,
#' as if cultivated in pure crops, and one row per host genotype 
#' (yields are expressed in weight or volume units / ha / cropping season),
#' \item planting_cost_perHa = a vector of planting costs (in monetary units / ha / cropping season),
#' \item market_value = a vector of market values of the production (in monetary units / weight or volume unit).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
64
#' }
65
#' @param GLAnoDis the value of absolute GLA in absence of disease (required to compute economic outputs).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
66
#' @param ylim_param a list of graphical parameters for each required output: bounds for y-axes for
67
#' audpc, gla, gla_rel, eco_cost, eco_yield, eco_product, eco_margin.
68
#' @param writeTXT a logical indicating if the output is written in a text file (TRUE) or not (FALSE).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
69
#' @param graphic a logical indicating if a tiff graphic of the output is generated (only if more than one year is simulated).
70
#' @param path path of text file (if writeTXT = TRUE) and tiff graphic (if graphic = TRUE) to be generated.
Loup Rimbaud's avatar
Loup Rimbaud committed
71
#' @details Outputs are computed every year for every cultivar as well as for the whole landscape. \describe{
Jean-Francois Rey's avatar
Jean-Francois Rey committed
72
73
#' \item{\strong{Epidemiological outputs.}}{
#' The epidemiological impact of pathogen spread can be evaluated by different measures: \enumerate{
74
75
76
77
78
79
#' \item Area Under Disease Progress Curve (AUDPC): average number of diseased host individuals (status I + R)
#' per time step and square meter.
#' \item Relative Area Under Disease Progress Curve (AUDPCr): average proportion of diseased host individuals 
#' (status I + R) relative to the total number of existing hosts (H+L+I+R).
#' \item Green Leaf Area (GLA): average number of healthy host individuals (status H) per time step and per square meter.
#' \item Relative Green Leaf Area (GLAr): average proportion of healthy host individuals (status H) relative to the total number
Loup Rimbaud's avatar
Loup Rimbaud committed
80
#' of existing hosts (H+L+I+R).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
81
82
83
84
#' }
#'  }
#' \item{\strong{Economic outputs.}}{
#' The economic outcome of a simulation can be evaluated using: \enumerate{
85
86
87
88
89
90
91
#' \item Crop yield: yearly crop yield (e.g. grains, fruits, wine) in weight (or volume) units
#' per hectare (depends on the number of productive hosts and associated theoretical yield).
#' \item Crop products: yearly product generated from sales, in monetary units per hectare
#' (depends on crop yield and market value).
#' \item Operational crop costs: yearly costs associated with crop planting in monetary units per hectare 
#' (depends on initial host density and planting cost).
#' \item Crop margin, i.e. products - operational costs, in monetary units per hectare.
Jean-Francois Rey's avatar
Jean-Francois Rey committed
92
93
94
#' }
#'  }
#'  }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
95
#' @return A list containing, for each required type of output, a matrix summarising the output for each year and cultivar
Loup Rimbaud's avatar
Loup Rimbaud committed
96
#' (as well as the whole landscape).
97
#' Each matrix can be written in a txt file (if writeTXT=TRUE), and illustrated in a graphic (if graphic=TRUE).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
98
99
#' @references Rimbaud L., Papaïx J., Rey J.-F., Barrett L. G. and Thrall P. H. (2018). Assessing the durability and efficiency of
#' landscape-based strategies to deploy plant resistance to pathogens. \emph{PLoS Computational Biology} 14(4):e1006067.
Loup Rimbaud's avatar
Loup Rimbaud committed
100
#' @seealso \link{evol_output}
Jean-Francois Rey's avatar
Jean-Francois Rey committed
101
#' @examples
102
#' \dontrun{
Jean-Francois Rey's avatar
Jean-Francois Rey committed
103
104
#' demo_landsepi()
#' }
105
#' @include RcppExports.R Math-Functions.R graphics.R
Jean-Francois Rey's avatar
Jean-Francois Rey committed
106
107
108
109
#' @importFrom sf st_read
#' @importFrom grDevices colorRampPalette dev.off graphics.off gray png tiff
#' @importFrom utils write.table
#' @export
Jean-Francois Rey's avatar
Jean-Francois Rey committed
110
111
112
epid_output <- function(types = "all", time_param, Npatho, area, rotation, croptypes, cultivars_param, eco_param,
                        GLAnoDis = cultivars_param$max_density[1] ## true for non-growing host
                        , ylim_param = list(
Loup's avatar
Loup committed
113
                          audpc = c(0, 0.76), #0.38),
114
115
                          audpc_rel = c(0, 1),
                          gla = c(0, 1.48),
Jean-Francois Rey's avatar
Jean-Francois Rey committed
116
117
                          gla_rel = c(0, 1),
                          eco_cost = c(0, NA),
118
                          eco_yield = c(0, NA),
Jean-Francois Rey's avatar
Jean-Francois Rey committed
119
                          eco_product = c(0, NA),
120
                          eco_margin = c(NA, NA)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
121
                        ),
122
                        writeTXT = TRUE, graphic = TRUE, path = getwd()) {
123
  valid_outputs <- c("audpc", "audpc_rel", "gla", "gla_rel", "eco_yield", "eco_product", "eco_cost", "eco_margin", "HLIR_dynamics")
Loup Rimbaud's avatar
Loup Rimbaud committed
124
  if (types[1] == "all") {
125
126
    types <- valid_outputs
  }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
127
  valid_outputs <- c(valid_outputs, paste0(c("H", "L", "I", "R", "HL", "HI", "HR", "LI", "LR", "IR", "HLI", "HLR", "HIR", "LIR", "HLIR"), "_dynamics"))
128
  if (is.na(sum(match(types, valid_outputs)))) {
129
    stop(paste("Error: valid types of outputs are", paste(valid_outputs, collapse = ", ")))
130
131
132
133
134
135
  }

  ## Time parameters
  Nyears <- time_param$Nyears
  nTSpY <- time_param$nTSpY
  nTS <- Nyears * nTSpY ## Total number of time-steps
Jean-Francois Rey's avatar
Jean-Francois Rey committed
136

Loup Rimbaud's avatar
Loup Rimbaud committed
137
  ## Landscape
Jean-Francois Rey's avatar
Jean-Francois Rey committed
138
  rotation <- data.frame(rotation)
139
  areaTot <- sum(area)
140
141
142
143
  Npoly <- length(area)

  ## Host parameters
  cultivar_names <- cultivars_param$name
Loup Rimbaud's avatar
Loup Rimbaud committed
144
  initial_density <- cultivars_param$initial_density
145
146
  max_density <- cultivars_param$max_density
  cultivars_genes_list <- cultivars_param$cultivars_genes_list
147
  yield_perIndivPerTS <- eco_param$yield_perHa * 1E-4 / nTSpY / GLAnoDis
148
  planting_cost_perIndiv <- eco_param$planting_cost_perHa * 1E-4 / initial_density
Loup Rimbaud's avatar
Loup Rimbaud committed
149
150
  market_value <- eco_param$market_value
  Nhost <- length(initial_density)
151
152


153
154
155
  ## Calculation of the carrying capacity and number of exisiting hosts
  K <- array(dim = c(Npoly, Nhost, Nyears)) ## for audpc
  C <- array(dim = c(Npoly, Nhost, Nyears)) ## for cost
Jean-Francois Rey's avatar
Jean-Francois Rey committed
156
  area_host <- matrix(0, nrow = Nyears, ncol = Nhost + 1) ## for cost
157
158
159
160
161
162
163
  for (y in 1:Nyears) {
    if (ncol(rotation) == 1) {
      index_year <- 1
    } else {
      index_year <- y
    }
    ts_year <- ((y - 1) * nTSpY + 1):(y * nTSpY)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
164

165
166
167
168
169
170
    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)
171
172
        C[poly, index_host, y] <- area[poly] * initial_density[index_host] * prop
        area_host[y, index_host] <- area_host[y, index_host] + area[poly]
173
      }
174
175
176
177
178
    } ## for poly
  } ## for y
  K_host <- apply(K, c(2, 3), sum, na.rm = TRUE)
  C_host <- apply(C, c(2, 3), sum, na.rm = TRUE)
  C_host[C_host == 0] <- NA
Jean-Francois Rey's avatar
Jean-Francois Rey committed
179
  area_host[, Nhost + 1] <- areaTot
Jean-Francois Rey's avatar
Jean-Francois Rey committed
180

181
  ## IMPORTATION OF THE SIMULATION OUTPUT (only those required depending on parameter "types")
182
  requireH <- sum(is.element(substr(types, 1, 3), c("gla", "eco"))) > 0
183
184
  requireL <- sum(is.element(types, c("audpc_rel", "gla_rel", "eco_yield", "eco_product", "eco_cost", "eco_margin"))) > 0
  requireIR <- sum(is.element(types, c("audpc", "audpc_rel", "gla_rel", "eco_yield", "eco_product", "eco_cost", "eco_margin"))) > 0
Jean-Francois Rey's avatar
Jean-Francois Rey committed
185
186
187
  if (graphic & sum(substr(types, nchar(types) - 7, nchar(types)) == "dynamics") > 0) {
    requireH <- 1
    requireL <- 1
188
189
    requireIR <- 1
  }
190
191
192
193
194
195
196
197
  H <- as.list(1:nTS)
  # Hjuv <- as.list(1:nTS)
  # P <- as.list(1:nTS)
  L <- as.list(1:nTS)
  I <- as.list(1:nTS)
  R <- as.list(1:nTS)
  index <- 0

Loup Rimbaud's avatar
Loup Rimbaud committed
198
  # print("Reading binary files to compute epidemiological model outputs")
199
200
201
  for (year in 1:Nyears) {
    # print(paste("year", year, "/", Nyears))
    if (requireH) {
202
203
      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")
204
205
206
207
208
209
210
211
212
213
214
215
216
217
      close(binfileH)
    }
    # binfileHjuv = file(paste(path, sprintf("/Hjuv-%02d", year), ".bin",sep=""), "rb")
    #  Hjuv.tmp <- readBin(con=binfileHjuv, what="int", n=Npoly*Nhost*nTSpY, size = 4, signed=T,endian="little")
    # close(binfileHjuv)
    # binfileP <- file(paste(path, sprintf("/P-%02d", year), ".bin", sep = ""), "rb")
    # P.tmp <- readBin(con = binfileP,
    #                  what = "int",
    #                  n = Npoly * Npatho * nTSpY,
    #                  size = 4,
    #                  signed = T,
    #                  endian = "little")
    # close(binfileP)
    if (requireL) {
218
219
      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")
220
221
222
      close(binfileL)
    }
    if (requireIR) {
223
224
      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")
225
226
      close(binfileI)

227
228
      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")
229
230
231
232
      close(binfileR)
    }

    for (t in 1:nTSpY) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
233
      if (requireH) {
234
        H[[t + index]] <- matrix(H.tmp[((Nhost * Npoly) * (t - 1) + 1):(t * (Nhost * Npoly))], ncol = Nhost, byrow = T)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
235
      }
236
237
      # Hjuv[[t + index]] <- matrix(Hjuv.tmp[((Nhost*Npoly)*(t-1)+1):(t*(Nhost*Npoly))],ncol=Nhost,byrow=T)
      # P[[t + index]] <- matrix(P.tmp[((Npatho * Npoly) * (t - 1) + 1):(t * (Npatho * Npoly))], ncol = Npatho, byrow = T)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
238
      if (requireL) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
239
240
241
242
        L[[t + index]] <- array(
          data = L.tmp[((Npatho * Npoly * Nhost) * (t - 1) + 1):(t * (Npatho * Npoly * Nhost))],
          dim = c(Nhost, Npatho, Npoly)
        )
Jean-Francois Rey's avatar
Jean-Francois Rey committed
243
244
      }
      if (requireIR) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
245
246
247
248
249
250
251
252
        I[[t + index]] <- array(
          data = I.tmp[((Npatho * Npoly * Nhost) * (t - 1) + 1):(t * (Npatho * Npoly * Nhost))],
          dim = c(Nhost, Npatho, Npoly)
        )
        R[[t + index]] <- array(
          data = R.tmp[((Npatho * Npoly * Nhost) * (t - 1) + 1):(t * (Npatho * Npoly * Nhost))],
          dim = c(Nhost, Npatho, Npoly)
        )
Jean-Francois Rey's avatar
Jean-Francois Rey committed
253
      }
254
255
256
257
258
259
260
261
262
263
264
265
266
    } ## for t

    index <- index + nTSpY
  } ## for year

  #### HLIR DYNAMIC
  H_host <- NULL
  L_host <- NULL
  I_host <- NULL
  R_host <- NULL
  if (requireH) {
    for (t in 1:nTS) {
      H_host <- cbind(H_host, apply(H[[t]], 2, sum))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
267
    }
268
269
270
271
  }
  if (requireL) {
    for (t in 1:nTS) {
      L_host <- cbind(L_host, apply(L[[t]], 1, sum))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
272
    }
273
274
275
276
277
  }
  if (requireIR) {
    for (t in 1:nTS) {
      I_host <- cbind(I_host, apply(I[[t]], 1, sum))
      R_host <- cbind(R_host, apply(R[[t]], 1, sum))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
278
    }
279
280
  }
  N_host <- H_host + L_host + I_host + R_host
Jean-Francois Rey's avatar
Jean-Francois Rey committed
281
282


283
284
285
  #####       COMPUTATION OF OUTPUTS
  ## as.numeric to allow memory allocation of long integer
  res <- list()
Jean-Francois Rey's avatar
Jean-Francois Rey committed
286

287
288
289
  ## Colours (+1 to avoid picking white)
  # COL <- c("#FF5555","#4F94CD","darkolivegreen4","#CD950C","black") ## colors: red, blue, green, orange, black
  grad_red <- colorRampPalette(c("#FF5555", "white"))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
290
  RED <- grad_red(Nhost + 1)
291
  grad_blue <- colorRampPalette(c("#4F94CD", "white"))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
292
  BLUE <- grad_blue(Nhost + 1)
293
  grad_green <- colorRampPalette(c("darkolivegreen4", "white"))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
294
295
296
297
  GREEN <- grad_green(Nhost + 1)
  grad_grey <- colorRampPalette(c("black", "white")) # "gray30"
  GRAY <- grad_grey(Nhost + 1)

298
  for (type in types) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
299

300
    ## Epidemic dynamics (H, L, I, R)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
301
    if (graphic & substr(type, nchar(type) - 7, nchar(type)) == "dynamics") {
302
      varToPlot <- strsplit(type, "_")[[1]][1]
Jean-Francois Rey's avatar
Jean-Francois Rey committed
303
304
305
306
307
308
      varToLegend <- strsplit(varToPlot, "")[[1]]

      tiff(
        filename = paste(path, "/", type, ".tiff", sep = ""),
        width = 180, height = 110, units = "mm", compression = "lzw", res = 300
      )
309
      par(xpd = NA, cex = .9, mar = c(4, 4.3, 3, 9))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327

      plot(0, 0,
        type = "n", xaxt = "n", yaxt = "n", bty = "n", xlim = c(1, nTS), ylim = c(0, 1), main = "Epidemic dynamics",
        ylab = "Proportion of hosts relative to carrying capacity", xlab = ""
      )
      for (host in 1:Nhost) {
        if (grepl("H", varToPlot)) {
          lines(1:nTS, H_host[host, ] / rep(K_host[host, ], each = nTSpY), col = GREEN[host], lty = host)
        }
        if (grepl("L", varToPlot)) {
          lines(1:nTS, L_host[host, ] / rep(K_host[host, ], each = nTSpY), col = BLUE[host], lty = host)
        }
        if (grepl("I", varToPlot)) {
          lines(1:nTS, I_host[host, ] / rep(K_host[host, ], each = nTSpY), col = RED[host], lty = host)
        }
        if (grepl("R", varToPlot)) {
          lines(1:nTS, R_host[host, ] / rep(K_host[host, ], each = nTSpY), col = GRAY[host], lty = host)
        }
328
      }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
329
330
331
      if (Nyears == 1) {
        axis(1, at = round(seq(1, nTS, length.out = 8)))
        title(xlab = "Days")
332
      } else {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
333
334
        axis(1, at = seq(1, nTS + 1, nTSpY * ((Nyears - 1) %/% 10 + 1)), labels = seq(0, Nyears, ((Nyears - 1) %/% 10 + 1)))
        title(xlab = "Years")
335
      }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
336
337
338
339
340
341
342
343
344
345
346
347
      axis(2, at = seq(0, 1, .2), las = 1)
      legend(nTS * 1.05, .8,
        title.adj = -.01, bty = "n", title = "Status:", legend = varToLegend, border = NA,
        fill = c(GREEN[1], BLUE[1], RED[1], GRAY[1])[is.element(c("H", "L", "I", "R"), varToLegend)], cex = 0.9
      )
      legend(nTS * 1.05, .3,
        title.adj = -.01, bty = "n", title = "Cultivar:", legend = cultivar_names,
        lty = 1:Nhost, lwd = 2, col = GRAY[1:Nhost], cex = 0.9
      )
      dev.off()
    } else { ## Other epidemiological outputs

348
349
350
      output_matrix <- data.frame(matrix(NA, ncol = Nhost + 1, nrow = Nyears))
      rownames(output_matrix) <- paste("year_", 1:Nyears, sep = "")
      colnames(output_matrix) <- c(cultivar_names, "total")
Jean-Francois Rey's avatar
Jean-Francois Rey committed
351

352
353
354
355
      for (y in 1:Nyears) {
        ts_year <- ((y - 1) * nTSpY + 1):(y * nTSpY)
        ## metrics for each host
        for (host in 1:Nhost) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
356
357
358
          switch(type,
            "audpc" = {
              numerator <- sum(as.numeric(I_host[host, ts_year]), as.numeric(R_host[host, ts_year]))
Loup's avatar
Loup committed
359
360
361
362
363
              denominator <- nTSpY * areaTot   ## sum(as.numeric(K_host[host, y] * length(ts_year)))
            },
            "audpc_rel" = {
              numerator <- sum(as.numeric(I_host[host, ts_year]), as.numeric(R_host[host, ts_year]))
              denominator <- sum(as.numeric(N_host[host, ts_year]))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
364
            },
365
            "gla" = {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
              numerator <- sum(as.numeric(H_host[host, ts_year]))
              denominator <- nTSpY * areaTot
            },
            "gla_rel" = {
              numerator <- sum(as.numeric(H_host[host, ts_year]))
              denominator <- sum(as.numeric(N_host[host, ts_year]))
            },
            { ## Default instructions (i.e. if substr(type,1,3)=="eco")
              numerator <- sum(
                yield_perIndivPerTS[host, "H"] * as.numeric(H_host[host, ts_year]),
                yield_perIndivPerTS[host, "L"] * as.numeric(L_host[host, ts_year]),
                yield_perIndivPerTS[host, "I"] * as.numeric(I_host[host, ts_year]),
                yield_perIndivPerTS[host, "R"] * as.numeric(R_host[host, ts_year])
              )
              denominator <- 1
            }
          )

384
385
386
387
          if (denominator > 0 & K_host[host, y] > 0) {
            output_matrix[y, host] <- numerator / denominator
          } ## if >0
        } ## for host
Jean-Francois Rey's avatar
Jean-Francois Rey committed
388

389
        ## metrics for all landscape
Jean-Francois Rey's avatar
Jean-Francois Rey committed
390
        switch(type,
Loup's avatar
Loup committed
391
392
393
394
395
396
397
398
          # "audpc" = {
          #   numerator_tot <- sum(as.numeric(I_host[, ts_year]), as.numeric(R_host[, ts_year]))
          #   denominator_tot <- sum(as.numeric(K_host[, y] * length(ts_year)))
          #   if (denominator_tot > 0) {
          #     output_matrix[y, "total"] <- numerator_tot / denominator_tot
          #   }
          # },
          "audpc_rel" = {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
399
            numerator_tot <- sum(as.numeric(I_host[, ts_year]), as.numeric(R_host[, ts_year]))
Loup's avatar
Loup committed
400
            denominator_tot <- sum(as.numeric(N_host[, ts_year]))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
401
402
403
404
405
406
407
408
409
410
411
            if (denominator_tot > 0) {
              output_matrix[y, "total"] <- numerator_tot / denominator_tot
            }
          },
          "gla_rel" = {
            numerator_tot <- sum(as.numeric(H_host[, ts_year]))
            denominator_tot <- sum(as.numeric(N_host[, ts_year]))
            if (denominator_tot > 0) {
              output_matrix[y, "total"] <- numerator_tot / denominator_tot
            }
          },
Loup's avatar
Loup committed
412
          { ## i.e. if audpc | gla | substr(type,1,3)=="eco"
Jean-Francois Rey's avatar
Jean-Francois Rey committed
413
            output_matrix[y, "total"] <- sum(output_matrix[y, 1:Nhost], na.rm = TRUE)
414
          }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
415
        )
416
      } ## for y
Jean-Francois Rey's avatar
Jean-Francois Rey committed
417

418
419
      ## Economic outputs
      if (substr(type, 1, 3) == "eco") {
420
421
        eco <- list(eco_yield = output_matrix)
        eco[["eco_product"]] <- data.frame(rep(market_value, each = Nyears) * eco[["eco_yield"]][, 1:Nhost])
422
        ## data.frame to avoid pb when Nhost=1
423
424
425
        colnames(eco[["eco_product"]]) <- cultivar_names ## useful if Nhost=1
        eco[["eco_product"]]$total <- apply(eco[["eco_product"]], 1, sum, na.rm = TRUE)
        eco[["eco_cost"]] <- data.frame(rep(planting_cost_perIndiv, each = Nyears) * t(C_host))
426
427
        colnames(eco[["eco_cost"]]) <- cultivar_names ## useful if Nhost=1
        eco[["eco_cost"]]$total <- apply(eco[["eco_cost"]], 1, sum, na.rm = TRUE)
428
        eco[["eco_margin"]] <- eco[["eco_product"]] - eco[["eco_cost"]]
Jean-Francois Rey's avatar
Jean-Francois Rey committed
429
430

        output_matrix <- eco[[type]] / (area_host * 1E-4) ## normalization by area in Ha
431
      }
432
433
      # product <- t(output_matrix[,1:Nhost])
      # benefit <- market_value %*% product
434
435
      # cost <- planting_cost_perIndiv %*% C_host
      # margin <- benefit - cost
Jean-Francois Rey's avatar
Jean-Francois Rey committed
436

437
      if (writeTXT) {
438
        write.table(output_matrix, file = paste(path, "/", type, ".txt", sep = ""), sep = ",")
439
      }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
440

441
442
443
444
445
446
      if (graphic & Nyears > 1) {
        ## Graphical parameters
        PCH <- 15:(15 + Nhost - 1)
        PCH.tot <- 0
        LTY.tot <- 1
        COL.tot <- RED[1]
Jean-Francois Rey's avatar
Jean-Francois Rey committed
447
448

        switch(type, "audpc" = {
449
          main_output <- "AUDPC"
Loup's avatar
Loup committed
450
451
452
          ylab_output <- expression("Equivalent nb of diseased hosts / time step / m"^2) ## expression('title'^2)
          round_ylab_output <- 2
        }, "audpc_rel" = {
Loup's avatar
Loup committed
453
          main_output <- "Relative AUDPC"
Loup's avatar
Loup committed
454
          ylab_output <- "Proportion of diseased hosts: (I+R)/N"
455
          round_ylab_output <- 2
456
457
458
        }, "gla" = {
          main_output <- expression("Green leaf area (GLA)") ## expression('title'[2])
          ylab_output <- expression("Equivalent nb of healthy hosts / time step / m"^2) ## expression('title'^2)
459
460
461
462
463
          round_ylab_output <- 2
        }, "gla_rel" = {
          main_output <- expression("Relative green leaf area (GLA"[r] * ")")
          ylab_output <- "Proportion of healthy hosts: H/N"
          round_ylab_output <- 2
464
465
        }, "eco_yield" = {
          main_output <- "Crop yield"
Loup Rimbaud's avatar
Loup Rimbaud committed
466
          ylab_output <- "Weight (or volume) units / ha"
467
          round_ylab_output <- 0
468
469
        }, "eco_product" = {
          main_output <- "Crop products"
470
          ylab_output <- "Monetary units / ha"
471
472
          round_ylab_output <- 0
        }, "eco_cost" = {
473
          main_output <- "Operational crop costs"
474
          ylab_output <- "Monetary units / ha"
475
          round_ylab_output <- 0
476
477
        }, "eco_margin" = {
          main_output <- "Margin"
478
          ylab_output <- "Monetary units / ha"
479
480
481
482
          round_ylab_output <- 0
        })
        ymin_output <- ylim_param[[type]][1]
        ymax_output <- ylim_param[[type]][2]
Jean-Francois Rey's avatar
Jean-Francois Rey committed
483
        if (is.na(ymin_output)) {
484
485
          ymin_output <- floor(min(0, min(output_matrix, na.rm = TRUE)))
        }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
486
        if (is.na(ymax_output)) {
487
488
489
          ymax_output <- ceiling(max(output_matrix, na.rm = TRUE))
        }
        output_mean_host <- apply(output_matrix, 2, mean, na.rm = TRUE) ## Average output per cultivar
Jean-Francois Rey's avatar
Jean-Francois Rey committed
490

491
        graphics.off()
Jean-Francois Rey's avatar
Jean-Francois Rey committed
492
493
494
495
        tiff(
          filename = paste(path, "/", type, ".tiff", sep = ""),
          width = 180, height = 110, units = "mm", compression = "lzw", res = 300
        )
496
497
498
499
        # m <- matrix(c(rep(1,5),2),6,1)    # matrix(c(rep(1,5),3,rep(2,5),3),6,2)
        # layout(m)
        # par(xpd=F, cex=.9,mar=c(5,4.1,4.1,2.1))
        par(xpd = NA, cex = .9, mar = c(4, 4.3, 3, 9))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
500
501
502
503
        plot(0, 0,
          type = "n", bty = "n", xaxt = "n", yaxt = "n", xlim = c(1, Nyears), ylim = c(ymin_output, ymax_output),
          main = main_output, xlab = "Years", ylab = ylab_output
        )
504
505
506
507
508
509
510
511
        axis(1, at = seq(1, Nyears, by = (Nyears - 1) %/% 10 + 1))
        axis(2, at = round(seq(ymin_output, ymax_output, length.out = 5), round_ylab_output), las = 1)
        for (host in 1:Nhost) {
          lines(1:Nyears, output_matrix[, host], type = "o", lwd = 2, lty = host, pch = PCH[host], col = GRAY[host])
          points(par("usr")[1], output_mean_host[host], pch = PCH[host], col = GRAY[host], xpd = TRUE)
        }
        lines(1:Nyears, output_matrix[, "total"], type = "o", lwd = 2, lty = LTY.tot, pch = PCH.tot, col = COL.tot)
        points(par("usr")[1], output_mean_host["total"], pch = PCH.tot, col = COL.tot, cex = 1, xpd = TRUE)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
512

513
        if (ymin_output < 0) {
Loup Rimbaud's avatar
Loup Rimbaud committed
514
          abline(h = 0, lwd = 1, lty = 3, col = BLUE[1], xpd = FALSE)
515
        }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
516
517
518
519
520
521

        legend(Nyears * 1.05, .5 * ymax_output,
          title.adj = -.01, bty = "n", title = "Cultivar:",
          legend = c(cultivar_names, "Whole landscape"), cex = 0.9, lty = c(1:Nhost, LTY.tot),
          lwd = 2, pch = c(PCH, PCH.tot), pt.cex = c(rep(1, Nhost), 1), col = c(GRAY[1:Nhost], COL.tot), seg.len = 2.5
        )
522
        dev.off()
523
      }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
524

525
526
527
      res[[type]] <- output_matrix
    } ## for type
  } ## else if output != HLIR_dynamics
Jean-Francois Rey's avatar
Jean-Francois Rey committed
528

529
530
  return(res)
}
Jean-Francois Rey's avatar
Jean-Francois Rey committed
531
532


Jean-Francois Rey's avatar
Jean-Francois Rey committed
533
#' @title Switch from index of genotype to indices of agressiveness on different components
Loup Rimbaud's avatar
Loup Rimbaud committed
534
#' @name switch_patho_to_aggr
Loup Rimbaud's avatar
Loup Rimbaud committed
535
#' @description Finds the level of aggressiveness on different components (targeted by different resistance genes)
Loup Rimbaud's avatar
Loup Rimbaud committed
536
537
538
539
540
#' from the index of a given pathogen genotype
#' @param index_patho index of pathogen genotype
#' @param Ngenes number of resistance genes
#' @param Nlevels_aggressiveness vector of the number of adaptation levels related to each resistance gene
#' @return a vector containing the indices of aggressiveness on the different components targeter by the resistance genes
Jean-Francois Rey's avatar
Jean-Francois Rey committed
541
542
#' @examples
#' switch_patho_to_aggr(5, 3, c(2, 2, 3))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
543
#' @export
544
switch_patho_to_aggr <- function(index_patho, Ngenes, Nlevels_aggressiveness) {
545
546
547
548
549
550
551
552
553
554
  ## Index_patho starts at 0
  ## Index aggressiveness starts at 0
  aggr <- numeric()
  remainder <- index_patho
  for (g in 1:Ngenes) {
    prod <- 1
    k <- g
    while (k < Ngenes) {
      k <- k + 1
      prod <- prod * Nlevels_aggressiveness[k]
Jean-Francois Rey's avatar
Jean-Francois Rey committed
555
    }
556
557
    aggr[g] <- remainder %/% prod ## Quotient
    remainder <- remainder %% prod
Jean-Francois Rey's avatar
Jean-Francois Rey committed
558
  }
559
560
  return(aggr)
}
Jean-Francois Rey's avatar
Jean-Francois Rey committed
561
562
563
564
565



#' @title Generation of evolutionary model outputs
#' @name evol_output
Loup Rimbaud's avatar
Loup Rimbaud committed
566
#' @description Generates evolutionary outputs from model simulations.
567
568
#' @param types a character string (or a vector of character strings if several outputs are to be computed) specifying the
#' type of outputs to generate (see details):\itemize{
Loup Rimbaud's avatar
Loup Rimbaud committed
569
570
#' \item "evol_patho": Dynamics of pathogen genotype frequencies
#' \item "evol_aggr": Evolution of pathogen aggressiveness
571
572
573
#' \item "durability": Durability of resistance genes
#' \item "all": compute all these outputs (default)
#' }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
574
#' @param time_param list of simulation parameters:\itemize{
Loup Rimbaud's avatar
Loup Rimbaud committed
575
576
577
578
#' \item Nyears = number cropping seasons,
#' \item nTSpY = number of time-steps per cropping season.
#' }
#' @param Npoly number of fields in the landscape.
Jean-Francois Rey's avatar
Jean-Francois Rey committed
579
580
581
#' @param cultivars_param list of parameters associated with each host genotype (i.e. cultivars)
#' when cultivated in pure crops:\itemize{
#' \item name = vector of cultivar names,
Loup Rimbaud's avatar
Loup Rimbaud committed
582
583
#' \item cultivars_genes_list = a list containing, for each host genotype, the indices of carried resistance genes.
#' }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
584
585
586
#' @param genes_param list of parameters associated with each resistance gene and with the evolution of
#' each corresponding pathogenicity gene:\itemize{
#' \item name = vector of names of resistance genes,
Loup Rimbaud's avatar
Loup Rimbaud committed
587
#' \item Nlevels_aggressiveness = vector containing the number of adaptation levels related to each resistance gene (i.e. 1 + number
Jean-Francois Rey's avatar
Jean-Francois Rey committed
588
#' of required mutations for a pathogenicity gene to fully adapt to the corresponding resistance gene),
Loup Rimbaud's avatar
Loup Rimbaud committed
589
#' }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
590
591
#' @param thres_breakdown an integer (or vector of integers) giving the threshold (i.e. number of infections) above which a
#' pathogen genotype is unlikely to go extinct and resistance is considered broken down, used to characterise the time to
Loup Rimbaud's avatar
Loup Rimbaud committed
592
#' invasion of resistant hosts (several values are computed if several thresholds are given in a vector).
593
#' @param writeTXT a logical indicating if the output is written in a text file (TRUE) or not (FALSE).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
594
#' @param graphic a logical indicating if graphics must be generated (TRUE) or not (FALSE).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
595
#' @param path a character string indicating the path of the repository where simulation output files are located and
Loup Rimbaud's avatar
Loup Rimbaud committed
596
#' where .txt files and graphics will be generated.
Jean-Francois Rey's avatar
Jean-Francois Rey committed
597
#'
Jean-Francois Rey's avatar
Jean-Francois Rey committed
598
#' @details For each pathogen genotype, several computations are performed: \itemize{
Loup Rimbaud's avatar
Loup Rimbaud committed
599
600
601
#' \item appearance: time to first appearance (as propagule);
#' \item R_infection: time to first true infection of a resistant host;
#' \item R_invasion: time when the number of infections of resistant hosts reaches a threshold above which
Jean-Francois Rey's avatar
Jean-Francois Rey committed
602
#' the genotype is unlikely to go extinct.}
603
#' The value Nyears + 1 time step is used if the genotype never appeared/infected/invaded.
Jean-Francois Rey's avatar
Jean-Francois Rey committed
604
#' @return A list containing, for each required type of output, a matrix summarising the output.
605
#' Each matrix can be written in a txt file (if writeTXT=TRUE), and illustrated in a graphic (if graphic=TRUE).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
606
#' @references Rimbaud L., Papaïx J., Rey J.-F., Barrett L. G. and Thrall P. H. (2018). Assessing the durability and efficiency
Loup Rimbaud's avatar
Loup Rimbaud committed
607
#' of landscape-based strategies to deploy plant resistance to pathogens. \emph{PLoS Computational Biology} 14(4):e1006067.
Loup Rimbaud's avatar
Loup Rimbaud committed
608
#' @seealso \link{epid_output}
Jean-Francois Rey's avatar
Jean-Francois Rey committed
609
#' @examples
610
#' \dontrun{
Jean-Francois Rey's avatar
Jean-Francois Rey committed
611
612
#' demo_landsepi()
#' }
613
#' @include RcppExports.R Math-Functions.R graphics.R
Jean-Francois Rey's avatar
Jean-Francois Rey committed
614
615
616
#' @importFrom grDevices colorRampPalette dev.off graphics.off gray png tiff
#' @importFrom utils write.table
#' @export
Jean-Francois Rey's avatar
Jean-Francois Rey committed
617
evol_output <- function(types = "all", time_param, Npoly, cultivars_param, genes_param, thres_breakdown = 50000,
618
                        writeTXT = TRUE, graphic = TRUE, path = getwd()) {
Loup Rimbaud's avatar
Loup Rimbaud committed
619

Loup Rimbaud's avatar
Loup Rimbaud committed
620
  valid_outputs <- c("evol_patho", "evol_aggr", "durability")
Loup Rimbaud's avatar
Loup Rimbaud committed
621
  if (types[1] == "all") {
622
623
624
    types <- valid_outputs
  }
  if (is.na(sum(match(types, valid_outputs)))) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
625
    stop(paste("Error: valid types of outputs are", paste(valid_outputs, collapse = ", ")))
626
  }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
627

Jean-Francois Rey's avatar
Jean-Francois Rey committed
628
629
630
631
  ## Time parameters
  Nyears <- time_param$Nyears
  nTSpY <- time_param$nTSpY
  nTS <- Nyears * nTSpY ## Total number of time-steps
Jean-Francois Rey's avatar
Jean-Francois Rey committed
632

Jean-Francois Rey's avatar
Jean-Francois Rey committed
633
  ## Host parameters
634
  Nhost <- length(cultivars_param$name)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
635
  cultivars_genes_list <- cultivars_param$cultivars_genes_list
Jean-Francois Rey's avatar
Jean-Francois Rey committed
636

Jean-Francois Rey's avatar
Jean-Francois Rey committed
637
  # find the indices of resistant hosts
Jean-Francois Rey's avatar
Jean-Francois Rey committed
638
639
640
  Rhosts <- (1:Nhost)[as.logical(lapply(cultivars_genes_list, function(x) {
    length(x) > 0
  }))]
Jean-Francois Rey's avatar
Jean-Francois Rey committed
641

Jean-Francois Rey's avatar
Jean-Francois Rey committed
642
  ## Gene parameters
Loup Rimbaud's avatar
Loup Rimbaud committed
643
  Nlevels_aggressiveness <- genes_param$Nlevels_aggressiveness
Loup Rimbaud's avatar
Loup Rimbaud committed
644
  Ngenes <- length(Nlevels_aggressiveness)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
645
  Npatho <- prod(Nlevels_aggressiveness)
Loup Rimbaud's avatar
Loup Rimbaud committed
646
  
Jean-Francois Rey's avatar
Jean-Francois Rey committed
647

Jean-Francois Rey's avatar
Jean-Francois Rey committed
648
649
650
651
652
  #### IMPORTATION OF THE SIMULATION OUTPUT
  P <- as.list(1:nTS)
  L <- as.list(1:nTS)
  I <- as.list(1:nTS)
  index <- 0
Jean-Francois Rey's avatar
Jean-Francois Rey committed
653

Loup Rimbaud's avatar
Loup Rimbaud committed
654
  # print("Reading binary files to compute evolutionary model outputs")
Jean-Francois Rey's avatar
Jean-Francois Rey committed
655
  for (year in 1:Nyears) {
656
657
658
659
660
    binfileP <- file(paste(path, sprintf("/P-%02d", year), ".bin", sep = ""), "rb")
    P.tmp <- readBin(con = binfileP, what = "int", n = Npoly * Npatho * nTSpY, size = 4, signed = T, endian = "little")
    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")
    binfileL <- file(paste(path, sprintf("/L-%02d", year), ".bin", sep = ""), "rb")
661
    L.tmp <- readBin(con = binfileL, what = "int", n = Npoly * Npatho * Nhost * nTSpY, size = 4, signed = T, endian = "little")
Jean-Francois Rey's avatar
Jean-Francois Rey committed
662
663


Jean-Francois Rey's avatar
Jean-Francois Rey committed
664
    for (t in 1:nTSpY) {
665
      P[[t + index]] <- matrix(P.tmp[((Npatho * Npoly) * (t - 1) + 1):(t * (Npatho * Npoly))], ncol = Npatho, byrow = T)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
666
667
668
669
670
671
672
673
      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)
      )
Jean-Francois Rey's avatar
Jean-Francois Rey committed
674
    }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
675

Jean-Francois Rey's avatar
Jean-Francois Rey committed
676
677
678
679
680
    index <- index + nTSpY
    close(binfileP)
    close(binfileL)
    close(binfileI)
  }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
681

Jean-Francois Rey's avatar
Jean-Francois Rey committed
682
683
684
  ####
  ####          COMPUTATION OF OUTPUTS
  ####
Loup Rimbaud's avatar
Loup Rimbaud committed
685
  res <- list()
Jean-Francois Rey's avatar
Jean-Francois Rey committed
686

Jean-Francois Rey's avatar
Jean-Francois Rey committed
687
688
  ## Matrix to switch from patho to aggr
  pathoToAggr <- matrix(NA, nrow = Npatho, ncol = Ngenes)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
689
  for (p in 1:Npatho) {
690
    pathoToAggr[p, ] <- switch_patho_to_aggr(p - 1, Ngenes, Nlevels_aggressiveness) + 1
Jean-Francois Rey's avatar
Jean-Francois Rey committed
691
692
  } ## -1/+1 because indices start at 0 in function switch

Jean-Francois Rey's avatar
Jean-Francois Rey committed
693
694
695
696
697
698
699
700
701
  ####  Summary of P, L, I for each pathogen genotype
  P_patho <- NULL
  L_patho <- NULL
  I_patho <- NULL
  IL_patho_Rhost <- NULL ## sum of L and I on resistant hosts
  for (t in 1:nTS) {
    P_patho <- cbind(P_patho, apply(P[[t]], 2, sum))
    L_patho <- cbind(L_patho, apply(L[[t]], 2, sum))
    I_patho <- cbind(I_patho, apply(I[[t]], 2, sum))
702
703
    if (length(Rhosts) > 0 & Npatho > 1) {
      IL_patho_Rhost <- cbind(IL_patho_Rhost, apply(L[[t]][Rhosts, , ] + I[[t]][Rhosts, , ], 1 + (length(Rhosts) > 1), sum))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
704
705
    } ## length(Rhosts)=2 --> need to use the 2nd dimension
  }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
706

Jean-Francois Rey's avatar
Jean-Francois Rey committed
707
708
  ## Genotype frequencies
  I_pathoProp <- matrix(NA, nrow = Npatho, ncol = nTS)
709
  Itot <- apply(matrix(I_patho, ncol = nTS), 2, sum) ## matrix to avoid pb if only 1 row
Jean-Francois Rey's avatar
Jean-Francois Rey committed
710
711
  for (p in 1:Npatho) {
    for (t in 1:nTS) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
712
      if (Itot[t] > 0) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
713
        I_pathoProp[p, t] <- I_patho[p, t] / Itot[t]
Jean-Francois Rey's avatar
Jean-Francois Rey committed
714
      }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
715
716
    }
  }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
717

Jean-Francois Rey's avatar
Jean-Francois Rey committed
718
  ## Computed variables
Loup Rimbaud's avatar
Loup Rimbaud committed
719
720
  if (length(thres_breakdown) == 1) {
    names(thres_breakdown) <- "R_invasion"
Jean-Francois Rey's avatar
Jean-Francois Rey committed
721
  } else {
Loup Rimbaud's avatar
Loup Rimbaud committed
722
    names(thres_breakdown) <- paste("R_invasion_", 1:length(thres_breakdown), sep = "")
Jean-Francois Rey's avatar
Jean-Francois Rey committed
723
  }
Loup Rimbaud's avatar
Loup Rimbaud committed
724
  output_variables <- c("appearance", "R_infection", names(thres_breakdown))
Jean-Francois Rey's avatar
Jean-Francois Rey committed
725

726
  ## "time_to_first" variables for each pathogen (must be computed in any case)
Loup Rimbaud's avatar
Loup Rimbaud committed
727
  patho_evol <- matrix(NA, ncol = 3 + length(thres_breakdown), nrow = Npatho)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
728
  colnames(patho_evol) <- c(output_variables, "finalPopSize_LI")
Jean-Francois Rey's avatar
Jean-Francois Rey committed
729
  rownames(patho_evol) <- paste("patho_", 1:Npatho, sep = "")
Loup Rimbaud's avatar
Loup Rimbaud committed
730
  colnames(pathoToAggr) <- genes_param$name
Jean-Francois Rey's avatar
Jean-Francois Rey committed
731
  rownames(pathoToAggr) <- paste("patho_", 1:Npatho, sep = "")
Jean-Francois Rey's avatar
Jean-Francois Rey committed
732

Jean-Francois Rey's avatar
Jean-Francois Rey committed
733
  for (patho in 1:Npatho) {
734
735
736
    patho_evol[patho, "appearance"] <- min(which(P_patho[patho, ] > 0), nTS + 1, na.rm = TRUE)
    patho_evol[patho, "R_infection"] <- min(which(IL_patho_Rhost[patho, ] > 0), nTS + 1, na.rm = TRUE)
    patho_evol[patho, "finalPopSize_LI"] <- as.numeric(sum(L_patho[patho, nTS], I_patho[patho, nTS]))
Loup Rimbaud's avatar
Loup Rimbaud committed
737
    for (k in 1:length(thres_breakdown)) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
738
739
740
741
      patho_evol[patho, names(thres_breakdown)[k]] <- min(which(IL_patho_Rhost[patho, ] > thres_breakdown[k]),
        nTS + 1,
        na.rm = TRUE
      )
Jean-Francois Rey's avatar
Jean-Francois Rey committed
742
743
    }
  }
744
  if (writeTXT & sum(is.element(types, c("evol_patho"))) > 0) {
745
    write.table(patho_evol, file = paste(path, "/evol_patho", ".txt", sep = ""), sep = ",")
Loup Rimbaud's avatar
Loup Rimbaud committed
746
    write.table(pathoToAggr, file = paste(path, "/evol_patho2aggr", ".txt", sep = ""), sep = ",")
Jean-Francois Rey's avatar
Jean-Francois Rey committed
747
  }
Loup Rimbaud's avatar
Loup Rimbaud committed
748
  res[["evol_patho"]] <- patho_evol
Jean-Francois Rey's avatar
Jean-Francois Rey committed
749

750
  ## "Time_to_first" variables for each level of aggressiveness
Loup Rimbaud's avatar
Loup Rimbaud committed
751
  if (sum(is.element(types, c("evol_aggr", "durability"))) > 0) {
752
753
754
755
756
757
758
759
760
761
762
    aggr_evol <- list()
    for (g in 1:Ngenes) {
      aggr_evol[[g]] <- matrix(NA, nrow = Nlevels_aggressiveness[g], ncol = length(output_variables))
      colnames(aggr_evol[[g]]) <- output_variables
      rownames(aggr_evol[[g]]) <- paste("level_", 1:Nlevels_aggressiveness[g], sep = "")
      for (a in 1:nrow(aggr_evol[[g]])) {
        for (k in colnames(aggr_evol[[g]])) {
          index_patho <- pathoToAggr[, g] == a
          aggr_evol[[g]][a, k] <- min(patho_evol[index_patho, k])
        } ## for k
      } ## for a
763
      if (writeTXT & sum(is.element(types, c("evol_aggr"))) > 0) {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
764
        write.table(aggr_evol[[g]], file = paste(path, "/evol_aggr_", genes_param$name[g], ".txt", sep = ""), sep = ",")
Loup Rimbaud's avatar
Loup Rimbaud committed
765
      }
766
    } ## for g
767
    names(aggr_evol) <- genes_param$name
Loup Rimbaud's avatar
Loup Rimbaud committed
768
    res[["aggr_evol"]] <- aggr_evol
Jean-Francois Rey's avatar
Jean-Francois Rey committed
769

770
771
772
    ## Durability relative to a given criterion
    if (sum(is.element(types, c("durability"))) > 0) {
      durability <- matrix(NA, nrow = 1, ncol = Ngenes)
Loup Rimbaud's avatar
Loup Rimbaud committed
773
      colnames(durability) <- names(aggr_evol)
774
775
776
777
778
      for (g in 1:Ngenes) {
        finalLevel <- nrow(aggr_evol[[g]]) ## last line
        criterion <- ncol(aggr_evol[[g]]) ## last column
        durability[, g] <- aggr_evol[[g]][finalLevel, criterion]
      }
779
      if (writeTXT & sum(is.element(types, c("durability"))) > 0) {
780
781
        write.table(durability, file = paste(path, "/durability", ".txt", sep = ""), sep = ",")
      }
Loup Rimbaud's avatar
Loup Rimbaud committed
782
      res[["durability"]] <- durability
783
784
    } ## if durability
  } ## if aggr_evol
Jean-Francois Rey's avatar
Jean-Francois Rey committed
785
786


Jean-Francois Rey's avatar
Jean-Francois Rey committed
787
788
789
790
791
  ####
  ####               GRAPHICS
  ####
  if (graphic) {
    graphics.off()
Jean-Francois Rey's avatar
Jean-Francois Rey committed
792

Jean-Francois Rey's avatar
Jean-Francois Rey committed
793
794
795
    ## Dynamics of pathotype frequencies (gray scales)
    I_aggrProp <- list()
    for (g in 1:Ngenes) {
796
      I_aggrProp[[g]] <- matrix(0, nrow = Nlevels_aggressiveness[g], ncol = nTS)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
797
798
799
      for (t in 1:nTS) {
        for (p in 1:Npatho) {
          aggr <- pathoToAggr[p, g]
800
          I_aggrProp[[g]][aggr, t] <- I_aggrProp[[g]][aggr, t] + I_pathoProp[p, t]
Jean-Francois Rey's avatar
Jean-Francois Rey committed
801
802
        }
      }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
803

Jean-Francois Rey's avatar
Jean-Francois Rey committed
804
805
806
807
808
      if (Nlevels_aggressiveness[g] > 1) {
        tiff(
          filename = paste(path, "/freqPatho_gene", g, ".tiff", sep = ""), width = 100, height = 100,
          units = "mm", compression = "lzw", res = 300
        )
809
810
811
812
        par(xpd = F, mar = c(4, 4, 2, 2))
        plot_freqPatho(genes_param$name[g], Nlevels_aggressiveness[g], I_aggrProp[[g]], nTS, Nyears, nTSpY)
        dev.off()
      }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
813
    }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
814

Jean-Francois Rey's avatar
Jean-Francois Rey committed
815
816
    ## Dynamics of genotypes frequencies (curves)
    # COL <- c("#FF5555","#4F94CD","darkolivegreen4","#CD950C","black") ## colors: red, blue, green, orange, black
Jean-Francois Rey's avatar
Jean-Francois Rey committed
817
818
819
820
    tiff(
      filename = paste(path, "/freqPathoGenotypes.tiff", sep = ""), width = 180, height = 110,
      units = "mm", compression = "lzw", res = 300
    )
Jean-Francois Rey's avatar
Jean-Francois Rey committed
821
    par(xpd = FALSE, mar = c(4, 4, 0, 9))
822
    plot(0, 0, type = "n", xlim = c(1, nTS), bty = "n", las = 1, ylim = c(0, 1), xaxt = "n", xlab = "", ylab = "Frequency")
Jean-Francois Rey's avatar
Jean-Francois Rey committed
823
824
825
826
    if (Nyears == 1) {
      axis(1, at = round(seq(1, nTS, length.out = 8)), las = 1)
      title(xlab = "Evolutionnary time (days)")
    } else {
Jean-Francois Rey's avatar
Jean-Francois Rey committed
827
828
829
830
      axis(1,
        las = 1, at = seq(1, nTS + 1, nTSpY * ((Nyears - 1) %/% 10 + 1)),
        labels = seq(0, Nyears, ((Nyears - 1) %/% 10 + 1))
      )
Jean-Francois Rey's avatar
Jean-Francois Rey committed
831
832
      title(xlab = "Evolutionnary time (years)")
    }
Jean-Francois Rey's avatar
Jean-Francois Rey committed
833
834

    for (patho in 1:Npatho) {
835
      lines(I_pathoProp[patho, ], col = patho, lty = patho, lwd = 1.5)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
836
837
    }

Jean-Francois Rey's avatar
Jean-Francois Rey committed
838
    par(xpd = TRUE)
Jean-Francois Rey's avatar
Jean-Francois Rey committed
839
840
841
842
    legend(nTS * 1.05, .5,
      title.adj = -.01, bty = "n", title = "Pathogen genotype:",
      legend = 1:Npatho, lty = 1:Npatho, lwd = 2, col = 1:Npatho
    )
Jean-Francois Rey's avatar
Jean-Francois Rey committed
843
844
    dev.off()
  } ## if graphic
Jean-Francois Rey's avatar
Jean-Francois Rey committed
845

Loup Rimbaud's avatar
Loup Rimbaud committed
846
  return(res[match(types, valid_outputs)])
Jean-Francois Rey's avatar
Jean-Francois Rey committed
847
848
849
850
851
}


#' @title Generation of a video
#' @name video
Loup Rimbaud's avatar
Loup Rimbaud committed
852
#' @description Generates a video showing the epidemic dynamics on a map representing the cropping landscape.
Loup Rimbaud's avatar
Loup Rimbaud committed
853
#' (requires ffmpeg library).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
854
855
#' @param audpc A dataframe containing audpc outputs (generated through epid_output). 1 year per line and
#' 1 column per cultivar, with an additional column for the average audpc in the landscape.
Jean-Francois Rey's avatar
Jean-Francois Rey committed
856
#' @param time_param list of simulation parameters:\itemize{
Loup Rimbaud's avatar
Loup Rimbaud committed
857
858
859
860
#' \item Nyears = number cropping seasons,
#' \item nTSpY = number of time-steps per cropping season.
#' }
#' @param Npatho number of pathogen genotypes.
Loup Rimbaud's avatar
Loup Rimbaud committed
861
862
#' @param landscape a sp object containing the landscape.
#' @param area a vector containing polygon areas (must be in square meters).
Jean-Francois Rey's avatar
Jean-Francois Rey committed
863
864
865
#' @param rotation a dataframe containing for each field (rows) and year (columns, named "year_1", "year_2", etc.),
#' the index of the cultivated croptype. Importantly, the matrix must contain 1 more column than the real number
#' of simulated years.
Loup Rimbaud's avatar
Loup Rimbaud committed
866
867
#' @param croptypes a dataframe with three columns named 'croptypeID' for croptype index,
#' 'cultivarID' for cultivar index and 'proportion' for the proportion of the cultivar within the croptype.
868
#' @param croptype_names a vector of croptype names (for legend).