diff --git a/DESCRIPTION b/DESCRIPTION
index 9fe0c773c9e3a45cf1c57a9847f5e228a0a07299..85ae92d43d3ea77ba6541300778e5170458ac2ce 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: hydroPSO
Type: Package
Title: Model-Independent Particle Swarm Optimisation for Environmental Models
-Version: 0.1-58-23
-Date: 2012-11-12
+Version: 0.1-58-24
+Date: 2012-11-15
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>
Description: This package implements a state-of-the-art version of the Particle Swarm Optimisation (PSO) algorithm (SPSO-2011 and SPSO-2007 capable), with a special focus on the calibration of environmental models. hydroPSO is model-independent, allowing the user to easily interface any model code with the calibration engine (PSO). It includes a series of controlling options and PSO variants to fine-tune the performance of the calibration engine to different calibration problems. An advanced sensitivity analysis function together with user-friendly plotting summaries facilitate the interpretation and assessment of the calibration results. Bugs reports/comments/questions are very welcomed.
diff --git a/NEWS b/NEWS
index 36d2a4e8b1b968817fe985cdd86f7b0ead3d535c..b8b0e091f10cec84767c02fffe7ef5185210843a 100755
--- a/NEWS
+++ b/NEWS
@@ -1,109 +1,111 @@
NEWS/ChangeLog for hydroPSO
--------------------------
0.1-59 (under-development & not working yet)
- o 'hydroPSO' : -) towards SPSO 2011 capable, which includes (check the latest \url{clerc.maurice.free.fr/pso/SPSO_descriptions.pdf}):
- * new default value for swarm size (npart=40),
- * new equation for updating the velocity ( V[t+1] = w*V + Gr -x + r*alea_sphere(Gr, ||Gr-x||) ),
- * new initialization of the velocity ( Vini = U(lower-Xini, upper-Xini) )
- * new confinement of the velocity ( V[t+1] = -0.5 * V[t], when x[t+1] > x_max | x[t+1] < x_min )
- * optional normalisation of parameter values
+ o 'hydroPSO' : -) towards SPSO 2011 capable, which includes (check the latest \url{clerc.maurice.free.fr/pso/SPSO_descriptions.pdf}):
+ * new default value for swarm size (npart=40),
+ * new equation for updating the velocity ( V[t+1] = w*V + Gr -x + r*alea_sphere(Gr, ||Gr-x||) ),
+ * new initialization of the velocity ( Vini = U(lower-Xini, upper-Xini) )
+ * new confinement of the velocity ( V[t+1] = -0.5 * V[t], when x[t+1] > x_max | x[t+1] < x_min )
+ * optional normalisation of parameter values
- -) improved performance: ~ 33% faster for 32-bit machines and 38% faster for 64-bit machines.
- Tested on 10-, 20- and 30-D benchmark functions with 'write2disk=FALSE'.
- -) 'fn' argument now can be any R function or a character. In the latter case, it can be "hydromod" or the name of a
- valid R function. In previous versions of 'hydroPSO' only a character type was accepted.
- -) now it handles models with sub-daily time step and with sub-daily observations (thanks to O. Rakovec !).
- Suggested dependence on hydroGOF >= 0.3-5 and hydroTSM >= 0.3-6
- -) new output file 'BestModel_out.txt', with the model outputs corresponding to the best parameter set.
- Only available when 'fn.name=="hydromod"'
- -) new 'normalise' parameter for the 'control' variable, in order to improve the performance when the search space is not
- a hypercube
- -) normalised swarm radius now it is computed using the median distance of all the particles the global best, instead of using the maximum distance
- (as proposed in Evers and Ghalia 2009). This was done in order to make easier the identification of the stagnation point for activating regrouping
- -) when 'use.RG=TRUE', the default values for 'RG.thr', 'RG.r', and 'RG.miniter' were changed (related to the change in the computation of swarm
- radius), from 1.1E-4, 0.8 and 5 to 1E-5, 2 and 100, respectively
- -) argument 'method' now allows the following 3 new values: 'spso2007', 'spso2011' and 'canonical', which set the values of the PSO
- engine to the ones of the SPSO 2007, SPSO 2011 and the canonical one, respectively.
- -) argument 'method'. The old (default) value 'pso' was replaced by the new (default) value 'spso2011'
- -) argument 'Vini.type' now allows 2 new values for setting the initial velocity of the particles according to the equation
- specified in the SPSO 2011: 'random2011' and 'lhs2011', which use a uniform distribution or with a Latin-hypercube
- sampling, respectively.
- -) argument 'Vini.type'. The old values 'random' and 'lhs' were replaced by 'random2007' and 'lhs2007', in order to make
- clear that it follows the equation described in SPSO 2007.
- -) argument 'Vini.type', when missing its value depends on the value of the 'method' argument:
- method == 'spso2007' => Vini.type='random2007'
- method != 'spso2007' => Vini.type='random2011'
- -) argument 'Xini.type', now its default value is 'random' instead of 'lhs', in order to make it compatible with the default
- value for 'method' (method='spso2011')
- -) argument 'npart', when missing its value depends on the value of the 'method' argument:
- method == 'spso2007' => npart=10+2*[sqrt(n)]
- method != 'spso2007' => npart=40
- -) argument 'boundary.wall' now allows the following new values: 'absorbing2011', which set 'boundary.wall' to the
- absorbing condition specified in the SPSO 2011
- -) argument 'boundary.wall'. The old value 'absorbing' was replaced by 'absorbing2007', in order to make clear that is
- follows the equation described in SPSO 2007.
- -) argument 'boundary.wall', when missing its value depends on the value of the 'method' argument:
- method == 'spso2007' => boundary.wall='absorbing2007'
- method != 'spso2007' => boundary.wall='absorbing2011'
- -) in the documentation, default values are now mentioned for each argument of 'control'.
- -) the 'PSO_logfile.txt' output file now contains information about: 'best.update', 'random.update', 'Xini.type', 'Vini.type'.
- In addition, information about 'IW.type' and 'IW.exp' is now written only when 'length(IW.w) > 1'
- -) output file 'BestParamPerParticle.txt' now has the GoF as first column (previously it was the last one)
- -) fixed velocity equation for 'boundary.wall="reflecting"' :
- ver <= 0.1-58 -> ver >= 0.1-59
- v[t+1] : v[t] -> -v[t] (when x[t+1] > x_max | x[t+1] < x_min )
- -) changed velocity equation for 'boundary.wall="damping"' :
- ver <= 0.1-58 -> ver >= 0.1-59
- v[t+1] : v[t] -> -v[t] (when x[t+1] > x_max | x[t+1] < x_min )
- -) when the control argument 'out.with.fit.iter' (not used so far) is set to TRUE, the number of iterations returned
- correspond to the effective number of iterations carried out, not the maximum defined by 'maxit'
- -) observed values are now correctly written to disk for sub-daily models
- -) the package now depends on R >= 2.13.0, not 2.10.0 as in previous releases, in order to be consistent with the
- byte-compiler setting of the package
- -) fixed some (very unlikely) error when 'IW.type="aiwf"'
+ -) improved performance: ~ 33% faster for 32-bit machines and 38% faster for 64-bit machines.
+ Tested on 10-, 20- and 30-D benchmark functions with 'write2disk=FALSE'.
+ -) 'fn' argument now can be any R function or a character. In the latter case, it can be "hydromod" or the name of a
+ valid R function. In previous versions of 'hydroPSO' only a character type was accepted.
+ -) now it handles models with sub-daily time step and with sub-daily observations (thanks to O. Rakovec !).
+ Suggested dependence on hydroGOF >= 0.3-5 and hydroTSM >= 0.3-6
+ -) new output file 'BestModel_out.txt', with the model outputs corresponding to the best parameter set.
+ Only available when 'fn.name=="hydromod"'
+ -) new 'normalise' parameter for the 'control' variable, in order to improve the performance when the search space is not
+ a hypercube
+ -) normalised swarm radius now it is computed using the median distance of all the particles the global best, instead of
+ using the maximum distance (as proposed in Evers and Ghalia 2009). This was done in order to make easier the
+ identification of the stagnation point for activating regrouping
+ -) when 'use.RG=TRUE', the default values for 'RG.thr', 'RG.r', and 'RG.miniter' were changed (related to the change in the
+ computation of swarm radius), from 1.1E-4, 0.8 and 5 to 1E-5, 2 and 100, respectively
+ -) argument 'method' now allows the following 3 new values: 'spso2007', 'spso2011' and 'canonical', which set the values of the PSO
+ engine to the ones of the SPSO 2007, SPSO 2011 and the canonical one, respectively.
+ -) argument 'method'. The old (default) value 'pso' was replaced by the new (default) value 'spso2011'
+ -) argument 'Vini.type' now allows 2 new values for setting the initial velocity of the particles according to the equation
+ specified in the SPSO 2011: 'random2011' and 'lhs2011', which use a uniform distribution or with a Latin-hypercube
+ sampling, respectively.
+ -) argument 'Vini.type'. The old values 'random' and 'lhs' were replaced by 'random2007' and 'lhs2007', in order to make
+ clear that it follows the equation described in SPSO 2007.
+ -) argument 'Vini.type', when missing its value depends on the value of the 'method' argument:
+ method == 'spso2007' => Vini.type='random2007'
+ method != 'spso2007' => Vini.type='random2011'
+ -) argument 'Xini.type', now its default value is 'random' instead of 'lhs', in order to make it compatible with the default
+ value for 'method' (method='spso2011')
+ -) argument 'npart', when missing its value depends on the value of the 'method' argument:
+ method == 'spso2007' => npart=10+2*[sqrt(n)]
+ method != 'spso2007' => npart=40
+ -) argument 'boundary.wall' now allows the following new values: 'absorbing2011', which set 'boundary.wall' to the
+ absorbing condition specified in the SPSO 2011
+ -) argument 'boundary.wall'. The old value 'absorbing' was replaced by 'absorbing2007', in order to make clear that is
+ follows the equation described in SPSO 2007.
+ -) argument 'boundary.wall', when missing its value depends on the value of the 'method' argument:
+ method == 'spso2007' => boundary.wall='absorbing2007'
+ method != 'spso2007' => boundary.wall='absorbing2011'
+ -) in the documentation, default values are now mentioned for each argument of 'control'.
+ -) the 'PSO_logfile.txt' output file now contains information about: 'best.update', 'random.update', 'Xini.type', 'Vini.type'.
+ In addition, information about 'IW.type' and 'IW.exp' is now written only when 'length(IW.w) > 1'
+ -) output file 'BestParamPerParticle.txt' now has the GoF as first column (previously it was the last one)
+ -) fixed velocity equation for 'boundary.wall="reflecting"' :
+ ver <= 0.1-58 -> ver >= 0.1-59
+ v[t+1] : v[t] -> -v[t] (when x[t+1] > x_max | x[t+1] < x_min )
+ -) changed velocity equation for 'boundary.wall="damping"' :
+ ver <= 0.1-58 -> ver >= 0.1-59
+ v[t+1] : v[t] -> -v[t] (when x[t+1] > x_max | x[t+1] < x_min )
+ -) when the control argument 'out.with.fit.iter' (not used so far) is set to TRUE, the number of iterations returned
+ correspond to the effective number of iterations carried out, not the maximum defined by 'maxit'
+ -) observed values are now correctly written to disk for sub-daily models
+ -) the package now depends on R >= 2.13.0, not 2.10.0 as in previous releases, in order to be consistent with the
+ byte-compiler setting of the package
+ -) fixed some (very unlikely) error when 'IW.type="aiwf"'
o Running 'hydroPSO' >= 0.1-59 with default settings will produce DIFFERENT RESULTS from those obtained with 'hydroPSO' <= 0.1-58, due to the
following changes in default values:
- ver <= 0.1-58 -> ver >= 0.1-59
- -) npart : 10+2*[sqrt(n)] -> 40
- -) Xini.type : 'lhs' -> 'random'
- -) Vini.type : 'lhs2007' -> 'random2011'
- -) boundary.wall: 'reflecting' -> 'absorbing2011'
+ ver <= 0.1-58 -> ver >= 0.1-59
+ -) npart : 10+2*[sqrt(n)] -> 40
+ -) Xini.type : 'lhs' -> 'random'
+ -) Vini.type : 'lhs2007' -> 'random2011'
+ -) boundary.wall: 'reflecting' -> 'absorbing2011'
- -) TVc1.type : 'non-linear' -> 'linear' (no effect, because 'use.TVc1=FALSE' by default)
- -) TVc2.type : 'non-linear' -> 'linear' (no effect, because 'use.TVc2=FALSE' by default)
- -) TVlambda.type: 'non-linear' -> 'linear' (no effect, because 'use.TVlambda=FALSE' by default)
+ -) TVc1.type : 'non-linear' -> 'linear' (no effect, because 'use.TVc1=FALSE' by default)
+ -) TVc2.type : 'non-linear' -> 'linear' (no effect, because 'use.TVc2=FALSE' by default)
+ -) TVlambda.type: 'non-linear' -> 'linear' (no effect, because 'use.TVlambda=FALSE' by default)
- -) boundary.wall: 'reflecting' in hydroPSO ver <= 0.1-58 is not longer equivalent to 'reflecting' in
- hydroPSO ver >= 0.1-58
+ -) 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 '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),
- 'srosenbrock' : shifted Rosenbrock (CEC 2005),
- '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 'test_functions' : -) new benchmark functions:
+ 'schwefel' : Schwefel function
+ 'ssphere' : shifted Sphere (CEC 2005),
+ 'srosenbrock' : shifted Rosenbrock (CEC 2005),
+ '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_convergence': -) the label 'Gbest' was replaced by "Global Optimum' in the title of the plot and in the label of the 'y' axis, in
+ order to make it more intuitive to people non-familiar with PSO
+ 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 '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_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_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 '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
diff --git a/R/hydroPSO2pest.R b/R/hydroPSO2pest.R
index 571cda4324e645b74be6f424ecd12c42e2f78e13..1e51ed2641bd3445c6990c7f507f9814d4d7ae21 100644
--- a/R/hydroPSO2pest.R
+++ b/R/hydroPSO2pest.R
@@ -9,17 +9,14 @@
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
-# Updates: 12-Nov-2012 ##
+# Updates: 15-Nov-2012 ##
################################################################################
-# Purpose : To read the 'ParamRanges.txt' hydroPSO input file ##
+# Purpose : To read the 'ParamRanges.txt' hydroPSO input file ##
################################################################################
-
-# 'drty.model': character, with the name of the directory with the input and exe
-# files of the model
-
-.read_paramranges <- function(drty.model, fname="ParamRanges.txt") {
-
- setwd(drty.model)
+.read_paramranges <- function(fname="ParamRanges.txt") {
+
+ if (!file.exists(fname))
+ stop("Invalid argument: file '", fname, "' does not exists !")
# Reading the file
x <- read.table(fname, header=TRUE, stringsAsFactors=FALSE)
@@ -36,15 +33,16 @@
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
-# Updates: 12-Nov-2012 ##
+# Updates: 15-Nov-2012 ##
################################################################################
# Purpose : To read the 'ParamFiles.txt' hydroPSO input file ##
################################################################################
-
-# 'drty.model': character, with the name of the directory with the input and exe
-# files of the model
.read_paramfiles <- function(fname="ParamFiles.txt") {
+ if (!file.exists(fname))
+ stop("Invalid argument: file '", fname, "' does not exists !")
+
+ # Reading the file
x <- read.table(fname, header=TRUE, stringsAsFactors=FALSE)
} # '.read_paramfiles' END
@@ -56,15 +54,15 @@
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
-# Updates: 12-Nov-2012 ##
+# Updates: 15-Nov-2012 ##
################################################################################
# Purpose : To read the 'Observations.txt' hydroPSO output file ##
################################################################################
-
-# 'drty.model': character, with the name of the directory with the input and exe
-# files of the model
.read_observations <- function(fname="Observations.txt") {
+ if (!file.exists(fname))
+ stop("Invalid argument: file '", fname, "' does not exists !")
+
# reading the first line of the file
x <- read.table(fname, header=FALSE, nrows=1, stringsAsFactors=FALSE)
@@ -78,16 +76,23 @@
################################################################################
-## .write.tpl ##
+## .modify.tpl ##
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
-# Updates: 12-Nov-2012 ##
+# Updates: 15-Nov-2012 ##
################################################################################
-# Purpose : To write a .tpl PEST input file ##
+# Purpose : To write a .tpl PEST input file based on the 'modelin.fname' ##
+# model input file ##
+# It identifies all the parameters that have to be changed within ##
+# the 'modelin.fname' file, and then it replaces by the name of ##
+# each parameter enclosed between user-defined tokens ##
################################################################################
-.write.tpl <- function(tpl.fname, modelin.fname, paramfiles, token="#") {
+.modify.tpl <- function(tpl.fname, modelin.fname, paramfiles, token="#", verbose=TRUE) {
+
+ if (!file.exists(tpl.fname))
+ stop("Invalid argument: file '", tpl.fname, "' does not exists !")
# Getting the line for the input file
row.index <- which(paramfiles[,3] == basename(modelin.fname))
@@ -116,8 +121,10 @@
# writing the .tpl to disk
writeLines(x, con=tpl.fname)
+
+ if (verbose) message("[ File '", basename(tpl.fname), "' successfully modified ]")
-} # '.write.tpl' END
+} # '.modify.tpl' END
@@ -127,28 +134,56 @@
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
-# Updates: 12-Nov-2012 ##
+# Updates: 15-Nov-2012 ##
################################################################################
# Purpose : To export the content of the hydroPSO files 'ParamRanges.txt' ##
# 'ParamFiles.txt' to PEST, as a single .pst files with its corres-##
# ponding .tpl and .ins files ##
################################################################################
+# 'paramfile.fname' : character, file name (full path) of the hydroPSO ##
+# input file storing the location and names of the model##
+# files that have to be modified for each parameter ##
+# subject to calibration. ##
+# By default this file is called 'ParamFiles.txt' and ##
+# -if the full path it is not specified- it is searched ##
+# for within the 'PSO.in' subdirectory of 'drty.model' ##
+# 'param.ranges' : character with the name of the hydroPSO input file ##
+# defining the minimum and maximum boundary values for ##
+# each one of the parameters to be calibrated. ##
+# By default this file is called 'ParamRanges.txt' and ##
+# -if the full path it is not specified- it is searched ##
+# for within the 'PSO.in' subdirectory of 'drty.model' ##
+# 'observations.fname': character with the name of the hydroPSO output file ##
+# storing the observed values used during the optimisa- ##
+# tion. ##
+# By default this file is called 'Observations.txt' and ##
+# -if the full path it is not specified- it is searched ##
+# for within the 'PSO.out' ##
+# subdirectory of 'drty.model'. ##
+# 'exe.fname' : character, model command line arguments to be entered ##
+# through a prompted string to execute the user-defined ##
+# model ##
+# 'drty.model' : character, path to the executable file of the model ##
+# specified in \code{exe.fname}. ALL the input files ##
+# required to run the model have to be located within ##
+# this directory ##
+# 'pst.fname' : character, with the name of the output .pst file ##
+################################################################################
hydroPSO2pest <- function(
- paramfiles.fname="ParamFiles.txt",
- paramranges.fname="ParamRanges.txt",
+ param.files="ParamFiles.txt",
+ param.ranges="ParamRanges.txt",
observations.fname="Observations.txt",
exe.fname,
drty.model=getwd(),
- #rscript.fname,
- pst.fname="hydroPSO2PEST.pst"
-
+ pst.fname="hydroPSO2PEST.pst",
+ verbose=TRUE
) {
- if (missing(paramfiles.fname))
- stop("Missing argument: 'paramfiles.fname'")
+ if (missing(param.files))
+ stop("Missing argument: 'param.files'")
- if (missing(paramranges.fname))
- stop("Missing argument: 'paramranges.fname'")
+ if (missing(param.ranges))
+ stop("Missing argument: 'param.ranges'")
if (missing(exe.fname))
stop("Missing argument: 'exe.fname'")
@@ -160,11 +195,11 @@ hydroPSO2pest <- function(
##############################################################################
- if (basename(paramfiles.fname)==paramfiles.fname)
- paramfiles.fname <- paste(drty.model, "/PSO.in/", paramfiles.fname, sep="")
+ if (basename(param.files)==param.files)
+ param.files <- paste(drty.model, "/PSO.in/", param.files, sep="")
- if (basename(paramranges.fname)==paramranges.fname)
- paramranges.fname <- paste(drty.model, "/PSO.in/", paramranges.fname, sep="")
+ if (basename(param.ranges)==param.ranges)
+ param.ranges <- paste(drty.model, "/PSO.in/", param.ranges, sep="")
if (basename(observations.fname)==observations.fname)
observations.fname <- paste(drty.model, "/PSO.out/", observations.fname, sep="")
@@ -172,14 +207,19 @@ hydroPSO2pest <- function(
if (basename(pst.fname)==pst.fname)
pst.fname <- paste(drty.model, "/", pst.fname, sep="")
- if (!file.exists(paramfiles.fname)) stop("Invalid argument: 'paramfiles.fname' does not exist !")
- if (!file.exists(paramranges.fname)) stop("Invalid argument: 'paramranges.fname' does not exist !")
- if (!file.exists(observations.fname)) stop("Invalid argument: 'observations.fname' does not exist in '", dirname(observations.fname), "' !")
+ if (!file.exists(param.files))
+ stop("Invalid argument: The file '", param.files, "' does not exist !")
+ if (!file.exists(param.ranges))
+ stop("Invalid argument: The file '", param.ranges, "' does not exist !")
+ if (!file.exists(observations.fname))
+ stop("Invalid argument: The file '", observations.fname, "' does not exist !")
##############################################################################
# 1) Reading ParamRanges.txt
- x <- .read_paramranges(drty.model, fname=paramranges.fname)
+ if (verbose) message(" ]")
+ if (verbose) message("[ 1) Reading '", basename(param.ranges), "' ... ]")
+ x <- .read_paramranges(fname=param.ranges)
# Param numbers
param.nmbrs <- x[, 1]
@@ -191,19 +231,25 @@ hydroPSO2pest <- function(
# Number of parameters
nparam <- length(param.names)
+ if (verbose) message("[ Number of parameters found :", nparam, " ]")
##############################################################################
# 2) Reading Observations.txt
+ if (verbose) message(" ]")
+ if (verbose) message("[ 2) Reading observations ('", basename(param.ranges), "') ... ]")
obs <- .read_observations(fname=observations.fname)
# Number of observations
nobs <- length(obs)
+ if (verbose) message("[ Number of observations found:", nobs, " ]")
##############################################################################
# 3) Reading ParamFiles.txt
- x <- .read_paramfiles(fname=paramfiles.fname)
+ if (verbose) message(" ]")
+ if (verbose) message("[ 3) Reading '", basename(param.files), "' ... ]")
+ x <- .read_paramfiles(fname=param.files)
# Param numbers
ParameterNmbr <- x[, 1]
@@ -224,18 +270,22 @@ hydroPSO2pest <- function(
# Number of tpl files
ntpls <- length(tpl.fnames)
+ if (verbose) message("[ Number of .tpl to be created:", ntpls, " ]")
##############################################################################
# 4) tpl creation
-
- for (t in 1:ntpls) {
-
+ if (verbose) message(" ]")
+ if (verbose) message("[ 4) Creating .tpl files ... ]")
+
+ for (t in 1:ntpls) {
tmp <- t
l <- nchar(t)
if (l < 3) tmp <- paste(paste(rep("0", 3-l), collapse=""), t, sep="")
tpl.fname <- paste(drty.model, "/", "file", tmp, ".tpl", sep="")
- fname.in <- paste(drty.model, "/", tpl.fnames[t], sep="")
- file.copy(fname.in, tpl.fname, overwrite=TRUE)
+ fname.in <- paste(drty.model, "/", tpl.fnames[t], sep="")
+ if (!file.exists(fname.in)) {
+ stop("Invalid argument: The file '", fname.in, "' does not exist !")
+ } else file.copy(fname.in, tpl.fname, overwrite=TRUE)
# Reading the new .tpl
tmp <- readLines(tpl.fname)
@@ -245,15 +295,19 @@ hydroPSO2pest <- function(
first.line <- paste("ptf ", token, sep="")
tmp <- c(first.line, tmp)
- # writing the .tpl to disk
+ # writing the .tpl to disk, with exactly the same content as the original input file
writeLines(tmp, con=tpl.fname)
- .write.tpl(tpl.fname=tpl.fname, modelin.fname=fname.in, paramfiles=x, token=token)
+ # writing parameter names enclosed by tokens within the .tpl
+ .modify.tpl(tpl.fname=tpl.fname, modelin.fname=fname.in, paramfiles=x,
+ token=token, verbose=verbose)
} # FOR 't' end
##############################################################################
# 5) ins creation
+ if (verbose) message(" ]")
+ if (verbose) message("[ 5) Creating a basic .ins file ... ]")
x <- c(NA, NA, NA, NA, NA)
x[1] <- "***THIS FILE IS TO BE CREATED BY THE USER BEFORE RUNNING PEST. SEE TEMPLATE BELOW***"
@@ -269,35 +323,49 @@ hydroPSO2pest <- function(
##############################################################################
# 6) pst creation
-
- ##pst.template <- system.file("hydroPSO2pest.pst", package="hydroPSO")
- pst.template <- "/dataMZB/2012/mzb_functions/R-mzb_packages/SVN/hydroPSO/trunk/inst/hydroPSO2pest.pst"
- if (!file.exists(dirname(pst.fname))) dir.create(dirname(pst.fname), recursive=TRUE)
+ if (verbose) message(" ]")
+ if (verbose) message("[ 6) Creating the .tpl file from template ... ]")
+
+ # Copying the .pst template provided with the package
+ pst.template <- system.file("hydroPSO2pest.pst", package="hydroPSO")
+ #pst.template <- "/dataMZB/2012/mzb_functions/R-mzb_packages/SVN/hydroPSO/trunk/inst/hydroPSO2pest.pst"
+ if (!file.exists(dirname(pst.template)))
+ stop("Invalid argument: The file '", pst.template, "' does not exist !")
+ if (!file.exists(dirname(pst.fname)))
+ dir.create(dirname(pst.fname), recursive=TRUE)
file.copy(pst.template, pst.fname, overwrite=TRUE)
- # modifying the .pst
+ ##############################################################################
+ # 7) modifying the .pst
+ if (verbose) message(" ]")
+ if (verbose) message("[ 7) Modifying the .tpl file ... ]")
+
x <- readLines(pst.fname)
L <- length(x)
- # Number of parameters
+ # 7.1) Number of parameters
+ if (verbose) message("[ 7.1) Writing number of parameters ... ]")
tmp <- as.character(nparam)
l <- nchar(tmp)
if (l < 4) tmp <- paste(paste(rep(" ", 4-l), collapse=""), nparam, sep="")
substr(x[4], start=2, stop=5) <- tmp
- # Number of observations
+ # 7.2) Number of observations
+ if (verbose) message("[ 7.2) Writing number of observations ... ]")
tmp <- as.character(nobs)
l <- nchar(tmp)
if (l < 6) tmp <- paste(paste(rep(" ", 6-l), collapse=""), nobs, sep="")
substr(x[4], start=7, stop=12) <- tmp
- # Number of parameter groups
+ # 7.3) Number of parameter groups
+ if (verbose) message("[ 7.3) Writing number of parametr groups ... ]")
tmp <- as.character(nparam)
l <- nchar(tmp)
if (l < 5) tmp <- paste(paste(rep(" ", 5-l), collapse=""), nparam, sep="")
substr(x[4], start=14, stop=18) <- tmp
- # Number of .tpl
+ # 7.4) Number of .tpl
+ if (verbose) message("[ 7.4) Writing number of .tpl files ... ]")
tmp <- as.character(ntpls)
l <- nchar(tmp)
if (l < 4) tmp <- paste(paste(rep(" ", 4-l), collapse=""), ntpls, sep="")
@@ -305,11 +373,10 @@ hydroPSO2pest <- function(
# Number of .ins
#substr(x[5], start=7, stop=12) <- "1"
- # it should be changed by the user
-
-
+ # it should be changed by the user
- # * parameter groups
+ # 7.5) * parameter groups
+ if (verbose) message("[ 7.5) Writing '* parameter groups' section ... ]")
x <- c(x[1:11], rep(x[12], nparam), x[13:L])
L <- L + nparam - 1
for (i in 1:nparam) {
@@ -320,7 +387,8 @@ hydroPSO2pest <- function(
} # FOR end
- # * parameter data
+ # 7.6) * parameter data
+ if (verbose) message("[ 7.6) Writing '* parameter data' section ... ]")
x <- c(x[1:(13+nparam-1)], rep(x[14+nparam-1], nparam), x[(15+nparam-1):L])
L <- L + nparam - 1
for (i in 1:nparam) {
@@ -354,13 +422,12 @@ hydroPSO2pest <- function(
} # FOR end
- # * observation data
- x <- c(x[1:(17+2*(nparam-1))], rep(x[17+2*(nparam-1)+1], nobs), x[(17+2*(nparam-1)+2):L])
-
+ # 7.7) * observation data
+ if (verbose) message("[ 7.7) Writing '* observation data' section ... ]")
+ x <- c(x[1:(17+2*(nparam-1))], rep(x[17+2*(nparam-1)+1], nobs), x[(17+2*(nparam-1)+2):L])
L <- L + nobs - 1
- for (i in 1:nobs) {
-
+ for (i in 1:nobs) {
# i-th observation name
l <- nchar(i)
if (l < 5) tmp <- paste("obs", paste(rep("0", 5-l), collapse=""), i, sep="")
@@ -371,21 +438,20 @@ hydroPSO2pest <- function(
tmp <- as.character(tmp)
l <- nchar(tmp)
if (l < 23) tmp <- paste(paste(rep(" ", 23-l), collapse=""), tmp, sep="")
- substr(x[17+2*(nparam-1)+i], start=14, stop=36) <- tmp
-
+ substr(x[17+2*(nparam-1)+i], start=14, stop=36) <- tmp
} # FOR end
-
- # * model command line
+ # 7.8) * model command line
+ if (verbose) message("[ 7.8) Writing '* model command line' section ... ]")
x[20+2*(nparam-1)+(nobs-1)] <- exe.fname
- # * model input/output
+ # 7.9) * model input/output
+ if (verbose) message("[ 798) Writing '* model input/output' section ... ]")
x <- c(x[1:(21+2*(nparam-1)+(nobs-1))], rep(x[21+2*(nparam-1)+(nobs-1)+1], ntpls), x[(21+2*(nparam-1)+(nobs-1)+2):L])
L <- L + ntpls - 1
- for (i in 1:ntpls) {
-
+ for (i in 1:ntpls) {
# i-th .tpl filename
l <- nchar(i)
if (l < 3) tmp <- paste("file", paste(rep("0", 3-l), collapse=""), i, sep="")
@@ -396,8 +462,7 @@ hydroPSO2pest <- function(
l <- nchar(tmp)
if (l < 12) tmp <- paste(paste(rep(" ", 12-l), collapse=""), tpl.fnames[i], sep="")
- substr(x[21+2*(nparam-1)+(nobs-1)+i], start=13, stop=13+l+1) <- tmp
-
+ substr(x[21+2*(nparam-1)+(nobs-1)+i], start=13, stop=13+l+1) <- tmp
} # FOR end
# writing to disk
@@ -407,21 +472,9 @@ hydroPSO2pest <- function(
##############################################################################
# 8) Output
message("[ ]")
- message("[ hydroPSO2PEST finished !! ]")
+ message("[ hydroPSO2PEST finished !! ]")
message("[ PEST input files available in: '", drty.model, "']")
+ message("[ ]")
message("[ Before running PEST, you MUST check *.pst and *.ins files ]")
} # 'pest2hydroPSO' END
-
-
-
-#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"
-#
-# )
diff --git a/R/pest2hydroPSO.R b/R/pest2hydroPSO.R
index 7e3601f7e9335e26974db9c4bfe07e748b516df2..8bb860aaaf5a169342777cc566267d2c88ad8dc9 100644
--- a/R/pest2hydroPSO.R
+++ b/R/pest2hydroPSO.R
@@ -69,7 +69,7 @@
# Purpose : To write the 'ParamFiles.txt' hydroPSO input file ##
################################################################################
.pst2paramfiles <- function(drty.model, tpls, inputs, param.names,
- fname.out="ParamFiles.txt", DecimalPlaces=5) {
+ fname.out="ParamFiles.txt", decimals=5) {
drty.bak <- getwd()
setwd(drty.model)
@@ -82,7 +82,7 @@
stop("Invalid argument: dimensions do not match !")
field.names <- c("ParameterNmbr", "ParameterName", "Filename",
- "Row.Number", "Col.Start", "Col.End", "DecimalPlaces")
+ "Row.Number", "Col.Start", "Col.End", "decimals")
#tmp <- matrix(NA, ncol=length(filed.names), nrow=nparam)
@@ -113,10 +113,10 @@
token.pos <- which(strsplit(x[l], '', useBytes=TRUE)[[1]]==token)
- out[[p]] <- rbind(out[[p]], c(p, param.names[p], inputs[f], l, token.pos[1], token.pos[2], DecimalPlaces) )
+ out[[p]] <- rbind(out[[p]], c(p, param.names[p], inputs[f], l, token.pos[1], token.pos[2], decimals) )
if (length(token.pos) >2) {
for ( t in seq(3, length(token.pos), by=2) )
- out[[p]] <- rbind(out[[p]], c(p, param.names[p], inputs[f], l, token.pos[t], token.pos[t+1], DecimalPlaces) )
+ out[[p]] <- rbind(out[[p]], c(p, param.names[p], inputs[f], l, token.pos[t], token.pos[t+1], decimals) )
} # IF end
} # IF end
@@ -135,7 +135,7 @@
tmp <- tmp[-1, ]
# identifying possible errors
- error.index <- which ( (as.numeric(tmp[,6]) - as.numeric(tmp[,5]) + 1 ) <= DecimalPlaces)
+ error.index <- which ( (as.numeric(tmp[,6]) - as.numeric(tmp[,5]) + 1 ) <= decimals)
if (length(error.index) > 0) {
warning("In ParamFiles.txt, decimal places have to be manually corrected:", paste(error.index) )
} # IF
@@ -154,18 +154,57 @@
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 08-Nov-2012 ##
-# Updates: 09-Nov-2012 ##
+# Updates: 09-Nov-2012 ; 15-Nov-2012 ##
################################################################################
# Purpose : To import the PEST input files (.pst, .tpl) to be used within ##
# hydroPSO (ParamFiles.txt, ParamRanges.txt, hydroPSO_Rscript.R) ##
################################################################################
+# 'pst.fname' : character, with name of the PEST input file (.pst), ##
+# which contains all the information regarding ##
+# parameters, observations and template files (.tpl and ##
+# .ins) used by PEST ##
+# 'drty.pest' : character, path to the executable file of PEST. ALL ##
+# the input files required to run PEST with the model ##
+# have to be located within this directory (.tpl and ##
+# .ins) ##
+# 'drty.model' : character, path to the executable file of the model ##
+# specified in \code{exe.fname}. ALL the input files ##
+# required to run the model have to be located within ##
+# this directory ##
+# 'drty.out' : character, name of the directory that will store all ##
+# the output files produced by this function ##
+# 'param.files' : character, file name (full path) of the hydroPSO ##
+# input file storing the location and names of the model##
+# files that have to be modified for each parameter ##
+# subject to calibration. ##
+# By default this file is called 'ParamFiles.txt' and ##
+# -if the full path it is not specified- it is searched ##
+# for within the 'PSO.in' subdirectory of 'drty.model' ##
+# 'param.ranges' : character with the name of the hydroPSO input file ##
+# defining the minimum and maximum boundary values for ##
+# each one of the parameters to be calibrated. ##
+# By default this file is called 'ParamRanges.txt' and ##
+# -if the full path it is not specified- it is searched ##
+# for within the 'PSO.in' subdirectory of 'drty.model' ##
+# 'observations.fname': character with the name of the hydroPSO output file ##
+# storing the observed values used during the optimisa- ##
+# tion. ##
+# By default this file is called 'Observations.txt' and ##
+# -if the full path it is not specified- it is searched ##
+# for within the 'PSO.out' ##
+# subdirectory of 'drty.model'. ##
+# 'exe.fname' : character, model command line arguments to be entered ##
+# through a prompted string to execute the user-defined ##
+# model ##
+# 'decimals' : ##
+################################################################################
pest2hydroPSO <- function(pst.fname,
drty.pest=NULL,
drty.model=NULL,
drty.out="PSO.in",
- paramfiles.fname="ParamFiles.txt",
- paramranges.fname="ParamRanges.txt",
- DecimalPlaces=5,
+ param.files="ParamFiles.txt",
+ param.ranges="ParamRanges.txt",
+ decimals=5,
verbose=TRUE) {
if (missing(pst.fname)) stop("PEST control file is missing ('pst.fname')")
@@ -175,11 +214,11 @@ pest2hydroPSO <- function(pst.fname,
if (drty.out=="PSO.in") drty.out <- paste(dirname(pst.fname), "/PSO.in", sep="")
if (!file.exists(drty.out)) dir.create(drty.out, recursive=TRUE)
- if (basename(paramfiles.fname)==paramfiles.fname)
- paramfiles.fname <- paste(drty.out, "/", paramfiles.fname, sep="")
+ if (basename(param.files)==param.files)
+ param.files <- paste(drty.out, "/", param.files, sep="")
- if (basename(paramranges.fname)==paramranges.fname)
- paramranges.fname <- paste(drty.out, "/", paramranges.fname, sep="")
+ if (basename(param.ranges)==param.ranges)
+ param.ranges <- paste(drty.out, "/", param.ranges, sep="")
##############################################################################
# 1) Reading .pst file
@@ -230,7 +269,7 @@ pest2hydroPSO <- function(pst.fname,
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)
+ fname.out=param.ranges)
##############################################################################
@@ -290,8 +329,8 @@ pest2hydroPSO <- function(pst.fname,
# Writing ParamFiles.txt
.pst2paramfiles(drty.model=drty.model, tpls=tpls, inputs=inputs,
- param.names=param.names, fname.out=paramfiles.fname,
- DecimalPlaces=DecimalPlaces)
+ param.names=param.names, fname.out=param.files,
+ decimals=decimals)
##############################################################################
diff --git a/R/plot_convergence.R b/R/plot_convergence.R
index 55e33c17e463234a4da0f0591d4dc7a4a42c4a0b..1d68d4bb819f58fc62c0cb68d50a249e7019b61b 100755
--- a/R/plot_convergence.R
+++ b/R/plot_convergence.R
@@ -8,7 +8,7 @@
# 'plot_convergence' #
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas #
# Started: 08-Nov-2011, #
-# Updates: 13-Ene-2012 #
+# Updates: 13-Ene-2012 ; 14-Nov-2012 #
################################################################################
@@ -17,9 +17,9 @@ plot_convergence <- function(x,
col=c("black", "darkolivegreen"),
lty=c(1,3),
lwd=c(2,2),
- main="Gbest & Normalized Swarm Radius vs Iteration Number",
+ main="Global Optimum & Normalized Swarm Radius vs Iteration Number",
xlab="Iteration Number",
- ylab=c("Gbest", expression(delta[norm]) ),
+ ylab=c("Global Optimum", expression(delta[norm]) ),
pch=c(15, 18),
cex=1,
cex.main=1.4, cex.axis=1.2, cex.lab=1.2,
diff --git a/R/read_convergence.R b/R/read_convergence.R
index 9cddfa3125f427d34b9341cbeaa26ffc7ca71a67..e57885a0119a474fa26e89cd56ccc30a9e91cf2b 100755
--- a/R/read_convergence.R
+++ b/R/read_convergence.R
@@ -9,7 +9,7 @@
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas #
# Started: 08-Nov-2011, #
-# Updates: 13-Ene-2012 ; 29-Ene-2012 ; 01-Feb-2012 #
+# Updates: 13-Ene-2012 ; 29-Ene-2012 ; 01-Feb-2012 ; 14-Nov-2012 #
################################################################################
# Purpose: To read the output file 'ConvergenceMeasures.txt' #
################################################################################
@@ -37,9 +37,9 @@ read_convergence <- function(file="ConvergenceMeasures.txt",
col=c("black", "darkolivegreen"),
lty=c(1,3),
lwd=c(2,2),
- main="Gbest & Normalized Swarm Radius vs Iteration Number",
+ main="Global Optimum & Normalized Swarm Radius vs Iteration Number",
xlab="Iteration Number",
- ylab=c("Gbest", expression(delta[norm]) ),
+ ylab=c("Global Optimum", expression(delta[norm]) ),
pch=c(15, 18),
cex=1,
cex.main=1.4,
diff --git a/man/ReadPlot_convergence.Rd b/man/ReadPlot_convergence.Rd
index 6200b1d2b175e02cafeb18a27b340fe13ddf9804..8a1663d9fa64b060f907ab3a8c32bcc712df52b8 100755
--- a/man/ReadPlot_convergence.Rd
+++ b/man/ReadPlot_convergence.Rd
@@ -18,15 +18,15 @@ This function reads a file containing different parameter sets and ther correspo
\usage{
read_convergence(file="ConvergenceMeasures.txt", MinMax=NULL, beh.thr=NA,
verbose=TRUE, plot=TRUE, col=c("black", "darkolivegreen"), lty=c(1,3),
- lwd=c(2,2), main="Gbest & Normalized Swarm Radius vs Iteration Number",
- xlab="Iteration Number", ylab=c("Gbest", expression(delta[norm])),
+ lwd=c(2,2), main="Global Optimum & Normalized Swarm Radius vs Iteration Number",
+ xlab="Iteration Number", ylab=c("Global Optimum", expression(delta[norm])),
pch=c(15, 18), cex=1, cex.main=1.4, cex.axis=1.2, cex.lab=1.2,
legend.pos="topright", ..., do.png=FALSE, png.width=1500, png.height=900,
png.res=90,png.fname="ConvergenceMeasures.png")
plot_convergence(x, verbose=TRUE, col=c("black", "darkolivegreen"), lty=c(1,3),
- lwd=c(2,2), main="Gbest & Normalized Swarm Radius vs Iteration Number",
- xlab="Iteration Number", ylab=c("Gbest", expression(delta[norm])),
+ lwd=c(2,2), main="Global Optimum & Normalized Swarm Radius vs Iteration Number",
+ xlab="Iteration Number", ylab=c("Global Optimum", expression(delta[norm])),
pch=c(15, 18), cex=1, cex.main=1.4, cex.axis=1.2, cex.lab=1.2,
legend.pos="topright", ..., do.png=FALSE, png.width=1500, png.height=900,
png.res=90, png.fname="ConvergenceMeasures.png")
diff --git a/man/hydroPSO.Rd b/man/hydroPSO.Rd
index 5575fea45742f08c7c39ff02593e4aa10f022c9e..c369f6f546b96eb0d342b32d53f4aa2d960b6229 100755
--- a/man/hydroPSO.Rd
+++ b/man/hydroPSO.Rd
@@ -87,7 +87,7 @@ character, path to the directory storing the output files generated by hydroPSO
}
\item{param.ranges}{
OPTIONAL. Used only when \code{fn='hydromod'} \cr
-character, name of the file storing the desired range of variation of each parameter
+character, name of the file defining the minimum and maximum boundary values for each one of the parameters to be calibrated
}
\item{digits}{
OPTIONAL. Used only when \code{write2disk=TRUE} \cr
diff --git a/man/hydromod.Rd b/man/hydromod.Rd
index 80ae4a623937c54fb43c2e670c8fa473cde32029..fd08addef9613d21e5bd9af4051c656695feae63 100755
--- a/man/hydromod.Rd
+++ b/man/hydromod.Rd
@@ -34,7 +34,8 @@ character, file name (full path) storing location and names of the files that ha
character, path to the executable file of the model specified in \code{exe.fname}. ALL the input files required to run the model have to be located within this directory
}
\item{exe.fname}{
-character, name and extension of the file used to run the model specified in \code{exe.fname}
+character, model command line arguments to be entered through a prompted string to execute the user-defined model
+
}
\item{stdout, stderr}{
where output to \sQuote{stdout} or \sQuote{stderr} should be sent. Possible values are \code{FALSE} (discard output, the default), \code{""}, to the R console. See \code{\link[base]{system2}}