Skip to content
GitLab
Select Git revision
  • 12-erreur-apiplatform-core-metadata-resource-resourcenamecollection
  • 2.x
  • 3.x default
  • 2-typage-des-champs-json
  • evolution-tests
  • refactoring
  • typescript
  • 1.0.x
  • 1.x
  • api-client
  • 0.1.x
  • 0.x
  • 3.0.0
  • v3.0.0-alpha-2
  • 3.0.0-alpha
  • 2.1.1
  • 2.1.0
  • 2.0.4
  • 2.0.3
  • 2.0.2
  • 2.0.1
  • 2.0.0
  • 1.0.2
  • 1.0.1
  • 1.0.0
  • 0.0.1
Find file BlameHistoryPermalink
ErrorCrit_KGE2.R 5.19 KiB
1 
ErrorCrit_KGE2 <- function(InputsCrit,OutputsModel, warnings = TRUE, verbose = TRUE){
2
3
4 
##Arguments_check________________________________
5 
  if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }
6 
  if(inherits(OutputsModel,"OutputsModel")==FALSE){ stop("OutputsModel must be of class 'OutputsModel' \n"); return(NULL); }
7
8
9 
##Initialisation_________________________________
10 
  CritName <- NA;
11 
  if(InputsCrit$transfo==""    ){ CritName <- "KGE'[Q]"      ; }
12 
  if(InputsCrit$transfo=="sqrt"){ CritName <- "KGE'[sqrt(Q)]"; }
13 
  if(InputsCrit$transfo=="log" ){ CritName <- "KGE'[log(Q)]" ; }
14 
  if(InputsCrit$transfo=="inv" ){ CritName <- "KGE'[1/Q]"    ; }
15 
  if(InputsCrit$transfo=="sort"){ CritName <- "KGE'[sort(Q)]"; }
16 
  CritValue       <- NA;
17 
  CritBestValue   <- +1;
18 
  Multiplier      <- -1; ### must be equal to -1 or +1 only
19
20
21 
##Data_preparation_______________________________
22 
  VarObs <- InputsCrit$Qobs  ; VarObs[!InputsCrit$BoolCrit] <- NA;
23 
  VarSim <- OutputsModel$Qsim; VarSim[!InputsCrit$BoolCrit] <- NA;
24 
  ##Data_transformation
25 
  if("Ind_zeroes" %in% names(InputsCrit) & "epsilon" %in% names(InputsCrit)){ if(length(InputsCrit$Ind_zeroes)>0){
26 
    VarObs <- VarObs + InputsCrit$epsilon;
27 
    VarSim <- VarSim + InputsCrit$epsilon;
28 
  } }
29 
  if(InputsCrit$transfo=="sqrt"){ VarObs <- sqrt(VarObs); VarSim <- sqrt(VarSim); }
30 
  if(InputsCrit$transfo=="log" ){ VarObs <- log(VarObs) ; VarSim <- log(VarSim) ; VarSim[VarSim      < -1E100] <- NA; }
31 
  if(InputsCrit$transfo=="inv" ){ VarObs <- 1/VarObs    ; VarSim <- 1/VarSim    ; VarSim[abs(VarSim) > 1E+100] <- NA; }
32 
  if(InputsCrit$transfo=="sort"){
33 
	VarSim[is.na(VarObs)] <- NA
34 
	VarSim <- sort(VarSim, na.last = TRUE)
35 
	VarObs <- sort(VarObs, na.last = TRUE)
36 
	InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
37 
  }
38 
  ##TS_ignore
39 
  TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit ;
40 
  Ind_TS_ignore <- which(TS_ignore); if(length(Ind_TS_ignore)==0){ Ind_TS_ignore <- NULL; }
41 
  if(sum(!TS_ignore)==0){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); }
42 
  if(sum(!TS_ignore)==1){ OutputsCrit <- list(NA); names(OutputsCrit) <- c("CritValue"); return(OutputsCrit); } ### to avoid a problem in standard deviation computation
43 
  if(inherits(OutputsModel,"hourly" )){ WarningTS <- 365; }
44 
  if(inherits(OutputsModel,"daily"  )){ WarningTS <- 365; }
45 
  if(inherits(OutputsModel,"monthly")){ WarningTS <-  12; }
46 
  if(inherits(OutputsModel,"yearly" )){ WarningTS <-   3; }
47 
  if(sum(!TS_ignore)<WarningTS & warnings){ warning("\t criterion computed on less than ", WarningTS, " time-steps") }
48 
  ##Other_variables_preparation
49 
  meanVarObs <- mean(VarObs[!TS_ignore]);
50 
  meanVarSim <- mean(VarSim[!TS_ignore]);
51 
  iCrit           <- 0
52 
  SubCritPrint    <- NULL
53 
  SubCritNames    <- NULL
54 
  SubCritValues   <- NULL
55
56
57 
##SubErrorCrit_____KGE_rPearson__________________
58 
  iCrit <- iCrit+1;
59 
  SubCritPrint[iCrit]  <- paste(CritName," cor(sim, obs, \"pearson\") =", sep = "")
60 
  SubCritNames[iCrit]  <- "r"
61 
  SubCritValues[iCrit] <- NA;
62 
  Numer <- sum( (VarObs[!TS_ignore]-meanVarObs)*(VarSim[!TS_ignore]-meanVarSim) );
63 
  Deno1 <- sqrt( sum((VarObs[!TS_ignore]-meanVarObs)^2) );
64 
  Deno2 <- sqrt( sum((VarSim[!TS_ignore]-meanVarSim)^2) );
65 
  if(Numer==0){ if(Deno1==0 & Deno2==0){ Crit <- 1; } else { Crit <- 0; }
66 
  } else { Crit  <- Numer/(Deno1*Deno2); }
67 
  if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; }
68
69
70 
##SubErrorCrit_____KGE_gamma______________________
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112
iCrit <- iCrit+1; SubCritPrint[iCrit] <- paste(CritName," cv(sim)/cv(obs) =", sep = "") SubCritNames[iCrit] <- "gamma" SubCritValues[iCrit] <- NA; if(meanVarSim==0){ if(sd(VarSim[!TS_ignore])==0){ CVsim <- 1; } else { CVsim <- 99999; } } else { CVsim <- sd(VarSim[!TS_ignore])/meanVarSim; } if(meanVarObs==0){ if(sd(VarObs[!TS_ignore])==0){ CVobs <- 1; } else { CVobs <- 99999; } } else { CVobs <- sd(VarObs[!TS_ignore])/meanVarObs; } if(CVsim==0 & CVobs==0){ Crit <- 1; } else { Crit <- CVsim/CVobs ; } if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; } ##SubErrorCrit_____KGE_beta______________________ iCrit <- iCrit+1; SubCritPrint[iCrit] <- paste(CritName," mean(sim)/mean(obs) =", sep = "") SubCritNames[iCrit] <- "beta" SubCritValues[iCrit] <- NA; if(meanVarSim==0 & meanVarObs==0){ Crit <- 1; } else { Crit <- meanVarSim/meanVarObs ; } if(is.numeric(Crit) & is.finite(Crit)){ SubCritValues[iCrit] <- Crit; } ##ErrorCrit______________________________________ if(sum(is.na(SubCritValues))==0){ CritValue <- ( 1 - sqrt( (SubCritValues[1]-1)^2 + (SubCritValues[2]-1)^2 + (SubCritValues[3]-1)^2 ) ); } ##Verbose______________________________________ if(verbose) { message("Crit. ", CritName, " = ", sprintf("%.4f", 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) }