Commit f7323d49 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

v1.0.13.9 minor syntax revisions in Calibration_Michel

parent 1f275390
Package: airGR Package: airGR
Type: Package Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.0.13.8 Version: 1.0.13.9
Date: 2018-08-30 Date: 2018-08-31
Authors@R: c( Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")), person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")), person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")),
......
...@@ -14,7 +14,7 @@ output: ...@@ -14,7 +14,7 @@ output:
### 1.0.13.8 Release Notes (2018-08-30) ### 1.0.13.8 Release Notes (2018-08-31)
#### Deprectated and defunct #### Deprectated and defunct
......
Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE) { Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions,
FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE) {
##_____Arguments_check_____________________________________________________________________
if (!inherits(InputsModel, "InputsModel")) { ##_____Arguments_check_____________________________________________________________________
stop("InputsModel must be of class 'InputsModel' \n") if (!inherits(InputsModel, "InputsModel")) {
return(NULL) stop("InputsModel must be of class 'InputsModel' \n")
} return(NULL)
if (!inherits(RunOptions, "RunOptions")) { }
stop("RunOptions must be of class 'RunOptions' \n") if (!inherits(RunOptions, "RunOptions")) {
return(NULL) stop("RunOptions must be of class 'RunOptions' \n")
} return(NULL)
if (!inherits(InputsCrit, "InputsCrit")) { }
stop("InputsCrit must be of class 'InputsCrit' \n") if (!inherits(InputsCrit, "InputsCrit")) {
return(NULL) stop("InputsCrit must be of class 'InputsCrit' \n")
} return(NULL)
if (!inherits(CalibOptions, "CalibOptions")) { }
stop("CalibOptions must be of class 'CalibOptions' \n") if (!inherits(CalibOptions, "CalibOptions")) {
return(NULL) stop("CalibOptions must be of class 'CalibOptions' \n")
} return(NULL)
if (!inherits(CalibOptions, "HBAN")) { }
stop("CalibOptions must be of class 'HBAN' if Calibration_Michel is used \n") if (!inherits(CalibOptions, "HBAN")) {
return(NULL) stop("CalibOptions must be of class 'HBAN' if Calibration_Michel is used \n")
} return(NULL)
}
##_check_FUN_TRANSFO
if (is.null(FUN_TRANSFO)) { ##_check_FUN_TRANSFO
if (identical(FUN_MOD, RunModel_GR4H )) { if (is.null(FUN_TRANSFO)) {
FUN_TRANSFO <- TransfoParam_GR4H if (identical(FUN_MOD, RunModel_GR4H )) {
} FUN_TRANSFO <- TransfoParam_GR4H
if (identical(FUN_MOD, RunModel_GR4J )) { }
FUN_TRANSFO <- TransfoParam_GR4J if (identical(FUN_MOD, RunModel_GR4J )) {
} FUN_TRANSFO <- TransfoParam_GR4J
if (identical(FUN_MOD, RunModel_GR5J )) { }
FUN_TRANSFO <- TransfoParam_GR5J if (identical(FUN_MOD, RunModel_GR5J )) {
} FUN_TRANSFO <- TransfoParam_GR5J
if (identical(FUN_MOD, RunModel_GR6J )) { }
FUN_TRANSFO <- TransfoParam_GR6J if (identical(FUN_MOD, RunModel_GR6J )) {
} FUN_TRANSFO <- TransfoParam_GR6J
if (identical(FUN_MOD, RunModel_GR2M )) { }
FUN_TRANSFO <- TransfoParam_GR2M if (identical(FUN_MOD, RunModel_GR2M )) {
FUN_TRANSFO <- TransfoParam_GR2M
}
if (identical(FUN_MOD, RunModel_GR1A )) {
FUN_TRANSFO <- TransfoParam_GR1A
}
if (identical(FUN_MOD, RunModel_CemaNeige )) {
FUN_TRANSFO <- TransfoParam_CemaNeige
}
if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
if (identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
FUN1 <- TransfoParam_GR4J
FUN2 <- TransfoParam_CemaNeige
} }
if (identical(FUN_MOD, RunModel_GR1A )) { if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
FUN_TRANSFO <- TransfoParam_GR1A FUN1 <- TransfoParam_GR5J
FUN2 <- TransfoParam_CemaNeige
} }
if (identical(FUN_MOD, RunModel_CemaNeige )) { if (identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
FUN_TRANSFO <- TransfoParam_CemaNeige FUN1 <- TransfoParam_GR6J
FUN2 <- TransfoParam_CemaNeige
} }
if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) { FUN_TRANSFO <- function(ParamIn, Direction) {
if (identical(FUN_MOD, RunModel_CemaNeigeGR4J)) { Bool <- is.matrix(ParamIn)
FUN1 <- TransfoParam_GR4J if (Bool == FALSE) {
FUN2 <- TransfoParam_CemaNeige ParamIn <- rbind(ParamIn)
}
if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
FUN1 <- TransfoParam_GR5J
FUN2 <- TransfoParam_CemaNeige
} }
if (identical(FUN_MOD,RunModel_CemaNeigeGR6J)) { ParamOut <- NA * ParamIn
FUN1 <- TransfoParam_GR6J NParam <- ncol(ParamIn)
FUN2 <- TransfoParam_CemaNeige ParamOut[, 1:(NParam-2)] <- FUN1(ParamIn[, 1:(NParam-2)], Direction)
} ParamOut[, (NParam-1):NParam ] <- FUN2(ParamIn[, (NParam-1):NParam ], Direction)
FUN_TRANSFO <- function(ParamIn, Direction) { if (Bool == FALSE) {
Bool <- is.matrix(ParamIn) ParamOut <- ParamOut[1, ]
if (Bool == FALSE) {
ParamIn <- rbind(ParamIn)
}
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 == FALSE) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
} }
return(ParamOut)
} }
if (is.null(FUN_TRANSFO)) {
stop("FUN_TRANSFO was not found (in Calibration function) \n")
return(NULL)
}
}
##_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) { if (is.null(FUN_TRANSFO)) {
stop("Calibration_Michel can handle a maximum of 20 parameters \n") stop("FUN_TRANSFO was not found (in Calibration function) \n")
return(NULL) return(NULL)
} }
HistParamR <- matrix(NA, nrow = 500*NParam, ncol = NParam) }
HistParamT <- matrix(NA, nrow = 500*NParam, ncol = NParam)
HistCrit <- matrix(NA, nrow = 500*NParam, ncol = 1) ##_variables_initialisation
CritName <- NULL ParamFinalR <- NULL
CritBestValue <- NULL ParamFinalT <- NULL
Multiplier <- NULL CritFinal <- NULL
CritOptim <- +1E100 NRuns <- 0
##_temporary_change_of_Outputs_Sim NIter <- 0
RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal ### this reduces the size of the matrix exchange with fortran and therefore speeds the calibration if ("StartParamDistrib" %in% names(CalibOptions)) {
PrefilteringType <- 2
} else {
PrefilteringType <- 1
##_____Parameter_Grid_Screening____________________________________________________________ }
if (PrefilteringType == 1) {
NParam <- ncol(CalibOptions$StartParamList)
##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter }
ProposeCandidatesGrid <- function(DistribParam) { if (PrefilteringType == 2) {
NewCandidates <- expand.grid(lapply(seq_len(ncol(DistribParamR)), function(x) DistribParam[, x])) NParam <- ncol(CalibOptions$StartParamDistrib)
NewCandidates <- unique(NewCandidates) # to avoid duplicates when a parameter is set }
Output <- list(NewCandidates = NewCandidates) if (NParam > 20) {
} stop("Calibration_Michel can handle a maximum of 20 parameters \n")
return(NULL)
}
##Creation_of_new_candidates_______________________________________________ HistParamR <- matrix(NA, nrow = 500 * NParam, ncol = NParam)
OptimParam <- is.na(CalibOptions$FixedParam) 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) { if (PrefilteringType == 1) {
CandidatesParamR <- CalibOptions$StartParamList message("List-Screening in progress (", appendLF = FALSE)
} }
if (PrefilteringType == 2) { if (PrefilteringType == 2) {
DistribParamR <- CalibOptions$StartParamDistrib message("Grid-Screening in progress (", appendLF = FALSE)
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)
} }
message("0%", appendLF = FALSE)
##Loop_to_test_the_various_candidates______________________________________ }
iNewOptim <- 0 for (iNew in 1:nrow(CandidatesParamR)) {
Ncandidates <- nrow(CandidatesParamR)
if (verbose & Ncandidates > 1) { if (verbose & Ncandidates > 1) {
if (PrefilteringType == 1) { for (k in c(2, 4, 6, 8)) {
message("List-Screening in progress (", appendLF = FALSE) if (iNew == round(k / 10 * Ncandidates)) {
} message(" ", 10 * k, "%", 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 <- FUN_MOD(InputsModel, RunOptions, Param)
##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit, OutputsModel, verbose = FALSE)
if (!is.na(OutputsCrit$CritValue)) {
if (OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim) {
CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier
iNewOptim <- iNew
} }
} }
##Storage_of_crit_info }
if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) { ##Model_run
CritName <- OutputsCrit$CritName Param <- CandidatesParamR[iNew, ]
CritBestValue <- OutputsCrit$CritBestValue OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
Multiplier <- OutputsCrit$Multiplier ##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit, OutputsModel, verbose = FALSE)
if (!is.na(OutputsCrit$CritValue)) {
if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
iNewOptim <- iNew
} }
} }
if (verbose & Ncandidates > 1) { ##Storage_of_crit_info
message(" 100%)\n", appendLF = FALSE) if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
CritName <- OutputsCrit$CritName
CritBestValue <- OutputsCrit$CritBestValue
Multiplier <- OutputsCrit$Multiplier
} }
}
if (verbose & Ncandidates > 1) {
##End_of_first_step_Parameter_Screening____________________________________ message(" 100%)\n", appendLF = FALSE)
ParamStartR <- CandidatesParamR[iNewOptim, ] }
if (!is.matrix(ParamStartR)) {
ParamStartR <- matrix(ParamStartR, nrow = 1)
##End_of_first_step_Parameter_Screening____________________________________
ParamStartR <- CandidatesParamR[iNewOptim, ]
if (!is.matrix(ParamStartR)) {
ParamStartR <- matrix(ParamStartR, nrow = 1)
}
ParamStartT <- FUN_TRANSFO(ParamStartR, "RT")
CritStart <- CritOptim
NRuns <- NRuns+nrow(CandidatesParamR)
if (verbose) {
if (Ncandidates > 1) {
message(sprintf("\t Screening completed (%s runs)", NRuns))
} }
ParamStartT <- FUN_TRANSFO(ParamStartR, "RT") if (Ncandidates == 1) {
CritStart <- CritOptim message("\t Starting point for steepest-descent local search:")
NRuns <- NRuns+nrow(CandidatesParamR)
if (verbose) {
if (Ncandidates > 1) {
message(sprintf("\t Screening completed (%s runs)", NRuns))
}
if (Ncandidates == 1) {
message("\t Starting point for steepest-descent local search:")
}
message("\t Param = ", paste(sprintf("%8.3f", ParamStartR), collapse = " , "))
message(sprintf("\t Crit %-12s = %.4f", CritName, CritStart*Multiplier))
} }
##Results_archiving________________________________________________________ message("\t Param = ", paste(sprintf("%8.3f", ParamStartR), collapse = " , "))
HistParamR[1, ] <- ParamStartR message(sprintf("\t Crit %-12s = %.4f", CritName, CritStart * Multiplier))
HistParamT[1, ] <- ParamStartT }
HistCrit[1, ] <- CritStart ##Results_archiving________________________________________________________
HistParamR[1, ] <- ParamStartR
HistParamT[1, ] <- ParamStartT
HistCrit[1, ] <- CritStart
##_____Steepest_Descent_Local_Search_______________________________________________________
##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure ##_____Steepest_Descent_Local_Search_______________________________________________________
ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam,Pace) {
##Format_checking
if (nrow(NewParamOptimT) != 1 | nrow(OldParamOptimT) != 1) { ##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure
stop("each input set must be a matrix of one single line \n") ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam,Pace) {
return(NULL) ##Format_checking
} if (nrow(NewParamOptimT) != 1 | nrow(OldParamOptimT) != 1) {
if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) { stop("each input set must be a matrix of one single line \n")
stop("each input set must have the same number of values \n") return(NULL)
return(NULL) }
} if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
##Proposal_of_new_parameter_sets ###(local search providing 2*NParam-1 new sets) stop("each input set must have the same number of values \n")
NParam <- ncol(NewParamOptimT) return(NULL)
VECT <- NULL }
for (I in 1:NParam) { ##Proposal_of_new_parameter_sets ###(local search providing 2 * NParam-1 new sets)
##We_check_that_the_current_parameter_should_indeed_be_optimised NParam <- ncol(NewParamOptimT)
if (OptimParam[I] == TRUE) { VECT <- NULL
for (J in 1:2) { for (I in 1:NParam) {
Sign <- 2 * J - 3 #Sign can be equal to -1 or +1 ##We_check_that_the_current_parameter_should_indeed_be_optimised
##We_define_the_new_potential_candidate if (OptimParam[I] == TRUE) {
Add <- TRUE for (J in 1:2) {
PotentialCandidateT <- NewParamOptimT Sign <- 2 * J - 3 #Sign can be equal to -1 or +1
PotentialCandidateT[1, I] <- NewParamOptimT[I] + Sign * Pace ##We_define_the_new_potential_candidate
##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary Add <- TRUE
if (PotentialCandidateT[1, I] < RangesT[1, I] ) { PotentialCandidateT <- NewParamOptimT
PotentialCandidateT[1,I] <- RangesT[1, I] PotentialCandidateT[1, I] <- NewParamOptimT[I] + Sign * Pace
} ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
if (PotentialCandidateT[1, I] > RangesT[2, I]) { if (PotentialCandidateT[1, I] < RangesT[1, I] ) {
PotentialCandidateT[1,I] <- RangesT[2,I] PotentialCandidateT[1,I] <- RangesT[1, I]
} }
##We_check_the_set_is_not_outside_the_range_of_possible_values if (PotentialCandidateT[1, I] > RangesT[2, I]) {
if (NewParamOptimT[I] == RangesT[1, I] & Sign < 0) { PotentialCandidateT[1,I] <- RangesT[2,I]
Add <- FALSE }
} ##We_check_the_set_is_not_outside_the_range_of_possible_values
if (NewParamOptimT[I] == RangesT[2, I] & Sign > 0) { if (NewParamOptimT[I] == RangesT[1, I] & Sign < 0) {
Add <- FALSE Add <- FALSE
} }
##We_check_that_this_set_has_not_been_tested_during_the_last_iteration if (NewParamOptimT[I] == RangesT[2, I] & Sign > 0) {
if (identical(PotentialCandidateT, OldParamOptimT)) { Add <- FALSE
Add <- FALSE }
} ##We_check_that_this_set_has_not_been_tested_during_the_last_iteration
##We_add_the_candidate_to_our_list if (identical(PotentialCandidateT, OldParamOptimT)) {
if (Add == TRUE) { Add <- FALSE
VECT <- c(VECT, PotentialCandidateT) }
} ##We_add_the_candidate_to_our_list
if (Add) {
VECT <- c(VECT, PotentialCandidateT)
} }
} }
} }
Output <- NULL
Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
return(Output)
}
##Initialisation_of_variables
if (verbose) {
message("Steepest-descent local search in progress")
} }
Pace <- 0.64 Output <- NULL
PaceDiag <- rep(0, NParam) Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
CLG <- 0.7^(1/NParam) return(Output)
Compt <- 0 }
CritOptim <- CritStart
##Conversion_of_real_parameter_values
RangesR <- CalibOptions$SearchRanges ##Initialisation_of_variables
RangesT <- FUN_TRANSFO(RangesR, "RT") if (verbose) {
NewParamOptimT <- ParamStartT message("Steepest-descent local search in progress")
OldParamOptimT <- ParamStartT }
Pace <- 0.64
PaceDiag <- rep(0, NParam)
##START_LOOP_ITER_________________________________________________________ CLG <- 0.7^(1 / NParam)
for (ITER in 1:(100*NParam)) { Compt <- 0
CritOptim <- CritStart
##Conversion_of_real_parameter_values
RangesR <- CalibOptions$SearchRanges
RangesT <- FUN_TRANSFO(RangesR, "RT")
NewParamOptimT <- ParamStartT
OldParamOptimT <- ParamStartT
##START_LOOP_ITER_________________________________________________________
for (ITER in 1:(100 * NParam)) {
##Exit_loop_when_Pace_becomes_too_small___________________________________ ##Exit_loop_when_Pace_becomes_too_small___________________________________
if (Pace < 0.01) { if (Pace < 0.01) {
break break
} }