## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE)

## -----------------------------------------------------------------------------
library(randomizr)

set.seed(20260824)
sims <- 2000

# An exact promise: the property either held on every draw or it did not.
exact <- function(draw, property) {
  ok <- replicate(sims, isTRUE(property(draw())))
  list(result = sprintf("holds on %d of %d draws", sum(ok), sims),
       holds = all(ok))
}

# A probability promise: no unit's realized rate may sit more than five Monte
# Carlo standard errors from the probability the package claims for it.
probability <- function(draw, target, condition = 1) {
  realized <- rowMeans(replicate(sims, draw()) == condition)
  gap <- max(abs(realized - target))
  tolerance <- 5 * sqrt(max(target * (1 - target)) / sims)
  list(result = sprintf("largest gap %.3f, tolerance %.3f", gap, tolerance),
       holds = gap < tolerance)
}

## -----------------------------------------------------------------------------
N <- 100
blocks <- rep(c("a", "b", "c", "d", "e"), each = 20)
clusters <- rep(1:20, each = 5)
regions <- rep(c("north", "south"), each = 50)
p_unit <- seq(0.05, 0.95, length.out = N)
x <- rnorm(N)

# helpers the properties below are written in terms of
by_group <- function(z, g) tapply(z, g, sum)
uniform_within <- function(z, g) all(tapply(z, g, function(x) length(unique(x))) == 1)

## -----------------------------------------------------------------------------
checks <- list(
  list("simple_ra()", "each unit is treated with probability 0.25",
       probability(function() simple_ra(N, prob = 0.25), 0.25)),

  list("simple_ra()", "the number treated is random, not fixed",
       local({
         m <- replicate(sims, sum(simple_ra(N, prob = 0.25)))
         list(result = sprintf("sd %.2f against a binomial %.2f, %d distinct counts",
                               sd(m), sqrt(N * 0.25 * 0.75), length(unique(m))),
              holds = abs(sd(m) / sqrt(N * 0.25 * 0.75) - 1) < 0.1)
       })),

  list("complete_ra()", "each unit is treated with probability 25/100",
       probability(function() complete_ra(N, m = 25), 0.25)),

  list("complete_ra()", "exactly 25 units are treated",
       exact(function() complete_ra(N, m = 25), function(z) sum(z) == 25)),

  list("complete_ra(num_arms = 3)", "the three arms get 33, 33 and 34 units",
       exact(function() complete_ra(N, num_arms = 3),
             function(z) all(sort(as.vector(table(z))) == c(33, 33, 34)))),

  list("complete_ra(num_arms = 3)", "each unit reaches arm T1 with probability 1/3",
       probability(function() complete_ra(N, num_arms = 3), 1 / 3, condition = "T1")),

  list("block_ra()", "exactly half of every block of 20 is treated",
       exact(function() block_ra(blocks = blocks, prob = 0.5),
             function(z) all(by_group(z, blocks) == 10))),

  list("block_ra()", "each unit is treated with its block's probability, 0.5",
       probability(function() block_ra(blocks = blocks, prob = 0.5), 0.5)),

  list("cluster_ra()", "no cluster is ever split across conditions",
       exact(function() cluster_ra(clusters = clusters),
             function(z) uniform_within(z, clusters))),

  list("cluster_ra()", "each of the 20 clusters is treated with probability 0.5",
       probability(function() tapply(cluster_ra(clusters = clusters), clusters, `[`, 1), 0.5)),

  list("block_and_cluster_ra()", "clusters stay intact and each region treats 25 units",
       exact(function() block_and_cluster_ra(clusters = clusters, blocks = regions),
             function(z) uniform_within(z, clusters) && all(by_group(z, regions) == 25))),

  list("balanced_ra()", "each unit is treated with its own probability",
       probability(function() balanced_ra(prob_unit = p_unit), p_unit)),

  list("balanced_ra()", "the number treated is the floor or ceiling of its target",
       exact(function() balanced_ra(prob_unit = p_unit),
             function(z) sum(z) %in% c(floor(sum(p_unit)), ceiling(sum(p_unit))))),

  list("balanced_ra(formula = ~ x)", "each unit is treated with its own probability",
       probability(function() balanced_ra(prob_unit = p_unit, formula = ~ x), p_unit)),

  list("balanced_ra(formula = ~ x)", "the number treated is the floor or ceiling of its target",
       exact(function() balanced_ra(prob_unit = p_unit, formula = ~ x),
             function(z) sum(z) %in% c(floor(sum(p_unit)), ceiling(sum(p_unit))))),

  list("declare_ra() + conduct_ra()", "realized rates match the declared probabilities",
       local({
         declaration <- declare_ra(
           blocks = blocks,
           block_m_each = rbind(c(15, 5), c(14, 6), c(13, 7), c(12, 8), c(11, 9)))
         probability(function() conduct_ra(declaration),
                     declaration$probabilities_matrix[, 2])
       })),

  list("complete_rs()", "exactly 25 units are sampled",
       exact(function() complete_rs(N, n = 25), function(s) sum(s) == 25)),

  list("strata_rs()", "exactly half of every stratum is sampled",
       exact(function() strata_rs(strata = blocks, prob = 0.5),
             function(s) all(by_group(s, blocks) == 10)))
)

guarantees <- data.frame(
  Function = vapply(checks, function(x) x[[1]], character(1)),
  Promise = vapply(checks, function(x) x[[2]], character(1)),
  Measured = vapply(checks, function(x) x[[3]]$result, character(1)),
  Holds = vapply(checks, function(x) x[[3]]$holds, logical(1))
)

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(guarantees, row.names = FALSE)

## -----------------------------------------------------------------------------
all(guarantees$Holds)