R/whv.R

In moocore: Core Mathematical Functions for Multi-Objective Optimization

#' Compute (total) weighted hypervolume given a set of rectangles
#'
#' Calculates the hypervolume weighted by a set of rectangles (with zero weight
#' outside the rectangles). The function [total_whv_rect()] calculates the
#' total weighted hypervolume as [hypervolume()]` + scalefactor *
#' abs(prod(reference - ideal)) * whv_rect()`. The details of the computation
#' are given by \citet{DiaLop2020ejor}.
#'
#' @inherit hypervolume params return
#'
#' @param rectangles `matrix()`\cr Weighted rectangles that will bias the
#'   computation of the hypervolume. Maybe generated by [eafdiff()] with
#'   `rectangles=TRUE` or by [choose_eafdiff()].
#'
#' @details
#'   TODO
#'
#' @seealso [read_datasets()], [eafdiff()], [choose_eafdiff()], [whv_hype()]
#'
#' @examples
#' rectangles <- as.matrix(read.table(header=FALSE, text='
#'  1.0  3.0  2.0  Inf    1
#'  2.0  3.5  2.5  Inf    2
#'  2.0  3.0  3.0  3.5    3
#' '))
#' whv_rect (matrix(2, ncol=2), rectangles, reference = 6)
#' whv_rect (matrix(c(2, 1), ncol=2), rectangles, reference = 6)
#' whv_rect (matrix(c(1, 2), ncol=2), rectangles, reference = 6)
#'
#' @references
#' \insertAllCited{}
#'
#'@concept metrics
#'@export
whv_rect <- function(x, rectangles, reference, maximise = FALSE)
{
  x <- as_double_matrix(x)
  nobjs <- ncol(x)
  if (nobjs != 2L) stop("sorry: only 2 objectives supported")
  if (ncol(rectangles) != 5L) stop("rectangles: invalid number of columns")
  if (is.null(reference)) stop("reference cannot be NULL")
  if (length(reference) == 1L) reference <- rep_len(reference, nobjs)
  # FIXME: This code does not handle maximisation yet.
  stopifnot(maximise == FALSE)
  ## if (any(maximise)) {
  ##   x <- transform_maximise(x, maximise)
  ##   if (length(maximise) == 1L) {
  ##     reference <- -reference
  ##     rectangles[,1L:4L] <- -rectangles[,1L:4L, drop = FALSE]
  ##   } else {
  ##     reference[maximise] <- -reference[maximise]
  ##     pos <- as.vector(matrix(1L:4L, nrow=2L)[, maximise, drop=FALSE])
  ##     rectangles[,pos] <- -rectangles[,pos, drop = FALSE]
  ##   }
  ## }
  .Call(rect_weighted_hv2d_C,
    t(x),
    t(rectangles),
    reference)
}


#' @inheritParams largest_eafdiff
#'
#' @param scalefactor `numeric(1)`\cr Real value within \eqn{(0,1]} that scales
#'   the overall weight of the differences. This is parameter psi (\eqn{\psi})
#'   in \citet{DiaLop2020ejor}.
#'
#' @examples
#' total_whv_rect (matrix(2, ncol=2), rectangles, reference = 6, ideal = c(1,1))
#' total_whv_rect (matrix(c(2, 1), ncol=2), rectangles, reference = 6, ideal = c(1,1))
#' total_whv_rect (matrix(c(1, 2), ncol=2), rectangles, reference = 6, ideal = c(1,1))
#'
#' @rdname whv_rect
#'
#' @concept metrics
#' @export
total_whv_rect <- function(x, rectangles, reference, maximise = FALSE, ideal = NULL, scalefactor = 0.1)
{
  x <- as.matrix(x)
  nobjs <- ncol(x)
  maximise <- as.logical(rep_len(maximise, nobjs))
  if (nobjs != 2L) stop("sorry: only 2 objectives supported")
  if (ncol(rectangles) != 5L) stop("invalid number of columns in rectangles (should be 5)")
  if (scalefactor <= 0 || scalefactor > 1) stop("scalefactor must be within (0,1]")

  hv <- hypervolume(x, reference, maximise = maximise)
  whv <- whv_rect(x, rectangles, reference, maximise = maximise)
  if (is.null(ideal)) {
    # FIXME: Should we include the range of the rectangles here?
    ideal <- get_ideal(x, maximise = maximise)
  }
  if (length(ideal) != nobjs) {
    stop("'ideal' should have same length as ncol(x)")
  }
  beta <- scalefactor * abs(prod(reference - ideal))
  #cat("beta: ", beta, "\n")
  hv + beta * whv
}

