emoptim: Evolutionary Multimodal Optimimization

#'  Neighborhood based crowding Differential Evolution (NCDE)  for multimodal optimization
#'
#'
#'  A niching differential evolution (DE) algorithm with neighborhood mutation strategy to
#'   solve multimodal optimization problems.
#'
#'
#'
#'
#' @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., Suganthan, P. N., & Liang, J. J. (2012). Differential evolution with neighborhood mutation for multimodal optimization. IEEE transactions on evolutionary computation, 16(5), 601-614.\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{iter}}{number of iterations. Defaults to \code{200}.}
#'   \item{\code{SF}}{F is a scale factor in \eqn{[0, 2]} used for scaling the differential vector. Defaults to \code{0.9}.}
#'   \item{\code{CR}}{crossover probability from \eqn{[0, 1]}. Defaults to \code{0.2}.}
#'   \item{\code{NP}}{Population size. Must be larger than 8. Defaults to \code{10 * length(lower)}. }

#'   \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}.}
#' }
#'
#' The DE strategy is 'DE/rand/1'.\cr
#' \code{fn} is maximized.\cr
#' The benchmarks show that \code{\link{ncde}} is more efficient than \code{\link{ferpsols}}
#' in finding all local maxima. \cr
#' \code{fn} must not return any \code{NaN}.\cr
#' The only stopping rule is the number of iterations.\cr
#'
#' The implemented algorithm here slightly differs from the original
#'  MATLAB code in a way that the off-springs for all \code{NP} members calculated by vectorization (\code{sapply})
#'   and not with \code{for} loop as the original MATLAB code.\cr
#' @return
#'  a list contains:
#' \describe{
#'   \item{\code{pop}}{a matrix; position of the population.}
#'   \item{\code{popval}}{a vector; corresponding population values.}
#'   \item{\code{nfeval}}{number of function evaluations.}
#'   \item{\code{maxima}}{position of particles after using \code{L-BFGS-B} when \code{control$hybrid = TRUE}. It should be local maxima. If \code{control$hybrid == FALSE}, then it is \code{NA}.}
#'   \item{\code{maximaval}}{fitness values of \code{maxima}.}
#' }
#'
#'
#' The concept of
#' niching is inspired by the way organisms evolve in nature.
#' When integrated with EAs, niching involves the formation
#' of subpopulations within a population. Each subpopulation
#' aims to locate one optimal solution and together the whole
#' population is expected to locate multiple optimal solutions in a single run (Qu et al., 2012).\cr
#'
#' @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)
#'
#'ncde(two_peak, 0, 20, control = list(seed = 66, NP = 50))
#'
#'# without local search
#'ncde(two_peak, 0, 20, control = list(seed = 66, NP = 50))
#'
#'# without refining and local search
#'ncde(two_peak, 0, 20, control = list(seed = 66, NP = 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 <- ncde(dmaxima, 0, 1, control = list(seed = 66, NP = 10))
#'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 <- ncde(Himmelblau, c(-6, -6), c(6, 6), control = list(seed = 66, NP = 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 <- ncde(Six_Hump, c(-1.9, -1.1), c(1.9, 1.1), control = list(seed = 66, NP = 10))
#'unique(round(res$maxima, 5))
#'
#'####################################################################################
#'## 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 <- ncde(Shubert, rep(-10, 2), rep(10, 2), control = list(seed = 66, NP = 70))
#'unique(round(res$maxima, 5)) ## all global maxima but not local maxima
#'
#'## 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))
#'}
#' @export


