CreateInputsCrit.R 12.2 KB
Newer Older
1
2
3
CreateInputsCrit <- function(FUN_CRIT,
                             InputsModel,
                             RunOptions,
4
                             Qobs,              # deprecated
5
                             Obs,
6
                             VarObs = "Q",
7
8
                             BoolCrit = NULL,
                             transfo = "",
9
                             weights = NULL,
10
                             Ind_zeroes = NULL, # deprecated
11
                             epsilon = NULL,
12
                             warnings = TRUE,
13
14
15
16
17
                             verbose = TRUE) {
  
  
  ObjectClass <- NULL
  
18
  
19
  ## ---------- check arguments
20
  
21
  if (!missing(Qobs)) {
22
    if (missing(Obs)) {
23
      if (warnings) {
24
        warning("argument 'Qobs' is deprecated. Please use 'Obs' and 'VarObs' instead")
25
      }
26
      Obs <- Qobs
27
      # VarObs <- "Qobs"
28
    } else {
29
      warning("argument 'Qobs' is deprecated. The values set in 'Obs' will be used instead")
30
    }
31
32
  }
  if (!missing(Ind_zeroes) & warnings) {
33
    warning("deprecated 'Ind_zeroes' argument")
34
35
  }
  if (!missing(verbose)) {
36
    warning("deprecated 'verbose' argument. Use 'warnings', instead")
37
38
39
  }
  
  
40
41
  ## check 'InputsModel'
  if (!inherits(InputsModel, "InputsModel")) {
42
    stop("'InputsModel' must be of class 'InputsModel'")
43
44
45
46
47
48
49
  }
  
  
  ## length of index of period to be used for the model run
  LLL <- length(InputsModel$DatesR[RunOptions$IndPeriod_Run])
  
  
50
51
  ## check 'Obs' and definition of idLayer
  vecObs <- unlist(Obs)
52
  if (length(vecObs) %% LLL != 0 | !is.numeric(vecObs)) {
53
    stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
54
  }
55
  if (!is.list(Obs)) {
56
    idLayer <- list(1L)
57
    Obs <- list(Obs)
58
  } else {
59
    idLayer <- lapply(Obs, function(i) {
60
61
62
63
64
65
        if (is.list(i)) {
          length(i)
        } else {
          1L
        }
      })
66
    Obs <- lapply(Obs, function(x) rowMeans(as.data.frame(x)))
67
  }
68

69
  
70
71
  ## create list of arguments
  listArgs <- list(FUN_CRIT   = FUN_CRIT,
72
                   Obs        = Obs,
73
                   VarObs     = VarObs,
74
                   BoolCrit   = BoolCrit,
75
                   idLayer    = idLayer,
76
77
78
79
80
81
82
83
84
85
86
                   transfo    = transfo,
                   weights    = weights,
                   epsilon    = epsilon)
  
  
  ## check lists lengths
  for (iArgs in names(listArgs)) {
    if (iArgs %in% c("weights", "BoolCrit", "epsilon")) {
      if (any(is.null(listArgs[[iArgs]]))) {
        listArgs[[iArgs]] <- lapply(seq_along(listArgs$FUN_CRIT), function(x) NULL)
      }
87
    }
88
    if (iArgs %in% c("FUN_CRIT", "VarObs", "transfo", "weights") & length(listArgs[[iArgs]]) > 1L) {
89
90
      listArgs[[iArgs]] <- as.list(listArgs[[iArgs]])
    }
91
92
    if (!is.list(listArgs[[iArgs]])) {
      listArgs[[iArgs]] <- list(listArgs[[iArgs]])
93
    }
94
95
  }
  
96
97
98
99
  ## check 'FUN_CRIT'
  listArgs$FUN_CRIT <- lapply(listArgs$FUN_CRIT, FUN = match.fun)
  
  
100
101
102
  ## check 'VarObs'
  if (missing(VarObs)) {
    listArgs$VarObs <- as.list(rep("Q", times = length(listArgs$Obs)))
103
    # if (warnings) {
104
    #   warning("'VarObs' automatically set to \"Q\"")
105
    # }
106
107
  }
  
108
  
109
110
111
  ## check 'VarObs' + 'RunOptions'
  if ("Q" %in% VarObs & !inherits(RunOptions, "GR")) {
    stop("'VarObs' cannot contain Q if a GR rainfall-runoff model is not used")
112
  }
113
114
  if (any(c("SCA", "SWE") %in% VarObs) & !inherits(RunOptions, "CemaNeige")) {
    stop("'VarObs' cannot contain SCA or SWE if CemaNeige is not used")
115
  }
116
  if ("SCA" %in% VarObs & inherits(RunOptions, "CemaNeige") & !"Gratio"   %in% RunOptions$Outputs_Sim) {
117
118
    stop("'Gratio' is missing in 'Outputs_Sim' of 'RunOptions', which is necessary to output SCA with CemaNeige")
  }
119
  if ("SWE" %in% VarObs & inherits(RunOptions, "CemaNeige") & !"SnowPack" %in% RunOptions$Outputs_Sim) {
120
121
122
123
    stop("'SnowPack' is missing in 'Outputs_Sim' of 'RunOptions', which is necessary to output SWE with CemaNeige")
  }

  
124
125
  ## check 'transfo'
  if (missing(transfo)) {
126
    listArgs$transfo <- as.list(rep("", times = length(listArgs$Obs)))
127
128
129
    # if (warnings) {
    #   warning("'transfo' automatically set to \"\"")
    # }
130
131
132
133
134
135
136
137
138
139
140
  }  
  
  ## check length of each args
  if (length(unique(sapply(listArgs, FUN = length))) != 1) {
    stopListArgs <- paste(sapply(names(listArgs), shQuote), collapse = ", ")
    stop(sprintf("arguments %s must have the same length", stopListArgs))
  }
  
  
  ## check 'RunOptions'
  if (!inherits(RunOptions , "RunOptions")) {
141
    stop("'RunOptions' must be of class 'RunOptions'")
142
143
  }
  
144
  
145
146
  ## check 'weights'
  if (length(listArgs$weights) > 1 & sum(unlist(listArgs$weights)) == 0 & !any(sapply(listArgs$weights, is.null))) {
147
    stop("sum of 'weights' cannot be equal to zero")
148
  }
149

150
  
151
152
  ## ---------- reformat
  
153
154
155
156
  ## reformat list of arguments
  listArgs2 <- lapply(seq_along(listArgs$FUN_CRIT), function(i) lapply(listArgs, "[[", i))
  
  ## preparation of warning messages
157
  inVarObs  <- c("Q", "SCA", "SWE")
158
  msgVarObs <- "'VarObs' must be a (list of) character vector(s) and one of %s"
159
  msgVarObs <- sprintf(msgVarObs, paste(sapply(inVarObs, shQuote), collapse = ", "))
160
  inTransfo  <- c("", "sqrt", "log", "inv", "sort")
161
  msgTransfo <- "'transfo' must be a (list of) character vector(s) and one of %s"
162
  msgTransfo <- sprintf(msgTransfo, paste(sapply(inTransfo, shQuote), collapse = ", "))
163
164
165
166
167
168
169
170
171
  
  
  ## ---------- loop on the list of inputs
  
  InputsCrit <- lapply(listArgs2, function(iListArgs2) {
    
    ## check 'FUN_CRIT'
    if (!(identical(iListArgs2$FUN_CRIT, ErrorCrit_NSE ) | identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE ) |
          identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE2) | identical(iListArgs2$FUN_CRIT, ErrorCrit_RMSE))) {
172
      stop("incorrect 'FUN_CRIT' for use in 'CreateInputsCrit'", call. = FALSE)
173
    }
174
    if (identical(iListArgs2$FUN_CRIT, ErrorCrit_RMSE) & length(listArgs$weights) > 1 & all(!is.null(unlist(listArgs$weights)))) {
175
      stop("calculating a composite criterion with the RMSE is not allowed since RMSE is not a dimensionless metric", call. = FALSE)
176
    }
177
    
178
179
180
    ## check 'Obs'
    if (!is.vector(iListArgs2$Obs) | length(iListArgs2$Obs) != LLL | !is.numeric(iListArgs2$Obs)) {
      stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
181
182
    }
    
183
184
    ## check 'BoolCrit'
    if (is.null(iListArgs2$BoolCrit)) {
185
      iListArgs2$BoolCrit <- rep(TRUE, length(iListArgs2$Obs))
186
    }
187
    if (!is.logical(iListArgs2$BoolCrit)) {
188
      stop("'BoolCrit' must be a (list of) vector(s) of boolean", call. = FALSE)
189
    }
190
    if (length(iListArgs2$BoolCrit) != LLL) {
191
      stop("'BoolCrit' and 'InputsModel' series must have the same length", call. = FALSE)
192
    }
193
    
194
195
    ## check 'VarObs'
    if (!is.vector(iListArgs2$VarObs) | length(iListArgs2$VarObs) != 1 | !is.character(iListArgs2$VarObs) | !all(iListArgs2$VarObs %in% inVarObs)) {
196
      stop(msgVarObs, call. = FALSE)
197
    }
198
    
199
200
201
    ## check 'VarObs' + 'Obs'
    if (any(iListArgs2$VarObs %in% "SCA")) {
      idSCA <- which(iListArgs2$VarObs == "SCA")
202
      if (length(idSCA) == 1L) {
203
        vecSCA <- iListArgs2$Obs
204
      } else {
205
        vecSCA <- unlist(iListArgs2$Obs[idSCA])
206
      }
207
      if (min(vecSCA, na.rm = TRUE) < 0 | max(vecSCA, na.rm = TRUE) > 1) {
208
        stop("'Obs' outside [0,1] for \"SCA\"", call. = FALSE)
209
210
      }
    } 
211
    inPosVarObs <- c("Q", "SWE")
212
213
    if (any(iListArgs2$VarObs %in% inPosVarObs)) {
      idQSS <- which(iListArgs2$VarObs %in% inPosVarObs)
214
      if (length(idQSS) == 1L) {
215
        vecQSS <- iListArgs2$Obs
216
      } else {
217
        vecQSS <- unlist(iListArgs2$Obs[idQSS])
218
      }
219
      if (min(vecQSS, na.rm = TRUE) < 0) {
220
        stop(sprintf("'Obs' outside [0,Inf[ for \"%s\"", iListArgs2$VarObs), call. = FALSE)
221
222
223
224
      }
    }
    
    
225
226
227
    ## check 'transfo'
    if (is.null(iListArgs2$transfo) | !is.vector(iListArgs2$transfo) | length(iListArgs2$transfo) != 1 | !is.character(iListArgs2$transfo) | !all(iListArgs2$transfo %in% inTransfo)) {
      stop(msgTransfo, call. = FALSE)
228
229
    }
    
230
231
232
    ## check 'weights'
    if (!is.null(iListArgs2$weights)) {
      if (!is.vector(iListArgs2$weights) | length(iListArgs2$weights) != 1 | !is.numeric(iListArgs2$weights) | any(iListArgs2$weights < 0)) {
233
        stop("'weights' must be a single (list of) positive or equal to zero value(s)", call. = FALSE)
234
235
236
      }
    }
    
237
238
239
    ## check 'epsilon'
    if (!is.null(iListArgs2$epsilon)) {
      if (!is.vector(iListArgs2$epsilon) | length(iListArgs2$epsilon) != 1 | !is.numeric(iListArgs2$epsilon) | any(iListArgs2$epsilon <= 0)) {
240
        stop("'epsilon' must be a single (list of) positive value(s)", call. = FALSE)
241
      }
242
243
    } else if (iListArgs2$transfo %in% c("log", "inv") & any(iListArgs2$Obs %in% 0) & warnings) {
      warning("zeroes detected in Obs: the corresponding time-steps will be excluded by the 'ErrorCrit*' functions as the epsilon argument was set to NULL", call. = FALSE)
244
245
    }
    
246
247
    ## check 'transfo' + 'FUN_CRIT'
    if (iListArgs2$transfo == "log" & warnings) {
248
      warn_log_kge <- "we do not advise using the %s with a log transformation on Obs (see the details section in the 'CreateInputsCrit' help)"
249
250
251
252
253
254
255
      if (identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE)) {
        warning(sprintf(warn_log_kge, "KGE"), call. = FALSE)
      }
      if (identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE2)) {
        warning(sprintf(warn_log_kge, "KGE'"), call. = FALSE)
      }
    }
