Ефективно индексиране/съединяване в data.table през множество зависими условия за алгоритъм за откриване на спиране

Редактиране: Истински набор от данни е наличен тук

С благодарност към

Уанг, Руи, Фанглин Чен, Женю Чен, Тиансин Ли, Габриела Харари, Стефани Тигнор, Ся Джоу, Дрор Бен-Зеев и Андрю Т. Кембъл. StudentLife: Оценка на психичното здраве, академичното представяне и поведенческите тенденции на студенти, използващи смартфони. В сборника на конференцията на ACM за повсеместните компютри. 2014 г.


Обяснение

Провеждам симулационно проучване, в което извършвам откриване на спиране на данни за местоположение (координати на ширина/дължина) въз основа на относително прости критерии.

Местоположение (A) е спирка, ако съществува друго местоположение (B) с клеймо за време от поне 180 секунди след A и ако всички местоположения между A и B включително имат разстояние от A по-малко или равно на 80 метра.

Опитах се да намаля данните така, че все още да работи, но не изисква действителни координати.

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) има разстояние в latlon от 90.

id 2 е спирка, защото всички местоположения между него и първото местоположение с разлика във времето › 180 (id 6) имат разстояние по-малко от 80 (0, 30, 40, 40, 50).

id 6 не е спирка, защото въпреки че id 10 е › 180 разлика във времето, id 7, което попада между, има разстояние, по-голямо от 80.

id 8 не е спирка, защото няма местоположение след поне 180 секунди след това.

В крайна сметка трябва да мога да присвоя идентификатор на спирка лакомо, ако открия например, че id 2 има точки, които отговарят на изискването за разстояние до id 7, местоположения с идентификатори 2:7 имат идентификатор на спирка 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 цикъл, за да задам идентификатора на спиране, както е подходящо.

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 пъти в моето симулационно проучване.

Моето предположение е, че най-добрият начин да се справите с това е с non-equi съединение в data.table.

Другото изпълнение, което изпълнявам в момента, функционира по-добре, когато броят на редовете в набора от данни прави масива твърде тежък за памет. Няма да превеждам това, за да работя върху данните, но е тук, в случай че даде на някой идеи. Сложих го в цикъл while, за да може да прескочи някои итерации, когато вече му е присвоен stop_id на определен брой точки. Ако всички точки 1:7 принадлежат на stop_id 1, те самите не се считат за кандидат спирки, просто преминаваме към тестване отново в точка 8. Технически връща различно решение, но спирките, които са достатъчно близки, се обединяват по-късно в този процес , така че крайният резултат едва ли ще се различава много.

За цикъл, без матрица

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)]


person Danielle McCool    schedule 21.04.2021    source източник


Отговори (1)


Използвайки възможностите за non-equi joins на data.table, бихте могли да присъедините данните към себе си, като същевременно избягвате декартово произведение, което би било твърде скъпо.

Тъй като data.table позволява само >, < или =, съединенията първо се правят на правоъгълни области, преди да се филтрират подходящите разстояния. На реалните данни, които предоставихте, това прави 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
person Waldi    schedule 26.04.2021
comment
Да, мисля, че това е моментът, в който опростяването на latlon в една променлива в името на простотата е проблематично. Разбира се, когато имате координатна двойка географска ширина и дължина, не можете да видите радиус, както можете с времето. Мога да прозорча широчина и дължина независимо, но това ще ми даде квадрат. Ще актуализирам въпроса си с реални данни. - person Danielle McCool; 26.04.2021
comment
с нетърпение за реални данни: прозорецът трябва да е възможен и върху x/y данни - person Waldi; 26.04.2021
comment
Актуализирах с реални данни. За съжаление, това го прави по-малко предсказуем, но изглежда не мога да измисля добър баланс поради естеството на този тип данни. - person Danielle McCool; 26.04.2021
comment
уведомете ме дали редакцията ми върви в правилната посока. Броят на откриванията изглежда важен, но мисля, че спазват правилата, които сте посочили. - person Waldi; 26.04.2021
comment
lat/long min/max са квадратни (imgur.com/JRNsxFI), но на теория, тъй като квадратът съдържа кръг, може да намали броя, който проверявам. Проверявате истинското разстояние независимо. Не мога наистина да кажа какво трябва да бъде стоп изходът. Ако е B.min, спиранията не са алчни. Ако id1 формира потенциален стоп с id2 и id3, а id2 формира потенциален стоп с id3, id1 грабва 2 и id2 грабва 3, което предпазва стоповете от агрегиране в достатъчна степен. Също така B.all страда от същия проблем, че може да стане твърде голям, защото се присъединява преди проверка на условията, което трябва да избягвам - person Danielle McCool; 26.04.2021
comment
Разделителната способност на спиране с параметрите от 80 метра и 180 секунди на данните връща 353 спирания със stopFinder и 356 спирания с по-ясния алгоритъм, поставен по-горе. - person Danielle McCool; 26.04.2021
comment
Вижте моята актуализация. Прав сте, това все още не е алчно, ще го проверя. - person Waldi; 26.04.2021
comment
Нека продължим тази дискусия в чата. - person Danielle McCool; 26.04.2021
comment
@DanielleMcCool, онлайн в чата - person Waldi; 26.04.2021