解决方案
这并不容易,但我想出了一种方法来完成这项任务,使用单个矢量化调用rep()
,加上一些脚手架代码:
XR <- 3;
YC <- 4;
X <- matrix(c(1:XR%%(YC+1),seq(10,by=10,length.out=XR),seq(0.1,by=0.1,length.out=XR)),XR,dimnames=list(NULL,c('rep','val','fill')));
X;
## rep val fill
## [1,] 1 10 0.1
## [2,] 2 20 0.2
## [3,] 3 30 0.3
Y <- matrix(rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep']))),XR,byrow=T);
Y;
## [,1] [,2] [,3] [,4]
## [1,] 10 0.1 0.1 0.1
## [2,] 20 20.0 0.2 0.2
## [3,] 30 30.0 30.0 0.3
(次要点:我选择将列名分配rep val fill
给X
,而不是a b c
问题中指定的,并且我在索引时在我的解决方案中使用了这些列名X
(而不是使用数字索引),因为我通常更喜欢最大化人类可读性在可能的情况下,但这个细节对于解决方案的正确性和性能可以忽略不计。)
表现
这实际上比@josilber 的解决方案具有显着的性能优势,因为他使用apply()
which 在矩阵的行上进行内部循环(在 R 语言中传统上称为“隐藏循环”),而我的解决方案的核心是对rep()
. 我这样说不是为了敲@josilber 的解决方案,这是一个很好的解决方案(我什至给了他一个赞成票!);这不是解决这个问题的最佳解决方案。
这是使用您在问题中指出的大量参数的性能优势演示:
XR <- 1e6;
YC <- 480;
X <- matrix(c(1:XR%%(YC+1),seq(10,by=10,length.out=XR),seq(0.1,by=0.1,length.out=XR)),XR,dimnames=list(NULL,c('rep','val','fill')));
X;
## rep val fill
## [1,] 1 10 0.1
## [2,] 2 20 0.2
## [3,] 3 30 0.3
## [4,] 4 40 0.4
## [5,] 5 50 0.5
## [6,] 6 60 0.6
## [7,] 7 70 0.7
## [8,] 8 80 0.8
## [9,] 9 90 0.9
## [10,] 10 100 1.0
## [11,] 11 110 1.1
## [12,] 12 120 1.2
## [13,] 13 130 1.3
##
## ... (snip) ...
##
## [477,] 477 4770 47.7
## [478,] 478 4780 47.8
## [479,] 479 4790 47.9
## [480,] 480 4800 48.0
## [481,] 0 4810 48.1
## [482,] 1 4820 48.2
## [483,] 2 4830 48.3
## [484,] 3 4840 48.4
## [485,] 4 4850 48.5
## [486,] 5 4860 48.6
## [487,] 6 4870 48.7
## [488,] 7 4880 48.8
## [489,] 8 4890 48.9
## [490,] 9 4900 49.0
## [491,] 10 4910 49.1
## [492,] 11 4920 49.2
##
## ... (snip) ...
##
## [999986,] 468 9999860 99998.6
## [999987,] 469 9999870 99998.7
## [999988,] 470 9999880 99998.8
## [999989,] 471 9999890 99998.9
## [999990,] 472 9999900 99999.0
## [999991,] 473 9999910 99999.1
## [999992,] 474 9999920 99999.2
## [999993,] 475 9999930 99999.3
## [999994,] 476 9999940 99999.4
## [999995,] 477 9999950 99999.5
## [999996,] 478 9999960 99999.6
## [999997,] 479 9999970 99999.7
## [999998,] 480 9999980 99999.8
## [999999,] 0 9999990 99999.9
## [1e+06,] 1 10000000 100000.0
josilber <- function() t(apply(X,1,function(x) rep(x[2:3],c(x[1],YC-x[1]))));
bgoldst <- function() matrix(rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep']))),XR,byrow=T);
system.time({ josilber(); });
## user system elapsed
## 65.719 3.828 71.623
system.time({ josilber(); });
## user system elapsed
## 60.375 2.609 66.724
system.time({ bgoldst(); });
## user system elapsed
## 5.422 0.593 6.033
system.time({ bgoldst(); });
## user system elapsed
## 5.203 0.797 6.002
只是为了证明@josilber 和我得到了完全相同的结果,即使对于这么大的输入:
identical(bgoldst(),josilber());
## [1] TRUE
解释
现在我将尝试解释解决方案的工作原理。为了解释,我将使用以下输入:
XR <- 6;
YC <- 4;
X <- matrix(c(1:XR%%(YC+1),seq(10,by=10,length.out=XR),seq(0.1,by=0.1,length.out=XR)),XR,dimnames=list(NULL,c('rep','val','fill')));
X;
## rep val fill
## [1,] 1 10 0.1
## [2,] 2 20 0.2
## [3,] 3 30 0.3
## [4,] 4 40 0.4
## [5,] 0 50 0.5
## [6,] 1 60 0.6
解决方案是:
Y <- matrix(rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep']))),XR,byrow=T);
Y;
## [,1] [,2] [,3] [,4]
## [1,] 10.0 0.1 0.1 0.1
## [2,] 20.0 20.0 0.2 0.2
## [3,] 30.0 30.0 30.0 0.3
## [4,] 40.0 40.0 40.0 40.0
## [5,] 0.5 0.5 0.5 0.5
## [6,] 60.0 0.6 0.6 0.6
在高层次上,解决方案是围绕形成一个组合val
和fill
向量的单个向量,然后以某种方式重复该组合向量,然后从结果中构建一个新矩阵。
重复步骤可以使用一次调用来完成,rep()
因为它支持向量化重复计数。换句话说,对于给定的向量输入x
,它可以采用一个向量输入,times
该向量输入指定重复 的每个元素的次数x
。因此,挑战就变成了构建适当的x
论点times
。
因此,解决方案首先提取 的val
和fill
列X
:
X[,c('val','fill')];
## val fill
## [1,] 10 0.1
## [2,] 20 0.2
## [3,] 30 0.3
## [4,] 40 0.4
## [5,] 50 0.5
## [6,] 60 0.6
如您所见,由于我们已经索引了两列,所以我们仍然有一个矩阵,即使我们没有指定drop=F
索引操作(参见R:提取或替换对象的部分)。正如将要看到的,这很方便。
在 R 中,矩阵的“矩阵角色”下面实际上只是一个普通的旧原子向量,矩阵的“向量角色”可以用于向量化操作。这就是我们如何将val
andfill
数据传递给rep()
并适当地重复这些元素的方式。
但是,在执行此操作时,重要的是要准确了解如何将矩阵视为向量。答案是向量是由跨行的元素形成的,然后才跨列的。(对于更高维的数组,随后的维数会跟随。IOW,向量的顺序是跨行,然后是列,然后是 z 切片等)
如果您仔细查看上面的矩阵,您会发现它不能用作我们的x
参数rep()
,因为val
首先会跟随 s,然后是fill
s。实际上,我们可以相当容易地构造一个times
参数来重复每个元素正确的次数,但结果向量将完全无序,并且无法将其重塑为所需的矩阵Y
。
实际上,在继续解释之前,我为什么不快速演示一下:
rep(X[,c('val','fill')],times=c(X[,'rep'],YC-X[,'rep']))
## [1] 10.0 20.0 20.0 30.0 30.0 30.0 40.0 40.0 40.0 40.0 60.0 0.1 0.1 0.1 0.2 0.2 0.3 0.5 0.5 0.5 0.5 0.6 0.6 0.6
尽管上面的向量在所有正确的重复中都有所有正确的元素,但顺序是这样的,它不能形成所需的输出矩阵Y
。
所以,我们可以通过首先转置提取来解决这个问题:
t(X[,c('val','fill')]);
## [,1] [,2] [,3] [,4] [,5] [,6]
## val 10.0 20.0 30.0 40.0 50.0 60.0
## fill 0.1 0.2 0.3 0.4 0.5 0.6
现在我们将val
andfill
向量相互交错,这样,当展平为向量时,当我们将它作为参数传递给内部使用它作为向量的函数时,就会发生这种情况,例如我们将使用rep()
'sx
参数,我们将以正确的顺序获取val
相应fill
的值,以便从中重建矩阵。让我通过将矩阵显式展平为一个向量来展示它的外观(如您所见,这种“展平”可以通过一个简单的c()
调用来完成):
c(t(X[,c('val','fill')]));
## [1] 10.0 0.1 20.0 0.2 30.0 0.3 40.0 0.4 50.0 0.5 60.0 0.6
所以,我们有我们的x
论点。现在我们只需要构造times
参数。
这实际上是相当棘手的。首先,我们可以认识到val
值的重复计数直接在 的rep
列中提供X
,因此我们在X[,'rep']
. 并且这些值的重复计数可以根据我在 中捕获的输出矩阵中的列数与前面提到的或 IOW的重复计数fill
之间的差异来计算。问题是,我们需要交错这两个向量以符合我们的论点。Y
YC
val
YC-X[,'rep']
x
我不知道在 R 中交错两个向量的任何“内置”方式;似乎没有任何功能可以做到这一点。在处理这个问题时,我为这项任务提出了两种不同的可能解决方案,其中一种在性能和简洁性方面似乎都更好。但是由于我编写了我的原始解决方案来使用“更差”的解决方案,并且只是后来(实际上在写这个解释时)才想到第二种和“更好”的解决方案,我将在这里解释这两种方法,从第一种和更糟糕的开始一。
交错解决方案#1
交错两个向量可以通过顺序组合向量来完成,然后用精心设计的索引向量索引该组合向量,该索引向量基本上从组合向量的前半部分来回跳跃到后半部分,依次拉出每个元素每一半以交替的方式。
为了构造这个索引向量,我从一个长度等于组合向量长度一半的顺序向量开始,每个元素重复一次:
rep(1:nrow(X),each=2);
## [1] 1 1 2 2 3 3 4 4 5 5 6 6
接下来,我添加一个二元素向量,该向量由0
组合向量的一半长度组成:
nrow(X)*0:1;
## [1] 0 6
第二个加数循环通过第一个加数,实现我们需要的交错:
rep(1:nrow(X),each=2)+nrow(X)*0:1;
## [1] 1 7 2 8 3 9 4 10 5 11 6 12
因此我们可以索引组合的重复向量来得到我们的times
论点:
c(X[,'rep'],YC-X[,'rep'])[rep(1:nrow(X),each=2)+nrow(X)*0:1];
## [1] 1 3 2 2 3 1 4 0 0 4 1 3
交错解决方案#2
交错两个向量也可以通过将这两个向量组合成一个矩阵然后再次展平它们来完成,以使它们自然地成为交错的。我相信最简单的方法是将rbind()
它们放在一起,然后立即将它们压平c()
:
c(rbind(X[,'rep'],YC-X[,'rep']));
## [1] 1 3 2 2 3 1 4 0 0 4 1 3
根据一些粗略的性能测试,似乎解决方案 #2 的性能更高,并且可以清楚地看到它更简洁。此外,可以很容易地将其他向量添加到rbind()
调用中,但是添加到解决方案 #1 会涉及更多内容(几个增量)。
性能测试(使用大型数据集):
il1 <- function() c(X[,'rep'],YC-X[,'rep'])[rep(1:nrow(X),each=2)+nrow(X)*0:1];
il2 <- function() c(rbind(X[,'rep'],YC-X[,'rep']));
identical(il1(),il2());
## [1] TRUE
system.time({ replicate(30,il1()); });
## user system elapsed
## 3.750 0.000 3.761
system.time({ replicate(30,il1()); });
## user system elapsed
## 3.810 0.000 3.815
system.time({ replicate(30,il2()); });
## user system elapsed
## 1.516 0.000 1.512
system.time({ replicate(30,il2()); });
## user system elapsed
## 1.500 0.000 1.503
所以完整的rep()
调用以正确的顺序为我们提供了我们的数据:
rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep'])));
## [1] 10.0 0.1 0.1 0.1 20.0 20.0 0.2 0.2 30.0 30.0 30.0 0.3 40.0 40.0 40.0 40.0 0.5 0.5 0.5 0.5 60.0 0.6 0.6 0.6
最后一步是使用 构建一个矩阵byrow=T
,因为这就是数据最终从rep()
. 而且我们还必须指定所需的行数,这与输入矩阵相同,XR
(或者,如果需要,我们可以指定列数YC
,或者甚至两者都指定):
Y <- matrix(rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep']))),XR,byrow=T);
Y;
## [,1] [,2] [,3] [,4]
## [1,] 10.0 0.1 0.1 0.1
## [2,] 20.0 20.0 0.2 0.2
## [3,] 30.0 30.0 30.0 0.3
## [4,] 40.0 40.0 40.0 40.0
## [5,] 0.5 0.5 0.5 0.5
## [6,] 60.0 0.6 0.6 0.6
我们完成了!