ncde <-  function(fn, lower, upper, control = NULL,...){
  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, ...)
  #fn1 <- fn # delet later



  # differential weighting factor from interval [0,2]
  if (is.null(control$SF))
    control$SF <- 0.9
  ### algorithms parameter
  #if (is.null(control$st))
  control$st <- 2 ## strategy in DE algorithm
  # crossover probability from interval [0,1]
  if (is.null(control$CR))
    control$CR <- 0.2
  if (is.null(control$NP)) ## humber of intial solutions
    control$NP <- 10 * length(lower)
  NP <- control$NP
  if (NP <= 8)
    stop("'NP' must be larger than 8")

  ### comomon param
  if (is.null(control$iter))
    control$iter <- 200
  maxiter <- control$iter
  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

  ##############################################################################
  ##initializing
  ## creating the lower matrix, required for vectorization
  lowermat <- matrix(rep(lower, NP), nrow = NP, byrow = TRUE)
  uppermat <- matrix(rep(upper, NP), nrow = NP, byrow = TRUE)
  D <- length(lower) ##dimension
  ## initialazing the population
  pop <- lowermat + (uppermat - lowermat)* matrix(runif(D * NP), ncol = D, nrow = control$NP)
  # adding the vertices
  vertices <- make_vertices(lower = lower, upper = upper)
  pop[1:dim(vertices)[1], ] <- vertices
  ## initial values
  popval <- apply(pop, 1, fn1)
  nfeval <- NP # umber of function evaluation up to now
  ibest <- which.max(popval) # best solution in population
  ##############################################################################


  for(ii in 2:maxiter){

    ############################################################################
    # neighborhood size
    ## 'm' in table 8
    # In neighborhood mutation, there is only one parameter
    # m which is the neighborhood size. This parameter controls
    # how many individuals are selected in each subpopulation
    # Generally, m should be chosen between 1/20 of the population
    # size and 1/5 of the population size. It can also be dynamically
    # set, from a relatively large value to a small value. Different
    # Fig. 1. Illustration of neighborhood mutation.
    # from other niching parameters, neighborhood size is easy to
    # choose as it can be made proportional to the population size.
    # Moreover, the performance of algorithm is not sensitive to the
    # change of the neighborhood size as evidenced in Table XXI (NCDE).


    if ((NP <= 8))
      neigh_size <- 2 + 2 * ((maxiter - ii)/maxiter) else
        if (NP <= 200)
          neigh_size <- 5 + 5 * ((maxiter - ii)/maxiter) else
            neigh_size <- 20 + 30 * ((maxiter- ii)/maxiter)
          # if neigh_size is less than 5 it produces an error in DE when we want to sample pm_id
          neigh_size <- floor(neigh_size)
          ## neigh_size is the number of elemenst in subpoulation
          ##########################################################################

          ##########################################################################
          # computing the euclidean distance
          euc_dist  <- as.matrix(dist(pop, method = "euclidean"))
          ##########################################################################

          ###########################################################################
          ## producing the off spring u_i
          ## producing the ofsprings u_i in each subpopulation_i that is pop[order(euc_dist[j, ]), ][1:neigh_size,]
          ui <- t(sapply(1:NP, function(j)DE(newpop = pop[order(euc_dist[j, ]), ,drop = FALSE][1:neigh_size, , drop = FALSE],
                                             bestpop = pop[j, , drop = FALSE],
                                             target = pop[j, , drop = FALSE],
                                             st = control$st, ## strategy
                                             SF = control$SF,
                                             CR = control$CR,
                                             lower = lower,
                                             upper = upper,
                                             dimension = D)))
          if (D == 1)
            ui <- t(ui)
          ##########################################################################

          ##########################################################################
          # compute the cost function for u_i
          tempval <- apply(ui, MARGIN = 1, fn1)
          nfeval <- nfeval + NP ## update number of function evaluations.
          ##########################################################################

          ##########################################################################
          for(i in 1:NP){

            ########################################################################
            ## randomly perumted the population
            pp <- sample(1:NP, size = NP, replace = FALSE)
            qq <- pop[pp, , drop = FALSE]
            ## calulate the distance of each ui to qq and find the smallest one
            min_id <- which.min(sqrt(apply(sweep(qq, 2, ui[i, , drop = FALSE])^2, 1, sum)))

            if (popval[pp[min_id]]< tempval[i]){
              pop[pp[min_id], ] <- ui[i, , drop = FALSE]
              popval[pp[min_id]] <- tempval[i]
            }
          }
          ##########################################################################

          ##########################################################################
          #### vectorization of comparing the ui with population


          ##########################################################################

          # update index of the global best
          ibest <- which.max(popval)

  }


  ##############################################################################
  # sorting the pop and its corresponding values
  pop <- pop[order(popval, decreasing = TRUE), , drop = FALSE]
  popval <- sort(popval, decreasing = TRUE)
  ##############################################################################

  if (control$hybrid){
    fn2 <- function(par)
      -fn1(par) ## because optim works with min
    #
    temp_optim <- t(sapply(1:NP, FUN =  function(j)optim(par = pop[j, ], fn = fn2, method = "L-BFGS-B", lower = lower, upper = upper,
                                                         control = control$contr_hybrid)))
    temp_optim <- t(sapply(1:NP, 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(pop = pop, popval = popval, nfeval = nfeval, maxima = maxima, maximaval = maximaval))
}