Skip to content
Snippets Groups Projects
Commit 47870e2f authored by Mauricio Zambrano-Bigiarini's avatar Mauricio Zambrano-Bigiarini
Browse files

hydroPSO: fixed small bug related to improper PSO evolution when... 

hydroPSO: fixed small bug related to improper PSO evolution when 'best.update==async'. source code was tidy up. Added '...'  parameter.
parent d392f96a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
+ 2
− 2
View file @ 47870e2f
Package: hydroPSO
Type: Package
Title: Model-Independent Particle Swarm Optimisation for Environmental Models
Version: 0.1-56
Date: 2012-06-14
Version: 0.1-56-1
Date: 2012-06-15
Author: Mauricio Zambrano-Bigiarini [aut, cre] and Rodrigo Rojas [ctb]
Author@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>
......
NEWS
+ 6
− 0
View file @ 47870e2f
NEWS/ChangeLog for hydroPSO
--------------------------
0.1-57 XX-Jun-2012
o 'hydroPSO' : -) added '...' parameter. It is only used when 'fn' is different from "hydromod". This is only done for 'optim' compatibility.
-) fixed small bug related to improper PSO evolution when 'best.update'=="async" (NOT the default option !)
AND (topology!="gbest" | method!="pso")
-) source code was tidy up
0.1-56 14-Jun-2012
o 'hydroPSO' : -) much less memory consumption for large values of 'maxit' and/or 'npart' (thanks to Paul Smith for reporting it !)
The files 'BestParamPerIter.txt', 'PbestPerIter.txt', 'LocalBestPerIter.txt', 'Velocities.txt' are now written at the end of
......
R/PSO_v2012.R
+ 120
− 85
View file @ 47870e2f
......@@ -1104,7 +1104,7 @@ Random.Topology.Generation <- function(npart, K,
# Updates: Dec-2010 #
# May-2011 ; 28-Oct-2011 ; 14-Nov-2011 ; 23-Nov-2011 ; #
# 15-Jan-2012 ; 23-Jan-2012 ; 30-Jan-2012 ; 23-Feb-2012 ; 23-Mar-2012 #
# 14-Jun-2012 #
# 14-Jun-2012 ; 15-Jun-2012 #
################################################################################
# 'lower' : minimum possible value for each parameter
# 'upper' : maximum possible value for each parameter
......@@ -1326,6 +1326,7 @@ Random.Topology.Generation <- function(npart, K,
hydroPSO <- function(
par,
fn="hydromod",
...,
method=c("pso", "ipso", "fips", "wfips"),
lower=-Inf,
upper=Inf,
......@@ -2045,13 +2046,48 @@ hydroPSO <- function(
X.neighbours <- Random.Topology.Generation(npart, K, drty.out, iter)
ModelOut <- vector("list", npart)
# IW: linear, non-linear, runif
if (!use.IW) {
w <- 1
} else {
if ( (IW.type == "linear") | (IW.type == "non-linear") ) {
w <- compute.value.with.iter(iter=iter.tv, niter=niter.tv,
nexp=IW.exp, val.ini=w.ini,
val.fin=w.fin)
} else if (IW.type == "runif") {
w <- runif(1, min=w.ini, max=w.fin)
} # ELSE end
} # ELSE end
# TVc1: linear, non-linear
if ( (use.TVc1) & (TVc1.type != "GLratio") ) {
if ( (TVc1.type == "linear") | (TVc1.type == "non-linear") )
c1 <- compute.value.with.iter(iter=iter.tv, niter=niter.tv,
nexp=TVc1.exp, val.ini=c1.ini,
val.fin=c1.fin)
} # If end
# TVc2
if (use.TVc2)
c2 <- compute.value.with.iter(iter=iter.tv, niter=niter.tv,
nexp=TVc2.exp, val.ini=c2.ini,
val.fin=c2.fin)
# lambda
if (use.TVlambda) {
lambda <- compute.value.with.iter(iter=iter.tv, niter=niter.tv,
nexp=TVlambda.exp, val.ini=vmax.ini,
val.fin=vmax.fin)
Vmax <- lambda*Lmax
} # IF end
##########################################################################
# 3.a) Evaluate the particles fitness
if ( fn.name != "hydromod" ) {
# Evaluating a Test Function
Xt.fitness[iter, 1:npart] <- apply(X, fn, MARGIN=1)
Xt.fitness[iter, 1:npart] <- apply(X, fn, MARGIN=1, ...)
GoF <- Xt.fitness[iter, 1:npart]
ModelOut[1:npart] <- GoF ###
......@@ -2077,9 +2113,9 @@ hydroPSO <- function(
model.FUN.args <- modifyList(model.FUN.args, list(param.values=X[part,]) )
hydromod.out <- do.call(model.FUN, as.list(model.FUN.args))
Xt.fitness[iter, part] <- as.numeric(hydromod.out[["GoF"]])
GoF <- Xt.fitness[iter, part]
ModelOut[[part]] <- hydromod.out[["sim"]]
Xt.fitness[iter, part] <- as.numeric(hydromod.out[["GoF"]])
GoF <- Xt.fitness[iter, part]
ModelOut[[part]] <- hydromod.out[["sim"]]
if(is.finite(GoF)) nfn <- nfn + 1
......@@ -2130,51 +2166,28 @@ hydroPSO <- function(
X.best.part <- tmp[["x.best"]]
gbest.fit <- tmp[["gbest.fit"]]
gbest.pos <- tmp[["gbest.pos"]]
} # IF end
suppressWarnings(ifelse(MinMax=="max", pbest.fit.iter <- max( Xt.fitness[iter, ], na.rm=TRUE ),
pbest.fit.iter <- min( Xt.fitness[iter, ], na.rm=TRUE) ) )
tmp <- UpdateLocalBest(pbest.fit=pbest.fit,
tmp <- UpdateLocalBest(pbest.fit=pbest.fit,
localBest.pos=LocalBest.pos,
localBest.fit=LocalBest.fit,
x.neighbours=X.neighbours,
MinMax=MinMax)
LocalBest.fit <- tmp[["localBest.fit"]]
LocalBest.pos <- tmp[["localBest.pos"]]
if ( method == "ipso" ) {
tmp <- UpdateNgbest(pbest.fit=pbest.fit,
ngbest=ngbest,
MinMax=MinMax)
ngbest.fit <- tmp[["ngbest.fit"]]
ngbest.pos <- tmp[["ngbest.pos"]]
} # IF end
GPbest.fit.rate <- mean(pbest.fit, na.rm=TRUE)
ifelse( (is.finite(GPbest.fit.rate) ) & ( GPbest.fit.rate !=0 ),
GPbest.fit.rate <- abs( ( gbest.fit - GPbest.fit.rate ) / GPbest.fit.rate ),
GPbest.fit.rate <- +Inf)
LocalBest.fit <- tmp[["localBest.fit"]]
LocalBest.pos <- tmp[["localBest.pos"]]
ifelse( (gbest.fit.prior != 0) & (is.finite(gbest.fit.prior) ) ,
gbest.fit.rate <- abs( ( gbest.fit - gbest.fit.prior ) / gbest.fit.prior ),
gbest.fit.rate <- +Inf)
out <- ComputeSwarmRadiusAndDiameter(x=X, gbest= X.best.part[gbest.pos, ], Lmax=Lmax,
MinMax=MinMax, pbest.fit=pbest.fit)
swarm.radius <- out[["swarm.radius"]]
swarm.diameter <- out[["swarm.diameter"]]
NormSwarmRadius <- swarm.radius/swarm.diameter
if ( method == "ipso" ) {
tmp <- UpdateNgbest(pbest.fit=pbest.fit,
ngbest=ngbest,
MinMax=MinMax)
ngbest.fit <- tmp[["ngbest.fit"]]
ngbest.pos <- tmp[["ngbest.pos"]]
} # IF end
if ( (verbose) & ( iter/REPORT == floor(iter/REPORT) ) )
message( "iter:", format(iter, width=6, justify="right"),
" Gbest:", formatC( gbest.fit, format="E", digits=digits, flag=" "),
" Gbest_rate:", format( round(gbest.fit.rate*100, 2), width=6, nsmall=2, justify="left"), "%",
" Iter_best_fit:", formatC(pbest.fit.iter, format="E", digits=digits, flag=" "),
" nSwarm_Radius:", formatC(NormSwarmRadius, format="E", digits=digits, flag=" "),
" |g-mean(p)|/mean(p):", format( round(GPbest.fit.rate*100, 2), width=6, nsmall=2, justify="left"), "%" )
} # IF end
# 'X.bak' is only needed to correctly compute the Normalised Swarm Radious
# for the current iteration
X.bak <- X
##########################################################################
################### Particles Loop (j) - Start ########################
......@@ -2183,8 +2196,9 @@ hydroPSO <- function(
if (write2disk) {
# File 'Model_Out.txt'
GoF <- Xt.fitness[iter, j]
# File 'Model_Out.txt'
if(is.finite(GoF)) {
writeLines(as.character(c(iter, j,
formatC(GoF, format="E", digits=digits, flag=" "),
......@@ -2194,9 +2208,9 @@ hydroPSO <- function(
writeLines("", OFout.Text.file)
# File 'Particles.txt' #
if(is.finite(Xt.fitness[iter, j])) {
if(is.finite(GoF)) {
writeLines(as.character( c(iter, j,
formatC(Xt.fitness[iter, j], format="E", digits=digits, flag=" "), #GoF
formatC(GoF, format="E", digits=digits, flag=" "), #GoF
formatC(X[j, ], format="E", digits=digits, flag=" ")
) ), Particles.TextFile, sep=" ")
} else writeLines(as.character( c(iter, j, "NA",
......@@ -2205,9 +2219,9 @@ hydroPSO <- function(
writeLines("", Particles.TextFile)
# File 'Velocities.txt' #
if(is.finite(Xt.fitness[iter, j])) {
if(is.finite(GoF)) {
writeLines( as.character( c(iter, j,
formatC(Xt.fitness[iter, j], format="E", digits=digits, flag=" "), # GoF
formatC(GoF, format="E", digits=digits, flag=" "), # GoF
formatC(V[j, ], format="E", digits=digits, flag=" ")
) ), Velocities.TextFile, sep=" ")
} else writeLines( as.character( c(iter, j, "NA",
......@@ -2233,16 +2247,30 @@ hydroPSO <- function(
X.best.part[j,] <- tmp[["x.best"]]
gbest.pos <- tmp[["gbest.pos"]]
gbest.fit <- tmp[["gbest.fit"]]
tmp <- UpdateLocalBest(pbest.fit=pbest.fit,
localBest.pos=LocalBest.pos,
localBest.fit=LocalBest.fit,
x.neighbours=X.neighbours,
MinMax=MinMax)
LocalBest.fit <- tmp[["localBest.fit"]]
LocalBest.pos <- tmp[["localBest.pos"]]
} # IF end
if ( method == "ipso" ) {
tmp <- UpdateNgbest(pbest.fit=pbest.fit,
ngbest=ngbest,
MinMax=MinMax)
ngbest.fit <- tmp[["ngbest.fit"]]
ngbest.pos <- tmp[["ngbest.pos"]]
} # IF end
if (use.IW) {
if ( (IW.type == "linear") | (IW.type == "non-linear") ) {
w <- compute.value.with.iter(iter=iter.tv, niter=niter.tv,
nexp=IW.exp, val.ini=w.ini,
val.fin=w.fin)
} else if (IW.type == "aiwf") {
} # IF end
### IW, TVc1, Tv2, lambda
# IW: aiwf, GLratio
if (use.IW) {
if (IW.type == "aiwf") {
w <- compute.w.aiwf(iter.fit= Xt.fitness[iter, ],
particle.pos =j,
gbest.fit=gbest.fit,
......@@ -2251,34 +2279,15 @@ hydroPSO <- function(
MinMax=MinMax
)
} else if (IW.type == "GLratio") {
} else if (IW.type == "GLratio") {
w <- compute.w.with.GLratio(MinMax, gbest.fit, pbest.fit)
} else if (IW.type == "runif") {
w <- runif(1, min=w.ini, max=w.fin)
} # ELSE end
} else w <- 1
if (use.TVc1) {
if ( (TVc1.type == "linear") | (TVc1.type == "non-linear") ) {
c1 <- compute.value.with.iter(iter=iter.tv, niter=niter.tv,
nexp=TVc1.exp, val.ini=c1.ini,
val.fin=c1.fin)
} else if (TVc1.type == "GLratio") {
c1 <- compute.c1.with.GLratio(MinMax, gbest.fit, pbest.fit[j])
} # ELSE IF end
} # If end
if (use.TVc2)
c2 <- compute.value.with.iter(iter=iter.tv, niter=niter.tv,
nexp=TVc2.exp, val.ini=c2.ini,
val.fin=c2.fin)
if (use.TVlambda) {
lambda <- compute.value.with.iter(iter=iter.tv, niter=niter.tv,
nexp=TVlambda.exp, val.ini=vmax.ini,
val.fin=vmax.fin)
Vmax <- lambda*Lmax
} # IF end
} # ELSE end
} # IF end
# TVc1: GLratio
if ( (use.TVc1) & (TVc1.type == "GLratio") )
c1 <- compute.c1.with.GLratio(MinMax, gbest.fit, pbest.fit[j])
########################################################################
# 3.b) Updating the velocity of all the particles
......@@ -2318,8 +2327,34 @@ hydroPSO <- function(
##########################################################################
################### Particles Loop (j) - End ##########################
##########################################################################
gbest.fit.iter[iter] <- gbest.fit
suppressWarnings(ifelse(MinMax=="max", pbest.fit.iter <- max( Xt.fitness[iter, ], na.rm=TRUE ),
pbest.fit.iter <- min( Xt.fitness[iter, ], na.rm=TRUE) ) )
GPbest.fit.rate <- mean(pbest.fit, na.rm=TRUE)
ifelse( (is.finite(GPbest.fit.rate) ) & ( GPbest.fit.rate !=0 ),
GPbest.fit.rate <- abs( ( gbest.fit - GPbest.fit.rate ) / GPbest.fit.rate ),
GPbest.fit.rate <- +Inf)
ifelse( (gbest.fit.prior != 0) & (is.finite(gbest.fit.prior) ) ,
gbest.fit.rate <- abs( ( gbest.fit - gbest.fit.prior ) / gbest.fit.prior ),
gbest.fit.rate <- +Inf)
out <- ComputeSwarmRadiusAndDiameter(x=X.bak, gbest= X.best.part[gbest.pos, ], Lmax=Lmax,
MinMax=MinMax, pbest.fit=pbest.fit)
swarm.radius <- out[["swarm.radius"]]
swarm.diameter <- out[["swarm.diameter"]]
NormSwarmRadius <- swarm.radius/swarm.diameter
if ( (verbose) & ( iter/REPORT == floor(iter/REPORT) ) )
message( "iter:", format(iter, width=6, justify="right"),
" Gbest:", formatC( gbest.fit, format="E", digits=digits, flag=" "),
" Gbest_rate:", format( round(gbest.fit.rate*100, 2), width=6, nsmall=2, justify="left"), "%",
" Iter_best_fit:", formatC(pbest.fit.iter, format="E", digits=digits, flag=" "),
" nSwarm_Radius:", formatC(NormSwarmRadius, format="E", digits=digits, flag=" "),
" |g-mean(p)|/mean(p):", format( round(GPbest.fit.rate*100, 2), width=6, nsmall=2, justify="left"), "%" )
##########################################################################
# Random Generation around gbest, if requested #
......
man/hydroPSO.Rd
+ 8
− 6
View file @ 47870e2f
......@@ -14,7 +14,7 @@ Enhanced Particle Swarm Optimisation algorithm
Particle Swarm Optimisation algorithm to calibrate environmental models. It includes a series of controlling options and PSO variants to improve the performance of the algorithm and customize it to different calibration problems
}
\usage{
hydroPSO(par, fn= "hydromod", method=c("pso", "ipso", "fips", "wfips"),
hydroPSO(par, fn= "hydromod", ..., method=c("pso", "ipso", "fips", "wfips"),
lower=-Inf, upper=Inf, control=list(),
model.FUN=NULL, model.FUN.args=list() )
}
......@@ -25,8 +25,13 @@ OPTIONAL. numeric with a first guess for the parameters to be optimised (\env{n}
All the particles are randomly initialised according to the value of \code{Xini.type}. If the user provides \env{m} parameter sets for \code{par}, they are used to overwrite the first \env{m} parameter sets randomly defined according to the value of \code{Xini.type}. If some elements in \code{par} are non finite (lower than \code{lower} or larger than \code{upper}) they are ignored \cr
}
\item{fn}{
character, name of a valid R function \code{c('sphere','ackley','griewank','rastrigin','rosenbrock','schafferF6')} to be optimised or character value \sQuote{'hydromod'}. When \code{fn='hydromod'} the algorithm uses \code{model.FUN} and \code{model.FUN.args} to extract the values simulated by the model and to compute its corresponding goodness-of-fit function. When \code{fn!='hydromod'} the algorithm uses the value(s) returned by \code{fn} as both model output and its corresponding goodness-of-fit \cr
When \code{fn='hydromod'} the algorithm will optimise the model defined by \code{model.FUN} and \code{model.args}
character, name of a valid R function to be optimised (minimized or maximized) or the character value \sQuote{'hydromod'}. \cr
-) When \code{fn!='hydromod'}, the first argument of \code{fn} has to be a vector of parameters over which optimisation is to take place. It should return a scalar result. When \code{fn!='hydromod'} the algorithm uses the value(s) returned by \code{fn} as both model output and its corresponding goodness-of-fit measure. \cr
-) When \code{fn=='hydromod'} the algorithm will optimise the model defined by \code{model.FUN} and \code{model.args}, which are used to extract the values simulated by the model and to compute its corresponding goodness-of-fit measure.
}
\item{\dots}{
OPTIONAL. Only used when \code{fn!='hydromod'}. \cr
further arguments to be passed to \code{fn}.
}
\item{method}{
character, variant of the PSO algorithm to be used. Valid values are in \code{c('pso', 'ipso', 'fips', 'wfips')}: \cr
......@@ -37,9 +42,6 @@ character, variant of the PSO algorithm to be used. Valid values are in \code{c(
\kbd{wfips}: same implementation as \kbd{fips} method, but the contribution of each particle is weighted according to their goodness-of-fit value (see Mendes et al., 2004) \cr
By default \code{method=pso}
}
%%\item{\dots}{
%%further arguments to be passed to \code{fn}
%%}
\item{lower}{
numeric, lower boundary for each parameter \cr
Note for \code{\link[stats]{optim}} users: in \kbd{hydroPSO} the length of \code{lower} and \code{upper} are used to defined the dimension of the solution space
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment