From 405b467882fcf5d9a60fbb171aef69679646111c Mon Sep 17 00:00:00 2001
From: jd_bunker <jbunker7@outlook.com>
Date: Mon, 20 Jul 2026 21:22:51 -0500
Subject: [PATCH 1/2] Overhaul allocate() - 20Jul2026

Overhauled allocate(), especially regarding the rounding and lbound/N.h/n.samp-based adjustments
* `outputs` Parameter: allows user to choose between raw, adjusted, and/or final (rounded) allocations
* c.h Check: no longer listing Neyman
* sum(lbound * c.h) > cost: bumped from warning to error after edits to adjustment/rounding step
* any(lbound > N.h): gives warning
* `outputs` check: every value must be 'raw', 'adjusted', and/or 'rounded'
* Adjustment Step: properly routed depending on allocation method
* Rounding: stochastic to deterministic
* Documentation: updated to reflect all of the above changes
* Tests: patched existing tests based on above changes (e.g., warning to error, adjustment/rounding handling); added a few new tests
---
 R/allocate.R                   | 1151 ++++++++++++++++++++++----------
 man/allocate.Rd                |   59 +-
 tests/testthat/test-allocate.R | 1071 ++++++++++++++++-------------
 3 files changed, 1401 insertions(+), 880 deletions(-)

