R/variation_diffmut.R

In fcampelo/MOEADr: Component-Wise MOEA/D Implementation

#' Differential Mutation
#'
#' Differential Mutation implementation for the MOEA/D
#'
#' This function generalizes many variations of the Differential Mutation
#' operator with general form:
#'
#' `u = x_basis + Phi(x_a - x_b)`
#'
#' Where u is the new candidate vector, `Phi != 0` is a real number,
#' and `x_basis`, `x_a` and `x_b` are distinct vectors.
#'
#' This routine is intended to be used internally by [perform_variation()],
#' and should not be called directly by the user.
#'
#' @param X Population matrix
#' @param P Matrix of selection probabilities (generated by
#' [define_neighborhood()])
#' @param B Matrix of neighborhoods (generated by [define_neighborhood()])
#' @param Phi Mutation parameter. Either a scalar numeric constant, or NULL for
#'            randomly chosen between `0` and `1` (independently sampled for
#'            each operation).
#' @param basis how to select the basis vector. Currently supported methods are:
#'        \itemize{
#'          \item `basis = "rand"`, for using a randomly sampled vector from the
#'                 population;
#'          \item `basis = "mean"`, for using the mean point of the
#'                 neighborhood;
#'          \item `basis = "wgi"`, for using the the weighted mean point of the
#'                 neighborhood.
#'        }
#' @param ... other parameters to be passed down to specific options of basis
#' vector generation (e.g., `Y`, `Yt`, `W`, `scaling` and `aggfun`, required
#' when `basis = "wgi"`).
#'
#' @return Matrix `X`' containing the mutated population
#'
#' @section References:
#' K. Price, R.M. Storn, J.A. Lampinen, "Differential Evolution: A
#' Practical Approach to Global Optimization", Springer 2005\cr
#'
#' F. Campelo, L.S. Batista, C. Aranha (2020): The {MOEADr} Package: A
#' Component-Based Framework for Multiobjective Evolutionary Algorithms Based on
#' Decomposition. Journal of Statistical Software \doi{10.18637/jss.v092.i06}\cr
#'
#' D. V. Arnold, “Weighted multirecombination evolution strategies,”
#' Theoretical Computer Science 361(1):18–37, 2006.
#'
#' @export

variation_diffmut <- function(X, P, B, Phi = NULL, basis = 'rand', ...){

  input.pars <- as.list(sys.call())[-1]

  # ========== Error catching and default value definitions
  assertthat::assert_that(
    is.numeric(X) && is.matrix(X),
    is.numeric(P) && is.matrix(P) && is_within(P, 0, 1, strict = FALSE),
    identical(nrow(X), nrow(P)),
    nrow(P) == ncol(P),
    is.numeric(B) && is.matrix(B),
    nrow(B) == nrow(X),
    is.null(Phi) || (is.numeric(Phi) && Phi != 0),
    is.element(basis, c('rand', 'mean', 'wgi')))
  # ==========

  dimX <- dim(X)
  # Generate replacement indexes for xbasis, x0, x1
  # (Basis is recreated if 'mean' or 'wgi')
  R <- t(sapply(1:dimX[1],
                FUN = function(i) {
                  sample.int(dimX[1],
                             size    = 3,
                             replace = FALSE,
                             prob    = P[, i]) }))

  if (is.null(Phi)) {
    Phi <- matrix(stats::runif(dimX[1]),
                  nrow  = dimX[1],
                  ncol  = dimX[2],
                  byrow = FALSE)
  }

  if (basis == "rand"){
    Xb <- X[R[, 1], , drop = FALSE]
  } else if (basis == "mean"){
    Xb <- t(sapply(1:nrow(X),
                   FUN = function(i) {
                     apply(X[B[i, , drop = FALSE], ],
                           MARGIN = 2,
                           FUN    = mean) }))
  } else if (basis == "wgi"){
    # Calculate scalarized function values for each neighborhood
    normYs <- scale_objectives(Y       = input.pars$Y,
                               Yt      = input.pars$Yt,
                               scaling = input.pars$scaling)
    bigZ   <- scalarize_values(normYs  = normYs,
                               W       = input.pars$W,
                               B       = B,
                               aggfun  = input.pars$aggfun)

    # Remove the last row, which refers to x_i^{(t-1)}
    bigZ <- t(bigZ[-nrow(bigZ), , drop = FALSE])

    # Calculate weights for WGI
    wgi.W <- log(ncol(bigZ) + 0.5) - log(1:ncol(bigZ))

    Xb <- t(sapply(1:dimX[1],
                   FUN = function(i){
                     indx <- order(bigZ[i, , drop = FALSE])
                     colSums(wgi.W * X[B[i, indx, drop = FALSE], ])
                   }))
  }

  # Perform mutations and return
  return(Xb + Phi * (X[R[, 2], , drop = FALSE] - X[R[, 3], , drop = FALSE]))
}