r - 跨多个相关条件的数据表中的高效索引/加入停止检测算法
问题描述
编辑：此处提供真实数据集
感谢
Wang, Rui, Fanglin Chen, Zhenyu Chen, Tianxing Li, Gabriella Harari, Stefanie Tignor, Xia Zhou, Dror Ben-Zeev, and Andrew T. Campbell。“StudentLife：使用智能手机评估大学生的心理健康、学习成绩和行为趋势。” 在 ACM 普适计算会议论文集上。2014 年。
解释
我正在进行一项模拟研究，其中我根据相对简单的标准对位置数据（纬度/经度坐标）执行停止检测。
如果在 A 之后存在时间戳至少为 180 秒的另一个位置 (B)，并且如果 A 和 B 之间的所有位置与 A 的距离小于或等于 80 米，则位置 (A) 是一个停靠点。
我试图减少数据，使其仍然有效，但不需要实际坐标。
data <- data.table(id = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
latlon = c(0, 50, 80, 90, 90, 100, 190, 110, 110, 110),
time = c(0, 60, 120, 180, 240, 300, 360, 420, 480, 520))
id 1不是停止，因为时间差 > 180 ( id 5) 的第一个位置在纬度中的距离为 90。
id 2是一个停靠点，因为其自身与第一个位置之间的时间差 > 180 ( id 6) 的所有位置的距离都小于 80 (0, 30, 40, 40, 50)。
id 6不是一个停止，因为即使id 10是 > 180 的时间差异，id 7它之间的距离大于 80。
id 8不是停止，因为之后至少 180 秒后没有位置。
最终，我需要能够贪婪地分配“stop id”，例如，如果我发现id 2有满足距离要求的点直到id 7，id 为 2:7 的位置的 stop id 为 2。
矩阵和for循环
如果我运行这个：
nrows <- nrow(data)
latlon_dist <- outer(data$latlon, data$latlon, `-`)
latlon_dist[upper.tri(latlon_dist)] <- NA
time_window <- outer(data$time, data$time, `-`)
time_window[upper.tri(time_window)] <- NA
foo <- function(x){
mindist <- min(which(x[, 1] > 80), nrows)
if (mindist >= min(which(x[, 2] > 180), nrows + 1)) mindist else NA
}
bar <- array(c(latlon_dist, time_window),
dim = c(nrows, nrows, 2))
apply(bar, 2, foo)
它给了> NA 7 7 NA NA NA NA NA NA NA我可以在 for 循环中使用的阈值来适当地设置停止 ID。
threshholds <- apply(bar, 2, foo) - 1
previous_threshhold <- 0
for (i in seq_along(threshholds)) {
current_threshhold <- threshholds[i]
if (!is.na(current_threshhold) && current_threshhold > previous_threshhold) {
data[i:current_threshhold, stop_id := i]
previous_threshhold <- current_threshhold
}
}
在这一点上，这是我能够保证准确性的唯一方法。我尝试过的所有其他事情我都认为是正确的，只是发现它的行为与这种情况不同。但正如您可能想象的那样，这是非常低效的，并且在我的模拟研究中运行了 116,000 次。
我的假设是处理此问题的最佳方法是在 data.table 中使用非 equi 连接。
当数据集中的行数使数组过于占用内存时，我目前正在运行的另一个实现更好。我不会将其翻译为处理数据，但它就在这里，以防它给任何人提供想法。我把它卡在了一个while循环中，所以当它已经为多个点分配了一个stop_id时，它可以跳过一些迭代。如果点 1:7 都属于 stop_id 1，则它们本身不被视为候选停靠点，我们只需在点 8 处再次进行测试。它在技术上返回不同的解决方案，但“足够接近”的停靠点稍后会合并这个过程，所以最终的结果不太可能相差太大。
For循环，没有矩阵
stopFinder <- function(dt){
nrows <- nrow(dt)
if (nrows < 20000){
return(quickStopFinder(dt))
}
i <- 1
remove_indices <- 0
while (i < nrows) {
this_ends <- dt[!remove_indices,
Position(
geodist_vec(rep(longitude[1], .N),
rep(latitude[1], .N),
longitude,
latitude,
paired = TRUE),
f = function(x) x > 80,
nomatch = .N + 1) ] + i - 1
# A) Do some number of points occur within the distance?
# B) If so, is it at least three minutes out?
if (this_ends > (i + 1) && dt[c(i, (this_ends - 1)), timestamp[.N] > time_window[1]]) {
# Last index is the one before the distance is broken
last_index_of_stop <- this_ends - 1
# Next run, we will remove all prior considerations
remove_indices <- c(1:last_index_of_stop)
# Set the point itself
dt[i,
`:=`(candidate_stop = TRUE,
stop_id = id,
within_stop = TRUE)]
# Set the attached points
dt[(i + 1):last_index_of_stop,
`:=`(within_stop = TRUE,
stop_id = i)]
# Start iterating again on the point that broke the distance
i <- this_ends
} else {
# If no stop, move on and leave out this point
remove_indices <- c(1:i)
i <- i + 1
}
}
dt[]
}
quickStopFinder 或多或少是我一开始分享的实现，它是内存密集型且速度慢，但比 stopFinder 慢一点。
以前，我有这样的东西作为基础，但它需要很多后续步骤，并且并不总是给我想要的结果，但我会为后代添加它。
res <- dt[dt,
on = .(timestamp >= timestamp_dup,
timestamp <= time_window)]
res[, dist := geodist_vec(x1 = longitude,
y1 = latitude,
x2 = i.longitude,
y2 = i.latitude,
paired = TRUE,
measure = "haversine")]
res[, candidate_stop := all(dist <= 80), i.id]
新的真实数据
使用来自真实数据的示例进行编辑：
这可以处理连接的情况，但会变得太大太快。数据量小的时候速度很快。
sm2 <- read.csv(file = "http://daniellemc.cool/sm.csv", row.names = NULL)
sm <- copy(sm2)
setDT(sm)
sm <- sm[, .(timestamp, longitude, latitude, id)]
sm[, timestamp := as.POSIXct(timestamp)]
sm[, id2 := id]
# This is problematic on my data because of how quickly it grows.
test <- sm[sm, on = .(id >= id)]
test[, i.id2 := NULL]
setnames(test, c("time.2", "longitude.2", "latitude.2", "id.1",
"id.2", "time.1", "longitude.1", "latitude.1"))
# Time and distance differences calculated between each pair
test[, distdiff := geodist_vec(longitude.1, latitude.1,
longitude.2, latitude.2,
paired = TRUE)]
test[, timediff := time.2 - time.1]
# Include the next distance to make sure there's at least one within distance and
# over 180 timediff.
test[, nextdistdiff := shift(distdiff, -1), id.1]
# Are all distances within 180 sec within 80, and is the next following also < 80
test[, dist_met := FALSE]
test[timediff < 180, dist_met := all(distdiff < 80 & nextdistdiff < 80), id.1]
test[, dist_met := any(dist_met), id.1]
# Test how many occur consecutively
# This keeps us from having > 80 dist but then coming back within 80
test[, consecutive := FALSE]
test[distdiff < 80, consecutive := c(TRUE, cummin(diff(id.2) == 1) == 1), id.1]
test[consecutive == TRUE & dist_met == TRUE, stop_id := min(id.1), id.2]
test[test[consecutive == TRUE & dist_met == TRUE], stop_id := i.stop_id, on = .(id.1 = id.2)]
test <- unique(test[, .(stop_id, id.1)])
# Join it back to the data.
sm[test, stop_id := stop_id, on = .(id = id.1)]
解决方案
使用 的非等值连接功能data.table，您可以将数据连接到自身，同时避免过于昂贵的笛卡尔积。
As data.tableonly allowed或，首先在矩形区域上进行连接，然后过滤出适当的距离>。根据您提供的真实数据，计算量减少了 7 倍。<=
library(data.table)
library(geosphere)
data <- copy(sm)
minduration <- 180
maxdistance <- 80
data[,latmin := destPoint(cbind(longitude,latitude),b = 180, d=maxdistance)[,2]]
data[,latmax := destPoint(cbind(longitude,latitude),b = 0 , d=maxdistance)[,2]]
data[,lonmin := destPoint(cbind(longitude,latitude),b = 270, d=maxdistance)[,1]]
data[,lonmax := destPoint(cbind(longitude,latitude),b = 90, d=maxdistance)[,1]]
data[,timestampmin := timestamp+minduration]
# Cross product with space and time windows
cross <- data[data,.(i.id,x.id,i.latitude,i.longitude,x.latitude,x.longitude,dist = distGeo(cbind(x.longitude,x.latitude),cbind(i.longitude,i.latitude)) ,i.timestamp,x.timestamp)
,on=.(timestamp>timestampmin,
longitude >= lonmin,
longitude<=lonmax,
latitude >= latmin,
latitude <= latmax)][
dist<maxdistance]
# Summarizing the results
cross[,.(keep=cumsum(fifelse(diff(x.id-i.id)==1,1,NA_integer_))),by=i.id][
!is.na(keep),.(startid = i.id,nextid = i.id+keep)][
!(startid %in% nextid)][
,.(maxid=max(nextid)),by=startid][
,.(stopid = min(startid)),by=maxid]
maxid stopid
1: 6 1
2: 18 10
3: 26 22
4: 33 28
5: 48 40
---
162: 4273 4269
163: 4276 4274
164: 4295 4294
165: 4303 4301
166: 4306 4305
推荐阅读
- sql - 拒绝特定用户对一组表的选择权限
- ios - 在 iOS 中使用 SwipeGestureRecognizer 时如何在视图/页面之间添加转换
- outlook - python在打开带有附件的Outlook客户端时重命名文件名
- python - 如何找到一列的最早时间戳和最新时间戳之间的时间窗口并将其按另一列分组？
- ffmpeg - 为什么 ffmpeg 刻录中文字幕(ass) 不换行？
- angular - Angular Material：弹出窗口：在窗口打开并单击按钮时将数据发送给父级
- javascript - 在 JavaScript 中捕获 true
- scroll - Angular8滚动到底部不适用于Safari 12
- html - 我正在尝试使用添加图像，但它在节点 js 中显示 null 的图像
- r - 解析时间并转换为分钟