R/ErrorCrit.R

-ErrorCrit <- function(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbose = TRUE) {
-  
+ErrorCrit <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
+
+
   ## ---------- Arguments check
-  
+
   if (!inherits(InputsCrit, "InputsCrit")) {
-    stop("InputsCrit must be of class 'InputsCrit' \n")
-    return(NULL)
-  }  
+    stop("InputsCrit must be of class 'InputsCrit'")
+  }
   if (!inherits(OutputsModel, "OutputsModel")) {
-    stop("OutputsModel must be of class 'OutputsModel' \n")
-    return(NULL)
-  }  
-  if (!missing(FUN_CRIT)) {
-    warning("argument 'FUN_CRIT' is deprecated. The error criterion function is now automatically get from the 'InputsCrit' object", call. = FALSE)
+    stop("OutputsModel must be of class 'OutputsModel'")
   }
-  
-  
-  
+
+
   ## ---------- Criterion computation
-  
+
   ## ----- Single criterion
   if (inherits(InputsCrit, "Single")) {
-    OutputsCrit <- InputsCrit$FUN_CRIT(InputsCrit,
-                                       OutputsModel,
-                                       warnings = warnings,
-                                       verbose = verbose)
+    FUN_CRIT <- match.fun(InputsCrit$FUN_CRIT)
+    OutputsCrit <- FUN_CRIT(InputsCrit = InputsCrit,
+                            OutputsModel = OutputsModel,
+                            warnings = warnings,
+                            verbose = verbose)
   }
-  
-  
+
+
   ## ----- Multiple criteria or Composite criterion
-  
+
   if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
     listOutputsCrit <- lapply(InputsCrit, FUN = function(iInputsCrit) {
-      iInputsCrit$FUN_CRIT(iInputsCrit, OutputsModel, warnings = warnings, verbose = verbose)
+      FUN_CRIT <- match.fun(iInputsCrit$FUN_CRIT)
+      FUN_CRIT(InputsCrit = iInputsCrit,
+               OutputsModel = OutputsModel,
+               warnings = warnings,
+               verbose = verbose)
     })
-    
+
     listValCrit  <- sapply(listOutputsCrit, function(x) x[["CritValue"]])
     listNameCrit <- sapply(listOutputsCrit, function(x) x[["CritName"]])
-    listweights  <- unlist(lapply(InputsCrit, function(x) x[["weights"]]))
-    
-    
+    listweights  <- unlist(lapply(InputsCrit, function(x) x[["Weights"]]))
+    listweights  <- listweights / sum(listweights)
+
     if (inherits(InputsCrit, "Compo")) {
-      CritValue <- sum(listValCrit * (listweights / sum(listweights)))
+      CritValue <- sum(listValCrit * listweights)
       OutputsCritCompo <- list(MultiCritValues  = listValCrit,
                                MultiCritNames   = listNameCrit,
                                MultiCritWeights = listweights)
@@ -57,19 +57,19 @@ ErrorCrit <- function(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbo
         message("Crit. Composite = ", sprintf("%.4f", CritValue))
         msgForm <- paste(sprintf("%.2f", listweights), listNameCrit, sep = " * ", collapse = ", ")
         msgForm <- unlist(strsplit(msgForm, split = ","))
-        msgFormSep <- rep(c(",", ",", ",\n\t\t     "), times = ceiling(length(msgForm)/3))[1: length(msgForm)]
+        msgFormSep <- rep(c(",", ",", ",\n\t\t    "), times = ceiling(length(msgForm)/3))[1:length(msgForm)]
         msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
-        msgForm <- gsub("\\,\\\n\\\t\\\t     $|\\,$", "", msgForm)
-        message("\tFormula: mean(", msgForm, ")\n")
+        msgForm <- gsub("\\,\\\n\\\t\\\t    $|\\,$", "", msgForm)
+        message("\tFormula: sum(", msgForm, ")\n")
       }
     } else {
       OutputsCrit <- listOutputsCrit
       class(OutputsCrit) <- c("Multi", "ErrorCrit")
     }
-    
+
   }
 
   return(OutputsCrit)
-  
+
 }
 

R/ErrorCrit_KGE.R

 ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
