Skip to content
GitLab
Commit 336d5e9b authored by unknown's avatar unknown
Browse files

v1.0.4.1 function Calibration_Michel() cleaned

parent 4b4c4aee
dev 115-add-pcict-date-format-managment 118-plot-outputsmodel-remove-case-sensitivity-to-the-which-argument 136-decreased-performance-of-calibration-execution-time 60-migrate-transfoparam-functions-and-createcaliboptions-to-s3-methods 96-parallelize-the-grid-screening-step-during-calibration airGRplus cleanAirGRplus master CRAN_v1.7.0 CRAN_v1.6.12 CRAN_v1.6.11 CRAN_v1.6.10.4 CRAN_v1.6.9.27 CRAN_v1.6.9.24 CRAN_v1.6.9.23 CRAN_v1.6.9.21 CRAN_v1.4.3.65 CRAN_v1.4.3.60 CRAN_v1.4.3.52 CRAN_v1.3.2.42 CRAN_v1.3.2.41 CRAN_v1.3.2.40 CRAN_v1.3.2.39 CRAN_v1.3.2.23 CRAN_v1.3.2.17 CRAN_v1.2.13.16 CRAN_v1.2.13.13 CRAN_v1.0.14.1 CRAN_v1.0.13.19 CRAN_v1.0.10.11 CRAN_v1.0.9.64 CRAN_v1.0.5.12 CRAN_v1.0.5.11 CRAN_v1.0.5.10
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 346 additions and 208 deletions
+346 -208
DESCRIPTION
+ 2
− 2
View file @ 336d5e9b
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.3 Version: 1.0.4.1
Date: 2016-12-09 Date: 2017-01-18
Authors@R: c( Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl")), person("Laurent", "Coron", role = c("aut", "trl")),
person("Charles", "Perrin", role = c("aut", "ths")), person("Charles", "Perrin", role = c("aut", "ths")),
......
R/Calibration_Michel.R
+ 340
− 202
View file @ 336d5e9b
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_____________________________________________________________________ ##_____Arguments_check_____________________________________________________________________
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); } if (!inherits(InputsModel, "InputsModel")) {
if(inherits(RunOptions,"RunOptions")==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); } stop("InputsModel must be of class 'InputsModel' \n")
if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); } return(NULL)
if(inherits(CalibOptions,"CalibOptions")==FALSE){ stop("CalibOptions must be of class 'CalibOptions' \n"); return(NULL); } }
if(inherits(CalibOptions,"HBAN")==FALSE){ stop("CalibOptions must be of class 'HBAN' if Calibration_Michel is used \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 ##_check_FUN_TRANSFO
if(is.null(FUN_TRANSFO)){ if (is.null(FUN_TRANSFO)) {
if(identical(FUN_MOD,RunModel_GR4H )){ FUN_TRANSFO <- TransfoParam_GR4H ; } if (identical(FUN_MOD, RunModel_GR4H )) {
if(identical(FUN_MOD,RunModel_GR4J )){ FUN_TRANSFO <- TransfoParam_GR4J ; } FUN_TRANSFO <- TransfoParam_GR4H
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_GR4J )) {
if(identical(FUN_MOD,RunModel_GR2M )){ FUN_TRANSFO <- TransfoParam_GR2M ; } FUN_TRANSFO <- TransfoParam_GR4J
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_GR5J )) {
if(identical(FUN_MOD,RunModel_CemaNeigeGR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){ FUN_TRANSFO <- TransfoParam_GR5J
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_GR6J )) {
if(identical(FUN_MOD,RunModel_CemaNeigeGR6J)){ FUN1 <- TransfoParam_GR6J; FUN2 <- TransfoParam_CemaNeige; } FUN_TRANSFO <- TransfoParam_GR6J
FUN_TRANSFO <- function(ParamIn,Direction){ }
Bool <- is.matrix(ParamIn); if (identical(FUN_MOD, RunModel_GR2M )) {
if(Bool==FALSE){ ParamIn <- rbind(ParamIn); } FUN_TRANSFO <- TransfoParam_GR2M
ParamOut <- NA*ParamIn; }
NParam <- ncol(ParamIn); if (identical(FUN_MOD, RunModel_GR1A )) {
ParamOut[, 1:(NParam-2)] <- FUN1(ParamIn[, 1:(NParam-2)],Direction); FUN_TRANSFO <- TransfoParam_GR1A
ParamOut[,(NParam-1):NParam ] <- FUN2(ParamIn[,(NParam-1):NParam ],Direction); }
if(Bool==FALSE){ ParamOut <- ParamOut[1,]; } if (identical(FUN_MOD, RunModel_CemaNeige )) {
return(ParamOut); 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)
} }
if(is.null(FUN_TRANSFO)){ stop("FUN_TRANSFO was not found (in Calibration function) \n"); return(NULL); }
} }
##_variables_initialisation ##_variables_initialisation
ParamFinalR <- NULL; ParamFinalT <- NULL; CritFinal <- NULL; ParamFinalR <- NULL
NRuns <- 0; NIter <- 0; ParamFinalT <- NULL
if("StartParamDistrib" %in% names(CalibOptions)){ PrefilteringType <- 2; } else { PrefilteringType <- 1; } CritFinal <- NULL
if(PrefilteringType==1){ NParam <- ncol(CalibOptions$StartParamList); } NRuns <- 0
if(PrefilteringType==2){ NParam <- ncol(CalibOptions$StartParamDistrib); } NIter <- 0
if(NParam>20){ stop("Calibration_Michel can handle a maximum of 20 parameters \n"); return(NULL); } if ("StartParamDistrib" %in% names(CalibOptions)) {
HistParamR <- matrix(NA,nrow=500*NParam,ncol=NParam); PrefilteringType <- 2
HistParamT <- matrix(NA,nrow=500*NParam,ncol=NParam); } else {
HistCrit <- matrix(NA,nrow=500*NParam,ncol=1); PrefilteringType <- 1
CritName <- NULL; }
CritBestValue <- NULL; if (PrefilteringType == 1) {
Multiplier <- NULL; NParam <- ncol(CalibOptions$StartParamList)
CritOptim <- +1E100; }
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 ##_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 RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal ### this reduces the size of the matrix exchange with fortran and therefore speeds the calibration
...@@ -59,98 +118,128 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU ...@@ -59,98 +118,128 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU
##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter ##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter
ProposeCandidatesGrid <- function(DistribParam){ ProposeCandidatesGrid <- function(DistribParam) {
##Managing_matrix_sizes ##Managing_matrix_sizes
Nvalmax <- nrow(DistribParam); Nvalmax <- nrow(DistribParam)
NParam <- ncol(DistribParam); NParam <- ncol(DistribParam)
##we_add_columns_to_MatDistrib_until_it_has_20_columns ##we_add_columns_to_MatDistrib_until_it_has_20_columns
DistribParam2 <- matrix(NA,nrow=Nvalmax,ncol=20); DistribParam2 <- matrix(NA, nrow = Nvalmax, ncol = 20)
DistribParam2[1:Nvalmax,1:NParam] <- DistribParam; DistribParam2[1:Nvalmax, 1:NParam] <- DistribParam
##we_check_the_number_of_values_to_test_for_each_param ##we_check_the_number_of_values_to_test_for_each_param
NbDistrib <- rep(1,20); NbDistrib <- rep(1, 20)
for(iC in 1:20){ NbDistrib[iC] <- max( 1 , Nvalmax-sum(is.na(DistribParam2[,iC])) ); } 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) ##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 ##NB_we_always_do_20_loops ###which_is_here_the_max_number_of_param_that_can_be_optimised
VECT <- NULL; 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 (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 (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 (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]){ 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, VECT <- c(VECT,
DistribParam2[iL01,01],DistribParam2[iL02,02],DistribParam2[iL03,03],DistribParam2[iL04,04],DistribParam2[iL05,05], 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[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[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]); DistribParam2[iL16,16],DistribParam2[iL17,17],DistribParam2[iL18,18],DistribParam2[iL19,19],DistribParam2[iL20,20])
} } } } } } } } } }
} } } } } } } } } }
} } } } } } } } } }
} } } } } } } } } }
MAT <- matrix(VECT,ncol=20,byrow=TRUE)[,1:NParam]; MAT <- matrix(VECT, ncol = 20, byrow = TRUE)[, 1:NParam]
if(is.matrix(MAT)==FALSE){ MAT <- cbind(MAT); } if (!is.matrix(MAT)) {
Output <- NULL; MAT <- cbind(MAT)
Output$NewCandidates <- MAT; }
return(Output); Output <- NULL
Output$NewCandidates <- MAT
return(Output)
} }
##Creation_of_new_candidates_______________________________________________ ##Creation_of_new_candidates_______________________________________________
OptimParam <- is.na(CalibOptions$FixedParam) OptimParam <- is.na(CalibOptions$FixedParam)
if(PrefilteringType==1){ CandidatesParamR <- CalibOptions$StartParamList; } if (PrefilteringType == 1) {
if(PrefilteringType==2){ DistribParamR <- CalibOptions$StartParamDistrib; DistribParamR[,!OptimParam] <- NA; CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates; } CandidatesParamR <- CalibOptions$StartParamList
}
if (PrefilteringType == 2) {
DistribParamR <- CalibOptions$StartParamDistrib
DistribParamR[,!OptimParam] <- NA
CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates
}
##Remplacement_of_non_optimised_values_____________________________________ ##Remplacement_of_non_optimised_values_____________________________________
CandidatesParamR <- apply(CandidatesParamR,1,function(x){ x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]; return(x); }); CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
if(NParam>1){ CandidatesParamR <- t(CandidatesParamR); } else { CandidatesParamR <- cbind(CandidatesParamR); } x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
return(x)})
if (NParam>1) {
CandidatesParamR <- t(CandidatesParamR)
} else { CandidatesParamR <- cbind(CandidatesParamR)
}
##Loop_to_test_the_various_candidates______________________________________ ##Loop_to_test_the_various_candidates______________________________________
iNewOptim <- 0; iNewOptim <- 0
Ncandidates <- nrow(CandidatesParamR); Ncandidates <- nrow(CandidatesParamR)
if(verbose & Ncandidates>1){ if (verbose & Ncandidates > 1) {
if(PrefilteringType==1){ message("List-Screening in progress (", appendLF = FALSE) } if (PrefilteringType == 1) {
if(PrefilteringType==2){ message("Grid-Screening in progress (", appendLF = FALSE) } message("List-Screening in progress (", appendLF = FALSE)
}
if (PrefilteringType == 2) {
message("Grid-Screening in progress (", appendLF = FALSE)
}
message("0%", appendLF = FALSE) message("0%", appendLF = FALSE)
} }
for(iNew in 1:nrow(CandidatesParamR)){ for (iNew in 1:nrow(CandidatesParamR)) {
if(verbose & Ncandidates>1){ if (verbose & Ncandidates > 1) {
for(k in c(2,4,6,8)){ if(iNew==round(k/10*Ncandidates)){ message(" ", 10*k, "%",appendLF = FALSE) } } for (k in c(2, 4, 6, 8)) {
if (iNew == round(k/10*Ncandidates)) {
message(" ", 10*k, "%", appendLF = FALSE)
}
}
} }
##Model_run ##Model_run
Param <- CandidatesParamR[iNew,]; Param <- CandidatesParamR[iNew, ]
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param); OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
##Calibration_criterion_computation ##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel, verbose = FALSE); OutputsCrit <- FUN_CRIT(InputsCrit, OutputsModel, verbose = FALSE)
if(!is.na(OutputsCrit$CritValue)){ if(OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim){ if (!is.na(OutputsCrit$CritValue)) {
CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier; if (OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim) {
iNewOptim <- iNew; CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier
} } iNewOptim <- iNew
}
}
##Storage_of_crit_info ##Storage_of_crit_info
if(is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)){ if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
CritName <- OutputsCrit$CritName; CritName <- OutputsCrit$CritName
CritBestValue <- OutputsCrit$CritBestValue; CritBestValue <- OutputsCrit$CritBestValue
Multiplier <- OutputsCrit$Multiplier; Multiplier <- OutputsCrit$Multiplier
} }
} }
if(verbose & Ncandidates>1){ message(" 100%)\n", appendLF = FALSE) } if (verbose & Ncandidates > 1) {
message(" 100%)\n", appendLF = FALSE)
}
##End_of_first_step_Parameter_Screening____________________________________ ##End_of_first_step_Parameter_Screening____________________________________
ParamStartR <- CandidatesParamR[iNewOptim,]; if(!is.matrix(ParamStartR)){ ParamStartR <- matrix(ParamStartR,nrow=1); } ParamStartR <- CandidatesParamR[iNewOptim, ]
ParamStartT <- FUN_TRANSFO(ParamStartR,"RT"); if (!is.matrix(ParamStartR)) {
CritStart <- CritOptim; ParamStartR <- matrix(ParamStartR, nrow = 1)
NRuns <- NRuns+nrow(CandidatesParamR); }
if(verbose){ ParamStartT <- FUN_TRANSFO(ParamStartR, "RT")
if(Ncandidates> 1){ CritStart <- CritOptim
NRuns <- NRuns+nrow(CandidatesParamR)
if (verbose) {
if (Ncandidates > 1) {
message("\t Screening completed (", NRuns, " runs):") message("\t Screening completed (", NRuns, " runs):")
} }
if(Ncandidates==1){ if (Ncandidates == 1) {
message("\t Starting point for steepest-descent local search:") message("\t Starting point for steepest-descent local search:")
} }
message("\t Param = ", paste(formatC(ParamStartR, format = "f", width = 8, digits = 3), collapse = " , ")) message("\t Param = ", paste(formatC(ParamStartR, format = "f", width = 8, digits = 3), collapse = " , "))
message("\t Crit ", format(CritName, width = 12, justify = "left"), " = ", formatC(CritStart*Multiplier, format = "f", digits = 4)) message("\t Crit ", format(CritName, width = 12, justify = "left"), " = ", formatC(CritStart*Multiplier, format = "f", digits = 4))
} }
##Results_archiving________________________________________________________ ##Results_archiving________________________________________________________
HistParamR[1,] <- ParamStartR; HistParamR[1, ] <- ParamStartR
HistParamT[1,] <- ParamStartT; HistParamT[1, ] <- ParamStartT
HistCrit[1,] <- CritStart; HistCrit[1, ] <- CritStart
...@@ -159,186 +248,235 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU ...@@ -159,186 +248,235 @@ Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FU
##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure ##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure
ProposeCandidatesLoc <- function(NewParamOptimT,OldParamOptimT,RangesT,OptimParam,Pace){ ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam,Pace) {
##Format_checking ##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 (nrow(NewParamOptimT) != 1 | nrow(OldParamOptimT) != 1) {
if(ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT)!=length(OptimParam)){ stop("each input set must have the same number of values \n"); return(NULL); } 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) ##Proposal_of_new_parameter_sets ###(local search providing 2*NParam-1 new sets)
NParam <- ncol(NewParamOptimT); NParam <- ncol(NewParamOptimT)
VECT <- NULL; VECT <- NULL
for(I in 1:NParam){ for (I in 1:NParam) {
##We_check_that_the_current_parameter_should_indeed_be_optimised ##We_check_that_the_current_parameter_should_indeed_be_optimised
if(OptimParam[I]==TRUE){ if (OptimParam[I] == TRUE) {
for(J in 1:2){ for (J in 1:2) {
Sign <- 2*J-3; #Sign can be equal to -1 or +1 Sign <- 2 * J - 3 #Sign can be equal to -1 or +1
##We_define_the_new_potential_candidate ##We_define_the_new_potential_candidate
Add <- TRUE; Add <- TRUE
PotentialCandidateT <- NewParamOptimT; PotentialCandidateT <- NewParamOptimT
PotentialCandidateT[1,I] <- NewParamOptimT[I]+Sign*Pace; PotentialCandidateT[1, I] <- NewParamOptimT[I] + Sign * Pace
##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary ##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[1, I] ) {
if(PotentialCandidateT[1,I]>RangesT[2,I]){ PotentialCandidateT[1,I] <- RangesT[2,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 ##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[1, I] & Sign < 0) {
if( NewParamOptimT[I]==RangesT[2,I] & Sign>0 ){ Add <- FALSE; } 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 ##We_check_that_this_set_has_not_been_tested_during_the_last_iteration
if(identical(PotentialCandidateT,OldParamOptimT)){ Add <- FALSE; } if (identical(PotentialCandidateT, OldParamOptimT)) {
Add <- FALSE
}
##We_add_the_candidate_to_our_list ##We_add_the_candidate_to_our_list
if(Add==TRUE){ VECT <- c(VECT,PotentialCandidateT); } if (Add == TRUE) {
VECT <- c(VECT, PotentialCandidateT)
}
} }
} }
} }
Output <- NULL; Output <- NULL
Output$NewCandidatesT <- matrix(VECT,ncol=NParam,byrow=TRUE); Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
return(Output); return(Output)
} }
##Initialisation_of_variables ##Initialisation_of_variables
if(verbose){ if (verbose) {
message("Steepest-descent local search in progress") message("Steepest-descent local search in progress")
} }
Pace <- 0.64; Pace <- 0.64
PaceDiag <- rep(0,NParam); PaceDiag <- rep(0, NParam)
CLG <- 0.7^(1/NParam); CLG <- 0.7^(1/NParam)
Compt <- 0; Compt <- 0
CritOptim <- CritStart; CritOptim <- CritStart
##Conversion_of_real_parameter_values ##Conversion_of_real_parameter_values
RangesR <- CalibOptions$SearchRanges; RangesR <- CalibOptions$SearchRanges
RangesT <- FUN_TRANSFO(RangesR,"RT"); RangesT <- FUN_TRANSFO(RangesR, "RT")
NewParamOptimT <- ParamStartT; NewParamOptimT <- ParamStartT
OldParamOptimT <- ParamStartT; OldParamOptimT <- ParamStartT
##START_LOOP_ITER_________________________________________________________ ##START_LOOP_ITER_________________________________________________________
for(ITER in 1:(100*NParam)){ for (ITER in 1:(100*NParam)) {
##Exit_loop_when_Pace_becomes_too_small___________________________________ ##Exit_loop_when_Pace_becomes_too_small___________________________________
if(Pace<0.01){ break; } if (Pace < 0.01) {
break
}
##Creation_of_new_candidates______________________________________________ ##Creation_of_new_candidates______________________________________________
CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT,OldParamOptimT,RangesT,OptimParam,Pace)$NewCandidatesT; CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace)$NewCandidatesT
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR"); CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
##Remplacement_of_non_optimised_values_____________________________________ ##Remplacement_of_non_optimised_values_____________________________________
CandidatesParamR <- apply(CandidatesParamR,1,function(x){ x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]; return(x); }); CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
if(NParam>1){ CandidatesParamR <- t(CandidatesParamR); } else { CandidatesParamR <- cbind(CandidatesParamR); } x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
return(x)})
if (NParam > 1) {
CandidatesParamR <- t(CandidatesParamR)
} else {
CandidatesParamR <- cbind(CandidatesParamR)
}
##Loop_to_test_the_various_candidates_____________________________________ ##Loop_to_test_the_various_candidates_____________________________________
iNewOptim <- 0; iNewOptim <- 0
for(iNew in 1:nrow(CandidatesParamR)){ for (iNew in 1:nrow(CandidatesParamR)) {
##Model_run ##Model_run
Param <- CandidatesParamR[iNew,]; Param <- CandidatesParamR[iNew, ]
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param); OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
##Calibration_criterion_computation ##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel, verbose = FALSE); OutputsCrit <- FUN_CRIT(InputsCrit, OutputsModel, verbose = FALSE)
if(!is.na(OutputsCrit$CritValue)){ if(OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim){ if (!is.na(OutputsCrit$CritValue)) {
CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier; if (OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim) {
iNewOptim <- iNew; CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier
} } iNewOptim <- iNew
}
}
} }
NRuns <- NRuns+nrow(CandidatesParamR); NRuns <- NRuns + nrow(CandidatesParamR)
##When_a_progress_has_been_achieved_______________________________________ ##When_a_progress_has_been_achieved_______________________________________
if(iNewOptim!=0){ if (iNewOptim != 0) {
##We_store_the_optimal_set ##We_store_the_optimal_set
OldParamOptimT <- NewParamOptimT; OldParamOptimT <- NewParamOptimT
NewParamOptimT <- matrix(CandidatesParamT[iNewOptim,1:NParam],nrow=1); 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 ##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
if(Compt>2*NParam){ if (Compt > 2*NParam) {
Pace <- Pace*2; Pace <- Pace * 2
Compt <- 0; Compt <- 0
} }
##We_update_PaceDiag ##We_update_PaceDiag
VectPace <- NewParamOptimT-OldParamOptimT; VectPace <- NewParamOptimT-OldParamOptimT
for(iC in 1:NParam){ if(OptimParam[iC]){ for (iC in 1:NParam) {
if(VectPace[iC]!=0){ PaceDiag[iC] <- CLG*PaceDiag[iC]+(1-CLG)*VectPace[iC]; } if (OptimParam[iC]) {
if(VectPace[iC]==0){ PaceDiag[iC] <- CLG*PaceDiag[iC]; } if (VectPace[iC] != 0) {
} } PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC]
}
if (VectPace[iC] == 0) {
PaceDiag[iC] <- CLG*PaceDiag[iC]
}
}
}
} else { } else {
##When_no_progress_has_been_achieved_we_decrease_the_pace_________________ ##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
Pace <- Pace/2; Pace <- Pace / 2
Compt <- 0; Compt <- 0
} }
##Test_of_an_additional_candidate_using_diagonal_progress_________________ ##Test_of_an_additional_candidate_using_diagonal_progress_________________
if(ITER>4*NParam){ if (ITER > 4*NParam) {
NRuns <- NRuns + 1
NRuns <- NRuns+1; iNewOptim <- 0
iNewOptim <- 0; iNew <- 1; iNew <- 1
CandidatesParamT <- NewParamOptimT+PaceDiag; if(!is.matrix(CandidatesParamT)){ CandidatesParamT <- matrix(CandidatesParamT,nrow=1); } CandidatesParamT <- NewParamOptimT+PaceDiag
##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary if (!is.matrix(CandidatesParamT)) {
for(iC in 1:NParam){ if(OptimParam[iC]){ CandidatesParamT <- matrix(CandidatesParamT, nrow = 1)
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]; } ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
} } for (iC in 1:NParam) {
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR"); if (OptimParam[iC]) {
##Model_run if (CandidatesParamT[iNew, iC] < RangesT[1, iC]) {
Param <- CandidatesParamR[iNew,]; CandidatesParamT[iNew, iC] <- RangesT[1, iC]
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param); }
##Calibration_criterion_computation if (CandidatesParamT[iNew, iC] > RangesT[2, iC]) {
OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel, verbose = FALSE); CandidatesParamT[iNew, iC] <- RangesT[2, iC]
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);
} }
}
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_______________________________________________________ ##Results_archiving_______________________________________________________
NewParamOptimR <- FUN_TRANSFO(NewParamOptimT,"TR"); NewParamOptimR <- FUN_TRANSFO(NewParamOptimT, "TR")
HistParamR[ITER+1,] <- NewParamOptimR; HistParamR[ITER+1, ] <- NewParamOptimR
HistParamT[ITER+1,] <- NewParamOptimT; HistParamT[ITER+1, ] <- NewParamOptimT
HistCrit[ITER+1,] <- CritOptim; 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="")); } ### 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_________________________________________________________ } ##END_LOOP_ITER_________________________________________________________
ITER <- ITER-1; ITER <- ITER-1
##Case_when_the_starting_parameter_set_remains_the_best_solution__________ ##Case_when_the_starting_parameter_set_remains_the_best_solution__________
if(CritOptim==CritStart & verbose){ if (CritOptim == CritStart & verbose) {
message("\t No progress achieved"); message("\t No progress achieved")
} }
##End_of_Steepest_Descent_Local_Search____________________________________ ##End_of_Steepest_Descent_Local_Search____________________________________
ParamFinalR <- NewParamOptimR; ParamFinalR <- NewParamOptimR
ParamFinalT <- NewParamOptimT; ParamFinalT <- NewParamOptimT
CritFinal <- CritOptim; CritFinal <- CritOptim
NIter <- 1+ITER; NIter <- 1 + ITER
if(verbose){ if (verbose) {
message("\t Calibration completed (", NIter, " iterations, ", NRuns, " runs)") message("\t Calibration completed (", NIter, " iterations, ", NRuns, " runs)")
message("\t Param = ", paste(formatC(ParamFinalR, format = "f", width = 8, digits = 3), collapse = " , ")) message("\t Param = ", paste(formatC(ParamFinalR, format = "f", width = 8, digits = 3), collapse = " , "))
message("\t Crit ", format(CritName, width = 12, justify = "left"), " = ", formatC(CritFinal*Multiplier, format = "f", digits = 4)) message("\t Crit ", format(CritName, width = 12, justify = "left"), " = ", formatC(CritFinal*Multiplier, format = "f", digits = 4))
} }
##Results_archiving_______________________________________________________ ##Results_archiving_______________________________________________________
HistParamR <- cbind(HistParamR[1:NIter,]); colnames(HistParamR) <- paste("Param",1:NParam,sep=""); HistParamR <- cbind(HistParamR[1:NIter, ])
HistParamT <- cbind(HistParamT[1:NIter,]); colnames(HistParamT) <- paste("Param",1:NParam,sep=""); colnames(HistParamR) <- paste0("Param", 1:NParam)
HistCrit <- cbind(HistCrit[1:NIter,]); ###colnames(HistCrit) <- paste("HistCrit"); 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; BoolCrit_Actual <- InputsCrit$BoolCrit
MatBoolCrit <- cbind( InputsCrit$BoolCrit , BoolCrit_Actual ); BoolCrit_Actual[OutputsCrit$Ind_notcomputed] <- FALSE
colnames(MatBoolCrit) <- c("BoolCrit_Requested","BoolCrit_Actual"); MatBoolCrit <- cbind(InputsCrit$BoolCrit, BoolCrit_Actual)
colnames(MatBoolCrit) <- c("BoolCrit_Requested", "BoolCrit_Actual")
##_____Output______________________________________________________________________________ ##_____Output______________________________________________________________________________
OutputsCalib <- list(as.double(ParamFinalR),CritFinal*Multiplier,NIter,NRuns,HistParamR,HistCrit*Multiplier,MatBoolCrit,CritName,CritBestValue); OutputsCalib <- list(ParamFinalR = as.double(ParamFinalR), CritFinal = CritFinal*Multiplier,
names(OutputsCalib) <- c("ParamFinalR","CritFinal","NIter","NRuns","HistParamR","HistCrit","MatBoolCrit","CritName","CritBestValue"); NIter = NIter, NRuns = NRuns,
class(OutputsCalib) <- c("OutputsCalib","HBAN"); HistParamR = HistParamR, HistCrit = HistCrit*Multiplier, MatBoolCrit = MatBoolCrit,
return(OutputsCalib); CritName = CritName, CritBestValue = CritBestValue)
class(OutputsCalib) <- c("OutputsCalib", "HBAN")
return(OutputsCalib)
......
inst/CITATION
+ 4
− 4
View file @ 336d5e9b
...@@ -17,14 +17,14 @@ citEntry(entry="Manual", ...@@ -17,14 +17,14 @@ citEntry(entry="Manual",
title = "airGR: Suite of GR hydrological models for precipitation-runoff modelling", title = "airGR: Suite of GR hydrological models for precipitation-runoff modelling",
author = personList(as.person("L. Coron"), as.person("C. Perrin"), as.person("C. Michel")), author = personList(as.person("L. Coron"), as.person("C. Perrin"), as.person("C. Michel")),
journal = "R News", journal = "R News",
year = "2016", year = "2017",
note = "R package version 1.0.3", note = "R package version 1.0.4",
url = "https://webgr.irstea.fr/airGR/?lang=en", url = "https://webgr.irstea.fr/airGR/?lang=en",
textVersion = textVersion =
paste("Coron, L., Perrin, C. and Michel, C.", paste("Coron, L., Perrin, C. and Michel, C.",
"(2016).", "(2017).",
"airGR: Suite of GR hydrological models for precipitation-runoff modelling.", "airGR: Suite of GR hydrological models for precipitation-runoff modelling.",
"R package version 1.0.3.", "R package version 1.0.4.",
"https://webgr.irstea.fr/airGR/?lang=en.", "https://webgr.irstea.fr/airGR/?lang=en.",
sep = " ") sep = " ")
) )
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