Skip to content
GitLab
Commit bb03e2fb authored by Dorchies David's avatar Dorchies David
Browse files

feat(SD): Debug LLR propagation model 

Saved work 11 May 2022

Refs HYCAR-Hydro/airgr#152, HYCAR-Hydro/airgr#153
parent a45508d7
dev
No related merge requests found
Pipeline #36194 failed with stages
in 6 minutes and 12 seconds
Hide whitespace changes
Inline Side-by-side
Showing
with 160 additions and 183 deletions
+160 -183
NAMESPACE
+ 81
− 79
View file @ bb03e2fb
#####################################
## Load DLL ##
#####################################
useDynLib(airGR, .registration = TRUE)
#####################################
## S3 methods ##
#####################################
S3method('[', InputsModel)
#S3method('[', OutputsModel) ### to add in version 2.0
S3method(plot, OutputsModel)
S3method(SeriesAggreg, data.frame)
S3method(SeriesAggreg, list)
S3method(SeriesAggreg, InputsModel)
S3method(SeriesAggreg, OutputsModel)
#####################################
## Export ##
#####################################
export(Calibration)
export(Calibration_Michel)
export(CreateCalibOptions)
export(CreateErrorCrit_GAPX)
export(CreateIniStates)
export(CreateInputsCrit)
export(CreateInputsCrit_Lavenne)
export(CreateInputsModel)
export(CreateRunOptions)
export(DataAltiExtrapolation_Valery)
export(ErrorCrit)
export(ErrorCrit_KGE)
export(ErrorCrit_KGE2)
export(ErrorCrit_NSE)
export(ErrorCrit_RMSE)
export(Imax)
export(PE_Oudin)
export(plot.OutputsModel) ### to remove from version 2.0
export(RunModel)
export(RunModel_CemaNeige)
export(RunModel_CemaNeigeGR4H)
export(RunModel_CemaNeigeGR5H)
export(RunModel_CemaNeigeGR4J)
export(RunModel_CemaNeigeGR5J)
export(RunModel_CemaNeigeGR6J)
export(RunModel_GR1A)
export(RunModel_GR2M)
export(RunModel_GR4H)
export(RunModel_GR5H)
export(RunModel_GR4J)
export(RunModel_GR5J)
export(RunModel_GR6J)
export(RunModel_Lag)
export(SeriesAggreg)
export(TransfoParam)
export(TransfoParam_CemaNeige)
export(TransfoParam_CemaNeigeHyst)
export(TransfoParam_GR1A)
export(TransfoParam_GR2M)
export(TransfoParam_GR4H)
export(TransfoParam_GR5H)
export(TransfoParam_GR4J)
export(TransfoParam_GR5J)
export(TransfoParam_GR6J)
export(TransfoParam_Lag)
export(.ErrorCrit)
export(.FeatModels)
#####################################
## Import ##
#####################################
import(stats)
import(graphics)
import(grDevices)
import(utils)
#####################################
## Load DLL ##
#####################################
useDynLib(airGR, .registration = TRUE)
#####################################
## S3 methods ##
#####################################
S3method('[', InputsModel)
#S3method('[', OutputsModel) ### to add in version 2.0
S3method(plot, OutputsModel)
S3method(SeriesAggreg, data.frame)
S3method(SeriesAggreg, list)
S3method(SeriesAggreg, InputsModel)
S3method(SeriesAggreg, OutputsModel)
#####################################
## Export ##
#####################################
export(Calibration)
export(Calibration_Michel)
export(CreateCalibOptions)
export(CreateErrorCrit_GAPX)
export(CreateIniStates)
export(CreateInputsCrit)
export(CreateInputsCrit_Lavenne)
export(CreateInputsModel)
export(CreateRunOptions)
export(DataAltiExtrapolation_Valery)
export(ErrorCrit)
export(ErrorCrit_KGE)
export(ErrorCrit_KGE2)
export(ErrorCrit_NSE)
export(ErrorCrit_RMSE)
export(Imax)
export(PE_Oudin)
export(plot.OutputsModel) ### to remove from version 2.0
export(RunModel)
export(RunModel_CemaNeige)
export(RunModel_CemaNeigeGR4H)
export(RunModel_CemaNeigeGR5H)
export(RunModel_CemaNeigeGR4J)
export(RunModel_CemaNeigeGR5J)
export(RunModel_CemaNeigeGR6J)
export(RunModel_GR1A)
export(RunModel_GR2M)
export(RunModel_GR4H)
export(RunModel_GR5H)
export(RunModel_GR4J)
export(RunModel_GR5J)
export(RunModel_GR6J)
export(RunModel_Lag)
export(RunModel_LLR)
export(SeriesAggreg)
export(TransfoParam)
export(TransfoParam_CemaNeige)
export(TransfoParam_CemaNeigeHyst)
export(TransfoParam_GR1A)
export(TransfoParam_GR2M)
export(TransfoParam_GR4H)
export(TransfoParam_GR5H)
export(TransfoParam_GR4J)
export(TransfoParam_GR5J)
export(TransfoParam_GR6J)
export(TransfoParam_Lag)
export(TransfoParam_LLR)
export(.ErrorCrit)
export(.FeatModels)
#####################################
## Import ##
#####################################
import(stats)
import(graphics)
import(grDevices)
import(utils)
R/Calibration_Michel.R
+ 2
− 2
View file @ bb03e2fb
......@@ -296,7 +296,7 @@ Calibration_Michel <- function(InputsModel,
for (iNew in 1:nrow(CandidatesParamR)) {
##Model_run
Param <- CandidatesParamR[iNew, ]
OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, FUN_SD = RunOptions$FUN_SD, ...)
##Calibration_criterion_computation
OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
if (!is.na(OutputsCrit$CritValue)) {
......@@ -357,7 +357,7 @@ Calibration_Michel <- function(InputsModel,
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
##Model_run
Param <- CandidatesParamR[iNew, ]
OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, FUN_SD = RunOptions$FUN_SD, ...)
##Calibration_criterion_computation
OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
......
R/CreateCalibOptions.R
+ 1
− 1
View file @ bb03e2fb
......@@ -205,7 +205,7 @@ CreateCalibOptions <- function(FUN_MOD,
ParamT <- cbind(ParamTSD, ParamT)
}
if (identical(RunModel_LLR, FUN_SD)){
ParamTSD <- matrix(c(+1.25, +1.25
ParamTSD <- matrix(c(+1.25, +1.25,
+2.50, +2.50,
+5.00, +5.00), ncol = 2, byrow = TRUE)
ParamT <- cbind(ParamTSD, ParamT)
......
R/RunModel.R
+ 1
− 1
View file @ bb03e2fb
......@@ -14,7 +14,7 @@ RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, FUN_SD = RunModel_
Param = Param[iFirstParamRunOffModel:length(Param)], ...)
if (inherits(InputsModel, "SD") && !identical(FUN_MOD, FUN_SD)) {
OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1:RunOptions$FeatFUN_MOD$NParamSD],
OutputsModel <- FUN_SD(InputsModel, RunOptions, Param[1:RunOptions$FeatFUN_MOD$NParamSD],
OutputsModel)
}
return(OutputsModel)
......
R/RunModel_LLR.R
+ 20
− 53
View file @ bb03e2fb
......@@ -63,71 +63,38 @@ RunModel_LLR <- function(InputsModel, RunOptions, Param, QcontribDown) {
QcontribDown)
}
## Parameters set up
TParam <- Param[1L]
KParam <- Param[2L]
C0 <- exp(-1/KParam)
C1 <- ((KParam/1) * (1-C0))-C0
C2 <- 1 - (KParam/1) * (1-C0)
## propagation time from upstream meshes to outlet
PT <- InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep
print(PT)
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 = ", "))
PT <- round(InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep)
PK <- sqrt(InputsModel$LengthHydro/max(InputsModel$LengthHydro)) * KParam
C0 <- exp(-1/PK)
C1 <- ((PK/1) * (1-C0))-C0
C2 <- 1 - (PK/1) * (1-C0)
## Lag model computation
Qsim_m3 <- QsimDown *
InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
Qsim_m3 <- matrix(QsimDown *
InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
, nrow = length(IndPeriod1), ncol = NbUpBasins)
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]
Qupstream_m3 <- InputsModel$Qupstream[IndPeriod1, upstream_basin]
Qroute <- Qllr <- Qupstream_m3
for (q_upstream in seq_along(Qupstream_m3)[2:max(seq_along(Qupstream_m3))]){
Qroute[q_upstream] <- C0 * Qroute[q_upstream - 1] + C1 * Qupstream_m3[q_upstream - 1] + C2 * Qupstream_m3[q_upstream]
Qllr[q_upstream + PT[upstream_basin]] <- Qroute[q_upstream]
}
Qsim_m3[, upstream_basin] <- Qllr[IndPeriod1]
}
Qsim_m3 <- rowSums(Qsim_m3)
## OutputsModel
if ("Qsim_m3" %in% RunOptions$Outputs_Sim) {
......
R/TransfoParam_LLR.R
+ 3
− 0
View file @ bb03e2fb
......@@ -2,6 +2,7 @@ TransfoParam_LLR <- function(ParamIn, Direction) {
## number of model parameters
NParam <- 2L
## check arguments
isVecParamIn <- is.vector(ParamIn)
if (isVecParamIn) {
......@@ -20,10 +21,12 @@ TransfoParam_LLR <- function(ParamIn, Direction) {
## transformation
if (Direction == "TR") {
ParamOut <- ParamIn
ParamOut[, 1] <- 20 * (ParamIn[, 1] + 10) / 20.0
ParamOut[, 2] <- 20 * (ParamIn[, 2] + 10) / 20.0
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- ParamIn[, 1] * 20.0 / 20 - 10
ParamOut[, 2] <- 20 * (ParamIn[, 2] + 10) / 20.0
}
......
R/TransfoParam_Lag.R
+ 38
− 39
View file @ bb03e2fb
TransfoParam_Lag <- function(ParamIn, Direction) {
## number of model parameters
NParam <- 1L
## 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 Lag model requires %i parameters", NParam))
}
## transformation
if (Direction == "TR") {
ParamOut <- 20 * (ParamIn + 10) / 20.0
}
if (Direction == "RT") {
ParamOut <- ParamIn * 20.0 / 20 - 10
}
## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
}
TransfoParam_Lag <- function(ParamIn, Direction) {
## number of model parameters
NParam <- 1L
## 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 Lag model requires %i parameters", NParam))
}
## transformation
if (Direction == "TR") {
ParamOut <- 20 * (ParamIn + 10) / 20.0
}
if (Direction == "RT") {
ParamOut <- ParamIn * 20.0 / 20 - 10
}
## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
}
R/Utils.R
+ 3
− 2
View file @ bb03e2fb
......@@ -76,9 +76,10 @@
if (!is.null(FUN_SD)){
if (identical(RunModel_Lag, FUN_SD)){
res$NParamSD <- FeatMod$NbParam[14]
}
if(identical(RunModel_LLR, FUN_SD)){
res$FUN_SD <- "Lag"
}else if(identical(RunModel_LLR, FUN_SD)){
res$NParamSD <- FeatMod$NbParam[15]
res$FUN_SD <- "LLR"
}
}
return(res)
......
R/UtilsCalibOptions.R
+ 11
− 6
View file @ bb03e2fb
......@@ -29,13 +29,11 @@
if (IsSD) {
if (identical(RunModel_Lag, FUN_SD)){
FUN_ROUT <- TransfoParam_Lag
}
if (identical(RunModel_LLR, FUN_SD)){
} else if (identical(RunModel_LLR, FUN_SD)){
FUN_ROUT <- TransfoParam_LLR
}
}
## set FUN_TRANSFO
if (! "CemaNeige" %in% FeatFUN_MOD$Class) {
if (!IsSD) {
......@@ -46,10 +44,17 @@
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
ParamOut[, (FeatFUN_MOD$NParamSD + 1):NParam] <- FUN_GR(ParamIn[, (FeatFUN_MOD$NParamSD + 1):NParam], Direction)
ParamOut[, FeatFUN_MOD$NParamSD ] <- FUN_ROUT(as.matrix(ParamIn[, FeatFUN_MOD$NParamSD]), Direction)
if (identical(RunModel_Lag, FUN_SD)){
ParamOut[, 1 ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
} else if (identical(RunModel_LLR, FUN_SD)){
ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
}
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
......@@ -103,7 +108,7 @@
NParam <- ncol(ParamIn)
ParamOut[, 2:(NParam - 4) ] <- FUN_GR(ParamIn[, 2:(NParam - 4)], Direction)
ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
ParamOut[, 1 ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
......@@ -124,7 +129,7 @@
ParamOut[, 2:(NParam - 2)] <- FUN_GR(ParamIn[, 2:(NParam - 2)], Direction)
}
ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
ParamOut[, 1 ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment