R/read_out.R

In hzambran/hydroPSO: Particle Swarm Optimisation, with Focus on Environmental Models

# File read_out.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 2010-2020 Mauricio Zambrano-Bigiarini & Rodrigo Rojas
# Distributed under GPL 2 or later

################################################################################
#                            'read_out'                                        # 
################################################################################
# Purpose: This function reads the values of the objective function / model for#
#          each particle and iteration                                         #
################################################################################
# Output : list with two elements:                                             #
#          -) model.values: numeric or matrix/data.frame with the values of the#
#                           objective function / model for each particle and   #
#                           iteration                                          #
#          -) model.gofs  : numeric vector with the goodness-of-fit value for  #
#                           each particle and iteration                        #
#          -) model.best  : numeric with the best model output                 #
#          -) model.obs   : numeric with the observed values used for hydroPSO,#
#                           or NA if they were not provided                    #
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas                         #  
# Started: 08-Nov-2011,                                                        #
# Updates: 27-Jan-2012 ; 02-Feb-2012 ; 15-Feb-2012 ; 15-Feb-2012 ; 22-Feb-2012 #  
#          23-Mar-2012 ; 06-Dec-2012                                           # 
#          29-May-2013                                                         #  
#          09-Abr-2014                                                         # 
#          30-Jul-2015                                                         # 
#          28-Feb-2020 ; 07-Mar-2020 ; 12-Mar-2020                             #
#          22-Nov-2023                                                         #
################################################################################
# Columns in 'of_out' are:
# Iter         : integer, with the iteration number for each row of the file
# Part         : integer, with the particle number for each row of the file
# GoF          : numeric, with the goodness of fit value for each particle and iteration
# Model_Output : numeric vector (one or more values for each row), with the model output for each particle and iteration

