对于给定长度，总和为 1 的十进制数（百分之一）的所有可能排列
All possible permutations of decimal numbers (hundredths) that sum up to 1 for a given length
考虑向量 s 如下：
s=seq(0.01, 0.99, 0.01)
> s
[1] 0.01 0.02 0.03 0.04 0.05 0.06 0.07
0.08 0.09 .......... 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99
现在给定 s 和固定长度 m，我想要一个矩阵 所有可能的排列 长度 m这样矩阵的每一行总和为 1（不包括蛮力方法）。
例如，如果 m=4（即列数），所需的矩阵将如下所示：
0.01 0.01 0.01 0.97
0.02 0.01 0.01 0.96
0.03 0.01 0.01 0.95
0.04 0.01 0.01 0.94
0.05 0.01 0.01 0.93
0.06 0.01 0.01 0.92
.
.
.
0.53 0.12 0.30 0.05
.
.
.
0.96 0.02 0.01 0.01
0.97 0.01 0.01 0.01
.
.
.
0.01 0.97 0.01 0.01
.
.
.
以m=4为例，内存密集型方法是：
raw <- data.table::CJ(s,s,s,s)
result <- raw[rowSums(raw) == 1, ]
head(result)
V1 V2 V3 V4
1: 0.01 0.01 0.01 0.97
2: 0.01 0.01 0.02 0.96
3: 0.01 0.01 0.03 0.95
4: 0.01 0.01 0.04 0.94
5: 0.01 0.01 0.05 0.93
6: 0.01 0.01 0.06 0.92
以下是使用递归执行此操作的方法：
permsum <- function(s,m) if (m==1L) matrix(s) else do.call(rbind,lapply(seq_len(s-m+1L),function(x) unname(cbind(x,permsum(s-x,m-1L)))));
res <- permsum(100L,4L);
head(res);
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 97
## [2,] 1 1 2 96
## [3,] 1 1 3 95
## [4,] 1 1 4 94
## [5,] 1 1 5 93
## [6,] 1 1 6 92
tail(res);
## [,1] [,2] [,3] [,4]
## [156844,] 95 2 2 1
## [156845,] 95 3 1 1
## [156846,] 96 1 1 2
## [156847,] 96 1 2 1
## [156848,] 96 2 1 1
## [156849,] 97 1 1 1
您可以除以 100 得到分数，而不是整数：
head(res)/100;
## [,1] [,2] [,3] [,4]
## [1,] 0.01 0.01 0.01 0.97
## [2,] 0.01 0.01 0.02 0.96
## [3,] 0.01 0.01 0.03 0.95
## [4,] 0.01 0.01 0.04 0.94
## [5,] 0.01 0.01 0.05 0.93
## [6,] 0.01 0.01 0.06 0.92
说明
首先让我们定义输入：
s 这是输出矩阵中每一行应相加的目标值。m 这是要在输出矩阵中生成的列数。
与浮点运算相比，使用整数运算计算结果更有效、更可靠，因此我将我的解决方案设计为仅使用整数。因此 s 是表示目标整数和的标量整数。
现在让我们检查 seq_len() 为非基本情况生成的序列：
seq_len(s-m+1L)
这会生成一个从 1 到最大可能值的序列，该值可能是 s 总和的一部分，剩余 m 列。比如想一下s=100,m=4的情况：我们可以用的最大数是97，参与的和是97+1+1+1。剩余的每一列将最大可能值减 1，这就是为什么我们在计算序列长度时必须从 s 中减去 m。
生成序列的每个元素都应被视为求和中加数的一个可能 "selection"。
do.call(rbind,lapply(seq_len(s-m+1L),function(x) ...))
对于每个选择，我们必须递归。我们可以使用 lapply() 来做到这一点。
为了向前跳转，lambda 将对 permsum() 进行一次递归调用，然后 cbind() 使用当前选择的 return 值。这将产生一个矩阵，对于这个递归级别，它总是具有相同的宽度。因此，lapply() 调用将 return 一个矩阵列表，所有矩阵的宽度都相同。然后我们必须将它们行绑定在一起，这就是为什么我们必须在这里使用 do.call(rbind,...) 技巧。
unname(cbind(x,permsum(s-x,m-1L)))
lambda 的主体相当简单；我们 cbind() 当前选择 x 与递归调用的 return 值，完成此子矩阵的求和。不幸的是，我们必须调用 unname()，否则最终从 x 参数设置的每一列将具有列名 x.
这里最重要的细节是递归调用的参数选择。首先，因为 lambda 参数 x 刚刚在当前递归计算期间被选中，我们必须从 s 中减去它以获得新的求和目标，即将到来的递归调用将负责实现该目标。因此第一个参数变为 s-x。其次，因为x的选择占了一列，所以我们必须从m中减去1，这样递归调用在其输出矩阵中产生的列就会少一列。
if (m==1L) matrix(s) else ...
最后，让我们检查一下基本情况。在递归函数的每次评估中，我们必须检查 m 是否达到 1，在这种情况下，我们可以简单地 return 所需的总和 s 本身。
浮点数差异
我调查了我的结果和 psidom 的结果之间的差异。例如：
library(data.table);
bgoldst <- function(s,m) permsum(s,m)/s;
psidom <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[rowSums(raw)==1,]; };
## helper function to sort a matrix by columns
smp <- function(m) m[do.call(order,as.data.frame(m)),];
s <- 100L; m <- 3L; ss <- seq_len(s-1L)/s;
x <- smp(bgoldst(s,m));
y <- smp(unname(as.matrix(psidom(ss,m))));
nrow(x);
## [1] 4851
nrow(y);
## [1] 4809
所以我们的两个结果之间有 42 行的差异。我决定尝试通过以下代码行找出究竟省略了哪些排列。基本上，它比较两个矩阵的每个元素并将比较结果打印为逻辑矩阵。我们可以向下扫描回滚以找到第一个不同的行。以下是摘录的输出：
x==do.call(rbind,c(list(y),rep(list(NA),nrow(x)-nrow(y))));
## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE
## [5,] TRUE TRUE TRUE
##
## ... snip ...
##
## [24,] TRUE TRUE TRUE
## [25,] TRUE TRUE TRUE
## [26,] TRUE TRUE TRUE
## [27,] TRUE TRUE TRUE
## [28,] TRUE TRUE TRUE
## [29,] TRUE FALSE FALSE
## [30,] TRUE FALSE FALSE
## [31,] TRUE FALSE FALSE
## [32,] TRUE FALSE FALSE
## [33,] TRUE FALSE FALSE
##
## ... snip ...
第 29 行出现了第一个差异。这是每个置换矩阵中该行周围的 window：
win <- 27:31;
x[win,]; y[win,];
## [,1] [,2] [,3]
## [1,] 0.01 0.27 0.72
## [2,] 0.01 0.28 0.71
## [3,] 0.01 0.29 0.70 (missing from y)
## [4,] 0.01 0.30 0.69 (missing from y)
## [5,] 0.01 0.31 0.68
## [,1] [,2] [,3]
## [1,] 0.01 0.27 0.72
## [2,] 0.01 0.28 0.71
## [3,] 0.01 0.31 0.68
## [4,] 0.01 0.32 0.67
## [5,] 0.01 0.33 0.66
有趣的是，当您手动计算总和时，缺失的排列通常总和正好为 1。起初我以为是 data.table 的 CJ() 函数对浮点数做了一些奇怪的事情，但进一步的测试似乎表明它是 rowSums() 正在做的事情：
0.01+0.29+0.70==1;
## [1] TRUE
ss[1L]+ss[29L]+ss[70L]==1;
## [1] TRUE
rowSums(CJ(0.01,0.29,0.70))==1; ## looks like CJ()'s fault, but wait...
## [1] FALSE
cj <- CJ(0.01,0.29,0.70);
cj$V1+cj$V2+cj$V3==1; ## not CJ()'s fault
## [1] TRUE
rowSums(matrix(c(0.01,0.29,0.70),1L,byrow=T))==1; ## rowSums()'s fault
## [1] FALSE
我们可以通过在浮点比较中应用手动（并且有点武断）容差来解决这个 rowSums() 怪癖。为此，我们需要取绝对差，然后与公差进行小于比较：
abs(rowSums(CJ(0.01,0.29,0.70))-1)<1e-10;
## [1] TRUE
因此：
psidom2 <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[abs(rowSums(raw)-1)<1e-10,]; };
y <- smp(unname(as.matrix(psidom2(ss,m))));
nrow(y);
## [1] 4851
identical(x,y);
## [1] TRUE
组合
感谢 Joseph Wood 指出这实际上是 排列。我最初将我的函数命名为 combsum()，但我将其重命名为 permsum() 以反映这一启示。而且，正如 Joseph 所建议的那样，可以修改算法以生成组合，这可以按如下方式完成，现在名副其实 combsum():
combsum <- function(s,m,l=s) if (m==1L) matrix(s) else do.call(rbind,lapply(seq((s+m-1L)%/%m,min(l,s-m+1L)),function(x) unname(cbind(x,combsum(s-x,m-1L,x)))));
res <- combsum(100L,4L);
head(res);
## [,1] [,2] [,3] [,4]
## [1,] 25 25 25 25
## [2,] 26 25 25 24
## [3,] 26 26 24 24
## [4,] 26 26 25 23
## [5,] 26 26 26 22
## [6,] 27 25 24 24
tail(res);
## [,1] [,2] [,3] [,4]
## [7148,] 94 3 2 1
## [7149,] 94 4 1 1
## [7150,] 95 2 2 1
## [7151,] 95 3 1 1
## [7152,] 96 2 1 1
## [7153,] 97 1 1 1
这需要 3 处更改。
首先，我添加了一个新参数l，代表"limit"。基本上，为了保证每个递归生成唯一的组合，我强制每个选择必须 小于或等于 当前组合中的任何先前选择。这需要将当前上限作为参数l。在顶层调用 l 可以默认为 s，这对于 m>1 的情况来说实际上太高了，但这不是问题，因为它只是两个上层之一将在序列生成期间应用的限制。
第二个变化当然是在 lapply() lambda 中进行递归调用时将最新选择 x 作为参数传递给 l。
最后的改变是最棘手的。现在必须按如下方式计算选择序列：
seq((s+m-1L)%/%m,min(l,s-m+1L))
必须将下限从 permsum() 中使用的 1 提高到仍然允许降幅组合的最低可能选择。当然，最低可能的选择取决于尚未生产的色谱柱数量；列越多，我们为以后的选择留下的 "room" 就越多。公式是在 m 上进行 s 的整数除法，但我们还必须有效地 "round up"，这就是为什么我在进行除法之前添加 m-1L。我也考虑过做浮点除法然后调用 as.integer(ceiling(...))，但我认为全整数方法要好得多。
例如，考虑 s=10,m=3 的情况。要在剩余 3 列的情况下产生 10 的总和，我们不能选择小于 4 的值，因为如果不沿组合上升，我们将没有足够的数量来产生 10。在这种情况下，公式将 12 除以 3 得到 4。
可以使用 permsum() 中使用的相同公式计算上限，但我们还必须通过调用 min().[=83 来应用当前限制 l =]
我已经使用以下代码验证了我的新 combsum() 与 Joseph 的 IntegerPartitionsOfLength() 函数对于许多随机测试用例的行为相同：
## helper function to sort a matrix within each row and then by columns
smc <- function(m) smp(t(apply(m,1L,sort)));
## test loop
for (i in seq_len(1000L)) {
repeat {
s <- sample(1:100,1L);
m <- sample(2:5,1L);
if (s>=m) break;
};
x <- combsum(s,m);
y <- IntegerPartitionsOfLength(s,m);
cat(paste0(s,',',m,'\n'));
if (!identical(smc(x),smc(y))) stop('bad.');
};
基准测试
常用自包含测试代码：
library(microbenchmark);
library(data.table);
library(partitions);
library(gtools);
permsum <- function(s,m) if (m==1L) matrix(s) else do.call(rbind,lapply(seq_len(s-m+1L),function(x) unname(cbind(x,permsum(s-x,m-1L)))));
combsum <- function(s,m,l=s) if (m==1L) matrix(s) else do.call(rbind,lapply(seq((s+m-1L)%/%m,min(l,s-m+1L)),function(x) unname(cbind(x,combsum(s-x,m-1L,x)))));
IntegerPartitionsOfLength <- function(n, Lim, combsOnly = TRUE) { a <- 0L:n; k <- 2L; a[2L] <- n; MyParts <- vector("list", length=P(n)); count <- 0L; while (!(k==1L) && k <= Lim + 1L) { x <- a[k-1L]+1L; y <- a[k]-1L; k <- k-1L; while (x<=y && k <= Lim) {a[k] <- x; y <- y-x; k <- k+1L}; a[k] <- x+y; if (k==Lim) { count <- count+1L; MyParts[[count]] <- a[1L:k]; }; }; MyParts <- MyParts[1:count]; if (combsOnly) {do.call(rbind, MyParts)} else {MyParts}; };
GetDecimalReps <- function(s,m) { myPerms <- permutations(m,m); lim <- nrow(myPerms); intParts <- IntegerPartitionsOfLength(s,m,FALSE); do.call(rbind, lapply(intParts, function(x) { unique(t(sapply(1L:lim, function(y) x[myPerms[y, ]]))); })); };
smp <- function(m) m[do.call(order,as.data.frame(m)),];
smc <- function(m) smp(t(apply(m,1L,sort)));
bgoldst.perm <- function(s,m) permsum(s,m)/s;
psidom2 <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[abs(rowSums(raw)-1)<1e-10,]; };
joseph.perm <- function(s,m) GetDecimalReps(s,m)/s;
bgoldst.comb <- function(s,m) combsum(s,m)/s;
joseph.comb <- function(s,m) IntegerPartitionsOfLength(s,m)/s;
排列
## small scale
s <- 10L; m <- 3L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## [1] TRUE
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),psidom2(ss,m),joseph.perm(s,m));
## Unit: microseconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 347.254 389.5920 469.1011 420.383 478.7575 1869.697 100
## psidom2(ss, m) 702.206 830.5015 1007.5111 907.265 1038.3405 2618.089 100
## joseph.perm(s, m) 1225.225 1392.8640 1722.0070 1506.833 1860.0745 4411.234 100
## large scale
s <- 100L; m <- 4L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## [1] TRUE
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),psidom2(ss,m),joseph.perm(s,m),times=5L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 1.286856 1.304177 1.426376 1.374411 1.399850 1.766585 5
## psidom2(ss, m) 6.673545 7.046951 7.416161 7.115375 7.629177 8.615757 5
## joseph.perm(s, m) 5.299452 10.499891 13.769363 12.680607 15.107748 25.259117 5
## very large scale
s <- 100L; m <- 5L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## Error: cannot allocate vector of size 70.9 Gb
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),joseph.perm(s,m),times=1L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 28.58359 28.58359 28.58359 28.58359 28.58359 28.58359 1
## joseph.perm(s, m) 50.51965 50.51965 50.51965 50.51965 50.51965 50.51965 1
组合
## small-scale
s <- 10L; m <- 3L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m));
## Unit: microseconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 161.225 179.6145 205.0898 187.3120 199.5005 1310.328 100
## joseph.comb(s, m) 172.344 191.8025 204.5681 197.7895 205.2735 437.489 100
## large-scale
s <- 100L; m <- 4L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m),times=5L);
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 409.0708 485.9739 556.4792 591.4774 627.419 668.4548 5
## joseph.comb(s, m) 2164.2134 3315.0138 3317.9725 3540.6240 3713.732 3856.2793 5
## very large scale
s <- 100L; m <- 6L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m),times=1L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 2.498588 2.498588 2.498588 2.498588 2.498588 2.498588 1
## joseph.comb(s, m) 12.344261 12.344261 12.344261 12.344261 12.344261 12.344261 1
这是一个 return 纯 combinations 的算法（顺序无关紧要）。它基于 Jerome Kelleher (link) 构建的整数分区算法。
library(partitions)
IntegerPartitionsOfLength <- function(n, Lim, combsOnly = TRUE) {
a <- 0L:n
k <- 2L
a[2L] <- n
MyParts <- vector("list", length=P(n))
count <- 0L
while (!(k==1L) && k <= Lim + 1L) {
x <- a[k-1L]+1L
y <- a[k]-1L
k <- k-1L
while (x<=y && k <= Lim) {a[k] <- x; y <- y-x; k <- k+1L}
a[k] <- x+y
if (k==Lim) {
count <- count+1L
MyParts[[count]] <- a[1L:k]
}
}
MyParts <- MyParts[1:count]
if (combsOnly) {do.call(rbind, MyParts)} else {MyParts}
}
system.time(res <- combsum(100L,5L))
user system elapsed
0.75 0.00 0.77
system.time(a <- IntegerPartitionsOfLength(100, 5))
user system elapsed
1.36 0.37 1.76
identical(smc(a),smc(res))
[1] TRUE
head(a)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 1 1 1 96
[2,] 1 1 1 2 95
[3,] 1 1 1 3 94
[4,] 1 1 1 4 93
[5,] 1 1 1 5 92
[6,] 1 1 1 6 91
一个非常大的例子（N.B。使用@bgoldst 创建的 smc 函数）：
system.time(a <- IntegerPartitionsOfLength(100L,6L))
user system elapsed
4.57 0.36 4.93
system.time(res <- combsum(100L,6L))
user system elapsed
3.69 0.00 3.71
identical(smc(a),smc(res))
[1] TRUE
## this would take a very long time with GetDecimalReps below
注意：IntegerPartitionsOfLength 只是 return 一组特定数字的 combinations 而不是一组数字的 permutations（顺序很重要）。例如。对于集合s = (1, 1, 3)，s的组合正好是s，而s的排列是：(1, 1, 3), (1, 3, 1), (3, 1, 1).
如果你想要像 OP 要求的那样的答案，你将不得不做这样的事情（这绝不是最好的方法，而且效率不如上面的@bgoldst permsum ):
library(gtools)
GetDecimalReps <- function(n) {
myPerms <- permutations(n,n); lim <- nrow(myPerms)
intParts <- IntegerPartitionsOfLength(100,n,FALSE)
do.call(rbind, lapply(intParts, function(x) {
unique(t(sapply(1L:lim, function(y) x[myPerms[y, ]])))
}))
}
system.time(a <- GetDecimalReps(4L))
user system elapsed
2.85 0.42 3.28
system.time(res <- combsum(100L,4L))
user system elapsed
1.35 0.00 1.34
head(a/100)
[,1] [,2] [,3] [,4]
[1,] 0.01 0.01 0.01 0.97
[2,] 0.01 0.01 0.97 0.01
[3,] 0.01 0.97 0.01 0.01
[4,] 0.97 0.01 0.01 0.01
[5,] 0.01 0.01 0.02 0.96
[6,] 0.01 0.01 0.96 0.02
tail(a/100)
[,1] [,2] [,3] [,4]
[156844,] 0.25 0.26 0.24 0.25
[156845,] 0.25 0.26 0.25 0.24
[156846,] 0.26 0.24 0.25 0.25
[156847,] 0.26 0.25 0.24 0.25
[156848,] 0.26 0.25 0.25 0.24
[156849,] 0.25 0.25 0.25 0.25
identical(smp(a),smp(res)) ## using the smp function created by @bgoldst
[1] TRUE
@bgoldst 上面的算法对于两种 return 类型（即 combinations/permutations）都是优越的。另请参阅上面@bgoldst 的优秀基准测试。作为结束语，您可以轻松修改 IntegerPartionsOfLength 以获得 1:100 的所有组合，只需将 k==Lim 更改为 k <= Lim 并设置 combsOnly = FALSE 以便 return 列表。干杯！
考虑向量 s 如下：
s=seq(0.01, 0.99, 0.01)
> s
[1] 0.01 0.02 0.03 0.04 0.05 0.06 0.07
0.08 0.09 .......... 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99
现在给定 s 和固定长度 m，我想要一个矩阵 所有可能的排列 长度 m这样矩阵的每一行总和为 1（不包括蛮力方法）。
例如，如果 m=4（即列数），所需的矩阵将如下所示：
0.01 0.01 0.01 0.97
0.02 0.01 0.01 0.96
0.03 0.01 0.01 0.95
0.04 0.01 0.01 0.94
0.05 0.01 0.01 0.93
0.06 0.01 0.01 0.92
.
.
.
0.53 0.12 0.30 0.05
.
.
.
0.96 0.02 0.01 0.01
0.97 0.01 0.01 0.01
.
.
.
0.01 0.97 0.01 0.01
.
.
.
以m=4为例，内存密集型方法是：
raw <- data.table::CJ(s,s,s,s)
result <- raw[rowSums(raw) == 1, ]
head(result)
V1 V2 V3 V4
1: 0.01 0.01 0.01 0.97
2: 0.01 0.01 0.02 0.96
3: 0.01 0.01 0.03 0.95
4: 0.01 0.01 0.04 0.94
5: 0.01 0.01 0.05 0.93
6: 0.01 0.01 0.06 0.92
以下是使用递归执行此操作的方法：
permsum <- function(s,m) if (m==1L) matrix(s) else do.call(rbind,lapply(seq_len(s-m+1L),function(x) unname(cbind(x,permsum(s-x,m-1L)))));
res <- permsum(100L,4L);
head(res);
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 97
## [2,] 1 1 2 96
## [3,] 1 1 3 95
## [4,] 1 1 4 94
## [5,] 1 1 5 93
## [6,] 1 1 6 92
tail(res);
## [,1] [,2] [,3] [,4]
## [156844,] 95 2 2 1
## [156845,] 95 3 1 1
## [156846,] 96 1 1 2
## [156847,] 96 1 2 1
## [156848,] 96 2 1 1
## [156849,] 97 1 1 1
您可以除以 100 得到分数，而不是整数：
head(res)/100;
## [,1] [,2] [,3] [,4]
## [1,] 0.01 0.01 0.01 0.97
## [2,] 0.01 0.01 0.02 0.96
## [3,] 0.01 0.01 0.03 0.95
## [4,] 0.01 0.01 0.04 0.94
## [5,] 0.01 0.01 0.05 0.93
## [6,] 0.01 0.01 0.06 0.92
说明
首先让我们定义输入：
s这是输出矩阵中每一行应相加的目标值。m这是要在输出矩阵中生成的列数。
与浮点运算相比，使用整数运算计算结果更有效、更可靠，因此我将我的解决方案设计为仅使用整数。因此 s 是表示目标整数和的标量整数。
现在让我们检查 seq_len() 为非基本情况生成的序列：
seq_len(s-m+1L)
这会生成一个从 1 到最大可能值的序列，该值可能是 s 总和的一部分，剩余 m 列。比如想一下s=100,m=4的情况：我们可以用的最大数是97，参与的和是97+1+1+1。剩余的每一列将最大可能值减 1，这就是为什么我们在计算序列长度时必须从 s 中减去 m。
生成序列的每个元素都应被视为求和中加数的一个可能 "selection"。
do.call(rbind,lapply(seq_len(s-m+1L),function(x) ...))
对于每个选择，我们必须递归。我们可以使用 lapply() 来做到这一点。
为了向前跳转，lambda 将对 permsum() 进行一次递归调用，然后 cbind() 使用当前选择的 return 值。这将产生一个矩阵，对于这个递归级别，它总是具有相同的宽度。因此，lapply() 调用将 return 一个矩阵列表，所有矩阵的宽度都相同。然后我们必须将它们行绑定在一起，这就是为什么我们必须在这里使用 do.call(rbind,...) 技巧。
unname(cbind(x,permsum(s-x,m-1L)))
lambda 的主体相当简单；我们 cbind() 当前选择 x 与递归调用的 return 值，完成此子矩阵的求和。不幸的是，我们必须调用 unname()，否则最终从 x 参数设置的每一列将具有列名 x.
这里最重要的细节是递归调用的参数选择。首先，因为 lambda 参数 x 刚刚在当前递归计算期间被选中，我们必须从 s 中减去它以获得新的求和目标，即将到来的递归调用将负责实现该目标。因此第一个参数变为 s-x。其次，因为x的选择占了一列，所以我们必须从m中减去1，这样递归调用在其输出矩阵中产生的列就会少一列。
if (m==1L) matrix(s) else ...
最后，让我们检查一下基本情况。在递归函数的每次评估中，我们必须检查 m 是否达到 1，在这种情况下，我们可以简单地 return 所需的总和 s 本身。
浮点数差异
我调查了我的结果和 psidom 的结果之间的差异。例如：
library(data.table);
bgoldst <- function(s,m) permsum(s,m)/s;
psidom <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[rowSums(raw)==1,]; };
## helper function to sort a matrix by columns
smp <- function(m) m[do.call(order,as.data.frame(m)),];
s <- 100L; m <- 3L; ss <- seq_len(s-1L)/s;
x <- smp(bgoldst(s,m));
y <- smp(unname(as.matrix(psidom(ss,m))));
nrow(x);
## [1] 4851
nrow(y);
## [1] 4809
所以我们的两个结果之间有 42 行的差异。我决定尝试通过以下代码行找出究竟省略了哪些排列。基本上，它比较两个矩阵的每个元素并将比较结果打印为逻辑矩阵。我们可以向下扫描回滚以找到第一个不同的行。以下是摘录的输出：
x==do.call(rbind,c(list(y),rep(list(NA),nrow(x)-nrow(y))));
## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE
## [5,] TRUE TRUE TRUE
##
## ... snip ...
##
## [24,] TRUE TRUE TRUE
## [25,] TRUE TRUE TRUE
## [26,] TRUE TRUE TRUE
## [27,] TRUE TRUE TRUE
## [28,] TRUE TRUE TRUE
## [29,] TRUE FALSE FALSE
## [30,] TRUE FALSE FALSE
## [31,] TRUE FALSE FALSE
## [32,] TRUE FALSE FALSE
## [33,] TRUE FALSE FALSE
##
## ... snip ...
第 29 行出现了第一个差异。这是每个置换矩阵中该行周围的 window：
win <- 27:31;
x[win,]; y[win,];
## [,1] [,2] [,3]
## [1,] 0.01 0.27 0.72
## [2,] 0.01 0.28 0.71
## [3,] 0.01 0.29 0.70 (missing from y)
## [4,] 0.01 0.30 0.69 (missing from y)
## [5,] 0.01 0.31 0.68
## [,1] [,2] [,3]
## [1,] 0.01 0.27 0.72
## [2,] 0.01 0.28 0.71
## [3,] 0.01 0.31 0.68
## [4,] 0.01 0.32 0.67
## [5,] 0.01 0.33 0.66
有趣的是，当您手动计算总和时，缺失的排列通常总和正好为 1。起初我以为是 data.table 的 CJ() 函数对浮点数做了一些奇怪的事情，但进一步的测试似乎表明它是 rowSums() 正在做的事情：
0.01+0.29+0.70==1;
## [1] TRUE
ss[1L]+ss[29L]+ss[70L]==1;
## [1] TRUE
rowSums(CJ(0.01,0.29,0.70))==1; ## looks like CJ()'s fault, but wait...
## [1] FALSE
cj <- CJ(0.01,0.29,0.70);
cj$V1+cj$V2+cj$V3==1; ## not CJ()'s fault
## [1] TRUE
rowSums(matrix(c(0.01,0.29,0.70),1L,byrow=T))==1; ## rowSums()'s fault
## [1] FALSE
我们可以通过在浮点比较中应用手动（并且有点武断）容差来解决这个 rowSums() 怪癖。为此，我们需要取绝对差，然后与公差进行小于比较：
abs(rowSums(CJ(0.01,0.29,0.70))-1)<1e-10;
## [1] TRUE
因此：
psidom2 <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[abs(rowSums(raw)-1)<1e-10,]; };
y <- smp(unname(as.matrix(psidom2(ss,m))));
nrow(y);
## [1] 4851
identical(x,y);
## [1] TRUE
组合
感谢 Joseph Wood 指出这实际上是 排列。我最初将我的函数命名为 combsum()，但我将其重命名为 permsum() 以反映这一启示。而且，正如 Joseph 所建议的那样，可以修改算法以生成组合，这可以按如下方式完成，现在名副其实 combsum():
combsum <- function(s,m,l=s) if (m==1L) matrix(s) else do.call(rbind,lapply(seq((s+m-1L)%/%m,min(l,s-m+1L)),function(x) unname(cbind(x,combsum(s-x,m-1L,x)))));
res <- combsum(100L,4L);
head(res);
## [,1] [,2] [,3] [,4]
## [1,] 25 25 25 25
## [2,] 26 25 25 24
## [3,] 26 26 24 24
## [4,] 26 26 25 23
## [5,] 26 26 26 22
## [6,] 27 25 24 24
tail(res);
## [,1] [,2] [,3] [,4]
## [7148,] 94 3 2 1
## [7149,] 94 4 1 1
## [7150,] 95 2 2 1
## [7151,] 95 3 1 1
## [7152,] 96 2 1 1
## [7153,] 97 1 1 1
这需要 3 处更改。
首先，我添加了一个新参数l，代表"limit"。基本上，为了保证每个递归生成唯一的组合，我强制每个选择必须 小于或等于 当前组合中的任何先前选择。这需要将当前上限作为参数l。在顶层调用 l 可以默认为 s，这对于 m>1 的情况来说实际上太高了，但这不是问题，因为它只是两个上层之一将在序列生成期间应用的限制。
第二个变化当然是在 lapply() lambda 中进行递归调用时将最新选择 x 作为参数传递给 l。
最后的改变是最棘手的。现在必须按如下方式计算选择序列：
seq((s+m-1L)%/%m,min(l,s-m+1L))
必须将下限从 permsum() 中使用的 1 提高到仍然允许降幅组合的最低可能选择。当然，最低可能的选择取决于尚未生产的色谱柱数量；列越多，我们为以后的选择留下的 "room" 就越多。公式是在 m 上进行 s 的整数除法，但我们还必须有效地 "round up"，这就是为什么我在进行除法之前添加 m-1L。我也考虑过做浮点除法然后调用 as.integer(ceiling(...))，但我认为全整数方法要好得多。
例如，考虑 s=10,m=3 的情况。要在剩余 3 列的情况下产生 10 的总和，我们不能选择小于 4 的值，因为如果不沿组合上升，我们将没有足够的数量来产生 10。在这种情况下，公式将 12 除以 3 得到 4。
可以使用 permsum() 中使用的相同公式计算上限，但我们还必须通过调用 min().[=83 来应用当前限制 l =]
我已经使用以下代码验证了我的新 combsum() 与 Joseph 的 IntegerPartitionsOfLength() 函数对于许多随机测试用例的行为相同：
## helper function to sort a matrix within each row and then by columns
smc <- function(m) smp(t(apply(m,1L,sort)));
## test loop
for (i in seq_len(1000L)) {
repeat {
s <- sample(1:100,1L);
m <- sample(2:5,1L);
if (s>=m) break;
};
x <- combsum(s,m);
y <- IntegerPartitionsOfLength(s,m);
cat(paste0(s,',',m,'\n'));
if (!identical(smc(x),smc(y))) stop('bad.');
};
基准测试
常用自包含测试代码：
library(microbenchmark);
library(data.table);
library(partitions);
library(gtools);
permsum <- function(s,m) if (m==1L) matrix(s) else do.call(rbind,lapply(seq_len(s-m+1L),function(x) unname(cbind(x,permsum(s-x,m-1L)))));
combsum <- function(s,m,l=s) if (m==1L) matrix(s) else do.call(rbind,lapply(seq((s+m-1L)%/%m,min(l,s-m+1L)),function(x) unname(cbind(x,combsum(s-x,m-1L,x)))));
IntegerPartitionsOfLength <- function(n, Lim, combsOnly = TRUE) { a <- 0L:n; k <- 2L; a[2L] <- n; MyParts <- vector("list", length=P(n)); count <- 0L; while (!(k==1L) && k <= Lim + 1L) { x <- a[k-1L]+1L; y <- a[k]-1L; k <- k-1L; while (x<=y && k <= Lim) {a[k] <- x; y <- y-x; k <- k+1L}; a[k] <- x+y; if (k==Lim) { count <- count+1L; MyParts[[count]] <- a[1L:k]; }; }; MyParts <- MyParts[1:count]; if (combsOnly) {do.call(rbind, MyParts)} else {MyParts}; };
GetDecimalReps <- function(s,m) { myPerms <- permutations(m,m); lim <- nrow(myPerms); intParts <- IntegerPartitionsOfLength(s,m,FALSE); do.call(rbind, lapply(intParts, function(x) { unique(t(sapply(1L:lim, function(y) x[myPerms[y, ]]))); })); };
smp <- function(m) m[do.call(order,as.data.frame(m)),];
smc <- function(m) smp(t(apply(m,1L,sort)));
bgoldst.perm <- function(s,m) permsum(s,m)/s;
psidom2 <- function(ss,m) { raw <- do.call(data.table::CJ,rep(list(ss),m)); raw[abs(rowSums(raw)-1)<1e-10,]; };
joseph.perm <- function(s,m) GetDecimalReps(s,m)/s;
bgoldst.comb <- function(s,m) combsum(s,m)/s;
joseph.comb <- function(s,m) IntegerPartitionsOfLength(s,m)/s;
排列
## small scale
s <- 10L; m <- 3L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## [1] TRUE
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),psidom2(ss,m),joseph.perm(s,m));
## Unit: microseconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 347.254 389.5920 469.1011 420.383 478.7575 1869.697 100
## psidom2(ss, m) 702.206 830.5015 1007.5111 907.265 1038.3405 2618.089 100
## joseph.perm(s, m) 1225.225 1392.8640 1722.0070 1506.833 1860.0745 4411.234 100
## large scale
s <- 100L; m <- 4L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## [1] TRUE
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),psidom2(ss,m),joseph.perm(s,m),times=5L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 1.286856 1.304177 1.426376 1.374411 1.399850 1.766585 5
## psidom2(ss, m) 6.673545 7.046951 7.416161 7.115375 7.629177 8.615757 5
## joseph.perm(s, m) 5.299452 10.499891 13.769363 12.680607 15.107748 25.259117 5
## very large scale
s <- 100L; m <- 5L; ss <- seq_len(s-1L)/s;
ex <- smp(bgoldst.perm(s,m));
identical(ex,smp(unname(as.matrix(psidom2(ss,m)))));
## Error: cannot allocate vector of size 70.9 Gb
identical(ex,smp(joseph.perm(s,m)));
## [1] TRUE
microbenchmark(bgoldst.perm(s,m),joseph.perm(s,m),times=1L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.perm(s, m) 28.58359 28.58359 28.58359 28.58359 28.58359 28.58359 1
## joseph.perm(s, m) 50.51965 50.51965 50.51965 50.51965 50.51965 50.51965 1
组合
## small-scale
s <- 10L; m <- 3L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m));
## Unit: microseconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 161.225 179.6145 205.0898 187.3120 199.5005 1310.328 100
## joseph.comb(s, m) 172.344 191.8025 204.5681 197.7895 205.2735 437.489 100
## large-scale
s <- 100L; m <- 4L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m),times=5L);
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 409.0708 485.9739 556.4792 591.4774 627.419 668.4548 5
## joseph.comb(s, m) 2164.2134 3315.0138 3317.9725 3540.6240 3713.732 3856.2793 5
## very large scale
s <- 100L; m <- 6L;
ex <- smc(bgoldst.comb(s,m));
identical(ex,smc(joseph.comb(s,m)));
## [1] TRUE
microbenchmark(bgoldst.comb(s,m),joseph.comb(s,m),times=1L);
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst.comb(s, m) 2.498588 2.498588 2.498588 2.498588 2.498588 2.498588 1
## joseph.comb(s, m) 12.344261 12.344261 12.344261 12.344261 12.344261 12.344261 1
这是一个 return 纯 combinations 的算法（顺序无关紧要）。它基于 Jerome Kelleher (link) 构建的整数分区算法。
library(partitions)
IntegerPartitionsOfLength <- function(n, Lim, combsOnly = TRUE) {
a <- 0L:n
k <- 2L
a[2L] <- n
MyParts <- vector("list", length=P(n))
count <- 0L
while (!(k==1L) && k <= Lim + 1L) {
x <- a[k-1L]+1L
y <- a[k]-1L
k <- k-1L
while (x<=y && k <= Lim) {a[k] <- x; y <- y-x; k <- k+1L}
a[k] <- x+y
if (k==Lim) {
count <- count+1L
MyParts[[count]] <- a[1L:k]
}
}
MyParts <- MyParts[1:count]
if (combsOnly) {do.call(rbind, MyParts)} else {MyParts}
}
system.time(res <- combsum(100L,5L))
user system elapsed
0.75 0.00 0.77
system.time(a <- IntegerPartitionsOfLength(100, 5))
user system elapsed
1.36 0.37 1.76
identical(smc(a),smc(res))
[1] TRUE
head(a)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 1 1 1 96
[2,] 1 1 1 2 95
[3,] 1 1 1 3 94
[4,] 1 1 1 4 93
[5,] 1 1 1 5 92
[6,] 1 1 1 6 91
一个非常大的例子（N.B。使用@bgoldst 创建的 smc 函数）：
system.time(a <- IntegerPartitionsOfLength(100L,6L))
user system elapsed
4.57 0.36 4.93
system.time(res <- combsum(100L,6L))
user system elapsed
3.69 0.00 3.71
identical(smc(a),smc(res))
[1] TRUE
## this would take a very long time with GetDecimalReps below
注意：IntegerPartitionsOfLength 只是 return 一组特定数字的 combinations 而不是一组数字的 permutations（顺序很重要）。例如。对于集合s = (1, 1, 3)，s的组合正好是s，而s的排列是：(1, 1, 3), (1, 3, 1), (3, 1, 1).
如果你想要像 OP 要求的那样的答案，你将不得不做这样的事情（这绝不是最好的方法，而且效率不如上面的@bgoldst permsum ):
library(gtools)
GetDecimalReps <- function(n) {
myPerms <- permutations(n,n); lim <- nrow(myPerms)
intParts <- IntegerPartitionsOfLength(100,n,FALSE)
do.call(rbind, lapply(intParts, function(x) {
unique(t(sapply(1L:lim, function(y) x[myPerms[y, ]])))
}))
}
system.time(a <- GetDecimalReps(4L))
user system elapsed
2.85 0.42 3.28
system.time(res <- combsum(100L,4L))
user system elapsed
1.35 0.00 1.34
head(a/100)
[,1] [,2] [,3] [,4]
[1,] 0.01 0.01 0.01 0.97
[2,] 0.01 0.01 0.97 0.01
[3,] 0.01 0.97 0.01 0.01
[4,] 0.97 0.01 0.01 0.01
[5,] 0.01 0.01 0.02 0.96
[6,] 0.01 0.01 0.96 0.02
tail(a/100)
[,1] [,2] [,3] [,4]
[156844,] 0.25 0.26 0.24 0.25
[156845,] 0.25 0.26 0.25 0.24
[156846,] 0.26 0.24 0.25 0.25
[156847,] 0.26 0.25 0.24 0.25
[156848,] 0.26 0.25 0.25 0.24
[156849,] 0.25 0.25 0.25 0.25
identical(smp(a),smp(res)) ## using the smp function created by @bgoldst
[1] TRUE
@bgoldst 上面的算法对于两种 return 类型（即 combinations/permutations）都是优越的。另请参阅上面@bgoldst 的优秀基准测试。作为结束语，您可以轻松修改 IntegerPartionsOfLength 以获得 1:100 的所有组合，只需将 k==Lim 更改为 k <= Lim 并设置 combsOnly = FALSE 以便 return 列表。干杯！