R/updt_best.R

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

#' Best Neighborhood Replacement Update for MOEA/D
#'
#' Population update using the best neighborhood replacement method for the
#' MOEADr package.
#'
#' The Best Neighborhood Replacement method consists of three steps:
#'
#'  - For each subproblem `i`, the best candidate solution `x_j` from the
#'           entire population is determined.
#'  - The neighborhood of subproblem `i` is replaced by the neighborhood
#'     of subproblem j. The size of this neighborhood is given by a parameter
#'     `Tr`.
#'  - The Restricted replacement (see [updt_restricted()]) is then
#'     applied using this new neighborhood.
#'
#' This update routine is intended to be used internally by the main [moead()]
#' function, and should not be called directly by the user.
#'
#'
#' @param update List containing the population update parameters. See
#' Section `Update Strategies` of the [moead()] documentation for
#' details. `update` must have the following key-value pairs:
#' - `update$Tr`: positive integer, neighborhood size for the update
#'   operation
#' - `update$nr`: positive integer, maximum number of copies of a given
#'   candidate solution.
#' @param X Matrix of candidate solutions
#' @param Xt Matrix of incumbent solutions
#' @param Y Matrix of objective function values of `X`
#' @param Yt Matrix of objective function values of `Xt`
#' @param V List object containing information about the constraint violations
#' of the candidate solutions, generated by [evaluate_population()]
#' @param Vt List object containing information about the constraint violations
#' of the incumbent solutions, generated by [evaluate_population()]
#' @param normYs List generated by [scale_objectives()], containing two matrices
#' of scaled objective values (`normYs$Y` and `normYs$Yt`) and two vectors,
#' containing the current estimates of the ideal (`normYs$minP`) and nadir
#' (`normYs$maxP`) points. See [scale_objectives()] for details.
#' @param W matrix of weights, generated by [generate_weights()].
#' @param aggfun List containing the aggregation function parameters. See
#' Section `Scalar Aggregation Functions` of the [moead()] documentation for
#' details.
#' @param BP Neighborhood list, generated by [define_neighborhood()].
#' @param constraint list containing the parameters defining the constraint
#' handling method. See Section `Constraint Handling` of the [moead()]
#' documentation for details.
#' @param ... other parameters (included for compatibility with generic call)
#'
#' @return List object containing the update population matrix (`X`),
#' and its corresponding matrix of objective function values (`Y`) and
#' constraint value list (`V`).
#'
#' @export
#'
#' @section References:
#' 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

updt_best <- function(update, X, Xt, Y, Yt, V, Vt,
                      normYs, W, BP, constraint, aggfun, ...){

  ## Verify that the necessary parameters exist.
  assertthat::assert_that(
    all(assertthat::has_name(update, c("nr", "Tr"))),
    assertthat::is.count(update$nr),
    assertthat::is.count(update$Tr))

  nr <- update$nr
  Tr <- update$Tr

  # Calculate scalarized performance of all individuals for all subproblems
  fullZ   <- scalarize_values(normYs = normYs,
                              W      = W,
                              B      = BP$fullB,
                              aggfun = aggfun)

  # Find the problem in which each CANDIDATE solution (not incumbent) performs
  # best
  best.indx <- apply(X      = fullZ[1:(nrow(fullZ) - 1), , drop = FALSE],
                     MARGIN = 1,
                     FUN    = which.min)

  best.subprob <- mapply(FUN      = function(i, j, B){B[i, j]},
                         i        = 1:nrow(BP$fullB),
                         j        = best.indx,
                         MoreArgs = list(B = BP$fullB))

  # Define restricted neighborhoods for best update (that is, the update
  # neighborhood of subproblem i is set as the neighborhood of best.subprob[i])
  bestB    <- BP$fullB[best.subprob, 1:Tr, drop = FALSE]

  # Assemble bigZ matrix according to neighborhood bestB
  bestZ <- scalarize_values(normYs = normYs,
                            W      = W,
                            B      = bestB,
                            aggfun = aggfun)

  best.sel.indx <- order_neighborhood(bigZ       = bestZ,
                                      B          = bestB,
                                      V          = V,
                                      Vt         = Vt,
                                      constraint = constraint)

  # ========= Code below here should be identical to updt_restricted =========#

  # Function for returning the selected solution (variable or objectives space)
  # for a subproblem:
  # - i: subproblem index
  # - sel.indx: matrix of selection indices
  # - XY: matrix of candidate solutions (in variable or objective space)
  # - XYt: matrix of incumbent solutions (in variable or objective space)
  # - B: matrix of neighborhoods
  do.update <- function(i, sel.indx, XY, XYt, B){
    for (j in sel.indx[i,]) {               #each element in b_i, in fitness order
      if (j > ncol(B)) return(XYt[i, , drop = FALSE])     # last row = incumbent solution
      else if (used[B[i, j]] < nr)          # tests if the current element is still available
      {
        used[B[i, j]] <<- used[B[i, j]] + 1 # modifies count matrix in parent env
        return(XY[B[i, j], , drop = FALSE])
      }
    }
  }

  # Vector of indices (random permutation), and deshuffling vector
  I  <- sample.int(nrow(X))
  I2 <- order(I)

  # Counter of how many time each solution has been used
  used <- rep(0, nrow(X))

  # Update matrix of candidate solutions
  Xnext <- t(vapply(X         = I,
                    FUN       = do.update,
                    FUN.VALUE = numeric(ncol(X)),
                    sel.indx  = best.sel.indx,
                    XY        = X,
                    XYt       = Xt,
                    B         = bestB,
                    USE.NAMES = FALSE))
  Xnext <- Xnext[I2, ]

  # Update matrix of function values
  used  <- rep(0, nrow(Y))
  Ynext <- t(vapply(X         = I,
                    FUN       = do.update,
                    FUN.VALUE = numeric(ncol(Y)),
                    sel.indx  = best.sel.indx,
                    XY        = Y,
                    XYt       = Yt,
                    B         = bestB,
                    USE.NAMES = FALSE))
  Ynext <- Ynext[I2, ]

  # Update list of constraint values
  if(is.null(V)){
    Vnext <- NULL
  } else{
    Vnext <- list(Cmatrix = NULL, Vmatrix = NULL, v = NULL)

    ## 1: Cmatrix
    used <- rep(0, nrow(Y))
    Vnext$Cmatrix <- t(vapply(X         = I,
                              FUN       = do.update,
                              FUN.VALUE = numeric(ncol(V$Cmatrix)),
                              sel.indx  = best.sel.indx,
                              XY        = V$Cmatrix,
                              XYt       = Vt$Cmatrix,
                              B         = bestB,
                              USE.NAMES = FALSE))
    ## 2: Vmatrix
    used <- rep(0, nrow(Y))
    Vnext$Vmatrix <- t(vapply(X         = I,
                              FUN       = do.update,
                              FUN.VALUE = numeric(ncol(V$Vmatrix)),
                              sel.indx  = best.sel.indx,
                              XY        = V$Vmatrix,
                              XYt       = Vt$Vmatrix,
                              B         = bestB,
                              USE.NAMES = FALSE))

    ## 3: v
    Vnext$v <- rowSums(Vnext$Vmatrix)
  }

  # Output
  return(list(X = Xnext,
              Y = Ynext,
              V = Vnext))
}