From dbe6bd3eacda0c89502ce7a109578af9a325bb6e Mon Sep 17 00:00:00 2001
From: Mauricio Zambrano-Bigiarini <hzambran@users.noreply.github.com>
Date: Fri, 9 Nov 2012 16:45:22 +0000
Subject: [PATCH] new function 'pest2hydroPSO.R
---
DESCRIPTION | 4 +-
NAMESPACE | 2 +-
NEWS | 2 +
R/pest2hydroPSO.R | 260 +++++++++++++++++++++++++++++++
inst/Rscripts/hydroPSO-Rscript.R | 92 +++++++++++
5 files changed, 357 insertions(+), 3 deletions(-)
create mode 100644 R/pest2hydroPSO.R
create mode 100644 inst/Rscripts/hydroPSO-Rscript.R
diff --git a/DESCRIPTION b/DESCRIPTION
index 3577344..fb198c9 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 fa9a593..60531c8 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 eac9a5d..80e3391 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 0000000..7a7f557
--- /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 0000000..b61cfc5
--- /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)
--
GitLab