Menu
I am trying to take the results of a which(..., arr.ind = TRUE) function and remove the rows that are not the first to "connect" with one another.
Examples:
#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
These examples can be read in using:
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))
The purpose of this is to take a dist matrix of Euclidean distances and turn it into a sequence of point-to-point distances that skip distances below a certain threshold. While there may be other ways to solve this problem, I am very interested in figuring out the best way to do this by filtering out rows from the which-matrix.
Reproducible example of my intended use:
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])
874
This is a perfect problem for Rcpp. Observe:
#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
Here are the benchmarks on example3 provided by the OP:
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
Here are some benchmarks on a larger example (didn't include get_graph as it was taking a very long time):
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
Test for equality:
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
UPDATE
If Rcpp isn't your flavor, here is a Base R translation of the above code that is still faster than the other solutions:
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
And for the larger example:
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
