R/ErrorCrit.R

-ErrorCrit <- function(InputsCrit,OutputsModel,FUN_CRIT, warnings = TRUE, verbose = TRUE){
-    return( FUN_CRIT(InputsCrit,OutputsModel, warnings = warnings, verbose = verbose) )
+ErrorCrit <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
+
+
+  ## ---------- Arguments check
+
+  if (!inherits(InputsCrit, "InputsCrit")) {
+    stop("InputsCrit must be of class 'InputsCrit'")
+  }
+  if (!inherits(OutputsModel, "OutputsModel")) {
+    stop("OutputsModel must be of class 'OutputsModel'")
+  }
+
+
+  ## ---------- Criterion computation
+
+  ## ----- Single criterion
+  if (inherits(InputsCrit, "Single")) {
+    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) {
+      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  <- listweights / sum(listweights)
+
+    if (inherits(InputsCrit, "Compo")) {
+      CritValue <- sum(listValCrit * listweights)
+      OutputsCritCompo <- list(MultiCritValues  = listValCrit,
+                               MultiCritNames   = listNameCrit,
+                               MultiCritWeights = listweights)
+      OutputsCrit <- list(CritValue       = CritValue,
+                          CritName        = "Composite",
+                          CritBestValue   = +1,
+                          Multiplier      = -1,
+                          Ind_notcomputed = NULL,
+                          CritCompo       = OutputsCritCompo,
+                          MultiCrit       = listOutputsCrit)
+      class(OutputsCrit) <- c("Compo", "ErrorCrit")
+      if (verbose) {
+        message("------------------------------------\n")
+        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)]
+        msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
+        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(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$Qobs
-    VarObs[!InputsCrit$BoolCrit] <- NA
-    VarSim <- OutputsModel$Qsim
-    VarSim[!InputsCrit$BoolCrit] <- NA
-    
-    ##Data_transformation
-    if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & verbose) {
-      if (any(VarObs %in% 0)) {
-        warning("zeroes detected in Qobs: the corresponding time-steps will be exclude from the criteria computation if the epsilon agrument of 'CreateInputsCrit' = NULL")
-      }
-      if (any(VarSim %in% 0)) {
-        warning("zeroes detected in Qsim: the corresponding time-steps will be exclude from the criteria computation if the epsilon agrument 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
@@ -134,70 +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
-    )
-    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(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$Qobs
-    VarObs[!InputsCrit$BoolCrit] <- NA
-    VarSim <- OutputsModel$Qsim
-    VarSim[!InputsCrit$BoolCrit] <- NA
-    
-    ##Data_transformation
-    if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & verbose) {
-      if (any(VarObs %in% 0)) {
-        warning("zeroes detected in Qobs: the corresponding time-steps will be exclude from the criteria computation if the epsilon agrument of 'CreateInputsCrit' = NULL")
-      }
-      if (any(VarSim %in% 0)) {
-        warning("zeroes detected in Qsim: the corresponding time-steps will be exclude from the criteria computation if the epsilon agrument 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 {
@@ -120,19 +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))
-    }
-    
-    
-    ##Output_________________________________________
-    OutputsCrit <- list(CritValue       = CritValue,
-                        CritName        = CritName,
-                        CritBestValue   = CritBestValue,
-                        Multiplier      = Multiplier,
-                        Ind_notcomputed = Ind_TS_ignore)
-    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/Imax.R

+
+Imax <- function(InputsModel,
+                 IndPeriod_Run,
+                 TestedValues = seq(from = 0.1, to = 3, by = 0.1)) {
+
+
+  ## ---------- check arguments
+
+  ##  InputsModel
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+  }
+  if (!inherits(InputsModel, "hourly")) {
+    stop("'InputsModel' must be of class 'hourly'")
+  }
+
+  ## IndPeriod_Run
+  if (!is.vector(IndPeriod_Run)) {
+    stop("'IndPeriod_Run' must be a vector of numeric values")
+  }
+  if (!inherits(IndPeriod_Run, "integer")) {
+    stop("'IndPeriod_Run' must be of type integer")
+  }
+  if (!identical(as.integer(IndPeriod_Run), IndPeriod_Run[1]:IndPeriod_Run[length(IndPeriod_Run)])) {
+    stop("'IndPeriod_Run' must be a continuous sequence of integers")
+  }
+
+  ## TestedValues
+  if (!(is.numeric(TestedValues))) {
+    stop("'TestedValues' must be 'numeric'")
+  }
+
+
+  ## ---------- hourly inputs aggregation
+
+  ## aggregate data at the daily time step
+  daily_data <- SeriesAggreg(InputsModel[IndPeriod_Run], Format = "%Y%m%d")
+
+
+  ## ---------- calculate interception
+
+  ## calculate total interception of daily GR models on the period
+  cum_daily <- sum(pmin(daily_data$Precip, daily_data$PotEvap))
+  if (anyNA(cum_daily)) {
+    stop("'IndPeriod_Run' must be set to select 24 hours by day")
+  }
+
+  ## calculate total interception of the GR5H interception store on the period
+  ## and compute difference with daily values
+  differences <- array(NA, c(length(TestedValues)))
+  for (Imax in TestedValues) {
+    C0 <- 0
+    cum_hourly <- 0
+    for (i in IndPeriod_Run) {
+      Ec  <- min(InputsModel$PotEvap[i], InputsModel$Precip[i] + C0)
+      Pth <- max(0, InputsModel$Precip[i] - (Imax-C0) - Ec)
+      C0  <- C0 + InputsModel$Precip[i] - Ec - Pth
+      cum_hourly <- cum_hourly + Ec
+    }
+    differences[which(Imax == TestedValues)] <- abs(cum_hourly - cum_daily)
+  }
+
+  ## return the Imax value that minimises the difference
+  return(TestedValues[which.min(differences)])
+
+}

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)
 }