Reading/Plotting all the output files generated by ‘hydroPSO’

Description

The function read_results reads the following output files of hydroPSO:

1) ‘BestParameterSet.txt’: best parameter set and its corresponding goodness-of-fit found during the optimisation
2) ‘Particles.txt’: parameter values and their corresponding goodness-of-fit value for all particles and iterations
3) ‘Velocities.txt’: velocity values and their corresponding goodness-of-fit value for all particles and iterations
4) ‘Model_out.txt’: values of the objective function/model output for each particle and iteration
5) ‘ConvergenceMeasures.txt’: convergence measures summarizing performance of hydroPSO
6) ‘Particles_GofPerIter.txt’: goodness-of-fit only for all the particles during all the iterations

The function plot_results takes the outputs of the read_results function and then produces the following plots:

1) Dotty plots of parameter values
2) Histograms of parameter values
3) Boxplots of parameter values
4) Correlation matrix among parameter values (optional)
5) Empirical CDFs of parameter values
6) Parameter values vs Number of Model Evaluations
7) (pseudo) 3D dotty plots of (selected) parameter values
8) GoF for each particle against Number of Model Evaluations
9) Velocity values vs Number of Model Evaluations
10a) Scatterplot between Best Simulated values and Observations (OPTIONAL, only if MinMax is provided)
10b) Empirical CDFs for model's output (only produced if obs is NOT a zoo object)
10b) ggof (See ggof) between Best Simulated values and Observations (OPTIONAL, only if obs is a zoo object)
10d) Empirical CDFs for selected quantiles of model's output (OPTIONAL, only if obs is a zoo object)
11) Convergence Measures (Gbest and normSwarmRadius) vs Iteration Number

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
read_results(drty.out="PSO.out", MinMax=NULL, beh.thr=NA, 
             modelout.cols=NULL, nsim=NULL, verbose=TRUE)
           
plot_results(drty.out="PSO.out", param.names, gof.name="GoF", MinMax=NULL, 
     beh.thr=NA, beh.col="red", beh.lty=1, beh.lwd=2, nrows="auto", 
     col="black", ylab=gof.name, main=NULL, pch=19, cex=0.5, cex.main=1.7, 
     cex.axis=1.3, cex.lab=1.5, breaks="Scott", freq=TRUE, do.pairs=FALSE,
     weights=NULL, byrow=FALSE, leg.cex=1.2,
     
     
     dp3D.names="auto", GOFcuts="auto", 
     colorRamp= colorRampPalette(c("darkred", "red", "orange", "yellow", 
     "green", "darkgreen", "cyan")), alpha=0.65, points.cex=0.7, 
     
     ptype="one",
     
     nsim=NULL,
     
     modelout.cols=NULL,
     ftype="o", FUN=mean, 
     quantiles.desired= c(0.05,0.5,0.95), 
     quantiles.labels= c("Q5","Q50","Q95"),
     
     legend.pos="topright", 
     
     do.png=FALSE, png.width=1500, png.height=900, png.res=90, 
     dotty.png.fname="Params_DottyPlots.png", 
     hist.png.fname ="Params_Histograms.png",
     bxp.png.fname="Params_Boxplots.png",
     ecdf.png.fname ="Params_ECDFs.png", 
     pruns.png.fname="Params_ValuesPerRun.png",
     dp3d.png.fname ="Params_dp3d.png", 
     pairs.png.fname="Params_Pairs.png",
     part.png.fname ="Particles_GofPerIter.png",
     vruns.png.fname="Velocities_ValuePerRun.png", 
     modelout.best.png.fname="ModelOut_BestSim_vs_Obs.png",
     modelout.quant.png.fname="ModelOut_Quantiles.png",
     conv.png.fname ="ConvergenceMeasures.png", verbose=TRUE)

Arguments

drty.out

character, path to the directory storing the output files generated by hydroPSO

param.names

character, names for the parameters in params that have to be plotted (param.names can be a subset of params).
Names for each parameter are taken from the first row of the ‘Particles.txt’ file

verbose

logical, if TRUE, progress messages are printed

gof.name

character, name of the goodness-of-fit variable in all plots

MinMax

OPTIONAL. character, indicates whether the optimum value in x corresponds to the minimum or maximum of the objective function. It is only used to identify the optimum on the plots
Valid values are in: c('min', 'max')

beh.thr

OPTIONAL. numeric, threshold to filter out parameter sets and model outputs with a non-acceptable performance (non behavioural parameter sets)

nsim

OPTIONAL. number simulated equivalent values of the model / objective function to be compared against observations.
It is only useful when the model to be calibrated returns NA instead of the simulated values for some parameter set(s) (e.g., MODFLOW). See read_out

modelout.cols

numeric, column number in file that store the outputs that have to be read/plotted, without counting the first three that correspond to iteration, particle and GoF. If modelout.cols=NULL, all the columns in will be read, but the first trhee that contains the iteration number, the particle number and the corresponding goodness-of-fit. See read_out

