Forked from HYCAR-Hydro / airGR
Source project has a limited visibility.
Calibration_HBAN.R 22.26 KiB
#*************************************************************************************************
#' Calibration algorithm which minimises the error criterion. \cr
#' \cr
#' The algorithm is based on a local search procedure. 
#' First, a screening is performed using either a rough predefined grid or a list of parameter sets
#' and then a simple steepest descent local search algorithm is performed.
#' 
#' A screening is first performed either from a rough predefined grid (considering various initial 
#' values for each paramete) or from a list of initial parameter sets. \cr
#' The best set identified in this screening is then used as a starting point for the steepest 
#' descent local search algorithm. \cr
#' For this search, the parameters are used in a transformed version, to obtain uniform           
#' variation ranges (and thus a similar pace), while the true ranges might be quite different. \cr
#' At each iteration, we start from a parameter set of NParam values (NParam being the number of 
#' free parameters of the chosen hydrological model) and we determine the 2*NParam-1 new candidates 
#' by changing one by one the different parameters (+/- pace). \cr
#' All these candidates are tested and the best one kept to be the starting point for the next    
#' iteration. At the end of each iteration, the pace is either increased or decreased to adapt    
#' the progression speed. A diagonal progress can occasionally be done.                        \cr
#' The calibration algorithm stops when the pace becomes too small.                            \cr
#' To optimise the exploration of the parameter space, transformation functions are used to convert
#' the model parameters. This is done using the TransfoParam functions.
#*************************************************************************************************
#' @title  Calibration algorithm which minimises the error criterion using the Irstea-HBAN procedure
#' @author Laurent Coron (August 2013)
#' @references
#'   Michel, C. (1991),
#'        Hydrologie appliquée aux petits bassins ruraux, Hydrology handout (in French), Cemagref, Antony, France.
#' @example tests/example_Calibration_HBAN.R
#' @seealso \code{\link{Calibration}}, \code{\link{Calibration_optim}},
#'          \code{\link{RunModel_GR4J}}, \code{\link{TransfoParam_GR4J}}, \code{\link{ErrorCrit_RMSE}},
#'          \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, 
#'          \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
#' @encoding UTF-8
#' @export
#_FunctionInputs__________________________________________________________________________________
#' @param  InputsModel         [object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details
#' @param  RunOptions          [object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details
#' @param  InputsCrit          [object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details
#' @param  CalibOptions        [object of class \emph{CalibOptions}] see \code{\link{CreateCalibOptions}} for details
#' @param  FUN_MOD             [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
#' @param  FUN_CRIT            [function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)
#' @param  FUN_TRANSFO         (optional) [function] model parameters transformation function, if the FUN_MOD used is native in the package FUN_TRANSFO is automatically defined
#' @param  quiet               (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
#_FunctionOutputs_________________________________________________________________________________
#' @return  [list] list containing the function outputs organised as follows:
#'          \tabular{ll}{
#'          \emph{$ParamFinalR  }  \tab   [numeric] parameter set obtained at the end of the calibration \cr
#'          \emph{$CritFinal    }  \tab   [numeric] error criterion obtained at the end of the calibration \cr
#'          \emph{$NIter        }  \tab   [numeric] number of iterations during the calibration \cr
#'          \emph{$NRuns        }  \tab   [numeric] number of model runs done during the calibration \cr
#'          \emph{$HistParamR   }  \tab   [numeric] table showing the progression steps in the search for optimal set: parameter values \cr
#'          \emph{$HistCrit     }  \tab   [numeric] table showing the progression steps in the search for optimal set: criterion values \cr
#'          \emph{$MatBoolCrit  }  \tab   [boolean] table giving the requested and actual time steps when the model is calibrated \cr
#'          \emph{$CritName     }  \tab   [character] name of the calibration criterion \cr
#'          \emph{$CritBestValue}  \tab   [numeric] theoretical best criterion value \cr
#'          }
#**************************************************************************************************
Calibration_HBAN <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO=NULL,quiet=FALSE){
##_____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_HBAN is used \n"); return(NULL); }  
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
##_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_HBAN 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],
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
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_______________________________________________ if(PrefilteringType==1){ CandidatesParamR <- CalibOptions$StartParamList; } if(PrefilteringType==2){ DistribParamR <- CalibOptions$StartParamDistrib; DistribParamR[,!CalibOptions$OptimParam] <- NA; CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates; } ##Remplacement_of_non_optimised_values_____________________________________ CandidatesParamR <- apply(CandidatesParamR,1,function(x){ x[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$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(!quiet & Ncandidates>1){ if(PrefilteringType==1){ cat(paste("\t List-Screening in progress (",sep="")); } if(PrefilteringType==2){ cat(paste("\t Grid-Screening in progress (",sep="")); } cat("0%"); } for(iNew in 1:nrow(CandidatesParamR)){ if(!quiet & Ncandidates>1){ for(k in c(2,4,6,8)){ if(iNew==round(k/10*Ncandidates)){ cat(paste(" ",10*k,"%",sep="")); } } } ##Model_run Param <- CandidatesParamR[iNew,]; OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param); ##Calibration_criterion_computation OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel); 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(!quiet & Ncandidates>1){ cat(" 100%) \n"); } ##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(!quiet){ if(Ncandidates> 1){ cat(paste("\t Screening completed (",NRuns," runs): \n",sep="")); } if(Ncandidates==1){ cat(paste("\t Starting point for steepest-descent local search: \n",sep="")); } cat(paste("\t Param = ",paste(formatC(ParamStartR,format="f",width=8,digits=3),collapse=" , "),"\n",sep="")); cat(paste("\t Crit ",format(CritName,width=12,justify="left")," = ",formatC(CritStart*Multiplier,format="f",digits=4),"\n",sep="")); } ##Results_archiving________________________________________________________ HistParamR[1,] <- ParamStartR; HistParamT[1,] <- ParamStartT; HistCrit[1,] <- CritStart;
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
##_____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(!quiet){ cat("\t Steepest-descent local search in progress \n"); } 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,CalibOptions$OptimParam,Pace)$NewCandidatesT; CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR"); ##Remplacement_of_non_optimised_values_____________________________________ CandidatesParamR <- apply(CandidatesParamR,1,function(x){ x[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam]; return(x); }); if(NParam>1){ CandidatesParamR <- t(CandidatesParamR); } else { CandidatesParamR <- cbind(CandidatesParamR); }
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
##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); 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(CalibOptions$OptimParam[iC]==TRUE){ 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(CalibOptions$OptimParam[iC]==TRUE){ 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); 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); } }
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405
##Results_archiving_______________________________________________________ NewParamOptimR <- FUN_TRANSFO(NewParamOptimT,"TR"); HistParamR[ITER+1,] <- NewParamOptimR; HistParamT[ITER+1,] <- NewParamOptimT; HistCrit[ITER+1,] <- CritOptim; ### if(!quiet){ 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 & !quiet){ cat("\t No progress achieved \n"); } ##End_of_Steepest_Descent_Local_Search____________________________________ ParamFinalR <- NewParamOptimR; ParamFinalT <- NewParamOptimT; CritFinal <- CritOptim; NIter <- 1+ITER; if(!quiet){ cat(paste("\t Calibration completed (",NIter," iterations, ",NRuns," runs): \n",sep="")); cat(paste("\t Param = ",paste(formatC(ParamFinalR,format="f",width=8,digits=3),collapse=" , "),"\n",sep="")); cat(paste("\t Crit ",format(CritName,width=12,justify="left")," = ",formatC(CritFinal*Multiplier,format="f",digits=4),"\n",sep="")); } ##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); }