尽可能高效地在矩阵中找到一个序列
Find a sequence in a matrix as efficiently as possible
要求不高。
请在发布您的答案之前！！
1) 确保你的函数没有给出其他数据的错误，模拟几个相似的矩阵。 （关闭种子）
2)确保你的函数比我的快
3)确保你的函数和我的完全一样，在不同的矩阵上模拟（关闭种子）
例如
for(i in 1:500){
m <- matrix(sample(c(F,T),30,T),ncol = 3) ; colnames(m) <- paste0("x",1:ncol(m))
res <- c(my_fun(m),your_function(m))
print(res)
if(sum(res)==1) break
}
m
4) 该函数应该适用于具有任意数量的行和列的矩阵
======================================== ==================
该函数在逻辑矩阵的第一列中查找 true，如果找到，则转到第 2 列和新行，依此类推。
如果找到序列 return true 如果没有 false
set.seed(15)
m <- matrix(sample(c(F,T),30,T),ncol = 3) ; colnames(m) <- paste0("x",1:ncol(m))
m
x1 x2 x3
[1,] FALSE TRUE TRUE
[2,] FALSE FALSE FALSE
[3,] TRUE TRUE TRUE
[4,] TRUE TRUE TRUE
[5,] FALSE FALSE FALSE
[6,] TRUE TRUE FALSE
[7,] FALSE TRUE FALSE
[8,] FALSE FALSE FALSE
[9,] FALSE FALSE TRUE
[10,] FALSE FALSE TRUE
我的慢示例函数
find_seq <- function(m){
colum <- 1
res <- rep(FALSE,ncol(m))
for(i in 1:nrow(m)){
if(m[i,colum]==TRUE){
res[colum] <- TRUE
print(c(row=i,col=colum))
colum <- colum+1}
if(colum>ncol(m)) break
}
all(res)
}
find_seq(m)
row col
3 1
row col
4 2
row col
9 3
[1] TRUE
如何让它尽可能快？
UPD=========================
microbenchmark::microbenchmark(Jean_Claude_Arbaut_fun(m),
+ ThomasIsCoding_fun(m),
+ my_fun(m))
Unit: microseconds
expr min lq mean median uq max neval cld
Jean_Claude_Arbaut_fun(m) 2.850 3.421 4.36179 3.9915 4.5615 27.938 100 a
ThomasIsCoding_fun(m) 14.824 15.965 17.92030 16.5350 17.1050 101.489 100 b
my_fun(m) 23.946 24.517 25.59461 25.0880 25.6580 42.192 100 c
更新
追求速度的可以试试下面的base R方案
TIC_fun <- function(m) {
p <- k <- 1
nr <- nrow(m)
nc <- ncol(m)
repeat {
if (p > nr) {
return(FALSE)
}
found <- FALSE
for (i in p:nr) {
if (m[i, k]) {
# print(c(row = i, col = k))
p <- i + 1
k <- k + 1
found <- TRUE
break
}
}
if (!found) {
return(FALSE)
}
if (k > nc) {
return(TRUE)
}
}
}
你会看到
Unit: microseconds
expr min lq mean median uq max neval
my_fun(m) 18.600 26.3010 41.46795 41.5510 44.3010 121.302 100
TIC_fun(m) 10.201 14.1515 409.89394 22.6505 24.4005 38906.601 100
上一个答案
您可以试试下面的代码
lst <- with(as.data.frame(which(m, arr.ind = TRUE)), split(row, col))
# lst <- apply(m, 2, which)
setNames(
stack(
setNames(
Reduce(function(x, y) y[y > x][1],
lst,
init = -Inf,
accumulate = TRUE
)[-1],
names(lst)
)
),
c("row", "col")
)
这给出了
row col
1 3 1
2 4 2
3 9 3
一个更有趣的实现可能是使用递归（只是为了好玩，但由于效率低下而不推荐）
f <- function(k) {
if (k == 1) {
return(data.frame(row = which(m[, k])[1], col = k))
}
s <- f(k - 1)
for (i in (tail(s, 1)$row + 1):nrow(m)) {
if (m[i, k]) {
return(rbind(s, data.frame(row = i, col = k)))
}
}
}
并给出
> f(ncol(m))
row col
1 3 1
2 4 2
3 9 3
与accumulate:
purrr::accumulate(apply(m, 2, which), .init = -Inf, ~ min(.y[.y > min(.x)]))[-1]
# or
purrr::accumulate(apply(m, 2, which), .init = -Inf, ~ .y[.y > .x][1])[-1]
# x1 x2 x3
# 3 4 9
有点难看（<<- 的原因），但它会完成工作..
tempval <- 0
lapply(split(m, col(m)), function(x) {
value <- which(x)[which(x) > tempval][1]
tempval <<- value
return(value)
})
# $`1`
# [1] 3
#
# $`2`
# [1] 4
#
# $`3`
# [1] 9
如果您的示例具有代表性，我们假设 nrow(m) >> ncol(m)。在这种情况下，将交互从行移动到列会更有效：
ff = function(m)
{
i1 = 1
for(j in 1:ncol(m)) {
if(i1 > nrow(m)) return(FALSE)
i1 = match(TRUE, m[i1:nrow(m), j]) + i1
#print(i1)
if(is.na(i1)) return(FALSE)
}
return(TRUE)
}
这里是一个专注于案件处理的功能。它比所有的都快，希望它是正确的:)
f <- \(m) {
stopifnot(dim(m)[2] == 3L)
e <- nrow(m)
x1 <- if (any(xx1 <- m[, 1])) {
which.max(xx1)
} else {
NA_integer_
}
x2 <- if (is.na(x1)) {
NA_integer_
}
else if (any(xx2 <- m[(x1 + 1):e, 2])) {
which.max(xx2) + x1
} else {
NA_integer_
}
x3 <- if (is.na(x2)) {
NA_integer_
}
else if (any(xx3 <- m[(x2 + 1):e, 3])) {
which.max(xx3) + x2
} else {
NA_integer_
}
!anyNA(c(x1, x2, x3))
}
f(m)
# [1] TRUE
m2 <- m
m2[, 3] <- FALSE
f(m2)
# [1] FALSE
数据：
set.seed(15)
m <- matrix(sample(c(FALSE, TRUE), 30, TRUE), ncol=3)
如果我对问题的理解正确，那么单次循环遍历行就足够了。这是使用 Rcpp 执行此操作的一种方法。我这里只return true/false 的答案，如果你需要索引，也可以。
library(Rcpp)
cppFunction('
bool hasSequence(LogicalMatrix m) {
int nrow = m.nrow(), ncol = m.ncol();
if (nrow > 0 && ncol > 0) {
int j = 0;
for (int i = 0; i < nrow; i++) {
if (m(i, j)) {
if (++j >= ncol) {
return true;
}
}
}
}
return false;
}')
a <- matrix(c(F, F, T, T, F, T, F, F, F, F,
T, F, T, T, F, T, T, F, F, F,
T, F, T, T, F, F, F, F, T, T), ncol = 3)
a
hasSequence(a)
为了同时获取索引，以下函数 return 是一个列表，其中至少包含一个元素（名为 'found'，true 或 false），如果找到 = true，则为另一个元素，名为 'indices':
cppFunction('
List findSequence(LogicalMatrix m) {
int nrow = m.nrow(), ncol = m.ncol();
IntegerVector indices(ncol);
if (nrow > 0 && ncol > 0) {
int j = 0;
for (int i = 0; i < nrow; i++) {
if (m(i, j)) {
indices(j) = i + 1;
if (++j >= ncol) {
return List::create(Named("found") = true,
Named("indices") = indices);
}
}
}
}
return List::create(Named("found") = false);
}')
findSequence(a)
一些 link 了解 Rcpp：
- High performance functions with Rcpp、哈德莉·威克姆
- Rcpp for everyone、津田正树
- Interfacing R with C/C++，马泰奥·法西奥罗
- Rcpp Gallery - Rcpp 包的文章和代码示例
您必须至少了解一点 C 语言（最好是 C++，但对于基本用法，您可以将 Rcpp 视为 C 以及一些针对 R 数据类型的神奇语法）。第一个 link 解释了 Rcpp 类型的基础知识（向量、矩阵和列表，如何分配、使用和 return 它们）。其他 link 是很好的补充。
要求不高。
请在发布您的答案之前！！
1) 确保你的函数没有给出其他数据的错误，模拟几个相似的矩阵。 （关闭种子）
2)确保你的函数比我的快
3)确保你的函数和我的完全一样，在不同的矩阵上模拟（关闭种子）
例如
for(i in 1:500){
m <- matrix(sample(c(F,T),30,T),ncol = 3) ; colnames(m) <- paste0("x",1:ncol(m))
res <- c(my_fun(m),your_function(m))
print(res)
if(sum(res)==1) break
}
m
4) 该函数应该适用于具有任意数量的行和列的矩阵
======================================== ==================
该函数在逻辑矩阵的第一列中查找 true，如果找到，则转到第 2 列和新行，依此类推。
如果找到序列 return true 如果没有 false
set.seed(15)
m <- matrix(sample(c(F,T),30,T),ncol = 3) ; colnames(m) <- paste0("x",1:ncol(m))
m
x1 x2 x3
[1,] FALSE TRUE TRUE
[2,] FALSE FALSE FALSE
[3,] TRUE TRUE TRUE
[4,] TRUE TRUE TRUE
[5,] FALSE FALSE FALSE
[6,] TRUE TRUE FALSE
[7,] FALSE TRUE FALSE
[8,] FALSE FALSE FALSE
[9,] FALSE FALSE TRUE
[10,] FALSE FALSE TRUE
我的慢示例函数
find_seq <- function(m){
colum <- 1
res <- rep(FALSE,ncol(m))
for(i in 1:nrow(m)){
if(m[i,colum]==TRUE){
res[colum] <- TRUE
print(c(row=i,col=colum))
colum <- colum+1}
if(colum>ncol(m)) break
}
all(res)
}
find_seq(m)
row col
3 1
row col
4 2
row col
9 3
[1] TRUE
如何让它尽可能快？
UPD=========================
microbenchmark::microbenchmark(Jean_Claude_Arbaut_fun(m),
+ ThomasIsCoding_fun(m),
+ my_fun(m))
Unit: microseconds
expr min lq mean median uq max neval cld
Jean_Claude_Arbaut_fun(m) 2.850 3.421 4.36179 3.9915 4.5615 27.938 100 a
ThomasIsCoding_fun(m) 14.824 15.965 17.92030 16.5350 17.1050 101.489 100 b
my_fun(m) 23.946 24.517 25.59461 25.0880 25.6580 42.192 100 c
更新
追求速度的可以试试下面的base R方案
TIC_fun <- function(m) {
p <- k <- 1
nr <- nrow(m)
nc <- ncol(m)
repeat {
if (p > nr) {
return(FALSE)
}
found <- FALSE
for (i in p:nr) {
if (m[i, k]) {
# print(c(row = i, col = k))
p <- i + 1
k <- k + 1
found <- TRUE
break
}
}
if (!found) {
return(FALSE)
}
if (k > nc) {
return(TRUE)
}
}
}
你会看到
Unit: microseconds
expr min lq mean median uq max neval
my_fun(m) 18.600 26.3010 41.46795 41.5510 44.3010 121.302 100
TIC_fun(m) 10.201 14.1515 409.89394 22.6505 24.4005 38906.601 100
上一个答案
您可以试试下面的代码
lst <- with(as.data.frame(which(m, arr.ind = TRUE)), split(row, col))
# lst <- apply(m, 2, which)
setNames(
stack(
setNames(
Reduce(function(x, y) y[y > x][1],
lst,
init = -Inf,
accumulate = TRUE
)[-1],
names(lst)
)
),
c("row", "col")
)
这给出了
row col
1 3 1
2 4 2
3 9 3
一个更有趣的实现可能是使用递归（只是为了好玩，但由于效率低下而不推荐）
f <- function(k) {
if (k == 1) {
return(data.frame(row = which(m[, k])[1], col = k))
}
s <- f(k - 1)
for (i in (tail(s, 1)$row + 1):nrow(m)) {
if (m[i, k]) {
return(rbind(s, data.frame(row = i, col = k)))
}
}
}
并给出
> f(ncol(m))
row col
1 3 1
2 4 2
3 9 3
与accumulate:
purrr::accumulate(apply(m, 2, which), .init = -Inf, ~ min(.y[.y > min(.x)]))[-1]
# or
purrr::accumulate(apply(m, 2, which), .init = -Inf, ~ .y[.y > .x][1])[-1]
# x1 x2 x3
# 3 4 9
有点难看（<<- 的原因），但它会完成工作..
tempval <- 0
lapply(split(m, col(m)), function(x) {
value <- which(x)[which(x) > tempval][1]
tempval <<- value
return(value)
})
# $`1`
# [1] 3
#
# $`2`
# [1] 4
#
# $`3`
# [1] 9
如果您的示例具有代表性，我们假设 nrow(m) >> ncol(m)。在这种情况下，将交互从行移动到列会更有效：
ff = function(m)
{
i1 = 1
for(j in 1:ncol(m)) {
if(i1 > nrow(m)) return(FALSE)
i1 = match(TRUE, m[i1:nrow(m), j]) + i1
#print(i1)
if(is.na(i1)) return(FALSE)
}
return(TRUE)
}
这里是一个专注于案件处理的功能。它比所有的都快，希望它是正确的:)
f <- \(m) {
stopifnot(dim(m)[2] == 3L)
e <- nrow(m)
x1 <- if (any(xx1 <- m[, 1])) {
which.max(xx1)
} else {
NA_integer_
}
x2 <- if (is.na(x1)) {
NA_integer_
}
else if (any(xx2 <- m[(x1 + 1):e, 2])) {
which.max(xx2) + x1
} else {
NA_integer_
}
x3 <- if (is.na(x2)) {
NA_integer_
}
else if (any(xx3 <- m[(x2 + 1):e, 3])) {
which.max(xx3) + x2
} else {
NA_integer_
}
!anyNA(c(x1, x2, x3))
}
f(m)
# [1] TRUE
m2 <- m
m2[, 3] <- FALSE
f(m2)
# [1] FALSE
数据：
set.seed(15)
m <- matrix(sample(c(FALSE, TRUE), 30, TRUE), ncol=3)
如果我对问题的理解正确，那么单次循环遍历行就足够了。这是使用 Rcpp 执行此操作的一种方法。我这里只return true/false 的答案，如果你需要索引，也可以。
library(Rcpp)
cppFunction('
bool hasSequence(LogicalMatrix m) {
int nrow = m.nrow(), ncol = m.ncol();
if (nrow > 0 && ncol > 0) {
int j = 0;
for (int i = 0; i < nrow; i++) {
if (m(i, j)) {
if (++j >= ncol) {
return true;
}
}
}
}
return false;
}')
a <- matrix(c(F, F, T, T, F, T, F, F, F, F,
T, F, T, T, F, T, T, F, F, F,
T, F, T, T, F, F, F, F, T, T), ncol = 3)
a
hasSequence(a)
为了同时获取索引，以下函数 return 是一个列表，其中至少包含一个元素（名为 'found'，true 或 false），如果找到 = true，则为另一个元素，名为 'indices':
cppFunction('
List findSequence(LogicalMatrix m) {
int nrow = m.nrow(), ncol = m.ncol();
IntegerVector indices(ncol);
if (nrow > 0 && ncol > 0) {
int j = 0;
for (int i = 0; i < nrow; i++) {
if (m(i, j)) {
indices(j) = i + 1;
if (++j >= ncol) {
return List::create(Named("found") = true,
Named("indices") = indices);
}
}
}
}
return List::create(Named("found") = false);
}')
findSequence(a)
一些 link 了解 Rcpp：
- High performance functions with Rcpp、哈德莉·威克姆
- Rcpp for everyone、津田正树
- Interfacing R with C/C++，马泰奥·法西奥罗
- Rcpp Gallery - Rcpp 包的文章和代码示例
您必须至少了解一点 C 语言（最好是 C++，但对于基本用法，您可以将 Rcpp 视为 C 以及一些针对 R 数据类型的神奇语法）。第一个 link 解释了 Rcpp 类型的基础知识（向量、矩阵和列表，如何分配、使用和 return 它们）。其他 link 是很好的补充。