diff --git a/DESCRIPTION b/DESCRIPTION
index fb198c9a02e42255e115d6b592fc2bda3c436955..b2ca4d5251c9cfe1c56eb76867ec33e6f4711c59 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-21
-Date: 2012-11-09
+Version: 0.1-58-22
+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>
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/NAMESPACE b/NAMESPACE
index 60531c8d77ad910c1ba25e9df8feb26e9f48778b..a384594234741b9b04f5118536dc0cb131dc37a4 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -56,7 +56,9 @@ export(hydroPSO,
srastrigin,
srosenbrock,
sschwefel1_2,
- ssphere
+ ssphere,
+ hydroPSO2pest,
+ pest2hydroPSO
)
S3method(params2ecdf, matrix)
diff --git a/NEWS b/NEWS
index 80e339106a7eabf1a959152e0a6152529075a151..7454d8625c5bb6396284278997ca23bfcd83db88 100755
--- a/NEWS
+++ b/NEWS
@@ -18,8 +18,6 @@ NEWS/ChangeLog for hydroPSO
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
- -) new function 'pest2hydroPSO' for importing PEST input files
-
-) 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
@@ -95,7 +93,10 @@ NEWS/ChangeLog for hydroPSO
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
+ -) improved documentation
+
+ o new function 'pest2hydroPSO' for importing PEST input files into hydroPSO
+ o new function 'hydroPSO2pest' for exporting hydroPSO input files to PEST
0.1-58 14-Sep-2012
diff --git a/R/hydroPSO2pest.R b/R/hydroPSO2pest.R
new file mode 100644
index 0000000000000000000000000000000000000000..8bcfa44df5619a69d326ddfef60d694c410f31bb
--- /dev/null
+++ b/R/hydroPSO2pest.R
@@ -0,0 +1,427 @@
+# File hydroPSO2pest.R
+# Part of the hydroPSO package, http://www.rforge.net/hydroPSO/
+# Copyright 2012-2012 Mauricio Zambrano-Bigiarini
+# Distributed under GPL 2 or later
+
+################################################################################
+## .read_paramrange ##
+################################################################################
+# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
+################################################################################
+# Created: 12-Nov-2012 ##
+# Updates: 12-Nov-2012 ##
+################################################################################
+# 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)
+
+ # Reading the file
+ x <- read.table(fname, header=TRUE, stringsAsFactors=FALSE)
+
+ return(x)
+
+} # '.read_paramranges' END
+
+
+
+################################################################################
+## .read_paramfiles ##
+################################################################################
+# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
+################################################################################
+# Created: 12-Nov-2012 ##
+# Updates: 12-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") {
+
+ x <- read.table(fname, header=TRUE, stringsAsFactors=FALSE)
+
+} # '.read_paramfiles' END
+
+
+################################################################################
+## .read_observations ##
+################################################################################
+# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
+################################################################################
+# Created: 12-Nov-2012 ##
+# Updates: 12-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") {
+
+ # reading the first line of the file
+ x <- read.table(fname, header=FALSE, nrows=1, stringsAsFactors=FALSE)
+
+ if (NCOL(x) == 1) {
+ x <- read.table(fname, header=FALSE)
+ } else x <- coredata(read.zoo(fname, header=FALSE))
+
+ return(as.numeric(x[,1]))
+
+} # '.read_observations' END
+
+
+################################################################################
+## .write.tpl ##
+################################################################################
+# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
+################################################################################
+# Created: 12-Nov-2012 ##
+# Updates: 12-Nov-2012 ##
+################################################################################
+# Purpose : To write a .tpl PEST input file ##
+################################################################################
+.write.tpl <- function(tpl.fname, modelin.fname, paramfiles, token="#") {
+
+ # Getting the line for the input file
+ row.index <- which(paramfiles[,3] == basename(modelin.fname))
+
+ # Reading the model input file
+ x <- readLines(tpl.fname)
+
+ # Writing the tpl file
+ for (r in row.index ) {
+ # parameter name
+ ParameterName <- as.character(paramfiles[r, 2])
+ # Row location
+ RowNumber <- paramfiles[r, 4]
+ # Columns location
+ ColStart <- paramfiles[r, 5]
+ ColEnd <- paramfiles[r, 6]
+
+ L <- ColEnd - ColStart + 1
+ nspaces <- L - 2 - nchar(ParameterName)
+ spaces <- paste(rep(" ", nspaces), collapse="")
+ new.stg <- paste(token, ParameterName, spaces, token, sep="")
+
+ # Replacing the line with the token and parameter name
+ substr(x[RowNumber], start=ColStart, stop=ColEnd) <- new.stg
+ } # FOR end
+
+ # writing the .tpl to disk
+ writeLines(x, con=tpl.fname)
+
+} # '.write.tpl' END
+
+
+
+################################################################################
+## hydroPSO2pest ##
+################################################################################
+# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
+################################################################################
+# Created: 12-Nov-2012 ##
+# Updates: 12-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 ##
+################################################################################
+hydroPSO2pest <- function(
+ paramfiles.fname="ParamFiles.txt",
+ paramranges.fname="ParamRanges.txt",
+ observations.fname="Observations.txt",
+ exe.fname,
+ drty.model=getwd(),
+ #rscript.fname,
+ pst.fname="hydroPSO2PEST.pst"
+
+ ) {
+
+ if (missing(paramfiles.fname))
+ stop("Missing argument: 'paramfiles.fname'")
+
+ if (missing(paramranges.fname))
+ stop("Missing argument: 'paramranges.fname'")
+
+ if (missing(exe.fname))
+ stop("Missing argument: 'exe.fname'")
+
+ if (missing(observations.fname))
+ stop("Missing argument: 'observations.fname'")
+
+ if (!file.exists(drty.model)) stop("Invalid argument: 'drty.model' does not exist !")
+
+ ##############################################################################
+
+ if (basename(paramfiles.fname)==paramfiles.fname)
+ paramfiles.fname <- paste(drty.model, "/PSO.in/", paramfiles.fname, sep="")
+
+ if (basename(paramranges.fname)==paramranges.fname)
+ paramranges.fname <- paste(drty.model, "/PSO.in/", paramranges.fname, sep="")
+
+ if (basename(observations.fname)==observations.fname)
+ observations.fname <- paste(drty.model, "/PSO.out/", observations.fname, sep="")
+
+ 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), "' !")
+
+ ##############################################################################
+
+ # 1) Reading ParamRanges.txt
+ x <- .read_paramranges(drty.model, fname=paramranges.fname)
+
+ # Param numbers
+ param.nmbrs <- x[, 1]
+ # Param names
+ param.names <- x[, 2]
+ # Param boundaries
+ param.min <- x[, 3]
+ param.max <- x[, 4]
+
+ # Number of parameters
+ nparam <- length(param.names)
+
+ ##############################################################################
+
+ # 2) Reading Observations.txt
+ obs <- .read_observations(fname=observations.fname)
+
+ # Number of observations
+ nobs <- length(obs)
+
+ ##############################################################################
+
+ # 3) Reading ParamFiles.txt
+ x <- .read_paramfiles(fname=paramfiles.fname)
+
+ # Param numbers
+ ParameterNmbr <- x[, 1]
+ # Param names
+ ParameterName <- x[, 2]
+ # Model input names
+ Filename <- as.factor(x[, 3])
+ # Row location
+ RowNumber <- x[, 4]
+ # Columns location
+ ColStart <- x[, 5]
+ ColEnd <- x[, 6]
+ # Number of decimal places
+ DecimalPlaces <- x[, 7]
+
+ # Name of single-occurrence input files to be modified
+ tpl.fnames <- as.character(levels(Filename))
+
+ # Number of tpl files
+ ntpls <- length(tpl.fnames)
+
+ ##############################################################################
+ # 4) tpl creation
+
+ 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)
+
+ # Reading the new .tpl
+ tmp <- readLines(tpl.fname)
+
+ # Adding a 1st line with the token
+ token <- "#"
+ first.line <- paste("ptf ", token, sep="")
+ tmp <- c(first.line, tmp)
+
+ # writing the .tpl to disk
+ writeLines(tmp, con=tpl.fname)
+
+ .write.tpl(tpl.fname=tpl.fname, modelin.fname=fname.in, paramfiles=x, token=token)
+ } # FOR 't' end
+
+
+ ##############################################################################
+ # 5) ins creation
+
+ x <- c(NA, NA, NA, NA, NA)
+ x[1] <- "***THIS FILE IS TO BE CREATED BY THE USER BEFORE RUNNING PEST. SEE TEMPLATE BELOW***"
+ x[2] <- "pif @ "
+ x[3] <- "@***OUTPUT KEYWORD***@"
+ x[4] <- "l1 (obs1)StartCol:EndCol"
+ x[5] <- "..."
+
+ # writing to disk
+ fname.out <- paste(drty.model, "/ModelOut.ins", sep="")
+ writeLines(x, con=fname.out)
+
+
+ ##############################################################################
+ # 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)
+ file.copy(pst.template, pst.fname, overwrite=TRUE)
+
+ # modifying the .pst
+ x <- readLines(pst.fname)
+ L <- length(x)
+
+ # 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
+ 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
+ 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
+ tmp <- as.character(ntpls)
+ l <- nchar(tmp)
+ if (l < 4) tmp <- paste(paste(rep(" ", 4-l), collapse=""), ntpls, sep="")
+ substr(x[5], start=2, stop=5) <- tmp
+
+ # Number of .ins
+ #substr(x[5], start=7, stop=12) <- "1"
+ # it should be changed by the user
+
+
+
+ # * parameter groups
+ x <- c(x[1:11], rep(x[12], nparam), x[13:L])
+ L <- L + nparam - 1
+ for (i in 1:nparam) {
+ tmp <- param.names[i]
+ l <- nchar(tmp)
+ if (l < 17) tmp <- paste(paste(rep(" ", 17-l), collapse=""), param.names[i], sep="")
+ substr(x[11+i], start=1, stop=17) <- tmp
+ } # FOR end
+
+
+ # * parameter data
+ 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) {
+ # param names
+ tmp <- param.names[i]
+ l <- nchar(tmp)
+ if (l < 17) tmp <- paste(paste(rep(" ", 17-l), collapse=""), param.names[i], sep="")
+ substr(x[13+nparam-1+i], start=1, stop=17) <- tmp
+
+ # param ini
+ tmp <- format((param.max[i] + param.min[i])/2, scientific=TRUE, digits=7)
+ tmp <- as.character(tmp)
+ l <- nchar(tmp)
+ if (l < 13) tmp <- paste(paste(rep(" ", 13-l), collapse=""), tmp, sep="")
+ substr(x[13+nparam-1+i], start=33, stop=45) <- tmp
+
+ # param min
+ tmp <- format(param.min[i], scientific=TRUE, digits=7)
+ tmp <- as.character(tmp)
+ l <- nchar(tmp)
+ if (l < 13) tmp <- paste(paste(rep(" ", 13-l), collapse=""), tmp, sep="")
+ substr(x[13+nparam-1+i], start=47, stop=59) <- tmp
+
+ # param max
+ tmp <- format(param.max[i], scientific=TRUE, digits=7)
+ tmp <- as.character(tmp)
+ l <- nchar(tmp)
+ if (l < 15) tmp <- paste(paste(rep(" ", 15-l), collapse=""), tmp, sep="")
+ substr(x[13+nparam-1+i], start=61, stop=76) <- tmp
+
+ } # 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])
+
+
+ L <- L + nobs - 1
+ 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="")
+ substr(x[17+2*(nparam-1)+i], start=1, stop=8) <- tmp
+
+ # i-th observation value
+ tmp <- format(obs[i], scientific=TRUE, digits=18)
+ 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
+
+ } # FOR end
+
+
+
+ # * model command line
+ x[20+2*(nparam-1)+(nobs-1)] <- exe.fname
+
+ # * model input/output
+ 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) {
+
+ # i-th .tpl filename
+ l <- nchar(i)
+ if (l < 3) tmp <- paste("file", paste(rep("0", 3-l), collapse=""), i, sep="")
+ substr(x[21+2*(nparam-1)+(nobs-1)+i], start=1, stop=11) <- tmp
+
+ # i-th model input file
+ tmp <- as.character(tpl.fnames[i])
+ 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
+
+ } # FOR end
+
+ # writing to disk
+ writeLines(x, con=pst.fname)
+
+
+ ##############################################################################
+ # 8) Output
+ message("[ ]")
+ message("[ hydroPSO2PEST finished !! ]")
+ message("[ PEST input files available in: '", drty.model, "']")
+ 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 7a7f55785ed0ad698e46d49f3ea084a1ffa17ae4..f16a4e7fe18ba1f36746add0735d34994d57ab28 100644
--- a/R/pest2hydroPSO.R
+++ b/R/pest2hydroPSO.R
@@ -1,4 +1,34 @@
-.pst2paramranges <- function(drty.model, names, ini, min, max, fname.out="ParamRanges.txt") {
+# File pest2hydroPSO.R
+# Part of the hydroPSO package, http://www.rforge.net/hydroPSO/
+# Copyright 2012-2012 Mauricio Zambrano-Bigiarini
+# Distributed under GPL 2 or later
+
+################################################################################
+## .pst2paramranges ##
+################################################################################
+# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
+################################################################################
+# Created: 08-Nov-2012 ##
+# Updates: 09-Nov-2012 ##
+################################################################################
+# Purpose : To write the 'ParamRanges.txt' hydroPSO input file ##
+################################################################################
+
+# 'drty.model': character, with the name of the directory with the input and exe
+# files of the model
+# 'names' : character vector, with the names of the parameters to be calibrated
+# 'ini' : numeric, with the initial guess used in PEST for each parameter
+# to be calibrated (not used in hydroPSO)
+# 'min' : numeric, with the lowest boundaries for the the parameters to be
+# calibrated
+# 'max' : numeric, with the highest boundaries for the the parameters to be
+# calibrated
+# 'fname.out' : character, with the name of the output text file with the boudaries
+# for each parameter to be calibrated
+
+
+.pst2paramranges <- function(drty.model, names, ini, min, max,
+ fname.out="ParamRanges.txt") {
drty.bak <- getwd()
setwd(drty.model)
@@ -28,7 +58,18 @@
-.pst2paramfiles <- function(drty.model, tpls, inputs, param.names, fname.out="ParamFiles.txt", DecimalPlaces=5) {
+################################################################################
+## .pst2paramfiles ##
+################################################################################
+# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
+################################################################################
+# Created: 08-Nov-2012 ##
+# Updates: 09-Nov-2012 ##
+################################################################################
+# Purpose : To write the 'ParamFiles.txt' hydroPSO input file ##
+################################################################################
+.pst2paramfiles <- function(drty.model, tpls, inputs, param.names,
+ fname.out="ParamFiles.txt", DecimalPlaces=5) {
drty.bak <- getwd()
setwd(drty.model)
@@ -107,9 +148,17 @@
} # '.pst2paramfiles' END
-
-
-
+################################################################################
+## pest2hydroPSO ##
+################################################################################
+# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
+################################################################################
+# Created: 08-Nov-2012 ##
+# Updates: 09-Nov-2012 ##
+################################################################################
+# Purpose : To import the PEST input files (.pst, .tpl) to be used within ##
+# hydroPSO (ParamFiles.txt, ParamRanges.txt, hydroPSO_Rscript.R) ##
+################################################################################
pest2hydroPSO <- function(pst.fname,
drty.pest=NULL,
drty.model=NULL,
diff --git a/inst/hydroPSO2pest.pst b/inst/hydroPSO2pest.pst
new file mode 100644
index 0000000000000000000000000000000000000000..20cfbe07d7cfd50f289ea5e4ff8b23dafeaf59ce
--- /dev/null
+++ b/inst/hydroPSO2pest.pst
@@ -0,0 +1,24 @@
+pcf
+* control data
+restart estimation
+ npar nobs npar 0 1
+ ntpl nins single point 1 0 0
+ 5.0 2.0 0.3 0.03 10 999
+ 3.0 3.0 0.001 0
+ 0.1
+ 30 0.01 3 3 0.01 3
+ 1 1 1
+* parameter groups
+param.name relative 0.01 0.0 switch 2.0 parabolic
+* parameter data
+param.name none relative param.ini param.min param.max param.name 1.0000 0.0000 1
+* observation groups
+obsgroup
+* observation data
+obs00001 obs.value 1.0 obsgroup
+* model command line
+exe.fname
+* model input/output
+file001.tpl input001.txt
+ModelOut.ins output.txt
+* prior information