Skip to content
Snippets Groups Projects
Commit 269eae6c authored by Mauricio Zambrano-Bigiarini's avatar Mauricio Zambrano-Bigiarini
Browse files

new argument 'do.pairs' for functions: 'read\_particles', 'plot\_particles', 'plot\_results'

parent 91539b52
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
+ 1
− 1
View file @ 269eae6c
Package: hydroPSO
Type: Package
Title: Model-Independent Particle Swarm Optimisation for Environmental Models
Version: 0.1-58-22
Version: 0.1-58-23
Date: 2012-11-12
Authors@R: c(person("Mauricio", "Zambrano-Bigiarini", email="mzb.devel@gmail.com", role=c("aut","cre")), person("Rodrigo", "Rojas", email="Rodrigo.RojasMujica@gmail.com", role=c("ctb")) )
Maintainer: Mauricio Zambrano-Bigiarini <mzb.devel@gmail.com>
......
NEWS
+ 19
− 9
View file @ 269eae6c
......@@ -78,6 +78,11 @@ NEWS/ChangeLog for hydroPSO
-) boundary.wall: 'reflecting' in hydroPSO ver <= 0.1-58 is not longer equivalent to 'reflecting' in
hydroPSO ver >= 0.1-58
o 'hydromod' : -) now it is able to handle sub-daily models
-) if the 'hydroGOF' package is not available, a simple correlation plot is produced between observations and
the best simulation
-) improved documentation
o 'test_functions': -) new benchmark functions:
'schwefel' : Schwefel function
'ssphere' : shifted Sphere (CEC 2005),
......@@ -85,15 +90,20 @@ NEWS/ChangeLog for hydroPSO
'srastrigin' : shifted Rastrigin (CEC 2005),
'sgriewank' : shifted Griewank (CEC 2005),
'sackley' : shifted Ackley (CEC 2005),
'sschwefel1_2': shifted Schwefel's Problem 1.2 (CEC 2005)
o 'plot_out' : -) argument 'MinMax' is not required any more
-) argument 'sim' may also be 'integer' (before it had to be 'numeric' != 'integer')
o 'hydromod' : -) now it is able to handle sub-daily models
-) if the 'hydroGOF' package is not available, a simple correlation plot is produced between observations and
the best simulation
-) improved documentation
'sschwefel1_2': shifted Schwefel's Problem 1.2 (CEC 2005)
o 'read_particles': -) new argument 'do.pairs' (by default 'do.pairs=FALSE'), to control if the correlation matrix among parameters has to be plotted or not
(up to hydroPSO <=0.1.58 it was always plotted)
o 'plot_particles': -) new argument 'do.pairs' (by default 'do.pairs=FALSE'), to control if the correlation matrix among parameters has to be plotted or not
(up to hydroPSO <=0.1.58 it was always plotted)
o 'plot_results' : -) new argument 'do.pairs' (by default 'do.pairs=FALSE'), to control if the correlation matrix among parameters has to be plotted or not
(up to hydroPSO <=0.1.58 it was always plotted)
o 'plot_out' : -) argument 'MinMax' is not required any more!
-) argument 'sim' may also be 'integer' (before it had to be 'numeric', which does not allow 'integer' objects !)
o new function 'pest2hydroPSO' for importing PEST input files into hydroPSO
o new function 'hydroPSO2pest' for exporting hydroPSO input files to PEST
......
R/hydroPSO2pest.R
+ 10
− 10
View file @ 269eae6c
......@@ -415,13 +415,13 @@ hydroPSO2pest <- function(
setwd("/mnt/netapp1/nahaUsers/rojasro/test_functions/SWAT2005_hydroPSO")
hydroPSO2pest( paramfiles.fname="ParamFiles.txt",
paramranges.fname="ParamRanges.txt",
observations.fname="PEST2hydroPSO_OBS.txt",
exe.fname="./swat2005.out",
drty.model="/mnt/netapp1/nahaUsers/rojasro/test_functions/SWAT2005_hydroPSO",
#rscript.fname,
pst.fname="hydroPSO2PEST_test1.pst"
)
#setwd("/mnt/netapp1/nahaUsers/rojasro/test_functions/SWAT2005_hydroPSO")
#hydroPSO2pest( paramfiles.fname="ParamFiles.txt",
# paramranges.fname="ParamRanges.txt",
# observations.fname="PEST2hydroPSO_OBS.txt",
# exe.fname="./swat2005.out",
# drty.model="/mnt/netapp1/nahaUsers/rojasro/test_functions/SWAT2005_hydroPSO",
# #rscript.fname,
# pst.fname="hydroPSO2PEST_test1.pst"
#
# )
R/pest2hydroPSO.R
+ 39
− 16
View file @ 269eae6c
......@@ -59,14 +59,14 @@
################################################################################
## .pst2paramfiles ##
## .pst2paramfiles ##
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 08-Nov-2012 ##
# Updates: 09-Nov-2012 ##
################################################################################
# Purpose : To write the 'ParamFiles.txt' hydroPSO input file ##
# Purpose : To write the 'ParamFiles.txt' hydroPSO input file ##
################################################################################
.pst2paramfiles <- function(drty.model, tpls, inputs, param.names,
fname.out="ParamFiles.txt", DecimalPlaces=5) {
......@@ -165,7 +165,8 @@ pest2hydroPSO <- function(pst.fname,
drty.out="PSO.in",
paramfiles.fname="ParamFiles.txt",
paramranges.fname="ParamRanges.txt",
DecimalPlaces=5) {
DecimalPlaces=5,
verbose=TRUE) {
if (missing(pst.fname)) stop("PEST control file is missing ('pst.fname')")
if (is.null(drty.pest)) drty.pest <- dirname(pst.fname)
......@@ -181,27 +182,33 @@ pest2hydroPSO <- function(pst.fname,
paramranges.fname <- paste(drty.out, "/", paramranges.fname, sep="")
##############################################################################
# Reading .pst file
# 1) Reading .pst file
if (verbose) message(" ]")
if (verbose) message("[ 1) Reading the .pst file '", pst.fname, "' ... ]")
x <- readLines(pst.fname)
##############################################################################
# 1.a) Getting the number of parameters and observations
# Getting the number of parameters and observations
values <- as.numeric(strsplit(x[4], " ")[[1]])
nna.index <- which(!is.na(values))
nparam <- values[nna.index][1]
nobs <- values[nna.index][2]
if (verbose) message("[ Number of parameters found :", nparam, " ]")
if (verbose) message("[ Number of observations found:", nobs, " ]")
# 1.b) Getting the number of input files (.tpl) and output files (.ins)
# Getting the number of input files (.tpl) and output files (.ins)
files <- strsplit(x[5], " ")[[1]]
spc.index <- which(files=="")
files <- files[-spc.index]
ntpl <- as.numeric(files[1])
nins <- as.numeric(files[2])
if (verbose) message("[ Number of .tpl found :", ntpls, " ]")
if (verbose) message("[ Number of .ins found :", nins, " ]")
##############################################################################
# 2) Getting Param names and Ranges
# Getting Param names and Ranges
L <- 0
params.stg <- "* parameter data"
ini.index <- which(x==params.stg)
......@@ -217,14 +224,19 @@ pest2hydroPSO <- function(pst.fname,
param.max <- as.numeric(tmp[,4])
} else stop("Invalid pst file: ", params.stg, " does not exist !")
# writing 'ParamRanges.txt'
##############################################################################
# 2) Writing 'ParamRanges.txt'
if (verbose) message(" ]")
if (verbose) message("[ 2) Writing 'ParamRanges.txt' into '", drty.model, "' ... ]")
.pst2paramranges(drty.model=drty.model, names=param.names, ini=param.ini,
min=param.min, max=param.max,
fname.out=paramranges.fname)
##############################################################################
# 3) Getting observations
# 3) Writing file with observations
if (verbose) message(" ]")
if (verbose) message("[ 3) Writing 'PEST2hydroPSO_OBS.txt' into '", drty.out, "' ... ]")
L <- 0
obs.stg <- "* observation data"
ini.index <- which(x==obs.stg)
......@@ -244,6 +256,8 @@ pest2hydroPSO <- function(pst.fname,
##############################################################################
# 4) Model command line
if (verbose) message(" ]")
if (verbose) message("[ 4) Searching model command line ... ]")
L <- 0
model.stg <- "* model command line"
ini.index <- which(x==model.stg)
......@@ -252,10 +266,14 @@ pest2hydroPSO <- function(pst.fname,
if (L > 0) {
model.exe <- x[ini.index+1]
} else stop("Invalid pst file: ", model.stg, " does not exist !")
if (verbose) message("[ Model command line : '", model.exe, "' ... ]")
##############################################################################
# 5) Getting Param filenames and locations (.tpl's and .ins's)
# 5) Getting Parameter filenames and locations (.tpl's and .ins's)
if (verbose) message(" ]")
if (verbose) message("[ 5) Writing model input/output ... ]")
L <- 0
io.stg <- "* model input/output"
ini.index <- which(x==io.stg)
......@@ -278,6 +296,8 @@ pest2hydroPSO <- function(pst.fname,
##############################################################################
# 6) Creating the Rscript used for running hydroPSO
if (verbose) message(" ]")
if (verbose) message("[ 6) Copying R script for running hydroPSO ... ]")
rscript.fname <- system.file("Rscripts/hydroPSO-Rscript.R", package="hydroPSO")
dst.fname <- paste(drty.model, "/Rscripts/hydroPSO-Rscript.R", sep="")
if (!file.exists(dirname(dst.fname))) dir.create(dirname(dst.fname), recursive=TRUE)
......@@ -285,6 +305,8 @@ pest2hydroPSO <- function(pst.fname,
##############################################################################
# 7) Modifying the Rscript used for running hydroPSO
if (verbose) message(" ]")
if (verbose) message("[ 6) Modifying R script for running hydroPSO ... ]")
# rading the Rscript
x <- readLines(dst.fname)
......@@ -301,9 +323,10 @@ pest2hydroPSO <- function(pst.fname,
##############################################################################
# 8) Output
message("[ PEST2hydroPSO finished !! ]")
message("[R script to run hydroPSO available in: '", dst.fname, "']")
message("[Before running hydroPSO, you MUST modify the section: ")
message(" 'User-defined variables' in '", basename(dst.fname), "']")
if (verbose) message(" ]")
if (verbose) message("[ PEST2hydroPSO finished !! ]")
if (verbose) message("[ R script to run hydroPSO available in: '", dst.fname, "']")
if (verbose) message("[ Before running hydroPSO, you MUST modify the section: ")
if (verbose) message(" 'User-defined variables' in '", basename(dst.fname), "']")
} # 'pest2hydroPSO' END
R/plot_particles.R
+ 50
− 46
View file @ 269eae6c
......@@ -8,7 +8,7 @@
# 'plot_particles' #
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas #
# Started: 08-Nov-2011, #
# Updates: 02-Feb-2012 ; 17-Feb-2012 ; 09-Mar-2012 #
# Updates: 02-Feb-2012 ; 17-Feb-2012 ; 09-Mar-2012 ; 12-Nov-2012 #
################################################################################
# This function plots the contents of the 'Particles.txt' ouput file of #
# hydroPSO, with the position and fitness value of all the particles in the #
......@@ -16,9 +16,11 @@
# The following plots are produced: #
# 1) Dotty Plots of Parameter Values #
# 2) Histograms of Parameter Values #
# 3) Empirical CDFs of Parameter Values #
# 4) Parameter Values Against Number of Model Evaluations #
# 5) Plotting pseudo-3D dotty plots of Parameter Values #
# 3) Boxplots of Parameter Values #
# 4) Correlation matrix among Parameter Values (optional) #
# 5) Empirical CDFs of Parameter Values #
# 6) Parameter Values Against Number of Model Evaluations #
# 7) Plotting pseudo-3D dotty plots of Parameter Values #
################################################################################
plot_particles <- function(#####################################################
......@@ -42,8 +44,9 @@ plot_particles <- function(#####################################################
cex.lab=1.5,
#...,
breaks="Scott",
freq=TRUE,
freq=TRUE,
do.pairs=FALSE,
#####################################################
# For ECDFs of parameter values ('params2ecdf')
weights=NULL,
......@@ -176,7 +179,7 @@ plot_particles <- function(#####################################################
)
#############################################################################
# 3) Plotting Histograms of Parameter Values
# 3) Plotting boxplots of Parameter Values
msg <- "[ Plotting boxplots for parameter values"
if (do.png) msg <- paste(msg, " into '", basename(bxp.png.fname), sep="")
msg <- paste(msg, "' ... ]", sep="")
......@@ -218,46 +221,47 @@ plot_particles <- function(#####################################################
#############################################################################
# 4) Plotting Correlation Matrix of Parameter Values (with hydroTSM::hydropairs)
if ( require(hydroTSM) ) {
msg <- "[ Plotting correlation matrix for parameter values"
if (do.png) msg <- paste(msg, " into '", basename(pairs.png.fname), sep="")
msg <- paste(msg, "' ... ]", sep="")
if (verbose) message(msg)
if (!do.png) x11()
if (do.pairs)
if ( require(hydroTSM) ) {
msg <- "[ Plotting correlation matrix for parameter values"
if (do.png) msg <- paste(msg, " into '", basename(pairs.png.fname), sep="")
msg <- paste(msg, "' ... ]", sep="")
if (verbose) message(msg)
if (!do.png) x11()
plot_params(params=params,
gofs=gofs,
ptype="pairs",
param.cols=1:nparam,
param.names=param.names,
#of.col=ncol(z),
of.name=gof.name,
MinMax=MinMax,
beh.thr=beh.thr,
beh.col=beh.col,
beh.lty=beh.lty,
beh.lwd=beh.lwd,
nrows=nrows,
col=col,
ylab=ylab,
main=main,
pch=pch,
cex=cex,
cex.main=cex.main,
cex.axis=cex.axis,
cex.lab=cex.lab,
#...,
breaks=breaks,
freq=freq,
verbose=FALSE,
# PNG options
do.png=do.png,
png.width=png.width,
png.height=png.height,
png.res=png.res,
png.fname=pairs.png.fname
)
} # IF end
plot_params(params=params,
gofs=gofs,
ptype="pairs",
param.cols=1:nparam,
param.names=param.names,
#of.col=ncol(z),
of.name=gof.name,
MinMax=MinMax,
beh.thr=beh.thr,
beh.col=beh.col,
beh.lty=beh.lty,
beh.lwd=beh.lwd,
nrows=nrows,
col=col,
ylab=ylab,
main=main,
pch=pch,
cex=cex,
cex.main=cex.main,
cex.axis=cex.axis,
cex.lab=cex.lab,
#...,
breaks=breaks,
freq=freq,
verbose=FALSE,
# PNG options
do.png=do.png,
png.width=png.width,
png.height=png.height,
png.res=png.res,
png.fname=pairs.png.fname
)
} else warning("'hydroTSM' package is missing: Correlation among parameters was not plotted !")
#############################################################################
# 5) Empirical CDFs of Parameter Values
......
R/plot_results.R
+ 5
− 3
View file @ 269eae6c
......@@ -10,7 +10,7 @@
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas #
# Started: 10-Nov-2011, #
# Updates: 13-Ene-2012 ; 15-Feb-2012 ; 21-Feb-2012 ; 09-Mar-2012 ; 23-Mar-2012 #
# 11-Jun-2012 #
# 11-Jun-2012 ; 12-Nov-2012 #
################################################################################
plot_results <- function(drty.out="PSO.out",
......@@ -36,6 +36,7 @@ plot_results <- function(drty.out="PSO.out",
#...,
breaks="Scott",
freq=TRUE,
do.pairs=FALSE,
#######################################################
#### For ECDFs of parameter values ('params2ecdf') ####
......@@ -164,7 +165,7 @@ plot_results <- function(drty.out="PSO.out",
# Checking 'param.names'
for ( i in 1:npar) {
if ( !(param.names[i] %in% colnames(params)) )
stop("Invalid argument: The field '", param.names[i], "' doesn't exist in 'params'")
stop("Invalid argument: The field '", param.names[i], "' does not exist in 'params' !")
} # IF end
# Subsetting
......@@ -192,7 +193,7 @@ plot_results <- function(drty.out="PSO.out",
# 1) Dotty Plots,
# 2) Histograms,
# 3) Boxplots
# 4) Correlation Matrix
# 4) Correlation Matrix (optional)
# 5) Empirical CDFs
# 6) Parameter Values Against Number of Model Evaluations
# 7) (pseudo)3D dotty plots
......@@ -217,6 +218,7 @@ plot_results <- function(drty.out="PSO.out",
#...,
breaks=breaks,
freq=freq,
do.pairs=do.pairs,
#####################################################
# For ECDFs of parameter values ('params2ecdf')
......
R/read_particles.R
+ 3
− 1
View file @ 269eae6c
......@@ -9,7 +9,7 @@
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas #
# Started: 08-Nov-2011, #
# Updates: 17-Feb-2012 ; 09-Mar-2012 #
# Updates: 17-Feb-2012 ; 09-Mar-2012 ; 12-Nov-2012 #
################################################################################
# This function reads the 'Particles.txt' ouput file of hydroPSO and returns #
# a data.frame with the position and fitness value of all the particles in the #
......@@ -41,6 +41,7 @@ read_particles <- function(file="Particles.txt",
#...,
breaks="Scott",
freq=TRUE,
do.pairs=FALSE,
# Parameters for the 3D dotty plots
dp3D.names="auto",
GOFcuts="auto",
......@@ -165,6 +166,7 @@ read_particles <- function(file="Particles.txt",
#...,
breaks=breaks,
freq=freq,
do.pairs=do.pairs,
# Parameters for the 3D dotty plots
dp3D.names=dp3D.names,
GOFcuts=GOFcuts,
......
man/ReadPlot_particles.Rd
+ 8
− 3
View file @ 269eae6c
......@@ -20,7 +20,7 @@ The function \code{plot_particles} takes the parameter sets and their correspond
1) Dotty plots \cr
2) Histograms \cr
3) Boxplots \cr
4) Correlation matrix \cr
4) Correlation matrix (optional) \cr
5) Empirical CDFs \cr
6) Parameter values vs Number of Model Evaluations \cr
7) (pseudo) 3D dotty plots
......@@ -31,7 +31,8 @@ read_particles(file="Particles.txt", verbose=TRUE, plot=TRUE,
beh.lwd=2, nrows="auto", col="black", ylab=gof.name, main=NULL, pch=19,
cex=0.5, cex.main=1.5, cex.axis=1.5, cex.lab=1.5,
%%...,
breaks="Scott", freq=TRUE, dp3D.names="auto", GOFcuts="auto",
breaks="Scott", freq=TRUE, do.pairs=FALSE,
dp3D.names="auto", GOFcuts="auto",
colorRamp= colorRampPalette(c("darkred", "red", "orange", "yellow",
"green", "darkgreen", "cyan")), alpha=1, points.cex=0.7,legend.pos="topleft",
do.png=FALSE, png.width=1500, png.height=900, png.res=90,
......@@ -48,7 +49,7 @@ plot_particles(params, gofs, gof.name="GoF", MinMax=NULL, beh.thr=NA,
ylab=gof.name, main=NULL, pch=19, cex=0.5, cex.main=1.5,
cex.axis=1.5, cex.lab=1.5,
%%...,
breaks="Scott", freq=TRUE,
breaks="Scott", freq=TRUE, do.pairs=FALSE,
weights=NULL, byrow=FALSE, leg.cex=1.5,
dp3D.names="auto", GOFcuts="auto",
colorRamp= colorRampPalette(c("darkred", "red", "orange", "yellow",
......@@ -153,6 +154,10 @@ breaks for plotting the histograms of the parameter sets. See \code{\link[graphi
\item{freq}{
OPTIONAL. Only used when \code{plot=TRUE} \cr
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 \code{\link[graphics]{hist}}
}
\item{do.pairs}{
OPTIONAL. Only used when \code{plot=TRUE} \cr
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.
}
\item{weights}{
OPTIONAL. Only used when \code{plot=TRUE} \cr
......
man/ReadPlot_results.Rd
+ 7
− 3
View file @ 269eae6c
......@@ -27,7 +27,7 @@ The function \code{plot_results} takes the outputs of the \code{read_results} fu
1) Dotty plots of parameter values \cr
2) Histograms of parameter values \cr
3) Boxplots of parameter values \cr
4) Correlation matrix among parameter values \cr
4) Correlation matrix among parameter values (optional) \cr
5) Empirical CDFs of parameter values\cr
6) Parameter values vs Number of Model Evaluations \cr
7) (pseudo) 3D dotty plots of (selected) parameter values \cr
......@@ -45,7 +45,7 @@ read_results(drty.out="PSO.out", MinMax=NULL, beh.thr=NA, modelout.cols=NULL, ve
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,
cex.axis=1.3, cex.lab=1.5, breaks="Scott", freq=TRUE, do.pairs=FALSE,
weights=NULL, byrow=FALSE, leg.cex=1.2,
%%...
%% Parameters for the 3D dotty plots
......@@ -160,6 +160,10 @@ breaks for plotting the histograms of the parameter sets. See \code{\link[graphi
\item{freq}{
OPTIONAL. Only used when \code{plot=TRUE} \cr
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 \code{\link[graphics]{hist}}
}
\item{do.pairs}{
OPTIONAL. Only used when \code{plot=TRUE} \cr
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.
}
\item{weights}{
OPTIONAL. Only used when \code{plot=TRUE} \cr
......@@ -286,7 +290,7 @@ The function \code{read_results} returns a list with the following elements:
\item{model.best}{numeric with the best model / objective function value. In order to be computed, the user has to provide a valid value for \code{MinMax}. See \code{\link{read_out}} }
\item{model.obs}{numeric with the observed values used during the optimisation. See \code{obs}}
\item{convergence.measures}{matrix/data.frame with the convergence measures. See \code{\link{read_convergence}} function}
\item{part.GofPerIter}{matrix/data.frame with the goodness-of-fit only for all the particles during all the iterations}
\item{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}
}
%%\references{
%%
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment