#' Fitness-Euclidean distance ratio particle swarm optimization (FER-PSO) with local search for multimodal optimization
#'
#'
#' In FER-PSO, Fitness and Euclidean distance Ratio (FER)
#' is calculated based on the fitness difference and
#' the Euclidean distance between a particle's personal best and
#' other personal bests of the particles in the population. The
#' key advantage is that FER-PSO removes the need of prespecifying
#' niching parameters that are commonly required in
#' existing niching evolutiobary algorithms for multimodal optimization.\cr
#' Furthermore, a local search technique has been used to enhance the ability to
#' locate most global or local optima.
#'
#' @param fn objective function that should be maximized. Should not return \code{NaN}.
#' @param lower lower bound.
#' @param upper upper bound.
#' @param control control parameters for the algorithm. See "Details".
#' @param ... extra arguments are passed to \code{fn}.
#'
#' @references
#' Qu, B. Y., Liang, J. J., & Suganthan, P. N. (2012). Niching particle swarm optimization with local search for multi-modal optimization. Information Sciences, 197, 131-143.\cr
#'
#' Li, X. (2007, July). A multimodal particle swarm optimizer based on fitness Euclidean-distance ratio. In Proceedings of the 9th annual conference on Genetic and evolutionary computation (pp. 78-85). ACM. \cr
#'
#' Based on MATLAB code that can be found in Suganthan's home page.
#'
#'
#' @details
#'
#' The \code{control} argument is a list that can supply any of the following components:
#' \describe{
#' \item{\code{swarm}}{swarm size. Defaults to \code{50}.}
#' \item{\code{iter}}{number of iterations. Defaults to \code{200}.}
#' \item{\code{w}}{inertia weight. Defaults to \code{0.729843788}.}
#' \item{\code{c1}}{acceleration factor. Defaults to \code{2.05}.}
#' \item{\code{c2}}{acceleration factor. Defaults to \code{2.05}.}
#' \item{\code{local}}{logical; local search should be performed? Defaults to \code{TRUE}}
#' \item{\code{vectorize_local}}{logical; vectorization for local search? Defaults to \code{TRUE}.}
#' \item{\code{seed}}{random seed.}
#' \item{\code{hybrid}}{logical; if true, before quiting the algorithm, an \code{L-BFGS-B}
#' search with the provided position as initial guess is done to improve the accuracy of the results.
#' Defaults to \code{TRUE}. Note that no attempt is done to control the maximal number of
#' function evaluations within the local search step (this can be done separately through \code{hybrid.control})
#' }
#' \item{\code{contr_hybrid}}{List with any additional control parameters to pass on to \code{\link{stats}{optim}}
#' when using \code{L-BFGS-B} for the local search. Defaults to \code{NULL}.}
#' }
#'
#' \code{fn} is maximized.\cr
#' This function is efficient when the purpose is finding all global maxima.
#' Please use \code{\link{ncde}} to find all local maxima.\cr
#' \code{fn} must not return any \code{NaN}.\cr
#' The only stopping rule is the number of iterations.
#' @return
#' a list contains:
#' \describe{
#' \item{\code{pbest}}{a matrix; position of the particles.}
#' \item{\code{pbestval}}{a vector; corresponding fitness value of each particle.}
#' \item{\code{nfeval}}{number of function evaluations.}
#' \item{\code{maxima}}{position of particles after using \code{L-BFGS-B}. It should be local maxima. If \code{control$ == FALSE}, then it is \code{NA}.}
#' \item{\code{maximaval}}{fitness values of \code{maxima}.}
#' }
#' @examples
#'####################################################################################
#'## Two-Peak Trap: global maximum on x = 20 and local maximum on x = 0
#'two_peak <- function(x)
#' y <- (160/15) * (15 - x) * (x < 15) + 40 * (x - 15) * (x >= 15)
#'
#'ferpsols(two_peak, 0, 20, control = list(seed = 66, swarm = 50))
#'
#'# without local search
#'ferpsols(two_peak, 0, 20, control = list(seed = 66, swarm = 50))
#'
#'# without refining and local search
#'ferpsols(two_peak, 0, 20, control = list(seed = 66, swarm = 50, hybrid = FALSE))
#'\dontrun{
#' ####################################################################################
#'# Decreasing Maxima: one global on x = 0.1 and four local maxima
#'dmaxima <- function(x)
#' y <- exp(-2 * log(2) * ((x - 0.1)/0.8)^2) * (sin(5 * pi * x))^6
#'
#'res <- ferpsols(dmaxima, 0, 1, control = list(seed = 66, swarm = 100))
#'unique(round(res$maxima, 5)) ## ferpsols can not find the local maxima
#'
#'## plot
#'x <- seq(0, 1, length.out = 400)
#'plot(x, dmaxima(x), type = "l")
#'
#'####################################################################################
#'# Himmelblau's function: four global optima on
#'# x = c(-2.80512, 3.13131), c(3.00000, 2.00000), c(-3.77931, -3.28319) and c(3.58443, -1.84813)
#'Himmelblau <- function(x){
#' y <- - (x[1]^2 + x[2] - 11)^2 - (x[1] + x[2]^2 - 7)^2
#' return(y)
#'}
#'
#'
#'res <- ferpsols(Himmelblau, c(-6, -6), c(6, 6), control = list(seed = 66, swarm = 50))
#'unique(round(res$maxima, 5))
#'
#'Himmelblau_plot <- function(x, y)
#' Himmelblau(x = c(x, y))
#'Himmelblau_plot <- Vectorize(Himmelblau_plot)
#'x <- y <- seq(-6, 6, length.out = 100)
#'persp(x, y, z = outer(X = x, Y = y, FUN = Himmelblau_plot))
#'
#'####################################################################################
#'# Six-Hump Camel Back: two global and two local maxima
#'Six_Hump <- function(x){
#' factor1 <- (4 - 2.1 * (x[1]^2) + (x[1]^4)/3) * (x[1]^2) + x[1] * x[2]
#' factor2 <- (-4 + 4 * (x[2]^2)) * (x[2]^2)
#' y <- -4 * (factor1 + factor2)
#' return(y)
#'}
#'res <- ferpsols(Six_Hump, c(-1.9, -1.1), c(1.9, 1.1), control = list(seed = 66, swarm = 200))
#'unique(round(res$maxima, 5)) ## can not find the local maxima
#'
#'####################################################################################
#'## 2D Inverted Shubert function :
#'# The global minima: 18 global minima f(x*) = -186.7309.
#'# the local maxima: sevral
#'
#'Shubert <- function(x){
#' j <- 1:5
#' out <- -(sum(j * cos((j + 1) * x[1] + j)) * sum(j * cos((j + 1) * x[2] + j)))
#' return(out)
#'}
#'
#'res <- ferpsols(Shubert, rep(-10, 2), rep(10, 2), control = list(seed = 66, swarm = 200))
#'unique(round(res$maxima, 5)) ## only the global maxima are found
#'
#'## plotting
#'Shubert_plot <- function(x, y)
#' Shubert(x = c(x, y))
#'Shubert_plot <- Vectorize(Shubert_plot)
#'y <- x <- seq(-10, 10, length.out = 40)
#'persp(x, y, z = outer(X = x, Y = y, FUN =Shubert_plot))
#'}
#' @importFrom lhs randomLHS augmentLHS
#' @importFrom stats dist optim runif
#' @export
ferpsols <- function(fn, lower, upper, control = list(), ...){
if (missing(fn))
stop("'fn' is missing")
if (missing(fn))
stop("'lower' is missing")
if (missing(fn))
stop("'upper' is missing")
fn1 <- function(par)
fn(par, ...)
#swarm size
if (is.null(control$swarm))
control$swarm <- 50
if (is.null(control$w)) ## inertia weight
control$w<- 0.729843788
if (is.null(control$c1)) ##accelaration factor accroding to equation 4
control$c1 <- 2.05
if (is.null(control$c2))
control$c2 <- 2.05
if (is.null(control$local))
control$local <- TRUE
if (is.null(control$vectorize_local))
control$vectorize_local <- TRUE
## common param
if (is.null(control$iter))
control$iter <- 200
if (!is.null(control$seed))
set.seed(control$seed)
if (is.null(control$hybrid))
control$hybrid <- TRUE
if (is.null(control$contr_hybrid))
control$contr_hybrid <- NULL
# parameters fixed during the run
# max velocity and minus min veloicty
maxvel <- 0.5 * (upper-lower)
swarm <- control$swarm
####################################################################
## max velocity as a matrix, needed for vectorization
maxvelmat <- matrix(rep(maxvel, swarm), nrow = swarm, byrow = TRUE)
#lower matrix and upper matrix
lowermat <- matrix(rep(lower, swarm), nrow = swarm, byrow = TRUE)
uppermat <- matrix(rep(upper, swarm), nrow = swarm, byrow = TRUE)
#####################################################################
##size of the search space Eq. 7. sos is ||s|| in the numerator of alpha Li (2007)
sos <- sqrt(sum((upper - lower)^2))
D <- length(lower) ##dimension
#####################################################################
## initializing
### we want to add the vertices to the inital positions later
vertices <- make_vertices(lower = lower, upper = upper)
A <- lhs::randomLHS(floor((swarm - dim(vertices)[1])/2), length(lower))
pos <- lhs::augmentLHS(A, ceiling((swarm - dim(vertices)[1])/2))
vrunif <- Vectorize(runif)
# pos <- vrunif (swarm, lower, upper)
#pos[1:dim(vertices)[1], ] <- vertices
pos <- rbind(pos, vertices)
# fitness
fit <- apply(pos, 1, fn1)
##initializing the fitness count
nfeval <- swarm
pbest <- pos
pbestval <- fit
#initialize the velocity of the particles
vel <- lowermat + (2 * uppermat * vrunif(swarm, rep(0, D), rep(1, D)))
#####################################################################
for(iter in 2:(control$iter)){
# global best fitness value and id
gbestval <- max(fit)
gbestid <- which.max(fit)
gworstval <- min(fit)
gworstid <- which.min(fit)
# compute the scaling factor (alpha in equation 7)
sf <- sos/(gbestval-gworstval)
# inital the nbest position
nbest <- pbest
# inital the nbest value
nbestval <- pbestval
# update the nbest and nbestval
##caluclate the euclidean distance for all pbest
eucdist <- dist(pbest, method = "euclidean")
### update_nbest update the nbest and nbest value
nbesttemp <- update_nbest(eucdist = as.matrix(eucdist) , pbest = pbest,
nbest = nbest, pbestval = pbestval, nbestval = nbestval, sf = sf)
nbest <- nbesttemp$nbest
nbestval <- nbesttemp$nbestval
## update velocity
temp1 <- control$c1 * vrunif(swarm, rep(0, D), rep(1, D)) * (pbest - pos) + control$c2 * vrunif(swarm, rep(0, D), rep(1, D)) * (nbest - pos)
vel <- control$w * vel + temp1
temp1 <- NA ## cleaning
#######################################################################################################
# velocity
##check the upper bound
vel <- (vel > maxvelmat) * maxvelmat + (vel <= maxvelmat) * vel
##check the lower bound
vel <- (vel < -maxvelmat) * -maxvelmat + (vel >= -maxvelmat) * vel
########################################################################################################
########################################################################################################
##update position
pos <- pos + vel
## correct out of bounds value
pos <- ((pos >= lowermat) & (pos <= uppermat)) * pos + ##if inside the box return pos
(pos < lowermat) * (lowermat + 0.25 * (uppermat - lowermat) * vrunif(swarm, rep(0, D), rep(1, D))) + #lower bound
(pos > uppermat) * (uppermat - 0.25 * (uppermat - lowermat) * vrunif(swarm, rep(0, D), rep(1, D))) ##upper bound
########################################################################################################
########################################################################################################
# computing the fitness values and updating the pbes and pbestval
fit <- apply(pos, 1, fn1)
nfeval <- swarm + nfeval
temp2 <- (pbestval > fit)
pbestval <- temp2 * pbestval + (1 - temp2) * fit
pbest <- temp2 * pbest + (1 - temp2) * pos
########################################################################################################
########################################################################################################
## local search
if (control$local && control$vectorize_local){
temp_euc <- as.matrix(dist(pbest, method = "euclidean"))
temp_near_id <- sapply(1:swarm, FUN = function(j)which(order(temp_euc[j, ]) == 2))
near <- pbest - pbest[temp_near_id , , drop = FALSE] * (pbestval >= pbestval[temp_near_id]) +
(pbest[temp_near_id , , drop = FALSE] - pbest) * (pbestval < pbestval[temp_near_id])
temp_pbest <- pbest + 2 * matrix(runif(D * swarm), nrow = swarm, ncol = D) * near
temp_pbest <- ((temp_pbest >= lowermat) & (temp_pbest <= uppermat)) * temp_pbest +
(temp_pbest < lowermat) * (lowermat + 0.25 * (uppermat - lowermat) * vrunif(swarm, rep(0, D), rep(1, D))) + #lower bound
(temp_pbest > uppermat) * (uppermat - 0.25 * (uppermat - lowermat) * vrunif(swarm, rep(0, D), rep(1, D)))
temp_pbestval <- apply(temp_pbest, 1, fn1)
nfeval <- nfeval + swarm ## updating the number of fitness evaluations
temp_improve <- temp_pbestval > pbestval
## updating pbest and pbestval if there is improvement
pbest <- pbest * (1-temp_improve) + temp_pbest * temp_improve
pbestval <- pbestval * (1-temp_improve) + temp_pbestval * temp_improve
}
if (control$local && !control$vectorize_local){
##local search
for (k in 1:swarm){
temp_euc <- order(sqrt(apply(sweep(pbest, 2, pbest[k, , drop = FALSE])^2, 1, sum)))
temp_near_id <- which(temp_euc == 2)
## producing a point near the pbest[k]
near <- pbest[k, , drop = FALSE] - pbest[temp_near_id , , drop = FALSE] * (pbestval[k] >= pbestval[temp_near_id]) +
(pbest[temp_near_id , , drop = FALSE] - pbest[k, , drop = FALSE]) * (pbestval[k] < pbestval[temp_near_id])
temp_pbest <- pbest[k, , drop = FALSE] + 2 * runif(D) * near
temp_pbest <- ((temp_pbest >= lower) & (temp_pbest <= upper)) * temp_pbest +
(temp_pbest < lower) * (lower + 0.25 * (upper - lower)*runif(D)) + #lower bound
(temp_pbest > upper) * (upper - 0.25 * (upper - lower)*runif(D)) # upper bound
temp_pbestval <- fn1(temp_pbest)
nfeval <- nfeval + 1
temp_improve <- temp_pbestval > pbestval[k]
## updating the number of fitness evaluations
pbest[k, ] <- pbest[k, , drop = FALSE] * (1-temp_improve) + temp_pbest * temp_improve
pbestval[k] <- pbestval[k] * (1-temp_improve) + temp_pbestval * temp_improve
}
}
} ## end of local search
########################################################################################################
pbest <- pbest[order(pbestval, decreasing = TRUE), , drop = FALSE]
pbestval <- sort(pbestval, decreasing = TRUE)
########################################################################################################
if (control$hybrid){
##### now we imrove all the particles with optim
fn2 <- function(par)
-fn1(par) ## because optim works with min
temp_optim <- t(sapply(1:swarm, FUN = function(j)optim(par = pbest[j, ], fn = fn2, method = "L-BFGS-B",
lower = lower, upper = upper,
control= control$contr_hybrid)))
temp_optim <- t(sapply(1:swarm, function(j) c(temp_optim[j, ]$par, temp_optim[j, ]$value)))
#temp_optim <- round(temp_optim, control$digits)
#temp_optim <- unique(temp_optim)
maxima <- temp_optim[, -(D+1), drop = FALSE]
maximaval <- -temp_optim[, (D+1)]
} else
maxima <- maximaval <- NA
########################################################################################################
return(list(pbest = pbest, pbestval = pbestval, nfeval = nfeval, maxima = maxima, maximaval = maximaval))
}
###########################################################################################
###########################################################################################
# fn <- function(xx)
# {
#
#
# factor1=(4-2.1*(xx[1]^2)+(xx[1]^4)/3)*(xx[1]^2)+xx[1]*xx[2]
# factor2=(-4+4*(xx[2]^2))*(xx[2]^2)
# y=-4*(factor1+factor2)
#
# return(y)
# }
#
# lower <- c(-1.9, -1.1)
# upper <- c(1.9, 1.1)
#
#
# #two maxima
# fn(c(0.089842, -0.712656))
# fn(c(-0.089842, 0.712656))
#
#
# control = list()
#
# time1 <- proc.time()
# Rprof(filename = "ferpso_rproof.out", line.profiling = T)
# test <- ferpso(fn, lower, upper, control = list(local = F, seed = 900, swarm = 500))
# Rprof(NULL)
# summaryRprof("ferpso_rproof.out", lines = "show")
# time1 <- proc.time() - time1
###########################################################################################
###########################################################################################
# update_nbest(eucdist = eucdist,
# pbest = pbest, nbest = nbest, pbestval = pbestval, nbestval = nbestval, sf = sf)
# Two-Peak Trap
#ferpso(fn1, 0, 20, control = list(local = TRUE, seed = 66, swarm = 50))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.