R中N个元素与q个元素的组合

Question

我有N=6 元素和q=3 元素符号为0,1,2。

我想在所有可能的位置创建N=6元素的所有向量，其中2元素等于0、2元素等于1和2元素等于2 .

这些向量的数量等于combn(6,2)*combn(4,2)*combn(2,2)=90。

这是在矩阵F 中构造这些90 向量的代码：

N=6
x<-c(1:N)
#########################################
A<-combn(x,2)
B<-matrix(0,ncol(A),length(x)) 
for( i in 1:ncol(A) ) 
{ 
y<-rep(0,N) 
y[A[1:nrow(A),i]]<-1 
B[i,]<-y 
}
######################################
E<-matrix(0,nrow(B),length(x)-nrow(A))  
for( i in 1:nrow(B) ) 
{ 
q=0 
for( j in 1:ncol(B) ) 
{ 
if( B[i,j]!=1 )  
{ 
q=q+1 
E[i,q]<-j 
}}}
########################################
ASD<-combn(E[i,],2) 
F<-matrix(0,nrow(B)*ncol(ASD),length(x)) 
q=0 
for( i in 1:nrow(B) ) 
{ 
ASD<-combn(E[i,],2) 
for( j in 1:ncol(ASD) ) 
{ 
B[i,ASD[1:nrow(ASD),j]]<-2 
q=q+1 
F[q,]<-B[i,] 
B[i,ASD[1:nrow(ASD),j]]<-0 
}}

还有其他不那么复杂的方法吗？

Answer 1

这是来自包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 和 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 正在寻找具有特定频率的数字的特定组合，我们可以使用 Fundamental theorem of arithmetic，它指出每个数字都可以写为质数的唯一组合的乘积。我们得到集合0, 1, 2，加1得到集合1, 2, 3。我们这样做是为了避免在取产品时出现很多零。

现在，我们的任务是找到所有组合，以使每个元素恰好出现两次。这意味着在我们将产品应用于我们的目标组合后，我们得到1*1*2*2*3*3 = 36（注意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 值。以下是一些示例：

1. v = c(10,13,267,1) 和 numReps = c(3,1,2,5) -->> myV = c(1,2,3,5) -->> target = 1^3 * 2^1 * 3^2 * 5^5 = 56250
2. v = 0:5 和 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
3. OP 示例：v = c(0,1,2) 和 numReps = c(2,2,2) -->> myV = c(1,2,3) -->> target = 1^2 * 2^2 * 3^2 = 36

在我们找到乘积等于target 值的所有排列后，我们只需使用索引将原始向量v 的内容映射到生成的矩阵。

例如，如果您在 OP 的示例中设置 N = 8，您会得到 c(0,1,2) 的所有排列，其中 0 完全重复 4 次，并且 1 和 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 的排列总数是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 解决方案，您甚至不必考虑将获得多少结果。很不错！！

更新

RcppAlgos（我是作者）的新版本刚刚在 CRAN 上发布。 permuteGeneral 现在有一个额外的参数，称为 freqs，它允许多重集的排列，这正是 OP 正在寻找的。​​p>

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

更新 2

正如@StéphaneLaurent 所指出的，包arrangements 已作为iterpc 的替代品发布（参见@RandyLai 的cmets）。它效率更高，并且能够处理更广泛的组合问题（例如分区）。以下是较大示例的基准：

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 中此类问题的更多信息，我在问题中写了extensive overview：R: Permutations and combinations with/without replacement and for distinct/non-distinct items/multiset。

Answer 2

您可以使用iterpc 包：

