R/verification.R

Defines functions verification

Documented in verification

## File verification.R
# Part of the hydroPSO R package, https://github.com/hzambran/hydroPSO
#                                 http://cran.r-project.org/web/packages/hydroPSO
#                                 http://www.rforge.net/hydroPSO/
## Copyright 2011-2020 Mauricio Zambrano-Bigiarini & Rodrigo Rojas
## Distributed under GPL 2 or later

################################################################################
#                           verification                                       #
################################################################################
# Purpose    : Run a hydrological model with different parameter sets (obtained#
#              during calibration) and to obtain the simulated time series for #
#              each one of the parameter sets                                  #
################################################################################
# Output     : A list of two elements:                                         #
#              ##1) sim: simulated values obtained by running the hydrological #
#                      model                                                   #
#              2) gofs    : goodness-of fit of the simulated values against    #
#                           observed ones, by using THE USER-DEFINED 'gof'     #
#                           measure                                            #
#              3) best.gof: goodness-of fit of the "best" parameter set        #
#              4) best.par: parameter values of the "best" paraemter set       # 
################################################################################
# Author  : Mauricio Zambrano-Bigiarini                                        #
# Started : 18-Jan-2011 at JRC Ispra                                           #
# Updates : 12-May-2011 ; 13-Feb-2012  ; 23-Feb-2012                           #
#           09-Abr-2014                                                        #
#           09-Mar-2020 ; 12-Mar-2020 ; 15-Mar-2020                            #
################################################################################
verification <- function(
                         fn="hydromod",  
                         par,
                         ...,
                         control=list(),
                         model.FUN=NULL,
                         model.FUN.args=list()                                
                        ) {
                     
  
  ##############################################################################
  #                            INPUT CHECKING                                  #
  ##############################################################################
        
  # Checking the name of the objective function
  if (missing(fn)) {
    stop("Missing argument: 'fn' must be provided")
  } else 
      if ( is.character(fn) | is.function(fn) )  {
        if (is.character(fn)) {
          if (fn=="hydromod") {
            fn.name <- fn
	    fn      <- match.fun(fn)
          } else if (fn=="hydromodInR") {
              fn.name <- fn
              fn      <- match.fun(model.FUN)
            } else stop("Invalid argument: valid character values for 'fn' are only: c('hydromod', 'hydromodInR')")
	} else if (is.function(fn)) {
	    fn.name <- as.character(substitute(fn))
	    fn      <- fn
	  } # ELSE end
      } else stop("Missing argument: 'class(fn)' must be in c('function', 'character')")      
        
  # Checking 'par'
  if (missing(par)) stop("Missing argument: 'par' must be provided")
                 
  ########################################################################        
  con <- list(
                
          drty.in=getwd(),
          drty.out="verification", # Character, with the name of the directory that will store the results of the LH-OAT. 
          digits=7,
                
          gof.name="GoF",          # Character, only used for identifying the goodness-of-fit of each model run
          MinMax=c("min", "max"),            # Character, indicating if PSO have to find a minimum or a maximum for the objective function. \cr
                                   # Valid values are in: \code{c('min', 'max')} \cr
          do.plots=FALSE,
          write2disk=TRUE,
          verbose= TRUE,           # logical, indicating if progress messages have to be printed
          
          parallel=c("none", "multicore", "parallel", "parallelWin"),
          par.nnodes=NA,
	  par.pkgs= c()
             )
             
  MinMax        <- match.arg(control[["MinMax"]], con[["MinMax"]])     
  parallel      <- match.arg(control[["parallel"]], con[["parallel"]])    

  nmsC <- names(con)
  con[(namc <- names(control))] <- control
  
  if (length(noNms <- namc[!namc %in% nmsC])) 
    warning("[Unknown names in control: ", paste(noNms, collapse = ", "), " (not used) !]")
    
  drty.in        <- con[["drty.in"]]
  drty.out       <- con[["drty.out"]]
  digits         <- con[["digits"]]

  gof.name       <- con[["gof.name"]]
  MinMax         <- con[["MinMax"]]
  do.plots       <- as.logical(con[["do.plots"]])  
  write2disk     <- as.logical(con[["write2disk"]])
  verbose        <- as.logical(con[["verbose"]])
  par.nnodes     <- con[["par.nnodes"]]
  par.pkgs       <- con[["par.pkgs"]] 
        
  ########################################################################
  ##################### Dummy checkings ##################################
     
  # 'hydromod' checkings
  if (fn.name=="hydromod") {

    # checking that 'model.FUN' is a valid function          
    if ( is.null(model.FUN) ) {
      stop( "'model.FUN' has to be defined !" )
    } else  {
        model.FUN.name <- model.FUN
        model.FUN      <- match.fun(model.FUN)
      } # ELSE end
  
    # checking 'model.FUN.args'
    if ( length(model.FUN.args)==0 ) {
      warning( "'model.FUN.args' is an empty list. Are you sure your model doesn't have any argument(s) ?" )
    } else {
        # Assigning the numeric values provided by the user in ' model.FUN.args' 
        # to the arguments of the hydrological model
        model.FUN.argsDefaults <- formals(model.FUN)
        model.FUN.args         <- modifyList(model.FUN.argsDefaults, model.FUN.args) 
      } # ELSe end
             
  } # IF end 

  if (fn.name=="hydromodInR") {
    if ( is.null(model.FUN) ) {
      stop( "'model.FUN' has to be defined !" )
    } else  {
        model.FUN.name <- as.character(substitute(model.FUN))
        model.FUN      <- match.fun(model.FUN)   
      } # ELSE end

    if (!("param.values" %in% names(formals(model.FUN)) ))
      stop("[ Invalid argument: 'param.values' must be the first argument of the 'model.FUN' function! ]")

    if (!("obs" %in% names(formals(model.FUN)) )) 
      stop("[ Invalid argument: 'obs' must be an argument of the 'model.FUN' function! ]")
   
    model.FUN.argsDefaults <- formals(model.FUN)
    if ( length(model.FUN.args) > 0 ) {
      model.FUN.args <- modifyList(model.FUN.argsDefaults, model.FUN.args) 
    } else model.FUN.args <- model.FUN.argsDefaults

  } # IF end 
          
  if ( is.matrix(par) | is.data.frame(par) ) {
    # Computing 'P', the Dimension of the Solution Space
    P <- ncol(par)

    # Computing the number of parameter sets
    nparamsets <- nrow(par)

    # Meaningful name of each one of the parameters
    if (is.null(colnames(par))) {
      Parameter.names <- paste("Param", 1:nparamsets, sep="")
    } else Parameter.names <- colnames(par)  
  
  } else if (is.numeric(par)) {
      # Computing 'P', the Dimension of the Solution Space
      P <- length(par)
 
      # Computing the number of parameter sets
      nparamsets <- 1

      # Meaningful name of each one of the parameters
      if (is.null(names(par))) {
        Parameter.names <- paste("Param", 1:P, sep="")
      } else Parameter.names <- names(par)    

    } else stop("Invalid argument: 'class(par)' must be in c('numeric', 'matrix', 'data.frame')")
   
  if (verbose) message( "[ Number of parameter sets read: ", nparamsets, "            ]" )   
  if (verbose) message( "[ Parameter names              : ", paste(Parameter.names, collapse = ", "), " ]")
  
  # If the user only provided a single parameter set, it is transformed into matrix
  if (nparamsets==1) par <- matrix(par, nrow=1)
        
  # Adding the parent path of 'drty.out', if it doesn't have it
  if (drty.out == basename(drty.out) )
    drty.out <- paste( getwd(), "/", drty.out, sep="")
        
  # Verifying that 'drty.out' directory exists. IF not, it is created
  if (!file.exists(file.path(drty.out))) {
    if (write2disk) {
      dir.create(file.path(drty.out))
      if (verbose) message("                                            ")
      if (verbose) message("[ Output directory '", basename(drty.out), "' was created on: '", dirname(drty.out), "' ]") 
      if (verbose) message("                                            ")
    } # IF end
  } # IF end  
  
  
  ##############################################################################
  #                            Writing Info File
  ##############################################################################  
     
  # File 'Verification-logfile.txt' #        
  InfoTXT.fname <- paste(file.path(drty.out), "/", "Verification-logfile.txt", sep="")
  InfoTXT.TextFile  <- file(InfoTXT.fname , "w+")
  #c(isOpen(Tfile, "r"), isOpen(Tfile, "w")) # both TRUE
  writeLines("================================================================================", InfoTXT.TextFile) # writing a separation line with a carriage return
    writeLines(c("Platform             :", sessionInfo()[[1]]$platform), InfoTXT.TextFile, sep="  ")
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines(c("R version            :", sessionInfo()[[1]]$version.string), InfoTXT.TextFile, sep="  ")
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines(c("hydroPSO version     :", sessionInfo()$otherPkgs$hydroPSO$Version), InfoTXT.TextFile, sep="  ")
    writeLines("", InfoTXT.TextFile) #writing a blank line with a carriage return
    writeLines(c("hydroPSO Built       :", sessionInfo()$otherPkgs$hydroPSO$Built), InfoTXT.TextFile, sep="  ")
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines("================================================================================", InfoTXT.TextFile) # writing a separation line with a carriage return
    Time.Ini <- Sys.time()
    writeLines(c("Starting Time        :", date()), InfoTXT.TextFile, sep="  ")
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines("================================================================================", InfoTXT.TextFile) # writing a separation line with a carriage return
    writeLines(c("drty.in              :", drty.in), InfoTXT.TextFile, sep=" ") 
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines(c("drty.out             :", drty.out), InfoTXT.TextFile, sep=" ") 
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines(c("digits               :", digits), InfoTXT.TextFile, sep=" ") 
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines(c("gof.name             :", gof.name), InfoTXT.TextFile, sep=" ") 
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines(c("do.plots             :", do.plots), InfoTXT.TextFile, sep=" ") 
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines(c("write2disk           :", write2disk), InfoTXT.TextFile, sep=" ") 
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines(c("verbose              :", verbose), InfoTXT.TextFile, sep=" ") 
    writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
    writeLines(c("parallel          :", parallel), InfoTXT.TextFile, sep=" ")  
    writeLines("", InfoTXT.TextFile)  
    if (parallel!="none") {
      writeLines(c("par.nnodes        :", par.nnodes), InfoTXT.TextFile, sep=" ") 
      writeLines("", InfoTXT.TextFile)
      writeLines(c("par.pkgs          :", par.pkgs), InfoTXT.TextFile, sep=" ") 
      writeLines("", InfoTXT.TextFile)     
    } # IF end
    if (fn.name=="hydromod") {
      try(writeLines(c("hydromod function    :", model.FUN.name), InfoTXT.TextFile, sep=" ") , TRUE)
      writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
      writeLines(c("hydromod args        :"), InfoTXT.TextFile, sep=" ") 
      writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
      for ( i in 1:length(model.FUN.args) ) {              
        arg.name  <- names(model.FUN.args)[i]
        arg.name  <- format(paste("  ", arg.name, sep=""), width=20, justify="left" )
        arg.value <- ""
        arg.value <- try(as.character( as.character(model.FUN.args[i])), TRUE)
             
        writeLines(c(arg.name, ":", arg.value), InfoTXT.TextFile, sep=" ") 
        writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return             
      } # FOR end
  } # IF end
  # Closing the text file
  close(InfoTXT.TextFile) 
  
  ########################################################################  
  #                        Text Files initialization                     #
  ########################################################################  

  # File 'Verification-ModelOut.txt' #
  model.out.text.fname <- paste(file.path(drty.out), "/", "Verification-ModelOut.txt", sep="")
  OFout.Text.file  <- file(model.out.text.fname, "w+")
  #c(isOpen(Tfile, "r"), isOpen(Tfile, "w")) # both TRUE
  #writeLines( paste("t", 1:length(obs)), OFout.Text.file, sep=" ")   
  writeLines(c("Param", "GoF", "Model_Output"), OFout.Text.file, sep="  ") 
  writeLines("", OFout.Text.file) # writing a blank line with a carriage return
  close(OFout.Text.file)  
        

  # File 'Verification-ParamValues.txt' #
  # with the parameters values for each partcile in each iteration
  gof.text.fname <- paste(file.path(drty.out), "/", "Verification-ParamValues.txt", sep="")
  Params.Text.file  <- file(gof.text.fname, "w+")           
  #c(isOpen(Tfile, "r"), isOpen(Tfile, "w")) # both TRUE
  writeLines(c("ParamNmbr", gof.name, Parameter.names), Params.Text.file, sep=" ")
  writeLines("", Params.Text.file) # writing a blank line with a carriage return
  close(Params.Text.file)          
        

  ############################################################################
  ##                                parallel                                 #
  ############################################################################
  if (parallel != "none") {
    
  #  if ( ( (parallel=="multicore") | (parallel=="parallel") ) & 
  #     ( (R.version$os=="mingw32") | (R.version$os=="mingw64") ) )
  #     stop("[ Fork clusters are not supported on Windows =>  'parallel' can not be set to '", parallel, "' ]")

  if (parallel=="multicore") {
     warning("[ Package 'parallel' is not available anymore in CRAN. It was changed to 'parallel='parallel' ]")
     parallel <- "parallel"
  } # IF end
    
  ifelse(parallel=="parallelWin", parallel.pkg <- "parallel",  parallel.pkg <- parallel) 
  if ( length(find.package(parallel.pkg, quiet=TRUE)) == 0 ) {               
    warning("[ Package '", parallel.pkg, "' is not installed =>  parallel='none' ]")
    parallel <- "none"
  }  else { 
      
       if (verbose) message("                                                 ")
       if (verbose) message("[ Parallel initialization ...                   ]")
      
       fn1 <- function(i, x) fn(x[i,])

       #require(parallel)           
       nnodes.pc <- parallel::detectCores()
       if (verbose) message("[ Number of cores/nodes detected: ", nnodes.pc, "             ]")
           
       if ( (parallel=="parallel") | (parallel=="parallelWin") ) {             
          logfile.fname <- paste(file.path(drty.out), "/", "parallel_logfile.txt", sep="") 
          if (file.exists(logfile.fname)) file.remove(logfile.fname)
       } # IF end
             
       if (is.na(par.nnodes)) {
         par.nnodes <- nnodes.pc
       } else if (par.nnodes > nnodes.pc) {
           warning("[ 'nnodes' > number of detected cores (", par.nnodes, ">", nnodes.pc, ") =>  par.nnodes=", nnodes.pc, "            ] !",)
           par.nnodes <- nnodes.pc
         } # ELSE end
 
       if (verbose) message("[ Number of cores/nodes used    : ", par.nnodes, "             ]")                 
               
       if (parallel=="parallel") {
           ifelse(write2disk, 
                  cl <- parallel::makeForkCluster(nnodes = par.nnodes, outfile=logfile.fname),
                  cl <- parallel::makeForkCluster(nnodes = par.nnodes) )         
       } else if (parallel=="parallelWin") {      
           ifelse(write2disk,
                  cl <- parallel::makeCluster(par.nnodes, outfile=logfile.fname),
                  cl <- parallel::makeCluster(par.nnodes) )
           pckgFn <- function(packages) {
             for(i in packages) library(i, character.only = TRUE)
           } # 'packFn' END
           parallel::clusterCall(cl, pckgFn, par.pkgs)
           parallel::clusterExport(cl, ls.str(mode="function",envir=.GlobalEnv) )
           if (fn.name=="hydromod") {
             parallel::clusterExport(cl, model.FUN.args$out.FUN)
             parallel::clusterExport(cl, model.FUN.args$gof.FUN)
           } # IF end                   
         } # ELSE end                   
                            
         if (fn.name=="hydromod") {
           if (!("model.drty" %in% names(formals(hydromod)) )) {
              stop("[ Invalid argument: 'model.drty' has to be an argument of the 'hydromod' function! ]")
           } else { # Copying the files in 'model.drty' as many times as the number of cores
             
               model.drty <- path.expand(model.FUN.args$model.drty)
                 
               files <- list.files(model.drty, full.names=TRUE, include.dirs=TRUE) 
               tmp <- which(basename(files)=="parallel")
               if (length(tmp) > 0) files <- files[-tmp]
               parallel.drty <- paste(file.path(model.drty), "/parallel", sep="")

               if (file.exists(parallel.drty)) {                      
                 if (verbose) message("[ Removing the 'parallel' directory ... ]")    
                 try(unlink(parallel.drty, recursive=TRUE, force=TRUE))
               } # IF end 
               dir.create(parallel.drty)

               mc.dirs <- character(par.nnodes)
               if (verbose) message("                                                     ")
               for (i in 1:par.nnodes) {
                 mc.dirs[i] <- paste(parallel.drty, "/", i, "/", sep="")
                 dir.create(mc.dirs[i])
                 if (verbose) message("[ Copying model input files to directory '", mc.dirs[i], "' ... ]")
                 file.copy(from=files, to=mc.dirs[i], overwrite=TRUE, recursive=TRUE)
               } # FOR end
                 
               tmp       <- ceiling(nparamsets/par.nnodes)        
               part.dirs <- rep(mc.dirs, tmp)[1:nparamsets]  
             } # ELSE end                 
         } # IF end
           
       } # ELSE end  
  
  }  # IF end    
  ############################################################################## 

  ##############################################################################
  #                            MAIN BODY
  ##############################################################################
  
  # Opening the file 'Verification-ModelOut.txt' for appending
  OFout.Text.file <- file(model.out.text.fname, "a")   

  # Opening the file 'Verification-ParamValues.txt' for appending
  Params.Text.file <- file(gof.text.fname, "a")

  gof.all <- numeric(nparamsets)


  # Evaluating an R Function
  if ( (fn.name != "hydromod") & (fn.name != "hydromodInR") ) {          
    if (parallel=="none") {
      hydromod.out <- apply(par, fn, MARGIN=1, ...)
    } else             
        if ( (parallel=="parallel") | (parallel=="parallelWin") ) #{
           hydromod.out <- parallel::parRapply(cl= cl, x=par, FUN=fn, ...)
        #} else if (parallel=="multicore")
        #    hydromod.out <- unlist(parallel::mclapply(1:npart, FUN=fn1, x=par, ..., mc.cores=par.nnodes, mc.silent=TRUE)) 
  } # IF end

  # Evaluating an R-based model
  if (fn.name == "hydromodInR") {
    if (verbose) message("                                                 ")
    if (verbose) message("=================================================")
    if (verbose) message("[ Running the model ...                         ]") 
    if (parallel=="none") {
      hydromod.out <- apply(par, model.FUN, MARGIN=1)
    } else             
        if ( (parallel=="parallel") | (parallel=="parallelWin") ) #{
          hydromod.out <- parallel::parRapply(cl= cl, x=par, FUN=model.FUN)
        #} else if (parallel=="multicore")
        #  hydromod.out <- unlist(parallel::mclapply(1:npart, FUN=fn1, x=par, mc.cores=par.nnodes, mc.silent=TRUE)) )
  } # IF end	 


  # Evaluating a system-console-based model
  for ( p in 1:nparamsets) {  

    # Getting the parameter set
    param.values.p <- as.numeric(par[p,])
    
    ##########################################################################
    # 2)                 Running the hydrological model                      #
    ##########################################################################
    
    # If the user wants to create a plot with obs vs sim
    if (do.plots) {
      do.png         <- TRUE
      png.fname      <- paste(drty.out, "/ParameterSet_", p, ".png", sep="")
      main           <- paste("Parameter Set:", p, sep="")
      model.FUN.args <- modifyList(model.FUN.args, list(do.png=do.png, png.fname=png.fname, main=main)) 
    } # IF end
    
    # Model evaluation 
    if (fn.name == "hydromod") {

      if (verbose) message("                    |                      ")
      if (verbose) message("                    |                      ") 
      if (verbose) message("==============================================================")
      if (verbose) message( paste("[ Running parameter set ", 
                                  format( p, width=4, justify="left" ), 
                                  "/", nparamsets, " => ", 
                                  format( round(100*p/nparamsets,2), width=7, justify="left" ), "%",
                                  ". Starting... ]", sep="") )
      if (verbose) message("==============================================================")

      # Running the console-based model
      model.FUN.args <- modifyList(model.FUN.args, list(param.values=param.values.p)) 
      hydromod.out   <- do.call(model.FUN, as.list(model.FUN.args)) 

      if (verbose) message("                              |                               ")
      if (verbose) message("                              |                               ") 
      if (verbose) message("==============================================================")
      if (verbose) message("[================ Verification finished ! ===================]")
      if (verbose) message("==============================================================")
    } # IF end
     
        
    ############################################################################
    # 3)                     Extracting simulated values                       #                                 
    ############################################################################
                  
    # Extracting the simulated values and the goodness of fit
    if ( fn.name == "hydromod" ) {
      sims <- as.numeric(hydromod.out[["sim"]])
      GoF  <- as.numeric(hydromod.out[["GoF"]])
    } else if ( fn.name == "hydromodInR" ) {
        sims  <- hydromod.out[[p]][["sim"]]
        GoF   <- hydromod.out[[p]][["GoF"]] 
      } else {
          sims <- as.numeric(hydromod.out[p])
          GoF  <- as.numeric(hydromod.out[p])
        } # ELSE end     
   
    gof.all[p] <- GoF

    # Writing to the 'Verification-ModelOut.txt' file
    suppressWarnings( tmp <- formatC(GoF, format="E", digits=digits, flag=" ") )
    if ( (fn.name != "hydromod") & (fn.name != "hydromodInR") ) {
      writeLines(as.character(c(p, tmp, tmp)), OFout.Text.file, sep="  ") 
    } else suppressWarnings( writeLines(as.character(c(p, tmp, formatC(sims, format="E", digits=digits, flag=" ") )), OFout.Text.file, sep="  ") )
    writeLines("", OFout.Text.file) # writing a blank line with a carriage return  
    
    # Writing to the 'Verification-ParamValues.txt' file
    suppressWarnings( writeLines(as.character(c(p, tmp, formatC(param.values.p, format="E", digits=digits, flag=" ") )), Params.Text.file, sep="  ") )
    writeLines("", Params.Text.file) # writing a blank line with a carriage return 
    
  } # FOR end

  # Closing the output text files
  close(OFout.Text.file)
  close(Params.Text.file)
  
  ##############################################################################
  # 4)                    Writing Ending and Elapsed Time                      #                                 
  ##############################################################################
  InfoTXT.TextFile <- file(InfoTXT.fname, "a")    
  #c(isOpen(Tfile, "r"), isOpen(Tfile, "w")) # both TRUE
  writeLines("================================================================================", InfoTXT.TextFile) # writing a separation line with a carriage return
  writeLines(c("Ending Time            :", date()), InfoTXT.TextFile, sep=" ")
  writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
  writeLines("================================================================================", InfoTXT.TextFile) # writing a separation line with a carriage return
  # Total time of the simulations
  Time.Fin <- Sys.time()
  writeLines(c("Elapsed Time           :", format(round(Time.Fin - Time.Ini, 2))), InfoTXT.TextFile, sep=" ")
  writeLines("", InfoTXT.TextFile) # writing a blank line with a carriage return
  writeLines("================================================================================", InfoTXT.TextFile) # writing a separation line with a carriage return
  close(InfoTXT.TextFile)

  ##############################################################################
  ##                                parallel                                   #
  ##############################################################################
  if (parallel!="none") {
    if ( (parallel=="parallel") | (parallel=="parallelWin") )   
         parallel::stopCluster(cl)   
    if (fn.name=="hydromod") {
      if (verbose) message("                                         ")
      if (verbose) message("[ Removing the 'parallel' directory ... ]")    
      unlink(dirname(mc.dirs[1]), recursive=TRUE)
    } # IF end
         
  } # IF end
  ##############################################################################


  
  ##############################################################################
  # 5)                    Creating the output                                  #                                 
  ##############################################################################

  if (verbose) message("=================================================")
  if (verbose) message("[ Reading 'Verification-ModelOut.txt' file ...  ]") 
  sims <- data.table::fread(file=model.out.text.fname, skip=1, data.table=FALSE)
  nc   <- ncol(sims)
  # Removing the first 2 columns in 'sims': ParameterSetNmbr, GoF
  sims <- sims[, 3:nc]  
  colnames(sims) <- paste0("sim", 1:(nc-2))
  colnames(sims) <- paste0("par", 1:nparamsets)


  ifelse(MinMax=="min", best.rowindex <- which.min(gof.all),
                        best.rowindex <- which.max(gof.all)) 

  best.gof <- gof.all[best.rowindex]                       
  best.par <- par[best.rowindex,]

  out <- list(gofs=gof.all, model.values=sims, best.gof=best.gof, best.param=best.par )

  if (verbose) message("=================================================")
  if (verbose) message("[                Finished !                     ]") 
  if (verbose) message("=================================================")

  return(out)

} # 'verification' END

Try the hydroPSO package in your browser

Any scripts or data that you put into this service are public.

hydroPSO documentation built on April 29, 2020, 9:37 a.m.