r - 保持连接的结果(...，arr.ind = TRUE)

Question

我正在尝试获取which(..., arr.ind = TRUE)函数的结果，并删除不是第一个“彼此连接”的行。
示例:

#example 1      example 2      example 3
   row col        row col        row col
     1   4          2   3          1   3
     2   4          2   4          2   5
     4   5          3   5          3   5
     3   6          2   7          4   6
     4   6          3   7          5   6
     3   7          4   7          6   8
     4   7          5   7          9  10
# should become (trimmed.mtx)
   row col        row col        row col
     1   4          2   3          1   3
     4   5          3   5          3   5
                    5   7          5   6
                                   6   8

可以使用以下示例读取这些示例:

example1 <- structure(list(row = c(1L, 2L, 4L, 3L, 4L, 3L, 4L), col = c(4L, 4L, 5L, 6L, 6L, 7L, 7L)), .Names = c("row", "col"), class = "data.frame", row.names = c(NA, -7L))
example2 <- structure(list(row = c(2L, 2L, 3L, 2L, 3L, 4L, 5L), col = c(3L, 4L, 5L, 7L, 7L, 7L, 7L)), .Names = c("row", "col"), class = "data.frame", row.names = c(NA, -7L))
example3 <- structure(list(row = c(1L, 2L, 3L, 4L, 5L, 6L, 9L), col = c(3L, 5L, 5L, 6L, 6L, 8L, 10L)), .Names = c("row", "col"), class = "data.frame", row.names = c(NA, -7L))

这样做的目的是获取一个欧几里得距离的距离矩阵，并将其转换为一系列点对点距离，从而跳过距离低于某个阈值的距离。尽管可能还有其他方法可以解决此问题，但我对通过从which-matrix过滤掉行来找出解决此问题的最佳方法非常感兴趣。
我的预期用途的可复制示例:

set.seed(81417) # Aug 14th, 2017
# Generate fake location data (temporally sequential)

x <- as.matrix(cbind(x = rnorm(10, 10, 3), y =  rnorm(10, 10, 3)))

# Find euclidean point-to-point distances and remove distances that are less than:
value = 5
# I attempted to do so by calculating an entire Euclidean distance matrix (dist())
# and then finding a path from point-to-nearest-point 
# using distances that are greater than the value 
d <- as.matrix(dist(x[,c("x","y")]))
d[lower.tri(d)] <- 0
mtx <- which(d > value, arr.ind = T)
mtx

# Change from EVERY point-to-point distance (mtx) > value
# to only the "connecting" points that exceed the skipping value
trimmed.mtx <- {?}

# final result
cbind(x[unique(c(trimmed.mtx)),],d[trimmed.mtx])

Answer 1

对于Rcpp来说，这是一个完美的问题。观察:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerMatrix findConnections(IntegerMatrix m) {
    int i = 0, j = 0, k = 1, n = m.nrow();

    // initialize matrix with same dimensions as m
    IntegerMatrix myConnections(n, 2);

    while (i < n) {
        // Populate with "connected" row
        myConnections(j,_) = m(i,_);

        // Search for next connection
        while (k < n && m(i, 1) != m(k, 0)) {k++;}
        i = k;
        j++;
    }

    // Subset matrix and output result
    IntegerMatrix subMatrix(j, 2);
    for (i = 0; i < j; i++) {subMatrix(i,_) = myConnections(i,_);}

    return subMatrix;
}


findConnections(as.matrix(example3))
     [,1] [,2]
[1,]    1    3
[2,]    3    5
[3,]    5    6
[4,]    6    8

以下是OP提供的 example3基准测试:

microbenchmark(get_path(example3),
               foo(example3),
               f(example3),
               findConnections(as.matrix(example3)))
