向量化包含循环和if子句的搜索函数

Question

我给出了两个非常大的数据集，并且我一直在尝试构建一个函数，该函数将从一个集合中找出某些关于其他数据的if子句的某些坐标组。 我的问题是，我写的功能很慢虽然我一直在阅读某些问题类似的问题的答案，但我还没有设法使其工作。
所以，如果我给出：

>head(CTSS)  
    V1  V2  V3 
1 chr1 564563 564598 
2 chr1 564620 564649 
3 chr1 565369 565404 
4 chr1 565463 565541 
5 chr1 565653 565697 
6 chr1 565861 565922 

和

> head(href) 
    chr  region start  end strand nu gene_id transcript_id 
1 chr1 start_codon 67000042 67000044  + . NM_032291  NM_032291 
2 chr1   CDS 67000042 67000051  + 0 NM_032291  NM_032291 
3 chr1  exon 66999825 67000051  + . NM_032291  NM_032291 
4 chr1   CDS 67091530 67091593  + 2 NM_032291  NM_032291 
5 chr1  exon 67091530 67091593  + . NM_032291  NM_032291 
6 chr1   CDS 67098753 67098777  + 1 NM_032291  NM_032291 

对于从HREF数据集我想找到的前两个值在起始列每个值CTSS数据集的第3列小于或等于，并将其保留在新的数据框中。
环路我写道：

y <- CTSS[order(-CTSS$V3), ]  
find_CTSS <- function(x, y) { 
    n <- length(x$start) 
    foo <- data.frame(matrix(0, n, 6)) 
    for (i in 1:n) 
    { 
     a <- which(y$V3 <= x$start[i]) 
     foo[i, ] = c(x$start[i], x$stop[i], y$V2[a[1]], y$V3[a[1]] , y$V2[a[2]], y$V3[a[2]]) 
    } 

print(foo) 

} 

Answer 1

您提供很少的数据（but see here），所以这是一个有点难以基准您的解决方案。看看下面的解决方案是否满足您的需求。

#make some fake data 
href <- data.frame(start = runif(10), stop = runif(10), other_col = sample(letters, 10)) 
CTSS <- data.frame(col1 = runif(100), col2 = runif(100)) 

# for each row in href (but extract only stop and start columns) 
result <- apply(X = href[, c("start", "stop")], MARGIN = 1, FUN = function(x, ctss) { 
      criterion <- x["start"] #make a criterion 
      #see which values are smaller or equal to this criterion (and sort them) 
      extracted <- sort(ctss[ctss$col2 <= criterion, "col2"]) 
      #extract last and one to last value 
      get.values <- extracted[c(length(extracted) - 1, length(extracted))] 
      #put values in data frame 
      out <- as.data.frame(matrix(get.values, ncol = 2)) 
      return(out) 
     }, ctss = CTSS) 

#pancake a list into a data.frame 
result <- do.call("rbind", result) 

Answer 2

我不知道我会花多少时间专注于此，所以我会向前迈进。当这类问题在APL杂志上收到单行答案时，我是一名APL家伙。后来，我成为了一名C++/STL人，并以新的着装代码学习了所有相同的东西。有时R会让我认为APL与PHP配合。

在这个问题中，数据框是一个分心。这是一个简单的矢量搜索，有些粘贴在一起。

对于性能关键向量搜索，您需要findInterval。搜索范围需要订购。 search-fors可以以任意顺序，但对于大型列表，您需要订购。

V <- sort (runif(10*1000*1000)) 
    U <- sort (runif(10*1000*1000)) 
    W <- findInterval (U, V) 

这运行在三个羊羔尾巴的摇晃。现在你有一对整数。第一列是1:length(U)，第二列的值是W内的整数索引。

sum(u==sort(u)[sort.int (sort.int (u, index.return=TRUE)$ix, index.return=TRUE)$ix]) 

好的，我的APL脑干有贡献。答案是长度（u），并演示了“粘合在一起”所需的逆向排序。

令人兴奋的事实：只有R中sort函数的特殊情况返回索引向量。在APL中，这是你从成绩上升/成绩下降的唯一答案。但是，嘿，这不像他们第一次做对了。

您必须修改findInterval的结果才能在匹配位置的小于一边选择两个元素，并且必须撤消两种类型才能粘合在一起。我怀疑你的运行时将会被这两种类型（非常长的列表）所占据，或者组装最终的数据框（用于较小的列表）。在我的系统上，对长度为100 * 1000 * 1000的数字列表进行排序开始出现问题。

夹在中间的findInterval的运行时间将是一片薄薄的生菜，这让我想起为什么我不打算闲逛。

Answer 3

我看到你想要的主要是加速。借用Roman Lustrik的代码，我看不出有什么优势可以在应用程序中进行排序。这真的会减慢速度。事实上，你想尽可能多地从apply（循环）中获得。所以下面的运行速度要快得多。

#all code using Roman Lustrik's made up data 

CTSSs <- CTSS[order(CTSS$col2),] #presort CTSS 
result <- lapply(X = href$start, FUN = function(x, ctss) { 
    extracted <- ctss$col2[ctss$col2 <= x] 
    get.values <- tail(extracted,2) 
    out <- matrix(get.values, ncol = 2) 
    return(out)}, ctss = CTSSs) 
#pancake a list into a data.frame 
result <- as.data.frame(do.call("rbind", result)) 

或者，我可以进一步遵循向量化的精神，真正地获得尽可能小的功能。

CTSSs <- CTSS[order(CTSS$col2),] #presort CTSS 
extracted <- lapply(href$start, function(x, ctss) { 
    ctss$col2[ctss$col2 <= x]}, ctss = CTSSs) 
get.values <- lapply(extracted, tail, n = 2) 
result <- t(sapply(get.values, matrix, ncol = 2)) 

#convert to a data.frame 
result <- as.data.frame(result) 

这可能会更快，或者你的情况也许没有，但是，你应该需要添加一个中间步骤，这也可能采取真正的矢量内置函数的优势，说，如果你想在做数学在将它们放入数据框之前的值，然后您可以轻松地做到这一点。另外，你会注意到现在我可以在矩阵阶段使用一个sapply/transpose，它比lapply/rbind更快。这通常是加速向量化你的代码的地方，而不是仅仅围绕它做一个大循环。 （顺便说一句，它可以更容易地在你的思维的每一步检查错误...也许这不是一个人呢。）

修订：

在进一步的思考，我意识到这是可以完全矢量。下面的代码将生成你想要比任何前面的例子快得多的东西。诀窍是使用cut（）和aggregate（）命令。

href <- href[order(href$start),] #just sorted so that the 0 at the beginning makes sense and the labels then match 
margin <- cut(CTSS$col2, breaks = c(0,href$start), labels = href$start, right = TRUE) 
result <- aggregate(col2 ~ margin, data = CTSS, FUN = function(x) tail(x,2)) 

可以重新格式化，结果如你所愿得到你想要什么，但应该做它的肉。您可能需要将margin列更改为numeric，以便它与href $ start匹配，并在上面的中间示例中使用与sapply类似的代码将上面的项目对列表转换为两个单独的列。这是循环或应用语句中的if（）语句，之前放慢了你的速度，而cut（）消除了这一点。