diff --git a/R/allocate.R b/R/allocate.R
index c0936fa..b0bbc53 100644
--- a/R/allocate.R
+++ b/R/allocate.R
@@ -1,370 +1,781 @@
-#' Sample allocation
-#'
-#' Compute the proportional, power, Neyman, and optimal sample allocations.
-#'
-#' @section Method:
-#' The \emph{allocate} function allocates a sample size \emph{n} on \emph{H} strata using one of the following allocation methods:
-#' \enumerate{
-#'   \item Proportional allocation \[\code{n.samp, N.h, allocation = "proportional"}\]
-#'         \deqn{n_h = n \times \frac{N_h}{\sum\limits_{h=1}^H N_h}}
-#'         where \cr
-#'         \eqn{n}: total sample size to be allocated (function input is \code{n.samp}), and \cr
-#'         \eqn{N_h}: population size of stratum \emph{h} (function input is \code{N.h}).
-#'   \item Power allocation \[\code{n.samp, N.h, power, allocation = "power"}\]
-#'         \deqn{n_h = n \times \frac{N_h^\alpha}{\sum\limits_{h=1}^H N_h^\alpha}}
-#'          where \cr
-#'          \eqn{\alpha}: a power value to control over-under-sampling with \eqn{0 \le \alpha \le 1} (function input is \code{power}).
-#'   \item Neyman allocation \[\code{n.samp, N.h, S.h, allocation = "neyman"}\]
-#'         \deqn{n_h = n \times \frac{N_h S_h}{\sum\limits_{h=1}^H N_h S_h}}
-#'         where \cr
-#'         \eqn{S_h}:  standard deviation of stratum \emph{h} (function input is \code{S.h}).
-#'    \item Optimal allocation
-#'           \itemize{
-#'             \item cost-constrained \[\code{N.h, S.h, c.h, cost, allocation = "optimal"}\]
-#'                   \deqn{n_h = (C-c_0) \times \frac{N_h S_h / \sqrt{c_h}}{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}}
-#'             where \cr
-#'             \eqn{c_h}: cost per unit in stratum \emph{h} (function input is \code{c.h}), and \cr
-#'             \eqn{(C - c_0)}: total variable cost (function input is \code{cost})
-#'             \item precision-constrained \[\code{N.h, S.h, c.h, variance, allocation = "optimal"}\]
-#'                   \deqn{n_h = N_h S_h / \sqrt{c_h} \times \frac{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}{V_0 \left(\sum\limits_{h=1}^H  N_h \right)^2 + \sum\limits_{h=1}^H N_h S_h^2}}
-#'             where \cr
-#'             \eqn{V_0}: fixed variance target for estimated mean (function input is \code{variance})
-#'           }
-#' }
-#'
-#' The table below presents the relevant inputs for each type; when irrelevant inputs are entered, an error message will be displayed.
-#'
-#' \tabular{lllllllll}{
-#'   \strong{allocation} \tab \strong{N.h} \tab \strong{n.samp} \tab \strong{S.h} \tab \strong{c.h} \tab \strong{cost} \tab \strong{variance} \tab \strong{lbound} \tab \strong{power} \cr
-#'   proportional                   \tab X \tab X \tab  \tab  \tab  \tab  \tab X \tab \cr
-#'   power                          \tab X \tab X \tab  \tab  \tab  \tab  \tab X \tab X\cr
-#'   neyman                         \tab X \tab X \tab X \tab X \tab  \tab  \tab X \tab \cr
-#'   optimal: cost-constrained      \tab X \tab  \tab X \tab X \tab X \tab  \tab X \tab \cr
-#'   optimal: precision-constrained \tab X \tab  \tab X \tab X \tab  \tab X \tab X \tab
-#' }
-#'
-#' @param allocation type of allocation, must be one of \code{"proportional"}, \code{"power"}, \code{"neyman"}, or \code{"optimal"}.
-#' @param n.samp total sample size to be allocated (positive integer of length 1). \cr\cr
-#'        required for the following allocation types: proportional, power, and Neyman, and \code{NULL} otherwise.
-#' @param N.h vector of population stratum sizes (\eqn{N_h}, all positive values), for example \code{c(150, 600, 250)}. \cr\cr
-#' required for all allocation types.
-#' @param S.h vector of stratum unit standard deviations (positive values same length as \code{N.h}) (\eqn{S_h}). \cr\cr
-#' required for the following allocation types: Neyman, and optimal, and \code{NULL} otherwise.
-#' @param c.h vector of cost per unit in stratum h (positive values same length as \code{N.h}) (\eqn{c_h}). \cr\cr
-#' required for the optimal allocation only, and \code{NULL} otherwise.
-#' @param cost total variable cost (positive value) \eqn{(C - c_0)}. \cr\cr
-#' required for the cost-constrained optimal allocation only, and \code{NULL} otherwise.
-#' @param variance fixed variance target for estimated mean (positive value) (\eqn{V_0}). \cr\cr
-#' required for the precision-constrained optimal allocation only, and \code{NULL} otherwise.
-#' @param power power value for power allocation (\eqn{0 \le \alpha \le 1}). \cr\cr
-#' required for the power allocation only, and \code{NULL} otherwise.
-#' @param lbound minimum stratum-level (positive integer of length 1). Default value is 2.
-#' @return Integer vector same length of N.h final allocation
-#' @export
-#' @examples
-#' # The first step is getting a frame summary
-#' #  Summarize the IPEDS dataset by OBEREG
-#' # - N: number of universities per region
-#' # - SD_ENRTOT: standard deviation of total enrollment per region
-#' # - Filter out rows with missing ENRTOT to ensure accurate variance estimates
-#'
-#' ipeds_summary <- ipeds |>
-#'   tidytable::filter(!is.na(ENRTOT)) |>
-#'   tidytable::group_by(OBEREG) |>
-#'   tidytable::summarize(
-#'     N = tidytable::n(),
-#'     SD_ENRTOT = stats::sd(ENRTOT)
-#'   ) |>
-#'   tidytable::ungroup()
-#'
-#' # Example of proportional allocation
-#' ipeds_summary |>
-#'   tidytable::mutate(
-#'     n = allocate("proportional", N.h = N, n.samp = 500)
-#'   )
-#'
-#' # Example of power allocation
-#' ipeds_summary |>
-#'   tidytable::mutate(
-#'     n = allocate("power", N.h = N, power = 0.5, n.samp = 500)
-#'   )
-#'
-#' # Example of Neyman allocation
-#' ipeds_summary |>
-#'   tidytable::mutate(
-#'     n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT)
-#'   )
-#'
-#' # Example of Neyman allocation with a lower bound of 5
-#' ipeds_summary |>
-#'   tidytable::mutate(
-#'     n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT, lbound = 5)
-#'   )
-allocate <- function(allocation, N.h, n.samp = NULL, S.h = NULL, c.h = NULL, cost = NULL, variance = NULL, power = NULL, lbound = 2) {
-  allocation <- match.arg(allocation, c("proportional", "power", "neyman", "optimal"))
-
-  ######
-  # Check inputs
-  .problems <- NULL # Initialize list of problems found with inputs
-  .addProblem <- function(parameter, condition, problems = .problems, allocation = NULL) { # Function to simplify addition to problems found with inputs to our running list (.problems)
-    if (is.null(parameter)) { # No parameter given
-      problem <- condition
-    } else if (parameter %in% c("allocation", "n.samp", "N.h", "S.h", "c.h", "cost", "variance", "power", "lbound")) {
-      problem <- paste0("The ", parameter, " parameter ", condition)
-    }
-    problemsNew <- c(problems, problem)
-    return(problemsNew)
-  }
-
-  .condition <- paste0('must be specified for allocation=="', allocation, '"')
-  # n.samp parameter
-  if (allocation %in% c("proportional", "power", "neyman") & is.null(n.samp)) {
-    .problems <- .addProblem(parameter = "n.samp", condition = .condition)
-  } else if (!allocation %in% c("proportional", "power", "neyman") & !is.null(n.samp)) {
-    warning("The n.samp parameter should be NULL", call. = FALSE)
-  }
-  # N.h parameter
-  if (allocation %in% c("proportional", "power", "neyman", "optimal") & is.null(N.h)) {
-    .problems <- .addProblem(parameter = "N.h", condition = .condition)
-  }
-  # S.h parameter
-  if (allocation %in% c("neyman", "optimal") & is.null(S.h)) {
-    .problems <- .addProblem(parameter = "S.h", condition = .condition)
-  } else if (!allocation %in% c("neyman", "optimal") & !is.null(S.h)) {
-    warning("The S.h parameter should be NULL", call. = FALSE)
-  }
-  # c.h parameter
-  if (allocation %in% c("optimal") & is.null(c.h)) {
-    .problems <- .addProblem(parameter = "c.h", condition = .condition)
-  } else if (!allocation %in% c("optimal") & !is.null(c.h)) {
-    warning("The c.h parameter should be NULL", call. = FALSE)
-  }
-  # power parameter
-  if (allocation %in% c("power") & is.null(power)) {
-    .problems <- .addProblem(parameter = "power", condition = .condition)
-  } else if (!allocation %in% c("power") & !is.null(power)) {
-    warning("The power parameter should be NULL", call. = FALSE)
-  }
-  # optimal allocation: only 1 of cost or variance should be provided
-  if (allocation == "optimal" & sum(is.null(cost), is.null(variance)) != 1) {
-    .problems <- .addProblem(parameter = NULL, condition = paste0('Exactly one of the cost and variance parameters should be supplied for allocation=="', allocation, '"'))
-  } else if (!allocation %in% c("optimal") & !(is.null(cost) & is.null(variance))) {
-    if (!is.null(cost)) {
-      warning("The cost parameter should be NULL", call. = FALSE)
-    }
-    if (!is.null(variance)) {
-      warning("The variance parameter should be NULL", call. = FALSE)
-    }
-  }
-
-  ###
-  # Check for miscellaneous unexpected parameter values
-
-  # n.samp
-  .condition <- "must be a positive integer of length 1"
-  if (!is.null(n.samp) & allocation %in% c("proportional", "power", "neyman")) {
-    if (!all(length(n.samp) == 1 & (typeof(n.samp) %in% c("integer") | (typeof(n.samp) == "double" & round(n.samp) == n.samp)) & n.samp > 0)) {
-      .problems <- .addProblem(parameter = "n.samp", condition = .condition)
-    }
-  }
-  # lbound
-  if (!(length(lbound) == 1 & (typeof(lbound) %in% c("integer") | (typeof(lbound) == "double" & round(lbound) == lbound)) & lbound > 0)) {
-    .problems <- .addProblem(parameter = "lbound", condition = .condition)
-  }
-
-  .condition <- "must be a vector of positive values (integers or non-integers)"
-  # N.h
-  if (!is.null(N.h)) {
-    if (!(length(N.h) >= 1 & typeof(N.h) %in% c("integer", "double") & all(N.h > 0))) {
-      .problems <- .addProblem(parameter = "N.h", condition = .condition)
-    }
-  }
-  .condition <- paste0(.condition, " that are the same length as N.h")
-  # S.h
-  if (!is.null(S.h) & allocation %in% c("neyman", "optimal")) {
-    if (!(length(S.h) >= 1 & typeof(S.h) %in% c("integer", "double") & all(S.h > 0) & length(S.h) == length(N.h))) {
-      .problems <- .addProblem(parameter = "S.h", condition = .condition)
-    }
-  }
-  # c.h
-  if (!is.null(c.h) & allocation %in% c("neyman", "optimal")) {
-    if (!(length(c.h) >= 1 & typeof(c.h) %in% c("integer", "double") & all(c.h > 0) & length(c.h) == length(N.h))) {
-      .problems <- .addProblem(parameter = "c.h", condition = .condition)
-    }
-  }
-  # cost
-  .condition <- "must be a positive value (integer or non-integer)"
-  if (!is.null(cost) & allocation %in% "optimal") {
-    if (!all(length(cost) == 1 & typeof(cost) %in% c("integer", "double") & cost > 0)) {
-      .problems <- .addProblem(parameter = "cost", condition = .condition)
-    }
-  }
-  # variance
-  if (allocation == "optimal" & !is.null(variance)) {
-    if (!all(length(variance) == 1 & typeof(variance) %in% c("integer", "double") & variance > 0 & length(variance) == 1)) {
-      .problems <- .addProblem(parameter = "variance", condition = .condition)
-    }
-  }
-
-  # power parameter
-  .condition <- "must be a positive value between 0 and 1, inclusive"
-  if (allocation == "power" & !is.null(power)) {
-    if (!(length(power) == 1 & typeof(power) %in% c("integer", "double") & 0 <= power & power <= 1)) {
-      .problems <- .addProblem(parameter = "power", condition = .condition)
-    }
-  }
-
-
-  if (allocation %in% c("proportional", "power", "neyman")) {
-    if (!(is.null(lbound) | is.null(N.h) | is.null(n.samp))) {
-      if (!all(lbound * length(N.h) <= n.samp)) {
-        .problems <- c(.problems, "lbound*length(N.h) must be less than or equal to n.samp")
-      }
-    }
-  }
-
-  if (!(is.null(N.h) | is.null(n.samp))) {
-    if (length(.problems) == 0) {
-      if (sum(N.h) < n.samp) {
-        warning("sum(N.h) is less than n.samp")
-      }
-    }
-  }
-
-  if (allocation %in% c("optimal")) {
-    if (!is.null(cost) & length(.problems) == 0) {
-      if (sum(lbound * c.h) > cost) {
-        warning("No solution exists for specified total cost (sum(lbound*c.h) > cost)")
-      }
-    }
-  }
-  ###
-  # Aggregate problems and stop the program if necessary
-  if (length(.problems) > 0) {
-    if (length(.problems) > 1) {
-      .problems <- paste0(
-        1:length(.problems),
-        ": ",
-        .problems
-      )
-    }
-    stop(
-      "\n",
-      paste0(.problems, collapse = "\n")
-    )
-  }
-
-  ######
-  # Moving onto the actual allocation
-  N <- sum(N.h)
-  if (allocation == "proportional") {
-    allocations <- n.samp * N.h / N
-  } else if (allocation == "power") {
-    N.h.powered <- N.h**power
-    N.powered <- sum(N.h.powered)
-    allocations <- n.samp * N.h.powered / N.powered
-  } else if (allocation == "neyman") {
-    propNum <- N.h * S.h # Numerator
-    propDen <- sum(propNum) # Denominator
-    allocations <- n.samp * propNum / propDen
-  } else if (allocation == "optimal") {
-    if (!is.null(cost)) { # Cost-constrained
-      propNum <- N.h * S.h / sqrt(c.h)
-      propDen <- sum(N.h * S.h * sqrt(c.h))
-      allocations <- cost * propNum / propDen
-    } else if (!is.null(variance)) { # Precision-constrained
-      propNum <- sum(N.h * S.h * sqrt(c.h))
-      propDen <- variance * sum(N.h)**2 + sum(N.h * S.h**2)
-      allocations <- N.h * S.h / sqrt(c.h) * propNum / propDen
-    }
-  }
-  # Calculate the total (raw) sample size
-
-
-  final_n <- max(ceiling(sum(allocations)), length(N.h) * lbound)
-
-
-  sizes <- allocations
-
-  num_groups <- length(N.h)
-
-  # Ensure each allocation is at least lbound
-  adjusted_allocations <- pmax(allocations, lbound)
-
-  # Calculate the sum of the adjusted allocations
-  total_allocated <- sum(adjusted_allocations)
-
-  # Calculate the difference from the total n
-  difference <- final_n - total_allocated
-
-  # If difference is positive, distribute it proportionally
-  if (difference != 0) {
-    remaining_sizes <- sizes
-    remaining_sizes <- adjusted_allocations - lbound
-    remaining_total <- sum(remaining_sizes)
-
-    if (remaining_total > 0) {
-      additional_allocations <- difference * remaining_sizes / remaining_total
-    } else {
-      additional_allocations <- rep(0, num_groups)
-    }
-
-    adjusted_allocations <- adjusted_allocations + additional_allocations
-  }
-
-  # Round the adjusted allocations
-  rounded_allocations <- round(adjusted_allocations)
-
-  # Step 5: Adjust the rounded allocations to ensure the sum equals n
-  total_allocated <- sum(rounded_allocations)
-  difference <- final_n - total_allocated
-  # Adjust the allocations by adding/subtracting the difference
-  while (difference != 0) {
-    i <- sample(1:num_groups, 1) # Randomly select an index
-    if (difference > 0) {
-      rounded_allocations[i] <- rounded_allocations[i] + 1
-      difference <- difference - 1
-    } else {
-      if (rounded_allocations[i] > lbound) {
-        rounded_allocations[i] <- rounded_allocations[i] - 1
-        difference <- difference + 1
-      }
-    }
-  }
-
-  # Prep for outputting
-  if (allocation == "proportional") {
-    inputs <- list("N.h" = N.h)
-  } else if (allocation == "power") {
-    inputs <- list("N.h" = N.h, "power" = power)
-  } else if (allocation == "neyman") {
-    inputs <- list("N.h" = N.h, "S.h" = S.h)
-  } else if (allocation == "optimal" & !is.null(cost)) {
-    inputs <- list("N.h" = N.h, "S.h" = S.h, "c.h" = c.h, "cost" = cost)
-  } else if (allocation == "optimal" & !is.null(variance)) {
-    inputs <- list("N.h" = N.h, "S.h" = S.h, "c.h" = c.h, "variance" = variance)
-  } else {
-    inputs <- list(N.h, S.h, c.h, cost, variance, power)
-  }
-  outputs <- as.integer(rounded_allocations)
-  if (allocation == "optimal") {
-    n.print <- final_n
-  } else {
-    n.print <- n.samp
-  }
-  message(paste0("Sample allocation of ", n.print, " using ", allocation, " with the relevant inputs:"))
-  for (i in 1:length(inputs)) {
-    message(paste0("  ",
-      names(inputs)[i],
-      " = ",
-      paste0(inputs[[i]], collapse = ", "),
-      collapse = ""
-    ))
-  }
-  message()
-  message("Output:")
-  message(paste0(outputs,
-    collapse = ", "
-  ))
-  return(outputs)
-}
+#' Sample allocation
+#'
+#' Compute the proportional, power, Neyman, and optimal sample allocations.
+#'
+#' @section Method:
+#' The \emph{allocate} function allocates a sample size \emph{n} on \emph{H} strata using one of the following allocation methods:
+#' \enumerate{
+#'   \item Proportional allocation \[\code{n.samp, N.h, allocation = "proportional"}\]
+#'         \deqn{n_h = n \times \frac{N_h}{\sum\limits_{h=1}^H N_h}}
+#'         where \cr
+#'         \eqn{n}: total sample size to be allocated (function input is \code{n.samp}), and \cr
+#'         \eqn{N_h}: population size of stratum \emph{h} (function input is \code{N.h}).
+#'   \item Power allocation \[\code{n.samp, N.h, power, allocation = "power"}\]
+#'         \deqn{n_h = n \times \frac{N_h^\alpha}{\sum\limits_{h=1}^H N_h^\alpha}}
+#'          where \cr
+#'          \eqn{\alpha}: a power value to control over-under-sampling with \eqn{0 \le \alpha \le 1} (function input is \code{power}).
+#'   \item Neyman allocation \[\code{n.samp, N.h, S.h, allocation = "neyman"}\]
+#'         \deqn{n_h = n \times \frac{N_h S_h}{\sum\limits_{h=1}^H N_h S_h}}
+#'         where \cr
+#'         \eqn{S_h}:  standard deviation of stratum \emph{h} (function input is \code{S.h}).
+#'    \item Optimal allocation
+#'           \itemize{
+#'             \item cost-constrained \[\code{N.h, S.h, c.h, cost, allocation = "optimal"}\]
+#'                   \deqn{n_h = (C−c_0) \times \frac{N_h S_h / \sqrt{c_h}}{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}}
+#'             where \cr
+#'             \eqn{c_h}: cost per unit in stratum \emph{h} (function input is \code{c.h}), and \cr
+#'             \eqn{(C – c_0)}: total variable cost (function input is \code{cost})
+#'             \item precision-constrained \[\code{N.h, S.h, c.h, variance, allocation = "optimal"}\]
+#'                   \deqn{n_h = N_h S_h / \sqrt{c_h} \times \frac{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}{V_0 \left(\sum\limits_{h=1}^H  N_h \right)^2 + \sum\limits_{h=1}^H N_h S_h^2}}
+#'             where \cr
+#'             \eqn{V_0}: fixed variance target for estimated mean (function input is \code{variance})
+#'           }
+#' }
+#'
+#' The table below presents the relevant inputs for each type; when irrelevant inputs are entered, an error message will be displayed.
+#'
+#' \tabular{lllllllll}{
+#'   \strong{allocation} \tab \strong{N.h} \tab \strong{n.samp} \tab \strong{S.h} \tab \strong{c.h} \tab \strong{cost} \tab \strong{variance} \tab \strong{lbound} \tab \strong{power} \cr
+#'   proportional                   \tab ✓ \tab ✓ \tab  \tab  \tab  \tab  \tab ✓ \tab \cr
+#'   power                          \tab ✓ \tab ✓ \tab  \tab  \tab  \tab  \tab ✓ \tab ✓\cr
+#'   neyman                         \tab ✓ \tab ✓ \tab ✓ \tab ✓ \tab  \tab  \tab ✓ \tab \cr
+#'   optimal: cost-constrained      \tab ✓ \tab  \tab ✓ \tab ✓ \tab ✓ \tab  \tab ✓ \tab \cr
+#'   optimal: precision-constrained \tab ✓ \tab  \tab ✓ \tab ✓ \tab  \tab ✓ \tab ✓ \tab
+#' }
+#'
+#' @param allocation type of allocation, must be one of \code{"proportional"}, \code{"power"}, \code{"neyman"}, or \code{"optimal"}.
+#' @param n.samp total sample size to be allocated (positive integer of length 1). \cr\cr
+#'        required for the following allocation types: proportional, power, and Neyman, and \code{NULL} otherwise.
+#' @param N.h vector of population stratum sizes (\eqn{N_h}, all positive values), for example \code{c(150, 600, 250)}. \cr\cr
+#' required for all allocation types.
+#' @param S.h vector of stratum unit standard deviations (positive values same length as \code{N.h}) (\eqn{S_h}). \cr\cr
+#' required for the following allocation types: Neyman, and optimal, and \code{NULL} otherwise.
+#' @param c.h vector of cost per unit in stratum h (positive values same length as \code{N.h}) (\eqn{c_h}). \cr\cr
+#' required for the optimal allocation only, and \code{NULL} otherwise.
+#' @param cost total variable cost (positive value) \eqn{(C – c_0)}. \cr\cr
+#' required for the cost-constrained optimal allocation only, and \code{NULL} otherwise.
+#' @param variance fixed variance target for estimated mean (positive value) (\eqn{V_0}). \cr\cr
+#' required for the precision-constrained optimal allocation only, and \code{NULL} otherwise.
+#' @param power power value for power allocation (\eqn{0 \le \alpha \le 1}). \cr\cr
+#' required for the power allocation only, and \code{NULL} otherwise.
+#' @param lbound minimum stratum-level (positive integer of length 1). Default value is 2.
+#' @param outputs
+#' character vector representing whether to output:\cr
+#'  \enumerate{
+#'    \item the raw allocations before accounting for N.h, lbound, and n.samp if needed \[\code{"raw"}\],
+#'    \item the continuous version after accounting for the above \[\code{"adjusted"}\], and/or
+#'    \item the rounded version of the above \[\code{"rounded"}\]
+#'  } \cr
+#'  Default is to only return the final rounded version \[\code{"rounded"}\]. \cr
+#'  If one version is requested, the result will be a numeric vector. Otherwise, the result will be a named list matching the requested outputs.
+#' @return If one output type is requested, a numeric vector of allocations.
+#' If multiple output types are requested, a named list containing the requested
+#' allocation vectors.
+#' @export
+#' @examples
+#' # The first step is getting a frame summary
+#' #  Summarize the IPEDS dataset by OBEREG
+#' # - N: number of universities per region
+#' # - SD_ENRTOT: standard deviation of total enrollment per region
+#' # - Filter out rows with missing ENRTOT to ensure accurate variance estimates
+#'
+#' ipeds_summary <- ipeds |>
+#'   tidytable::filter(!is.na(ENRTOT)) |>
+#'   tidytable::group_by(OBEREG) |>
+#'   tidytable::summarize(
+#'     N = tidytable::n(),
+#'     SD_ENRTOT = stats::sd(ENRTOT)
+#'   ) |>
+#'   tidytable::ungroup()
+#'
+#' # Example of proportional allocation
+#' ipeds_summary |>
+#'   tidytable::mutate(
+#'     n = allocate("proportional", N.h = N, n.samp = 500)
+#'   )
+#'
+#' # Example of power allocation
+#' ipeds_summary |>
+#'   tidytable::mutate(
+#'     n = allocate("power", N.h = N, power = 0.5, n.samp = 500)
+#'   )
+#'
+#' # Example of Neyman allocation
+#' ipeds_summary |>
+#'   tidytable::mutate(
+#'     n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT)
+#'   )
+#'
+#' # Example of Neyman allocation with a lower bound of 5
+#' ipeds_summary |>
+#'   tidytable::mutate(
+#'     n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT, lbound = 5)
+#'   )
+
+allocate <- function(allocation, N.h, n.samp = NULL, S.h = NULL, c.h = NULL, cost = NULL, variance = NULL, power = NULL, lbound = 2, outputs = "rounded") {
+  allocation <- match.arg(allocation, c("proportional", "power", "neyman", "optimal"))
+  outputs <- match.arg(outputs, c("raw", "adjusted", "rounded"), several.ok = TRUE)
+
+  ######
+  # Check inputs
+  .problems <- NULL # Initialize list of problems found with inputs
+  .addProblem <- function(parameter, condition, problems = .problems, allocation = NULL) { # Function to simplify addition to problems found with inputs to our running list (.problems)
+    if (is.null(parameter)) { # No parameter given
+      problem <- condition
+    } else if (parameter %in% c("allocation", "n.samp", "N.h", "S.h", "c.h", "cost", "variance", "power", "lbound")) {
+      problem <- paste0("The ", parameter, " parameter ", condition)
+    }
+    problemsNew <- c(problems, problem)
+    return(problemsNew)
+  }
+
+  .condition <- paste0('must be specified for allocation=="', allocation, '"')
+  # n.samp parameter
+  if (allocation %in% c("proportional", "power", "neyman") & is.null(n.samp)) {
+    .problems <- .addProblem(parameter = "n.samp", condition = .condition)
+  } else if (!allocation %in% c("proportional", "power", "neyman") & !is.null(n.samp)) {
+    warning("The n.samp parameter should be NULL", call. = FALSE)
+  }
+  # N.h parameter
+  if (allocation %in% c("proportional", "power", "neyman", "optimal") & is.null(N.h)) {
+    .problems <- .addProblem(parameter = "N.h", condition = .condition)
+  }
+  # S.h parameter
+  if (allocation %in% c("neyman", "optimal") & is.null(S.h)) {
+    .problems <- .addProblem(parameter = "S.h", condition = .condition)
+  } else if (!allocation %in% c("neyman", "optimal") & !is.null(S.h)) {
+    warning("The S.h parameter should be NULL", call. = FALSE)
+  }
+  # c.h parameter
+  if (allocation %in% c("optimal") & is.null(c.h)) {
+    .problems <- .addProblem(parameter = "c.h", condition = .condition)
+  } else if (!allocation %in% c("optimal") & !is.null(c.h)) {
+    warning("The c.h parameter should be NULL", call. = FALSE)
+  }
+  # power parameter
+  if (allocation %in% c("power") & is.null(power)) {
+    .problems <- .addProblem(parameter = "power", condition = .condition)
+  } else if (!allocation %in% c("power") & !is.null(power)) {
+    warning("The power parameter should be NULL", call. = FALSE)
+  }
+  # optimal allocation: only 1 of cost or variance should be provided
+  if (allocation == "optimal" & sum(is.null(cost), is.null(variance)) != 1) {
+    .problems <- .addProblem(parameter = NULL, condition = paste0('Exactly one of the cost and variance parameters should be supplied for allocation=="', allocation, '"'))
+  } else if (!allocation %in% c("optimal") & !(is.null(cost) & is.null(variance))) {
+    if (!is.null(cost)) {
+      warning("The cost parameter should be NULL", call. = FALSE)
+    }
+    if (!is.null(variance)) {
+      warning("The variance parameter should be NULL", call. = FALSE)
+    }
+  }
+
+  ###
+  # Check for miscellaneous unexpected parameter values
+
+  # n.samp
+  .condition <- "must be a positive integer of length 1"
+  if (!is.null(n.samp) & allocation %in% c("proportional", "power", "neyman")) {
+    if (!all(length(n.samp) == 1 & (typeof(n.samp) %in% c("integer") | (typeof(n.samp) == "double" & round(n.samp) == n.samp)) & n.samp > 0)) {
+      .problems <- .addProblem(parameter = "n.samp", condition = .condition)
+    }
+  }
+  # lbound
+  if (!(length(lbound) == 1 & (typeof(lbound) %in% c("integer") | (typeof(lbound) == "double" & round(lbound) == lbound)) & lbound > 0)) {
+    .problems <- .addProblem(parameter = "lbound", condition = .condition)
+  }
+
+  .condition <- "must be a vector of positive values (integers or non-integers)"
+  # N.h
+  if (!is.null(N.h)) {
+    if (!(length(N.h) >= 1 & typeof(N.h) %in% c("integer", "double") & all(N.h > 0))) {
+      .problems <- .addProblem(parameter = "N.h", condition = .condition)
+    }
+  }
+  .condition <- paste0(.condition, " that are the same length as N.h")
+  # S.h
+  if (!is.null(S.h) & allocation %in% c("neyman", "optimal")) {
+    if (!(length(S.h) >= 1 & typeof(S.h) %in% c("integer", "double") & all(S.h > 0) & length(S.h) == length(N.h))) {
+      .problems <- .addProblem(parameter = "S.h", condition = .condition)
+    }
+  }
+  # c.h
+  if (!is.null(c.h) & allocation %in% c("optimal")) {
+    if (!(length(c.h) >= 1 & typeof(c.h) %in% c("integer", "double") & all(c.h > 0) & length(c.h) == length(N.h))) {
+      .problems <- .addProblem(parameter = "c.h", condition = .condition)
+    }
+  }
+  # cost
+  .condition <- "must be a positive value (integer or non-integer)"
+  if (!is.null(cost) & allocation %in% "optimal") {
+    if (!all(length(cost) == 1 & typeof(cost) %in% c("integer", "double") & cost > 0)) {
+      .problems <- .addProblem(parameter = "cost", condition = .condition)
+    }
+  }
+  # variance
+  if (allocation == "optimal" & !is.null(variance)) {
+    if (!all(length(variance) == 1 & typeof(variance) %in% c("integer", "double") & variance > 0 & length(variance) == 1)) {
+      .problems <- .addProblem(parameter = "variance", condition = .condition)
+    }
+  }
+
+  # power parameter
+  .condition <- "must be a positive value between 0 and 1, inclusive"
+  if (allocation == "power" & !is.null(power)) {
+    if (!(length(power) == 1 & typeof(power) %in% c("integer", "double") & 0 <= power & power <= 1)) {
+      .problems <- .addProblem(parameter = "power", condition = .condition)
+    }
+  }
+
+
+  if (allocation %in% c("proportional", "power", "neyman")) {
+    if (!(is.null(lbound) | is.null(N.h) | is.null(n.samp))) {
+      if (!all(lbound * length(N.h) <= n.samp)) {
+        .problems <- c(.problems, "lbound*length(N.h) must be less than or equal to n.samp")
+      }
+    }
+  }
+
+  if (!(is.null(N.h) | is.null(n.samp))) {
+    if (length(.problems) == 0) {
+      if (sum(N.h) < n.samp) {
+        warning("sum(N.h) is less than n.samp")
+      }
+    }
+  }
+
+  if (allocation %in% c("optimal")) {
+    if (!is.null(cost) & length(.problems) == 0) {
+      if (sum(lbound * c.h) > cost) {
+        .problems <- c(.problems, "sum(lbound*c.h) must be less than or equal to cost")
+      }
+    }
+  }
+
+  if (allocation %in% c("proportional", "power", "neyman","optimal")){
+    if (!is.null(N.h)){
+      if (any(lbound > N.h)){
+        warning("lbound > N.h for at least one stratum")
+      }
+    }
+  }
+
+  # outputs parameter
+  .condition <- 'must be one or more of: "raw", "adjusted", "rounded"'
+  if (!all(outputs %in% c("raw", "adjusted", "rounded"))) {
+    .problems <- .addProblem(parameter = "outputs", condition = .condition)
+  }
+
+  ###
+  # Aggregate problems and stop the program if necessary
+  if (length(.problems) > 0) {
+    if (length(.problems) > 1) {
+      .problems <- paste0(
+        1:length(.problems),
+        ": ",
+        .problems
+      )
+    }
+    stop(
+      "\n",
+      paste0(.problems, collapse = "\n")
+    )
+  }
+
+  ######
+  #Fixed-total continuous adjustment:
+  # * Respects lower bound (lbound)
+  # * Respects upper bound (N.h)
+  # * Sums to sample size (n.samp)
+  # * Unrounded sample sizes
+  .adjust_fixed_total <- function(weights, n.samp, lbound, N.h) {
+    h <- length(weights)
+
+    lbound.h <- rep(lbound, h) #Extending lbound to a vector
+
+    # Raw proportional target
+    raw <- n.samp * weights / sum(weights)
+
+    #If raw already respects both bounds, use it directly.
+    if (all(raw >= lbound.h) && all(raw <= N.h)) {
+      return(as.numeric(raw))
+    }
+
+    #Otherwise fall back to iterative capacity-and-floor enforcement
+    alloc <- lbound.h
+    remaining <- n.samp - sum(alloc)
+    capacity <- N.h - alloc #Amount left in N.h after allocating lbound
+
+    if (remaining == 0) {
+      return(as.numeric(alloc))
+    }
+
+    #Initialize strata variables
+    active <- rep(TRUE, h) #Does stratum still need to be processed
+    extra <- rep(0, h) #How much to allocate to stratum above raw_allocations
+
+
+    repeat {
+      w <- weights
+      w[!active] <- 0
+
+      if (sum(w) <= 0) {
+        stop("No feasible solution: cannot distribute remaining sample.")
+      }
+
+      proposed <- remaining * w / sum(w) #Proposed additional allocation
+      hit <- active & (proposed > capacity) #Assigned all of stratum
+
+      #If didn't hit any statum's upper bound, break here
+      if (!any(hit)) {
+        extra[active] <- remaining * weights[active] / sum(weights[active])
+        break
+      }
+
+      #If hit a stratum's N, allocate its remaining capacity and set inactive
+      alloc[hit] <- alloc[hit] + capacity[hit]
+      remaining <- n.samp - sum(alloc)
+      active[hit] <- FALSE
+      capacity[hit] <- 0
+    }
+
+
+    adjusted_allocations <- alloc + extra
+    adjusted_allocations
+  }
+
+  ######
+  # Calculate stratified mean variance
+
+  .stratified_mean_variance <- function(n.h, N.h, S.h, N = sum(N.h)) {
+    W.h <- N.h / N
+    sum(W.h^2 * (1 - n.h / N.h) * S.h^2 / n.h)
+  }
+
+  ######
+  # Precision-constrained continuous adjustment:
+  # * Respects lower bound (lbound)
+  # * Respects upper bound (N.h)
+  # * Hits variance target
+  # * Unrounded sample sizes
+
+  .adjust_precision_constrained <- function(N.h, S.h, c.h, variance, lbound){
+    h <- length(N.h)
+    N <- sum(N.h)
+
+    a.h <- N.h * S.h / sqrt(c.h)
+
+    lbound.h <- rep(lbound, h)
+
+    # Initialize strata variables
+    adjusted_allocations <- rep(NA_real_, h)
+    fixed <- rep(FALSE, h) #Treat stratum fixed (TRUE) vs. free (FALSE)
+
+
+    V.lower <- .stratified_mean_variance(lbound.h, N.h, S.h)
+    V.upper <- .stratified_mean_variance(N.h, N.h, S.h)
+
+    if (variance >= V.lower) {
+      warning("The lower-bound allocation already satisfies the variance target; returning lower-bound allocation.")
+      return(lbound.h)
+    }
+
+    if (variance < V.upper) {
+      stop("No feasible allocation: even a census of all strata does not achieve the variance target.")
+    }
+
+
+    repeat {
+      # Compute variance contribution from fixed strata
+      if (any(fixed)) {
+        V.fixed <- .stratified_mean_variance(
+          n.h = adjusted_allocations[fixed],
+          N.h = N.h[fixed],
+          S.h = S.h[fixed],
+          N = N
+        )
+      } else {
+        V.fixed <- 0
+      }
+
+      # If all strata are fixed, break
+      if (all(fixed)) {
+        V.total <- .stratified_mean_variance(adjusted_allocations, N.h, S.h)
+        if (abs(V.total - variance) > 1e-10) {
+          warning("No feasible allocation exactly meets the variance target under the imposed bounds.")
+        }
+        break
+      }
+
+      # Recompute scaling constant over free strata
+      free <- !fixed
+      A.free <- sum(N.h[free] * S.h[free] * sqrt(c.h[free]))
+      B.free <- sum(N.h[free] * S.h[free]^2)
+
+      denom <- N^2 * (variance - V.fixed) + B.free
+
+      if (denom <= 0) {
+        stop("No feasible allocation solution for the requested variance after applying bounds.")
+      }
+
+      k.free <- A.free / denom
+
+      # Candidate allocations for free strata
+      candidate <- k.free * a.h[free]
+
+      # Check for violations
+      below <- candidate < lbound.h[free]
+      above <- candidate > N.h[free]
+
+      if (!any(below | above)){
+        adjusted_allocations[free] <- candidate
+        break
+      }
+
+      # Fix violating strata
+      free_idx <- which(free)
+
+      if (any(below)) {
+        idx <- free_idx[below]
+        adjusted_allocations[idx] <- lbound.h[idx]
+        fixed[idx] <- TRUE
+      }
+
+      if (any(above)) {
+        idx <- free_idx[above]
+        adjusted_allocations[idx] <- N.h[idx]
+        fixed[idx] <- TRUE
+      }
+
+    }
+    return(adjusted_allocations)
+  }
+
+  ######
+  # Precision-constrained rounding
+
+  .round_precision_constrained <- function(adjusted_allocations, N.h, S.h, c.h, variance, lbound) {
+    rounded_allocations <- floor(adjusted_allocations + 1e-9)
+    rounded_allocations <- pmax(rounded_allocations, lbound)
+    rounded_allocations <- pmin(rounded_allocations, N.h)
+    rounded_allocations <- as.integer(rounded_allocations)
+
+    frac <- adjusted_allocations - floor(adjusted_allocations + 1e-9)
+
+    while (.stratified_mean_variance(rounded_allocations, N.h, S.h) > variance) {
+      candidates <- which(rounded_allocations < N.h)
+
+      if (length(candidates) == 0) {
+        warning("No feasible integer allocation found to satisfy variance target.")
+        break
+      }
+
+      current_var <- .stratified_mean_variance(rounded_allocations, N.h, S.h)
+      scores <- rep(-Inf, length(N.h))
+
+      for (h in candidates) {
+        trial <- rounded_allocations
+        trial[h] <- trial[h] + 1L
+        new_var <- .stratified_mean_variance(trial, N.h, S.h)
+        var_gain <- current_var - new_var
+
+        #leftover relative to adjusted continuous allocation
+        leftover <- pmax(N.h[h] - adjusted_allocations[h], 0)
+
+        #normalize leftover
+        leftover_factor <- 1 + leftover / N.h[h]
+
+        #fractional preference
+        frac_factor <- 1 + frac[h]
+
+        #combined score
+        scores[h] <- (var_gain / c.h[h]) * leftover_factor * frac_factor
+      }
+
+      best_h <- which.max(scores)
+      if (!is.finite(scores[best_h]) || scores[best_h] <= 0) {
+        warning("Could not improve enough to satisfy variance constraint.")
+        break
+      }
+
+      rounded_allocations[best_h] <- rounded_allocations[best_h] + 1L
+    }
+
+    return(rounded_allocations)
+  }
+
+
+  ######
+  # Rounding fixed-totals
+
+  .round_fixed_total <- function(adjusted_allocations, n.samp, N.h, lbound, score) {
+    rounded_allocations <- floor(adjusted_allocations + 1e-9)
+    rounded_allocations <- pmax(rounded_allocations, lbound)
+    rounded_allocations <- pmin(rounded_allocations, N.h)
+    rounded_allocations <- as.integer(rounded_allocations)
+
+    deficit <- n.samp - sum(rounded_allocations)
+
+    if (deficit < 0) {
+      stop("Initial rounded allocation exceeds the target total.")
+    }
+    while (deficit > 0) {
+      candidates <- which(rounded_allocations < N.h)
+      if (length(candidates) == 0) {
+        stop("No feasible way to allocate remaining units.")
+      }
+
+
+      scores <- score(rounded_allocations)
+      scores[-candidates] <- -Inf
+
+
+      j <- candidates[which.max(scores[candidates])]
+
+
+      if (!is.finite(scores[j])) {
+        stop("No feasible scoring candidate for additional allocation.")
+      }
+
+      rounded_allocations[j] <- rounded_allocations[j] + 1L
+      deficit <- deficit - 1L
+    }
+
+    return(rounded_allocations)
+
+  }
+
+  .score_proportional <- function(adjusted_allocations, rounded_allocations, N.h) {
+    gap <- pmax(adjusted_allocations - rounded_allocations, 0)
+    leftover <- pmax(N.h - adjusted_allocations, 0)
+    weight <- N.h / sum(N.h)
+
+    gap * (1 + leftover / N.h) * (1 + weight)
+  }
+
+  .score_power <- function(adjusted_allocations, rounded_allocations, N.h, power) {
+    gap <- pmax(adjusted_allocations - rounded_allocations, 0)
+    leftover <- pmax(N.h - adjusted_allocations, 0)
+    weight_raw <- N.h^power
+    weight <- weight_raw / sum(weight_raw)
+
+    gap * (1 + leftover / N.h) * (1 + weight)
+  }
+
+  .score_neyman <- function(adjusted_allocations, rounded_allocations, N.h, S.h) {
+    gap <- pmax(adjusted_allocations - rounded_allocations, 0)
+    leftover <- pmax(N.h - adjusted_allocations, 0)
+    weight_raw <- N.h * S.h
+    weight <- weight_raw / sum(weight_raw)
+
+    gap * (1 + leftover / N.h) * (1 + weight)
+  }
+
+
+
+  ######
+  # Moving onto the actual allocation
+  N <- sum(N.h)
+  if (allocation == "proportional") {
+    raw_allocations <- n.samp * N.h / N
+
+    weights <- N.h
+
+    #Respect lower-bound, upper bound, and n.samp
+    adjusted_allocations <- .adjust_fixed_total(
+      weights,
+      n.samp,
+      lbound,
+      N.h
+    )
+
+    rounded_allocations <- .round_fixed_total(
+      adjusted_allocations,
+      n.samp,
+      N.h,
+      lbound,
+      score = function(current) .score_proportional(adjusted_allocations, current, N.h)
+    )
+  } else if (allocation == "power") {
+    N.h.powered <- N.h^power
+    raw_allocations <- n.samp * N.h.powered / sum(N.h.powered)
+
+    weights <- N.h.powered
+
+    #Respect lower-bound, upper bound, and n.samp
+    adjusted_allocations <- .adjust_fixed_total(
+      weights,
+      n.samp,
+      lbound,
+      N.h
+    )
+
+    rounded_allocations <- .round_fixed_total(
+      adjusted_allocations,
+      n.samp,
+      N.h,
+      lbound,
+      score = function(current) .score_power(adjusted_allocations, current, N.h, power)
+    )
+
+
+  } else if (allocation == "neyman") {
+    propNum <- N.h * S.h # Numerator
+    propDen <- sum(propNum) # Denominator
+    raw_allocations <- n.samp * propNum / propDen
+
+    weights <- propNum
+
+    #Respect lower-bound, upper bound, and n.samp
+    adjusted_allocations <- .adjust_fixed_total(
+      weights,
+      n.samp,
+      lbound,
+      N.h
+    )
+
+    rounded_allocations <- .round_fixed_total(
+      adjusted_allocations,
+      n.samp,
+      N.h,
+      lbound,
+      score = function(current) .score_neyman(adjusted_allocations, current, N.h, S.h)
+    )
+  } else if (allocation == "optimal") {
+    if (!is.null(cost)) { # Cost-constrained
+      propNum <- N.h * S.h / sqrt(c.h)
+      propDen <- sum(N.h * S.h * sqrt(c.h))
+      raw_allocations <- cost * propNum / propDen
+
+      #Lower-bound adjusted version
+      lbound.h <- rep(lbound, length(N.h))
+      baseline_cost <- sum(lbound.h * c.h)
+      adjusted_target <- pmax(raw_allocations, lbound.h)
+      adjusted_target <- pmin(adjusted_target, N.h)
+
+      adjusted_cost <- sum(adjusted_target * c.h)
+
+      if (adjusted_cost > cost){
+        #Rescale only the excess above lbound to stay on budget
+        excess <- adjusted_target - lbound.h
+        excess_cost <- sum(excess * c.h)
+        rho <- if (excess_cost > 0) (cost - baseline_cost) / excess_cost else 0
+        rho <- max(min(rho, 1), 0)
+        adjusted_allocations <- lbound.h + rho * excess
+      } else {
+        adjusted_allocations <- adjusted_target
+      }
+
+      #Deterministic rounded version under budget
+      rounded_allocations <- floor(adjusted_allocations + 1e-9)
+      rounded_allocations <- pmax(rounded_allocations, lbound.h)
+      rounded_allocations <- pmin(rounded_allocations, N.h)
+      rounded_allocations <- as.integer(rounded_allocations)
+
+      remaining_budget <- cost - sum(rounded_allocations * c.h)
+      frac <- adjusted_allocations - floor(adjusted_allocations + 1e-9)
+
+      repeat {
+        candidates <- which(
+          rounded_allocations < N.h &
+            (c.h <= remaining_budget)
+        )
+        if (length(candidates) == 0) break
+
+        # prioritize largest fractional parts per unit cost
+        gap <- pmax(adjusted_allocations[candidates] - rounded_allocations[candidates], 0)
+        leftover <- pmax(N.h[candidates] - adjusted_allocations[candidates], 0)
+
+        score <- (gap / c.h[candidates]) *
+          (1 + leftover / N.h[candidates])
+
+        j <- candidates[which.max(score)]
+
+        if (score[which.max(score)] <= 0) break
+
+        rounded_allocations[j] <- rounded_allocations[j] + 1L
+        remaining_budget <- cost - sum(rounded_allocations * c.h)
+      }
+
+    } else if (!is.null(variance)) { # Precision-constrained
+      propNum <- sum(N.h * S.h * sqrt(c.h))
+      propDen <- variance * sum(N.h)**2 + sum(N.h * S.h**2)
+      raw_allocations <- N.h * S.h / sqrt(c.h) * propNum / propDen
+
+      adjusted_allocations <- .adjust_precision_constrained(N.h, S.h, c.h, variance, lbound)
+
+      rounded_allocations <- .round_precision_constrained(adjusted_allocations, N.h, S.h, c.h, variance, lbound)
+    }
+  }
+
+
+  #Note: use the largest remainder (of extra sample left) [NOT NECESSARILY HOW MUCH IS CURRENTLY ALLOCATED)
+
+  # Prep for outputting
+  if (allocation == "proportional") {
+    inputs <- list("N.h" = N.h)
+  } else if (allocation == "power") {
+    inputs <- list("N.h" = N.h, "power" = power)
+  } else if (allocation == "neyman") {
+    inputs <- list("N.h" = N.h, "S.h" = S.h)
+  } else if (allocation == "optimal" & !is.null(cost)) {
+    inputs <- list("N.h" = N.h, "S.h" = S.h, "c.h" = c.h, "cost" = cost)
+  } else if (allocation == "optimal" & !is.null(variance)) {
+    inputs <- list("N.h" = N.h, "S.h" = S.h, "c.h" = c.h, "variance" = variance)
+  }
+
+  output <- list()
+  out.length <- 0
+  if (any(outputs == "raw")){
+    out.length <- out.length + 1
+    output[out.length] <- list(raw_allocations)
+    names(output)[out.length] <- "raw"
+  }
+  if (any(outputs == "adjusted")){
+    out.length <- out.length + 1
+    output[out.length] <- list(adjusted_allocations)
+    names(output)[out.length] <- "adjusted"
+  }
+  if (any(outputs == "rounded")){
+    out.length <- out.length + 1
+    output[out.length] <- list(rounded_allocations)
+    names(output)[out.length] <- "rounded"
+  }
+
+  #output <- as.integer(rounded_allocations)
+  if (allocation == "optimal") {
+    n.print <- sum(rounded_allocations)
+    if (!is.null(c.h)){
+      actual_cost <- sum(rounded_allocations * c.h)
+      n.print <- n.print |>
+        paste0(" (sample cost: ",round(actual_cost,digits=1),")")
+    }
+  } else {
+    n.print <- n.samp
+  }
+  message(paste0("Sample allocation of ", n.print, " using ", allocation, " with the relevant inputs:"))
+  for (i in 1:length(inputs)) {
+    message(paste0("  ",
+                   names(inputs)[i],
+                   " = ",
+                   paste0(inputs[[i]], collapse = ", "),
+                   collapse = ""
+    ))
+  }
+  message()
+  message("Output:")
+  for (i in 1:length(output)) {
+    message(paste0(" ",
+                   names(output)[i],
+                   " = ",
+                   paste0(output[[i]], collapse = ", "),
+                   collapse = ""
+    ))
+  }
+
+  if (length(output) == 1){
+    output <- output[[1]]
+  }
+
+  return(output)
+}
diff --git a/man/allocate.Rd b/man/allocate.Rd
index 9823997..a7426c6 100644
--- a/man/allocate.Rd
+++ b/man/allocate.Rd
@@ -17,7 +17,8 @@ allocate(
 )
 }
 \arguments{
-\item{allocation}{type of allocation, must be one of \code{"proportional"}, \code{"power"}, \code{"neyman"}, or \code{"optimal"}.}
+\item{allocation}{type of allocation,
+must be one of \code{"proportional"}, \code{"power"}, \code{"neyman"}, or \code{"optimal"}.}
 
 \item{N.h}{vector of population stratum sizes (\eqn{N_h}, all positive values), for example \code{c(150, 600, 250)}. \cr\cr
 required for all allocation types.}
@@ -25,14 +26,16 @@ required for all allocation types.}
 \item{n.samp}{total sample size to be allocated (positive integer of length 1). \cr\cr
 required for the following allocation types: proportional, power, and Neyman, and \code{NULL} otherwise.}
 
