diff --git a/R/allocate.R b/R/allocate.R
index c0936fa..649c2e1 100644
--- a/R/allocate.R
+++ b/R/allocate.R
@@ -21,7 +21,7 @@
 #'    \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}}}
+#'                   \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})
@@ -59,7 +59,18 @@
 #' @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
+#' @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
@@ -100,8 +111,10 @@
 #'   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) {
+
+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
@@ -187,7 +200,7 @@ allocate <- function(allocation, N.h, n.samp = NULL, S.h = NULL, c.h = NULL, cos
     }
   }
   # c.h
-  if (!is.null(c.h) & allocation %in% c("neyman", "optimal")) {
+  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)
     }
@@ -234,10 +247,25 @@ allocate <- function(allocation, N.h, n.samp = NULL, S.h = NULL, c.h = NULL, cos
   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)")
+        .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) {
@@ -254,84 +282,435 @@ allocate <- function(allocation, N.h, n.samp = NULL, S.h = NULL, c.h = NULL, cos
     )
   }
 
+  ######
+  #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") {
-    allocations <- n.samp * N.h / N
+    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
-    N.powered <- sum(N.h.powered)
-    allocations <- n.samp * N.h.powered / N.powered
+    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
-    allocations <- n.samp * propNum / propDen
+    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))
-      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
+      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)
 
-  final_n <- max(ceiling(sum(allocations)), length(N.h) * lbound)
+      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
+      }
 
-  sizes <- allocations
+      #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)
 
-  num_groups <- length(N.h)
+      remaining_budget <- cost - sum(rounded_allocations * c.h)
+      frac <- adjusted_allocations - floor(adjusted_allocations + 1e-9)
 
-  # Ensure each allocation is at least lbound
-  adjusted_allocations <- pmax(allocations, lbound)
+      repeat {
+        candidates <- which(
+          rounded_allocations < N.h &
+            (c.h <= remaining_budget)
+        )
+        if (length(candidates) == 0) break
 
-  # Calculate the sum of the adjusted allocations
-  total_allocated <- sum(adjusted_allocations)
+        # 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)
 
-  # Calculate the difference from the total n
-  difference <- final_n - total_allocated
+        score <- (gap / c.h[candidates]) *
+          (1 + leftover / N.h[candidates])
 
-  # If difference is positive, distribute it proportionally
-  if (difference != 0) {
-    remaining_sizes <- sizes
-    remaining_sizes <- adjusted_allocations - lbound
-    remaining_total <- sum(remaining_sizes)
+        j <- candidates[which.max(score)]
 
-    if (remaining_total > 0) {
-      additional_allocations <- difference * remaining_sizes / remaining_total
-    } else {
-      additional_allocations <- rep(0, num_groups)
-    }
+        if (score[which.max(score)] <= 0) break
 
-    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
+        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)
@@ -343,28 +722,60 @@ allocate <- function(allocation, N.h, n.samp = NULL, S.h = NULL, c.h = NULL, cos
     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)
+
+  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 <- final_n
+    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 = ""
+                   names(inputs)[i],
+                   " = ",
+                   paste0(inputs[[i]], collapse = ", "),
+                   collapse = ""
     ))
   }
   message()
   message("Output:")
-  message(paste0(outputs,
-    collapse = ", "
-  ))
-  return(outputs)
+  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..fafb361 100644
--- a/man/allocate.Rd
+++ b/man/allocate.Rd
@@ -13,7 +13,8 @@ allocate(
   cost = NULL,
   variance = NULL,
   power = NULL,
-  lbound = 2
+  lbound = 2,
+  outputs = "rounded"
 )
 }
 \arguments{
@@ -41,9 +42,20 @@ required for the precision-constrained optimal allocation only, and \code{NULL}
 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{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.
@@ -68,8 +80,8 @@ 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})
 \item precision-constrained [\code{N.h, S.h, c.h, variance, allocation = "optimal"}]
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"
+  )
+})