From ca415b7f711bc2df64c9a52c69363e96956941eb Mon Sep 17 00:00:00 2001
From: Mauricio Zambrano-Bigiarini <hzambran@users.noreply.github.com>
Date: Fri, 9 Nov 2012 00:56:52 +0000
Subject: [PATCH] hydroPSO.R: new regrouping ready (tested with random and
 gbest topologies)

---
 DESCRIPTION     |   4 +-
 NEWS            |   6 ++-
 R/PSO_v2012.R   | 121 +++++++++++++++++++++++-------------------------
 man/hydroPSO.Rd |  15 ++++--
 4 files changed, 74 insertions(+), 72 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index b9ee9f8..3577344 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-19
-Date: 2012-11-06
+Version: 0.1-58-20
+Date: 2012-11-08
 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/NEWS b/NEWS
index 5e13086..eac9a5d 100755
--- a/NEWS
+++ b/NEWS
@@ -18,6 +18,10 @@ 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
+                            -) 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 
+                               radius), from 1.1E-4, 0.8 and 5 to 1E-5, 2 and 100, respectively
                             -) argument 'method' now allows the following 3 new values: 'spso2007', 'spso2011' and 'canonical', which set the values of the PSO 
                                engine to the ones of the SPSO 2007, SPSO 2011 and the canonical one, respectively. 
                             -) argument 'method'. The old (default) value 'pso' was replaced by the new (default) value 'spso2011'
@@ -53,7 +57,7 @@ NEWS/ChangeLog for hydroPSO
                                  v[t+1] :   v[t]           ->      -v[t]                 (when x[t+1] > x_max | x[t+1] < x_min )
                             -) when the control argument 'out.with.fit.iter' (not used so far) is set to TRUE, the number of iterations returned
                                correspond to the effective number of iterations carried out, not the maximum defined by 'maxit'
-                            -) now it is able to correctly write the observed values for sub-daily models
+                            -) observed values are now correctly written to disk for sub-daily models
                             -) the package now depends on R >= 2.13.0, not 2.10.0 as in previous releases, in order to be consistent with the 
                                byte-compiler setting of the package
                             -) fixed some (very unlikely) error when 'IW.type="aiwf"'
