尽可能高效地查找矩阵中的序列

Question

要求很少。

在发布您的答案之前请！！

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)该函数应适用于具有任意行数和列数的矩阵

=================================================== ======== 该函数在逻辑矩阵的第一列中查找 a true，如果找到 true，则转到第 2 列和新行，依此类推。如果找到序列，如果没有找到则返回truefalse

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

Answer 1

更新

如果你追求速度，可以尝试下面的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

Answer 2

如果你的例子具有代表性，我们假设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)
}

Answer 3

如果我正确地理解了这个问题，那么对行进行一次循环就足够了。这是使用 Rcpp 执行此操作的一种方法。这里我只返回正确/错误的答案，如果你需要索引，也是可行的。

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)

为了获取索引，以下函数返回一个列表，其中至少有一个元素（名为“found”，true 或 false），如果 find = 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)

了解 Rcpp 的一些链接：

• Rcpp 的高性能函数，Hadley Wickham
• 适合所有人的 Rcpp，Masaki E. Tsuda
• R 与 C/C++ 接口，Matteo Fasiolo
• Rcpp Gallery - Rcpp 包的文章和代码示例

您必须至少了解一点 C 语言（最好是 C++，但对于基本用法，您可以将 Rcpp 视为具有 R 数据类型的一些神奇语法的 C）。第一个链接解释了 Rcpp 类型的基础知识（向量、矩阵和列表，如何分配、使用和返回它们）。其他链接是很好的补充。