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);
        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)==FALSE){ 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("\t Screening completed (", NRuns, " runs):")
      }
      if(Ncandidates==1){
        message("\t Starting point for steepest-descent local search:")
      }
      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))
    }
    ##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("\t Calibration completed (", NIter, " iterations, ", NRuns, " runs)")
      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))
    }
    ##Results_archiving_______________________________________________________
  	HistParamR <- cbind(HistParamR[1:NIter,]); colnames(HistParamR) <- paste("Param",1:NParam,sep="");
  	HistParamT <- cbind(HistParamT[1:NIter,]); colnames(HistParamT) <- paste("Param",1:NParam,sep="");
	  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(as.double(ParamFinalR),CritFinal*Multiplier,NIter,NRuns,HistParamR,HistCrit*Multiplier,MatBoolCrit,CritName,CritBestValue);
    names(OutputsCalib) <- c("ParamFinalR","CritFinal","NIter","NRuns","HistParamR","HistCrit","MatBoolCrit","CritName","CritBestValue");
    class(OutputsCalib) <- c("OutputsCalib","HBAN");
    return(OutputsCalib);



}