-\item{S.h}{vector of stratum unit standard deviations (positive values same length as \code{N.h}) (\eqn{S_h}). \cr\cr
+\item{S.h}{vector of stratum unit standard deviations
+(positive values same length as \code{N.h}) (\eqn{S_h}). \cr\cr
 required for the following allocation types: Neyman, and optimal, and \code{NULL} otherwise.}
 
 \item{c.h}{vector of cost per unit in stratum h (positive values same length as \code{N.h}) (\eqn{c_h}). \cr\cr
 required for the optimal allocation only, and \code{NULL} otherwise.}
 
-\item{cost}{total variable cost (positive value) \eqn{(C - c_0)}. \cr\cr
-required for the cost-constrained optimal allocation only, and \code{NULL} otherwise.}
+\item{cost}{total variable cost (positive value) \eqn{(C – c_0)}. \cr\cr
+required for the cost-constrained optimal
+allocation only, and \code{NULL} otherwise.}
 
 \item{variance}{fixed variance target for estimated mean (positive value) (\eqn{V_0}). \cr\cr
 required for the precision-constrained optimal allocation only, and \code{NULL} otherwise.}
@@ -40,7 +43,8 @@ required for the precision-constrained optimal allocation only, and \code{NULL}
 \item{power}{power value for power allocation (\eqn{0 \le \alpha \le 1}). \cr\cr
 required for the power allocation only, and \code{NULL} otherwise.}
 
-\item{lbound}{minimum stratum-level (positive integer of length 1). Default value is 2.}
+\item{lbound}{minimum stratum-level
+(positive integer of length 1). Default value is 2.}
 }
 \value{
 Integer vector same length of N.h final allocation
@@ -68,10 +72,10 @@ where \cr
 \item Optimal allocation
 \itemize{
 \item cost-constrained [\code{N.h, S.h, c.h, cost, allocation = "optimal"}]
-\deqn{n_h = (C-c_0) \times \frac{N_h S_h / \sqrt{c_h}}{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}}
-where \cr
+\deqn{n_h = (C−c_0) \times \frac{N_h S_h / \sqrt{c_h}}{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}}
+            where \cr
 \eqn{c_h}: cost per unit in stratum \emph{h} (function input is \code{c.h}), and \cr
-\eqn{(C - c_0)}: total variable cost (function input is \code{cost})
+\eqn{(C – c_0)}: total variable cost (function input is \code{cost})
 \item precision-constrained [\code{N.h, S.h, c.h, variance, allocation = "optimal"}]
 \deqn{n_h = N_h S_h / \sqrt{c_h} \times \frac{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}{V_0 \left(\sum\limits_{h=1}^H  N_h \right)^2 + \sum\limits_{h=1}^H N_h S_h^2}}
 where \cr
@@ -81,53 +85,22 @@ where \cr
 
 The table below presents the relevant inputs for each type; when irrelevant inputs are entered, an error message will be displayed.
 
-\tabular{lllllllll}{
-\strong{allocation} \tab \strong{N.h} \tab \strong{n.samp} \tab \strong{S.h} \tab \strong{c.h} \tab \strong{cost} \tab \strong{variance} \tab \strong{lbound} \tab \strong{power} \cr
-proportional                   \tab X \tab X \tab  \tab  \tab  \tab  \tab X \tab \cr
-power                          \tab X \tab X \tab  \tab  \tab  \tab  \tab X \tab X\cr
-neyman                         \tab X \tab X \tab X \tab X \tab  \tab  \tab X \tab \cr
-optimal: cost-constrained      \tab X \tab  \tab X \tab X \tab X \tab  \tab X \tab \cr
-optimal: precision-constrained \tab X \tab  \tab X \tab X \tab  \tab X \tab X \tab
+\tabular{lllll}{
+\strong{allocation} \tab \strong{N.h} \tab \strong{n.h} \tab \strong{power} \tab \strong{lbound} \cr
+proportional \tab ✓ \tab ✓  \tab   \tab ✓\cr
+power        \tab ✓ \tab ✓  \tab ✓ \tab ✓
 }
 }
 
 \examples{
-# The first step is getting a frame summary
-#  Summarize the IPEDS dataset by OBEREG
-# - N: number of universities per region
-# - SD_ENRTOT: standard deviation of total enrollment per region
-# - Filter out rows with missing ENRTOT to ensure accurate variance estimates
-
-ipeds_summary <- ipeds |>
-  tidytable::filter(!is.na(ENRTOT)) |>
-  tidytable::group_by(OBEREG) |>
-  tidytable::summarize(
-    N = tidytable::n(),
-    SD_ENRTOT = stats::sd(ENRTOT)
-  ) |>
-  tidytable::ungroup()
-
 # Example of proportional allocation
 ipeds_summary |>
   tidytable::mutate(
     n = allocate("proportional", N.h = N, n.samp = 500)
   )
-
 # Example of power allocation
 ipeds_summary |>
   tidytable::mutate(
     n = allocate("power", N.h = N, power = 0.5, n.samp = 500)
   )
-
-# Example of Neyman allocation
-ipeds_summary |>
-  tidytable::mutate(
-    n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT)
-  )
-
-# Example of Neyman allocation with a lower bound of 5
-ipeds_summary |>
-  tidytable::mutate(
-    n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT, lbound = 5)
-  )
 }
