Calibration_Michel.R 16.26 KiB
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)
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
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) { Output <- list(NewCandidates = expand.grid(lapply(seq_len(ncol(DistribParamR)), function(x) DistribParam[, x]))) } ##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) {