Calibration_Michel.R

Calibration_Michel <- function(InputsModel, 
                               RunOptions, 
                               InputsCrit, 
                               CalibOptions,
                               FUN_MOD, 
                               FUN_CRIT,           # deprecated
                               FUN_TRANSFO = NULL, 
                               verbose = TRUE) {
  
  
  FUN_MOD  <- match.fun(FUN_MOD)
  if (!missing(FUN_CRIT)) {
    FUN_CRIT <- match.fun(FUN_CRIT)
  }
  if (!is.null(FUN_TRANSFO)) {
    FUN_TRANSFO <- match.fun(FUN_TRANSFO)
  }
  
  
  ##_____Arguments_check_____________________________________________________________________
  if (!inherits(InputsModel, "InputsModel")) {
    stop("'InputsModel' must be of class 'InputsModel'")
  }  
  if (!inherits(RunOptions, "RunOptions")) {
    stop("'RunOptions' must be of class 'RunOptions'")
  }  
  if (!inherits(InputsCrit, "InputsCrit")) {
    stop("'InputsCrit' must be of class 'InputsCrit'")
  }
  if (inherits(InputsCrit, "Multi")) {
    stop("'InputsCrit' must be of class 'Single' or 'Compo'")
  }
  if (inherits(InputsCrit, "Single")) {
    listVarObs <- InputsCrit$VarObs
  }
  if (inherits(InputsCrit, "Compo")) {
    listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
  }
  if ("SCA" %in% listVarObs & !"Gratio" %in% RunOptions$Outputs_Cal) {
    warning("Missing 'Gratio' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SCA")
    RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "Gratio")
  }
  if ("SWE" %in% listVarObs & !"SnowPack" %in% RunOptions$Outputs_Cal) {
    warning("Missing 'SnowPack' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SWE")
    RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "SnowPack")
  }
  if (!inherits(CalibOptions, "CalibOptions")) {
    stop("'CalibOptions' must be of class 'CalibOptions'")
  }  
  if (!inherits(CalibOptions, "HBAN")) {
    stop("'CalibOptions' must be of class 'HBAN' if 'Calibration_Michel' is used")
  }
  if (!missing(FUN_CRIT)) {
    warning("argument 'FUN_CRIT' is deprecated. The error criterion function is now automatically get from the 'InputsCrit' object")
  }  
  
  
  ##_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    )) {
      if (inherits(CalibOptions, "hysteresis")) {
        FUN_TRANSFO <- TransfoParam_CemaNeigeHyst
      } else {
        FUN_TRANSFO <- TransfoParam_CemaNeige
      }
    }
    if (identical(FUN_MOD, RunModel_CemaNeigeGR4H) |
        identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
      if (identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
        FUN1 <- TransfoParam_GR4H
      }
      if (identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
        FUN1 <- TransfoParam_GR4J
      }
      if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
        FUN1 <- TransfoParam_GR5J
      }
      if (identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
        FUN1 <- TransfoParam_GR6J
      }
      if (inherits(CalibOptions, "hysteresis")) {
        FUN2 <- TransfoParam_CemaNeigeHyst
      } else {
        FUN2 <- TransfoParam_CemaNeige
      }
      if (inherits(CalibOptions, "hysteresis")) {
        FUN_TRANSFO <- function(ParamIn, Direction) {
          Bool <- is.matrix(ParamIn)
          if (!Bool) {
            ParamIn <- rbind(ParamIn)
          }
          ParamOut <- NA * ParamIn
          NParam   <- ncol(ParamIn)
          ParamOut[,          1:(NParam-4)] <- FUN1(ParamIn[,          1:(NParam-4)], Direction)
          ParamOut[, (NParam-3):NParam    ] <- FUN2(ParamIn[, (NParam-3):NParam    ], Direction)
          if (!Bool) {
            ParamOut <- ParamOut[1, ]
          }
          return(ParamOut)
        }
      } else {
        FUN_TRANSFO <- function(ParamIn, Direction) {
          Bool <- is.matrix(ParamIn)
          if (!Bool) {
            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) {
            ParamOut <- ParamOut[1, ]
          }
          return(ParamOut)
        }
      }
      
    }
    if (is.null(FUN_TRANSFO)) {
      stop("'FUN_TRANSFO' was not found (in 'Calibration' function)")
    }
  }
  
  ##_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")
  }
  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) {
      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 <- ErrorCrit(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")
    }
    if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
      stop("each input set must have the same number of values")
    }
    ##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]) {
        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)) {
    
    
    ##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 <- ErrorCrit(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]) { 
          PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[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 <- ErrorCrit(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))
    if (inherits(InputsCrit, "Compo")) {
      listweights  <- OutputsCrit$CritCompo$MultiCritWeights
      listNameCrit <- OutputsCrit$CritCompo$MultiCritNames
      msgForm <- paste(sprintf("%.2f", listweights), listNameCrit, sep = " * ", collapse = ", ")
      msgForm <- unlist(strsplit(msgForm, split = ","))
      msgFormSep <- rep(c(",", ",", ",\n\t\t    "), times = ceiling(length(msgForm)/3))[1: length(msgForm)]
      msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
      msgForm <- gsub("\\,\\\n\\\t\\\t    $|\\,$", "", msgForm)
      message("\tFormula: sum(", msgForm, ")")
    }
  }
  ##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)
  
  
  
}