diff --git a/DESCRIPTION b/DESCRIPTION
index 4b22ecb57e4cc000e8476ab3b6b0f0f7c6bc3c9d..4bfc5252534c79d2cc076a4612225772bc3a43ec 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.0.13.8
-Date: 2018-08-30
+Version: 1.0.13.9
+Date: 2018-08-31
Authors@R: c(
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")),
diff --git a/NEWS.rmd b/NEWS.rmd
index 0c75d1dd44fb35b423e2c3830b6e469610be6c92..9d066e0a2872eaafa4cb71ccffb64fc59b62568e 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -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
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index 691d113401d5038da98f1ccb2ef6368bfc7d11c7..d1f0dcddf5422a7136b50b7513fa371c085e7bd2 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -1,315 +1,319 @@
-Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE) {
-
-
-##_____Arguments_check_____________________________________________________________________
- if (!inherits(InputsModel, "InputsModel")) {
- stop("InputsModel must be of class 'InputsModel' \n")
- return(NULL)
- }
- if (!inherits(RunOptions, "RunOptions")) {
- stop("RunOptions must be of class 'RunOptions' \n")
- return(NULL)
- }
- if (!inherits(InputsCrit, "InputsCrit")) {
- stop("InputsCrit must be of class 'InputsCrit' \n")
- return(NULL)
- }
- if (!inherits(CalibOptions, "CalibOptions")) {
- 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")
- return(NULL)
- }
-
-
- ##_check_FUN_TRANSFO
- if (is.null(FUN_TRANSFO)) {
- 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_GR5J )) {
- FUN_TRANSFO <- TransfoParam_GR5J
- }
- if (identical(FUN_MOD, RunModel_GR6J )) {
- FUN_TRANSFO <- TransfoParam_GR6J
- }
- if (identical(FUN_MOD, RunModel_GR2M )) {
- FUN_TRANSFO <- TransfoParam_GR2M
+Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions,
+ FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE) {
+
+
+ ##_____Arguments_check_____________________________________________________________________
+ if (!inherits(InputsModel, "InputsModel")) {
+ stop("InputsModel must be of class 'InputsModel' \n")
+ return(NULL)
+ }
+ if (!inherits(RunOptions, "RunOptions")) {
+ stop("RunOptions must be of class 'RunOptions' \n")
+ return(NULL)
+ }
+ if (!inherits(InputsCrit, "InputsCrit")) {
+ stop("InputsCrit must be of class 'InputsCrit' \n")
+ return(NULL)
+ }
+ if (!inherits(CalibOptions, "CalibOptions")) {
+ 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")
+ return(NULL)
+ }
+
+
+ ##_check_FUN_TRANSFO
+ if (is.null(FUN_TRANSFO)) {
+ 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_GR5J )) {
+ FUN_TRANSFO <- TransfoParam_GR5J
+ }
+ 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_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 )) {
- FUN_TRANSFO <- TransfoParam_GR1A
+ if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
+ FUN1 <- TransfoParam_GR5J
+ FUN2 <- TransfoParam_CemaNeige
}
- if (identical(FUN_MOD, RunModel_CemaNeige )) {
- FUN_TRANSFO <- TransfoParam_CemaNeige
+ if (identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
+ FUN1 <- TransfoParam_GR6J
+ FUN2 <- 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_CemaNeigeGR5J)) {
- FUN1 <- TransfoParam_GR5J
- FUN2 <- TransfoParam_CemaNeige
+ FUN_TRANSFO <- function(ParamIn, Direction) {
+ Bool <- is.matrix(ParamIn)
+ if (Bool == FALSE) {
+ ParamIn <- rbind(ParamIn)
}
- if (identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
- FUN1 <- TransfoParam_GR6J
- FUN2 <- TransfoParam_CemaNeige
- }
- FUN_TRANSFO <- function(ParamIn, Direction) {
- Bool <- is.matrix(ParamIn)
- 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)
+ 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)
}
- 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) {
- stop("Calibration_Michel can handle a maximum of 20 parameters \n")
+ if (is.null(FUN_TRANSFO)) {
+ stop("FUN_TRANSFO was not found (in Calibration function) \n")
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)
- 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)
+ }
+
+ ##_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) {
+ stop("Calibration_Michel can handle a maximum of 20 parameters \n")
+ 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)
+ 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) {
- CandidatesParamR <- CalibOptions$StartParamList
+ message("List-Screening in progress (", appendLF = FALSE)
}
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)
+ message("Grid-Screening in progress (", appendLF = FALSE)
}
-
- ##Loop_to_test_the_various_candidates______________________________________
- iNewOptim <- 0
- Ncandidates <- nrow(CandidatesParamR)
+ message("0%", appendLF = FALSE)
+ }
+ for (iNew in 1: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 <- 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
+ for (k in c(2, 4, 6, 8)) {
+ if (iNew == round(k / 10 * Ncandidates)) {
+ message(" ", 10 * k, "%", appendLF = FALSE)
}
- }
- ##Storage_of_crit_info
- if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
- CritName <- OutputsCrit$CritName
- CritBestValue <- OutputsCrit$CritBestValue
- Multiplier <- OutputsCrit$Multiplier
+ }
+ }
+ ##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
}
}
- if (verbose & Ncandidates > 1) {
- message(" 100%)\n", appendLF = FALSE)
+ ##Storage_of_crit_info
+ if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
+ CritName <- OutputsCrit$CritName
+ CritBestValue <- OutputsCrit$CritBestValue
+ Multiplier <- OutputsCrit$Multiplier
}
-
-
- ##End_of_first_step_Parameter_Screening____________________________________
- ParamStartR <- CandidatesParamR[iNewOptim, ]
- if (!is.matrix(ParamStartR)) {
- ParamStartR <- matrix(ParamStartR, nrow = 1)
+ }
+ if (verbose & Ncandidates > 1) {
+ message(" 100%)\n", appendLF = FALSE)
+ }
+
+
+ ##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")
- CritStart <- CritOptim
- 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))
+ if (Ncandidates == 1) {
+ message("\t Starting point for steepest-descent local search:")
}
- ##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
- ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam,Pace) {
- ##Format_checking
- if (nrow(NewParamOptimT) != 1 | nrow(OldParamOptimT) != 1) {
- stop("each input set must be a matrix of one single line \n")
- return(NULL)
- }
- if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
- stop("each input set must have the same number of values \n")
- return(NULL)
- }
- ##Proposal_of_new_parameter_sets ###(local search providing 2*NParam-1 new sets)
- NParam <- ncol(NewParamOptimT)
- VECT <- NULL
- for (I in 1:NParam) {
- ##We_check_that_the_current_parameter_should_indeed_be_optimised
- if (OptimParam[I] == TRUE) {
- for (J in 1:2) {
- Sign <- 2 * J - 3 #Sign can be equal to -1 or +1
- ##We_define_the_new_potential_candidate
- Add <- TRUE
- PotentialCandidateT <- NewParamOptimT
- 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[1, I] ) {
- PotentialCandidateT[1,I] <- RangesT[1, I]
- }
- if (PotentialCandidateT[1, I] > RangesT[2, I]) {
- PotentialCandidateT[1,I] <- RangesT[2,I]
- }
- ##We_check_the_set_is_not_outside_the_range_of_possible_values
- if (NewParamOptimT[I] == RangesT[1, I] & Sign < 0) {
- Add <- FALSE
- }
- if (NewParamOptimT[I] == RangesT[2, I] & Sign > 0) {
- Add <- FALSE
- }
- ##We_check_that_this_set_has_not_been_tested_during_the_last_iteration
- if (identical(PotentialCandidateT, OldParamOptimT)) {
- Add <- FALSE
- }
- ##We_add_the_candidate_to_our_list
- if (Add == TRUE) {
- VECT <- c(VECT, PotentialCandidateT)
- }
+ message("\t Param = ", paste(sprintf("%8.3f", ParamStartR), collapse = " , "))
+ message(sprintf("\t Crit %-12s = %.4f", CritName, CritStart * Multiplier))
+ }
+ ##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
+ ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam,Pace) {
+ ##Format_checking
+ if (nrow(NewParamOptimT) != 1 | nrow(OldParamOptimT) != 1) {
+ stop("each input set must be a matrix of one single line \n")
+ return(NULL)
+ }
+ if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
+ stop("each input set must have the same number of values \n")
+ return(NULL)
+ }
+ ##Proposal_of_new_parameter_sets ###(local search providing 2 * NParam-1 new sets)
+ NParam <- ncol(NewParamOptimT)
+ VECT <- NULL
+ for (I in 1:NParam) {
+ ##We_check_that_the_current_parameter_should_indeed_be_optimised
+ if (OptimParam[I] == TRUE) {
+ for (J in 1:2) {
+ Sign <- 2 * J - 3 #Sign can be equal to -1 or +1
+ ##We_define_the_new_potential_candidate
+ Add <- TRUE
+ PotentialCandidateT <- NewParamOptimT
+ 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[1, I] ) {
+ PotentialCandidateT[1,I] <- RangesT[1, I]
+ }
+ if (PotentialCandidateT[1, I] > RangesT[2, I]) {
+ PotentialCandidateT[1,I] <- RangesT[2,I]
+ }
+ ##We_check_the_set_is_not_outside_the_range_of_possible_values
+ if (NewParamOptimT[I] == RangesT[1, I] & Sign < 0) {
+ Add <- FALSE
+ }
+ if (NewParamOptimT[I] == RangesT[2, I] & Sign > 0) {
+ Add <- FALSE
+ }
+ ##We_check_that_this_set_has_not_been_tested_during_the_last_iteration
+ if (identical(PotentialCandidateT, OldParamOptimT)) {
+ Add <- FALSE
+ }
+ ##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
- PaceDiag <- rep(0, NParam)
- CLG <- 0.7^(1/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)) {
-
-
+ 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
+ PaceDiag <- rep(0, NParam)
+ CLG <- 0.7^(1 / 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___________________________________
if (Pace < 0.01) {
break
}
-
-
+
+
##Creation_of_new_candidates______________________________________________
CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace)$NewCandidatesT
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
##Remplacement_of_non_optimised_values_____________________________________
CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
- return(x)})
+ return(x)
+ })
if (NParam > 1) {
CandidatesParamR <- t(CandidatesParamR)
} else {
CandidatesParamR <- cbind(CandidatesParamR)
}
-
-
+
+
##Loop_to_test_the_various_candidates_____________________________________
iNewOptim <- 0
for (iNew in 1:nrow(CandidatesParamR)) {
@@ -319,23 +323,23 @@ Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions
##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
+ if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
+ CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
iNewOptim <- iNew
}
}
}
NRuns <- NRuns + nrow(CandidatesParamR)
-
-
+
+
##When_a_progress_has_been_achieved_______________________________________
if (iNewOptim != 0) {
##We_store_the_optimal_set
OldParamOptimT <- NewParamOptimT
NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
- Compt <- Compt+1
+ Compt <- Compt + 1
##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
- if (Compt > 2*NParam) {
+ if (Compt > 2 * NParam) {
Pace <- Pace * 2
Compt <- 0
}
@@ -347,14 +351,14 @@ Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions
}
}
} else {
- ##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
+ ##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
Pace <- Pace / 2
Compt <- 0
}
-
-
+
+
##Test_of_an_additional_candidate_using_diagonal_progress_________________
- if (ITER > 4*NParam) {
+ if (ITER > 4 * NParam) {
NRuns <- NRuns + 1
iNewOptim <- 0
iNew <- 1
@@ -379,8 +383,8 @@ Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions
OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit, OutputsModel, verbose = FALSE)
- if (OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim) {
- CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier
+ if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
+ CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
iNewOptim <- iNew
}
##When_a_progress_has_been_achieved
@@ -388,65 +392,66 @@ Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions
OldParamOptimT <- NewParamOptimT
NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
}
-
+
}
-
+
##Results_archiving_______________________________________________________
NewParamOptimR <- FUN_TRANSFO(NewParamOptimT, "TR")
HistParamR[ITER+1, ] <- NewParamOptimR
HistParamT[ITER+1, ] <- NewParamOptimT
HistCrit[ITER+1, ] <- CritOptim
### if (verbose) { cat(paste("\t Iter ",formatC(ITER,format="d",width=3), " Crit ",formatC(CritOptim,format="f",digits=4), " Pace ",formatC(Pace,format="f",digits=4), "\n",sep=""))}
-
-
-
- } ##END_LOOP_ITER_________________________________________________________
- ITER <- ITER-1
-
- ##Case_when_the_starting_parameter_set_remains_the_best_solution__________
- if (CritOptim == CritStart & verbose) {
- message("\t No progress achieved")
- }
- ##End_of_Steepest_Descent_Local_Search____________________________________
- ParamFinalR <- NewParamOptimR
- ParamFinalT <- NewParamOptimT
- CritFinal <- CritOptim
- NIter <- 1 + ITER
- if (verbose) {
- message(sprintf("\t Calibration completed (%s iterations, %s runs)", NIter, NRuns))
- message("\t Param = ", paste(sprintf("%8.3f", ParamFinalR), collapse = " , "))
- message(sprintf("\t Crit %-12s = %.4f", CritName, CritFinal*Multiplier))
- }
- ##Results_archiving_______________________________________________________
- HistParamR <- cbind(HistParamR[1:NIter, ])
- colnames(HistParamR) <- paste0("Param", 1:NParam)
- HistParamT <- cbind(HistParamT[1:NIter, ])
- colnames(HistParamT) <- paste0("Param", 1:NParam)
- HistCrit <- cbind(HistCrit[1:NIter, ])
- ###colnames(HistCrit) <- paste("HistCrit")
-
- BoolCrit_Actual <- InputsCrit$BoolCrit
- BoolCrit_Actual[OutputsCrit$Ind_notcomputed] <- FALSE
- MatBoolCrit <- cbind(InputsCrit$BoolCrit, BoolCrit_Actual)
- colnames(MatBoolCrit) <- c("BoolCrit_Requested", "BoolCrit_Actual")
-
-
-##_____Output______________________________________________________________________________
- OutputsCalib <- list(ParamFinalR = as.double(ParamFinalR), CritFinal = CritFinal*Multiplier,
- NIter = NIter, NRuns = NRuns,
- HistParamR = HistParamR, HistCrit = HistCrit*Multiplier, MatBoolCrit = MatBoolCrit,
- CritName = CritName, CritBestValue = CritBestValue)
- class(OutputsCalib) <- c("OutputsCalib", "HBAN")
- return(OutputsCalib)
-
-
-
+
+ } ##END_LOOP_ITER_________________________________________________________
+ ITER <- ITER - 1
+
+
+ ##Case_when_the_starting_parameter_set_remains_the_best_solution__________
+ if (CritOptim == CritStart & verbose) {
+ message("\t No progress achieved")
+ }
+
+ ##End_of_Steepest_Descent_Local_Search____________________________________
+ ParamFinalR <- NewParamOptimR
+ ParamFinalT <- NewParamOptimT
+ CritFinal <- CritOptim
+ NIter <- 1 + ITER
+ if (verbose) {
+ message(sprintf("\t Calibration completed (%s iterations, %s runs)", NIter, NRuns))
+ message("\t Param = ", paste(sprintf("%8.3f", ParamFinalR), collapse = " , "))
+ message(sprintf("\t Crit %-12s = %.4f", CritName, CritFinal * Multiplier))
+ }
+ ##Results_archiving_______________________________________________________
+ HistParamR <- cbind(HistParamR[1:NIter, ])
+ colnames(HistParamR) <- paste0("Param", 1:NParam)
+ HistParamT <- cbind(HistParamT[1:NIter, ])
+ colnames(HistParamT) <- paste0("Param", 1:NParam)
+ HistCrit <- cbind(HistCrit[1:NIter, ])
+ ###colnames(HistCrit) <- paste("HistCrit")
+
+ BoolCrit_Actual <- InputsCrit$BoolCrit
+ BoolCrit_Actual[OutputsCrit$Ind_notcomputed] <- FALSE
+ MatBoolCrit <- cbind(InputsCrit$BoolCrit, BoolCrit_Actual)
+ colnames(MatBoolCrit) <- c("BoolCrit_Requested", "BoolCrit_Actual")
+
+
+ ##_____Output______________________________________________________________________________
+ OutputsCalib <- list(ParamFinalR = as.double(ParamFinalR), CritFinal = CritFinal * Multiplier,
+ NIter = NIter, NRuns = NRuns,
+ HistParamR = HistParamR, HistCrit = HistCrit * Multiplier,
+ MatBoolCrit = MatBoolCrit,
+ CritName = CritName, CritBestValue = CritBestValue)
+ class(OutputsCalib) <- c("OutputsCalib", "HBAN")
+ return(OutputsCalib)
+
+
+
}
-
+