read_out <- function(file="Model_out.txt", 
                     modelout.cols=NULL, 
                     nsim=NULL,
                     obs, 
                     MinMax=NULL, 
                     beh.thr=NA, 
                     verbose=TRUE,
                     #### Plotting arguments ####
                     plot=TRUE,
                     ptype=c("corr", "ts", "ecdf", "quant2ecdf"), 
                     ftype="dm", # OPTIONAL, only used when 'ptype=="ts"'.See [hydroGOF]{ggof}
                     FUN=mean,   # OPTIONAL, only used when 'ptype=="ts"'.See [hydroGOF]{ggof}                     
                     #### OPTIONS for ptype="quant2ecdf' #####
                     weights=NULL,
                     byrow=TRUE, 
                     quantiles.desired= c(0.05, 0.5, 0.95),
                     quantiles.labels= c("Q5", "Q50", "Q95"),
                     main=NULL,
                     ylab="Probability",
                     col="blue",
                     leg.cex=1.2,
                     leg.pos="bottomright",
                     cex.axis=1.2, 
                     cex.main=1.2, 
                     cex.lab=1.2,                     
                     #### PNG options ### 
                     do.png=FALSE, png.width=1500, png.height=900, png.res=90,
                     png.fname="ModelOut_vs_Obs.png" ) {                         
                         
  
  ##############################################################################
  # 1)                            Checkings                                    #
  ##############################################################################
  # Checking that 'file' exists
  if ( !file.exists(file) )
     stop( "Invalid argument value: The file '", basename(file), "' doesn't exist !!")

  # Checking 'MinMax'
  if ( !is.null(MinMax) ) {
     if ( !(MinMax %in% c("min", "max")) )
        stop("Invalid argument: 'MinMax' must be in c('min', 'max')")
  } else { # MinMax was not provided: it is read from the output 'PSO_logfile.txt' file
           # It assumes that 'PSO_logfile.txt' is located in the same directory than the 'Model_out.txt' file
      PSOlog <- read.table(file="PSO_logfile.txt", skip=9, nrows=1)
      MinMax <- as.character(PSOlog[1, 3])
      if (verbose) message("[ 'MinMax' was read from the 'PSO_logfile.txt' file, and set to '", MinMax, "' ]")
    } # ELSE end

  # Checking 'beh.thr'
  if ( !is.na(beh.thr) ) {
     if ( is.null(MinMax) ) {
        warning("Missing argument: 'MinMax' has to be provided before using 'beh.thr' !!")  
     } # IF end
  } # IF end
  
  # Checking 'plot' and 'MinMax'
  valid.types <- c("corr", "ts") 
  if (plot &  (length(which(!is.na(match(ptype, valid.types )))) <= 0) ) {
    if (is.null(MinMax)) {
      warning("Missing argument: 'MinMax' has to be provided before plotting when 'ptype' is in c('corr', 'ts') => 'plot=FALSE'")
      plot <- FALSE
    } # IF end
  } # IF end

  ##############################################################################
  # 2)                            Reading                                      #
  ##############################################################################

  # Reading the file
  if (verbose) message( "                                                     ")  
  if (verbose) message( "[ Reading the file '", basename(file), "' ... ]" )  
  if (!missing(nsim)) {
    cnames <- paste("sim", 1:nsim, sep="") 
    cnames <- c("Iter", "Part", "GoF", cnames)      
    sim  <- data.table::fread(file=file, header=FALSE, skip=1, fill=TRUE,  col.names=cnames, data.table=FALSE)
  } else sim  <- data.table::fread(file=file, header=FALSE, skip=1, fill=TRUE, data.table=FALSE)

  # Amount of total model outputs / parameter sets 
  nouts <- nrow(sim)
  if (verbose) message( "[ Total number of parameter sets: ", nouts, " ]" )       
  
  # computing the number of columns in 'file'
  ncols <- ncol(sim)  
  
  # Getting the GoF for each particle and iteration
  gofs <- sim[, 3]

  # Getting the Model Outputs for each particle and iteration
  outputs <- sim[, 4:ncols]
  
  # Computing the number of values in each model output
  lnsim  <- NCOL(outputs)
  if (verbose) message( "[ Number of model outputs for each parameter set ('nsim'): ", lnsim, " ]" )    
  
  # giving more meaningful names to the model outputs
  if (lnsim > 1)
    colnames(outputs) <- paste("Sim", 1:lnsim, sep="")
  
  # If the user only wants some columns of the model output file
  if (!is.null(modelout.cols)) {
    outputs <- outputs[, modelout.cols]
    
    lnsim  <- NCOL(outputs)
    if (verbose) message( "[ Number of extracted model outputs for each parameter set: ", lnsim, " ]" ) 
  } # IF end   
  
  # Filtering out those parameter sets above/below a certain threshold
  if (!is.na(beh.thr)) {  
     # Checking 'beh.thr'
     mx <- max(gofs, na.rm=TRUE)
      if (beh.thr > mx)
       stop("Invalid argument: 'beh.thr' must be lower than ", mx ,"!!")
    
    # Computing the row index of the behavioural parameter sets
    ifelse(MinMax=="min", beh.row.index <- which(gofs <= beh.thr), 
                          beh.row.index <- which(gofs >= beh.thr) )
    
    # Removing non-behavioural parameter sets and gofs
    if ( is.matrix(outputs) | is.data.frame(outputs) ) {
      outputs <- outputs[beh.row.index, ]
    } else outputs <- outputs[beh.row.index]
    gofs <- gofs[beh.row.index]
   
    # Amount of behavioural parameter sets 
    nbeh <- nrow(outputs)
    if (verbose) message( "[ Number of behavioural model outputs           : ", nbeh, " ]" ) 
    
  } # IF end
  
  # Selecting best model output
  if ( !is.null(MinMax) ) {
    ifelse(MinMax=="min", best.row.index <- which.min(gofs), 
                          best.row.index <- which.max(gofs) )
    if ( is.matrix(outputs) | is.data.frame(outputs) ) {
      best <- as.numeric(outputs[best.row.index, ]) 
    } else best <- as.numeric(outputs[best.row.index])
    
    # If the user only wants some columns of the model output file
    if (!is.null(modelout.cols)) best <- best[modelout.cols]
    
  } else {
      best <- NA
      plot <- FALSE
      if (verbose) warning("[ Note: You didn't provide 'MinMax' => model.best=NA & some plots will be skipped !! ]")
    } # ELSE end

  ##############################################################################
  # 3)                        Getting Observed Values                          #
  ##############################################################################

  if (missing(obs)) {
    fname <- "Observations.txt"
    if (file.exists(fname)) {
       if  ( length(find.package("zoo", quiet=TRUE)) != 0 ) {
        obs   <- read.zoo(fname) # zoo::read.zoo
        dates <- time(obs)
        # If the observed data do not have dates, they are transformed into numeric
        if ( TRUE && !( is(dates, "POSIXct") | is(dates, "POSIXt") | is(dates, "Date") ) ) obs <- coredata(obs)
      } else {
        obs <- read.table(fname, header=FALSE, skip=0)  
        # skippping the column with the index
        obs <- obs[,1]
      } # ELSE end 
      
      # If the user only wants some columns of the model output file
      if (!is.null(modelout.cols)) {
        obs   <- obs[modelout.cols]
        dates <- dates[modelout.cols] 
      } # IF end
  
    } else obs <- NA
  } # IF end

  # Checking length(obs)
  if ( !is.null(obs) & is.numeric(obs) ) {
    if ( length(obs) != lnsim ) 
      stop("Invalid argument: 'length(obs) != ncol(sims)' ", length(obs), "!=", lnsim, " !!")
  } # IF end
  
  ##############################################################################
  # 4)                            Output                                       #
  ##############################################################################  
  # creating the output
  out <- list(model.values=outputs, model.gofs=gofs, model.best=best, model.obs=obs)

  ##############################################################################
  # 5)                            Plotting                                     #
  ##############################################################################
  if ( plot & is.numeric(obs) ) {
  
       plot_out(sim=best, 
                obs=obs, 
                dates=dates, 
                ptype=ptype, 
                MinMax=MinMax, 
                ftype=ftype, # OPTIONAL, only used when 'ptype=="ts"'.See [hydroGOF]{ggof}
                FUN=FUN,     # OPTIONAL, only used when 'ptype=="ts"'.See [hydroGOF]{ggof}
                verbose=verbose,
                #### OPTIONS for ptype="quant2ecdf' #####
                weights=weights,
                byrow=byrow, 
                quantiles.desired= quantiles.desired,
                quantiles.labels= quantiles.labels,
                main=main,
                ylab=ylab,
                col=col,
                leg.cex=leg.cex,
                leg.pos=leg.pos,
                cex.axis=cex.axis, 
                cex.main=cex.main, 
                cex.lab=cex.lab,                     
                #### PNG options ### 
                do.png=do.png, 
                png.width=png.width, 
                png.height=png.height, 
                png.res=png.res,
                png.fname=png.fname)
  } # IF end
  
  return(out)
  
}  # 'read_out' END