R/estimate_probabilities.R

In PrzeChoj/gips: Gaussian Model Invariant by Permutation Symmetry

#' Estimate posterior probabilities from Metropolis-Hastings run
#'
#' We use the "second approach" from the paper.
#'
#' @param perms A list of `gips_perm` objects. Visited groups during the Metropolis-Hastings run.
#' @param show_progress_bar A boolean. Indicate whether or not to show the progress bar.
#'
#' @returns A named numeric vector. Names: character representations of permutations.
#' Elements: estimated posterior probabilities of permutations.
#' @noRd
estimate_probabilities <- function(perms, show_progress_bar = FALSE) {
  if (show_progress_bar) {
    progressBar <- utils::txtProgressBar(min = 0, max = length(perms), initial = 1)
  }
  group_representatives <- sapply(1:length(perms), function(i) {
    if (show_progress_bar) {
      utils::setTxtProgressBar(progressBar, i)
    }
    as.character(get_group_representative(perms[[i]]))
  })
  repr_counts <- table(group_representatives)
  if (show_progress_bar) {
    close(progressBar)
  }

  repr_weights <- sapply(names(repr_counts), function(p_str) {
    perm <- permutations::char2cycle(p_str)
    p_order <- permutations::permorder(perm)
    1 / numbers::eulersPhi(p_order)
  })
  unnormalized_probabilities <- repr_counts * repr_weights / length(perms)
  probabilities <- unnormalized_probabilities / sum(unnormalized_probabilities)
  probabilities <- as.numeric(probabilities)
  names(probabilities) <- names(unnormalized_probabilities)

  sort(probabilities, decreasing = TRUE)
}


#' Get a representative of a cyclic permutation group
#'
#' Essentially a "nu" function from paper (beginning of section 4.1.1)
#' `get_representative(perm) == nu(< perm >)`, where "`nu`" is from the paper,
#' and "`< perm >`" is a cyclic group generated by the permutation.
#'
#' @param perm `gips_perm`
#'
#' @returns Object of a `gips_perm` class.
#' @noRd
get_group_representative <- function(perm) {
  size <- attr(perm, "size")
  if (size == 0) {
    return(perm)
  }
  perm <- permutations::as.cycle(perm)
  p_order <- permutations::permorder(perm)
  coprimes <- get_coprimes(p_order)
  all_perms <- lapply(coprimes, function(cp) {
    permutations::cycle_power(perm, cp)
  })
  if (length(all_perms) == 0) {
    return(gips_perm(permutations::nullword, 0))
  }
  all_perms_as_vectors <- lapply(all_perms, function(p) {
    as.integer(permutations::cycle2word(p))
  })
  all_perms_as_strings <- sapply(all_perms_as_vectors, function(v) {
    paste(v, collapse = " ")
  })
  min_index <- stringi::stri_order(all_perms_as_strings)[1]
  gips_perm(all_perms[[min_index]], size)
}

#' Get coprime numbers
#'
#' @param n A single integer.
#'
#' @returns All integers smaller than n, that are coprime to n. Exception: n = 1,
#' in which case 1 is returned.
#' @noRd
get_coprimes <- function(n) {
  if (n == 1) {
    return(1)
  }
  smaller_ints <- 1:(n - 1)
  are_coprime <- sapply(smaller_ints, function(m) numbers::coprime(n, m))
  smaller_ints[are_coprime]
}


#' Calculate exact posterior probabilities
#'
#' We use the "second approach" from the paper.
#'
#' @param perms An output of `permutations::allperms()`.
#' @param log_posteriories A vector of all values of log posteriories of all `perms`.
#' @param show_progress_bar A boolean. Indicate whether or not to show two progress bars.
#'
#' @returns A named numeric vector. Names: character representations of permutations.
#' Elements: estimated posterior probabilities of permutations.
#' @noRd
calculate_probabilities <- function(perms, log_posteriories, show_progress_bar = FALSE) {
  if (show_progress_bar) {
    progressBar <- utils::txtProgressBar(min = 0, max = length(perms), initial = 1)
  }

  for (i in 1:19) {
    perms_size <- i
    if (OEIS_A051625[i] == length(perms) || OEIS_A000142[i] == length(perms)) {
      break
    }
  }

  if (perms_size == 19) {
    rlang::abort("There is sth wrong with this sequence of permutations!",
      "i" = "The length of the permutation vector has to be an element of OEIS sequence A051625 or A000142",
      "x" = paste0("You have the length of permutation vector = ", length(perms))
    )
  }

  group_representatives <- character(0)
  for (i in 1:length(perms)) {
    if (show_progress_bar) {
      utils::setTxtProgressBar(progressBar, i)
    }

    # We should use `g_perm <- gips_perm(perms[i], perms_size)`,
    # but this is faster, because it lacks safe checks:
    g_perm <- gips_perm_no_checks(perms[i], perms_size)

    group_representatives[i] <- as.character(get_group_representative(g_perm))
  }
  
  if (show_progress_bar) {
    close(progressBar)
  }


  # get rid of the repeated permutations:
  if (show_progress_bar) {
    progressBar <- utils::txtProgressBar(min = 0, max = length(group_representatives), initial = 1)
  }
  groups_unrepeated_log_posteriories <- log_posteriories[1]
  names(groups_unrepeated_log_posteriories)[1] <- group_representatives[1]
  for (i in 2:length(group_representatives)) {
    if (show_progress_bar) {
      utils::setTxtProgressBar(progressBar, i)
    }
    if (!(group_representatives[i] %in% names(groups_unrepeated_log_posteriories))) {
      groups_unrepeated_log_posteriories[length(groups_unrepeated_log_posteriories) + 1] <- log_posteriories[i]
      names(groups_unrepeated_log_posteriories)[length(groups_unrepeated_log_posteriories)] <- group_representatives[i]
    }
  }

  if (show_progress_bar) {
    close(progressBar)
  }

  sort(change_log_probabilities_unnorm_to_probabilities(groups_unrepeated_log_posteriories), decreasing = TRUE)
}