R/de.R

In globpso: Swarm Intelligence Optimization

#' Differential Evolution Algorithms for Minimization Problems
#'
#' The general-purpose implementation of differential evolution algorithms for minimizing an user-defined objective function.
#'
#' @param objFunc The R or Rcpp compiled objective function. See the example.
#' @param lower The vector of finite lower bounds of the search domain.
#' Its length should be equal to the dimension of the domain space.
#' @param upper The vector of finite upper bounds of the search domain.
#' Its length should be equal to the dimension of the domain space.
#' @param init The vector of initial population.
#' Its length should be equal to the dimension of the domain space.
#' When there are more than one initial vectors, specify \code{init} as
#' a matrix. Each row vector represents one initial point.
#' The default for \code{init} is \code{NULL}.
#' @param fixed The vector of real values and NA values that controls DE to search only 
#' for the NA-valued components. 
#' @param DE_INFO The list of DE parameters generated by \code{getDEInfo()}.
#' @param seed The random seed that controls initial population of DE. The default is \code{NULL}.
#' @param verbose The logical value controls if DE would reports the updating progress. The default is \code{TRUE}.
#' @param ... Further arguments to be passed to \code{objFunc}.
#' @return An List.
#' \describe{
#' \item{par}{ the global best particle.}
#' \item{val}{ the objective function value of the global best particle.}
#' \item{history}{ a vector of objective function values of the global best particle in DE search history.}
#' \item{cputime}{ the computational time in seconds.}
#' }
# @details
# TBD
#' @examples
#' # three-dimensional function
#' objf <- function(x, loc) {
#'   val <- 0
#'   for (i in 1:length(x)) {
#'     val <- val + (x[i] - loc)^2
#'   }
#'   return(val)
#' }
#' 
#' upp_bound <- rep(5, 3)
#' low_bound <- rep(-5, 3)
#' loc_shift <- 1
#'
#' alg_setting <- getDEInfo(nPop = 32, maxIter = 100, deType = "rand-1", 
#'                          sf = 0.5, cr = 0.1)
#' res <- diffevo(objFunc = objf, lower = low_bound, upper = upp_bound, 
#'                DE_INFO = alg_setting, loc = loc_shift)
#' res$par
#' res$val
#' 
#' # C++ function example
#' \donttest{
#' library(Rcpp)
#' library(RcppArmadillo)
#' objf_c <- cppFunction('double objf_c(SEXP x, SEXP loc) {
#'     double val = 0;
#'     double loc_c = (double)Rcpp::as<double>(loc);
#'    arma::rowvec x_c = (arma::rowvec)Rcpp::as<arma::rowvec>(x);
#'    for (arma::uword i = 0; i < x_c.n_elem; i++) {
#'      val += (x_c(i) - loc_c)*(x_c(i) - loc_c);
#'    }
#'     return val;
#'  }', depends = "RcppArmadillo")
#'
#' alg_setting <- getDEInfo(nPop = 32, maxIter = 100, deType = "rand-1", 
#'                          sf = 0.5, cr = 0.1)
#' res_c <- diffevo(objFunc = objf_c, lower = low_bound, upper = upp_bound, 
#'                        DE_INFO = alg_setting, loc = 1)
#' res_c$par
#' res_c$val
#' }
#' @references Storn, R., & Price, K. (1997). Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces. Journal of global optimization, 11, 341-359.
#' @name diffevo
#' @rdname diffevo
#' @importFrom methods hasArg
#' @importFrom Rcpp evalCpp cppFunction sourceCpp
#' @export
diffevo <- function(objFunc, lower, upper, init = NULL, fixed = NULL, 
                    DE_INFO = NULL, seed = NULL, verbose = TRUE, ...) {
  
  stopifnot(all(is.finite(lower)), all(is.finite(upper)), 
            length(lower) == length(upper), all(upper >= lower)
  )
  hasInitPop <- 0
  if (!is.null(init)) {
    if (is.vector(init)) {
      stopifnot(length(lower) == length(init), all(upper >= init), 
                all(is.finite(init)), all(init >= lower))
      init <- matrix(init, 1, length(init))
      hasInitPop <- 1
    } else if (is.matrix(init)) {
      for (i in 1:nrow(init)) {
        stopifnot(length(lower) == length(init[i,]), all(upper >= init[i,]), 
                  all(is.finite(init[i,])), all(init[i,] >= lower))
      }
      hasInitPop <- 1
    }
  }
  if (is.null(fixed)) { 
    fixed <- rep(NA, length(lower))
  } else {
    stopifnot(length(fixed) == length(lower))
  }
  
  if (is.null(DE_INFO)) {
    nPop <- 32
    if (hasInitPop) { nPop = max(32, nrow(init)) }
    DE_INFO <- getDEInfo(nPop = nPop, maxIter = 100)
    if (verbose) message(paste0("Use the default settings for DE See '?getDEInfo'."))
  }
  stopifnot(all(names(DE_INFO) == names(getDEInfo())))
  if (hasInitPop) { stopifnot(DE_INFO$nPop >= nrow(init)) }
  
  environment <- new.env()
  
  # Fill the rest of DE options in DE_INFO
  DE_INFO$varUpper <- upper
  DE_INFO$varLower <- lower
  DE_INFO$dPop <- length(upper)
  DE_INFO$hasInitPop <- hasInitPop
  DE_INFO$initPop <- init
  DE_INFO$fixedDims <- fixed
  # Start
  set.seed(seed)
  cputime <- system.time(
    deOut <- cppDE(objFunc, DE_INFO, environment, FALSE, verbose)
  )[3]
  if (verbose) message(paste0("CPU time: ", round(cputime, 2), " seconds."))
  
  return(
    list(par = deOut$GBest, val = deOut$fGBest, 
         pbest = deOut$PBest,
         fpbest = deOut$fPBest,
         history = deOut$fGBestHist, cputime = cputime)
  )
}


