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_NSE.R

-ErrorCrit_NSE <- function(InputsCrit,OutputsModel, warnings = TRUE, verbose = TRUE){
+ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
 
+  ## Arguments check
+  if (!inherits(OutputsModel, "OutputsModel")) {
+    stop("'OutputsModel' must be of class 'OutputsModel'")
+  }
 
-##Arguments_check________________________________
-  if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }  
-  if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); }  
+  EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "NSE", OutputsModel = OutputsModel, warnings = warnings)
 
+  CritValue <- NA
 
-##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
+  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 {
+      Crit <- (1 - Emod / Eref)
+    }
+    if (is.numeric(Crit) & is.finite(Crit)) {
+      CritValue <- Crit
+    }
 
-##Data_preparation_______________________________
-  VarObs <- InputsCrit$Qobs  ; VarObs[!InputsCrit$BoolCrit] <- NA; 
-  VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA;  
-  ##Data_transformation
-  if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
-    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)
+    ## Verbose
+    if (verbose) {
+      message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
+    }
   }
-  ##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);
-  if(Emod==0 & Eref==0){ Crit <- 0; } else { Crit <- (1-Emod/Eref); }
-  if(is.numeric(Crit) & is.finite(Crit)){ CritValue <- Crit; }
 
 
-##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)
+  ## 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/PEdaily_Oudin.R

-PEdaily_Oudin <- function(JD,Temp,LatRad){
-
-    PE_Oudin_D <- rep(NA,length(Temp));
-    for(k in 1:length(Temp)){
-
-      FI <- LatRad  ### latitude in rad
-      ### FI <- LatDeg/(180/pi)  ### conversion from deg to rad
-      COSFI <- cos(FI)
-      AFI <- abs(LatRad/42.)
-     
-      TETA <- 0.4093*sin(JD[k]/58.1-1.405)
-      COSTETA <- cos(TETA)
-      COSGZ <- max(0.001,cos(FI-TETA))
-      GZ <- acos(COSGZ)
-      COSGZ2 <- COSGZ*COSGZ
-      if(COSGZ2 >= 1){ SINGZ <- 0. } else { SINGZ <- sqrt(1.-COSGZ2) }
-      COSOM <- 1.-COSGZ/COSFI/COSTETA
-      if(COSOM < -1.){ COSOM <- -1. }
-      if(COSOM >  1.){ COSOM <-  1. }
-      COSOM2 <- COSOM*COSOM
-      if(COSOM2 >= 1.){ SINOM <- 0. } else { SINOM <- sqrt(1.-COSOM2) }
-      OM <- acos(COSOM)
-      COSPZ <- COSGZ+COSFI*COSTETA*(SINOM/OM-1.)
-      if(COSPZ < 0.001){ COSPZ <- 0.001 }
-      ETA <- 1.+cos(JD[k]/58.1)/30.
-      GE <- 446.*OM*COSPZ*ETA
-      
-      if(Temp[k] >= -5.0) { PE_Oudin_D[k] <- GE*(Temp[k]+5.)/100./28.5 } else { PE_Oudin_D[k] <- 0 }
-
-    }
-
-    return(PE_Oudin_D);
-
-}
-