256

257
258
    ## Create InputsCrit
    iInputsCrit <- list(FUN_CRIT   = iListArgs2$FUN_CRIT,
259
                        Obs        = iListArgs2$Obs,
260
                        VarObs     = iListArgs2$VarObs,
261
                        BoolCrit   = iListArgs2$BoolCrit,
262
                        idLayer    = iListArgs2$idLayer,
263
264
265
266
267
                        transfo    = iListArgs2$transfo,
                        epsilon    = iListArgs2$epsilon,
                        weights    = iListArgs2$weights)
    class(iInputsCrit) <- c("Single", "InputsCrit", ObjectClass)
    return(iInputsCrit)
268
    
269
270
271
  })
  names(InputsCrit) <- paste0("IC", seq_along(InputsCrit))
  
272
  
273
  listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
274
275
276
  inCnVarObs <- c("SCA", "SWE")
  if (!"ZLayers" %in% names(InputsModel)) {
    if(any(listVarObs %in% inCnVarObs)) {
277
      stop(sprintf("'VarObs' can not be equal to %i if CemaNeige is not used",
278
279
280
281
282
283
                   paste(sapply(inCnVarObs, shQuote), collapse = " or ")))
    }
  } else {
    listGroupLayer0 <- sapply(InputsCrit, FUN = "[[", "idLayer")
    listGroupLayer <- rep(listVarObs, times = listGroupLayer0)
    tabGroupLayer  <- as.data.frame(table(listGroupLayer))
284
    colnames(tabGroupLayer) <- c("VarObs", "freq")
285
286
287
    nLayers <- length(InputsModel$ZLayers)
    for (iInCnVarObs in inCnVarObs) {
      if (any(listVarObs %in% iInCnVarObs)) {
288
        if (tabGroupLayer[tabGroupLayer$VarObs %in% iInCnVarObs, "freq"] != nLayers) {
289
          stop(sprintf("'Obs' must contain %i vector(s) about %s", nLayers, iInCnVarObs))
290
291
292
293
294
295
296
297
298
299
300
301
302
        }
      }
    }
  }
  
  ## define idLayer as an index of the layer to use
  for (iInCnVarObs in unique(listVarObs)) {
    if (iInCnVarObs == "Q") {
      for (i in which(listVarObs == iInCnVarObs)) {
        InputsCrit[[i]]$idLayer <- NA
      }
    } else {
      aa <- listGroupLayer0[listVarObs == iInCnVarObs]
303
304
      aa <- unname(aa)
      bb <- cumsum(c(0, aa[-length(aa)]))
305
306
307
308
309
310
311
312
313
314
315
      cc <- lapply(seq_along(aa), function(x) seq_len(aa[x]) + bb[x])
      k <- 1
      for (i in which(listVarObs == iInCnVarObs)) {
        InputsCrit[[i]]$idLayer <- cc[[k]]
        k <- k + 1
      }
    }
  }

  
  ## if only one criterion --> not a list of InputsCrit but directly an InputsCrit
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
  if (length(InputsCrit) < 2) {
    InputsCrit <- InputsCrit[[1L]]
    InputsCrit["weights"] <- list(weights = NULL)
  } else {
    if (any(sapply(listArgs$weights, is.null))) {
      for (iListArgs in InputsCrit) {
        iListArgs$weights <- NULL
      }
      class(InputsCrit) <- c("Multi", "InputsCrit", ObjectClass)
    } else {
      class(InputsCrit) <- c("Compo", "InputsCrit", ObjectClass)
    }
    combInputsCrit <- combn(x = length(InputsCrit), m = 2)
    apply(combInputsCrit, MARGIN = 2, function(i) {
      equalInputsCrit <- identical(InputsCrit[[i[1]]], InputsCrit[[i[2]]])
      if(equalInputsCrit) {
332
        warning(sprintf("elements %i and %i of the criteria list are identical. This might not be necessary", i[1], i[2]), call. = FALSE)
333
334
      }
    })
335
  }
336
337
338
339
  
  return(InputsCrit)
  
}