#' Generation function of DE parameter settings
#' 
#' Create a list with DE parameters for Minimization.
#' 
#' @param nPop A integer number of population size in DE algorithm.
#' @param maxIter A integer number of maximal DE iterations.
# @param checkConv A logical value which controls whether DE checks the stopping criterion during updating procedure.
# Specify \code{TRUE} for DE to compute the stopping criterion \eqn{|f'-f|<\varepsilon}
# where \eqn{f'} and \eqn{f} are the objective function values in the previous and current iterations, respectively.
# The default is \code{FALSE}.
#' @param freeRun A number between \eqn{[0,1]} that controls the percentage of DE iterations which are free from examining the
#' stopping criterion, \eqn{|f'-f|<\varepsilon}
#' where \eqn{f'} and \eqn{f} are the objective function values in the previous and current iterations, respectively.
#' The default is 1.0 implying the DE will completely ignore the stopping criterion.
#' Otherwise, the DE checks the stopping criterion after free iterations.
#' @param tol A small value for the tolerance, \eqn{\varepsilon}, in the stopping criterion.
#' For \code{freeRun} smaller than 1.0, the default is \code{1e-6}. Otherwise, this value would not affect the algorithm.
#' @param deType string. The type of DE. This package current supports the following types:
#' \describe{
#' \item{"rand-1"}{ Mutation operation on the current position with one random direction}
#' \item{"rand-2"}{ Mutation operation on the current position with two random directions}
#' \item{"best-1"}{ Mutation operation on the best position with one random direction}
#' \item{"best-2"}{ Mutation operation on the best position with two random directions}
#' \item{"rand_to-best-1"}{ Mutation operation on the current position with direction to the best and one random direction}
#' \item{"rand-to-best-2"}{ Mutation operation on the current position with direction to the best and two random directions}
#' }
#' The default type is `rand-1`.
#' @param sf The value of scaling factor in DE updating procedure. The default is 0.5.
#' @param cr The value of crossover rate in DE updating procedure. The default is 0.1.
#' @return A list of DE parameter settings.
#' @examples
#' DE_INFO <- getDEInfo(nPop = 32, maxIter = 100)
#' @name getDEInfo
#' @rdname getDEInfo
#' @export
getDEInfo <- function(nPop = 32, maxIter = 100, deType = "rand-1",
                      #dPop = NULL, varUpper = NULL, varLower = NULL, checkConv = 0,
                      freeRun = 1.0, tol = 1e-6, sf = 0.5, cr = 0.1
) {
  
  stopifnot(length(nPop) == 1)
  if (!(deType %in% c("rand-1", "rand-2", "best-1", "best-2", "rand-to-best-1", "rand-to-best-2"))) {
    stop("Currently the function supports: \ndeType = 'rand-1' for DE/rand/1 algorithm.\ndeType = 'rand-2' for DE/rand/2 algorithm.\ndeType = 'best-1' for DE/best/1 algorithm.\ndeType = 'best-2' for DE/best/2 algorithm.\ndeType = 'rand-to-best-1' for DE/rand-to-best/1 algorithm.\ndeType = 'rand-to-best-2' for DE/rand-to-best/2 algorithm.
		")
  } else {
    if (deType == "rand-1") {
      typeDE <- 0
      if (nPop < 4) { stop("population size must > 4 for DE/rand/1 algorithm.") }
    } else if (deType == "rand-2") {
      typeDE <- 1
      if (nPop < 6) { stop("population size must > 6 for DE/rand/2 algorithm.") }
    } 
    else if (deType == "best-1") {
      typeDE <- 2
      if (nPop < 3) { stop("population size must > 5 for DE/best/1 algorithm.") }
    } 
    else if (deType == "best-2") {
      typeDE <- 3
      if (nPop < 5) { stop("population size must > 5 for DE/best/2 algorithm.") }
    } 
    else if (deType == "rand-to-best-1") {
      typeDE <- 4
      if (nPop < 3) { stop("population size must > 3 for DE/rand-to-best/1 algorithm.") }
    } 
    else if (deType == "rand-to-best-2") {
      typeDE <- 5
      if (nPop < 5) { stop("population size must > 5 for DE/rand-to-best/2 algorithm.") }
    } 
  }
  
  nLoop <- length(nPop)
  #if (length(dPop))     dPop     <- numeric(nLoop)
  #if (length(varUpper))   varUpper   <- matrix(0, nLoop)
  #if (length(varLower))   varLower   <- matrix(0, nLoop)
  #if (length(maxIter))    maxIter    <- rep(100   , nLoop)
  #if (length(checkConv) < nLoop)  checkConv  <- rep(checkConv, nLoop)
  if (length(typeDE) < nLoop)    typeDE    <- rep(typeDE, nLoop)
  if (length(freeRun) < nLoop)    freeRun    <- rep(freeRun, nLoop)
  if (length(tol) < nLoop)        tol        <- rep(tol, nLoop)
  if (length(sf) < nLoop)         sf         <- rep(sf, nLoop)
  if (length(cr) < nLoop)         cr         <- rep(cr, nLoop)
  
  list(
    nPop = nPop, dPop = "autogen", varUpper = "autogen", varLower = "autogen", initPop = "autogen", 
    hasInitPop = "autogen", fixedDims = "autogen", maxIter = maxIter, 
    typeDE = typeDE, #checkConv = checkConv,
    freeRun = freeRun, tol = tol, sf = sf, cr = cr
  )
}