diff --git a/tests/testthat/test-allocate.R b/tests/testthat/test-allocate.R
index 2ee4dd2..ad78b37 100644
--- a/tests/testthat/test-allocate.R
+++ b/tests/testthat/test-allocate.R
@@ -1,467 +1,604 @@
-#########
-# No/invalid allocation method
-
-test_that("allocate fails when no/invalid allocation method is provided", {
-  expect_error(allocate(),
-    regexp = 'argument "allocation" is missing, with no default'
-  )
-  expect_error(allocate(allocation = "random"),
-    regexp = "'arg' should be one of .proportional., .power., .neyman., .optimal."
-  )
-  expect_error(allocate(N.h = c(10, 20, 30)),
-    regexp = 'argument "allocation" is missing, with no default'
-  )
-  expect_error(allocate(n.h = c(10, 20, 30)),
-    regexp = "unused argument \\(n\\.h = c\\(10, 20, 30\\)\\)"
-  )
-})
-
-#########
-# Proportional allocation method
-test_that("Proportional - allocate throws error/warning when invalid inputs are provided", {
-  expect_error(allocate("proportional"),
-    regexp = 'argument "N\\.h" is missing, with no default'
-  )
-  expect_error(allocate("proportional", N.h = c(10, 20, 30)),
-    regexp = 'The n\\.samp parameter must be specified for allocation=="proportional"'
-  )
-  expect_warning(allocate("proportional", N.h = c(10, 20, 30), n.samp = 90),
-    regexp = "sum\\(N\\.h\\) is less than n.samp"
-  )
-  expect_error(allocate("proportional", N.h = c(10, 20, 30), n.samp = 5),
-    regexp = "lbound\\*length\\(N\\.h\\) must be less than or equal to n\\.samp"
-  )
-  expect_error(allocate("proportional", N.h = numeric(), n.samp = 6),
-    regexp = "The N\\.h parameter must be a vector of positive values \\(integers or non-integers\\)"
-  )
-  expect_error(allocate("proportional", N.h = c(0, 10, 20), n.samp = 6),
-    regexp = "The N\\.h parameter must be a vector of positive values \\(integers or non-integers\\)"
-  )
-  expect_error(allocate("proportional", N.h = c(-10, 10, 20), n.samp = 6),
-    regexp = "The N\\.h parameter must be a vector of positive values \\(integers or non-integers\\)"
-  )
-  expect_error(allocate("proportional", N.h = c(10, 20, 30), n.samp = 0),
-    regexp = "1: The n\\.samp parameter must be a positive integer of length 1\n2: lbound\\*length\\(N\\.h\\) must be less than or equal to n.samp"
-  )
-  expect_error(allocate("proportional", N.h = c(10, 20, 30), n.samp = -1),
-    regexp = "1: The n\\.samp parameter must be a positive integer of length 1\n2: lbound\\*length\\(N\\.h\\) must be less than or equal to n.samp"
-  )
-  expect_error(allocate("proportional", N.h = c(10, 20, 30), n.samp = c(1, 10)),
-    regexp = "1: The n\\.samp parameter must be a positive integer of length 1\n2: lbound\\*length\\(N\\.h\\) must be less than or equal to n.samp"
-  )
-})
-
-
-
-test_that("Proportional - proportions in returned sample match the relative proportions in N.h (when min(n.samp*N.h/sum(N.h)) > lbound and all((n.samp*N.h/sum(N.h)) %% 1 == 0)", {
-  # Formula: .n.samp*.N.h/sum(.N.h)
-  expect_equal(
-    allocate("proportional", N.h = c(2, 4), n.samp = 3, lbound = 1),
-    c(1, 2)
-  )
-  expect_equal(
-    allocate("proportional", N.h = c(5, 5), n.samp = 4, lbound = 1),
-    c(2, 2)
-  )
-  expect_equal(
-    allocate("proportional", N.h = c(20, 30, 40), n.samp = 9, lbound = 2),
-    c(2, 3, 4)
-  )
-  expect_equal(
-    allocate("proportional", N.h = c(25, 50, 100), n.samp = 21, lbound = 3),
-    c(3, 6, 12)
-  )
-})
-
-
-test_that("Proportional - does sum(rounded_allocations) = n.samp?", {
-  expect_equal(
-    sum(allocate("proportional", N.h = c(2, 4), n.samp = 3, lbound = 1)),
-    3
-  )
-  expect_equal(
-    sum(allocate("proportional", N.h = c(5, 5), n.samp = 4, lbound = 1)),
-    4
-  )
-  expect_equal(
-    sum(allocate("proportional", N.h = c(20, 30, 40), n.samp = 9)),
-    9
-  )
-  expect_equal(
-    sum(allocate("proportional", N.h = c(25, 50, 100), n.samp = 21, lbound = 3)),
-    21
-  )
-  expect_equal(
-    sum(allocate("proportional", N.h = c(25, 50, 100), n.samp = 9, lbound = 3)),
-    9
-  )
-})
-
-
-test_that("Proportional - does min(rounded_allocations) >= lbound?", {
-  expect_gte(
-    min(allocate("proportional", N.h = c(2, 4), n.samp = 3, lbound = 1)),
-    1
-  )
-  expect_gte(
-    min(allocate("proportional", N.h = c(5, 5), n.samp = 4, lbound = 1)),
-    1
-  )
-  expect_gte(
-    min(allocate("proportional", N.h = c(20, 30, 40), n.samp = 9)),
-    2
-  )
-  expect_gte(
-    min(allocate("proportional", N.h = c(25, 50, 100), n.samp = 21, lbound = 3)),
-    3
-  )
-  expect_gte(
-    min(allocate("proportional", N.h = c(25, 50, 100), n.samp = 9, lbound = 3)),
-    3
-  )
-})
-
-
-
-#########
-# Power allocation method
-
-test_that("Power - allocate throws error/warning when invalid inputs are provided", {
-  expect_error(allocate("power"),
-    regexp = 'argument "N\\.h" is missing, with no default'
-  )
-  expect_error(allocate("power", N.h = c(10, 20, 30)),
-    regexp = '1: The n\\.samp parameter must be specified for allocation=="power"\n2: The power parameter must be specified for allocation=="power"'
-  )
-  expect_error(allocate("power", N.h = c(10, 20, 30), n.samp = 90),
-    regexp = 'The power parameter must be specified for allocation=="power"'
-  )
-  expect_error(allocate("power", N.h = c(10, 20, 30), n.samp = 90, power = -0.5),
-    regexp = "The power parameter must be a positive value between 0 and 1, inclusive"
-  )
-  expect_error(allocate("power", N.h = c(10, 20, 30), n.samp = 90, power = 2),
-    regexp = "The power parameter must be a positive value between 0 and 1, inclusive"
-  )
-  expect_warning(allocate("power", N.h = c(10, 20, 30), n.samp = 90, power = 0.5),
-    regexp = "sum\\(N\\.h\\) is less than n.samp"
-  )
-  expect_error(allocate("power", N.h = c(10, 20, 30), n.samp = 5, power = 0.5),
-    regexp = "lbound\\*length\\(N\\.h\\) must be less than or equal to n\\.samp"
-  )
-})
-
-test_that("Power - proportions in returned sample match the unrounded sample (when min(n.samp*N.h/sum(N.h)) > lbound and expected unrounded sample are all integers)", {
-  # Formula:
-  # N.h.powered <- N.h**power
-  # N.powered <- sum(N.h.powered)
-  # allocations <- n.samp * N.h.powered / N.powered
-  expect_equal(
-    allocate("power", power = 0.5, N.h = c(4, 4, 16), n.samp = 20),
-    c(5, 5, 10)
-  )
-  expect_equal(
-    allocate("power", power = 0.25, N.h = c(16, 81), n.samp = 10),
-    c(4, 6)
-  )
-  expect_equal(
-    allocate("power", power = 1 / 3, N.h = c(8, 27), n.samp = 15),
-    c(6, 9)
-  )
-  expect_equal(
-    allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24),
-    c(6, 8, 10)
-  )
-})
-
-
-test_that("Power - does sum(rounded_allocations) = n.samp?", {
-  expect_equal(
-    sum(allocate("power", power = 0.5, N.h = c(4, 4, 16), n.samp = 20)),
-    20
-  )
-  expect_equal(
-    sum(allocate("power", power = 0.25, N.h = c(16, 81), n.samp = 10)),
-    10
-  )
-  expect_equal(
-    sum(allocate("power", power = 1 / 3, N.h = c(8, 27), n.samp = 15)),
-    15
-  )
-  expect_equal(
-    sum(allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24)),
-    24
-  )
-  expect_equal(
-    sum(allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24, lbound = 8)),
-    24
-  )
-})
-
-
-test_that("Power - does min(rounded_allocations) >= lbound?", {
-  expect_gte(
-    min(allocate("power", power = 0.5, N.h = c(4, 4, 16), n.samp = 20)),
-    2
-  )
-  expect_gte(
-    min(allocate("power", power = 0.25, N.h = c(16, 81), n.samp = 10)),
-    2
-  )
-  expect_gte(
-    min(allocate("power", power = 1 / 3, N.h = c(8, 27), n.samp = 15)),
-    2
-  )
-  expect_gte(
-    min(allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24)),
-    2
-  )
-  expect_gte(
-    min(allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24, lbound = 8)),
-    8
-  )
-})
-
-
-#########
-# Neyman allocation method
-
-test_that("Neyman - allocate throws error/warning when invalid inputs are provided", {
-  expect_error(allocate("neyman"),
-    regexp = 'argument "N\\.h" is missing, with no default'
-  )
-  expect_error(allocate("neyman", N.h = c(10, 20, 30)),
-    regexp = '1: The n\\.samp parameter must be specified for allocation=="neyman"\n2: The S\\.h parameter must be specified for allocation=="neyman"'
-  )
-  expect_error(allocate("neyman", N.h = c(10, 20, 30), n.samp = 90),
-    regexp = 'The S\\.h parameter must be specified for allocation=="neyman"'
-  )
-  expect_error(allocate("neyman", N.h = c(10, 20, 30), n.samp = 90, S.h = 0.5),
-    regexp = "The S\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N\\.h"
-  )
-  expect_error(allocate("neyman", N.h = c(10, 20), n.samp = 90, S.h = c(-0.5, 1.5)),
-    regexp = "The S\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N\\.h"
-  )
-  expect_error(allocate("neyman", N.h = c(10, 20, 30), n.samp = 90, S.h = c(0, 1, 2)),
-    regexp = "The S\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N\\.h"
-  )
-  expect_warning(allocate("neyman", N.h = c(10, 20, 30), n.samp = 90, S.h = c(0.5, 1.5, 2.5)),
-    regexp = "sum\\(N\\.h\\) is less than n.samp"
-  )
-  expect_error(allocate("neyman", N.h = c(20, 30), n.samp = 3, S.h = c(1, 2)),
-    regexp = "lbound\\*length\\(N\\.h\\) must be less than or equal to n\\.samp"
-  )
-})
-
-test_that("Neyman - proportions in returned sample match the unrounded sample (when min(n.samp*N.h/sum(N.h)) > lbound and expected unrounded sample are all integers)", {
-  # Formula:
-  # propNum <- N.h * S.h # Numerator
-  # propDen <- sum(propNum) # Denominator
-  # allocations <- n.samp * propNum / propDen
-  expect_equal(
-    allocate("neyman", N.h = c(8, 8, 32), n.samp = 20, S.h = c(1, 1, 2)),
-    c(2, 2, 16)
-  )
-  expect_equal(
-    allocate("neyman", N.h = c(10, 10, 10), n.samp = 20, S.h = c(1, 1, 2)),
-    c(5, 5, 10)
-  )
-  expect_equal(
-    allocate("neyman", N.h = c(30, 60), n.samp = 12, S.h = c(1, 1)),
-    c(4, 8)
-  )
-  expect_equal(
-    allocate("neyman", N.h = c(25, 50, 100), n.samp = 99, S.h = c(5, 15, 2.5), lbound = 1),
-    c(11, 66, 22)
-  )
-})
-
-
-test_that("Neyman - does sum(rounded_allocations) = n.samp?", {
-  expect_equal(
-    sum(allocate("neyman", N.h = c(8, 8, 32), n.samp = 20, S.h = c(1, 1, 2))),
-    20
-  )
-  expect_equal(
-    sum(allocate("neyman", N.h = c(10, 10, 10), n.samp = 20, S.h = c(1, 1, 2))),
-    20
-  )
-  expect_equal(
-    sum(allocate("neyman", N.h = c(30, 60), n.samp = 12, S.h = c(1, 1))),
-    12
-  )
-  expect_equal(
-    sum(allocate("neyman", N.h = c(25, 50, 100), n.samp = 99, S.h = c(5, 15, 2.5), lbound = 1)),
-    99
-  )
-  expect_equal(
-    sum(allocate("neyman", N.h = c(100, 50, 150), n.samp = 100, S.h = c(7.5, 10, 2.5), lbound = 25)),
-    100
-  )
-})
-
-
-test_that("Neyman - does min(rounded_allocations) >= lbound?", {
-  expect_gte(
-    min(allocate("neyman", N.h = c(8, 8, 32), n.samp = 20, S.h = c(1, 1, 2))),
-    2
-  )
-  expect_gte(
-    min(allocate("neyman", N.h = c(10, 10, 10), n.samp = 20, S.h = c(1, 1, 2))),
-    2
-  )
-  expect_gte(
-    min(allocate("neyman", N.h = c(30, 60), n.samp = 12, S.h = c(1, 1))),
-    2
-  )
-  expect_gte(
-    min(allocate("neyman", N.h = c(25, 50, 100), n.samp = 99, S.h = c(5, 15, 2.5), lbound = 1)),
-    1
-  )
-  expect_gte(
-    min(allocate("neyman", N.h = c(100, 50, 150), n.samp = 100, S.h = c(7.5, 10, 2.5), lbound = 25)),
-    25
-  )
-})
-
-
-#########
-# Optimal (cost) allocation method
-
-test_that("Optimal (cost) - allocate throws error/warning when invalid inputs are provided", {
-  expect_error(allocate("optimal"),
-    regexp = 'argument "N\\.h" is missing, with no default'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30)),
-    regexp = '1: The S\\.h parameter must be specified for allocation=="optimal"\n2: The c\\.h parameter must be specified for allocation=="optimal"\n3: Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1)),
-    regexp = '1: The S\\.h parameter must be specified for allocation=="optimal"\n2: Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5)),
-    regexp = 'Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), cost = 100, variance = 5),
-    regexp = 'Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), cost = 0),
-    regexp = "The cost parameter must be a positive value \\(integer or non-integer\\)"
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), cost = c(100, 50)),
-    regexp = "The cost parameter must be a positive value \\(integer or non-integer\\)"
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(0, 2, 1), S.h = c(0.5, 1, 1.5), cost = 100),
-    regexp = "The c\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N.h"
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2), S.h = c(0.5, 1, 1.5), cost = 100),
-    regexp = "The c\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N.h"
-  )
-})
-
-test_that("Optimal (cost) - proportions in returned sample match the unrounded sample (when min(n.samp*N.h/sum(N.h)) > lbound and expected unrounded sample are all integers)", {
-  # Formula:
-  # propNum <- N.h * S.h / sqrt(c.h)
-  # propDen <- sum(N.h * S.h * sqrt(c.h))
-  # allocations <- cost * propNum / propDen
-  expect_equal(
-    allocate("optimal", N.h = c(100, 150), S.h = c(1.5, 1), c.h = c(1, 4), cost = 24),
-    c(8, 4)
-  )
-  expect_equal(
-    allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5), c.h = c(1, 1, 4), cost = 100),
-    c(16, 12, 18)
-  )
-})
-
-
-test_that("Optimal (cost) - does the total cost of the allocated sample align with the cost parameter (when a solution exists (e.g., (sum(lbound*c.h) > cost)?", {
-  expect_equal(
-    sum(c(1, 4) * allocate("optimal", N.h = c(100, 150), S.h = c(1.5, 1), c.h = c(1, 4), cost = 24)),
-    24
-  )
-  expect_equal(
-    sum(c(1, 1, 4) * allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5), c.h = c(1, 1, 4), cost = 100)),
-    100
-  )
-})
-
-
-test_that("Optimal (cost) - does min(rounded_allocations) >= lbound?", {
-  expect_gte(
-    min(allocate("optimal", N.h = c(100, 150), S.h = c(1.5, 1), c.h = c(1, 4), cost = 24, lbound = 4)),
-    4
-  )
-  expect_gte(
-    min(allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5), c.h = c(1, 1, 4), cost = 100), lbound = 12),
-    12
-  )
-  expect_warning(expect_gte(
-    min(allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5), c.h = c(1, 4, 4), cost = 5)),
-    2
-  ))
-})
-
-
-#########
-# Optimal (variance) allocation method
-
-test_that("Optimal (variance) - allocate throws error/warning when invalid inputs are provided", {
-  expect_error(allocate("optimal"),
-    regexp = 'argument "N\\.h" is missing, with no default'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30)),
-    regexp = '1: The S\\.h parameter must be specified for allocation=="optimal"\n2: The c\\.h parameter must be specified for allocation=="optimal"\n3: Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1)),
-    regexp = '1: The S\\.h parameter must be specified for allocation=="optimal"\n2: Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5)),
-    regexp = 'Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), cost = 100, variance = 5),
-    regexp = 'Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), variance = 0),
-    regexp = "The variance parameter must be a positive value \\(integer or non-integer\\)"
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), variance = c(100, 50)),
-    regexp = "The variance parameter must be a positive value \\(integer or non-integer\\)"
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(0, 2, 1), S.h = c(0.5, 1, 1.5), variance = 100),
-    regexp = "The c\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N.h"
-  )
-  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2), S.h = c(0.5, 1, 1.5), variance = 100),
-    regexp = "The c\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N.h"
-  )
-})
-
-test_that("Optimal (variance) - proportions in returned sample match the unrounded sample (when min(n.samp*N.h/sum(N.h)) > lbound and expected unrounded sample are all integers)", {
-  # Formula:
-  # propNum <- sum(N.h * S.h * sqrt(c.h))
-  # propDen <- variance * sum(N.h)**2 + sum(N.h * S.h**2)
-  # allocations <- N.h * S.h / sqrt(c.h) * propNum / propDen
-  expect_equal(
-    allocate("optimal", N.h = c(36, 64), S.h = c(9, 9), c.h = c(4, 1), variance = 1),
-    c(11, 39)
-  )
-  expect_equal(
-    allocate("optimal", N.h = c(100, 64, 144), S.h = c(49, 36, 64), c.h = c(16, 4, 9), variance = 3.5),
-    c(51, 48, 128)
-  )
-})
-
-
-
-test_that("Optimal (variance) - does min(rounded_allocations) >= lbound?", {
-  expect_gte(
-    min(allocate("optimal", N.h = c(36, 64), S.h = c(9, 9), c.h = c(4, 1), variance = 1)),
-    2
-  )
-  expect_gte(
-    min(allocate("optimal", N.h = c(100, 64, 144), S.h = c(49, 36, 64), c.h = c(16, 4, 9), variance = 3.5)),
-    2
-  )
-  expect_gte(
-    min(allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5), c.h = c(1, 4, 4), variance = 5)),
-    2
-  )
-})
+#########
+# No/invalid allocation method
+
+test_that("allocate fails when no/invalid allocation method is provided", {
+  expect_error(allocate(),
+    regexp = 'argument "allocation" is missing, with no default'
+  )
+  expect_error(allocate(allocation = "random"),
+    regexp = "'arg' should be one of .proportional., .power., .neyman., .optimal."
+  )
+  expect_error(allocate(N.h = c(10, 20, 30)),
+    regexp = 'argument "allocation" is missing, with no default'
+  )
+  expect_error(allocate(n.h = c(10, 20, 30)),
+    regexp = "unused argument \\(n\\.h = c\\(10, 20, 30\\)\\)"
+  )
+})
+
+#########
+# Proportional allocation method
+test_that("Proportional - allocate throws error/warning when invalid inputs are provided", {
+  expect_error(allocate("proportional"),
+    regexp = 'argument "N\\.h" is missing, with no default'
+  )
+  expect_error(allocate("proportional", N.h = c(10, 20, 30)),
+    regexp = 'The n\\.samp parameter must be specified for allocation=="proportional"'
+  )
+
+  warnings <- capture_warnings(
+    expect_error(
+      allocate("proportional", N.h = c(10, 20, 30), n.samp = 90),
+      regexp = "No feasible solution: cannot distribute remaining sample\\."
+    )
+  )
+  expect_match(warnings, "sum\\(N\\.h\\) is less than n\\.samp", all = FALSE)
+
+  expect_error(allocate("proportional", N.h = c(10, 20, 30), n.samp = 5),
+    regexp = "lbound\\*length\\(N\\.h\\) must be less than or equal to n\\.samp"
+  )
+  expect_error(allocate("proportional", N.h = numeric(), n.samp = 6),
+    regexp = "The N\\.h parameter must be a vector of positive values \\(integers or non-integers\\)"
+  )
+  warnings <- capture_warnings(
+    expect_error(allocate("proportional", N.h = c(0, 10, 20), n.samp = 6),
+                 regexp = "The N\\.h parameter must be a vector of positive values \\(integers or non-integers\\)"
+    )
+  )
+  expect_match(warnings, "lbound > N\\.h for at least one stratum")
+
+  warnings <- capture_warnings(
+    expect_error(allocate("proportional", N.h = c(-10, 10, 20), n.samp = 6),
+                 regexp = "The N\\.h parameter must be a vector of positive values \\(integers or non-integers\\)"
+    )
+  )
+  expect_match(warnings, "lbound > N\\.h for at least one stratum")
+
+  expect_error(allocate("proportional", N.h = c(10, 20, 30), n.samp = 0),
+    regexp = "1: The n\\.samp parameter must be a positive integer of length 1\n2: lbound\\*length\\(N\\.h\\) must be less than or equal to n.samp"
+  )
+  expect_error(allocate("proportional", N.h = c(10, 20, 30), n.samp = -1),
+    regexp = "1: The n\\.samp parameter must be a positive integer of length 1\n2: lbound\\*length\\(N\\.h\\) must be less than or equal to n.samp"
+  )
+  expect_error(allocate("proportional", N.h = c(10, 20, 30), n.samp = c(1, 10)),
+    regexp = "1: The n\\.samp parameter must be a positive integer of length 1\n2: lbound\\*length\\(N\\.h\\) must be less than or equal to n.samp"
+  )
+})
+
+
+
+test_that("Proportional - proportions in returned sample match the relative proportions in N.h (when min(n.samp*N.h/sum(N.h)) > lbound and all((n.samp*N.h/sum(N.h)) %% 1 == 0)", {
+  # Formula: .n.samp*.N.h/sum(.N.h)
+  expect_equal(
+    allocate("proportional", N.h = c(2, 4), n.samp = 3, lbound = 1),
+    c(1, 2)
+  )
+  expect_equal(
+    allocate("proportional", N.h = c(5, 5), n.samp = 4, lbound = 1),
+    c(2, 2)
+  )
+  expect_equal(
+    allocate("proportional", N.h = c(20, 30, 40), n.samp = 9, lbound = 2),
+    c(2, 3, 4)
+  )
+  expect_equal(
+    allocate("proportional", N.h = c(25, 50, 100), n.samp = 21, lbound = 3),
+    c(3, 6, 12)
+  )
+})
+
+
+test_that("Proportional - does sum(rounded_allocations) = n.samp?", {
+  expect_equal(
+    sum(allocate("proportional", N.h = c(2, 4), n.samp = 3, lbound = 1)),
+    3
+  )
+  expect_equal(
+    sum(allocate("proportional", N.h = c(5, 5), n.samp = 4, lbound = 1)),
+    4
+  )
+  expect_equal(
+    sum(allocate("proportional", N.h = c(20, 30, 40), n.samp = 9)),
+    9
+  )
+  expect_equal(
+    sum(allocate("proportional", N.h = c(25, 50, 100), n.samp = 21, lbound = 3)),
+    21
+  )
+  expect_equal(
+    sum(allocate("proportional", N.h = c(25, 50, 100), n.samp = 9, lbound = 3)),
+    9
+  )
+})
+
+
+test_that("Proportional - does min(rounded_allocations) >= lbound?", {
+  expect_gte(
+    min(allocate("proportional", N.h = c(2, 4), n.samp = 3, lbound = 1)),
+    1
+  )
+  expect_gte(
+    min(allocate("proportional", N.h = c(5, 5), n.samp = 4, lbound = 1)),
+    1
+  )
+  expect_gte(
+    min(allocate("proportional", N.h = c(20, 30, 40), n.samp = 9)),
+    2
+  )
+  expect_gte(
+    min(allocate("proportional", N.h = c(25, 50, 100), n.samp = 21, lbound = 3)),
+    3
+  )
+  expect_gte(
+    min(allocate("proportional", N.h = c(25, 50, 100), n.samp = 9, lbound = 3)),
+    3
+  )
+})
+
+
+
+#########
+# Power allocation method
+
+test_that("Power - allocate throws error/warning when invalid inputs are provided", {
+  expect_error(allocate("power"),
+    regexp = 'argument "N\\.h" is missing, with no default'
+  )
+  expect_error(allocate("power", N.h = c(10, 20, 30)),
+    regexp = '1: The n\\.samp parameter must be specified for allocation=="power"\n2: The power parameter must be specified for allocation=="power"'
+  )
+  expect_error(allocate("power", N.h = c(10, 20, 30), n.samp = 90),
+    regexp = 'The power parameter must be specified for allocation=="power"'
+  )
+  expect_error(allocate("power", N.h = c(10, 20, 30), n.samp = 90, power = -0.5),
+    regexp = "The power parameter must be a positive value between 0 and 1, inclusive"
+  )
+  expect_error(allocate("power", N.h = c(10, 20, 30), n.samp = 90, power = 2),
+    regexp = "The power parameter must be a positive value between 0 and 1, inclusive"
+  )
+  warnings <- capture_warnings(
+    expect_error(allocate("power", N.h = c(10, 20, 30), n.samp = 90, power = 0.5),
+                   regexp = "No feasible solution: cannot distribute remaining sample."
+                   )
+  )
+
+  expect_match(warnings, "sum\\(N\\.h\\) is less than n\\.samp")
+
+
+  expect_error(allocate("power", N.h = c(10, 20, 30), n.samp = 5, power = 0.5),
+    regexp = "lbound\\*length\\(N\\.h\\) must be less than or equal to n\\.samp"
+  )
+})
+
+test_that("Power - proportions in returned sample match the unrounded sample", {
+  # Formula:
+  # N.h.powered <- N.h^power
+  # N.powered <- sum(N.h.powered)
+  # allocations <- n.samp * N.h.powered / N.powered
+  expect_equal(
+    allocate("power", power = 0.5, N.h = c(8, 8, 32), n.samp = 20),
+    c(5, 5, 10)
+  )
+  expect_equal(
+    allocate("power", power = 0.25, N.h = c(16, 81), n.samp = 10),
+    c(4, 6)
+  )
+  expect_equal(
+    allocate("power", power = 1 / 3, N.h = c(8, 27), n.samp = 15),
+    c(6, 9)
+  )
+  expect_equal(
+    allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24),
+    c(6, 8, 10)
+  )
+})
+
+
+test_that("Power - does sum(rounded_allocations) = n.samp?", {
+  expect_equal(
+    sum(allocate("power", power = 0.5, N.h = c(4, 4, 16), n.samp = 20)),
+    20
+  )
+  expect_equal(
+    sum(allocate("power", power = 0.25, N.h = c(16, 81), n.samp = 10)),
+    10
+  )
+  expect_equal(
+    sum(allocate("power", power = 1 / 3, N.h = c(8, 27), n.samp = 15)),
+    15
+  )
+  expect_equal(
+    sum(allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24)),
+    24
+  )
+  expect_equal(
+    sum(allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24, lbound = 8)),
+    24
+  )
+})
+
+
+test_that("Power - does min(rounded_allocations) >= lbound?", {
+  expect_gte(
+    min(allocate("power", power = 0.5, N.h = c(4, 4, 16), n.samp = 20)),
+    2
+  )
+  expect_gte(
+    min(allocate("power", power = 0.25, N.h = c(16, 81), n.samp = 10)),
+    2
+  )
+  expect_gte(
+    min(allocate("power", power = 1 / 3, N.h = c(8, 27), n.samp = 15)),
+    2
+  )
+  expect_gte(
+    min(allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24)),
+    2
+  )
+  expect_gte(
+    min(allocate("power", power = 1 / 2, N.h = c(9, 16, 25), n.samp = 24, lbound = 8)),
+    8
+  )
+})
+
+
+#########
+# Neyman allocation method
+
+test_that("Neyman - allocate throws error/warning when invalid inputs are provided", {
+  expect_error(allocate("neyman"),
+    regexp = 'argument "N\\.h" is missing, with no default'
+  )
+  expect_error(allocate("neyman", N.h = c(10, 20, 30)),
+    regexp = '1: The n\\.samp parameter must be specified for allocation=="neyman"\n2: The S\\.h parameter must be specified for allocation=="neyman"'
+  )
+  expect_error(allocate("neyman", N.h = c(10, 20, 30), n.samp = 90),
+    regexp = 'The S\\.h parameter must be specified for allocation=="neyman"'
+  )
+  expect_error(allocate("neyman", N.h = c(10, 20, 30), n.samp = 90, S.h = 0.5),
+    regexp = "The S\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N\\.h"
+  )
+  expect_error(allocate("neyman", N.h = c(10, 20), n.samp = 90, S.h = c(-0.5, 1.5)),
+    regexp = "The S\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N\\.h"
+  )
+  expect_error(allocate("neyman", N.h = c(10, 20, 30), n.samp = 90, S.h = c(0, 1, 2)),
+    regexp = "The S\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N\\.h"
+  )
+  warnings <- capture_warnings(
+    expect_error(allocate("neyman", N.h = c(10, 20, 30), n.samp = 90, S.h = c(0.5, 1.5, 2.5)),
+                 regexp = "No feasible solution: cannot distribute remaining sample."
+                 )
+  )
+  expect_match(warnings, "sum\\(N\\.h\\) is less than n\\.samp")
+
+  expect_error(allocate("neyman", N.h = c(20, 30), n.samp = 3, S.h = c(1, 2)),
+    regexp = "lbound\\*length\\(N\\.h\\) must be less than or equal to n\\.samp"
+  )
+})
+
+test_that("Neyman - proportions in returned sample match the unrounded sample", {
+  # Formula:
+  # propNum <- N.h * S.h # Numerator
+  # propDen <- sum(propNum) # Denominator
+  # allocations <- n.samp * propNum / propDen
+  expect_equal(
+    allocate("neyman", N.h = c(8, 8, 32), n.samp = 20, S.h = c(1, 1, 2)),
+    c(2, 2, 16)
+  )
+  expect_equal(
+    allocate("neyman", N.h = c(10, 10, 10), n.samp = 20, S.h = c(1, 1, 2)),
+    c(5, 5, 10)
+  )
+  expect_equal(
+    allocate("neyman", N.h = c(30, 60), n.samp = 12, S.h = c(1, 1)),
+    c(4, 8)
+  )
+  expect_equal(
+    allocate("neyman", N.h = c(50, 100, 200), n.samp = 99, S.h = c(5, 15, 2.5), lbound = 1),
+    c(11, 66, 22)
+  )
+})
+
+
+test_that("Neyman - does sum(rounded_allocations) = n.samp?", {
+  expect_equal(
+    sum(allocate("neyman", N.h = c(8, 8, 32), n.samp = 20, S.h = c(1, 1, 2))),
+    20
+  )
+  expect_equal(
+    sum(allocate("neyman", N.h = c(10, 10, 10), n.samp = 20, S.h = c(1, 1, 2))),
+    20
+  )
+  expect_equal(
+    sum(allocate("neyman", N.h = c(30, 60), n.samp = 12, S.h = c(1, 1))),
+    12
+  )
+  expect_equal(
+    sum(allocate("neyman", N.h = c(25, 50, 100), n.samp = 99, S.h = c(5, 15, 2.5), lbound = 1)),
+    99
+  )
+  expect_equal(
+    sum(allocate("neyman", N.h = c(100, 50, 150), n.samp = 100, S.h = c(7.5, 10, 2.5), lbound = 25)),
+    100
+  )
+})
+
+
+test_that("Neyman - does min(rounded_allocations) >= lbound?", {
+  expect_gte(
+    min(allocate("neyman", N.h = c(8, 8, 32), n.samp = 20, S.h = c(1, 1, 2))),
+    2
+  )
+  expect_gte(
+    min(allocate("neyman", N.h = c(10, 10, 10), n.samp = 20, S.h = c(1, 1, 2))),
+    2
+  )
+  expect_gte(
+    min(allocate("neyman", N.h = c(30, 60), n.samp = 12, S.h = c(1, 1))),
+    2
+  )
+  expect_gte(
+    min(allocate("neyman", N.h = c(25, 50, 100), n.samp = 99, S.h = c(5, 15, 2.5), lbound = 1)),
+    1
+  )
+  expect_gte(
+    min(allocate("neyman", N.h = c(100, 50, 150), n.samp = 100, S.h = c(7.5, 10, 2.5), lbound = 25)),
+    25
+  )
+})
+
+
+#########
+# Optimal (cost) allocation method
+
+test_that("Optimal (cost) - allocate throws error/warning when invalid inputs are provided", {
+  expect_error(allocate("optimal"),
+    regexp = 'argument "N\\.h" is missing, with no default'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30)),
+    regexp = '1: The S\\.h parameter must be specified for allocation=="optimal"\n2: The c\\.h parameter must be specified for allocation=="optimal"\n3: Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1)),
+    regexp = '1: The S\\.h parameter must be specified for allocation=="optimal"\n2: Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5)),
+    regexp = 'Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), cost = 100, variance = 5),
+    regexp = 'Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), cost = 0),
+    regexp = "The cost parameter must be a positive value \\(integer or non-integer\\)"
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), cost = c(100, 50)),
+    regexp = "The cost parameter must be a positive value \\(integer or non-integer\\)"
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(0, 2, 1), S.h = c(0.5, 1, 1.5), cost = 100),
+    regexp = "The c\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N.h"
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2), S.h = c(0.5, 1, 1.5), cost = 100),
+    regexp = "The c\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N.h"
+  )
+})
+
+test_that("Optimal (cost) - proportions in returned sample match the unrounded sample", {
+  # Formula:
+  # propNum <- N.h * S.h / sqrt(c.h)
+  # propDen <- sum(N.h * S.h * sqrt(c.h))
+  # allocations <- cost * propNum / propDen
+  expect_equal(
+    allocate("optimal", N.h = c(100, 150), S.h = c(1.5, 1), c.h = c(1, 4), cost = 24),
+    c(8, 4)
+  )
+  expect_equal(
+    allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5), c.h = c(1, 1, 4), cost = 100),
+    c(16, 12, 18)
+  )
+})
+
+
+test_that("Optimal (cost) - does the total cost of the allocated sample align with the cost parameter (when a solution exists (e.g., (sum(lbound*c.h) > cost)?", {
+  expect_equal(
+    sum(c(1, 4) * allocate("optimal", N.h = c(100, 150), S.h = c(1.5, 1), c.h = c(1, 4), cost = 24)),
+    24
+  )
+  expect_equal(
+    sum(c(1, 1, 4) * allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5), c.h = c(1, 1, 4), cost = 100)),
+    100
+  )
+})
+
+
+test_that("Optimal (cost) - does min(rounded_allocations) >= lbound?", {
+  expect_gte(
+    min(allocate("optimal", N.h = c(100, 150), S.h = c(1.5, 1), c.h = c(1, 4), cost = 24, lbound = 4)),
+    4
+  )
+  expect_gte(
+    min(allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5), c.h = c(1, 1, 4), cost = 100, lbound = 12)),
+    12
+  )
+
+  expect_error(
+    allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5),
+             c.h = c(1, 4, 4), cost = 5),
+    regexp = "sum\\(lbound\\*c\\.h\\) must be less than or equal to cost"
+  )
+
+})
+
+
+#########
+# Optimal (variance) allocation method
+
+test_that("Optimal (variance) - allocate throws error/warning when invalid inputs are provided", {
+  expect_error(allocate("optimal"),
+    regexp = 'argument "N\\.h" is missing, with no default'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30)),
+    regexp = '1: The S\\.h parameter must be specified for allocation=="optimal"\n2: The c\\.h parameter must be specified for allocation=="optimal"\n3: Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1)),
+    regexp = '1: The S\\.h parameter must be specified for allocation=="optimal"\n2: Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5)),
+    regexp = 'Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), cost = 100, variance = 5),
+    regexp = 'Exactly one of the cost and variance parameters should be supplied for allocation=="optimal"'
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), variance = 0),
+    regexp = "The variance parameter must be a positive value \\(integer or non-integer\\)"
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2, 1), S.h = c(0.5, 1, 1.5), variance = c(100, 50)),
+    regexp = "The variance parameter must be a positive value \\(integer or non-integer\\)"
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(0, 2, 1), S.h = c(0.5, 1, 1.5), variance = 100),
+    regexp = "The c\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N.h"
+  )
+  expect_error(allocate("optimal", N.h = c(10, 20, 30), c.h = c(3, 2), S.h = c(0.5, 1, 1.5), variance = 100),
+    regexp = "The c\\.h parameter must be a vector of positive values \\(integers or non-integers\\) that are the same length as N.h"
+  )
+})
+
+test_that("Optimal (variance) - rounded allocation meets variance target", {
+
+ smv <- function(n.h, N.h, S.h) {
+    N <- sum(N.h); W.h <- N.h / N
+    sum(W.h^2 * (1 - n.h / N.h) * S.h^2 / n.h)
+  }
+
+  r <- allocate("optimal", N.h = c(36, 64), S.h = c(9, 9), c.h = c(4, 1), variance = 1)
+  expect_lte(smv(r, c(36, 64), c(9, 9)), 1)
+  expect_true(all(r <= c(36, 64)))
+  expect_true(all(r >= 2))
+
+  r <- allocate("optimal", N.h = c(100, 64, 144), S.h = c(49, 36, 64), c.h = c(16, 4, 9), variance = 3.5)
+  expect_lte(smv(r, c(100, 64, 144), c(49, 36, 64)), 3.5)
+  expect_true(all(r <= c(100, 64, 144)))
+  expect_true(all(r >= 2))
+
+})
+
+
+
+test_that("Optimal (variance) - does min(rounded_allocations) >= lbound?", {
+  expect_gte(
+    min(allocate("optimal", N.h = c(36, 64), S.h = c(9, 9), c.h = c(4, 1), variance = 1)),
+    2
+  )
+  expect_gte(
+    min(allocate("optimal", N.h = c(100, 64, 144), S.h = c(49, 36, 64), c.h = c(16, 4, 9), variance = 3.5)),
+    2
+  )
+  expect_warning(
+    result <- allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5), c.h = c(1, 4, 4), variance = 5),
+    regexp = "lower-bound allocation already satisfies"
+  )
+  expect_gte(min(result), 2)
+
+})
+
+
+#####
+# outputs parameter tests
+test_that("outputs = 'rounded' (default) returns a bare integer vector", {
+  x <- allocate("proportional", N.h = c(20, 30, 40), n.samp = 9)
+  expect_false(is.list(x))
+  expect_length(x, 3)
+  expect_type(x, "integer")
+})
+
+test_that("outputs = 'raw' returns a bare numeric vector", {
+  x <- allocate("proportional", N.h = c(20, 30, 40), n.samp = 9, outputs = "raw")
+  expect_false(is.list(x))
+  expect_length(x, 3)
+  expect_type(x, "double")
+})
+
+test_that("outputs = c('raw','adjusted','rounded') returns a named list in order", {
+  x <- allocate("proportional", N.h = c(20, 30, 40), n.samp = 9,
+                outputs = c("raw", "adjusted", "rounded"))
+  expect_type(x, "list")
+  expect_named(x, c("raw", "adjusted", "rounded"))
+  expect_length(x$raw, 3)
+  expect_length(x$adjusted, 3)
+  expect_length(x$rounded, 3)
+})
+
+test_that("outputs = c('adjusted','rounded') returns a two-element named list", {
+  x <- allocate("neyman", N.h = c(30, 60), n.samp = 12, S.h = c(1, 1),
+                outputs = c("adjusted", "rounded"))
+  expect_type(x, "list")
+  expect_named(x, c("adjusted", "rounded"))
+})
+
+test_that("outputs with invalid value errors", {
+  expect_error(
+    allocate("proportional", N.h = c(20, 30, 40), n.samp = 9, outputs = "foo"),
+    regexp = "'arg' should be one of"
+  )
+})
+
+
+######
+# Upper-bound invariant across allocation modes
+test_that("Rounded allocations never exceed N.h", {
+  # Capacity-limited proportional (n.samp equals sum(N.h))
+  r <- allocate("proportional", N.h = c(4, 4, 16), n.samp = 20)
+  expect_true(all(r <= c(4, 4, 16)))
+  expect_equal(sum(r), 20)
+
+  # Capacity-limited Neyman with a stratum forced against its cap
+  r <- allocate("neyman", N.h = c(25, 50, 100), n.samp = 99,
+                S.h = c(5, 15, 2.5), lbound = 1)
+  expect_true(all(r <= c(25, 50, 100)))
+  expect_equal(sum(r), 99)
+
+  # Capacity-limited Power
+  r <- allocate("power", N.h = c(4, 4, 16), n.samp = 20, power = 0.5)
+  expect_true(all(r <= c(4, 4, 16)))
+  expect_equal(sum(r), 20)
+})
+
+
+
+######
+# Cost-constrained: budget invariant
+test_that("Cost-constrained allocation never exceeds budget", {
+  r <- allocate("optimal", N.h = c(100, 150), S.h = c(1.5, 1),
+                c.h = c(1, 4), cost = 24)
+  expect_lte(sum(r * c(1, 4)), 24)
+
+  r <- allocate("optimal", N.h = c(100, 75, 150), S.h = c(1, 1, 1.5),
+                c.h = c(1, 1, 4), cost = 100)
+  expect_lte(sum(r * c(1, 1, 4)), 100)
+})
+
+######
+# lbound == N.h edge case
+test_that("lbound == N.h in one stratum is handled", {
+  r <- allocate("proportional", N.h = c(5, 20, 30), n.samp = 15, lbound = 5)
+  expect_equal(r[1], 5L)
+  expect_equal(sum(r), 15)
+  expect_true(all(r <= c(5, 20, 30)))
+  expect_true(all(r >= 5))
+})
+
+
+#####
+# Irrelevant-parameter warnings
+test_that("Irrelevant parameters trigger warnings", {
+  expect_warning(
+    allocate("proportional", N.h = c(10, 20), n.samp = 5, S.h = c(1, 2)),
+    regexp = "S\\.h parameter should be NULL"
+  )
+  expect_warning(
+    allocate("neyman", N.h = c(10, 20), n.samp = 5, S.h = c(1, 2), power = 0.5),
+    regexp = "power parameter should be NULL"
+  )
+  expect_warning(
+    allocate("proportional", N.h = c(10, 20), n.samp = 5, cost = 100),
+    regexp = "cost parameter should be NULL"
+  )
+})