#' Approximation of the (weighted) hypervolume by Monte-Carlo sampling (2D only)
#'
#' Return an estimation of the hypervolume of the space dominated by the input
#' data following the procedure described by \citet{AugBadBroZit2009gecco}. A
#' weight distribution describing user preferences may be specified.
#'
#' @inheritParams hypervolume
#'
#' @param ideal `numeric()`\cr Ideal point as a vector of numerical values.
#'
#' @param nsamples `integer(1)`\cr Number of samples for Monte-Carlo sampling.
#'
#' @param seed `integer(1)`\cr Random seed.
#'
#' @param dist `character(1)`\cr Weight distribution type. See Details.
#'
#' @param mu `numeric()`\cr Parameter of the weight distribution. See Details.
#'
#' @details
#' The current implementation only supports 2 objectives.
#'
#' A weight distribution  \citep{AugBadBroZit2009gecco} can be provided via the `dist` argument. The ones currently supported are:
#'  * `"uniform"` corresponds to the default hypervolume (unweighted).
#'  * `"point"` describes a goal in the objective space, where the parameter `mu` gives the coordinates of the goal. The resulting weight distribution is a multivariate normal distribution centred at the goal.
#' * `"exponential"` describes an exponential distribution with rate parameter `1/mu`, i.e., \eqn{\lambda = \frac{1}{\mu}}.
#'
#' @return A single numerical value.
#'
#' @references
#' \insertAllCited{}
#'
#' @seealso    [read_datasets()], [eafdiff()], [whv_rect()]
#'
#' @doctest
#' @expect equal(3.99807)
#' whv_hype(matrix(2, ncol=2), reference = 4, ideal = 1, seed = 42)
#' @expect equal(3.00555)
#' whv_hype(matrix(c(3,1), ncol=2), reference = 4, ideal = 1, seed = 42)
#' @expect equal(1.14624)
#' whv_hype(matrix(2, ncol=2), reference = 4, ideal = 1, seed = 42,
#'          dist = "exponential", mu=0.2)
#' @expect equal(1.66815)
#' whv_hype(matrix(c(3,1), ncol=2), reference = 4, ideal = 1, seed = 42,
#'          dist = "exponential", mu=0.2)
#' @expect equal(0.64485)
#' whv_hype(matrix(2, ncol=2), reference = 4, ideal = 1, seed = 42,
#'          dist = "point", mu=c(2.9,0.9))
#' @expect equal(4.03632)
#' whv_hype(matrix(c(3,1), ncol=2), reference = 4, ideal = 1, seed = 42,
#'          dist = "point", mu=c(2.9,0.9))
#' @concept metrics
#' @export
whv_hype <- function(x, reference, ideal, maximise = FALSE,
                     nsamples = 1e5L, seed = NULL, dist = "uniform", mu = NULL)
{
  x <- as_double_matrix(x)
  nobjs <- ncol(x)

  if (is.null(reference)) stop("reference cannot be NULL")
  if (length(reference) == 1L) reference <- rep_len(reference, nobjs)
  if (is.null(ideal)) stop("ideal cannot be NULL")
  if (length(ideal) == 1L) ideal <- rep_len(ideal, nobjs)
  if (nobjs != 2L) stop("sorry: only 2 objectives supported")

  if (any(maximise)) {
    if (all(maximise)) {
      x <- -x
      reference <- -reference
      ideal <- -ideal
    } else if (length(maximise) != nobjs) {
      stop("length of maximise must be either 1 or ncol(x)")
    } else {
      x[,maximise] <- -x[,maximise]
      reference[maximise] <- -reference[maximise]
      ideal[maximise] <- -ideal[maximise]
    }
  }

  seed <- if (is.null(seed)) get_seed() else as_integer(seed)

  if (dist != "uniform" && is.null(mu)) {
    stop("parameter 'mu' is required when dist ='", dist, "'")
  }

  .Call(whv_hype_C,
    t(x),
    as.double(ideal),
    as.double(reference),
    as.integer(nsamples),
    dist,
    seed,
    mu)
}

get_seed <- function() sample.int(.Machine$integer.max, 1)