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_CemaNeigeGR5J)) { FUN1 <- TransfoParam_GR5J FUN2 <- TransfoParam_CemaNeige } 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) } } 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") 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) { ##Managing_matrix_sizes Nvalmax <- nrow(DistribParam) NParam <- ncol(DistribParam) ##we_add_columns_to_MatDistrib_until_it_has_20_columns DistribParam2 <- matrix(NA, nrow = Nvalmax, ncol = 20) DistribParam2[1:Nvalmax, 1:NParam] <- DistribParam ##we_check_the_number_of_values_to_test_for_each_param NbDistrib <- rep(1, 20) for (iC in 1:20) { NbDistrib[iC] <- max(1, Nvalmax-sum(is.na(DistribParam2[, iC]))) } ##Loop_on_the_various_values_to_test ###(if 4 param and 3 values for each => 3^4 sets) ##NB_we_always_do_20_loops ###which_is_here_the_max_number_of_param_that_can_be_optimised VECT <- NULL for (iL01 in 1:NbDistrib[01]) { for (iL02 in 1:NbDistrib[02]) { for (iL03 in 1:NbDistrib[03]) { for (iL04 in 1:NbDistrib[04]) { for (iL05 in 1:NbDistrib[05]) { for (iL06 in 1:NbDistrib[06]) { for (iL07 in 1:NbDistrib[07]) { for (iL08 in 1:NbDistrib[08]) { for (iL09 in 1:NbDistrib[09]) { for (iL10 in 1:NbDistrib[10]) { for (iL11 in 1:NbDistrib[11]) { for (iL12 in 1:NbDistrib[12]) { for (iL13 in 1:NbDistrib[13]) { for (iL14 in 1:NbDistrib[14]) { for (iL15 in 1:NbDistrib[15]) { for (iL16 in 1:NbDistrib[16]) { for (iL17 in 1:NbDistrib[17]) { for (iL18 in 1:NbDistrib[18]) { for (iL19 in 1:NbDistrib[19]) { for (iL20 in 1:NbDistrib[20]) { VECT <- c(VECT, DistribParam2[iL01,01],DistribParam2[iL02,02],DistribParam2[iL03,03],DistribParam2[iL04,04],DistribParam2[iL05,05], DistribParam2[iL06,06],DistribParam2[iL07,07],DistribParam2[iL08,08],DistribParam2[iL09,09],DistribParam2[iL10,10], DistribParam2[iL11,11],DistribParam2[iL12,12],DistribParam2[iL13,13],DistribParam2[iL14,14],DistribParam2[iL15,15], DistribParam2[iL16,16],DistribParam2[iL17,17],DistribParam2[iL18,18],DistribParam2[iL19,19],DistribParam2[iL20,20]) } } } } } } } } } } } } } } } } } } } } MAT <- matrix(VECT, ncol = 20, byrow = TRUE)[, 1:NParam] if (!is.matrix(MAT)) { MAT <- cbind(MAT) } Output <- NULL Output$NewCandidates <- MAT return(Output) } ##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 <- 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)) { CritName <- OutputsCrit$CritName CritBestValue <- OutputsCrit$CritBestValue Multiplier <- OutputsCrit$Multiplier } } 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)) } 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________________________________________________________ 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) } } } } 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)}) if (NParam > 1) { CandidatesParamR <- t(CandidatesParamR) } else { CandidatesParamR <- cbind(CandidatesParamR) } ##Loop_to_test_the_various_candidates_____________________________________ iNewOptim <- 0 for (iNew in 1:nrow(CandidatesParamR)) { ##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 } } } 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 ##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row if (Compt > 2*NParam) { Pace <- Pace * 2 Compt <- 0 } ##We_update_PaceDiag VectPace <- NewParamOptimT-OldParamOptimT for (iC in 1:NParam) { if (OptimParam[iC]) { if (VectPace[iC] != 0) { PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC] } if (VectPace[iC] == 0) { PaceDiag[iC] <- CLG*PaceDiag[iC] } } } } else { ##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) { NRuns <- NRuns + 1 iNewOptim <- 0 iNew <- 1 CandidatesParamT <- NewParamOptimT+PaceDiag if (!is.matrix(CandidatesParamT)) { CandidatesParamT <- matrix(CandidatesParamT, nrow = 1) } ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary for (iC in 1:NParam) { if (OptimParam[iC]) { if (CandidatesParamT[iNew, iC] < RangesT[1, iC]) { CandidatesParamT[iNew, iC] <- RangesT[1, iC] } if (CandidatesParamT[iNew, iC] > RangesT[2, iC]) { CandidatesParamT[iNew, iC] <- RangesT[2, iC] } } } CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR") ##Model_run Param <- CandidatesParamR[iNew, ] 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 iNewOptim <- iNew } ##When_a_progress_has_been_achieved if (iNewOptim != 0) { 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) }