Unit: microseconds
                                expr      min        lq       mean   median        uq        max neval cld
                  get_path(example3) 3345.999 3519.0255 6361.76978 3714.014 3892.9930 202511.942   100   b
                       foo(example3)  215.514  239.3230  360.81086  257.180  278.3200  10256.384   100  a 
                         f(example3)  936.355 1034.4645 1175.60323 1073.668 1142.4270   9676.755   100  a 
findConnections(as.matrix(example3))   52.135   60.3445   71.62075   67.528   80.4585    103.858   100  a 

以下是一个较大示例的一些基准测试(由于耗时很长，因此未包含 get_graph):

set.seed(6221)
x <- as.matrix(cbind(x = rnorm(1000, 10, 3), y =  rnorm(1000, 10, 3)))
value = 5
d <- as.matrix(dist(x[,c("x","y")]))
d[lower.tri(d)] <- 0
mtxLarge <- which(d > value, arr.ind = T)
mtxLargeFoo <- data.frame(mtxLarge, row.names = NULL) ## this is for the function foo
                                            ## as we don't want to include
                                            ## the time it takes to create
                                            ## a data.frame every time.

microbenchmark(foo(mtxLargeFoo),
               f(mtxLarge),
               findConnections(as.matrix(mtxLarge)), times = 10, unit = "relative")
Unit: relative
                                expr      min       lq     mean   median       uq      max neval cld
                    foo(mtxLargeFoo) 3168.479 3376.909 2660.377 3424.276 2319.434 1960.161    10  b 
                         f(mtxLarge) 8307.009 8436.569 6420.919 8319.151 5184.557 4610.922    10   c
findConnections(as.matrix(mtxLarge))    1.000    1.000    1.000    1.000    1.000    1.000    10 a  

测试是否相等:

a <- findConnections(as.matrix(mtxLarge))
b <- foo(mtxLargeFoo)
c <- f(mtxLarge)
sapply(1:2, function(x) identical(a[,x], b[,x], c[, x]))
[1] TRUE TRUE

更新
如果 Rcpp不是您的喜好，这是上述代码的Base R转换，它仍然比其他解决方案快:

findConnectionsBase <- function(m) {
    n <- nrow(m)
    myConnections <- matrix(integer(0), nrow = n, ncol = 2)
    i <- j <- 1L
    k <- 2L
    while (i <= n) {
        myConnections[j, ] <- m[i, ]
        while (k <= n && m[i, 2] != m[k, 1]) {k <- k + 1L}
        i <- k
        j <- j + 1L
    }
    myConnections[!is.na(myConnections[,1]), ]
}

microbenchmark(get_path(example3),
           foo(example3),
           f(example3),
           BaseR = findConnectionsBase(as.matrix(example3)),
           Rcpp = findConnections(as.matrix(example3)))
Unit: microseconds
              expr      min        lq       mean   median       uq        max neval cld
get_path(example3) 3128.844 3204.3765 6057.18995 3406.137 3849.274 188685.016   100   b
     foo(example3)  239.734  251.4325  399.71418  267.648  301.309  12455.441   100  a 
       f(example3)  899.409  961.3950 1145.72695 1014.555 1127.237   9583.982   100  a 
             BaseR   79.638   89.2850  103.63571   97.905  111.657    212.230   100  a 
              Rcpp   48.850   55.8290   64.24807   61.781   69.170    123.151   100  a 

对于较大的示例:

microbenchmark(foo(mtxLargeFoo),
               f(mtxLarge),
               BaseR = findConnectionsBase(as.matrix(mtxLarge)),
               Rcpp = findConnections(as.matrix(mtxLarge)), times = 10, unit = "relative")
               Unit: relative
            expr       min        lq      mean    median        uq       max neval  cld
foo(mtxLargeFoo) 2651.9626 2555.0515 1606.2785 1703.0256 1711.4850  671.9115    10   c 
     f(mtxLarge) 6812.7195 6433.2009 3976.6135 4218.1703 4105.1138 1642.2768    10    d
           BaseR  787.9947  733.4528  440.2043  478.9412  435.4744  167.7491    10  b  
            Rcpp    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    10 a