Calibration_Michel.R

Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO=NULL, verbose = TRUE){


##_____Arguments_check_____________________________________________________________________
    if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
    if(inherits(RunOptions,"RunOptions")==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }  
    if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); 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); }  


   ##_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);