beh.col

OPTIONAL. Only used when plot=TRUE
character, colour for drawing a horizontal line for separating behavioural from non behavioural parameter sets

beh.lty

OPTIONAL. Only used when plot=TRUE
numeric, line type for drawing a horizontal line for separating behavioural from non behavioural parameter sets

beh.lwd

OPTIONAL. Only used when plot=TRUE
numeric, width for drawing a horizontal line for separating behavioural from non behavioural parameter sets

nrows

OPTIONAL. Only used when plot=TRUE
numeric, number of rows to be used in the plotting window
If nrowsis set to auto, the number of rows is automatically computed depending on the number of columns of x

col

OPTIONAL. Only used when plot=TRUE
character, colour to be used for drawing the points of the dotty plots

ylab

OPTIONAL. Only used when plot=TRUE
character, label for the 'y' axis

main

OPTIONAL. Only used when plot=TRUE
character, title for the plot

pch

OPTIONAL. Only used when plot=TRUE
numeric, type of symbol to be used for drawing the points of the dotty plots. (e.g., 1: white circle)

cex

OPTIONAL. Only used when plot=TRUE
numeric, values controlling the size of text and points with respect to the default

cex.main

OPTIONAL. Only used when plot=TRUE
numeric, magnification for main titles relative to the current setting of cex

cex.axis

OPTIONAL. Only used when plot=TRUE
numeric, magnification for axis annotation relative to the current setting of cex

cex.lab

OPTIONAL. Only used when plot=TRUE
numeric, magnification for x and y labels relative to the current setting of cex

breaks

OPTIONAL. Only used when plot=TRUE
breaks for plotting the histograms of the parameter sets. See hist

freq

OPTIONAL. Only used when plot=TRUE
logical, if TRUE, the histogram graphic is a representation of frequencies, the counts component of the result; if FALSE, probability densities, component density, are plotted (so that the histogram has a total area of one). Defaults to TRUE if and only if breaks are equidistant (and probability is not specified). See hist

do.pairs

OPTIONAL. Only used when plot=TRUE
logical, indicates whether a correlation matrix among parameters has to be plotted. If the number of parameter sets tried during the optimisation is large, it may require some time.

weights

OPTIONAL. Only used when plot=TRUE
numeric vector, values of the weights to be used for computing the empirical CDFs. See params2ecdf

byrow

OPTIONAL. Only used when plot=TRUE
logical, indicates whether the computations have to be made for each column or for each row of x. See params2ecdf

leg.cex

OPTIONAL. Only used when plot=TRUE
character expansion factor *relative* to current 'par("cex")'. Used for text, and provides the default for 'pt.cex' and 'title.cex'. Default value = 1.2

dp3D.names

character, name for all the parameters (usually only the most sensitive ones) that will be used for plotting pseudo-3D dotty plots
If dp3D.names='auto' half the number of parameters in file are chosen randomly for plotting. See plot_NparOF

GOFcuts

numeric, specifies at which values of the objective function gof.name the colours of the plot have to change. See plot_NparOF

colorRamp

R function defining the colour ramp to be used for colouring the pseudo-3D dotty plots of Parameter Values, OR character representing those colours. See plot_NparOF

alpha

numeric between 0 and 1 representing the transparency level to apply to the colors of the pseudo-3D dotty plots. See plot_NparOF

points.cex

size of the points to be plotted

ptype

character, represents the type of plot. Valid values are: in c("one", "many"), for plotting all the particles in the same figure or in one windows per particle, respectively
See plot_GofPerParticle

ftype

OPTIONAL. Only used when plot=TRUE and the observed values provided by the user were zoo objects. See plot_out and ggof.

FUN

OPTIONAL. Only used when plot=TRUE and the observed values provided by the user were zoo objects. See plot_out and ggof

quantiles.desired

numeric vector, quantiles to be computed. Default values are c(.025, .5, .975) ( => 2.5%, 50%, 97.5% ). See plot_out

quantiles.labels

OPTIONAL. Only used when plot=TRUE
character vector, names to quantiles.desired. Default value is c("Q5", "Q50", "Q95"). See plot_out

legend.pos

See plot_convergence

do.png

logical, indicates if all the figures have to be saved into PNG files instead of the screen device

png.width

OPTIONAL. Only used when do.png=TRUE
numeric, width of the PNG device. See png

png.height

OPTIONAL. Only used when do.png=TRUE
numeric, height of the PNG device. See png

png.res

OPTIONAL. Only used when do.png=TRUE
numeric, nominal resolution in ppi which will be recorded in the PNG file, if a positive integer of the device. See png

dotty.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the dotty plots of the parameter values.

hist.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the histograms of the parameter values.

bxp.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the boxplots of the parameter values

ecdf.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the empirical CDFs of the parameter values.

pruns.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the parameter values vs the number of model evaluations

