From 47870e2f548938965eaf8419597be8ee6047a2d9 Mon Sep 17 00:00:00 2001
From: Mauricio Zambrano-Bigiarini <hzambran@users.noreply.github.com>
Date: Fri, 15 Jun 2012 16:23:02 +0000
Subject: [PATCH] hydroPSO: fixed small bug related to improper PSO evolution
when 'best.update==async'. source code was tidy up. Added '...' parameter.
---
DESCRIPTION | 4 +-
NEWS | 6 ++
R/PSO_v2012.R | 205 ++++++++++++++++++++++++++++--------------------
man/hydroPSO.Rd | 14 ++--
4 files changed, 136 insertions(+), 93 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index df13e4c..825808e 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-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>
diff --git a/NEWS b/NEWS
index a3b31b8..8ffc02a 100755
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,12 @@
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
diff --git a/R/PSO_v2012.R b/R/PSO_v2012.R
index 4be491e..5de6ac4 100755
--- a/R/PSO_v2012.R
+++ b/R/PSO_v2012.R
@@ -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 #
diff --git a/man/hydroPSO.Rd b/man/hydroPSO.Rd
index 56f2d1f..5ccd7cc 100755
--- a/man/hydroPSO.Rd
+++ b/man/hydroPSO.Rd
@@ -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
--
GitLab