From a1c68e603bae5ca5f6890efa80cf18b07d53c54c Mon Sep 17 00:00:00 2001
From: JD Bunker <jbunker7@outlook.com>
Date: Tue, 21 Jul 2026 14:45:12 -0500
Subject: [PATCH 2/2] Minor allocate() Tweaks - 21Jul2026
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Making minor tweaks to the allocate() function related to non-ASCII symbols (`–` to `-` and `✓` to `X`). I also regenerated the allocate.Rd file, which incorporates the new `outputs` argument and seemingly expands the table and examples section.
---
 R/allocate.R    | 1562 +++++++++++++++++++++++------------------------
 man/allocate.Rd |   71 ++-
 2 files changed, 836 insertions(+), 797 deletions(-)

diff --git a/R/allocate.R b/R/allocate.R
index b0bbc53..649c2e1 100644
--- a/R/allocate.R
+++ b/R/allocate.R
@@ -1,781 +1,781 @@
-#' Sample allocation
-#'
-#' Compute the proportional, power, Neyman, and optimal sample allocations.
-#'
-#' @section Method:
-#' The \emph{allocate} function allocates a sample size \emph{n} on \emph{H} strata using one of the following allocation methods:
-#' \enumerate{
-#'   \item Proportional allocation \[\code{n.samp, N.h, allocation = "proportional"}\]
-#'         \deqn{n_h = n \times \frac{N_h}{\sum\limits_{h=1}^H N_h}}
-#'         where \cr
-#'         \eqn{n}: total sample size to be allocated (function input is \code{n.samp}), and \cr
-#'         \eqn{N_h}: population size of stratum \emph{h} (function input is \code{N.h}).
-#'   \item Power allocation \[\code{n.samp, N.h, power, allocation = "power"}\]
-#'         \deqn{n_h = n \times \frac{N_h^\alpha}{\sum\limits_{h=1}^H N_h^\alpha}}
-#'          where \cr
-#'          \eqn{\alpha}: a power value to control over-under-sampling with \eqn{0 \le \alpha \le 1} (function input is \code{power}).
-#'   \item Neyman allocation \[\code{n.samp, N.h, S.h, allocation = "neyman"}\]
-#'         \deqn{n_h = n \times \frac{N_h S_h}{\sum\limits_{h=1}^H N_h S_h}}
-#'         where \cr
-#'         \eqn{S_h}:  standard deviation of stratum \emph{h} (function input is \code{S.h}).
-#'    \item Optimal allocation
-#'           \itemize{
-#'             \item cost-constrained \[\code{N.h, S.h, c.h, cost, allocation = "optimal"}\]
-#'                   \deqn{n_h = (C−c_0) \times \frac{N_h S_h / \sqrt{c_h}}{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}}
-#'             where \cr
-#'             \eqn{c_h}: cost per unit in stratum \emph{h} (function input is \code{c.h}), and \cr
-#'             \eqn{(C – c_0)}: total variable cost (function input is \code{cost})
-#'             \item precision-constrained \[\code{N.h, S.h, c.h, variance, allocation = "optimal"}\]
-#'                   \deqn{n_h = N_h S_h / \sqrt{c_h} \times \frac{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}{V_0 \left(\sum\limits_{h=1}^H  N_h \right)^2 + \sum\limits_{h=1}^H N_h S_h^2}}
-#'             where \cr
-#'             \eqn{V_0}: fixed variance target for estimated mean (function input is \code{variance})
-#'           }
-#' }
-#'
-#' The table below presents the relevant inputs for each type; when irrelevant inputs are entered, an error message will be displayed.
-#'
-#' \tabular{lllllllll}{
-#'   \strong{allocation} \tab \strong{N.h} \tab \strong{n.samp} \tab \strong{S.h} \tab \strong{c.h} \tab \strong{cost} \tab \strong{variance} \tab \strong{lbound} \tab \strong{power} \cr
-#'   proportional                   \tab ✓ \tab ✓ \tab  \tab  \tab  \tab  \tab ✓ \tab \cr
-#'   power                          \tab ✓ \tab ✓ \tab  \tab  \tab  \tab  \tab ✓ \tab ✓\cr
-#'   neyman                         \tab ✓ \tab ✓ \tab ✓ \tab ✓ \tab  \tab  \tab ✓ \tab \cr
-#'   optimal: cost-constrained      \tab ✓ \tab  \tab ✓ \tab ✓ \tab ✓ \tab  \tab ✓ \tab \cr
-#'   optimal: precision-constrained \tab ✓ \tab  \tab ✓ \tab ✓ \tab  \tab ✓ \tab ✓ \tab
-#' }
-#'
-#' @param allocation type of allocation, must be one of \code{"proportional"}, \code{"power"}, \code{"neyman"}, or \code{"optimal"}.
-#' @param n.samp total sample size to be allocated (positive integer of length 1). \cr\cr
-#'        required for the following allocation types: proportional, power, and Neyman, and \code{NULL} otherwise.
-#' @param N.h vector of population stratum sizes (\eqn{N_h}, all positive values), for example \code{c(150, 600, 250)}. \cr\cr
-#' required for all allocation types.
-#' @param S.h vector of stratum unit standard deviations (positive values same length as \code{N.h}) (\eqn{S_h}). \cr\cr
-#' required for the following allocation types: Neyman, and optimal, and \code{NULL} otherwise.
-#' @param c.h vector of cost per unit in stratum h (positive values same length as \code{N.h}) (\eqn{c_h}). \cr\cr
-#' required for the optimal allocation only, and \code{NULL} otherwise.
-#' @param cost total variable cost (positive value) \eqn{(C – c_0)}. \cr\cr
-#' required for the cost-constrained optimal allocation only, and \code{NULL} otherwise.
-#' @param variance fixed variance target for estimated mean (positive value) (\eqn{V_0}). \cr\cr
-#' required for the precision-constrained optimal allocation only, and \code{NULL} otherwise.
-#' @param power power value for power allocation (\eqn{0 \le \alpha \le 1}). \cr\cr
-#' required for the power allocation only, and \code{NULL} otherwise.
-#' @param lbound minimum stratum-level (positive integer of length 1). Default value is 2.
-#' @param outputs
-#' character vector representing whether to output:\cr
-#'  \enumerate{
-#'    \item the raw allocations before accounting for N.h, lbound, and n.samp if needed \[\code{"raw"}\],
-#'    \item the continuous version after accounting for the above \[\code{"adjusted"}\], and/or
-#'    \item the rounded version of the above \[\code{"rounded"}\]
-#'  } \cr
-#'  Default is to only return the final rounded version \[\code{"rounded"}\]. \cr
-#'  If one version is requested, the result will be a numeric vector. Otherwise, the result will be a named list matching the requested outputs.
-#' @return If one output type is requested, a numeric vector of allocations.
-#' If multiple output types are requested, a named list containing the requested
-#' allocation vectors.
-#' @export
-#' @examples
-#' # The first step is getting a frame summary
-#' #  Summarize the IPEDS dataset by OBEREG
-#' # - N: number of universities per region
-#' # - SD_ENRTOT: standard deviation of total enrollment per region
-#' # - Filter out rows with missing ENRTOT to ensure accurate variance estimates
-#'
-#' ipeds_summary <- ipeds |>
-#'   tidytable::filter(!is.na(ENRTOT)) |>
-#'   tidytable::group_by(OBEREG) |>
-#'   tidytable::summarize(
-#'     N = tidytable::n(),
-#'     SD_ENRTOT = stats::sd(ENRTOT)
-#'   ) |>
-#'   tidytable::ungroup()
-#'
-#' # Example of proportional allocation
-#' ipeds_summary |>
-#'   tidytable::mutate(
-#'     n = allocate("proportional", N.h = N, n.samp = 500)
-#'   )
-#'
-#' # Example of power allocation
-#' ipeds_summary |>
-#'   tidytable::mutate(
-#'     n = allocate("power", N.h = N, power = 0.5, n.samp = 500)
-#'   )
-#'
-#' # Example of Neyman allocation
-#' ipeds_summary |>
-#'   tidytable::mutate(
-#'     n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT)
-#'   )
-#'
-#' # Example of Neyman allocation with a lower bound of 5
-#' ipeds_summary |>
-#'   tidytable::mutate(
-#'     n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT, lbound = 5)
-#'   )
-
-allocate <- function(allocation, N.h, n.samp = NULL, S.h = NULL, c.h = NULL, cost = NULL, variance = NULL, power = NULL, lbound = 2, outputs = "rounded") {
-  allocation <- match.arg(allocation, c("proportional", "power", "neyman", "optimal"))
-  outputs <- match.arg(outputs, c("raw", "adjusted", "rounded"), several.ok = TRUE)
-
-  ######
-  # Check inputs
-  .problems <- NULL # Initialize list of problems found with inputs
-  .addProblem <- function(parameter, condition, problems = .problems, allocation = NULL) { # Function to simplify addition to problems found with inputs to our running list (.problems)
-    if (is.null(parameter)) { # No parameter given
-      problem <- condition
-    } else if (parameter %in% c("allocation", "n.samp", "N.h", "S.h", "c.h", "cost", "variance", "power", "lbound")) {
-      problem <- paste0("The ", parameter, " parameter ", condition)
-    }
-    problemsNew <- c(problems, problem)
-    return(problemsNew)
-  }
-
-  .condition <- paste0('must be specified for allocation=="', allocation, '"')
-  # n.samp parameter
-  if (allocation %in% c("proportional", "power", "neyman") & is.null(n.samp)) {
-    .problems <- .addProblem(parameter = "n.samp", condition = .condition)
-  } else if (!allocation %in% c("proportional", "power", "neyman") & !is.null(n.samp)) {
-    warning("The n.samp parameter should be NULL", call. = FALSE)
-  }
-  # N.h parameter
-  if (allocation %in% c("proportional", "power", "neyman", "optimal") & is.null(N.h)) {
-    .problems <- .addProblem(parameter = "N.h", condition = .condition)
-  }
-  # S.h parameter
-  if (allocation %in% c("neyman", "optimal") & is.null(S.h)) {
-    .problems <- .addProblem(parameter = "S.h", condition = .condition)
-  } else if (!allocation %in% c("neyman", "optimal") & !is.null(S.h)) {
-    warning("The S.h parameter should be NULL", call. = FALSE)
-  }
-  # c.h parameter
-  if (allocation %in% c("optimal") & is.null(c.h)) {
-    .problems <- .addProblem(parameter = "c.h", condition = .condition)
-  } else if (!allocation %in% c("optimal") & !is.null(c.h)) {
-    warning("The c.h parameter should be NULL", call. = FALSE)
-  }
-  # power parameter
-  if (allocation %in% c("power") & is.null(power)) {
-    .problems <- .addProblem(parameter = "power", condition = .condition)
-  } else if (!allocation %in% c("power") & !is.null(power)) {
-    warning("The power parameter should be NULL", call. = FALSE)
-  }
-  # optimal allocation: only 1 of cost or variance should be provided
-  if (allocation == "optimal" & sum(is.null(cost), is.null(variance)) != 1) {
-    .problems <- .addProblem(parameter = NULL, condition = paste0('Exactly one of the cost and variance parameters should be supplied for allocation=="', allocation, '"'))
-  } else if (!allocation %in% c("optimal") & !(is.null(cost) & is.null(variance))) {
-    if (!is.null(cost)) {
-      warning("The cost parameter should be NULL", call. = FALSE)
-    }
-    if (!is.null(variance)) {
-      warning("The variance parameter should be NULL", call. = FALSE)
-    }
-  }
-
-  ###
-  # Check for miscellaneous unexpected parameter values
-
-  # n.samp
-  .condition <- "must be a positive integer of length 1"
-  if (!is.null(n.samp) & allocation %in% c("proportional", "power", "neyman")) {
-    if (!all(length(n.samp) == 1 & (typeof(n.samp) %in% c("integer") | (typeof(n.samp) == "double" & round(n.samp) == n.samp)) & n.samp > 0)) {
-      .problems <- .addProblem(parameter = "n.samp", condition = .condition)
-    }
-  }
-  # lbound
-  if (!(length(lbound) == 1 & (typeof(lbound) %in% c("integer") | (typeof(lbound) == "double" & round(lbound) == lbound)) & lbound > 0)) {
-    .problems <- .addProblem(parameter = "lbound", condition = .condition)
-  }
-
-  .condition <- "must be a vector of positive values (integers or non-integers)"
-  # N.h
-  if (!is.null(N.h)) {
-    if (!(length(N.h) >= 1 & typeof(N.h) %in% c("integer", "double") & all(N.h > 0))) {
-      .problems <- .addProblem(parameter = "N.h", condition = .condition)
-    }
-  }
-  .condition <- paste0(.condition, " that are the same length as N.h")
-  # S.h
-  if (!is.null(S.h) & allocation %in% c("neyman", "optimal")) {
-    if (!(length(S.h) >= 1 & typeof(S.h) %in% c("integer", "double") & all(S.h > 0) & length(S.h) == length(N.h))) {
-      .problems <- .addProblem(parameter = "S.h", condition = .condition)
-    }
-  }
-  # c.h
-  if (!is.null(c.h) & allocation %in% c("optimal")) {
-    if (!(length(c.h) >= 1 & typeof(c.h) %in% c("integer", "double") & all(c.h > 0) & length(c.h) == length(N.h))) {
-      .problems <- .addProblem(parameter = "c.h", condition = .condition)
-    }
-  }
-  # cost
-  .condition <- "must be a positive value (integer or non-integer)"
-  if (!is.null(cost) & allocation %in% "optimal") {
-    if (!all(length(cost) == 1 & typeof(cost) %in% c("integer", "double") & cost > 0)) {
-      .problems <- .addProblem(parameter = "cost", condition = .condition)
-    }
-  }
-  # variance
-  if (allocation == "optimal" & !is.null(variance)) {
-    if (!all(length(variance) == 1 & typeof(variance) %in% c("integer", "double") & variance > 0 & length(variance) == 1)) {
-      .problems <- .addProblem(parameter = "variance", condition = .condition)
-    }
-  }
-
-  # power parameter
-  .condition <- "must be a positive value between 0 and 1, inclusive"
-  if (allocation == "power" & !is.null(power)) {
-    if (!(length(power) == 1 & typeof(power) %in% c("integer", "double") & 0 <= power & power <= 1)) {
-      .problems <- .addProblem(parameter = "power", condition = .condition)
-    }
-  }
-
-
-  if (allocation %in% c("proportional", "power", "neyman")) {
-    if (!(is.null(lbound) | is.null(N.h) | is.null(n.samp))) {
-      if (!all(lbound * length(N.h) <= n.samp)) {
-        .problems <- c(.problems, "lbound*length(N.h) must be less than or equal to n.samp")
-      }
-    }
-  }
-
-  if (!(is.null(N.h) | is.null(n.samp))) {
-    if (length(.problems) == 0) {
-      if (sum(N.h) < n.samp) {
-        warning("sum(N.h) is less than n.samp")
-      }
-    }
-  }
-
-  if (allocation %in% c("optimal")) {
-    if (!is.null(cost) & length(.problems) == 0) {
-      if (sum(lbound * c.h) > cost) {
-        .problems <- c(.problems, "sum(lbound*c.h) must be less than or equal to cost")
-      }
-    }
-  }
-
-  if (allocation %in% c("proportional", "power", "neyman","optimal")){
-    if (!is.null(N.h)){
-      if (any(lbound > N.h)){
-        warning("lbound > N.h for at least one stratum")
-      }
-    }
-  }
-
-  # outputs parameter
-  .condition <- 'must be one or more of: "raw", "adjusted", "rounded"'
-  if (!all(outputs %in% c("raw", "adjusted", "rounded"))) {
-    .problems <- .addProblem(parameter = "outputs", condition = .condition)
-  }
-
-  ###
-  # Aggregate problems and stop the program if necessary
-  if (length(.problems) > 0) {
-    if (length(.problems) > 1) {
-      .problems <- paste0(
-        1:length(.problems),
-        ": ",
-        .problems
-      )
-    }
-    stop(
-      "\n",
-      paste0(.problems, collapse = "\n")
-    )
-  }
-
-  ######
-  #Fixed-total continuous adjustment:
-  # * Respects lower bound (lbound)
-  # * Respects upper bound (N.h)
-  # * Sums to sample size (n.samp)
-  # * Unrounded sample sizes
-  .adjust_fixed_total <- function(weights, n.samp, lbound, N.h) {
-    h <- length(weights)
-
-    lbound.h <- rep(lbound, h) #Extending lbound to a vector
-
-    # Raw proportional target
-    raw <- n.samp * weights / sum(weights)
-
-    #If raw already respects both bounds, use it directly.
-    if (all(raw >= lbound.h) && all(raw <= N.h)) {
-      return(as.numeric(raw))
-    }
-
-    #Otherwise fall back to iterative capacity-and-floor enforcement
-    alloc <- lbound.h
-    remaining <- n.samp - sum(alloc)
-    capacity <- N.h - alloc #Amount left in N.h after allocating lbound
-
-    if (remaining == 0) {
-      return(as.numeric(alloc))
-    }
-
-    #Initialize strata variables
-    active <- rep(TRUE, h) #Does stratum still need to be processed
-    extra <- rep(0, h) #How much to allocate to stratum above raw_allocations
-
-
-    repeat {
-      w <- weights
-      w[!active] <- 0
-
-      if (sum(w) <= 0) {
-        stop("No feasible solution: cannot distribute remaining sample.")
-      }
-
-      proposed <- remaining * w / sum(w) #Proposed additional allocation
-      hit <- active & (proposed > capacity) #Assigned all of stratum
-
-      #If didn't hit any statum's upper bound, break here
-      if (!any(hit)) {
-        extra[active] <- remaining * weights[active] / sum(weights[active])
-        break
-      }
-
-      #If hit a stratum's N, allocate its remaining capacity and set inactive
-      alloc[hit] <- alloc[hit] + capacity[hit]
-      remaining <- n.samp - sum(alloc)
-      active[hit] <- FALSE
-      capacity[hit] <- 0
-    }
-
-
-    adjusted_allocations <- alloc + extra
-    adjusted_allocations
-  }
-
-  ######
-  # Calculate stratified mean variance
-
-  .stratified_mean_variance <- function(n.h, N.h, S.h, N = sum(N.h)) {
-    W.h <- N.h / N
-    sum(W.h^2 * (1 - n.h / N.h) * S.h^2 / n.h)
-  }
-
-  ######
-  # Precision-constrained continuous adjustment:
-  # * Respects lower bound (lbound)
-  # * Respects upper bound (N.h)
-  # * Hits variance target
-  # * Unrounded sample sizes
-
-  .adjust_precision_constrained <- function(N.h, S.h, c.h, variance, lbound){
-    h <- length(N.h)
-    N <- sum(N.h)
-
-    a.h <- N.h * S.h / sqrt(c.h)
-
-    lbound.h <- rep(lbound, h)
-
-    # Initialize strata variables
-    adjusted_allocations <- rep(NA_real_, h)
-    fixed <- rep(FALSE, h) #Treat stratum fixed (TRUE) vs. free (FALSE)
-
-
-    V.lower <- .stratified_mean_variance(lbound.h, N.h, S.h)
-    V.upper <- .stratified_mean_variance(N.h, N.h, S.h)
-
-    if (variance >= V.lower) {
-      warning("The lower-bound allocation already satisfies the variance target; returning lower-bound allocation.")
-      return(lbound.h)
-    }
-
-    if (variance < V.upper) {
-      stop("No feasible allocation: even a census of all strata does not achieve the variance target.")
-    }
-
-
-    repeat {
-      # Compute variance contribution from fixed strata
-      if (any(fixed)) {
-        V.fixed <- .stratified_mean_variance(
-          n.h = adjusted_allocations[fixed],
-          N.h = N.h[fixed],
-          S.h = S.h[fixed],
-          N = N
-        )
-      } else {
-        V.fixed <- 0
-      }
-
-      # If all strata are fixed, break
-      if (all(fixed)) {
-        V.total <- .stratified_mean_variance(adjusted_allocations, N.h, S.h)
-        if (abs(V.total - variance) > 1e-10) {
-          warning("No feasible allocation exactly meets the variance target under the imposed bounds.")
-        }
-        break
-      }
-
-      # Recompute scaling constant over free strata
-      free <- !fixed
-      A.free <- sum(N.h[free] * S.h[free] * sqrt(c.h[free]))
-      B.free <- sum(N.h[free] * S.h[free]^2)
-
-      denom <- N^2 * (variance - V.fixed) + B.free
-
-      if (denom <= 0) {
-        stop("No feasible allocation solution for the requested variance after applying bounds.")
-      }
-
-      k.free <- A.free / denom
-
-      # Candidate allocations for free strata
-      candidate <- k.free * a.h[free]
-
-      # Check for violations
-      below <- candidate < lbound.h[free]
-      above <- candidate > N.h[free]
-
-      if (!any(below | above)){
-        adjusted_allocations[free] <- candidate
-        break
-      }
-
-      # Fix violating strata
-      free_idx <- which(free)
-
-      if (any(below)) {
-        idx <- free_idx[below]
-        adjusted_allocations[idx] <- lbound.h[idx]
-        fixed[idx] <- TRUE
-      }
-
-      if (any(above)) {
-        idx <- free_idx[above]
-        adjusted_allocations[idx] <- N.h[idx]
-        fixed[idx] <- TRUE
-      }
-
-    }
-    return(adjusted_allocations)
-  }
-
-  ######
-  # Precision-constrained rounding
-
-  .round_precision_constrained <- function(adjusted_allocations, N.h, S.h, c.h, variance, lbound) {
-    rounded_allocations <- floor(adjusted_allocations + 1e-9)
-    rounded_allocations <- pmax(rounded_allocations, lbound)
-    rounded_allocations <- pmin(rounded_allocations, N.h)
-    rounded_allocations <- as.integer(rounded_allocations)
-
-    frac <- adjusted_allocations - floor(adjusted_allocations + 1e-9)
-
-    while (.stratified_mean_variance(rounded_allocations, N.h, S.h) > variance) {
-      candidates <- which(rounded_allocations < N.h)
-
-      if (length(candidates) == 0) {
-        warning("No feasible integer allocation found to satisfy variance target.")
-        break
-      }
-
-      current_var <- .stratified_mean_variance(rounded_allocations, N.h, S.h)
-      scores <- rep(-Inf, length(N.h))
-
-      for (h in candidates) {
-        trial <- rounded_allocations
-        trial[h] <- trial[h] + 1L
-        new_var <- .stratified_mean_variance(trial, N.h, S.h)
-        var_gain <- current_var - new_var
-
-        #leftover relative to adjusted continuous allocation
-        leftover <- pmax(N.h[h] - adjusted_allocations[h], 0)
-
-        #normalize leftover
-        leftover_factor <- 1 + leftover / N.h[h]
-
-        #fractional preference
-        frac_factor <- 1 + frac[h]
-
-        #combined score
-        scores[h] <- (var_gain / c.h[h]) * leftover_factor * frac_factor
-      }
-
-      best_h <- which.max(scores)
-      if (!is.finite(scores[best_h]) || scores[best_h] <= 0) {
-        warning("Could not improve enough to satisfy variance constraint.")
-        break
-      }
-
-      rounded_allocations[best_h] <- rounded_allocations[best_h] + 1L
-    }
-
-    return(rounded_allocations)
-  }
-
-
-  ######
-  # Rounding fixed-totals
-
-  .round_fixed_total <- function(adjusted_allocations, n.samp, N.h, lbound, score) {
-    rounded_allocations <- floor(adjusted_allocations + 1e-9)
-    rounded_allocations <- pmax(rounded_allocations, lbound)
-    rounded_allocations <- pmin(rounded_allocations, N.h)
-    rounded_allocations <- as.integer(rounded_allocations)
-
-    deficit <- n.samp - sum(rounded_allocations)
-
-    if (deficit < 0) {
-      stop("Initial rounded allocation exceeds the target total.")
-    }
-    while (deficit > 0) {
-      candidates <- which(rounded_allocations < N.h)
-      if (length(candidates) == 0) {
-        stop("No feasible way to allocate remaining units.")
-      }
-
-
-      scores <- score(rounded_allocations)
-      scores[-candidates] <- -Inf
-
-
-      j <- candidates[which.max(scores[candidates])]
-
-
-      if (!is.finite(scores[j])) {
-        stop("No feasible scoring candidate for additional allocation.")
-      }
-
-      rounded_allocations[j] <- rounded_allocations[j] + 1L
-      deficit <- deficit - 1L
-    }
-
-    return(rounded_allocations)
-
-  }
-
-  .score_proportional <- function(adjusted_allocations, rounded_allocations, N.h) {
-    gap <- pmax(adjusted_allocations - rounded_allocations, 0)
-    leftover <- pmax(N.h - adjusted_allocations, 0)
-    weight <- N.h / sum(N.h)
-
-    gap * (1 + leftover / N.h) * (1 + weight)
-  }
-
-  .score_power <- function(adjusted_allocations, rounded_allocations, N.h, power) {
-    gap <- pmax(adjusted_allocations - rounded_allocations, 0)
-    leftover <- pmax(N.h - adjusted_allocations, 0)
-    weight_raw <- N.h^power
-    weight <- weight_raw / sum(weight_raw)
-
-    gap * (1 + leftover / N.h) * (1 + weight)
-  }
-
-  .score_neyman <- function(adjusted_allocations, rounded_allocations, N.h, S.h) {
-    gap <- pmax(adjusted_allocations - rounded_allocations, 0)
-    leftover <- pmax(N.h - adjusted_allocations, 0)
-    weight_raw <- N.h * S.h
-    weight <- weight_raw / sum(weight_raw)
-
-    gap * (1 + leftover / N.h) * (1 + weight)
-  }
-
-
-
-  ######
-  # Moving onto the actual allocation
-  N <- sum(N.h)
-  if (allocation == "proportional") {
-    raw_allocations <- n.samp * N.h / N
-
-    weights <- N.h
-
-    #Respect lower-bound, upper bound, and n.samp
-    adjusted_allocations <- .adjust_fixed_total(
-      weights,
-      n.samp,
-      lbound,
-      N.h
-    )
-
-    rounded_allocations <- .round_fixed_total(
-      adjusted_allocations,
-      n.samp,
-      N.h,
-      lbound,
-      score = function(current) .score_proportional(adjusted_allocations, current, N.h)
-    )
-  } else if (allocation == "power") {
-    N.h.powered <- N.h^power
-    raw_allocations <- n.samp * N.h.powered / sum(N.h.powered)
-
-    weights <- N.h.powered
-
-    #Respect lower-bound, upper bound, and n.samp
-    adjusted_allocations <- .adjust_fixed_total(
-      weights,
-      n.samp,
-      lbound,
-      N.h
-    )
-
-    rounded_allocations <- .round_fixed_total(
-      adjusted_allocations,
-      n.samp,
-      N.h,
-      lbound,
-      score = function(current) .score_power(adjusted_allocations, current, N.h, power)
-    )
-
-
-  } else if (allocation == "neyman") {
-    propNum <- N.h * S.h # Numerator
-    propDen <- sum(propNum) # Denominator
-    raw_allocations <- n.samp * propNum / propDen
-
-    weights <- propNum
-
-    #Respect lower-bound, upper bound, and n.samp
-    adjusted_allocations <- .adjust_fixed_total(
-      weights,
-      n.samp,
-      lbound,
-      N.h
-    )
-
-    rounded_allocations <- .round_fixed_total(
-      adjusted_allocations,
-      n.samp,
-      N.h,
-      lbound,
-      score = function(current) .score_neyman(adjusted_allocations, current, N.h, S.h)
-    )
-  } else if (allocation == "optimal") {
-    if (!is.null(cost)) { # Cost-constrained
-      propNum <- N.h * S.h / sqrt(c.h)
-      propDen <- sum(N.h * S.h * sqrt(c.h))
-      raw_allocations <- cost * propNum / propDen
-
-      #Lower-bound adjusted version
-      lbound.h <- rep(lbound, length(N.h))
-      baseline_cost <- sum(lbound.h * c.h)
-      adjusted_target <- pmax(raw_allocations, lbound.h)
-      adjusted_target <- pmin(adjusted_target, N.h)
-
-      adjusted_cost <- sum(adjusted_target * c.h)
-
-      if (adjusted_cost > cost){
-        #Rescale only the excess above lbound to stay on budget
-        excess <- adjusted_target - lbound.h
-        excess_cost <- sum(excess * c.h)
-        rho <- if (excess_cost > 0) (cost - baseline_cost) / excess_cost else 0
-        rho <- max(min(rho, 1), 0)
-        adjusted_allocations <- lbound.h + rho * excess
-      } else {
-        adjusted_allocations <- adjusted_target
-      }
-
-      #Deterministic rounded version under budget
-      rounded_allocations <- floor(adjusted_allocations + 1e-9)
-      rounded_allocations <- pmax(rounded_allocations, lbound.h)
-      rounded_allocations <- pmin(rounded_allocations, N.h)
-      rounded_allocations <- as.integer(rounded_allocations)
-
-      remaining_budget <- cost - sum(rounded_allocations * c.h)
-      frac <- adjusted_allocations - floor(adjusted_allocations + 1e-9)
-
-      repeat {
-        candidates <- which(
-          rounded_allocations < N.h &
-            (c.h <= remaining_budget)
-        )
-        if (length(candidates) == 0) break
-
-        # prioritize largest fractional parts per unit cost
-        gap <- pmax(adjusted_allocations[candidates] - rounded_allocations[candidates], 0)
-        leftover <- pmax(N.h[candidates] - adjusted_allocations[candidates], 0)
-
-        score <- (gap / c.h[candidates]) *
-          (1 + leftover / N.h[candidates])
-
-        j <- candidates[which.max(score)]
-
-        if (score[which.max(score)] <= 0) break
-
-        rounded_allocations[j] <- rounded_allocations[j] + 1L
-        remaining_budget <- cost - sum(rounded_allocations * c.h)
-      }
-
-    } else if (!is.null(variance)) { # Precision-constrained
-      propNum <- sum(N.h * S.h * sqrt(c.h))
-      propDen <- variance * sum(N.h)**2 + sum(N.h * S.h**2)
-      raw_allocations <- N.h * S.h / sqrt(c.h) * propNum / propDen
-
-      adjusted_allocations <- .adjust_precision_constrained(N.h, S.h, c.h, variance, lbound)
-
-      rounded_allocations <- .round_precision_constrained(adjusted_allocations, N.h, S.h, c.h, variance, lbound)
-    }
-  }
-
-
-  #Note: use the largest remainder (of extra sample left) [NOT NECESSARILY HOW MUCH IS CURRENTLY ALLOCATED)
-
-  # Prep for outputting
-  if (allocation == "proportional") {
-    inputs <- list("N.h" = N.h)
-  } else if (allocation == "power") {
-    inputs <- list("N.h" = N.h, "power" = power)
-  } else if (allocation == "neyman") {
-    inputs <- list("N.h" = N.h, "S.h" = S.h)
-  } else if (allocation == "optimal" & !is.null(cost)) {
-    inputs <- list("N.h" = N.h, "S.h" = S.h, "c.h" = c.h, "cost" = cost)
-  } else if (allocation == "optimal" & !is.null(variance)) {
-    inputs <- list("N.h" = N.h, "S.h" = S.h, "c.h" = c.h, "variance" = variance)
-  }
-
-  output <- list()
-  out.length <- 0
-  if (any(outputs == "raw")){
-    out.length <- out.length + 1
-    output[out.length] <- list(raw_allocations)
-    names(output)[out.length] <- "raw"
-  }
-  if (any(outputs == "adjusted")){
-    out.length <- out.length + 1
-    output[out.length] <- list(adjusted_allocations)
-    names(output)[out.length] <- "adjusted"
-  }
-  if (any(outputs == "rounded")){
-    out.length <- out.length + 1
-    output[out.length] <- list(rounded_allocations)
-    names(output)[out.length] <- "rounded"
-  }
-
-  #output <- as.integer(rounded_allocations)
-  if (allocation == "optimal") {
-    n.print <- sum(rounded_allocations)
-    if (!is.null(c.h)){
-      actual_cost <- sum(rounded_allocations * c.h)
-      n.print <- n.print |>
-        paste0(" (sample cost: ",round(actual_cost,digits=1),")")
-    }
-  } else {
-    n.print <- n.samp
-  }
-  message(paste0("Sample allocation of ", n.print, " using ", allocation, " with the relevant inputs:"))
-  for (i in 1:length(inputs)) {
-    message(paste0("  ",
-                   names(inputs)[i],
-                   " = ",
-                   paste0(inputs[[i]], collapse = ", "),
-                   collapse = ""
-    ))
-  }
-  message()
-  message("Output:")
-  for (i in 1:length(output)) {
-    message(paste0(" ",
-                   names(output)[i],
-                   " = ",
-                   paste0(output[[i]], collapse = ", "),
-                   collapse = ""
-    ))
-  }
-
-  if (length(output) == 1){
-    output <- output[[1]]
-  }
-
-  return(output)
-}
+#' Sample allocation
+#'
+#' Compute the proportional, power, Neyman, and optimal sample allocations.
+#'
+#' @section Method:
+#' The \emph{allocate} function allocates a sample size \emph{n} on \emph{H} strata using one of the following allocation methods:
+#' \enumerate{
+#'   \item Proportional allocation \[\code{n.samp, N.h, allocation = "proportional"}\]
+#'         \deqn{n_h = n \times \frac{N_h}{\sum\limits_{h=1}^H N_h}}
+#'         where \cr
+#'         \eqn{n}: total sample size to be allocated (function input is \code{n.samp}), and \cr
+#'         \eqn{N_h}: population size of stratum \emph{h} (function input is \code{N.h}).
+#'   \item Power allocation \[\code{n.samp, N.h, power, allocation = "power"}\]
+#'         \deqn{n_h = n \times \frac{N_h^\alpha}{\sum\limits_{h=1}^H N_h^\alpha}}
+#'          where \cr
+#'          \eqn{\alpha}: a power value to control over-under-sampling with \eqn{0 \le \alpha \le 1} (function input is \code{power}).
+#'   \item Neyman allocation \[\code{n.samp, N.h, S.h, allocation = "neyman"}\]
+#'         \deqn{n_h = n \times \frac{N_h S_h}{\sum\limits_{h=1}^H N_h S_h}}
+#'         where \cr
+#'         \eqn{S_h}:  standard deviation of stratum \emph{h} (function input is \code{S.h}).
+#'    \item Optimal allocation
+#'           \itemize{
+#'             \item cost-constrained \[\code{N.h, S.h, c.h, cost, allocation = "optimal"}\]
+#'                   \deqn{n_h = (C−c_0) \times \frac{N_h S_h / \sqrt{c_h}}{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}}
+#'             where \cr
+#'             \eqn{c_h}: cost per unit in stratum \emph{h} (function input is \code{c.h}), and \cr
+#'             \eqn{(C - c_0)}: total variable cost (function input is \code{cost})
+#'             \item precision-constrained \[\code{N.h, S.h, c.h, variance, allocation = "optimal"}\]
+#'                   \deqn{n_h = N_h S_h / \sqrt{c_h} \times \frac{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}{V_0 \left(\sum\limits_{h=1}^H  N_h \right)^2 + \sum\limits_{h=1}^H N_h S_h^2}}
+#'             where \cr
+#'             \eqn{V_0}: fixed variance target for estimated mean (function input is \code{variance})
+#'           }
+#' }
+#'
+#' The table below presents the relevant inputs for each type; when irrelevant inputs are entered, an error message will be displayed.
+#'
+#' \tabular{lllllllll}{
+#'   \strong{allocation} \tab \strong{N.h} \tab \strong{n.samp} \tab \strong{S.h} \tab \strong{c.h} \tab \strong{cost} \tab \strong{variance} \tab \strong{lbound} \tab \strong{power} \cr
+#'   proportional                   \tab X \tab X \tab  \tab  \tab  \tab  \tab X \tab \cr
+#'   power                          \tab X \tab X \tab  \tab  \tab  \tab  \tab X \tab X\cr
+#'   neyman                         \tab X \tab X \tab X \tab X \tab  \tab  \tab X \tab \cr
+#'   optimal: cost-constrained      \tab X \tab  \tab X \tab X \tab X \tab  \tab X \tab \cr
+#'   optimal: precision-constrained \tab X \tab  \tab X \tab X \tab  \tab X \tab X \tab
+#' }
+#'
+#' @param allocation type of allocation, must be one of \code{"proportional"}, \code{"power"}, \code{"neyman"}, or \code{"optimal"}.
+#' @param n.samp total sample size to be allocated (positive integer of length 1). \cr\cr
+#'        required for the following allocation types: proportional, power, and Neyman, and \code{NULL} otherwise.
+#' @param N.h vector of population stratum sizes (\eqn{N_h}, all positive values), for example \code{c(150, 600, 250)}. \cr\cr
+#' required for all allocation types.
+#' @param S.h vector of stratum unit standard deviations (positive values same length as \code{N.h}) (\eqn{S_h}). \cr\cr
+#' required for the following allocation types: Neyman, and optimal, and \code{NULL} otherwise.
+#' @param c.h vector of cost per unit in stratum h (positive values same length as \code{N.h}) (\eqn{c_h}). \cr\cr
+#' required for the optimal allocation only, and \code{NULL} otherwise.
+#' @param cost total variable cost (positive value) \eqn{(C - c_0)}. \cr\cr
+#' required for the cost-constrained optimal allocation only, and \code{NULL} otherwise.
+#' @param variance fixed variance target for estimated mean (positive value) (\eqn{V_0}). \cr\cr
+#' required for the precision-constrained optimal allocation only, and \code{NULL} otherwise.
+#' @param power power value for power allocation (\eqn{0 \le \alpha \le 1}). \cr\cr
+#' required for the power allocation only, and \code{NULL} otherwise.
+#' @param lbound minimum stratum-level (positive integer of length 1). Default value is 2.
+#' @param outputs
+#' character vector representing whether to output:\cr
+#'  \enumerate{
+#'    \item the raw allocations before accounting for N.h, lbound, and n.samp if needed \[\code{"raw"}\],
+#'    \item the continuous version after accounting for the above \[\code{"adjusted"}\], and/or
+#'    \item the rounded version of the above \[\code{"rounded"}\]
+#'  } \cr
+#'  Default is to only return the final rounded version \[\code{"rounded"}\]. \cr
+#'  If one version is requested, the result will be a numeric vector. Otherwise, the result will be a named list matching the requested outputs.
+#' @return If one output type is requested, a numeric vector of allocations.
+#' If multiple output types are requested, a named list containing the requested
+#' allocation vectors.
+#' @export
+#' @examples
+#' # The first step is getting a frame summary
+#' #  Summarize the IPEDS dataset by OBEREG
+#' # - N: number of universities per region
+#' # - SD_ENRTOT: standard deviation of total enrollment per region
+#' # - Filter out rows with missing ENRTOT to ensure accurate variance estimates
+#'
+#' ipeds_summary <- ipeds |>
+#'   tidytable::filter(!is.na(ENRTOT)) |>
+#'   tidytable::group_by(OBEREG) |>
+#'   tidytable::summarize(
+#'     N = tidytable::n(),
+#'     SD_ENRTOT = stats::sd(ENRTOT)
+#'   ) |>
+#'   tidytable::ungroup()
+#'
+#' # Example of proportional allocation
+#' ipeds_summary |>
+#'   tidytable::mutate(
+#'     n = allocate("proportional", N.h = N, n.samp = 500)
+#'   )
+#'
+#' # Example of power allocation
+#' ipeds_summary |>
+#'   tidytable::mutate(
+#'     n = allocate("power", N.h = N, power = 0.5, n.samp = 500)
+#'   )
+#'
+#' # Example of Neyman allocation
+#' ipeds_summary |>
+#'   tidytable::mutate(
+#'     n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT)
+#'   )
+#'
+#' # Example of Neyman allocation with a lower bound of 5
+#' ipeds_summary |>
+#'   tidytable::mutate(
+#'     n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT, lbound = 5)
+#'   )
+
+allocate <- function(allocation, N.h, n.samp = NULL, S.h = NULL, c.h = NULL, cost = NULL, variance = NULL, power = NULL, lbound = 2, outputs = "rounded") {
+  allocation <- match.arg(allocation, c("proportional", "power", "neyman", "optimal"))
+  outputs <- match.arg(outputs, c("raw", "adjusted", "rounded"), several.ok = TRUE)
+
+  ######
+  # Check inputs
+  .problems <- NULL # Initialize list of problems found with inputs
+  .addProblem <- function(parameter, condition, problems = .problems, allocation = NULL) { # Function to simplify addition to problems found with inputs to our running list (.problems)
+    if (is.null(parameter)) { # No parameter given
+      problem <- condition
+    } else if (parameter %in% c("allocation", "n.samp", "N.h", "S.h", "c.h", "cost", "variance", "power", "lbound")) {
+      problem <- paste0("The ", parameter, " parameter ", condition)
+    }
+    problemsNew <- c(problems, problem)
+    return(problemsNew)
+  }
+
+  .condition <- paste0('must be specified for allocation=="', allocation, '"')
+  # n.samp parameter
+  if (allocation %in% c("proportional", "power", "neyman") & is.null(n.samp)) {
+    .problems <- .addProblem(parameter = "n.samp", condition = .condition)
+  } else if (!allocation %in% c("proportional", "power", "neyman") & !is.null(n.samp)) {
+    warning("The n.samp parameter should be NULL", call. = FALSE)
+  }
+  # N.h parameter
+  if (allocation %in% c("proportional", "power", "neyman", "optimal") & is.null(N.h)) {
+    .problems <- .addProblem(parameter = "N.h", condition = .condition)
+  }
+  # S.h parameter
+  if (allocation %in% c("neyman", "optimal") & is.null(S.h)) {
+    .problems <- .addProblem(parameter = "S.h", condition = .condition)
+  } else if (!allocation %in% c("neyman", "optimal") & !is.null(S.h)) {
+    warning("The S.h parameter should be NULL", call. = FALSE)
+  }
+  # c.h parameter
+  if (allocation %in% c("optimal") & is.null(c.h)) {
+    .problems <- .addProblem(parameter = "c.h", condition = .condition)
+  } else if (!allocation %in% c("optimal") & !is.null(c.h)) {
+    warning("The c.h parameter should be NULL", call. = FALSE)
+  }
+  # power parameter
+  if (allocation %in% c("power") & is.null(power)) {
+    .problems <- .addProblem(parameter = "power", condition = .condition)
+  } else if (!allocation %in% c("power") & !is.null(power)) {
+    warning("The power parameter should be NULL", call. = FALSE)
+  }
+  # optimal allocation: only 1 of cost or variance should be provided
+  if (allocation == "optimal" & sum(is.null(cost), is.null(variance)) != 1) {
+    .problems <- .addProblem(parameter = NULL, condition = paste0('Exactly one of the cost and variance parameters should be supplied for allocation=="', allocation, '"'))
+  } else if (!allocation %in% c("optimal") & !(is.null(cost) & is.null(variance))) {
+    if (!is.null(cost)) {
+      warning("The cost parameter should be NULL", call. = FALSE)
+    }
+    if (!is.null(variance)) {
+      warning("The variance parameter should be NULL", call. = FALSE)
+    }
+  }
+
+  ###
+  # Check for miscellaneous unexpected parameter values
+
+  # n.samp
+  .condition <- "must be a positive integer of length 1"
+  if (!is.null(n.samp) & allocation %in% c("proportional", "power", "neyman")) {
+    if (!all(length(n.samp) == 1 & (typeof(n.samp) %in% c("integer") | (typeof(n.samp) == "double" & round(n.samp) == n.samp)) & n.samp > 0)) {
+      .problems <- .addProblem(parameter = "n.samp", condition = .condition)
+    }
+  }
+  # lbound
+  if (!(length(lbound) == 1 & (typeof(lbound) %in% c("integer") | (typeof(lbound) == "double" & round(lbound) == lbound)) & lbound > 0)) {
+    .problems <- .addProblem(parameter = "lbound", condition = .condition)
+  }
+
+  .condition <- "must be a vector of positive values (integers or non-integers)"
+  # N.h
+  if (!is.null(N.h)) {
+    if (!(length(N.h) >= 1 & typeof(N.h) %in% c("integer", "double") & all(N.h > 0))) {
+      .problems <- .addProblem(parameter = "N.h", condition = .condition)
+    }
+  }
+  .condition <- paste0(.condition, " that are the same length as N.h")
+  # S.h
+  if (!is.null(S.h) & allocation %in% c("neyman", "optimal")) {
+    if (!(length(S.h) >= 1 & typeof(S.h) %in% c("integer", "double") & all(S.h > 0) & length(S.h) == length(N.h))) {
+      .problems <- .addProblem(parameter = "S.h", condition = .condition)
+    }
+  }
+  # c.h
+  if (!is.null(c.h) & allocation %in% c("optimal")) {
+    if (!(length(c.h) >= 1 & typeof(c.h) %in% c("integer", "double") & all(c.h > 0) & length(c.h) == length(N.h))) {
+      .problems <- .addProblem(parameter = "c.h", condition = .condition)
+    }
+  }
+  # cost
+  .condition <- "must be a positive value (integer or non-integer)"
+  if (!is.null(cost) & allocation %in% "optimal") {
+    if (!all(length(cost) == 1 & typeof(cost) %in% c("integer", "double") & cost > 0)) {
+      .problems <- .addProblem(parameter = "cost", condition = .condition)
+    }
+  }
+  # variance
+  if (allocation == "optimal" & !is.null(variance)) {
+    if (!all(length(variance) == 1 & typeof(variance) %in% c("integer", "double") & variance > 0 & length(variance) == 1)) {
+      .problems <- .addProblem(parameter = "variance", condition = .condition)
+    }
+  }
+
+  # power parameter
+  .condition <- "must be a positive value between 0 and 1, inclusive"
+  if (allocation == "power" & !is.null(power)) {
+    if (!(length(power) == 1 & typeof(power) %in% c("integer", "double") & 0 <= power & power <= 1)) {
+      .problems <- .addProblem(parameter = "power", condition = .condition)
+    }
+  }
+
+
+  if (allocation %in% c("proportional", "power", "neyman")) {
+    if (!(is.null(lbound) | is.null(N.h) | is.null(n.samp))) {
+      if (!all(lbound * length(N.h) <= n.samp)) {
+        .problems <- c(.problems, "lbound*length(N.h) must be less than or equal to n.samp")
+      }
+    }
+  }
+
+  if (!(is.null(N.h) | is.null(n.samp))) {
+    if (length(.problems) == 0) {
+      if (sum(N.h) < n.samp) {
+        warning("sum(N.h) is less than n.samp")
+      }
+    }
+  }
+
+  if (allocation %in% c("optimal")) {
+    if (!is.null(cost) & length(.problems) == 0) {
+      if (sum(lbound * c.h) > cost) {
+        .problems <- c(.problems, "sum(lbound*c.h) must be less than or equal to cost")
+      }
+    }
+  }
+
+  if (allocation %in% c("proportional", "power", "neyman","optimal")){
+    if (!is.null(N.h)){
+      if (any(lbound > N.h)){
+        warning("lbound > N.h for at least one stratum")
+      }
+    }
+  }
+
+  # outputs parameter
+  .condition <- 'must be one or more of: "raw", "adjusted", "rounded"'
+  if (!all(outputs %in% c("raw", "adjusted", "rounded"))) {
+    .problems <- .addProblem(parameter = "outputs", condition = .condition)
+  }
+
+  ###
+  # Aggregate problems and stop the program if necessary
+  if (length(.problems) > 0) {
+    if (length(.problems) > 1) {
+      .problems <- paste0(
+        1:length(.problems),
+        ": ",
+        .problems
+      )
+    }
+    stop(
+      "\n",
+      paste0(.problems, collapse = "\n")
+    )
+  }
+
+  ######
+  #Fixed-total continuous adjustment:
+  # * Respects lower bound (lbound)
+  # * Respects upper bound (N.h)
+  # * Sums to sample size (n.samp)
+  # * Unrounded sample sizes
+  .adjust_fixed_total <- function(weights, n.samp, lbound, N.h) {
+    h <- length(weights)
+
+    lbound.h <- rep(lbound, h) #Extending lbound to a vector
+
+    # Raw proportional target
+    raw <- n.samp * weights / sum(weights)
+
+    #If raw already respects both bounds, use it directly.
+    if (all(raw >= lbound.h) && all(raw <= N.h)) {
+      return(as.numeric(raw))
+    }
+
+    #Otherwise fall back to iterative capacity-and-floor enforcement
+    alloc <- lbound.h
+    remaining <- n.samp - sum(alloc)
+    capacity <- N.h - alloc #Amount left in N.h after allocating lbound
+
+    if (remaining == 0) {
+      return(as.numeric(alloc))
+    }
+
+    #Initialize strata variables
+    active <- rep(TRUE, h) #Does stratum still need to be processed
+    extra <- rep(0, h) #How much to allocate to stratum above raw_allocations
+
+
+    repeat {
+      w <- weights
+      w[!active] <- 0
+
+      if (sum(w) <= 0) {
+        stop("No feasible solution: cannot distribute remaining sample.")
+      }
+
+      proposed <- remaining * w / sum(w) #Proposed additional allocation
+      hit <- active & (proposed > capacity) #Assigned all of stratum
+
+      #If didn't hit any statum's upper bound, break here
+      if (!any(hit)) {
+        extra[active] <- remaining * weights[active] / sum(weights[active])
+        break
+      }
+
+      #If hit a stratum's N, allocate its remaining capacity and set inactive
+      alloc[hit] <- alloc[hit] + capacity[hit]
+      remaining <- n.samp - sum(alloc)
+      active[hit] <- FALSE
+      capacity[hit] <- 0
+    }
+
+
+    adjusted_allocations <- alloc + extra
+    adjusted_allocations
+  }
+
+  ######
+  # Calculate stratified mean variance
+
+  .stratified_mean_variance <- function(n.h, N.h, S.h, N = sum(N.h)) {
+    W.h <- N.h / N
+    sum(W.h^2 * (1 - n.h / N.h) * S.h^2 / n.h)
+  }
+
+  ######
+  # Precision-constrained continuous adjustment:
+  # * Respects lower bound (lbound)
+  # * Respects upper bound (N.h)
+  # * Hits variance target
+  # * Unrounded sample sizes
+
+  .adjust_precision_constrained <- function(N.h, S.h, c.h, variance, lbound){
+    h <- length(N.h)
+    N <- sum(N.h)
+
+    a.h <- N.h * S.h / sqrt(c.h)
+
+    lbound.h <- rep(lbound, h)
+
+    # Initialize strata variables
+    adjusted_allocations <- rep(NA_real_, h)
+    fixed <- rep(FALSE, h) #Treat stratum fixed (TRUE) vs. free (FALSE)
+
+
+    V.lower <- .stratified_mean_variance(lbound.h, N.h, S.h)
+    V.upper <- .stratified_mean_variance(N.h, N.h, S.h)
+
+    if (variance >= V.lower) {
+      warning("The lower-bound allocation already satisfies the variance target; returning lower-bound allocation.")
+      return(lbound.h)
+    }
+
+    if (variance < V.upper) {
+      stop("No feasible allocation: even a census of all strata does not achieve the variance target.")
+    }
+
+
+    repeat {
+      # Compute variance contribution from fixed strata
+      if (any(fixed)) {
+        V.fixed <- .stratified_mean_variance(
+          n.h = adjusted_allocations[fixed],
+          N.h = N.h[fixed],
+          S.h = S.h[fixed],
+          N = N
+        )
+      } else {
+        V.fixed <- 0
+      }
+
+      # If all strata are fixed, break
+      if (all(fixed)) {
+        V.total <- .stratified_mean_variance(adjusted_allocations, N.h, S.h)
+        if (abs(V.total - variance) > 1e-10) {
+          warning("No feasible allocation exactly meets the variance target under the imposed bounds.")
+        }
+        break
+      }
+
+      # Recompute scaling constant over free strata
+      free <- !fixed
+      A.free <- sum(N.h[free] * S.h[free] * sqrt(c.h[free]))
+      B.free <- sum(N.h[free] * S.h[free]^2)
+
+      denom <- N^2 * (variance - V.fixed) + B.free
+
+      if (denom <= 0) {
+        stop("No feasible allocation solution for the requested variance after applying bounds.")
+      }
+
+      k.free <- A.free / denom
+
+      # Candidate allocations for free strata
+      candidate <- k.free * a.h[free]
+
+      # Check for violations
+      below <- candidate < lbound.h[free]
+      above <- candidate > N.h[free]
+
+      if (!any(below | above)){
+        adjusted_allocations[free] <- candidate
+        break
+      }
+
+      # Fix violating strata
+      free_idx <- which(free)
+
+      if (any(below)) {
+        idx <- free_idx[below]
+        adjusted_allocations[idx] <- lbound.h[idx]
+        fixed[idx] <- TRUE
+      }
+
+      if (any(above)) {
+        idx <- free_idx[above]
+        adjusted_allocations[idx] <- N.h[idx]
+        fixed[idx] <- TRUE
+      }
+
+    }
+    return(adjusted_allocations)
+  }
+
+  ######
+  # Precision-constrained rounding
+
+  .round_precision_constrained <- function(adjusted_allocations, N.h, S.h, c.h, variance, lbound) {
+    rounded_allocations <- floor(adjusted_allocations + 1e-9)
+    rounded_allocations <- pmax(rounded_allocations, lbound)
+    rounded_allocations <- pmin(rounded_allocations, N.h)
+    rounded_allocations <- as.integer(rounded_allocations)
+
+    frac <- adjusted_allocations - floor(adjusted_allocations + 1e-9)
+
+    while (.stratified_mean_variance(rounded_allocations, N.h, S.h) > variance) {
+      candidates <- which(rounded_allocations < N.h)
+
+      if (length(candidates) == 0) {
+        warning("No feasible integer allocation found to satisfy variance target.")
+        break
+      }
+
+      current_var <- .stratified_mean_variance(rounded_allocations, N.h, S.h)
+      scores <- rep(-Inf, length(N.h))
+
+      for (h in candidates) {
+        trial <- rounded_allocations
+        trial[h] <- trial[h] + 1L
+        new_var <- .stratified_mean_variance(trial, N.h, S.h)
+        var_gain <- current_var - new_var
+
+        #leftover relative to adjusted continuous allocation
+        leftover <- pmax(N.h[h] - adjusted_allocations[h], 0)
+
+        #normalize leftover
+        leftover_factor <- 1 + leftover / N.h[h]
+
+        #fractional preference
+        frac_factor <- 1 + frac[h]
+
+        #combined score
+        scores[h] <- (var_gain / c.h[h]) * leftover_factor * frac_factor
+      }
+
+      best_h <- which.max(scores)
+      if (!is.finite(scores[best_h]) || scores[best_h] <= 0) {
+        warning("Could not improve enough to satisfy variance constraint.")
+        break
+      }
+
+      rounded_allocations[best_h] <- rounded_allocations[best_h] + 1L
+    }
+
+    return(rounded_allocations)
+  }
+
+
+  ######
+  # Rounding fixed-totals
+
+  .round_fixed_total <- function(adjusted_allocations, n.samp, N.h, lbound, score) {
+    rounded_allocations <- floor(adjusted_allocations + 1e-9)
+    rounded_allocations <- pmax(rounded_allocations, lbound)
+    rounded_allocations <- pmin(rounded_allocations, N.h)
+    rounded_allocations <- as.integer(rounded_allocations)
+
+    deficit <- n.samp - sum(rounded_allocations)
+
+    if (deficit < 0) {
+      stop("Initial rounded allocation exceeds the target total.")
+    }
+    while (deficit > 0) {
+      candidates <- which(rounded_allocations < N.h)
+      if (length(candidates) == 0) {
+        stop("No feasible way to allocate remaining units.")
+      }
+
+
+      scores <- score(rounded_allocations)
+      scores[-candidates] <- -Inf
+
+
+      j <- candidates[which.max(scores[candidates])]
+
+
+      if (!is.finite(scores[j])) {
+        stop("No feasible scoring candidate for additional allocation.")
+      }
+
+      rounded_allocations[j] <- rounded_allocations[j] + 1L
+      deficit <- deficit - 1L
+    }
+
+    return(rounded_allocations)
+
+  }
+
+  .score_proportional <- function(adjusted_allocations, rounded_allocations, N.h) {
+    gap <- pmax(adjusted_allocations - rounded_allocations, 0)
+    leftover <- pmax(N.h - adjusted_allocations, 0)
+    weight <- N.h / sum(N.h)
+
+    gap * (1 + leftover / N.h) * (1 + weight)
+  }
+
+  .score_power <- function(adjusted_allocations, rounded_allocations, N.h, power) {
+    gap <- pmax(adjusted_allocations - rounded_allocations, 0)
+    leftover <- pmax(N.h - adjusted_allocations, 0)
+    weight_raw <- N.h^power
+    weight <- weight_raw / sum(weight_raw)
+
+    gap * (1 + leftover / N.h) * (1 + weight)
+  }
+
+  .score_neyman <- function(adjusted_allocations, rounded_allocations, N.h, S.h) {
+    gap <- pmax(adjusted_allocations - rounded_allocations, 0)
+    leftover <- pmax(N.h - adjusted_allocations, 0)
+    weight_raw <- N.h * S.h
+    weight <- weight_raw / sum(weight_raw)
+
+    gap * (1 + leftover / N.h) * (1 + weight)
+  }
+
+
+
+  ######
+  # Moving onto the actual allocation
+  N <- sum(N.h)
+  if (allocation == "proportional") {
+    raw_allocations <- n.samp * N.h / N
+
+    weights <- N.h
+
+    #Respect lower-bound, upper bound, and n.samp
+    adjusted_allocations <- .adjust_fixed_total(
+      weights,
+      n.samp,
+      lbound,
+      N.h
+    )
+
+    rounded_allocations <- .round_fixed_total(
+      adjusted_allocations,
+      n.samp,
+      N.h,
+      lbound,
+      score = function(current) .score_proportional(adjusted_allocations, current, N.h)
+    )
+  } else if (allocation == "power") {
+    N.h.powered <- N.h^power
+    raw_allocations <- n.samp * N.h.powered / sum(N.h.powered)
+
+    weights <- N.h.powered
+
+    #Respect lower-bound, upper bound, and n.samp
+    adjusted_allocations <- .adjust_fixed_total(
+      weights,
+      n.samp,
+      lbound,
+      N.h
+    )
+
+    rounded_allocations <- .round_fixed_total(
+      adjusted_allocations,
+      n.samp,
+      N.h,
+      lbound,
+      score = function(current) .score_power(adjusted_allocations, current, N.h, power)
+    )
+
+
+  } else if (allocation == "neyman") {
+    propNum <- N.h * S.h # Numerator
+    propDen <- sum(propNum) # Denominator
+    raw_allocations <- n.samp * propNum / propDen
+
+    weights <- propNum
+
+    #Respect lower-bound, upper bound, and n.samp
+    adjusted_allocations <- .adjust_fixed_total(
+      weights,
+      n.samp,
+      lbound,
+      N.h
+    )
+
+    rounded_allocations <- .round_fixed_total(
+      adjusted_allocations,
+      n.samp,
+      N.h,
+      lbound,
+      score = function(current) .score_neyman(adjusted_allocations, current, N.h, S.h)
+    )
+  } else if (allocation == "optimal") {
+    if (!is.null(cost)) { # Cost-constrained
+      propNum <- N.h * S.h / sqrt(c.h)
+      propDen <- sum(N.h * S.h * sqrt(c.h))
+      raw_allocations <- cost * propNum / propDen
+
+      #Lower-bound adjusted version
+      lbound.h <- rep(lbound, length(N.h))
+      baseline_cost <- sum(lbound.h * c.h)
+      adjusted_target <- pmax(raw_allocations, lbound.h)
+      adjusted_target <- pmin(adjusted_target, N.h)
+
+      adjusted_cost <- sum(adjusted_target * c.h)
+
+      if (adjusted_cost > cost){
+        #Rescale only the excess above lbound to stay on budget
+        excess <- adjusted_target - lbound.h
+        excess_cost <- sum(excess * c.h)
+        rho <- if (excess_cost > 0) (cost - baseline_cost) / excess_cost else 0
+        rho <- max(min(rho, 1), 0)
+        adjusted_allocations <- lbound.h + rho * excess
+      } else {
+        adjusted_allocations <- adjusted_target
+      }
+
+      #Deterministic rounded version under budget
+      rounded_allocations <- floor(adjusted_allocations + 1e-9)
+      rounded_allocations <- pmax(rounded_allocations, lbound.h)
+      rounded_allocations <- pmin(rounded_allocations, N.h)
+      rounded_allocations <- as.integer(rounded_allocations)
+
+      remaining_budget <- cost - sum(rounded_allocations * c.h)
+      frac <- adjusted_allocations - floor(adjusted_allocations + 1e-9)
+
+      repeat {
+        candidates <- which(
+          rounded_allocations < N.h &
+            (c.h <= remaining_budget)
+        )
+        if (length(candidates) == 0) break
+
+        # prioritize largest fractional parts per unit cost
+        gap <- pmax(adjusted_allocations[candidates] - rounded_allocations[candidates], 0)
+        leftover <- pmax(N.h[candidates] - adjusted_allocations[candidates], 0)
+
+        score <- (gap / c.h[candidates]) *
+          (1 + leftover / N.h[candidates])
+
+        j <- candidates[which.max(score)]
+
+        if (score[which.max(score)] <= 0) break
+
+        rounded_allocations[j] <- rounded_allocations[j] + 1L
+        remaining_budget <- cost - sum(rounded_allocations * c.h)
+      }
+
+    } else if (!is.null(variance)) { # Precision-constrained
+      propNum <- sum(N.h * S.h * sqrt(c.h))
+      propDen <- variance * sum(N.h)**2 + sum(N.h * S.h**2)
+      raw_allocations <- N.h * S.h / sqrt(c.h) * propNum / propDen
+
+      adjusted_allocations <- .adjust_precision_constrained(N.h, S.h, c.h, variance, lbound)
+
+      rounded_allocations <- .round_precision_constrained(adjusted_allocations, N.h, S.h, c.h, variance, lbound)
+    }
+  }
+
+
+  #Note: use the largest remainder (of extra sample left) [NOT NECESSARILY HOW MUCH IS CURRENTLY ALLOCATED)
+
+  # Prep for outputting
+  if (allocation == "proportional") {
+    inputs <- list("N.h" = N.h)
+  } else if (allocation == "power") {
+    inputs <- list("N.h" = N.h, "power" = power)
+  } else if (allocation == "neyman") {
+    inputs <- list("N.h" = N.h, "S.h" = S.h)
+  } else if (allocation == "optimal" & !is.null(cost)) {
+    inputs <- list("N.h" = N.h, "S.h" = S.h, "c.h" = c.h, "cost" = cost)
+  } else if (allocation == "optimal" & !is.null(variance)) {
+    inputs <- list("N.h" = N.h, "S.h" = S.h, "c.h" = c.h, "variance" = variance)
+  }
+
+  output <- list()
+  out.length <- 0
+  if (any(outputs == "raw")){
+    out.length <- out.length + 1
+    output[out.length] <- list(raw_allocations)
+    names(output)[out.length] <- "raw"
+  }
+  if (any(outputs == "adjusted")){
+    out.length <- out.length + 1
+    output[out.length] <- list(adjusted_allocations)
+    names(output)[out.length] <- "adjusted"
+  }
+  if (any(outputs == "rounded")){
+    out.length <- out.length + 1
+    output[out.length] <- list(rounded_allocations)
+    names(output)[out.length] <- "rounded"
+  }
+
+  #output <- as.integer(rounded_allocations)
+  if (allocation == "optimal") {
+    n.print <- sum(rounded_allocations)
+    if (!is.null(c.h)){
+      actual_cost <- sum(rounded_allocations * c.h)
+      n.print <- n.print |>
+        paste0(" (sample cost: ",round(actual_cost,digits=1),")")
+    }
+  } else {
+    n.print <- n.samp
+  }
+  message(paste0("Sample allocation of ", n.print, " using ", allocation, " with the relevant inputs:"))
+  for (i in 1:length(inputs)) {
+    message(paste0("  ",
+                   names(inputs)[i],
+                   " = ",
+                   paste0(inputs[[i]], collapse = ", "),
+                   collapse = ""
+    ))
+  }
+  message()
+  message("Output:")
+  for (i in 1:length(output)) {
+    message(paste0(" ",
+                   names(output)[i],
+                   " = ",
+                   paste0(output[[i]], collapse = ", "),
+                   collapse = ""
+    ))
+  }
+
+  if (length(output) == 1){
+    output <- output[[1]]
+  }
+
+  return(output)
+}
diff --git a/man/allocate.Rd b/man/allocate.Rd
index a7426c6..fafb361 100644
--- a/man/allocate.Rd
+++ b/man/allocate.Rd
@@ -13,12 +13,12 @@ allocate(
   cost = NULL,
   variance = NULL,
   power = NULL,
-  lbound = 2
+  lbound = 2,
+  outputs = "rounded"
 )
 }
 \arguments{
-\item{allocation}{type of allocation,
-must be one of \code{"proportional"}, \code{"power"}, \code{"neyman"}, or \code{"optimal"}.}
+\item{allocation}{type of allocation, must be one of \code{"proportional"}, \code{"power"}, \code{"neyman"}, or \code{"optimal"}.}
 
 \item{N.h}{vector of population stratum sizes (\eqn{N_h}, all positive values), for example \code{c(150, 600, 250)}. \cr\cr
 required for all allocation types.}
@@ -26,16 +26,14 @@ required for all allocation types.}
 \item{n.samp}{total sample size to be allocated (positive integer of length 1). \cr\cr
 required for the following allocation types: proportional, power, and Neyman, and \code{NULL} otherwise.}
 
