diff --git a/DESCRIPTION b/DESCRIPTION
index 357734401960645faa31ac2c46baeb694cf2fd2a..fb198c9a02e42255e115d6b592fc2bda3c436955 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-20
-Date: 2012-11-08
+Version: 0.1-58-21
+Date: 2012-11-09
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 fa9a593f831766e5e73091a12684e30797e7850e..60531c8d77ad910c1ba25e9df8feb26e9f48778b 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,7 @@
importFrom("Hmisc", wtd.quantile, wtd.Ecdf)
importFrom("sp", coordinates, spplot)
importFrom("lattice", xyplot, axis.default)
-importFrom("zoo", is.zoo)
+importFrom("zoo", is.zoo, read.zoo, coredata)
import(grid)
#importFrom("hydroTSM", hydroplot, hydropairs, dip, mip, yip, vector2zoo)
#importFrom("hydroGOF", ggof)
diff --git a/NEWS b/NEWS
index eac9a5dfb75ec225c8f269734f3fd496a97518dd..80e339106a7eabf1a959152e0a6152529075a151 100755
--- a/NEWS
+++ b/NEWS
@@ -18,6 +18,8 @@ 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
diff --git a/R/pest2hydroPSO.R b/R/pest2hydroPSO.R
new file mode 100644
index 0000000000000000000000000000000000000000..7a7f55785ed0ad698e46d49f3ea084a1ffa17ae4
--- /dev/null
+++ b/R/pest2hydroPSO.R
@@ -0,0 +1,260 @@
+.pst2paramranges <- function(drty.model, names, ini, min, max, fname.out="ParamRanges.txt") {
+
+ drty.bak <- getwd()
+ setwd(drty.model)
+
+ nparam <- length(names)
+
+ if ( (nparam!=length(ini)) | (nparam!=length(min)) | (nparam!=length(max)) )
+ stop("Invalid argument: dimensions do not match !")
+
+ field.names <- c("ParameterNmbr", "ParameterName", "MinValue", "MaxValue", "IniPEST")
+
+ tmp <- matrix(NA, ncol=5, nrow=nparam)
+ tmp <- as.data.frame(tmp)
+ tmp[,1] <- 1:nparam
+ tmp[,2] <- names
+ tmp[,3] <- min
+ tmp[,4] <- max
+ tmp[,5] <- ini
+
+ colnames(tmp) <- field.names
+
+ write.table(tmp, file=fname.out, row.names=FALSE, quote=FALSE)
+
+ setwd(drty.bak)
+
+} # '.pst2paramranges' END
+
+
+
+.pst2paramfiles <- function(drty.model, tpls, inputs, param.names, fname.out="ParamFiles.txt", DecimalPlaces=5) {
+
+ drty.bak <- getwd()
+ setwd(drty.model)
+
+ ntpl <- length(tpls)
+ ninputs <- length(inputs)
+ nparam <- length(param.names)
+
+ if ( (ntpl!=ninputs) )
+ stop("Invalid argument: dimensions do not match !")
+
+ field.names <- c("ParameterNmbr", "ParameterName", "Filename",
+ "Row.Number", "Col.Start", "Col.End", "DecimalPlaces")
+
+ #tmp <- matrix(NA, ncol=length(filed.names), nrow=nparam)
+
+ # output creation
+ tmp <- matrix(NA, ncol=length(field.names), nrow=1)
+ tmp <- as.data.frame(tmp)
+ colnames(tmp) <- field.names
+
+ out <- vector("list", nparam)
+ for (i in 1:nparam) out[[i]] <- tmp
+
+
+ for (p in 1:nparam) {
+
+ for (f in 1:ntpl) {
+
+ if (!file.exists(tpls[f])) stop("Invalid argument: file '", tpls[f], " does not exist !!")
+ x <- readLines(tpls[f])
+
+ # getting the token
+ token <- strsplit(x[1], " ", useBytes=TRUE)[[1]][2]
+
+ for (l in 2:length(x)) {
+
+ exists <- grep(param.names[p], x[l], useBytes=TRUE)
+
+ if (length(exists) > 0) {
+
+ 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) )
+ 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) )
+ } # IF end
+
+ } # IF end
+
+ } # FOR l end
+
+ } # FOR f end
+
+ # removing dummy 1st row
+ out[[p]] <- out[[p]][-1,]
+ } # FOR 'p' end
+
+ for (p in 1:nparam) tmp <- rbind(tmp, out[[p]])
+
+ # removing dummy 1st row
+ tmp <- tmp[-1, ]
+
+ # identifying possible errors
+ error.index <- which ( (as.numeric(tmp[,6]) - as.numeric(tmp[,5]) + 1 ) <= DecimalPlaces)
+ if (length(error.index) > 0) {
+ warning("In ParamFiles.txt, decimal places have to be manually corrected:", paste(error.index) )
+ } # IF
+
+ # Writing the output file
+ write.table(tmp, file=fname.out, row.names=FALSE, quote=FALSE)
+
+ setwd(drty.bak)
+
+} # '.pst2paramfiles' END
+
+
+
+
+
+pest2hydroPSO <- function(pst.fname,
+ drty.pest=NULL,
+ drty.model=NULL,
+ drty.out="PSO.in",
+ paramfiles.fname="ParamFiles.txt",
+ paramranges.fname="ParamRanges.txt",
+ DecimalPlaces=5) {
+
+ if (missing(pst.fname)) stop("PEST control file is missing ('pst.fname')")
+ if (is.null(drty.pest)) drty.pest <- dirname(pst.fname)
+ if (is.null(drty.model)) drty.model <- dirname(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(paramranges.fname)==paramranges.fname)
+ paramranges.fname <- paste(drty.out, "/", paramranges.fname, sep="")
+
+ ##############################################################################
+ # Reading .pst file
+ x <- readLines(pst.fname)
+
+ ##############################################################################
+ # 1.a) 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]
+
+ # 1.b) 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])
+
+
+ ##############################################################################
+ # 2) Getting Param names and Ranges
+ L <- 0
+ params.stg <- "* parameter data"
+ ini.index <- which(x==params.stg)
+
+ L <- length(ini.index)
+ if (L > 0) {
+ suppressWarnings(tmp <- read.table(file=pst.fname,skip=ini.index,header=FALSE,nrows=nparam,colClasses
+ =c("character","NULL","NULL","numeric","numeric","numeric","NULL","NULL","NULL","NULL"),fill=TRUE,na.strings=""))
+
+ param.names <- tmp[,1]
+ param.ini <- as.numeric(tmp[,2])
+ param.min <- as.numeric(tmp[,3])
+ param.max <- as.numeric(tmp[,4])
+ } else stop("Invalid pst file: ", params.stg, " does not exist !")
+
+ # writing 'ParamRanges.txt'
+ .pst2paramranges(drty.model=drty.model, names=param.names, ini=param.ini,
+ min=param.min, max=param.max,
+ fname.out=paramranges.fname)
+
+
+ ##############################################################################
+ # 3) Getting observations
+ L <- 0
+ obs.stg <- "* observation data"
+ ini.index <- which(x==obs.stg)
+
+ L <- length(ini.index)
+ if (L > 0) {
+ suppressWarnings(tmp <- read.table(file=pst.fname,skip=ini.index,header=FALSE,nrows=nobs,colClasses
+ =c("NULL","numeric","NULL","NULL"),fill=TRUE,na.strings=""))
+
+ obs <- tmp[,1]
+ } else stop("Invalid pst file: ", obs.stg, " does not exist !")
+
+ # Writing the output file
+ obs.fname <- paste(drty.out, "/PEST2hydroPSO_OBS.txt", sep="")
+ write.table(obs, file=obs.fname, col.names=FALSE, row.names=FALSE, quote=FALSE)
+
+
+ ##############################################################################
+ # 4) Model command line
+ L <- 0
+ model.stg <- "* model command line"
+ ini.index <- which(x==model.stg)
+
+ L <- length(ini.index)
+ if (L > 0) {
+ model.exe <- x[ini.index+1]
+ } else stop("Invalid pst file: ", model.stg, " does not exist !")
+
+
+ ##############################################################################
+ # 5) Getting Param filenames and locations (.tpl's and .ins's)
+ L <- 0
+ io.stg <- "* model input/output"
+ ini.index <- which(x==io.stg)
+
+ L <- length(ini.index)
+ if (L > 0) {
+ suppressWarnings(tmp <- read.table(file=pst.fname,skip=ini.index,header=FALSE,nrows=ntpl+nins,colClasses
+ =c("character","character"),fill=TRUE,na.strings=""))
+
+ tpls <- tmp[1:ntpl,1]
+ inputs <- tmp[1:ntpl,2]
+ ins <- tmp[(ntpl+1):(ntpl+nins),]
+ } else stop("Invalid pst file: ", io.stg, " does not exist !")
+
+ # Writing ParamFiles.txt
+ .pst2paramfiles(drty.model=drty.model, tpls=tpls, inputs=inputs,
+ param.names=param.names, fname.out=paramfiles.fname,
+ DecimalPlaces=DecimalPlaces)
+
+
+ ##############################################################################
+ # 6) Creating the Rscript used 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)
+ file.copy(rscript.fname, dst.fname, overwrite=TRUE)
+
+ ##############################################################################
+ # 7) Modifying the Rscript used for running hydroPSO
+
+ # rading the Rscript
+ x <- readLines(dst.fname)
+
+ # Changing model directory
+ x[25] <- paste("model.drty <- \"", drty.model, "\"", sep="")
+
+ # Changing model exe name
+ x[64] <- paste("exe.fname= \"", model.exe, "\"", sep="")
+
+ # Writing the modified file
+ if (file.exists(dst.fname)) file.remove(dst.fname)
+ writeLines(x, con=dst.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), "']")
+
+} # 'pest2hydroPSO' END
diff --git a/inst/Rscripts/hydroPSO-Rscript.R b/inst/Rscripts/hydroPSO-Rscript.R
new file mode 100644
index 0000000000000000000000000000000000000000..b61cfc5ae9b2003a7674ba9c3297c607e92e4993
--- /dev/null
+++ b/inst/Rscripts/hydroPSO-Rscript.R
@@ -0,0 +1,92 @@
+################################################################################
+## File hydroPSO-Rscript.R #
+## Part of the hydroPSO R package: #
+## http://www.rforge.net/hydroPSO/ ; #
+## http://cran.r-project.org/web/packages/hydroPSO #
+## Copyright 2011-2012 Mauricio Zambrano-Bigiarini & Rodrigo Rojas #
+## Distributed under GPL 2 or later #
+## #
+## R script created with he PEST2hydroPSO function #
+## #
+## Created by Mauricio Zambrano-Bigiarini and Rodrigo Rojas. 08-Nov-2012 #
+## Last saved: 08-Nov-2012 #
+################################################################################
+
+###Loading required libraries
+library(hydroPSO)
+library(hydroGOF)
+# library(hydroTSM) # OPTIONAL package
+
+################################################################################
+################## User-defined variables - START ###########################
+################################################################################
+
+###Definition of working directory: input, output and model files paths
+model.drty <- "user.model.drty"
+setwd(model.drty)
+
+###Period of analysis.
+###Only required for out.FUN
+###To define subperiod of analysis to compute the goodness-of-fit (GoF) measures
+Sim.Ini <- "YYYY-MM-DD"
+Sim.Fin <- "YYYY-MM-DD"
+gof.Ini <- "YYYY-MM-DD"
+gof.Fin <- "YYYY-MM-DD"
+
+
+###Function for reading the simulated equivalents
+out.FUN="read_output", # name of user-defined function for reading model outputs
+out.FUN.args=list( ###START out.FUN.args
+ file="filename.out" # name of the output file
+ #,,, # additional arguments for 'out.FUN'
+ ) ###END out.FUN.args
+
+
+###Goodness-of-fit function, either pre-defined from hydroGOF (e.g., ssq) or
+###customized
+gof.FUN <- "ssq" # sum of squared residuals. PEST default.
+ # any other model performance measure could be used
+gof.FUN.args <- list()
+
+################################################################################
+################### User-defined variables - END ############################
+################################################################################
+
+###Getting the OBSERVATIONS
+obs.fname <- "PEST2hydroPSO_OBS.txt"
+obs.fname <- paste(model.drty, "/PSO.in/", obs.fname)
+obs <- read.table(obs.fname)
+
+###MAIN model function
+model.FUN.args=list(
+ model.drty=model.drty,
+ param.files=paste(model.drty,"/PSO.in/ParamFiles.txt",sep=""),
+ exe.fname="ModelCommandLine", #TODO
+ #stdout="",
+ out.FUN=out.FUN,
+ out.FUN.args=out.FUN.args,
+ gof.FUN=gof.FUN,
+ gof.FUN.args=gof.FUN.args,
+ #gof.Ini=gof.Ini, # un-comment if 'gof.Ini' is defined
+ #gof.Fin=gof.Fin, # un-comment if 'gof.Fin' is defined
+ obs=obs
+) ###END model.FUN.args
+
+### MAIN PSO ALGORITHM
+### For hydroPSO fine-tuning parameters, see Zambrano-Bigiarini and Rojas, 2012
+hydroPSO(
+ fn="hydromod",
+ model.FUN="hydromod",
+ model.FUN.args=model.FUN.args,
+ control=list( ###START control options
+ param.ranges="ParamRanges.txt",
+ MinMax="min", # minimisation of ssq
+ #maxit=1000, # case dependent
+ #npart=40, # case dependent
+ #,,, # additional arguments for hydroPSO
+ REPORT=5 # frequency of screen messages
+ ) ###END control options
+) ###END MAIN hydroPSO ALGORITHM
+
+# Plotting the results
+plot_results(MinMax="min", do.png=TRUE)