-  
-  
-    ##Arguments_check________________________________
-    if (!inherits(InputsCrit, "InputsCrit")) {
-      stop("InputsCrit must be of class 'InputsCrit' \n")
-      return(NULL)
-    }
-    if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
-      stop("InputsCrit must be of class 'Single' \n")
-      return(NULL)
-    }
-    if (!inherits(OutputsModel, "OutputsModel")) {
-      stop("OutputsModel must be of class 'OutputsModel' \n")
-      return(NULL)
-    }
-    
-    
-    ##Initialisation_________________________________
-    CritName <- NA
-    if (InputsCrit$transfo == "") {
-      CritName <- "KGE[Q]"
-    }
-    if (InputsCrit$transfo == "sqrt") {
-      CritName <- "KGE[sqrt(Q)]"
-    }
-    if (InputsCrit$transfo == "log") {
-      CritName <- "KGE[log(Q)]"
-    }
-    if (InputsCrit$transfo == "inv") {
-      CritName <- "KGE[1/Q]"
-    }
-    if (InputsCrit$transfo == "sort") {
-      CritName <- "KGE[sort(Q)]"
-    }
-    CritValue       <- NA
-    CritBestValue   <- +1
-    Multiplier      <- -1
-    ### must be equal to -1 or +1 only
-    
-    
-    ##Data_preparation_______________________________
-    VarObs <- InputsCrit$obs
-    VarObs[!InputsCrit$BoolCrit] <- NA
-    VarSim <- OutputsModel$Qsim
-    VarSim[!InputsCrit$BoolCrit] <- NA
-    
-    ##Data_transformation
-    if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) {
-      if (any(VarObs %in% 0)) {
-        warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL")
-      }
-      if (any(VarSim %in% 0)) {
-        warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL")
-      }  
-    }
-    if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) {
-      VarObs <- VarObs + InputsCrit$epsilon
-      VarSim <- VarSim + InputsCrit$epsilon
-    }
-    if (InputsCrit$transfo == "sqrt") {
-      VarObs <- sqrt(VarObs)
-      VarSim <- sqrt(VarSim)
-    }
-    if (InputsCrit$transfo == "log") {
-      VarObs <- log(VarObs)
-      VarSim <- log(VarSim)
-      VarSim[VarSim      < -1e100] <- NA
-    }
-    if (InputsCrit$transfo == "inv") {
-      VarObs <- 1 / VarObs
-      VarSim <- 1 / VarSim
-      VarSim[abs(VarSim) > 1e+100] <- NA
-    }
-    if (InputsCrit$transfo == "sort") {
-      VarSim[is.na(VarObs)] <- NA
-      VarSim <- sort(VarSim, na.last = TRUE)
-      VarObs <- sort(VarObs, na.last = TRUE)
-      InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
-    }
-    
-    ##TS_ignore
-    TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit
-    Ind_TS_ignore <- which(TS_ignore)
-    if (length(Ind_TS_ignore) == 0) {
-      Ind_TS_ignore <- NULL
-    }
-    if (sum(!TS_ignore) == 0) {
-      OutputsCrit <- list(NA)
-      names(OutputsCrit) <- c("CritValue")
-      return(OutputsCrit)
-    }
-    if (sum(!TS_ignore) == 1) {
-      OutputsCrit <- list(NA)
-      names(OutputsCrit) <- c("CritValue")
-      return(OutputsCrit)
-    } ### to avoid a problem in standard deviation computation
-    if (inherits(OutputsModel, "hourly")) {
-      WarningTS <- 365
-    }
-    if (inherits(OutputsModel, "daily")) {
-      WarningTS <- 365
-    }
-    if (inherits(OutputsModel, "monthly")) {
-      WarningTS <-  12
-    }
-    if (inherits(OutputsModel, "yearly")) {
-      WarningTS <-   3
-    }
-    if (sum(!TS_ignore) < WarningTS & warnings) {
-      warning("\t criterion computed on less than ", WarningTS, " time-steps ")
-    }
-    
-    ##Other_variables_preparation
-    meanVarObs <- mean(VarObs[!TS_ignore])
-    meanVarSim <- mean(VarSim[!TS_ignore])
-    
-    iCrit           <- 0
-    SubCritPrint    <- NULL
-    SubCritNames    <- NULL
-    SubCritValues   <- NULL
-    
-    
-    ##SubErrorCrit_____KGE_rPearson__________________
-    iCrit <- iCrit + 1
-    SubCritPrint[iCrit]  <- paste(CritName, " cor(sim, obs, \"pearson\") =", sep = "")
-    SubCritValues[iCrit] <- NA
-    SubCritNames[iCrit]  <- "r"
-    
-    Numer <- sum((VarObs[!TS_ignore] - meanVarObs) * (VarSim[!TS_ignore] - meanVarSim))
-    Deno1 <- sqrt(sum((VarObs[!TS_ignore] - meanVarObs) ^ 2))
-    Deno2 <- sqrt(sum((VarSim[!TS_ignore] - meanVarSim) ^ 2))
-    
+
+  ## Arguments check
+  if (!inherits(OutputsModel, "OutputsModel")) {
+    stop("'OutputsModel' must be of class 'OutputsModel'")
+  }
+
+  EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "KGE", OutputsModel = OutputsModel, warnings = warnings)
+
+  CritValue     <- NA
+  SubCritValues <- rep(NA, 3)
+  SubCritNames  <- c("r", "alpha", "beta")
+  SubCritPrint  <- rep(NA, 3)
+
+  if (EC$CritCompute) {
+    ## Other variables preparation
+    meanVarObs <- mean(EC$VarObs[!EC$TS_ignore])
+    meanVarSim <- mean(EC$VarSim[!EC$TS_ignore])
+
+    ## SubErrorCrit KGE rPearson
+    SubCritPrint[1L] <- paste0(EC$CritName, " cor(sim, obs, \"pearson\") =")
+
+    Numer <- sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) * (EC$VarSim[!EC$TS_ignore] - meanVarSim))
+    Deno1 <- sqrt(sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) ^ 2))
+    Deno2 <- sqrt(sum((EC$VarSim[!EC$TS_ignore] - meanVarSim) ^ 2))
+
     if (Numer == 0) {
       if (Deno1 == 0 & Deno2 == 0) {
         Crit <- 1
@@ -138,72 +31,65 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
         Crit <- 0
       }
     } else {
-      Crit  <- Numer / (Deno1 * Deno2)
+      Crit <- Numer / (Deno1 * Deno2)
     }
     if (is.numeric(Crit) & is.finite(Crit)) {
-      SubCritValues[iCrit] <- Crit
-    }
-    
-    
-    ##SubErrorCrit_____KGE_alpha_____________________
-    iCrit <- iCrit + 1
-    SubCritPrint[iCrit]  <- paste(CritName, " sd(sim)/sd(obs)          =", sep = "")
-    SubCritValues[iCrit] <- NA
-    SubCritNames[iCrit]  <- "alpha"
-    
-    Numer <- sd(VarSim[!TS_ignore])
-    Denom <- sd(VarObs[!TS_ignore])
-    
+      SubCritValues[1L] <- Crit
+    }
+
+    ## SubErrorCrit KGE alpha
+    SubCritPrint[2L] <- paste0(EC$CritName, " sd(sim)/sd(obs)          =")
+
+    Numer <- sd(EC$VarSim[!EC$TS_ignore])
+    Denom <- sd(EC$VarObs[!EC$TS_ignore])
+
     if (Numer == 0 & Denom == 0) {
       Crit <- 1
     } else {
       Crit <- Numer / Denom
     }
     if (is.numeric(Crit) & is.finite(Crit)) {
-      SubCritValues[iCrit] <- Crit
+      SubCritValues[2L] <- Crit
     }