-\item{S.h}{vector of stratum unit standard deviations
-(positive values same length as \code{N.h}) (\eqn{S_h}). \cr\cr
+\item{S.h}{vector of stratum unit standard deviations (positive values same length as \code{N.h}) (\eqn{S_h}). \cr\cr
 required for the following allocation types: Neyman, and optimal, and \code{NULL} otherwise.}
 
 \item{c.h}{vector of cost per unit in stratum h (positive values same length as \code{N.h}) (\eqn{c_h}). \cr\cr
 required for the optimal allocation only, and \code{NULL} otherwise.}
 
-\item{cost}{total variable cost (positive value) \eqn{(C – c_0)}. \cr\cr
-required for the cost-constrained optimal
-allocation only, and \code{NULL} otherwise.}
+\item{cost}{total variable cost (positive value) \eqn{(C - c_0)}. \cr\cr
+required for the cost-constrained optimal allocation only, and \code{NULL} otherwise.}
 
 \item{variance}{fixed variance target for estimated mean (positive value) (\eqn{V_0}). \cr\cr
 required for the precision-constrained optimal allocation only, and \code{NULL} otherwise.}
@@ -43,11 +41,21 @@ required for the precision-constrained optimal allocation only, and \code{NULL}
 \item{power}{power value for power allocation (\eqn{0 \le \alpha \le 1}). \cr\cr
 required for the power allocation only, and \code{NULL} otherwise.}
 
-\item{lbound}{minimum stratum-level
-(positive integer of length 1). Default value is 2.}
+\item{lbound}{minimum stratum-level (positive integer of length 1). Default value is 2.}
+
+\item{outputs}{character vector representing whether to output:\cr
+\enumerate{
+\item the raw allocations before accounting for N.h, lbound, and n.samp if needed [\code{"raw"}],
+\item the continuous version after accounting for the above [\code{"adjusted"}], and/or
+\item the rounded version of the above [\code{"rounded"}]
+} \cr
+Default is to only return the final rounded version [\code{"rounded"}]. \cr
+If one version is requested, the result will be a numeric vector. Otherwise, the result will be a named list matching the requested outputs.}
 }
 \value{
-Integer vector same length of N.h final allocation
+If one output type is requested, a numeric vector of allocations.
+If multiple output types are requested, a named list containing the requested
+allocation vectors.
 }
 \description{
 Compute the proportional, power, Neyman, and optimal sample allocations.
@@ -75,7 +83,7 @@ where \cr
 \deqn{n_h = (C−c_0) \times \frac{N_h S_h / \sqrt{c_h}}{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}}
             where \cr
 \eqn{c_h}: cost per unit in stratum \emph{h} (function input is \code{c.h}), and \cr
-\eqn{(C – c_0)}: total variable cost (function input is \code{cost})
+\eqn{(C - c_0)}: total variable cost (function input is \code{cost})
 \item precision-constrained [\code{N.h, S.h, c.h, variance, allocation = "optimal"}]
 \deqn{n_h = N_h S_h / \sqrt{c_h} \times \frac{\sum\limits_{h=1}^H N_h S_h \sqrt{c_h}}{V_0 \left(\sum\limits_{h=1}^H  N_h \right)^2 + \sum\limits_{h=1}^H N_h S_h^2}}
 where \cr
@@ -85,22 +93,53 @@ where \cr
 
 The table below presents the relevant inputs for each type; when irrelevant inputs are entered, an error message will be displayed.
 
-\tabular{lllll}{
-\strong{allocation} \tab \strong{N.h} \tab \strong{n.h} \tab \strong{power} \tab \strong{lbound} \cr
-proportional \tab ✓ \tab ✓  \tab   \tab ✓\cr
-power        \tab ✓ \tab ✓  \tab ✓ \tab ✓
+\tabular{lllllllll}{
+\strong{allocation} \tab \strong{N.h} \tab \strong{n.samp} \tab \strong{S.h} \tab \strong{c.h} \tab \strong{cost} \tab \strong{variance} \tab \strong{lbound} \tab \strong{power} \cr
+proportional                   \tab X \tab X \tab  \tab  \tab  \tab  \tab X \tab \cr
+power                          \tab X \tab X \tab  \tab  \tab  \tab  \tab X \tab X\cr
+neyman                         \tab X \tab X \tab X \tab X \tab  \tab  \tab X \tab \cr
+optimal: cost-constrained      \tab X \tab  \tab X \tab X \tab X \tab  \tab X \tab \cr
+optimal: precision-constrained \tab X \tab  \tab X \tab X \tab  \tab X \tab X \tab
 }
 }
 
 \examples{
+# The first step is getting a frame summary
+#  Summarize the IPEDS dataset by OBEREG
+# - N: number of universities per region
+# - SD_ENRTOT: standard deviation of total enrollment per region
+# - Filter out rows with missing ENRTOT to ensure accurate variance estimates
+
+ipeds_summary <- ipeds |>
+  tidytable::filter(!is.na(ENRTOT)) |>
+  tidytable::group_by(OBEREG) |>
+  tidytable::summarize(
+    N = tidytable::n(),
+    SD_ENRTOT = stats::sd(ENRTOT)
+  ) |>
+  tidytable::ungroup()
+
 # Example of proportional allocation
 ipeds_summary |>
   tidytable::mutate(
     n = allocate("proportional", N.h = N, n.samp = 500)
   )
+
 # Example of power allocation
 ipeds_summary |>
   tidytable::mutate(
     n = allocate("power", N.h = N, power = 0.5, n.samp = 500)
   )
+
+# Example of Neyman allocation
+ipeds_summary |>
+  tidytable::mutate(
+    n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT)
+  )
+
+# Example of Neyman allocation with a lower bound of 5
+ipeds_summary |>
+  tidytable::mutate(
+    n = allocate("neyman", N.h = N, n.samp = 500, S.h = SD_ENRTOT, lbound = 5)
+  )
 }