Nothing
# File hydroPSO2pest.R
# Part of the hydroPSO package, http://www.rforge.net/hydroPSO/
# Copyright 2012-2012 Mauricio Zambrano-Bigiarini & Rodrigo Rojas
# Distributed under GPL 2 or later
################################################################################
## .read_paramrange ##
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
# Updates: 15-Nov-2012 ##
################################################################################
# Purpose : To read the 'ParamRanges.txt' hydroPSO input file ##
################################################################################
.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)
return(x)
} # '.read_paramranges' END
################################################################################
## .read_paramfiles ##
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
# Updates: 15-Nov-2012 ##
################################################################################
# Purpose : To read the 'ParamFiles.txt' hydroPSO input file ##
################################################################################
.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
################################################################################
## .read_observations ##
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
# Updates: 15-Nov-2012 ; 28-Nov-2012 ##
################################################################################
# Purpose : To read the 'Observations.txt' hydroPSO output file ##
################################################################################
.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)
if (NCOL(x) == 1) {
x <- read.table(fname, header=FALSE)
} else x <- coredata(read.zoo(fname, header=FALSE)) # zoo::coredata ; zoo::read.zoo
return(as.numeric(x))
} # '.read_observations' END
################################################################################
## .modify.tpl ##
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
# Updates: 15-Nov-2012 ##
################################################################################
# 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 ##
################################################################################
.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))
# 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)
if (verbose) message("[ File '", basename(tpl.fname), "' successfully modified ]")
} # '.modify.tpl' END
################################################################################
## hydroPSO2pest ##
################################################################################
# Author : Mauricio Zambrano-Bigiarini & Rodrigo Rojas ##
################################################################################
# Created: 12-Nov-2012 ##
# Updates: 15-Nov-2012 ; 28-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(
param.files="ParamFiles.txt",
param.ranges="ParamRanges.txt",
observations.fname="Observations.txt",
exe.fname,
drty.model=getwd(),
pst.fname="hydroPSO2PEST.pst",
verbose=TRUE
) {
if (missing(exe.fname))
stop("Missing argument: 'exe.fname'")
if (!file.exists(drty.model)) stop("Invalid argument: 'drty.model' does not exist !")
##############################################################################
if (basename(param.files)==param.files)
param.files <- paste(drty.model, "/PSO.in/", param.files, 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="")
if (basename(pst.fname)==pst.fname)
pst.fname <- paste(drty.model, "/", pst.fname, sep="")
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
if (verbose) message(" ]")
if (verbose) message("[ 1) Reading '", basename(param.ranges), "' ... ]")
x <- .read_paramranges(fname=param.ranges)
# 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)
if (verbose) message("[ Number of parameters found :", nparam, " ]")
##############################################################################
# 2) Reading Observations.txt
if (verbose) message(" ]")
if (verbose) message("[ 2) Reading observations ('", basename(observations.fname), "') ... ]")
obs <- .read_observations(fname=observations.fname)
# Number of observations
nobs <- length(obs)
if (verbose) message("[ Number of observations found:", nobs, " ]")
##############################################################################
# 3) Reading ParamFiles.txt
if (verbose) message(" ]")
if (verbose) message("[ 3) Reading '", basename(param.files), "' ... ]")
x <- .read_paramfiles(fname=param.files)
# 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)
if (verbose) message("[ Number of .tpl to be created:", ntpls, " ]")
##############################################################################
# 4) tpl creation
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="")
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)
# Adding a 1st line with the token
token <- "#"
first.line <- paste("ptf ", token, sep="")
tmp <- c(first.line, tmp)
# writing the .tpl to disk, with exactly the same content as the original input file
writeLines(tmp, con=tpl.fname)
# 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***"
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
if (verbose) message(" ]")
if (verbose) message("[ 6) Creating the .pst 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)
##############################################################################
# 7) modifying the .pst
if (verbose) message(" ]")
if (verbose) message("[ 7) Writing the .pst file ... ]")
x <- readLines(pst.fname)
L <- length(x)
# 7.1) Number of parameters
if (verbose) message("[ 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
# 7.2) Number of observations
if (verbose) message("[ 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
# 7.3) Number of parameter groups
if (verbose) message("[ 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
# 7.4) Number of .tpl
if (verbose) message("[ Writing number of .tpl files ... ]")
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
# 7.5) * parameter groups
if (verbose) message("[ 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) {
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
# 7.6) * parameter data
if (verbose) message("[ 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) {
# 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
# 7.7) * observation data
if (verbose) message("[ 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) {
# 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
# 7.8) * model command line
if (verbose) message("[ Writing '* model command line' section ... ]")
x[20+2*(nparam-1)+(nobs-1)] <- exe.fname
# 7.9) * model input/output
if (verbose) message("[ 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) {
# 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("[ ]")
message("[ PEST input files available in: ]")
message("[ '", drty.model, "']")
message("[ ]")
message("[ Before running PEST, you MUST check *.pst and *.ins files ]")
} # 'pest2hydroPSO' END
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.