-    
-    
-    ##SubErrorCrit_____KGE_beta______________________
-    iCrit <- iCrit + 1
-    SubCritPrint[iCrit]  <- paste(CritName, " mean(sim)/mean(obs)      =", sep = "")
-    SubCritValues[iCrit] <- NA
-    SubCritNames[iCrit]  <- "beta"
-    
+
+    ## SubErrorCrit KGE beta
+    SubCritPrint[3L] <- paste0(EC$CritName, " mean(sim)/mean(obs)      =")
+
     if (meanVarSim == 0 & meanVarObs == 0) {
       Crit <- 1
     } else {
       Crit <- meanVarSim / meanVarObs
     }
     if (is.numeric(Crit) & is.finite(Crit)) {
-      SubCritValues[iCrit] <- Crit
+      SubCritValues[3L] <- Crit
     }
-    
-    
-    ##ErrorCrit______________________________________
+
+    ## ErrorCrit
     if (sum(is.na(SubCritValues)) == 0) {
-      CritValue <- (1 - sqrt((SubCritValues[1] - 1)^2 + (SubCritValues[2] - 1)^2 + (SubCritValues[3] - 1)^2))
+      CritValue <- (1 - sqrt((SubCritValues[1L] - 1)^2 + (SubCritValues[2L] - 1)^2 + (SubCritValues[3L] - 1)^2))
     }