dp3d.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the pseudo-3D plots of all the parameters defined in dp3D.names

pairs.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the correlation matrix among the parameters and goodness-of-fit values in params and gofs. See plot_particles and hydropairs

part.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the goodness-of-fit for all the particles along the iterations

vruns.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the velocity values vs the number of model evaluations

modelout.best.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the observed values against its best simulated counterpart. See plot_out

modelout.quant.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with some quantiles of simulated values against its observed counterparts. See plot_out

conv.png.fname

OPTIONAL. Only used when do.png=TRUE
character, filename used to store the PNG file with the convergence measures. See plot_convergence

Value

The function read_results returns a list with the following elements:

best.param

numeric with the best parameter set

best.gof

numeric with the best fitness value of the objective function

params

data.frame with all the parameter sets tested during the optimisation

gofs

numeric with all the fitness values computed during the optimisation (each element in gofs corresponds to one row of params)

model.values

numeric or matrix/data.frame with the values of the objective function / model for each particle and iteration. See read_out

model.best

numeric with the best model / objective function value. In order to be computed, the user has to provide a valid value for MinMax. See read_out

model.obs

numeric with the observed values used during the optimisation. See obs

convergence.measures

matrix/data.frame with the convergence measures. See read_convergence function

part.GofPerIter

matrix/data.frame with the goodness-of-fit values for all the particles during all the iterations. It has as many columns as parameters to be optimised and as many rows as the number of iterations effectively carried out

Author(s)

Mauricio Zambrano-Bigiarini, mzb.devel@gmail.com

See Also

hydroPSO, read_best, read_particles, read_velocities, read_out, read_convergence, read_GofPerParticle, plot_ParamsPerIter

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
## Not run: 
# Setting the user home directory as working directory
setwd("~")

# Number of dimensions to be optimised
D <- 5

# Boundaries of the search space (Ackley test function)
lower <- rep(-32, D)
upper <- rep(32, D)

# Setting the seed
set.seed(100)

# Running PSO with the 'ackley' test function, writing the results to text files
hydroPSO(fn=ackley, lower=lower, upper=upper)
  
# Reading all the results and storing them in a variable
res <- read_results()

# Plotting all the results with a goodness-of-fit value lower than 5
plot_results(MinMax="min", beh.thr=5)

## End(Not run)

## Not run: 
################################################################################
#######################  SPSO-2007 example START ###############################
################################################################################
# Number of dimensions to be optimised
D <- 10

# boundaries for the test function
lower <- rep(-100, D)      # sphere
#lower <- rep(-5.12, D)    # rastrigin
#lower <- rep(-32, D)      # ackley

fn <- sphere
#fn <-rastrigin
#fn <-ackley

#######################################
#####   SPSO-2007 parameters   ########
npart  <- 10+floor(2*sqrt(D))
c1     <- 0.5+log(2)
c2     <- 0.5+log(2)
abstol <- 1e-20 
reltol <- 1e-20
maxit  <- 1000

use.IW        <- TRUE
IW.w          <- 1/(2*log(2))
REPORT        <- 100
lambda        <- 1
boundary.wall <- "absorbing2007"
#######################################

# Setting the user home directory as working directory
setwd("~")

# Running PSO  and writing the results to text files
set.seed(100)

hydroPSO(fn= fn, method="spso2007", lower=lower, upper=-lower,
         control=list(MinMax="min", maxit=maxit, npart=npart,  
                      c1=c1, c2=c2, 
                      use.IW=use.IW, IW.w=IW.w, 
                      topology="random", lambda=lambda, K=3,
                      Xini.type="random", Vini.type="random2007",
                      best.update="sync",
                      boundary.wall=boundary.wall, 
                      write2disk=TRUE, plot=FALSE, REPORT=REPORT,
                      abstol=abstol, reltol=reltol 
                      )
         )

# Plotting all the results 
plot_results(MinMax="min")

################################################################################
#######################  SPSO-2007 example END   ###############################
################################################################################


################################################################################
############### recommended hydroPSO configuration - START #####################
################################################################################

# Running PSO  and writing the results to text files
set.seed(100)
hydroPSO(fn= fn, method="spso2011", lower=lower, upper=-lower,
         control=list(MinMax="min", maxit=maxit, npart=40, 
                      c1=2.05, c2=2.05, 
                      use.IW=FALSE, use.CF=TRUE, 
                      topology="random", K=11,
                      use.TVlambda=TRUE, TVlambda.rng=c(1, 0.5), 
                      Xini.type="lhs", Vini.type="lhs2011",
                      best.update="sync",
                      boundary.wall="absorbing2011", 
                      write2disk=FALSE, plot=FALSE, REPORT=REPORT,
                      abstol=abstol, reltol=reltol 
                      )
         )

# compare the final optimum value and the number of function calls with those
# obtained in the SPSO-2007 example

################################################################################
################ recommended hydroPSO configuration - END ######################
################################################################################

## End(Not run)