An error occurred while loading the file. Please try again.
-
Delaigue Olivier authored
v1.6.1.11 fix: FUN_LAG function now takes a matrix in input in CreateCalibOptions and Calibration_Michel #34
9cb63cfd
Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD, FUN_CRIT, # deprecated FUN_TRANSFO = NULL, verbose = TRUE) { FUN_MOD <- match.fun(FUN_MOD) if (!missing(FUN_CRIT)) { FUN_CRIT <- match.fun(FUN_CRIT) } if (!is.null(FUN_TRANSFO)) { FUN_TRANSFO <- match.fun(FUN_TRANSFO) } ##_____Arguments_check_____________________________________________________________________ if (!inherits(InputsModel, "InputsModel")) { stop("'InputsModel' must be of class 'InputsModel'") } if (!inherits(RunOptions, "RunOptions")) { stop("'RunOptions' must be of class 'RunOptions'") } if (!inherits(InputsCrit, "InputsCrit")) { stop("'InputsCrit' must be of class 'InputsCrit'") } if (inherits(InputsCrit, "Multi")) { stop("'InputsCrit' must be of class 'Single' or 'Compo'") } if (inherits(InputsCrit, "Single")) { listVarObs <- InputsCrit$VarObs } if (inherits(InputsCrit, "Compo")) { listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs") } if ("SCA" %in% listVarObs & !"Gratio" %in% RunOptions$Outputs_Cal) { warning("Missing 'Gratio' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SCA") RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "Gratio") } if ("SWE" %in% listVarObs & !"SnowPack" %in% RunOptions$Outputs_Cal) { warning("Missing 'SnowPack' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SWE") RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "SnowPack") } if (!inherits(CalibOptions, "CalibOptions")) { stop("'CalibOptions' must be of class 'CalibOptions'") } if (!inherits(CalibOptions, "HBAN")) { stop("'CalibOptions' must be of class 'HBAN' if 'Calibration_Michel' is used") } if (!missing(FUN_CRIT)) { warning("argument 'FUN_CRIT' is deprecated. The error criterion function is now automatically get from the 'InputsCrit' object") } ##_check_FUN_TRANSFO if (is.null(FUN_TRANSFO)) { if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR4H)) { FUN1 <- TransfoParam_GR4H } if (identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) { FUN1 <- TransfoParam_GR5H } if (identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J)) { FUN1 <- TransfoParam_GR4J } if (identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) { FUN1 <- TransfoParam_GR5J7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
}
if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
FUN1 <- TransfoParam_GR6J
}
if (identical(FUN_MOD, RunModel_GR2M )) {
FUN1 <- TransfoParam_GR2M
}
if (identical(FUN_MOD, RunModel_GR1A )) {
FUN1 <- TransfoParam_GR1A
}
##_set_FUN_LAG
if (inherits(CalibOptions, "SD")) {
FUN_LAG <- TransfoParam_LAG
}
if (identical(FUN_MOD, RunModel_CemaNeige )) {
if (inherits(CalibOptions, "hysteresis")) {
FUN_TRANSFO <- TransfoParam_CemaNeigeHyst
} else {
FUN_TRANSFO <- TransfoParam_CemaNeige
}
}
if (identical(FUN_MOD, RunModel_GR4H) | identical(FUN_MOD, RunModel_GR5H) |
identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_GR6J)) {
if (!inherits(CalibOptions, "SD")) {
FUN_TRANSFO <- FUN1
} else {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
ParamOut[, 1:(NParam-1)] <- FUN1(ParamIn[, 1:(NParam-1)], Direction)
ParamOut[, NParam ] <- FUN_LAG(as.matrix(ParamIn[, NParam ]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
}
if (identical(FUN_MOD, RunModel_CemaNeigeGR4H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H) |
identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
if (inherits(CalibOptions, "hysteresis")) {
FUN2 <- TransfoParam_CemaNeigeHyst
} else {
FUN2 <- TransfoParam_CemaNeige
}
if (inherits(CalibOptions, "hysteresis") & !inherits(CalibOptions, "SD")) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
ParamOut[, 1:(NParam-4)] <- FUN1(ParamIn[, 1:(NParam-4)], Direction)
ParamOut[, (NParam-3):NParam ] <- FUN2(ParamIn[, (NParam-3):NParam ], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
if (!inherits(CalibOptions, "hysteresis") & !inherits(CalibOptions, "SD")) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
ParamOut[, 1:(NParam-2)] <- FUN1(ParamIn[, 1:(NParam-2)], Direction)
ParamOut[, (NParam-1):NParam ] <- FUN2(ParamIn[, (NParam-1):NParam ], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
if (inherits(CalibOptions, "hysteresis") & inherits(CalibOptions, "SD")) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
ParamOut[, 1:(NParam-5)] <- FUN1( ParamIn[, 1:(NParam-5)], Direction)
ParamOut[, (NParam-4):(NParam-1)] <- FUN2( ParamIn[, (NParam-4):(NParam-1)], Direction)
ParamOut[, NParam ] <- FUN_LAG(as.matrix(ParamIn[, NParam ]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
if (!inherits(CalibOptions, "hysteresis") & inherits(CalibOptions, "SD")) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
ParamOut[, 1:(NParam-3)] <- FUN1( ParamIn[, 1:(NParam-3)], Direction)
ParamOut[, (NParam-2):(NParam-1)] <- FUN2( ParamIn[, (NParam-2):(NParam-1)], Direction)
ParamOut[, NParam ] <- FUN_LAG(as.matrix(ParamIn[, NParam ]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
}
if (is.null(FUN_TRANSFO)) {
stop("'FUN_TRANSFO' was not found (in 'Calibration' function)")
}
}
##_variables_initialisation
ParamFinalR <- NULL
ParamFinalT <- NULL
CritFinal <- NULL
NRuns <- 0
NIter <- 0
if ("StartParamDistrib" %in% names(CalibOptions)) {
PrefilteringType <- 2
} else {
PrefilteringType <- 1
}
if (PrefilteringType == 1) {
NParam <- ncol(CalibOptions$StartParamList)
}
if (PrefilteringType == 2) {
NParam <- ncol(CalibOptions$StartParamDistrib)
}
if (NParam > 20) {
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
stop("Calibration_Michel can handle a maximum of 20 parameters")
}
HistParamR <- matrix(NA, nrow = 500 * NParam, ncol = NParam)
HistParamT <- matrix(NA, nrow = 500 * NParam, ncol = NParam)
HistCrit <- matrix(NA, nrow = 500 * NParam, ncol = 1)
CritName <- NULL
CritBestValue <- NULL
Multiplier <- NULL
CritOptim <- +1e100
##_temporary_change_of_Outputs_Sim
RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal ### this reduces the size of the matrix exchange with fortran and therefore speeds the calibration
##_____Parameter_Grid_Screening____________________________________________________________
##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter
ProposeCandidatesGrid <- function(DistribParam) {
NewCandidates <- expand.grid(lapply(seq_len(ncol(DistribParamR)), function(x) DistribParam[, x]))
NewCandidates <- unique(NewCandidates) # to avoid duplicates when a parameter is set
Output <- list(NewCandidates = NewCandidates)
}
##Creation_of_new_candidates_______________________________________________
OptimParam <- is.na(CalibOptions$FixedParam)
if (PrefilteringType == 1) {
CandidatesParamR <- CalibOptions$StartParamList
}
if (PrefilteringType == 2) {
DistribParamR <- CalibOptions$StartParamDistrib
DistribParamR[, !OptimParam] <- NA
CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates
}
##Remplacement_of_non_optimised_values_____________________________________
CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
return(x)
})
if (NParam > 1) {
CandidatesParamR <- t(CandidatesParamR)
} else {
CandidatesParamR <- cbind(CandidatesParamR)
}
##Loop_to_test_the_various_candidates______________________________________
iNewOptim <- 0
Ncandidates <- nrow(CandidatesParamR)
if (verbose & Ncandidates > 1) {
if (PrefilteringType == 1) {
message("List-Screening in progress (", appendLF = FALSE)
}
if (PrefilteringType == 2) {
message("Grid-Screening in progress (", appendLF = FALSE)
}
message("0%", appendLF = FALSE)
}
for (iNew in 1:nrow(CandidatesParamR)) {
if (verbose & Ncandidates > 1) {
for (k in c(2, 4, 6, 8)) {
if (iNew == round(k / 10 * Ncandidates)) {
message(" ", 10 * k, "%", appendLF = FALSE)
}
}
}
##Model_run
Param <- CandidatesParamR[iNew, ]
OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD)