simulate_gaussian <- function(n, mean, cov) {
  z <- matrix(rnorm(n * length(mean)), ncol = length(mean))
  sweep(z %*% chol(cov), 2, mean, "+")
}

all_permutations <- function(x) {
  if (length(x) <= 1L) {
    return(list(x))
  }

  out <- vector("list", factorial(length(x)))
  idx <- 1L
  for (i in seq_along(x)) {
    for (perm in all_permutations(x[-i])) {
      out[[idx]] <- c(x[i], perm)
      idx <- idx + 1L
    }
  }
  out
}

normalize_labels <- function(x) {
  lev <- sort(unique(x))
  vapply(x, function(value) which(lev == value), integer(1))
}

best_label_accuracy <- function(truth, pred) {
  truth <- normalize_labels(truth)
  pred <- normalize_labels(pred)
  k <- max(c(truth, pred))
  perms <- all_permutations(seq_len(k))

  max(vapply(
    perms,
    function(mapping) mean(truth == mapping[pred]),
    numeric(1)
  ))
}