-    
-    
-    ##Verbose______________________________________
+
+    ## Verbose
     if (verbose) {
-      message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue))
+      message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
       message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " "))
     }
-    
-    
-    ##Output_________________________________________
-    OutputsCrit <- list(CritValue       = CritValue,
-                        CritName        = CritName,
-                        SubCritValues   = SubCritValues,
-                        SubCritNames    = SubCritNames,
-                        CritBestValue   = CritBestValue,
-                        Multiplier      = Multiplier,
-                        Ind_notcomputed = Ind_TS_ignore
-    )
-    
-    class(OutputsCrit) <- c("KGE", "ErrorCrit")
-    return(OutputsCrit)
-    
   }
+
+
+  ## Output
+  OutputsCrit <- list(CritValue       = CritValue,
+                      CritName        = EC$CritName,
+                      SubCritValues   = SubCritValues,
+                      SubCritNames    = SubCritNames,
+                      CritBestValue   = EC$CritBestValue,
+                      Multiplier      = EC$Multiplier,
+                      Ind_notcomputed = EC$Ind_TS_ignore)
+
+  class(OutputsCrit) <- c("KGE", "ErrorCrit")
+  return(OutputsCrit)
+
+}
+
+class(ErrorCrit_KGE) <- c("FUN_CRIT", class(ErrorCrit_KGE))
+

R/ErrorCrit_NSE.R

 ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