> library(iterpc)
> I <- iterpc(c(2,2,2), labels=c(0,1,2), ordered=TRUE)
> getall(I)
      [,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
 [7,]    0    1    0    1    2    2
 [8,]    0    1    0    2    1    2
 [9,]    0    1    0    2    2    1
[10,]    0    1    1    0    2    2
[11,]    0    1    1    2    0    2
[12,]    0    1    1    2    2    0
[13,]    0    1    2    0    1    2
[14,]    0    1    2    0    2    1
[15,]    0    1    2    1    0    2
[16,]    0    1    2    1    2    0
[17,]    0    1    2    2    0    1
[18,]    0    1    2    2    1    0
[19,]    0    2    0    1    1    2
[20,]    0    2    0    1    2    1
[21,]    0    2    0    2    1    1
[22,]    0    2    1    0    1    2
[23,]    0    2    1    0    2    1
[24,]    0    2    1    1    0    2
[25,]    0    2    1    1    2    0
[26,]    0    2    1    2    0    1
[27,]    0    2    1    2    1    0
[28,]    0    2    2    0    1    1
[29,]    0    2    2    1    0    1
[30,]    0    2    2    1    1    0
[31,]    1    0    0    1    2    2
[32,]    1    0    0    2    1    2
[33,]    1    0    0    2    2    1
[34,]    1    0    1    0    2    2
[35,]    1    0    1    2    0    2
[36,]    1    0    1    2    2    0
[37,]    1    0    2    0    1    2
[38,]    1    0    2    0    2    1
[39,]    1    0    2    1    0    2
[40,]    1    0    2    1    2    0
[41,]    1    0    2    2    0    1
[42,]    1    0    2    2    1    0
[43,]    1    1    0    0    2    2
[44,]    1    1    0    2    0    2
[45,]    1    1    0    2    2    0
[46,]    1    1    2    0    0    2
[47,]    1    1    2    0    2    0
[48,]    1    1    2    2    0    0
[49,]    1    2    0    0    1    2
[50,]    1    2    0    0    2    1
[51,]    1    2    0    1    0    2
[52,]    1    2    0    1    2    0
[53,]    1    2    0    2    0    1
[54,]    1    2    0    2    1    0
[55,]    1    2    1    0    0    2
[56,]    1    2    1    0    2    0
[57,]    1    2    1    2    0    0
[58,]    1    2    2    0    0    1
[59,]    1    2    2    0    1    0
[60,]    1    2    2    1    0    0
[61,]    2    0    0    1    1    2
[62,]    2    0    0    1    2    1
[63,]    2    0    0    2    1    1
[64,]    2    0    1    0    1    2
[65,]    2    0    1    0    2    1
[66,]    2    0    1    1    0    2
[67,]    2    0    1    1    2    0
[68,]    2    0    1    2    0    1
[69,]    2    0    1    2    1    0
[70,]    2    0    2    0    1    1
[71,]    2    0    2    1    0    1
[72,]    2    0    2    1    1    0
[73,]    2    1    0    0    1    2
[74,]    2    1    0    0    2    1
[75,]    2    1    0    1    0    2
[76,]    2    1    0    1    2    0
[77,]    2    1    0    2    0    1
[78,]    2    1    0    2    1    0
[79,]    2    1    1    0    0    2
[80,]    2    1    1    0    2    0
[81,]    2    1    1    2    0    0
[82,]    2    1    2    0    0    1
[83,]    2    1    2    0    1    0
[84,]    2    1    2    1    0    0
[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

编辑 2018-04-28

iterpc 现在已弃用，取而代之的是 arrangements。

Answer 3

建议了两种方法：一种效率低下，一种更有效但费力。 （在这种情况下，我将“效率”等同于缩放，而不是代码量或执行时间。也就是说，只要你只创建 90 行，那么你就可以了。如果这是一个简化的问题并且您确实需要扩展到更大的矩阵，那么permutations 可能会超出内存和/或 R 的容量。）

这两种解决方案都比您的代码短一些。第一个比较清晰，只有4行代码；第二个确实有点晦涩（并且似乎进入了索引间接“inception”），但实际上仍然只需要 13 行代码。第二个可能会有所减少，但我的“游戏”时间用完了:-)

效率低下

一种方法是创建所有排列并过滤掉重复。只要您的“N”不会太大，这就会起作用。

library(gtools)
v <- rep(0:2, 2)
p <- permutations(6, 6)
p[] <- v[p]
p <- p[!duplicated(p),]

head(p)
#      [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    0    1    2    0    1    2
# [2,]    0    1    2    0    2    1
# [3,]    0    1    2    1    0    2
# [4,]    0    1    2    1    2    0
# [5,]    0    1    2    2    0    1
# [6,]    0    1    2    2    1    0
tail(p)
#       [,1] [,2] [,3] [,4] [,5] [,6]
# [85,]    2    2    0    1    0    1
# [86,]    2    2    0    1    1    0
# [87,]    2    2    0    0    1    1
# [88,]    2    2    1    0    0    1
# [89,]    2    2    1    0    1    0
# [90,]    2    2    1    1    0    0

验证每一行的每个元素都恰好有两个：

all(apply(p, 1, table) == 2)
# [1] TRUE

低效

一种效率较低（因此更加努力）的方法：使用combn(6,2) 和combn(4,2) 创建列索引矩阵，然后适当地分配“因子”。 （稍后会更有意义。）

（注意：我经常在转置矩阵方面更好地考虑这些问题；您可以轻松地做到这一点，只需调整代码以交换列/行。）

我们需要的是类似expand.grid 的东西，一次用于两列。所以我们将从较小的问题开始：

left2 <- t(combn(6, 2))
mid2 <- t(combn(4, 2))
left2
#       [,1] [,2]
#  [1,]    1    2
#  [2,]    1    3
#  [3,]    1    4
#  [4,]    1    5
#  [5,]    1    6
#  [6,]    2    3
#  [7,]    2    4
#  [8,]    2    5
#  [9,]    2    6
# [10,]    3    4
# [11,]    3    5
# [12,]    3    6
# [13,]    4    5
# [14,]    4    6
# [15,]    5    6
mid2
#      [,1] [,2]
# [1,]    1    2
# [2,]    1    3
# [3,]    1    4
# [4,]    2    3
# [5,]    2    4
# [6,]    3    4

现在，网格将在这两个矩阵的行索引上展开。

eg <- expand.grid(a = 1:15, b = 1:6)
head(eg)
#   a b
# 1 1 1
# 2 2 1
# 3 3 1
# 4 4 1
# 5 5 1
# 6 6 1
inds <- cbind(left2[eg$a,], mid2[eg$b,])
head(inds)
#      [,1] [,2] [,3] [,4]
# [1,]    1    2    1    2
# [2,]    1    3    1    2
# [3,]    1    4    1    2
# [4,]    1    5    1    2
# [5,]    1    6    1    2
# [6,]    2    3    1    2
inds[25,,drop=FALSE]
#      [,1] [,2] [,3] [,4]
# [1,]    3    4    1    3

这意味着，对于第 25 行，我们应该将第 3 列和第 4 列替换为第一个因子（例如 0）。然后，其余列中 (1,2,5,6)，我们应该将第 1 列和第 3 列替换为第二个因子（例如 1）。再说一遍，c(1,2,5,6)[c(1,3)] 等于将第 1 列和第 5 列替换为第二个值 (1)。 （第三个值2 将进入所有剩余的插槽。）

所以，要从上面找到c(1,2,5,6)，我们可以使用setdiff(1:6,...)：

afterleft2 <- t(apply(left2[eg$a,], 1, function(a) setdiff(1:6, a)))
head( afterleft2 )
#      [,1] [,2] [,3] [,4]
# [1,]    3    4    5    6
# [2,]    2    4    5    6
# [3,]    2    3    5    6
# [4,]    2    3    4    6
# [5,]    2    3    4    5
# [6,]    1    4    5    6
afterleft2[25,,drop=FALSE]
#      [,1] [,2] [,3] [,4]
# [1,]    1    2    5    6

让我们开始修复inds 第三和第四列。

inds[,3] <- afterleft2[ cbind(1:90, mid2[eg$b,1]) ]
inds[,4] <- afterleft2[ cbind(1:90, mid2[eg$b,2]) ]
head(inds)
#      [,1] [,2] [,3] [,4]
# [1,]    1    2    3    4
# [2,]    1    3    2    4
# [3,]    1    4    2    3
# [4,]    1    5    2    3
# [5,]    1    6    2    3
# [6,]    2    3    1    4
inds[25,,drop=FALSE]
#      [,1] [,2] [,3] [,4]
# [1,]    3    4    1    5

由此我们看到第 25 行有我们期望的“1”和“5”。

现在结束：

nr <- nrow(inds)
out <- matrix(nrow = nr, ncol = 6L)
out[cbind(1:nr,inds[,1])] <- 0L
out[cbind(1:nr,inds[,2])] <- 0L
out[cbind(1:nr,inds[,3])] <- 1L
out[cbind(1:nr,inds[,4])] <- 1L
head(out)
#      [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    0    0    1    1   NA   NA
# [2,]    0    1    0    1   NA   NA
# [3,]    0    1    1    0   NA   NA
# [4,]    0    1    1   NA    0   NA
# [5,]    0    1    1   NA   NA    0
# [6,]    1    0    0    1   NA   NA
out[25,,drop=FALSE]
#      [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    1   NA    0    0    1   NA

我上面提到的“剩余插槽”（对于第三个值）都是NA，按设计。

out[is.na(out)] <- 2L
head(out)
#      [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    0    0    1    1    2    2
# [2,]    0    1    0    1    2    2
# [3,]    0    1    1    0    2    2
# [4,]    0    1    1    2    0    2
# [5,]    0    1    1    2    2    0
# [6,]    1    0    0    1    2    2
out[25,,drop=FALSE]
#      [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    1    2    0    0    1    2

现在进行快速健全性检查，以确保我们的 out 变量在每行中恰好有两个元素。

all(apply(out, 1, table) == 2)
# [1] TRUE