diff --git a/R/PSO_v2012.R b/R/PSO_v2012.R
index fb8a0ba..975e406 100755
--- a/R/PSO_v2012.R
+++ b/R/PSO_v2012.R
@@ -1060,13 +1060,13 @@ UpdateNgbest <- function(pbest.fit, ngbest, MinMax) {
 ################################################################################
 #                    ComputeSwarmRadiusAndDiameter                             #
 ################################################################################
-# Author : Mauricio Zambrano-Bigiarini
-# Started: 12-Jan-2011
-# Updates: 12-Jan-2011 ; 28-Oct-2011
-#          06-Nov-2012 ; 07-Nov-2012
+# Author : Mauricio Zambrano-Bigiarini                                         #
+# Started: 12-Jan-2011                                                         #
+# Updates: 12-Jan-2011 ; 28-Oct-2011                                           #
+#          06-Nov-2012 ; 07-Nov-2012                                           #
 ################################################################################
-# Purpose: Function for computing the swarm radius, for detecting premature 
-#          convergence, in order to avoid stagnation 
+# Purpose: Function for computing the swarm radius, for detecting premature    #
+#          convergence, in order to avoid stagnation                           #
 ################################################################################
 # X      : matrix, with the current parameter values for all the particles
 # gbest  : numeric, with the parameter values for the best particle in the swarm
@@ -1109,17 +1109,25 @@ ComputeSwarmRadiusAndDiameter <- function(x, gbest, Lmax) {
 # Updates: 18-Nov-2011                                                         #
 #          06-Nov-2012 ; 07-Nov-2012 ; 08-Nov-2012                             #
 ################################################################################
-# Purpose: Function for regrouping the swarm in a search space centred around 
-#          the global best, which is hoped to be both, small enough for  
-#          efficient search and large enough to allow the swarm to escape from 
-#          the current local best
-################################################################################
-# Reference: Evers, G.I.; Ben Ghalia, M. 2009. Regrouping particle swarm 
-#            optimization: A new global optimization algorithm with improved  
-#            performance consistency across benchmarks. 
-#            Systems, Man and Cybernetics, 2009. SMC 2009. 
-#            IEEE International Conference on, vol., no., pp.3901-3908,
-#            doi: 10.1109/ICSMC.2009.5346625
+# Purpose: Function for regrouping the swarm in a search space centered around #
+#          the global best, which is hoped to be both, small enough for        #
+#          efficient search and large enough to allow the swarm to escape from #
+#          the current local best.                                             #
+#          There are 4 differences wrt Evers and Ghalia 2009:                  #
+#          -) swarm radius: median is used instead of max                      #
+#          -) computation of the new range of parameter space, which           #
+#             corresponds to the maximum boundary of all the swarm, instead of #
+#             abs(x-Gbest)                                                     #
+#          -) regrouping factor: user-defined instead of 6/(5*ro)              #
+#          -) velocity is re-initialized using Vini.type instead of using the  #
+#             formula proposed by Evers and Ghalia 2009                        #
+################################################################################
+# Reference: Evers, G.I.; Ben Ghalia, M. 2009. Regrouping particle swarm       #
+#            optimization: A new global optimization algorithm with improved   #
+#            performance consistency across benchmarks.                        #
+#            Systems, Man and Cybernetics, 2009. SMC 2009.                     #
+#            IEEE International Conference on, vol., no., pp.3901-3908,        #
+#            DOI: 10.1109/ICSMC.2009.5346625                                   #
 ################################################################################
 RegroupingSwarm <- function(x, 
                             xini.type, 
@@ -1140,31 +1148,25 @@ RegroupingSwarm <- function(x,
   # Regrouping factor 
   rf <- RG.r          # user-defined
   #rf <- 6/(5*RG.thr) # Evers & Ghalia
-  #rf <- (1/RG.thr)/2 # MZB
   
   # name of each parameter  
   param.IDs <- row.names(x.Range)
 
   # Removing possible attributes
-  gbest      <- as.numeric( gbest ) 
   x.min.rng  <- as.numeric( x.Range[ ,1] )
   x.max.rng  <- as.numeric( x.Range[ ,2] )
   
+  # Computing current boundaries for the whole swarm
   xmin    <- apply(x, MARGIN=2, FUN=min) 
   xmax    <- apply(x, MARGIN=2, FUN=max) 
-  xMinMaxO <- cbind(xmin, xmax)
-  
-  message("Boundaries0  :")
-  print(xMinMaxO)
+  xMinMaxO <- cbind(xmin, xmax)  
+  #message("Boundaries0  :")
+  #print(xMinMaxO)
   
   # Maximum length of the parameter space in each dimension
-  RangeO <- xmax - xmin 
-  
-  message("RangeO  :")
-  print(RangeO)
-  
-  # Maximum length of the parameter space in each dimension
-  #RangeO <- x.max.rng - x.min.rng 
+  RangeO <- xmax - xmin   
+  #message("RangeO  :")
+  #print(RangeO)
 
   # Transforming the 'gbest' into a matrix, in order to make easier some 
   # further computations
@@ -1174,11 +1176,12 @@ RegroupingSwarm <- function(x,
   # Is equal to the product of the regrouping factor with the maximum distance of 
   # each particle to the global best, for each dimension
   #RangeNew <- rf * apply( abs(x-Gbest), MARGIN=2, FUN=max) ## Evers & Ghalia
-  RangeNew <- rf * (xmax - xmin)                            ## MZB
+  RangeNew <- rf * abs(xmax - xmin)                        ## MZB
   
   # Making sure that the new range for each dimension is no larger than the original one
-  RangeNew <- pmin(abs(x.max.rng - x.min.rng), RangeNew)
-  
+  RangeNew <- pmin(abs(x.max.rng - x.min.rng), RangeNew)  
+  #message("RangeNew:")
+  #print(RangeNew)
 
   # Re-initializing particle's positions around gbest
   for (part in 1:npart) {
@@ -1196,41 +1199,31 @@ RegroupingSwarm <- function(x,
   xmin    <- apply(x, MARGIN=2, FUN=min)  ## MZB
   xmax    <- apply(x, MARGIN=2, FUN=max)  ## MZB
   xMinMax <- cbind(xmin, xmax)            ## MZB  
-
-  message("Gbest:")
-  print(gbest)
-  message("BoundariesNew:")
-  print(xMinMax)
-  message("              ")  
-  
-  message("RangeNew:")
-  print(RangeNew)
-  
-  
-  vmin    <- apply(v, MARGIN=2, FUN=min) 
-  vmax    <- apply(v, MARGIN=2, FUN=max) 
-  vMinMax <- cbind(vmin, vmax)
-  message("OldBoundariesV:")
-  print(vMinMax)
-  
-  
-  #x <- InitializateX(npart=npart, x.MinMax=xMinMax, x.ini.type=xini.type)  
+  #message("Gbest:")
+  #print(gbest)
+  #message("BoundariesNew:")
+  #print(xMinMax)   
+  
+  # Printing old velocities
+  #vmin    <- apply(v, MARGIN=2, FUN=min) 
+  #vmax    <- apply(v, MARGIN=2, FUN=max) 
+  #vMinMax <- cbind(vmin, vmax)
+  #message("OldBoundariesV:")
+  #print(vMinMax)
+  
+  # Re-initializing velocities
   v <- InitializateV(npart=npart, x.MinMax=xMinMax, v.ini.type=vini.type, Xini=x)
-  #v <- InitializateV(npart=npart, x.MinMax=xMinMaxO, v.ini.type=vini.type, Xini=x) 
-  #v <- v
-  
-  vmin    <- apply(v, MARGIN=2, FUN=min) 
-  vmax    <- apply(v, MARGIN=2, FUN=max) 
-  vMinMax <- cbind(vmin, vmax)
   
-  message("NewBoundariesV:")
-  print(vMinMax)
-  
-  #v <- InitializateV(npart=npart, x.MinMax=xMinMax, v.ini.type=vini.type, Xini=x)
+  # Printing new velocities
+  #vmin    <- apply(v, MARGIN=2, FUN=min) 
+  #vmax    <- apply(v, MARGIN=2, FUN=max) 
+  #vMinMax <- cbind(vmin, vmax)  
+  #message("NewBoundariesV:")
+  #print(vMinMax)
   
   # Relative change achieved in each dimension
-  rel.change        <- (RangeNew-RangeO)/RangeO
-  names(rel.change) <- param.IDs 
+  #rel.change        <- (RangeNew-RangeO)/RangeO
+  #names(rel.change) <- param.IDs 
 
   out      <- list(3)
   out[[1]] <- x
diff --git a/man/hydroPSO.Rd b/man/hydroPSO.Rd
index fe916b4..76f2e94 100755
--- a/man/hydroPSO.Rd
+++ b/man/hydroPSO.Rd
@@ -279,23 +279,28 @@ numeric, non-linear modulation index \cr
 When \code{TVlambda.type='linear'} \code{TVlambda.exp} is automatically set to 1 
 }
   \item{use.RG}{
-logical, indicates if the swarm should be regrouped when premature convergence is detected. When \code{use.RG=TRUE} the swarm is regrouped in a search space centred around the current global best. This updated search space is hoped to be both small enough for efficient search and large enough to allow the swarm to escape from stagnation (see Evers and Ghalia, 2009)\cr
-By default \code{use.RG=FALSE}
+logical, indicates if the swarm should be regrouped when premature convergence is detected. By default \code{use.RG=FALSE} \cr
+When \code{use.RG=TRUE} the swarm is regrouped in a search space centred around the current global best. This updated search space is hoped to be both small enough for efficient search and large enough to allow the swarm to escape from stagnation (see Evers and Ghalia, 2009)\cr
+There are 4 differences wrt Evers and Ghalia 2009: \cr
+-) swarm radius: median is used instead of max \cr
+-) computation of the new range of parameter space, which corresponds to the maximum boundaries of the whle swarm, instead of \sQuote{abs(x-Gbest)} \cr
+-) regrouping factor: \code{RG.r} instead of \sQuote{6/(5*ro)} \cr
+-) velocity is re-initialized using \code{Vini.type} instead of using the formula proposed by Evers and Ghalia 2009 
 }
   \item{RG.thr}{
 ONLY required when \code{use.RG=TRUE} \cr
 numeric, positive number representing the \kbd{stagnation threshold} used to decide whether the swarm has to be regrouped or not. See Evers and Galia (2009) for further details \cr
 Regrouping occurs when the \kbd{normalised swarm radius} is less than \code{RG.thr}\cr
-By default \code{RG.thr=1.1e-4}
+By default \code{RG.thr=1E-5}
 }
   \item{RG.r}{
 ONLY required when \code{use.RG=TRUE}. \cr
 numeric, positive number representing the \kbd{regrouping factor}, which is used to regroup the swarm in a search space centred around the current global best (see Evers and Galia, 2009 for further details)\cr
-By default \code{RG.thr=0.8}
+By default \code{RG.thr=2}
 }
   \item{RG.miniter}{
 ONLY required when \code{use.RG=TRUE} \cr
-numeric, minimum number of iterations needed before regrouping. By default \code{RG.miniter=5}
+numeric, minimum number of iterations needed before each new regrouping. By default \code{RG.miniter=100}
 }
 %%  \item{use.DS}{
 %%CPSO
-- 
GitLab