Answer to: R code for a constrained randomization in a three-dimensional array?

Starting with this: set.seed(42) ary <- `dim<-`(sample(rep(c("A","B","C"), times=c(152,152,88))), c(4,7,14)) ### Each pipe (X by Y combination) should contain roughly the same number of samples apply(ary, 1:2, table) |> apply(1, mean) # A B C # 5.428571 5.428571 3.142857 ### Each depth level (Z coordinate) should also contain roughly the same number of samples apply(ary, 3, table) |> apply(1, mean) # A B C # 10.857143 10.857143 6.285714 However, I think you may be under-constraining here: while this meets your "average" constraints, it may not meet reasonable distributions: ### Each pipe ... apply(ary, 1:2, table) |> apply(1, range) # A B C # [1,] 1 1 1 # [2,] 10 9 6 ### Each depth level ... apply(ary, 3, table) |> apply(1, range) # A B C # [1,] 9 6 3 # [2,] 15 14 9 This means that while you do average around 5.4 for each of A and B, they both range from 1-9 within any given pipe. One can work around this by randomly iterating. Let's say you'd be willing to accept ranges of 3-8 for A and B; similarly, enforcing ranges of (say) 8-16 for the depth-level tables. Side note: the apply code worked above, but occasionally when one or more of the treatments are missing in a particular pipe, then it won't work. I'll create a table2 that is more immune to this problem. set.seed(43) ary <- `dim<-`(sample(rep(c("A","B","C"), times=c(152,152,88))), c(4,7,14)) apply(ary, 1:2, table) |> apply(1, mean) # [1] NA NA NA NA Warning in mean.default(newX[, i], ...) : # argument is not numeric or logical: returning NA # Warning in mean.default(newX[, i], ...) : # argument is not numeric or logical: returning NA # Warning in mean.default(newX[, i], ...) : # argument is not numeric or logical: returning NA # Warning in mean.default(newX[, i], ...) : # argument is not numeric or logical: returning NA table2 <- function(x, always = c("A", "B", "C")) table(c(x, always)) - 1L apply(ary, 1:2, table2) |> apply(1, mean) # A B C # 5.428571 5.428571 3.142857 Back to the loop: test1 <- function(ar) { x <- apply(ar, 1:2, table2) |> apply(1, range) all(x[1, c("A", "B")] >= 3 & x[2, c("A", "B")] <= 8) } test2 <- function(ar) { x <- apply(ar, 3, table) |> apply(1, range) all(x[1, c("A", "B")] >= 8 & x[2, c("A", "B")] <= 16) } while (TRUE) { ary <- `dim<-`(sample(rep(c("A","B","C"), times=c(152,152,88))), c(4,7,14)) if (!test1(ary)) next if (!test2(ary)) next break } apply(ary, 1:2, table2) |> apply(1, mean) # A B C # 5.428571 5.428571 3.142857 apply(ary, 1:2, table2) |> apply(1, range) # A B C # [1,] 3 3 1 # [2,] 8 8 6 apply(ary, 3, table2) |> apply(1, mean) # A B C # 10.857143 10.857143 6.285714 apply(ary, 3, table2) |> apply(1, range) # A B C # [1,] 8 8 2 # [2,] 15 16 9 Those ranges look better. Note that with this approach, the more constraints you add, the more likely you are to over-constrain into an unachievable combination. The loop above took 6-ish seconds, very achievable. Making the additional tests tighter will take a bit longer to achieve a feasible array. This template/approach may help you to achieve what you need: if you add constraints, add simple test functions to do a simple test, then add it to the sequence of if (!testn(ary)) next expressions before the break. At its worst, this busy-cycles the cpu looking for feasible solutions, it is easy to break out of this with C-c or the stop button in an IDE. If I want to do this programmatically and guard against infinite loops, I'll add a meta-constraint maxiters <- 1e6 (or such) and change the loop to be iter <- 0; while (iter < maxiters) { iter <- iter + 1; ... }, ensuring that it will break out of the while loop eventually. If you're curious, my resulting array: ary <- structure(c("B", "A", "A", "A", "B", "A", "C", "A", "A", "B", "A", "B", "A", "C", "A", "C", "A", "B", "B", "A", "B", "B", "B", "C", "A", "B", "C", "A", "C", "A", "C", "A", "B", "A", "C", "A", "B", "B", "A", "B", "B", "B", "C", "B", "B", "B", "B", "A", "B", "C", "C", "B", "A", "B", "C", "A", "B", "C", "A", "A", "A", "B", "A", "B", "C", "A", "B", "A", "B", "B", "A", "C", "B", "A", "B", "B", "B", "C", "C", "B", "C", "A", "C", "B", "A", "A", "A", "B", "B", "B", "B", "B", "B", "C", "B", "C", "C", "B", "B", "B", "B", "B", "A", "A", "C", "A", "C", "A", "A", "A", "B", "A", "A", "A", "A", "B", "C", "A", "A", "B", "B", "C", "B", "B", "B", "A", "B", "A", "A", "A", "B", "B", "C", "C", "B", "A", "A", "C", "C", "A", "C", "A", "C", "B", "A", "A", "B", "A", "A", "B", "A", "B", "A", "A", "A", "B", "B", "A", "C", "B", "B", "B", "C", "B", "C", "B", "C", "A", "C", "C", "B", "A", "A", "A", "B", "C", "A", "B", "A", "C", "B", "B", "C", "B", "C", "B", "A", "A", "A", "C", "B", "C", "A", "C", "B", "B", "A", "C", "C", "A", "B", "B", "B", "B", "B", "A", "C", "B", "B", "B", "A", "C", "A", "A", "A", "C", "A", "B", "C", "A", "B", "A", "A", "A", "A", "B", "A", "A", "B", "B", "A", "A", "A", "C", "B", "B", "B", "B", "B", "A", "B", "B", "B", "B", "C", "B", "A", "B", "A", "B", "A", "B", "A", "B", "A", "B", "C", "B", "A", "A", "C", "B", "B", "C", "B", "A", "A", "C", "C", "A", "C", "B", "A", "B", "A", "C", "A", "C", "A", "C", "C", "A", "B", "C", "A", "C", "A", "B", "A", "C", "A", "A", "B", "B", "A", "C", "C", "B", "C", "C", "A", "A", "B", "B", "B", "B", "B", "C", "C", "B", "B", "C", "B", "B", "B", "A", "A", "A", "B", "C", "C", "B", "C", "A", "A", "C", "A", "A", "B", "B", "A", "C", "B", "C", "A", "A", "B", "A", "A", "A", "A", "B", "B", "A", "B", "B", "A", "C", "A", "A", "A", "C", "B", "C", "C", "C", "A", "A", "A", "A", "B", "C", "B", "B", "A", "A", "A", "B", "A", "A", "B", "A", "A", "C", "C", "A", "A", "A", "B", "A", "A", "B", "B", "B", "B", "A", "A", "B", "B", "A", "B", "C"), dim = c(4L, 7L, 14L))