-  
-  
-    ##Arguments_check________________________________
-    if (!inherits(InputsCrit, "InputsCrit")) {
-      stop("InputsCrit must be of class 'InputsCrit' \n")
-      return(NULL)
-    }
-	if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
-      stop("InputsCrit must be of class 'Single' \n")
-      return(NULL)
-    }
-    if (!inherits(OutputsModel, "OutputsModel")) {
-      stop("OutputsModel must be of class 'OutputsModel' \n")
-      return(NULL)
-    }
-    
-    
-    ##Initialisation_________________________________
-    CritName <- NA
-    if (InputsCrit$transfo == "") {
-      CritName <- "NSE[Q]"
-    }
-    if (InputsCrit$transfo == "sqrt") {
-      CritName <- "NSE[sqrt(Q)]"
-    }
-    if (InputsCrit$transfo == "log") {
-      CritName <- "NSE[log(Q)]"
-    }
-    if (InputsCrit$transfo == "inv") {
-      CritName <- "NSE[1/Q]"
-    }
-    if (InputsCrit$transfo == "sort") {
-      CritName <- "NSE[sort(Q)]"
-    }
-    CritValue       <- NA
-    CritBestValue   <- +1
-    Multiplier      <- -1
-    ### must be equal to -1 or +1 only
-    
-    
-    ##Data_preparation_______________________________
-    VarObs <- InputsCrit$obs
-    VarObs[!InputsCrit$BoolCrit] <- NA
-    VarSim <- OutputsModel$Qsim
-    VarSim[!InputsCrit$BoolCrit] <- NA
-    
-    ##Data_transformation
-    if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) {
-      if (any(VarObs %in% 0)) {
-        warning("zeroes detected in Qobs: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL")
-      }
-      if (any(VarSim %in% 0)) {
-        warning("zeroes detected in Qsim: the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL")
-      }  
-    }
-    if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) {
-      VarObs <- VarObs + InputsCrit$epsilon
-      VarSim <- VarSim + InputsCrit$epsilon
-    }
-    if (InputsCrit$transfo == "sqrt") {
-      VarObs <- sqrt(VarObs)
-      VarSim <- sqrt(VarSim)
-    }
-    if (InputsCrit$transfo == "log") {
-      VarObs <- log(VarObs)
-      VarSim <- log(VarSim)
-      VarSim[VarSim      < -1e100] <- NA
-    }
-    if (InputsCrit$transfo == "inv") {
-      VarObs <- 1 / VarObs
-      VarSim <- 1 / VarSim
-      VarSim[abs(VarSim) > 1e+100] <- NA
-    }
-    if (InputsCrit$transfo == "sort") {
-      VarSim[is.na(VarObs)] <- NA
-      VarSim <- sort(VarSim, na.last = TRUE)
-      VarObs <- sort(VarObs, na.last = TRUE)
-      InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
-    }
-    
-    ##TS_ignore
-    TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit
-    Ind_TS_ignore <-  which(TS_ignore)
-    if (length(Ind_TS_ignore) == 0) {
-      Ind_TS_ignore <- NULL
-    }
-    if (sum(!TS_ignore) == 0) {
-      OutputsCrit <- list(NA)
-      names(OutputsCrit) <- c("CritValue")
-      return(OutputsCrit)
-    }
-    if (inherits(OutputsModel, "hourly")) {
-      WarningTS <- 365
-    }
-    if (inherits(OutputsModel, "daily")) {
-      WarningTS <- 365
-      
-    }
-    if (inherits(OutputsModel, "monthly")) {
-      WarningTS <-  12
-    }
-    if (inherits(OutputsModel, "yearly")) {
-      WarningTS <-   3
-    }
-    if (sum(!TS_ignore) < WarningTS & warnings) {
-      warning("\t criterion computed on less than ", WarningTS, " time-steps")
-    }
-    
-    ##Other_variables_preparation
-    meanVarObs <- mean(VarObs[!TS_ignore])
-    meanVarSim <- mean(VarSim[!TS_ignore])
-    
-    
-    ##ErrorCrit______________________________________
-    Emod <- sum((VarSim[!TS_ignore] - VarObs[!TS_ignore])^2)
-    Eref <- sum((VarObs[!TS_ignore] - mean(VarObs[!TS_ignore]))^2)
-    
+
+  ## Arguments check
+  if (!inherits(OutputsModel, "OutputsModel")) {
+    stop("'OutputsModel' must be of class 'OutputsModel'")
+  }
+
+  EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "NSE", OutputsModel = OutputsModel, warnings = warnings)
+
+  CritValue <- NA
+
+  if (EC$CritCompute) {
+    ## ErrorCrit
+    Emod <- sum((EC$VarSim[!EC$TS_ignore] - EC$VarObs[!EC$TS_ignore])^2)
+    Eref <- sum((EC$VarObs[!EC$TS_ignore] - mean(EC$VarObs[!EC$TS_ignore]))^2)
+
     if (Emod == 0 & Eref == 0) {
       Crit <- 0
     } else {
@@ -124,21 +22,24 @@ ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
     if (is.numeric(Crit) & is.finite(Crit)) {
       CritValue <- Crit
     }
-    
-    
-    ##Verbose______________________________________
+
+    ## Verbose
     if (verbose) {
-      message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue), "\n")
-    }
-    
-    
-    ##Output_________________________________________
-    OutputsCrit <- list(CritValue       = CritValue,
-                        CritName        = CritName,
-                        CritBestValue   = CritBestValue,
-                        Multiplier      = Multiplier,
-                        Ind_notcomputed = Ind_TS_ignore)
-    
-    class(OutputsCrit) <- c("NSE", "ErrorCrit")
-    return(OutputsCrit)
+      message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
+    }
+  }
+
+
+  ## Output
+  OutputsCrit <- list(CritValue       = CritValue,
+                      CritName        = EC$CritName,
+                      CritBestValue   = EC$CritBestValue,
+                      Multiplier      = EC$Multiplier,
+                      Ind_notcomputed = EC$Ind_TS_ignore)
+
+  class(OutputsCrit) <- c("NSE", "ErrorCrit")
+  return(OutputsCrit)
+
 }
+
+class(ErrorCrit_NSE) <- c("FUN_CRIT", class(ErrorCrit_NSE))

R/RunModel.R

-RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD) {
-  return(FUN_MOD(InputsModel, RunOptions, Param))
+RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, ...) {
+
+  FUN_MOD <- match.fun(FUN_MOD)
+
+  if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
+    # Lag model take one parameter at the beginning of the vector
+    iFirstParamRunOffModel <- 2
+    RunOptions$FeatFUN_MOD$NbParam <- RunOptions$FeatFUN_MOD$NbParam - 1
+  } else {
+    # All parameters
+    iFirstParamRunOffModel <- 1
+  }
+
+  OutputsModel <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions,
+                          Param = Param[iFirstParamRunOffModel:length(Param)], ...)
+
+  if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
+    OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel)
+  }
+  return(OutputsModel)
 }