Skip to content
GitLab

Source

Target

Hide whitespace changes
Inline Side-by-side
Showing
with 1213 additions and 541 deletions
+1213 -541
R/RunModel_GR4J.R
View file @ 4bccfee2
RunModel_GR4J <- function(InputsModel,RunOptions,Param){ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
NParam <- 4; .ArgumentsCheckGR(InputsModel, RunOptions, Param)
FortranOutputs <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
"AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
##Arguments_check Param <- as.double(Param)
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"daily" )==FALSE){ stop("InputsModel must be of class 'daily' \n"); return(NULL); }
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); }
if(!is.vector(Param) | !is.numeric(Param)){ stop("Param must be a numeric vector \n"); return(NULL); }
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); }
Param <- as.double(Param);
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Input_data_preparation Param_X1X3_threshold <- 1e-2
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; } Param_X4_threshold <- 0.5
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run); if (Param[1L] < Param_X1X3_threshold) {
LInputSeries <- as.integer(length(IndPeriod1)) warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs)); Param[1L] <- Param_X1X3_threshold
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); } }
##Input_data_preparation if (Param[3L] < Param_X1X3_threshold) {
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries; warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim; Param[3L] <- Param_X1X3_threshold
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim; }
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Use_of_IniResLevels ## Input data preparation
if(!is.null(RunOptions$IniResLevels)){ if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm) RunOptions$IndPeriod_WarmUp <- NULL
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm) }
} IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
##Call_fortan ## Use of IniResLevels
RESULTS <- .Fortran("frun_GR4J",PACKAGE="airGR", if (!is.null(RunOptions$IniResLevels)) {
##inputs RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
LInputs=LInputSeries, ### length of input and output series RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d] }
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
##Output_data_preparation ## Call GR model Fortan
##OutputsModel_only RESULTS <- .Fortran("frun_gr4j", PACKAGE = "airGR",
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){ ## inputs
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]); LInputs = LInputSeries, ### length of input and output series
names(OutputsModel) <- FortranOutputs[IndOutputs]; } InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
##DatesR_and_OutputsModel_only InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){ NParam = as.integer(length(Param)), ### number of model parameter
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), Param = Param, ### parameter set
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) ); NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); } StateStart = RunOptions$IniStates, ### state variables used when the model run starts
##OutputsModel_and_StateEnd_only NOutputs = as.integer(length(IndOutputs)), ### number of output series
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){ IndOutputs = IndOutputs, ### indices of output series
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), ## outputs
list(RESULTS$StateEnd) ); Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/d]
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); } StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
##DatesR_and_OutputsModel_and_StateEnd )
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]), if ("StateEnd" %in% RunOptions$Outputs_Sim) {
list(RESULTS$StateEnd) ); RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); } RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
##End UH1 = RESULTS$StateEnd[(1:20) + 7],
rm(RESULTS); UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
class(OutputsModel) <- c("OutputsModel","daily","GR"); GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
return(OutputsModel); verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
} }
R/RunModel_GR5H.R 0 → 100644
View file @ 4bccfee2
RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsIntStore <- inherits(RunOptions, "interception")
if (IsIntStore) {
Imax <- RunOptions$Imax
} else {
Imax <- -99
}
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Output data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
if (IsIntStore) {
RunOptions$IniStates[4] <- RunOptions$IniResLevels[4] * Imax ### interception store level (mm)
}
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr5h", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
Imax = Imax, ### maximal capacity of interception store
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/h]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5H, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
IntStore = RESULTS$StateEnd[4L],
UH1 = NULL,
UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
}
R/RunModel_GR5J.R
View file @ 4bccfee2
RunModel_GR5J <- function(InputsModel,RunOptions,Param){ RunModel_GR5J <- function(InputsModel, RunOptions, Param) {
NParam <- 5; .ArgumentsCheckGR(InputsModel, RunOptions, Param)
FortranOutputs <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
"AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
##Arguments_check Param <- as.double(Param)
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"daily" )==FALSE){ stop("InputsModel must be of class 'daily' \n"); return(NULL); }
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); }
if(!is.vector(Param) | !is.numeric(Param)){ stop("Param must be a numeric vector \n"); return(NULL); }
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); }
Param <- as.double(Param);
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs));
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); }
##Output_data_preparation Param_X1X3_threshold <- 1e-2
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries; Param_X4_threshold <- 0.5
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim; if (Param[1L] < Param_X1X3_threshold) {
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim; warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
##Use_of_IniResLevels }
if(!is.null(RunOptions$IniResLevels)){ if (Param[3L] < Param_X1X3_threshold) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm) warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm) Param[3L] <- Param_X1X3_threshold
} }
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Call_fortan ## Input data preparation
RESULTS <- .Fortran("frun_GR5J",PACKAGE="airGR", if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
##inputs RunOptions$IndPeriod_WarmUp <- NULL
LInputs=LInputSeries, ### length of input and output series }
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d] IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d] LInputSeries <- as.integer(length(IndPeriod1))
NParam=as.integer(length(Param)), ### number of model parameter if ("all" %in% RunOptions$Outputs_Sim) {
Param=Param, ### parameter set IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising } else {
StateStart=RunOptions$IniStates, ### state variables used when the model run starts IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
NOutputs=as.integer(length(IndOutputs)), ### number of output series }
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL, UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
##Output_data_preparation ## Output data preparation
##OutputsModel_only IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){ ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]); ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
##End ## Use of IniResLevels
rm(RESULTS); if (!is.null(RunOptions$IniResLevels)) {
class(OutputsModel) <- c("OutputsModel","daily","GR"); RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
return(OutputsModel); RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
}
} ## Call GR model Fortan
RESULTS <- .Fortran("frun_gr5j", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/d]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL,
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
}
R/RunModel_GR6J.R
View file @ 4bccfee2
RunModel_GR6J <- function(InputsModel,RunOptions,Param){ RunModel_GR6J <- function(InputsModel, RunOptions, Param) {
NParam <- 6; .ArgumentsCheckGR(InputsModel, RunOptions, Param)
FortranOutputs <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1",
"Rout", "Exch", "AExch1", "AExch2", "AExch", "QR", "QR1", "Exp", "QD", "Qsim");
##Arguments_check Param <- as.double(Param)
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(InputsModel,"daily" )==FALSE){ stop("InputsModel must be of class 'daily' \n"); return(NULL); }
if(inherits(InputsModel,"GR" )==FALSE){ stop("InputsModel must be of class 'GR' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }
if(inherits(RunOptions,"GR" )==FALSE){ stop("RunOptions must be of class 'GR' \n"); return(NULL); }
if(!is.vector(Param) | !is.numeric(Param)){ stop("Param must be a numeric vector \n"); return(NULL); }
if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); }
Param <- as.double(Param);
Param_X1X3X6_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[1L] <- Param_X1X3X6_threshold
}
if (Param[3L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[3L] <- Param_X1X3X6_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
if (Param[6L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[6L] <- Param_X1X3X6_threshold
}
##Input_data_preparation Param_X1X3X6_threshold <- 1e-2
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; } Param_X4_threshold <- 0.5
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run); if (Param[1L] < Param_X1X3X6_threshold) {
LInputSeries <- as.integer(length(IndPeriod1)) warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs)); Param[1L] <- Param_X1X3X6_threshold
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); } }
if (Param[3L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[3L] <- Param_X1X3X6_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
if (Param[6L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[6L] <- Param_X1X3X6_threshold
}
##Output_data_preparation ## Input data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries; if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim; RunOptions$IndPeriod_WarmUp <- NULL
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
##Use_of_IniResLevels LInputSeries <- as.integer(length(IndPeriod1))
if(!is.null(RunOptions$IniResLevels)){ if ("all" %in% RunOptions$Outputs_Sim) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm) IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm) } else {
RunOptions$IniStates[3] <- RunOptions$IniResLevels[3] ### exponential store level (mm) IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
} }
##Call_fortan ## Output data preparation
RESULTS <- .Fortran("frun_GR6J",PACKAGE="airGR", IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
##inputs ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
LInputs=LInputSeries, ### length of input and output series ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR6J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
##Output_data_preparation ## Use of IniResLevels
##OutputsModel_only if (!is.null(RunOptions$IniResLevels)) {
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){ RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]); RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
names(OutputsModel) <- FortranOutputs[IndOutputs]; } RunOptions$IniStates[3] <- RunOptions$IniResLevels[3] ### exponential store level (mm)
##DatesR_and_OutputsModel_only }
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
##End ## Call GR model Fortan
rm(RESULTS); RESULTS <- .Fortran("frun_gr6j", PACKAGE = "airGR",
class(OutputsModel) <- c("OutputsModel","daily","GR"); ## inputs
return(OutputsModel); LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/d]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR6J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
UH1 = RESULTS$StateEnd[(1:20) + 7],
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
} }
R/RunModel_Lag.R 0 → 100644
View file @ 4bccfee2
RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
NParam <- 1
## argument check
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
}
if (!inherits(InputsModel, "SD")) {
stop("'InputsModel' must be of class 'SD'")
}
if (!inherits(RunOptions, "RunOptions")) {
stop("'RunOptions' must be of class 'RunOptions'")
}
if (!is.vector(Param) | !is.numeric(Param)) {
stop("'Param' must be a numeric vector")
}
if (sum(!is.na(Param)) != NParam) {
stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
}
if (inherits(QcontribDown, "OutputsModel")) {
if (is.null(QcontribDown$Qsim)) {
stop("'QcontribDown' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
}
if (length(QcontribDown$Qsim) != length(RunOptions$IndPeriod_Run)) {
stop("Time series Qsim in 'QcontribDown' should have the same length as 'RunOptions$IndPeriod_Run'")
}
if (!identical(RunOptions$IndPeriod_WarmUp, 0L) && !identical(RunOptions$Outputs_Sim, RunOptions$Outputs_Cal)) {
# This test is not necessary during calibration but usefull in other cases because
# WarmUpQsim is then used for downstream sub-basins because of the delay in Qupstream
if (is.null(QcontribDown$RunOptions$WarmUpQsim) ||
length(QcontribDown$RunOptions$WarmUpQsim) != length(RunOptions$IndPeriod_WarmUp)) {
stop("Time series WarmUpQsim in 'QcontribDown' should have the same length as 'RunOptions$IndPeriod_WarmUp'")
}
}
} else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
if (length(QcontribDown) != length(RunOptions$IndPeriod_Run)) {
stop("'QcontribDown' should have the same length as 'RunOptions$IndPeriod_Run'")
}
} else {
stop("'QcontribDown' must be a numeric vector or a 'OutputsModel' object")
}
# data set up
NbUpBasins <- length(InputsModel$LengthHydro)
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
if (inherits(QcontribDown, "OutputsModel")) {
OutputsModel <- QcontribDown
if (is.null(OutputsModel$RunOptions$WarmUpQsim)) {
OutputsModel$RunOptions$WarmUpQsim <- rep(NA, length(RunOptions$IndPeriod_WarmUp))
}
QsimDown <- c(OutputsModel$RunOptions$WarmUpQsim, OutputsModel$Qsim)
} else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
OutputsModel <- list()
class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
QsimDown <- c(rep(NA, length(RunOptions$IndPeriod_WarmUp)),
QcontribDown)
}
## propagation time from upstream meshes to outlet
PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / RunOptions$FeatFUN_MOD$TimeStep
HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
## set up initial states
IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
if (length(IniSD) > 0) {
if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
stop(
sprintf(
"SD initial states has a length of %i and a length of %i is required",
length(IniSD),
sum(floor(PT)) + NbUpBasins
)
)
}
IniStates <- lapply(seq_len(NbUpBasins), function(x) {
iStart <- 1
if (x > 1) {
iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
}
as.vector(IniSD[iStart:(iStart + PT[x])])
})
} else {
IniStates <- lapply(
seq_len(NbUpBasins),
function(iUpBasins) {
iWarmUp <- seq.int(
from = max(1, IndPeriod1[1] - floor(PT[iUpBasins]) - 1),
to = max(1, IndPeriod1[1] - 1)
)
ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
if (length(ini) != floor(PT[iUpBasins] + 1)) {
# If warm-up period is not enough long complete beginning with first value
ini <- c(rep(ini[1], floor(PT[iUpBasins] + 1) - length(ini)), ini)
}
return(as.vector(ini))
}
)
}
# message("IniStates: ", paste(IniStates, collapse = ", "))
## Lag model computation
Qsim_m3 <- QsimDown *
InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
for (upstream_basin in seq_len(NbUpBasins)) {
Qupstream <- c(IniStates[[upstream_basin]],
InputsModel$Qupstream[IndPeriod1, upstream_basin])
# message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
Qsim_m3 <- Qsim_m3 +
Qupstream[2:(1 + LInputSeries)] * HUTRANS[1, upstream_basin] +
Qupstream[1:LInputSeries] * HUTRANS[2, upstream_basin]
}
## OutputsModel
if ("Qsim_m3" %in% RunOptions$Outputs_Sim) {
OutputsModel$Qsim_m3 <- Qsim_m3[IndPeriod2]
}
if ("Qsim" %in% RunOptions$Outputs_Sim) {
# Convert back Qsim to mm
OutputsModel$Qsim <- Qsim_m3[IndPeriod2] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
# message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
}
if ("QsimDown" %in% RunOptions$Outputs_Sim) {
# Convert back Qsim to mm
OutputsModel$QsimDown <- QsimDown[IndPeriod2]
}
# Warning for negative flows or NAs only in extended outputs
if (length(RunOptions$Outputs_Sim) > 2) {
if (any(OutputsModel$Qsim[!is.na(OutputsModel$Qsim)] < 0)) {
warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
}
# Warning for NAs
if (any(is.na(OutputsModel$Qsim))) {
warning(length(which(is.na(OutputsModel$Qsim))), " time steps with NA values")
}
}
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
SD <- lapply(seq(NbUpBasins), function(x) {
lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)]
InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x]
})
if (is.null(OutputsModel$StateEnd)) {
OutputsModel$StateEnd <- CreateIniStates(RunModel_Lag, InputsModel, SD = SD)
} else {
OutputsModel$StateEnd$SD <- SD
}
# message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", "))
}
if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
OutputsModel$RunOptions$WarmUpQsim <- Qsim_m3[seq_len(length(RunOptions$IndPeriod_WarmUp))] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
}
if ("Param" %in% RunOptions$Outputs_Sim) {
OutputsModel$RunOptions$Param <- c(Param, OutputsModel$RunOptions$Param)
}
class(OutputsModel) <- c(class(OutputsModel), "SD")
return(OutputsModel)
}
R/SeriesAggreg.InputsModel.R 0 → 100644
View file @ 4bccfee2
SeriesAggreg.InputsModel <- function(x, Format, ...) {
SeriesAggreg.list(x,
Format,
ConvertFun = NA,
except = c("ZLayers", "LengthHydro", "BasinAreas"),
...)
}
R/SeriesAggreg.OutputsModel.R 0 → 100644
View file @ 4bccfee2
SeriesAggreg.OutputsModel <- function(x, Format, ...) {
SeriesAggreg.list(x,
Format,
ConvertFun = NA,
except = c("RunOptions", "StateEnd"),
...)
}
R/SeriesAggreg.R
View file @ 4bccfee2
SeriesAggreg <- function(TabSeries, TimeFormat, NewTimeFormat, ConvertFun, YearFirstMonth = 1, TimeLag = 0, verbose = TRUE){ SeriesAggreg <- function(x, Format, ...) {
UseMethod("SeriesAggreg")
}
##_Arguments_check
##check_TabSeries
if(is.null(TabSeries) ){ stop("TabSeries must be a dataframe containing the dates and data to be converted \n"); return(NULL); }
if(!is.data.frame(TabSeries)){ stop("TabSeries must be a dataframe containing the dates and data to be converted \n"); return(NULL); }
if(ncol(TabSeries)<2){ stop("TabSeries must contain at least two columns (including the coulmn of dates \n"); return(NULL); }
##check_TimeFormat
if(!any(class(TabSeries[,1]) %in% "POSIXt")) {
stop("TabSeries first column must be a vector of class POSIXlt or POSIXct \n")
return(NULL)
}
if(any(class(TabSeries[,1]) %in% "POSIXlt")) {
TabSeries[,1] <- as.POSIXct(TabSeries[,1])
}
for(iCol in 2:ncol(TabSeries)){
if(!is.numeric(TabSeries[,iCol])){ stop("TabSeries columns (other than the first one) be of numeric class \n"); return(NULL); } }
if(is.null( TimeFormat)){ stop("TimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
if(!is.vector( TimeFormat)){ stop("TimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
if(!is.character(TimeFormat)){ stop("TimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
if(length(TimeFormat)!=1 ){ stop("TimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
if(TimeFormat %in% c("hourly","daily","monthly","yearly")==FALSE){
stop("TimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
##check_NewTimeFormat
if(is.null( NewTimeFormat)){ stop("NewTimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
if(!is.vector( NewTimeFormat)){ stop("NewTimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
if(!is.character(NewTimeFormat)){ stop("NewTimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
if(length(NewTimeFormat)!=1 ){ stop("NewTimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
if(NewTimeFormat %in% c("hourly","daily","monthly","yearly")==FALSE){
stop("NewTimeFormat must be 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
##check_ConvertFun
if(is.null( ConvertFun)){ stop("ConvertFun must be a vector of character \n"); return(NULL); }
if(!is.vector( ConvertFun)){ stop("ConvertFun must be a vector of character \n"); return(NULL); }
if(!is.character(ConvertFun)){ stop("ConvertFun must be a vector of character \n"); return(NULL); }
if(length(ConvertFun)!=(ncol(TabSeries)-1)){ stop(paste("ConvertFun must be of length ",ncol(TabSeries)-1," (length=ncol(TabSeries)-1) \n",sep="")); return(NULL); }
if(sum(ConvertFun %in% c("sum","mean")==FALSE)!=0){ stop("ConvertFun elements must be either 'sum' or 'mean' \n"); return(NULL); }
##check_YearFirstMonth
if(is.null( YearFirstMonth)){ stop("YearFirstMonth must be an integer between 1 and 12 \n"); return(NULL); }
if(!is.vector( YearFirstMonth)){ stop("YearFirstMonth must be an integer between 1 and 12 \n"); return(NULL); }
if(!is.numeric(YearFirstMonth)){ stop("YearFirstMonth must be an integer between 1 and 12 \n"); return(NULL); }
YearFirstMonth <- as.integer(YearFirstMonth);
if(length(YearFirstMonth)!=1){ stop(paste("YearFirstMonth must be only one integer between 1 and 12 \n",sep="")); return(NULL); }
if(YearFirstMonth %in% (1:12) == FALSE){ stop(paste("YearFirstMonth must be only one integer between 1 and 12 \n",sep="")); return(NULL); }
##check_DatesR_integrity
if(TimeFormat=="hourly" ){ by <- "hours" ; }
if(TimeFormat=="daily" ){ by <- "days" ; }
if(TimeFormat=="monthly"){ by <- "months"; }
if(TimeFormat=="yearly" ){ by <- "years" ; }
TmpDatesR <- seq(from=TabSeries[1,1],to=tail(TabSeries[,1],1),by=by)
if(!identical(TabSeries[,1],TmpDatesR)){ stop("Problem detected in TabSeries dates vector (in comparison with seq(from=TabSeries[1,1],to=tail(TabSeries[,1],1))) \n"); return(NULL); }
##check_conversion_direction
if((TimeFormat == "daily" & NewTimeFormat %in% c("hourly") ) |
(TimeFormat == "monthly" & NewTimeFormat %in% c("hourly","daily") ) |
(TimeFormat == "yearly" & NewTimeFormat %in% c("hourly","daily","monthly"))){
stop("only time aggregation can be performed \n"); return(NULL); }
##check_if_conversion_not_needed
if((TimeFormat == "hourly" & NewTimeFormat=="hourly" ) |
(TimeFormat == "daily" & NewTimeFormat=="daily" ) |
(TimeFormat == "monthly" & NewTimeFormat=="monthly") |
(TimeFormat == "yearly" & NewTimeFormat=="yearly" )){
if(verbose){ warning("\t The old and new format are identical \n\t -> no time-step conversion was performed \n"); return(TabSeries); } }
##_Time_step_conversion
##_Handle_conventional_difference_between_hourly_series_and_others
TmpDatesR <- TabSeries[,1];
#if(TimeFormat=="hourly"){ TmpDatesR <- TmpDatesR - 60*60; }
TmpDatesR <- TmpDatesR + TimeLag
Hmax <- "00"; if(TimeFormat=="hourly"){ Hmax <- "23" }
##_Identify_the_part_of_the_series_to_be_aggregated
NDaysInMonth <- list("31",c("28","29"),"31","30","31","30","31","31","30","31","30","31")
YearLastMonth <- YearFirstMonth+11; if(YearLastMonth>12){ YearLastMonth <- YearLastMonth-12; }
YearFirstMonthTxt <- formatC(YearFirstMonth,format="d",width=2,flag="0")
YearLastMonthTxt <- formatC(YearLastMonth,format="d",width=2,flag="0")
if(NewTimeFormat=="daily" ){ Ind1 <- which(format(TmpDatesR, "%H")=="00");
Ind2 <- which(format(TmpDatesR, "%H")==Hmax);
if(Ind2[1]<Ind1[1]){ Ind2 <- Ind2[2:length(Ind2)]; }
if(tail(Ind1,1)>tail(Ind2,1)){ Ind1 <- Ind1[1:(length(Ind1)-1)]; }
### Aggr <- NULL; iii <- 0; for(kkk in 1:length(Ind1)){
### iii <- iii+1; Aggr <- c(Aggr,rep(iii,length(Ind1[kkk]:Ind2[kkk]))); }
Aggr <- as.numeric(format(TmpDatesR[min(Ind1):max(Ind2)],"%Y%m%d")); ### more efficient
NewDatesR <- data.frame(seq(from=TmpDatesR[min(Ind1)],to=TmpDatesR[max(Ind2)],by="days"))
}
if(NewTimeFormat=="monthly"){ Ind1 <- which(format(TmpDatesR, "%d%H")=="0100");
Ind2 <- which(format(TmpDatesR,"%m%d%H") %in% paste(c("0131","0228","0229","0331","0430","0531","0630","0731","0831","0930","1031","1130","1231"),Hmax,sep=""));
Ind2[1:(length(Ind2)-1)][diff(Ind2)==1] <- NA; Ind2 <- Ind2[!is.na(Ind2)]; ### to keep only feb 29 if both feb 28 and feb 29 exists
if(Ind2[1]<Ind1[1]){ Ind2 <- Ind2[2:length(Ind2)]; }
if(tail(Ind1,1)>tail(Ind2,1)){ Ind1 <- Ind1[1:(length(Ind1)-1)]; }
### Aggr <- NULL; iii <- 0; for(kkk in 1:length(Ind1)){
### iii <- iii+1; Aggr <- c(Aggr,rep(iii,length(Ind1[kkk]:Ind2[kkk]))); }
Aggr <- as.numeric(format(TmpDatesR[min(Ind1):max(Ind2)],"%Y%m")); ### more efficient
NewDatesR <- data.frame(seq(from=TmpDatesR[min(Ind1)],to=TmpDatesR[max(Ind2)],by="months"))
}
if(NewTimeFormat=="yearly" ){ Ind1 <- which(format(TmpDatesR,"%m%d%H") %in% paste(YearFirstMonthTxt,"0100",sep=""));
Ind2 <- which(format(TmpDatesR,"%m%d%H") %in% paste(YearLastMonthTxt,NDaysInMonth[[YearLastMonth]],Hmax,sep=""));
Ind2[1:(length(Ind2)-1)][diff(Ind2)==1] <- NA; Ind2 <- Ind2[!is.na(Ind2)]; ### to keep only feb 29 if both feb 28 and feb 29 exists
if(Ind2[1]<Ind1[1]){ Ind2 <- Ind2[2:length(Ind2)]; }
if(tail(Ind1,1)>tail(Ind2,1)){ Ind1 <- Ind1[1:(length(Ind1)-1)]; }
Aggr <- NULL; iii <- 0; for(kkk in 1:length(Ind1)){
iii <- iii+1; Aggr <- c(Aggr,rep(iii,length(Ind1[kkk]:Ind2[kkk]))); }
### Aggr <- as.numeric(format(TmpDatesR[min(Ind1):max(Ind2)],"%Y")); ### not working if YearFirstMonth != 01
NewDatesR <- data.frame(seq(from=TmpDatesR[min(Ind1)],to=TmpDatesR[max(Ind2)],by="years"))
}
##_Aggreation_and_export
NewTabSeries <- data.frame(NewDatesR)
for(iCol in 2:ncol(TabSeries)){
AggregData <- aggregate(TabSeries[min(Ind1):max(Ind2),iCol],by=list(Aggr),FUN=ConvertFun[iCol-1],na.rm=F)[,2]
NewTabSeries <- data.frame(NewTabSeries,AggregData)
}
names(NewTabSeries) <- names(TabSeries)
return(NewTabSeries);
}
\ No newline at end of file
R/SeriesAggreg.data.frame.R 0 → 100644
View file @ 4bccfee2
SeriesAggreg.data.frame <- function(x,
Format,
ConvertFun,
TimeFormat = NULL, # deprecated
NewTimeFormat = NULL, # deprecated
YearFirstMonth = 1,
TimeLag = 0,
...) {
## Arguments checks
if (!is.null(TimeFormat)) {
warning("deprecated 'TimeFormat' argument", call. = FALSE)
}
if (missing(Format)) {
Format <- .GetSeriesAggregFormat(NewTimeFormat)
} else if (!is.null(NewTimeFormat)) {
warning("deprecated 'NewTimeFormat' argument: 'Format' argument is used instead",
call. = FALSE)
}
if (is.null(Format)) {
stop("argument 'Format' is missing")
}
## check x
if (!is.data.frame(x)) {
stop("'x' must be a data.frame containing the dates and data to be aggregated")
}
if (ncol(x) < 2) {
stop("'x' must contain at least two columns (including the column of dates)")
}
## check x date column
if (!inherits(x[[1L]], "POSIXt")) {
stop("'x' first column must be a vector of class 'POSIXlt' or 'POSIXct'")
}
if (inherits(x[[1L]], "POSIXlt")) {
x[[1L]] <- as.POSIXct(x[[1L]])
}
## check x other columns (boolean converted to numeric)
apply(x[, -1L, drop = FALSE],
MARGIN = 2,
FUN = function(iCol) {
if (!is.numeric(iCol)) {
stop("'x' columns (other than the first one) must be of numeric class")
}
})
## check Format
listFormat <- c("%Y%m%d", "%Y%m", "%Y", "%m", "%d")
Format <- gsub(pattern = "[[:punct:]]+", replacement = "%", Format)
Format <- match.arg(Format, choices = listFormat)
## check ConvertFun
if (length(ConvertFun) != (ncol(x) - 1)) {
stop(sprintf("'ConvertFun' must be of length %i (ncol(x)-1)", ncol(x) - 1))
}
listConvertFun <- lapply(unique(ConvertFun), function(y) {
if (!grepl("^q\\d+$", y, ignore.case = TRUE)) {
return(match.fun(y))
}
})
names(listConvertFun) <- unique(ConvertFun)
lapply(ConvertFun, function(y) {
if (!grepl("^q\\d+$", y, ignore.case = TRUE)) {
TestOutput <- listConvertFun[[y]](1:10)
if (!is.numeric(TestOutput)) {
stop(sprintf("Returned value of '%s' function should be numeric", y))
}
if (length(TestOutput) != 1) {
stop(sprintf("Returned value of '%s' function should be of length 1", y))
}
}
})
## check YearFirstMonth
msgYearFirstMonth <- "'YearFirstMonth' should be a single vector of numeric value between 1 and 12"
YearFirstMonth <- match(YearFirstMonth, 1:12)
if (anyNA(YearFirstMonth)) {
stop(msgYearFirstMonth)
}
if (length(YearFirstMonth) != 1) {
stop(msgYearFirstMonth)
}
if (YearFirstMonth != 1 & Format != "%Y") {
warning("'YearFirstMonth' is ignored because Format != '%Y'")
}
## check TimeLag
msgTimeLag <- "'TimeLag' should be a single vector of a positive numeric value"
if (!is.vector(TimeLag)) {
stop(msgTimeLag)
}
if (!is.numeric(TimeLag)) {
stop(msgTimeLag)
}
if (length(TimeLag) != 1 | !any(TimeLag >= 0)) {
stop(msgTimeLag)
}
TabSeries0 <- x
colnames(TabSeries0)[1L] <- "DatesR"
TabSeries0$DatesR <- TabSeries0$DatesR + TimeLag
TabSeries2 <- TabSeries0
if (!Format %in% c("%d", "%m")) {
start <- sprintf("%i-01-01 00:00:00",
as.numeric(format(TabSeries2$DatesR[1L], format = "%Y")) - 1)
stop <- sprintf("%i-12-31 00:00:00",
as.numeric(format(TabSeries2$DatesR[nrow(TabSeries2)], format = "%Y")) + 1)
unitTs <- format(diff(x[1:2, 1]))
if (gsub("[0-9]+ ", "", unitTs) == "hours") {
byTs <- "hours"
} else {
if (gsub(" days$", "", unitTs) == "1") {
byTs <- "days"
} else {
byTs <- "months"
}
}
fakeTs <- data.frame(DatesR = seq(from = as.POSIXct(start, tz = "UTC"),
to = as.POSIXct(stop , tz = "UTC"),
by = byTs) + TimeLag)
TabSeries2 <- merge(fakeTs, TabSeries2, by = "DatesR", all.x = TRUE)
}
TabSeries2$DatesRini <- TabSeries2$DatesR %in% TabSeries0$DatesR
TabSeries2$Selec2 <- format(TabSeries2$DatesR, Format)
if (nchar(Format) > 2 | Format == "%Y") {
# Compute aggregation
TabSeries2$Selec <- !duplicated(TabSeries2$Selec2)
if (all(TabSeries2$Selec)) {
warning("the requested time 'Format' is the same as the one in 'x'. No time-step conversion was performed")
return(x)
}
if (Format == "%Y") {
yfm <- sprintf("%02.f", YearFirstMonth)
spF1 <- "%m"
spF2 <- "%Y-%m"
TabSeries2$Selec1 <- format(TabSeries2$DatesR, spF1)
TabSeries2$Selec2 <- format(TabSeries2$DatesR, spF2)
TabSeries2$Selec <- !duplicated(TabSeries2$Selec2) & TabSeries2$Selec1 == yfm
}
TabSeries2$Fac2 <- cumsum(TabSeries2$Selec)
} else {
# Compute regime
if (Format == "%d") {
spF2 <- "%m-%d"
TabSeries2$Selec2 <- format(TabSeries2$DatesR, spF2)
}
TabSeries2$Fac2 <- TabSeries2$Selec2
TabSeries2$Selec <- !duplicated(TabSeries2$Selec2)
}
listTsAggreg <- lapply(names(listConvertFun), function(y) {
if (any(ConvertFun == y)) {
colTsAggreg <- c("Fac2", colnames(x)[-1L][ConvertFun == y])
if (grepl("^q\\d+$", y, ignore.case = TRUE)) {
probs <- as.numeric(gsub("^q", "", y, ignore.case = TRUE)) / 100
if (probs < 0 || probs > 1) {
stop("'Q...' format of argument 'ConvertFun' must be an integer between 0 and 100")
}
aggregate(. ~ Fac2,
data = TabSeries2[, colTsAggreg],
FUN = quantile,
na.action = na.pass,
probs = probs,
type = 8,
na.rm = TRUE)
} else {
aggregate(. ~ Fac2,
data = TabSeries2[, colTsAggreg],
FUN = listConvertFun[[y]],
na.action = na.pass)
}
} else {
NULL
}
})
listTsAggreg <- listTsAggreg[!sapply(listTsAggreg, is.null)]
tsAggreg <- do.call(cbind, listTsAggreg)
tsAggreg <- tsAggreg[, !duplicated(colnames(tsAggreg))]
tsAggreg <- merge(tsAggreg,
TabSeries2[, c("Fac2", "DatesR", "DatesRini", "Selec")],
by = "Fac2",
all.x = TRUE,
all.y = FALSE)
tsAggreg <- tsAggreg[tsAggreg$Selec & tsAggreg$DatesRini, ]
tsAggreg <- tsAggreg[, colnames(TabSeries0)]
tsAggreg <- tsAggreg[order(tsAggreg$DatesR), ] # reorder by date especially for regime time series
colnames(tsAggreg)[1L] <- colnames(x)[1L] # keep original column names
return(tsAggreg)
}
R/SeriesAggreg.list.R 0 → 100644
View file @ 4bccfee2
SeriesAggreg.list <- function(x,
Format,
ConvertFun,
NewTimeFormat = NULL,
simplify = FALSE,
except = NULL,
recursive = TRUE,
...) {
classIni <- class(x)
class(x) <- "list" # in order to avoid the use of '['.InputsModel' or '['.OutputsModel' when x[i] is used
if (missing(Format)) {
Format <- .GetSeriesAggregFormat(NewTimeFormat)
} else if (!is.null(NewTimeFormat)) {
warning("deprecated 'NewTimeFormat' argument: 'Format' argument is used instead",
call. = FALSE)
}
# Check ConvertFun
if (any(classIni %in% c("InputsModel", "OutputsModel"))) {
if (!all(is.na(ConvertFun))) {
warning("Argument 'ConvertFun' is ignored on 'InputsModel' and 'OutputsModel' objects")
}
} else if (length(ConvertFun) != 1) {
stop("Argument 'ConvertFun' must be of length 1 with 'list' object")
} else if (!is.character(ConvertFun)) {
stop("Argument 'ConvertFun' must be a character")
}
# Determination of DatesR
if (!is.null(x$DatesR)) {
if (!inherits(x$DatesR, "POSIXt")) {
stop("'x$DatesR' should be of class 'POSIXt'")
}
DatesR <- x$DatesR
} else {
# Auto-detection of POSIXt item in Tabseries
itemPOSIXt <- which(sapply(x, function(x) {
inherits(x, "POSIXt")
}, simplify = TRUE))[1]
if (is.na(itemPOSIXt)) {
stop("At least one item of argument 'x' should be of class 'POSIXt'")
}
warning("Item 'DatesR' not found in 'x' argument: the item ",
names(x)[itemPOSIXt],
" has been automatically chosen")
DatesR <- x[[names(x)[itemPOSIXt]]]
}
# Selection of numeric items for aggregation
numericCols <- names(which(sapply(x, inherits, "numeric")))
arrayCols <- names(which(sapply(x, inherits, "array")))
numericCols <- setdiff(numericCols, arrayCols)
if (!is.null(except)) {
if (!inherits(except, "character")) {
stop("Argument 'except' should be a 'character'")
}
numericCols <- setdiff(numericCols, except)
}
cols <- x[numericCols]
lengthCols <- sapply(cols, length, simplify = TRUE)
if (any(lengthCols != length(DatesR))) {
sErr <- paste0(names(lengthCols)[lengthCols != length(DatesR)],
" (", lengthCols[lengthCols != length(DatesR)], ")",
collapse = ", ")
warning("The length of the following `numeric` items in 'x' ",
"is different than the length of 'DatesR (",
length(DatesR),
")': they will be ignored in the aggregation: ",
sErr)
cols <- cols[lengthCols == length(DatesR)]
}
dfOut <- NULL
if (length(cols)) {
# Treating aggregation at root level
if (is.na(ConvertFun)) {
ConvertFun2 <- .GetAggregConvertFun(names(cols), Format)
} else {
ConvertFun2 <- rep(ConvertFun, length(cols))
}
dfOut <- SeriesAggreg(cbind(DatesR, as.data.frame(cols)),
Format,
...,
ConvertFun = ConvertFun2)
}
if (simplify) {
# Returns data.frame of numeric found in the first level of the list
return(dfOut)
} else {
res <- list()
# Convert aggregated data.frame into list
if (!is.null(dfOut)) {
res <- as.list(dfOut)
## To be consistent with InputsModel class and because plot.OutputsModel use the POSIXlt class
res$DatesR <- as.POSIXlt(res$DatesR)
}
# Exploration of embedded lists and data.frames
if (is.null(recursive) || recursive) {
listCols <- x[!names(x) %in% except]
listCols <- listCols[sapply(listCols, inherits, "list")]
dfCols <- x[sapply(x, inherits, "data.frame")]
dfCols <- c(dfCols, x[sapply(x, inherits, "matrix")])
listCols <- listCols[setdiff(names(listCols), names(dfCols))]
if (length(listCols) > 0) {
if (is.na(ConvertFun)) {
# Check for predefined ConvertFun for all sub-elements
listConvertFun <- .GetAggregConvertFun(names(listCols), Format)
}
# Run SeriesAggreg for each embedded list
listRes <- lapply(names(listCols), function(y) {
listCols[[y]]$DatesR <- DatesR
if (is.na(ConvertFun)) {
SeriesAggreg.list(listCols[[y]],
Format = Format,
recursive = NULL,
...,
ConvertFun = listConvertFun[y])
} else {
SeriesAggreg.list(listCols[[y]],
Format = Format,
recursive = NULL,
...)
}
})
names(listRes) <- names(listCols)
if (is.null(res$DatesR)) {
# Copy DatesR in top level list
res$DatesR <- listRes[[1]]$DatesR
}
# Remove DatesR in embedded lists
lapply(names(listRes), function(x) {
listRes[[x]]$DatesR <<- NULL
})
res <- c(res, listRes)
}
if (length(dfCols) > 0) {
# Processing matrix and dataframes
for (i in length(dfCols)) {
key <- names(dfCols)[i]
m <- dfCols[[i]]
if (nrow(m) != length(DatesR)) {
warning(
"The number of rows of the 'matrix' item ",
key, " (", nrow(m),
") is different than the length of 'DatesR ('", length(DatesR),
"), it will be ignored in the aggregation"
)
} else {
if (is.na(ConvertFun)) {
ConvertFun2 <- rep(.GetAggregConvertFun(key, Format), ncol(m))
} else {
ConvertFun2 <- rep(ConvertFun, ncol(m))
}
res[[key]] <- SeriesAggreg.data.frame(data.frame(DatesR, m),
Format = Format,
ConvertFun = ConvertFun2)
}
}
}
}
# Store elements that are not aggregated
res <- c(res, x[setdiff(names(x), names(res))])
class(res) <- gsub("hourly|daily|monthly|yearly",
.GetSeriesAggregClass(Format),
classIni)
return(res)
}
}
R/TransfoParam.R
View file @ 4bccfee2
TransfoParam <- function(ParamIn, Direction, FUN_TRANSFO) { TransfoParam <- function(ParamIn, Direction, FUN_TRANSFO) {
FUN_TRANSFO <- match.fun(FUN_TRANSFO)
return(FUN_TRANSFO(ParamIn, Direction)) return(FUN_TRANSFO(ParamIn, Direction))
} }
R/TransfoParam_CemaNeige.R
View file @ 4bccfee2
TransfoParam_CemaNeige <- function(ParamIn, Direction) { TransfoParam_CemaNeige <- function(ParamIn, Direction) {
NParam <- 2 ## number of model parameters
Bool <- is.matrix(ParamIn) NParam <- 2L
if (Bool == FALSE) {
ParamIn <- rbind(ParamIn)
}
## check arguments
isVecParamIn <- is.vector(ParamIn)
if (isVecParamIn) {
ParamIn <- matrix(ParamIn, nrow = 1)
}
if (!inherits(ParamIn, "matrix")) {
stop("'ParamIn' must be of class 'matrix'")
}
if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
}
if (ncol(ParamIn) != NParam) { if (ncol(ParamIn) != NParam) {
stop(paste( "the CemaNeige module requires ", NParam, " parameters \n", sep = "" )) stop(sprintf("the CemaNeige module requires %i parameters", NParam))
return(NULL)
} }
## transformation
if (Direction == "TR") { if (Direction == "TR") {
ParamOut <- ParamIn ParamOut <- ParamIn
ParamOut[, 1] <- (ParamIn[, 1] + 9.99) / 19.98 ### CemaNeige X1 (weighting coefficient for snow pack thermal state) ParamOut[, 1] <- (ParamIn[, 1] + 9.99) / 19.98 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- exp(ParamIn[, 2]) / 200 ### CemaNeige X2 (degree-day melt coefficient) ParamOut[, 2] <- exp(ParamIn[, 2]) / 200 ### CemaNeige X2 (degree-day melt coefficient)
} }
if (Direction == "RT") { if (Direction == "RT") {
ParamOut <- ParamIn ParamOut <- ParamIn
ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state) ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- log(ParamIn[, 2] * 200) ### CemaNeige X2 (degree-day melt coefficient) ParamOut[, 2] <- log(ParamIn[, 2] * 200) ### CemaNeige X2 (degree-day melt coefficient)
} }
if (Bool == FALSE) {
ParamOut <- ParamOut[1, ] ## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
} }
return(ParamOut) return(ParamOut)
......
R/TransfoParam_CemaNeigeHyst.R 0 → 100644
View file @ 4bccfee2
TransfoParam_CemaNeigeHyst <- function(ParamIn, Direction) {
## number of model parameters
NParam <- 4L
## check arguments
isVecParamIn <- is.vector(ParamIn)
if (isVecParamIn) {
ParamIn <- matrix(ParamIn, nrow = 1)
}
if (!inherits(ParamIn, "matrix")) {
stop("'ParamIn' must be of class 'matrix'")
}
if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
}
if (ncol(ParamIn) != NParam) {
stop(sprintf("the CemaNeige module with linear hysteresis requires %i parameters", NParam))
}
## transformation
if (Direction == "TR") {
ParamOut <- ParamIn
ParamOut[, 1] <- (ParamIn[, 1] + 9.99) / 19.98 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- exp(ParamIn[, 2]) / 200 ### CemaNeige X2 (degree-day melt coefficient)
ParamOut[, 3] <- (ParamIn[, 3] * 5) + 50 ### Hyst Gaccum
ParamOut[, 4] <- (ParamIn[, 4] / 19.98) + 0.5 ### Hyst CV
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- log(ParamIn[, 2] * 200) ### CemaNeige X2 (degree-day melt coefficient)
ParamOut[, 3] <- (ParamIn[, 3] - 50) / 5 ### Hyst Gaccum
ParamOut[, 4] <- (ParamIn[, 4] - 0.5) * 19.98 ### Hyst CV
}
## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
}
R/TransfoParam_GR1A.R
View file @ 4bccfee2
TransfoParam_GR1A <- function(ParamIn, Direction) { TransfoParam_GR1A <- function(ParamIn, Direction) {
NParam <- 1 ## number of model parameters
Bool <- is.matrix(ParamIn) NParam <- 1L
if (Bool == FALSE) {
ParamIn <- rbind(ParamIn)
}
## check arguments
isVecParamIn <- is.vector(ParamIn)
if (isVecParamIn) {
ParamIn <- matrix(ParamIn, nrow = 1)
}
if (!inherits(ParamIn, "matrix")) {
stop("'ParamIn' must be of class 'matrix'")
}
if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
}
if (ncol(ParamIn) != NParam) { if (ncol(ParamIn) != NParam) {
stop(paste("the GR1A model requires ", NParam, " parameters \n", sep = "")) stop(sprintf("the GR1A model requires %i parameters", NParam))
return(NULL)
} }
## transformation
if (Direction == "TR") { if (Direction == "TR") {
ParamOut <- (ParamIn + 10.0) / 8 ParamOut <- (ParamIn + 10.0) / 8
} }
if (Direction == "RT") { if (Direction == "RT") {
ParamOut <- ParamIn * 8 - 10.0 ParamOut <- ParamIn * 8 - 10.0
} }
if (Bool == FALSE) {
ParamOut <- ParamOut[1, ] ## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
} }
return(ParamOut) return(ParamOut)
} }
R/TransfoParam_GR2M.R
View file @ 4bccfee2
TransfoParam_GR2M <- function(ParamIn, Direction) { TransfoParam_GR2M <- function(ParamIn, Direction) {
NParam <- 2 ## number of model parameters
Bool <- is.matrix(ParamIn) NParam <- 2L
if (Bool == FALSE) {
ParamIn <- rbind(ParamIn) ## check arguments
isVecParamIn <- is.vector(ParamIn)
if (isVecParamIn) {
ParamIn <- matrix(ParamIn, nrow = 1)
}
if (!inherits(ParamIn, "matrix")) {
stop("'ParamIn' must be of class 'matrix'")
}
if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
} }
if (ncol(ParamIn) != NParam) { if (ncol(ParamIn) != NParam) {
stop(paste("the GR2M model requires ", NParam, " parameters \n", sep = "")) stop(sprintf("the GR2M model requires %i parameters", NParam))
return(NULL)
} }
## transformation
if (Direction == "TR") { if (Direction == "TR") {
ParamOut <- ParamIn ParamOut <- ParamIn
ParamOut[, 1] <- exp(ParamIn[, 1]) ParamOut[, 1] <- exp(ParamIn[, 1])
ParamOut[, 2] <- ParamIn[, 2] / 4 + 2.5 ParamOut[, 2] <- ParamIn[, 2] / 4 + 2.5
} }
if (Direction == "RT") { if (Direction == "RT") {
ParamOut <- ParamIn ParamOut <- ParamIn
ParamOut[, 1] <- log(ParamIn[, 1]) ParamOut[, 1] <- log(ParamIn[, 1])
ParamOut[, 2] <- (ParamIn[, 2] - 2.5) * 4 ParamOut[, 2] <- (ParamIn[, 2] - 2.5) * 4
} }
if (Bool == FALSE) {
ParamOut <- ParamOut[1, ] ## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
} }
return(ParamOut) return(ParamOut)
} }
R/TransfoParam_GR4H.R
View file @ 4bccfee2
TransfoParam_GR4H <- function(ParamIn,Direction){ TransfoParam_GR4H <- function(ParamIn, Direction) {
NParam <- 4; ## number of model parameters
Bool <- is.matrix(ParamIn); NParam <- 4L
if(Bool==FALSE){ ParamIn <- rbind(ParamIn); }
if(ncol(ParamIn)!=NParam){ stop(paste("the GR4H model requires ",NParam," parameters \n",sep="")); return(NULL); }
## check arguments
if(Direction=="TR"){ isVecParamIn <- is.vector(ParamIn)
ParamOut <- ParamIn; if (isVecParamIn) {
ParamOut[,1] <- exp(ParamIn[,1]); ### GR4H X1 (production store capacity) ParamIn <- matrix(ParamIn, nrow = 1)
ParamOut[,2] <- sinh(ParamIn[, 2] / 3); ### GR4H X2 (groundwater exchange coefficient) }
ParamOut[,3] <- exp(ParamIn[,3]); ### GR4H X3 (routing store capacity) if (!inherits(ParamIn, "matrix")) {
ParamOut[,4] <- 480 + (480 - 0.5) * (ParamIn[,4] - 9.99) / 19.98; ### GR4H X4 (unit hydrograph time constant) stop("'ParamIn' must be of class 'matrix'")
} }
if(Direction=="RT"){ if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
ParamOut <- ParamIn; stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
ParamOut[,1] <- log(ParamIn[,1]); ### GR4H X1 (production store capacity) }
ParamOut[,2] <- 3 * asinh(ParamIn[,2]); ### GR4H X2 (groundwater exchange coefficient) if (ncol(ParamIn) != NParam) {
ParamOut[,3] <- log(ParamIn[,3]); ### GR4H X3 (routing store capacity) stop(sprintf("the GR4H model requires %i parameters", NParam))
ParamOut[,4] <- (ParamIn[,4] - 480) * 19.98 / (480 - 0.5) + 9.99; ### GR4H X4 (unit hydrograph time constant) }
## transformation
if (Direction == "TR") {
ParamOut <- ParamIn
ParamOut[, 1] <- exp(ParamIn[, 1]) ### GR4H X1 (production store capacity)
ParamOut[, 2] <- sinh(ParamIn[, 2] / 3) ### GR4H X2 (groundwater exchange coefficient)
ParamOut[, 3] <- exp(ParamIn[, 3]) ### GR4H X3 (routing store capacity)
ParamOut[, 4] <- 480 + (480 - 0.5) * (ParamIn[, 4] - 9.99) / 19.98 ### GR4H X4 (unit hydrograph time constant)
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- log(ParamIn[, 1]) ### GR4H X1 (production store capacity)
ParamOut[, 2] <- 3 * asinh(ParamIn[, 2]) ### GR4H X2 (groundwater exchange coefficient)
ParamOut[, 3] <- log(ParamIn[, 3]) ### GR4H X3 (routing store capacity)
ParamOut[, 4] <- (ParamIn[, 4] - 480) * 19.98 / (480 - 0.5) + 9.99 ### GR4H X4 (unit hydrograph time constant)
} }
if(Bool==FALSE){ ParamOut <- ParamOut[1,]; }
return(ParamOut); ## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
} }
R/TransfoParam_GR4J.R
View file @ 4bccfee2
TransfoParam_GR4J <- function(ParamIn,Direction){ TransfoParam_GR4J <- function(ParamIn, Direction) {
NParam <- 4;
Bool <- is.matrix(ParamIn); ## number of model parameters
if(Bool==FALSE){ ParamIn <- rbind(ParamIn); } NParam <- 4L
if(ncol(ParamIn)!=NParam){ stop(paste("the GR4J model requires ",NParam," parameters \n",sep="")); return(NULL); }
if(Direction=="TR"){ ## check arguments
ParamOut <- ParamIn; isVecParamIn <- is.vector(ParamIn)
ParamOut[,1] <- exp(ParamIn[,1]); ### GR4J X1 (production store capacity) if (isVecParamIn) {
ParamOut[,2] <- sinh(ParamIn[,2]); ### GR4J X2 (groundwater exchange coefficient) ParamIn <- matrix(ParamIn, nrow = 1)
ParamOut[,3] <- exp(ParamIn[,3]); ### GR4J X3 (routing store capacity) }
ParamOut[,4] <- 20+19.5*(ParamIn[,4]-9.99)/19.98; ### GR4J X4 (unit hydrograph time constant) if (!inherits(ParamIn, "matrix")) {
} stop("'ParamIn' must be of class 'matrix'")
if(Direction=="RT"){ }
ParamOut <- ParamIn; if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
ParamOut[,1] <- log(ParamIn[,1]); ### GR4J X1 (production store capacity) stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
ParamOut[,2] <- asinh(ParamIn[,2]); ### GR4J X2 (groundwater exchange coefficient) }
ParamOut[,3] <- log(ParamIn[,3]); ### GR4J X3 (routing store capacity) if (ncol(ParamIn) != NParam) {
ParamOut[,4] <- 9.99+19.98*(ParamIn[,4]-20)/19.5; ### GR4J X4 (unit hydrograph time constant) stop(sprintf("the GR4J model requires %i parameters", NParam))
}
## transformation
if (Direction == "TR") {
ParamOut <- ParamIn
ParamOut[, 1] <- exp(ParamIn[, 1]) ### GR4J X1 (production store capacity)
ParamOut[, 2] <- sinh(ParamIn[, 2]) ### GR4J X2 (groundwater exchange coefficient)
ParamOut[, 3] <- exp(ParamIn[, 3]) ### GR4J X3 (routing store capacity)
ParamOut[, 4] <- 20 + 19.5 * (ParamIn[, 4] - 9.99) / 19.98 ### GR4J X4 (unit hydrograph time constant)
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- log(ParamIn[, 1]) ### GR4J X1 (production store capacity)
ParamOut[, 2] <- asinh(ParamIn[, 2]) ### GR4J X2 (groundwater exchange coefficient)
ParamOut[, 3] <- log(ParamIn[, 3]) ### GR4J X3 (routing store capacity)
ParamOut[, 4] <- 9.99 + 19.98 * (ParamIn[, 4] - 20) / 19.5 ### GR4J X4 (unit hydrograph time constant)
} }
if(Bool==FALSE){ ParamOut <- ParamOut[1,]; }
return(ParamOut); ## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
} }
R/TransfoParam_GR5H.R 0 → 100644
View file @ 4bccfee2
TransfoParam_GR5H <- function(ParamIn, Direction) {
## number of model parameters
NParam <- 5L
## check arguments
isVecParamIn <- is.vector(ParamIn)
if (isVecParamIn) {
ParamIn <- matrix(ParamIn, nrow = 1)
}
if (!inherits(ParamIn, "matrix")) {
stop("'ParamIn' must be of class 'matrix'")
}
if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
}
if (ncol(ParamIn) != NParam) {
stop(sprintf("the GR5H model requires %i parameters", NParam))
}
## transformation
if (Direction == "TR") {
ParamOut <- ParamIn
ParamOut[, 1] <- exp(ParamIn[, 1]) ### GR5H X1 (production store capacity)
ParamOut[, 2] <- sinh(ParamIn[, 2]) ### GR5H X2 (groundwater exchange coefficient)
ParamOut[, 3] <- exp(ParamIn[, 3]) ### GR5H X3 (routing store capacity)
ParamOut[, 4] <- 480 + (480 - 0.01) * (ParamIn[, 4] - 10) / 20 ### GR5H X4 (unit hydrograph time constant)
ParamOut[, 5] <- (ParamIn[, 5] + 10) / 20 ### GR5H X5 (groundwater exchange coefficient 2)
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- log(ParamIn[, 1]) ### GR5H X1 (production store capacity)
ParamOut[, 2] <- asinh(ParamIn[, 2]) ### GR5H X2 (groundwater exchange coefficient)
ParamOut[, 3] <- log(ParamIn[, 3]) ### GR5H X3 (routing store capacity)
ParamOut[, 4] <- (ParamIn[, 4] - 480) * 20 / (480 - 0.01) + 10 ### GR5H X4 (unit hydrograph time constant)
ParamOut[, 5] <- ParamIn[, 5] * 20 - 10 ### GR5H X5 (groundwater exchange coefficient 2)
}
## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
}
R/TransfoParam_GR5J.R
View file @ 4bccfee2
TransfoParam_GR5J <- function(ParamIn,Direction){ TransfoParam_GR5J <- function(ParamIn, Direction) {
NParam <- 5;
Bool <- is.matrix(ParamIn); ## number of model parameters
if(Bool==FALSE){ ParamIn <- rbind(ParamIn); } NParam <- 5L
if(ncol(ParamIn)!=NParam){ stop(paste("the GR5J model requires ",NParam," parameters \n",sep="")); return(NULL); }
if(Direction=="TR"){ ## check arguments
ParamOut <- ParamIn; isVecParamIn <- is.vector(ParamIn)
ParamOut[,1] <- exp(ParamIn[,1]); ### GR5J X1 (production store capacity) if (isVecParamIn) {
ParamOut[,2] <- sinh(ParamIn[,2]); ### GR5J X2 (groundwater exchange coefficient 1) ParamIn <- matrix(ParamIn, nrow = 1)
ParamOut[,3] <- exp(ParamIn[,3]); ### GR5J X3 (routing store capacity) }
ParamOut[,4] <- 20+19.5*(ParamIn[,4]-9.99)/19.98; ### GR5J X4 (unit hydrograph time constant) if (!inherits(ParamIn, "matrix")) {
ParamOut[,5] <- (ParamIn[,5] + 9.99) / 19.98; ### GR5J X5 (groundwater exchange coefficient 2) stop("'ParamIn' must be of class 'matrix'")
} }
if(Direction=="RT"){ if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
ParamOut <- ParamIn; stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
ParamOut[,1] <- log(ParamIn[,1]); ### GR5J X1 (production store capacity) }
ParamOut[,2] <- asinh(ParamIn[,2]); ### GR5J X2 (groundwater exchange coefficient 1) if (ncol(ParamIn) != NParam) {
ParamOut[,3] <- log(ParamIn[,3]); ### GR5J X3 (routing store capacity) stop(sprintf("the GR5J model requires %i parameters", NParam))
ParamOut[,4] <- 9.99+19.98*(ParamIn[,4]-20)/19.5; ### GR5J X4 (unit hydrograph time constant) }
ParamOut[,5] <- ParamIn[,5] * 19.98 - 9.99; ### GR5J X5 (groundwater exchange coefficient 2)
}
## transformation
if(Bool==FALSE){ ParamOut <- ParamOut[1,]; } if (Direction == "TR") {
return(ParamOut); ParamOut <- ParamIn
ParamOut[, 1] <- exp(ParamIn[, 1]) ### GR5J X1 (production store capacity)
ParamOut[, 2] <- sinh(ParamIn[, 2]) ### GR5J X2 (groundwater exchange coefficient 1)
ParamOut[, 3] <- exp(ParamIn[, 3]) ### GR5J X3 (routing store capacity)
ParamOut[, 4] <- 20 + 19.5 * (ParamIn[, 4] - 9.99) / 19.98 ### GR5J X4 (unit hydrograph time constant)
ParamOut[, 5] <- (ParamIn[, 5] + 9.99) / 19.98 ### GR5J X5 (groundwater exchange coefficient 2)
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- log(ParamIn[, 1]) ### GR5J X1 (production store capacity)
ParamOut[, 2] <- asinh(ParamIn[, 2]) ### GR5J X2 (groundwater exchange coefficient 1)
ParamOut[, 3] <- log(ParamIn[, 3]) ### GR5J X3 (routing store capacity)
ParamOut[, 4] <- 9.99 + 19.98 * (ParamIn[, 4] - 20) / 19.5 ### GR5J X4 (unit hydrograph time constant)
ParamOut[, 5] <- ParamIn[, 5] * 19.98 - 9.99 ### GR5J X5 (groundwater exchange coefficient 2)
}
## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
} }
R/TransfoParam_GR6J.R
View file @ 4bccfee2
TransfoParam_GR6J <- function(ParamIn,Direction){ TransfoParam_GR6J <- function(ParamIn, Direction) {
NParam <- 6;
Bool <- is.matrix(ParamIn); ## number of model parameters
if(Bool==FALSE){ ParamIn <- rbind(ParamIn); } NParam <- 6L
if(ncol(ParamIn)!=NParam){ stop(paste("the GR6J model requires ",NParam," parameters \n",sep="")); return(NULL); }
if(Direction=="TR"){ ## check arguments
ParamOut <- ParamIn; isVecParamIn <- is.vector(ParamIn)
ParamOut[,1] <- exp(ParamIn[,1]); ### GR6J X1 (production store capacity) if (isVecParamIn) {
ParamOut[,2] <- sinh(ParamIn[,2]); ### GR6J X2 (groundwater exchange coefficient 1) ParamIn <- matrix(ParamIn, nrow = 1)
ParamOut[,3] <- exp(ParamIn[,3]); ### GR6J X3 (routing store capacity) }
ParamOut[,4] <- 20+19.5*(ParamIn[,4]-9.99)/19.98; ### GR6J X4 (unit hydrograph time constant) if (!inherits(ParamIn, "matrix")) {
ParamOut[,5] <- ParamIn[,5]/5.; ### GR5J X5 (groundwater exchange coefficient 2) stop("'ParamIn' must be of class 'matrix'")
ParamOut[,6] <- exp(ParamIn[,6]); ### GR6J X6 (coefficient for emptying exponential store) }
} if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
if(Direction=="RT"){ stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
ParamOut <- ParamIn; }
ParamOut[,1] <- log(ParamIn[,1]); ### GR6J X1 (production store capacity) if (ncol(ParamIn) != NParam) {
ParamOut[,2] <- asinh(ParamIn[,2]); ### GR6J X2 (groundwater exchange coefficient 1) stop(sprintf("the GR6J model requires %i parameters", NParam))
ParamOut[,3] <- log(ParamIn[,3]); ### GR6J X3 (routing store capacity) }
ParamOut[,4] <- 9.99+19.98*(ParamIn[,4]-20)/19.5; ### GR6J X4 (unit hydrograph time constant)
ParamOut[,5] <- ParamIn[,5] * 5.; ### GR5J X5 (groundwater exchange coefficient 2)
ParamOut[,6] <- log(ParamIn[,6]); ### GR6J X6 (coefficient for emptying exponential store) ## transformation
} if (Direction == "TR") {
ParamOut <- ParamIn
if(Bool==FALSE){ ParamOut <- ParamOut[1,]; } ParamOut[, 1] <- exp(ParamIn[, 1]) ### GR6J X1 (production store capacity)
return(ParamOut); ParamOut[, 2] <- sinh(ParamIn[, 2]) ### GR6J X2 (groundwater exchange coefficient 1)
ParamOut[, 3] <- exp(ParamIn[, 3]) ### GR6J X3 (routing store capacity)
ParamOut[, 4] <- 20 + 19.5 * (ParamIn[, 4] - 9.99) / 19.98 ### GR6J X4 (unit hydrograph time constant)
ParamOut[, 5] <- ParamIn[, 5] / 5 ### GR6J X5 (groundwater exchange coefficient 2)
ParamOut[, 6] <- exp(ParamIn[, 6]) ### GR6J X6 (coefficient for emptying exponential store)
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- log(ParamIn[, 1]) ### GR6J X1 (production store capacity)
ParamOut[, 2] <- asinh(ParamIn[, 2]) ### GR6J X2 (groundwater exchange coefficient 1)
ParamOut[, 3] <- log(ParamIn[, 3]) ### GR6J X3 (routing store capacity)
ParamOut[, 4] <- 9.99 + 19.98 * (ParamIn[, 4] - 20) / 19.5 ### GR6J X4 (unit hydrograph time constant)
ParamOut[, 5] <- ParamIn[, 5] * 5 ### GR6J X5 (groundwater exchange coefficient 2)
ParamOut[, 6] <- log(ParamIn[, 6]) ### GR6J X6 (coefficient for emptying exponential store)
}
## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
} }