Nothing
#' @title Perform Nondominated Sorting
#'
#' @description
#' Assign elements of `fitnesses` to nondominated fronts.
#'
#' The first nondominated front is the set of individuals that is not dominated by any other
#' individual with respect to any fitness dimension, i.e. where no other individual exists that
#' has all fitness values greater or equal, with at least one fitness value strictly greater.
#'
#' The n'th nondominated front is the set of individuals that is not dominated by any other
#' individual that is not in any nondominated front with smaller n.
#'
#' Fitnesses are *maximized*, so the individuals in lower numbered nondominated fronts tend
#' to have higher fitness values.
#'
#' @template param_fitnesses
#' @param epsilon (`numeric`)\cr
#' Epsilon-vaue for non-dominance. A value is epsilon-dominated by another if it is at least `epsilon` smaller than
#' the other in all dimensions, and more than `epsilon` smaller than the other in one dimension. `epsilon` may
#' be a scalar, in which case it is used for all dimensions or a vector, in which case its length must match
#' the number of dimensions. Default 0.
#'
#' @return `list`: `$front`: Vector assigning each individual in `fitnesses` its nondominated front.
#' `$domcount`: Length N vector counting the number of individuals that dominate the given individual.
#'
#' @export
rank_nondominated = function(fitnesses, epsilon = 0) {
# this may or may not be similar to https://core.ac.uk/download/pdf/30341871.pdf ; haven't read the paper,
# but at first glance it looks like what I'm doing here so I don't think I am doing something dumb.
# They suggest binary search for fronts, which is probably a neat idea.
# For large matrices this seems to be a little faster than ecr's C code
assert_matrix(fitnesses, mode = "numeric", any.missing = FALSE, min.cols = 1, min.rows = 1)
assert(check_number(epsilon, lower = 0, finite = TRUE), check_numeric(epsilon, lower = 0, len = ncol(fitnesses), finite = TRUE, any.missing = FALSE))
fitnesses_t = t(fitnesses)
ordering = do.call(order, c(as.data.frame(unname(fitnesses)), list(decreasing = TRUE))) # ordering along first dimension, tie break 2nd dimension, tie break .. nth dimension
fronts = list() # matrices of nondominated fronts. The matrices here are transposed for more efficient comparison (so they have ncol(fitnesses) rows)
frontindices = list() # row indices of elements of nondominated fronts
assignedfronts = numeric(length(ordering)) # assigned front to each individual
domcount = numeric(length(ordering)) # number of individuals that dominate individuals
lastdomcount = 1 # number of individuals from the last considered front that dominate this individual ; init to 1 so in first loop we create a front
for (orow in ordering) { # go through ordered rows, with highest according to 1st dimension first
numdomd = 0 # number of individuals that dominate this individual
rowvals = fitnesses_t[, orow]
epsrowvals = rowvals + epsilon
for (ifront in seq_along(fronts)) { # go from first to last front, check if individual is nondominated there
# dominated: *all* fitnesses greater or equal than, at least one actually greater than
# The following line is the "hotspot" of this function and makes up ~90% of the runtime
# 'colSums(...) > 0' is used in place of 'any()'
lastdomcount = sum(colSums(fronts[[ifront]] < epsrowvals) == 0 & colSums(fronts[[ifront]] > epsrowvals) > 0)
numdomd = numdomd + lastdomcount
if (!lastdomcount) break
}
if (lastdomcount) {
# dominated by all fronts --> opening a new front
ifront = length(fronts) + 1
fronts[[ifront]] = matrix(rowvals)
frontindices[[ifront]] = orow
} else {
# nondominated by front <ifront> --> appending to this front
curfi = frontindices[[ifront]] # indices of individuals in curfi
curfi[[length(curfi) + 1]] = orow # add current index to it
frontindices[[ifront]] = curfi # save changed list
fronts[[ifront]] = fitnesses_t[, curfi, drop = FALSE] # also create the matrix
}
assignedfronts[[orow]] = ifront
domcount[[orow]] = numdomd
}
list(fronts = assignedfronts, domcount = domcount)
}
#' @title Calculate Crowding Distance
#'
#' @description
#' Takes a `matrix` of fitness values and calculates the crowding distance for individuals in that `matrix`.
#'
#' Individuals that are minimal or maximal with respect to at least one dimension are assigned infinite
#' crowding distance.
#'
#' Individuals are assumed to be in a (epsilon-) nondominated front.
#'
#' @template param_fitnesses
#'
#' @return `numeric`: Vector of crowding distances.
#'
#' @export
dist_crowding = function(fitnesses) {
assert_matrix(fitnesses, mode = "numeric", any.missing = FALSE, min.cols = 1, min.rows = 1)
result = numeric(nrow(fitnesses))
for (col in as.data.frame(unname(fitnesses))) {
o = order(col)
dists = diff(col[o])
result[o] = result[o] + c(Inf, dists) + c(dists, Inf)
}
result
}
#' @title Calculate Hypervolume Contribution
#'
#' @description
#' Takes a `matrix` of fitness values and calculates the hypervolume contributions of individuals in that `matrix`.
#'
#' Hypervolume contribution of an individual I is the difference between the dominated hypervolume of a set of
#' individuals including I, where the fitness of I is increased by `epsilon`, and the dominated hypervolume of
#' the same set but excluding I.
#'
#' Individuals that are less than another individual more than `epsilon` in any dimension have hypervolume contribution
#' of 0.
#'
#' @template param_fitnesses
#' @template param_nadir
#' @param epsilon (`numeric`)\cr
#' Added to each individual before calculating its particular hypervolume contribution. `epsilon` may
#' be a scalar, in which case it is used for all dimensions, or a vector, in which case its length must match
#' the number of dimensions. Default 0.
#' @return `numeric`: The vector of dominated hypervolume contributions for each individual in `fitnesses`.
#'
#' @examples
#' (fitnesses = matrix(c(1, 5, 2, 3, 0, 3, 1, 0, 10, 8), ncol = 2))
#'
#' # to see the fitness matrix, use:
#' ## plot(fitnesses, pch = as.character(1:5))
#'
#' domhv_contribution(fitnesses)
#' @export
domhv_contribution = function(fitnesses, nadir = 0, epsilon = 0) {
assert_matrix(fitnesses, mode = "numeric", any.missing = FALSE, min.cols = 1, min.rows = 1)
assert(check_number(nadir, lower = 0, finite = TRUE), check_numeric(epsilon, lower = 0, len = ncol(fitnesses), finite = TRUE, any.missing = FALSE))
assert(check_number(epsilon, lower = 0, finite = TRUE), check_numeric(epsilon, lower = 0, len = ncol(fitnesses), finite = TRUE, any.missing = FALSE))
has_epsilon = any(epsilon > 0)
nd = nondominated(fitnesses, epsilon = epsilon)
nonzeroes = nd$front
strongfront = fitnesses[nd$strong_front, , drop = FALSE]
result = numeric(nrow(fitnesses))
volume_baseline = domhv(strongfront, nadir, prefilter = FALSE)
result[nonzeroes] = map_dbl(nonzeroes, function(nz) {
replace = which(nz == nd$strong_front)
if (has_epsilon) {
if (length(replace)) {
strongfront[replace, ] = strongfront[replace, ] + epsilon
} else {
strongfront = rbind(strongfront, fitnesses[nz, ] + epsilon)
}
volume_with = domhv(strongfront, nadir, prefilter = FALSE)
} else {
volume_with = volume_baseline
}
if (length(replace)) {
volume_without = domhv(strongfront[-replace, , drop = FALSE], nadir, prefilter = FALSE)
} else {
volume_without = volume_baseline
}
volume_with - volume_without
})
result
}
#' @title Calculate Hypervolume Improvement
#'
#' @description
#' Takes a `matrix` of fitness values and calculates the hypervolume improvement of individuals in that `matrix`, one by one,
#' over the `baseline` individuals.
#'
#' The hypervolume improvement for each point is the measure of all points that have fitnesses that are
#' * greater than the respective value in `nadir` in all dimensions, and
#' * smaller than the respective value in the given point in all dimensions, and
#' * greater than all points in `baseline` in at least one dimension.
#'
#' Individuals in `fitnesses` are considered independently of each other. A possible speedup is achieved because
#' `baseline` individuals only need to be pre-filtered once.
#'
#' @template param_fitnesses
#' @param baseline (`matrix` | `NULL`)\cr
#' Fitness-matrix with one column per objective, giving a population over which the hypervolume improvement should be calculated.
#' If `NULL`, the hypervolume of each individual in `fitnesses` is calculated.
#' @template param_nadir
#' @return `numeric`: The vector of dominated hypervolume contributions for each individual in `fitnesses`.
#' @examples
#' (fitnesses = matrix(c(1, 5, 2, 3, 0, 3, 1, 0, 10, 8), ncol = 2))
#'
#' # to see the fitness matrix, use:
#' ## plot(fitnesses, pch = as.character(1:5))
#'
#' domhv_improvement(fitnesses)
#'
#' domhv_improvement(fitnesses, fitnesses[1, , drop = FALSE])
#' @export
domhv_improvement = function(fitnesses, baseline = NULL, nadir = 0) {
assert_matrix(fitnesses, mode = "numeric", any.missing = FALSE, min.cols = 1, min.rows = 1)
assert(check_number(nadir, finite = TRUE), check_numeric(nadir, finite = TRUE, any.missing = FALSE, len = ncol(fitnesses)))
if (is.null(baseline) || test_matrix(baseline, mode = "numeric", any.missing = FALSE, ncols = ncol(fitnesses), nrows = 0)) {
return(apply(fitnesses, 1, function(fi) {
prod(pmax(fi - nadir, 0))
}))
}
assert_matrix(baseline, mode = "numeric", any.missing = FALSE, ncols = ncol(fitnesses))
baseline = baseline[nondominated(baseline)$strong_front, , drop = FALSE]
blt = t(baseline)
blhv = domhv(baseline, nadir = nadir, prefilter = FALSE)
apply(fitnesses, 1, function(fi) {
# no HV improvement if:
# * there is an individual in baseline (i.e. blt) that has all values >= the fitness individual
# * fi is worse or equal nadir in any dimension
if (any(colSums(blt < fi) == 0) || any(nadir >= fi)) {
return(0)
}
domhv(rbind(baseline, fi), nadir = nadir, prefilter = FALSE) - blhv
})
}
#' @title Calculate Dominated Hypervolume
#'
#' @description
#' Use Chan's algorithm (`r format_bib("chan2013klee")`) to calculate dominated hypervolume.
#'
#' @template param_fitnesses
#' @template param_nadir
#' @param prefilter (`logical(1)`)\cr
#' Whether to make a first pass that filters out dominated individuals.
#' If it can be guaranteed that all individuals are non-dominated, setting this to `FALSE` improves performance a bit.
#' Otherwise the recommended value is the default `FALSE`.
#' @param on_worse_than_nadir (`character(1)`)
#' Action when individuals that do not dominate the nadir are found. One of `"quiet"` (ignore), `"warn"` (give warning, default), or `"stop"` (throw error).
#' @return `numeric(1)`: The dominated hypervolume of individuals in `fitnesses`.
#' @examples
#' (fitnesses = matrix(c(1, 5, 2, 3, 0, 3, 1, 0, 10, 8), ncol = 2))
#'
#' # to see the fitness matrix, use:
#' ## plot(fitnesses, pch = as.character(1:5))
#'
#' domhv(fitnesses)
#' @export
domhv = function(fitnesses, nadir = 0, prefilter = TRUE, on_worse_than_nadir = "warn") {
assert_matrix(fitnesses, mode = "numeric", any.missing = FALSE, min.cols = 1, min.rows = 1)
dim = ncol(fitnesses)
assert(check_number(nadir, finite = TRUE), check_numeric(nadir, len = dim, any.missing = FALSE, finite = TRUE))
assert_choice(on_worse_than_nadir, c("quiet", "warn", "stop"))
if (any(t(fitnesses) < nadir)) {
switch(on_worse_than_nadir, quiet = identity, warn = warning, stop = stop)("Found fitness worse than nadir")
fitnesses = fitnesses[colSums(t(fitnesses) < nadir) == 0, , drop = FALSE]
}
if (prefilter) {
fitnesses = fitnesses[nondominated(fitnesses)$strong_front, , drop = FALSE]
# cat("prefiltered:\n")
# print(fitnesses)
}
fitnesses_t = t(fitnesses)
zenith = apply(fitnesses, 2, max)
if (length(nadir) == 1) nadir = rep(nadir, dim)
weight_lut = 2^(seq_len(dim) / dim)
weight_lut[[1]] = 0 # first dimension is not regarded for weight
weight_lut = c(weight_lut, weight_lut)
domhv_recurse = function(fitnesses_t, nadir, zenith, dimension) {
# cat(sprintf("\nFrom %s to %s dim %s:\n", str_collapse(round(nadir, 4)), str_collapse(round(zenith, 4)), dimension))
# print(round(t(fitnesses_t), 4))
above = colSums(fitnesses_t >= zenith)
if (any(above == dim)) {
# one point completely dominates the current window --> no volume
# cat("completely dominated, returning 0\n")
return(0)
}
cutting = above == dim - 1L
cutcube = fitnesses_t[, cutting, drop = FALSE]
# see if we need to do "simplify". if no -> done with this loop
if (length(cutcube)) {
# the points that go beyond zenith stay at nadir (because they are not what we care about)
# but the points that do not indicate how much we cut away. From them we take the max, because
# there could be multiple orthants that cut away space.
# (We always cut away from below, because we have (-inf)^d-orthants)
cutcube[cutcube >= zenith] = -Inf
nadir = pmax(nadir, matrixStats::rowMaxs(cutcube)) # relatively lean dependency
#> nadir = pmax(nadir, apply(cutcube, 1L, max)) # base R
#> nadir = pmax(nadir, Rfast::rowMaxs(cutcube, value = TRUE))
fitnesses_t = fitnesses_t[, !cutting, drop = FALSE]
# for each simplify step, we need to check again if things fall out of the current window
below = colSums(fitnesses_t <= nadir) != 0
#> below = Rfast::colAny(fitnesses_t <= nadir) # doesn't seem to add much
# left out below the nadir --> no influence on current window
if (any(below)) {
fitnesses_t = fitnesses_t[, !below, drop = FALSE]
}
if (!length(fitnesses_t)) {
# cat("nothing left, returning full\n")
return(prod(zenith - nadir))
}
}
# one point left: subtract its volume from entire volume
if (ncol(fitnesses_t) == 1L) {
# cat("one point left, returning its volume\n")
return(prod(zenith - nadir) - prod(pmin(zenith, fitnesses_t) - nadir)) # I don't like what this does to numerics
}
repeat {
dimpoints = fitnesses_t[dimension, ]
which.cutpoints = which(dimpoints < zenith[dimension])
if (length(which.cutpoints)) break
dimension = (dimension %% dim) + 1L
}
cutpoints = dimpoints[which.cutpoints]
# Maybe a different cut is more appropriate here, because more points fall out on the upper half.
# - points "below" are eliminated 100% points "above" are eliminated only if they have one dimension left
# - on the other hand the majority case could be that most points have only one dimension overlap
#> cutat = Rfast::nth(cutpoints, (length(cutpoints) + 1L) / 2) # not the correct way of doing this, because it doesn't really count the (d-2)-faces *within* the bounds.
# interpolate = TRUE in the following because we never want to get the maximum element.
if (length(cutpoints) == 1) {
cutat = cutpoints
} else if (length(cutpoints) == 2) {
cutat = min(cutpoints)
} else {
# calculate weights; disregard dimension itself by subtracting its own weight again at the end.
cutpointweights = colSums((fitnesses_t[, which.cutpoints, drop = FALSE] < zenith) * weight_lut[seq.int(dim + 1 - dimension, length.out = dim)])
cutat = matrixStats::weightedMedian(cutpoints, cutpointweights, interpolate = FALSE, ties = "min") # n log(n), apparently does sorting internally, pathetic.
if (cutat == max(cutpoints)) {
otherpoints = cutpoints[cutpoints != cutat]
if (length(otherpoints)) {
cutat = max(otherpoints)
}
}
}
#> cutat = sort(cutpoints)[(length(cutpoints) + 1L) / 2]
## alternative: sort partially
#> cutrank = as.integer((length(cutpoints) + 1L) / 2)
#> cutat = sort(cutpoints, partial = cutrank)[[cutrank]]
cutnadir = nadir
cutnadir[dimension] = cutat
cutzenith = zenith
cutzenith[dimension] = cutat
dimension = (dimension %% dim) + 1L
# pre-emptively drop lower part for upper half
fitnesses_upper = fitnesses_t[, dimpoints > cutat, drop = FALSE]
if (length(fitnesses_upper)) {
# cat("recursing upper\n")
result_upper = domhv_recurse(fitnesses_upper, cutnadir, zenith, dimension)
} else {
result_upper = prod(zenith - cutnadir)
}
# cat("recursing lower\n")
result_upper + domhv_recurse(fitnesses_t, nadir, cutzenith, dimension)
}
# problem seems to be upper of first lower: doesn't exist for clash.
prod(zenith - nadir) - domhv_recurse(fitnesses_t, nadir, zenith, 1)
}
# similar to rank_nondominated, but return only the first front.
# @return named `list`: `front` (`integer`): indices of epsilon-nondominated front; `strong_front` (`integer`): indices of dominating front (without epsilon)
nondominated = function(fitnesses, epsilon = 0) {
assert_matrix(fitnesses, mode = "numeric", any.missing = FALSE, min.cols = 1, min.rows = 1)
assert(check_number(epsilon, lower = 0, finite = TRUE), check_numeric(epsilon, lower = 0, len = ncol(fitnesses), any.missing = FALSE, finite = TRUE))
fitnesses_t = t(fitnesses)
ordering = do.call(order, c(as.data.frame(unname(fitnesses)), list(decreasing = TRUE))) # ordering along first dimension, tie break 2nd dimension, tie break .. nth dimension
curfi = ordering[[1]] # row indices of elements of nondominated fronts
curfi_strong = ordering[[1]] # front of dominating individuals
front = fitnesses_t[, curfi, drop = FALSE] # matrix of nondominated front, for quick comparison
for (orow in ordering[-1]) { # go through ordered rows, with highest according to 1st dimension first. (But skip the actual highest because we already added it.)
rowvals = fitnesses_t[, orow]
epsrowvals = rowvals + epsilon
# dominated: *all* fitnesses greater or equal than, at least one actually greater than
# The following line is the "hotspot" of this function and makes up ~90% of the runtime
# 'colSums(...) > 0' is used in place of 'any()'
dominated = any(colSums(front < epsrowvals) == 0 & colSums(front > epsrowvals) > 0)
if (!dominated) {
# strong dominating front
if (all(colSums(front <= rowvals) > 0)) {
curfi_strong[[length(curfi_strong) + 1]] = orow
front = fitnesses_t[, curfi_strong, drop = FALSE] # also create the matrix
}
# nondominated by front --> appending to front
curfi[[length(curfi) + 1]] = orow # add current index to it
}
}
list(front = curfi, strong_front = curfi_strong)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.