Let's say I have 6 numbers:
a <- c(1,2,3,4,5,6)
I want to list all possible 3 digit combinations of these 6 numbers, including repetitions.
The desired result would look like this:
1 1 1
1 2 3
1 2 4
...
I do not want to include elements that have the same 3 numbers but in a different order:
for example
1 2 3
3 2 1
should exclude one of them
The combinations function from gtools can do this:
library(gtools)
combinations(n=6, r=3, v=a, repeats.allowed=TRUE)
[,1] [,2] [,3] [1,] 1 1 1 [2,] 1 1 2 [3,] 1 1 3 ... [54,] 5 5 6 [55,] 5 6 6 [56,] 6 6 6
expand.grid returns a data.frame of combinations, picking one from each set you supply it. If you 1 2 3 is not the same as 3 2 1, subset it to get just the rows you want.
df <- expand.grid(1:6, 1:6, 1:6)
df[df$Var1 <= df$Var2 & df$Var2 <= df$Var3,]
#
# Var1 Var2 Var3
# 1 1 1 1
# 37 1 1 2
# 43 1 2 2
# 44 2 2 2
# 73 1 1 3
# 79 1 2 3
# 80 2 2 3
# 85 1 3 3
# 86 2 3 3
# 87 3 3 3
# 109 1 1 4
# 115 1 2 4
# 116 2 2 4
# 121 1 3 4
# 122 2 3 4
# 123 3 3 4
# 127 1 4 4
# 128 2 4 4
# 129 3 4 4
# 130 4 4 4
# 145 1 1 5
# 151 1 2 5
# 152 2 2 5
# 157 1 3 5
# 158 2 3 5
# 159 3 3 5
# 163 1 4 5
# 164 2 4 5
# 165 3 4 5
# 166 4 4 5
# 169 1 5 5
# 170 2 5 5
# 171 3 5 5
# 172 4 5 5
# 173 5 5 5
# 181 1 1 6
# 187 1 2 6
# 188 2 2 6
# 193 1 3 6
# 194 2 3 6
# 195 3 3 6
# 199 1 4 6
# 200 2 4 6
# 201 3 4 6
# 202 4 4 6
# 205 1 5 6
# 206 2 5 6
# 207 3 5 6
# 208 4 5 6
# 209 5 5 6
# 211 1 6 6
# 212 2 6 6
# 213 3 6 6
# 214 4 6 6
# 215 5 6 6
# 216 6 6 6
Below, is a general function that allows you to specify constraints on the output. For example, I have had many situations to arise whereby I needed all n-tuples of a given set such that there product was less than a given bound. Before I wrote this function, I was forced to use combinations and search for those rows that met my condition. This took a lot of time and a lot of memory.
Combo <- function(n,r,v=1:n,li=10^8,fun1="prod",fun2="<",repeats.allowed=FALSE) {
## where fun1 is a general function such as "prod", "sum", "sd", etc.
## and fun2 is a comparison operator such as "<", "<=", ">", "==", etc.
myfun <- match.fun(FUN = fun1)
operator1 <- match.fun(FUN = fun2)
operator2 <- match.fun(FUN = fun2)
myv <- sort(v)
if (fun2 %in% c(">",">=")) {
myv <- rev(myv)
TheLim <- min(v)
} else {
TheLim <- max(v)
if (fun2 == "==") {
operator1 <- match.fun(FUN = "<=")
}
}
if (!repeats.allowed) {
m <- matrix(numeric(0),combinat::nCm(n,r),r)
v1 <- myv; n1 <- length(v); t <- TRUE; count <- 0L
while (t) {
t <- operator1(myfun(v1[1:r]),li)
while (t && length(v1)>=r) {
t_1 <- operator2(myfun(v1[1:r]),li)
if (t_1) {count <- count+1L; m[count,] <- v1[1:r]}
v1 <- v1[-r]
t <- operator1(myfun(v1[1:r],na.rm=TRUE),li)
}
if (t) {
s <- 1:length(v1)
mymax <- myv[n1-(r-s)]
t1 <- which(!v1==mymax)
if (length(t1)>0) {
e <- max(t1)
v1[e] <- myv[which(myv==v1[e])+1L]
v1 <- c(v1[1:e],myv[(which(myv==v1[e])+1L):n1])
} else {
return(m[!is.na(m[,1]),])
}
} else {
r1 <- r-1L
while (r1>=1L && !t) {
v1[r1] <- myv[which(myv==v1[r1])+1L]
if (v1[r1]==TheLim) {r1 <- r1-1L; next}
v1 <- c(v1[1:r1],myv[(which(myv==v1[r1])+1L):n1])
t <- operator1(myfun(v1[1:r],na.rm=TRUE),li) && length(v1)>=r
r1 <- r1-1L
}
if (!t) {return(m[!is.na(m[,1]),])}
}
}
} else {
MySet <- 1:n
for (i in 1:(r-1L)) {MySet <- sapply(1:n, function(x) sum(MySet[1:x]))}
m <- matrix(numeric(0),nrow=MySet[n],ncol=r)
v1 <- c(rep(myv[1], r),myv[2:n]); n1 <- length(v); t <- TRUE; count <- 0L
while (t) {
t <- operator1(myfun(v1[1:r]),li)
while (t && length(v1)>=r) {
t_1 <- operator2(myfun(v1[1:r]),li)
if (t_1) {count <- count+1L; m[count,] <- v1[1:r]}
v1 <- v1[-r]
t <- operator1(myfun(v1[1:r],na.rm=TRUE),li)
}
if (t) {
s <- 1:length(v1)
t1 <- which(!v1==TheLim)
if (length(t1)>0) {
e <- max(t1)
v1[e] <- myv[which(myv==v1[e])+1L]
tSize <- r - length(myv[1:e])
if (!v1[e]==TheLim) {
v1 <- c(v1[1:e],rep(v1[e],tSize),myv[(which(myv==v1[e])+1L):n1])
} else {
v1 <- c(v1[1:e],rep(v1[e],tSize))
}
} else {
return(m[!is.na(m[,1]),])
}
} else {
r1 <- r-1L
while (r1>=1L && !t) {
if (v1[r1]==TheLim) {r1 <- r1-1L; next}
v1[r1] <- myv[which(myv==v1[r1])+1L]
tSize <- r - length(myv[1:r1])
v1 <- c(v1[1:r1],rep(v1[r1],tSize),myv[(which(myv==v1[r1])+1L):n1])
t <- operator1(myfun(v1[1:r],na.rm=TRUE),li) && length(v1)>=r
r1 <- r1-1L
}
if (!t) {return(m[!is.na(m[,1]),])}
}
}
}
}
Below are some examples:
## return all 3-tuple combinations of 1 through 6 such
## that the PRODUCT is less than 10
> Combo(n=6, r=3, v=1:6, li=10, fun1="prod", fun2="<", repeats.allowed=TRUE)
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 1 2
. . .
[10,] 1 3 3
[11,] 2 2 2
## return all 3-tuple combinations of 1 through 6 such
## that the SUM is less than 10
> Combo(n=6, r=3, v=1:6, li=10, fun1="sum", fun2="<", repeats.allowed=TRUE)
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 1 2
[3,] 1 1 3
. . .
[20,] 2 3 3
[21,] 2 3 4
[22,] 3 3 3
Here are some cool example involving prime numbers:
> library(numbers)
> myps <- Primes(1000)
> system.time(t1 <- Combo(n=length(myps), r=3, v=myps, li=10^5, fun1="prod", fun2="<", repeats.allowed=TRUE))
user system elapsed
0.18 0.00 0.18
> nrow(t1)
[1] 13465
> set.seed(42)
> t1[sample(nrow(t1),5),]
[,1] [,2] [,3]
[1,] 13 31 197
[2,] 17 19 167
[3,] 2 131 227
[4,] 11 11 751
[5,] 5 31 151
> object.size(t1)
323360 bytes
> system.time(t2 <- combinations(n=length(myps), r=3, v=myps, repeats.allowed=TRUE))
user system elapsed
3.63 0.00 3.68
> nrow(t2)
[1] 804440
> system.time(t3 <- t2[which(sapply(1:nrow(t2), function(x) prod(t2[x,]) < 10^5)),])
user system elapsed
1.55 0.00 1.54
> nrow(t3)
[1] 13465
> object.size(t2)
19306760 bytes
As you can see, the Combo function is much faster and performed in one step, whereas the combinations/sapply duo is slow (more than 5 seconds) and two clunky steps. The Combo function also returns an object that is almost 60 times smaller.
Here is another cool example. Let's say you want to find all 3-tuples of the first 168 primes (i.e. primes < 1000) such that the standard deviation is less than 50. No problem (with the same setup as above):
> system.time(t1 <- Combo(n=length(myps), r=3, v=myps, li=50, fun1="sd", fun2="<", repeats.allowed=TRUE))
user system elapsed
1.49 0.00 1.48
> system.time(t3 <- t2[which(sapply(1:nrow(t2), function(x) sd(t2[x,]) < 50)),])
user system elapsed
19.89 0.00 19.89
> nrow(t1)
[1] 22906
> nrow(t3)
[1] 22906
> all(t3==t1)
[1] TRUE
It should be noted that all function combinations don't work. For example if you let fun1="sd" and fun2=">", the above code will return 0 matches. Cheers!
A simple way to do this is using three for loop. Does this result what you want?
for (x in 1:6) {
for (y in x:6) {
for (z in y:6) {
print(paste(x,y,z))
}
}
}
