这是开发版本包中的超快单衬RcppAlgos.
devtools::install_github("jwood000/RcppAlgos")
library(RcppAlgos)
myPerms <– permuteGeneral(3,6,TRUE,"prod","==",36) - 1L
myPerms
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 1 1 2 2
[2,] 0 0 1 2 1 2
[3,] 0 0 1 2 2 1
[4,] 0 0 2 1 1 2
[5,] 0 0 2 1 2 1
[6,] 0 0 2 2 1 1
.
.
.
[,1] [,2] [,3] [,4] [,5] [,6]
[85,] 2 2 0 0 1 1
[86,] 2 2 0 1 0 1
[87,] 2 2 0 1 1 0
[88,] 2 2 1 0 0 1
[89,] 2 2 1 0 1 0
[90,] 2 2 1 1 0 0
以下是一些基准rcppAlgo, r2eOne, r2eTwo, and OPFun是每个方法的代码的函数包装器。
microbenchmark(rcppAlgo(),r2eOne(),r2eTwo(),OPFun(N=6) unit = "relative")
Unit: relative
expr min lq mean median uq max neval
rcppAlgo() 1.00000 1.00000 1.00000 1.0000 1.00000 1.000000 100
r2eOne() 471.56007 473.26487 194.01669 267.9402 274.46604 8.373630 100
r2eTwo() 50.71091 48.84173 24.01617 27.8441 34.02326 2.044374 100
OPFun(N=6) 37.35899 24.38966 22.38029 19.7059 19.51935 31.18059 100
解释
由于OP正在寻找具有特定频率的数字的特定组合，我们可以利用算术基本定理 https://en.wikipedia.org/wiki/Fundamental_theorem_of_arithmetic，它指出每个数字都可以写成素数的唯一组合的乘积。我们得到了集合0, 1, 2加 1 给出集合1, 2, 3。我们这样做是为了避免在购买产品时出现很多零。
现在，我们的任务是找到所有组合，使每个元素恰好出现两次。这意味着当我们将产品应用于我们的目标组合后，我们得到1*1*2*2*3*3 = 36 (N.B. 1不是质数，但可以忽略，因为1*n = n for all n）。现在问题很简单了。
我们只需找到所有组合，使得产品等于36然后减去1回到我们原来的一组数字，瞧！
通用解决方案
下面，我们有一个通用解决方案，可用于查找给定向量的所有排列，其中每个元素重复特定次数。
library(RcppAlgos) ## for primeSieve and permuteGeneral
MakePerms <- function(v, numReps, myCap = NULL) {
m <- sum(numReps)
n <- length(v)
## Generate some primes using prime
## number theorem; fudging a bit to
## ensure we get n-1 prime numbers
myPs <- primeSieve(2*n*log(n))[1:(n-1)]
## Set up vector that will be tested
myV <- c(1L, myPs)
target <- prod(myV^numReps)
ps <- permuteGeneral(myV, m, TRUE, "prod", "==", target, myCap)
for (j in 1:n) {ps[ps == myV[j]] <- v[j]}
ps
}
它在很大程度上依赖于算术基本定理的素因式分解的唯一性和一点索引（不像上面的小例子那么简单，但仍然只有 7 行，而且仍然超级快）。
我们首先创建第一个向量n-1素数并粘上1去完成myV。然后我们提出每个元素myV每个元素所需的重复次数由下式给出numReps并拿走产品以获得我们的target价值。这里有些例子：
-
v = c(10,13,267,1) and numReps = c(3,1,2,5)-->>myV = c(1,2,3,5)-->>target = 1^3 * 2^1 * 3^2 * 5^5 = 56250
-
v = 0:5 and numReps = c(1,2,1,2,2,2)-->>myV = c(1,2,3,5,7,11)-->>target = 1^1 * 2^2 * 3^1 * 5^2 * 7^2 * 11^2 = 1778700
- 操作示例：
v = c(0,1,2) and numReps = c(2,2,2)-->>myV = c(1,2,3)-->>target = 1^2 * 2^2 * 3^2 = 36
在我们找到所有其乘积等于的排列之后target值，我们只需映射原始向量的内容v使用索引生成的矩阵。
例如，如果您设置N = 8在OP的例子中你得到了所有排列c(0,1,2) with 0准确地重复4次，以及1 and 2重复两次。
t1 <- OPFun(N=8)
t2 <- MakePerms(0:2, c(4,2,2))
all.equal(t1[do.call(order, as.data.frame(t1)), ],
t2[do.call(order, as.data.frame(t2)), ])
[1] TRUE
microbenchmark(fun2(8), MakePerms(0:2, c(4,2,2)), unit = "relative")
Unit: relative
expr min lq mean median uq max neval
OPFun(8) 23.25099 22.56178 18.64762 19.52436 18.37387 10.90934 100
MakePerms(0:2, c(4, 2, 2)) 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 100
应该注意的是，可能的排列数量快速增长，因此像这样的尝试MakePerms(0:5, rep(2, 6))将失败，因为 的排列总数0:5 12 times is 12^6 = 2,985,984 > 2^31 - 1（即矩阵的最大行数Rcpp）。然而，我们并不期望所有这些都符合我们的标准，所以如果我们设定一个上限，比如说10^7，我们就会成功。观察：
a <- MakePerms(0:5, rep(2, 6), 10^7)
nrow(a)
7484400
set.seed(17)
a[sample(nrow(a), 10), ]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 0 5 3 3 1 2 4 4 5 1 0 2
[2,] 5 4 2 1 1 0 3 4 5 2 3 0
[3,] 2 4 5 3 5 1 3 0 1 0 4 2
[4,] 4 3 3 1 2 5 0 5 4 1 0 2
[5,] 2 2 5 3 4 1 0 3 5 1 0 4
[6,] 3 1 1 5 0 3 2 0 2 4 4 5
[7,] 1 1 4 2 0 5 4 0 3 5 3 2
[8,] 1 0 4 2 4 2 5 1 3 0 5 3
[9,] 4 3 4 1 5 0 0 2 2 1 3 5
[10,] 1 0 5 3 2 0 1 4 3 4 2 5
指某东西的用途myCap也可以大大提高效率。
microbenchmark(withOutCap = MakePerms(0:5, c(1,2,1,2,1,2)),
withCap = MakePerms(0:5, c(1,2,1,2,1,2), 10^5),
times = 15)
Unit: milliseconds
expr min lq mean median uq max neval
withOutCap 219.64847 246.4718 275.04672 282.52829 299.33816 311.2031 15
withCap 22.56437 30.6904 33.30469 31.70443 37.50858 41.6095 15
identical(MakePerms(0:5, c(1,2,1,2,1,2)), MakePerms(0:5, c(1,2,1,2,1,2), 10^5))
[1] TRUE
iterpc解决方案
似乎到目前为止提供的答案都是严格学术性的，因为@StéphaneLaurent 提供的答案要优越得多。超级通用，一条线，而且超级快！！
microbenchmark(iter = getall(iterpc(c(2,2,2), labels=c(0,1,2), ordered=TRUE)),
rcppAlg = MakePerms(0:2, c(2,2,2)))
Unit: microseconds
expr min lq mean median uq max neval
iter 428.885 453.2975 592.53164 540.154 683.9585 1165.772 100
rcppAlg 62.418 74.5205 93.44926 81.749 108.4660 216.454 100
随着排列数量的增加，故事也会发生变化。观察：
microbenchmark(iter = getall(iterpc(c(2,2,2,2), labels=c(0,1,2,3), ordered=TRUE)),
rcppAlg = MakePerms(0:3, c(2,2,2,2)),
rcppAlgCap = MakePerms(0:3, c(2,2,2,2), 5000))
Unit: microseconds
expr min lq mean median uq max neval
iter 877.246 1052.7060 1394.636 1150.0895 1265.088 8914.980 100
rcppAlg 964.446 1449.7115 2084.944 1787.9350 1906.242 10921.156 100
如果您使用myCap, MakePerms是有点快。这并不重要，因为随着iterpc解决方案，您甚至不必考虑会得到多少结果。很不错！！
Update
新版本的RcppAlgos（我是该书的作者）刚刚在 CRAN 上发布。现在还有一个额外的论据permuteGeneral called freqs它允许多重集的排列，这正是OP正在寻找的。
microbenchmark(iter = getall(iterpc(c(2,2,2,2), labels=0:3, ordered=TRUE)),
newRcppAlgos = permuteGeneral(0:3, freqs = c(2,2,2,2)))
Unit: microseconds
expr min lq mean median uq max neval
iter 457.442 482.8365 609.98678 508.6150 572.581 4037.048 100
newRcppAlgos 33.159 43.3975 56.40026 48.5665 58.194 625.691 100
microbenchmark(iter = getall(iterpc(c(5,4,3,2), labels=0:3, ordered=TRUE)),
newRcppAlgos = permuteGeneral(0:3, freqs = c(5,4,3,2)))
Unit: milliseconds
expr min lq mean median uq max neval
iter 480.25976 552.54343 567.9155 565.23066 579.0258 751.8556 100
newRcppAlgos 83.41194 87.03957 104.6279 95.67596 107.3572 181.1119 100
identical(getall(iterpc(c(5,4,3,2), labels=0:3, ordered=TRUE)),
permuteGeneral(0:3, freqs = c(5,4,3,2)))
[1] TRUE
nrow(permuteGeneral(0:3, freqs = c(5,4,3,2)))
[1] 2522520
Update 2
正如@StéphaneLaurent 指出的那样，该包arrangements已作为替代品发布iterpc（参见@RandyLai 的评论）。它效率更高，能够处理更广泛的组合问题（例如分区）。以下是较大示例的基准：
microbenchmark(arrangements = permutations(x = 0:3, freq = c(5,4,3,2)),
RcppAlgos = permuteGeneral(0:3, freqs = c(5,4,3,2)))
Unit: milliseconds
expr min lq mean median uq max neval
arrangements 97.10078 98.67154 113.5953 100.56261 131.3244 163.8912 100
RcppAlgos 92.13122 93.84818 108.1845 95.72691 101.2647 165.7248 100
...几乎相同的结果。
一个巨大的好处是arrangements是一次获得一个排列（或成块）的能力getnext。这允许用户生成超过2^31 - 1结果并总体上提供了更大的灵活性。
有关此类问题的更多信息，请参见R，我写了一个广泛的概述 https://stackoverflow.com/a/47983855/4408538对于这个问题：R：有/无替换以及不同/非不同项目/多重集的排列和组合 https://stackoverflow.com/q/22569176/4408538.
