989*_*989 4 combinations r matrix
考虑向量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
.
.
.
以下是使用递归的方法:
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减少了可能的最高值,这就是为什么我们必须减去m从s计算序列长度时.
生成序列的每个元素应被视为求和中加数的一个可能的"选择".
do.call(rbind,lapply(seq_len(s-m+1L),function(x) ...))
对于每个选择,我们必须递归.我们可以用它lapply()来做到这一点.
为了向前跳,lambda将对当前选择进行单个递归调用permsum(),然后cbind()返回值.这将产生一个矩阵,对于这种递归级别总是具有相同的宽度.因此,lapply()调用将返回矩阵列表,所有矩阵都具有相同的宽度.然后我们必须将它们绑定在一起,这就是为什么我们必须在do.call(rbind,...)这里使用这个技巧.
unname(cbind(x,permsum(s-x,m-1L)))
lambda的身体相当简单; 我们使用递归调用的返回值进行cbind()当前选择x,完成此子矩阵的求和.不幸的是我们必须调用unname(),否则每个最终从x参数设置的列都将具有列名x.
这里最重要的细节是递归调用的参数选择.首先,因为拉姆达的说法x刚刚在当前递归评估过程中选择了,我们必须减去它s得到一个新的概括目标,这即将到来的递归调用将负责实现.因此第一个论点变成了s-x.其次,因为选择x占用一列,我们必须从中减去1 m,这样递归调用将在其输出矩阵中产生少一列.
if (m==1L) matrix(s) else ...
最后,让我们来看一下基本情况.在每次递归函数的评估中,我们必须检查是否m已达到1,在这种情况下,我们可以简单地返回所需的和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行我们有第一个差异.这是每个排列矩阵中该行周围的窗口:
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()以反映这个启示.并且,正如约瑟夫所建议的,可以修改算法以产生组合,这可以按照以下方式完成,现在符合名称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,代表"限制".基本上,为了保证每个递归生成唯一组合,我强制每个选择必须小于或等于当前组合中的任何先前选择.这需要将当前上限作为参数l.在上层呼叫l可以直接被默认为s,这实际上是太高反正情况m>1,但是这不是一个问题,因为这只是将序列生成过程中应用了两个上限之一.
第二个变化当然是将最新的选择x作为参数传递l给lapply()lambda中的递归调用.
最后的改变是最棘手的.现在必须按如下方式计算选择顺序:
seq((s+m-1L)%/%m,min(l,s-m+1L))
必须将下限从使用的1提高permsum()到最低可能的选择,这仍然允许下降幅度组合.当然,最低的选择取决于尚未生产的柱数; 列越多,我们必须留下更多的"空间"以供将来选择.其计算公式为取的整数除法s的m,但我们也必须有效地"全面上涨",这就是为什么我添加m-1L之前采取分工.我还考虑过进行浮点除法然后调用as.integer(ceiling(...)),但我认为全整数方法要好得多.
例如,考虑一下s=10,m=3.为了产生10个剩余3列的总和,我们不能选择小于4,因为那样我们就没有足够的数量来产生10而没有沿着组合上升.在这种情况下,公式将12除以3得到4.
上限可以使用相同的公式计算permsum(),除了我们还必须l使用调用来应用当前限制min().
我已经验证了我的新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