Feat/directional correlation delay

R/edge_delay.R

@@ -0,0 +1,206 @@
+#' Directional correlation delay based edge lists
+#'
+#' \code{edge_delay} returns edge lists defined by the directional correlation
+#' delay between individuals. The function expects a \code{data.table} with
+#' relocation data, distance based edge lists, individual identifiers and a window argument. The
+#' window argument is used to specify the temporal window within which to consider
+#' the directional correlation delay. Relocation data should be in two columns
+#' representing the X and Y coordinates.
+#'
+#' The \code{edges} and \code{DT} must be \code{data.table}s. If your data is a
+#' \code{data.frame}, you can convert it by reference using
+#' \code{\link[data.table:setDT]{data.table::setDT}}.
+#'
+#' The \code{edges} and \code{DT} are internally matched in this function using
+#' the columns \code{timegroup} (from \code{group_times}) and \code{ID1} and
+#' \code{ID2} (in \code{edges}, from \code{dyad_id}) with \code{id} (in
+#' \code{DT}). This function expects a \code{fusionID} present, generated with
+#' the \code{fusion_id} function.
+#' The \code{id}, and \code{direction} arguments expect the names
+#' of a column in \code{DT} which correspond to the id, and direction columns.
+#'
+#' @inheritParams centroid_fusion
+#' @inheritParams direction_group
+#' @param window temporal window in unit of timegroup column generated with
+#'   \code{group_times}, eg. \code{window = 4} corresponds to the 4 timegroups
+#'   before and after the focal observation
+#'
+#' @return \code{edge_delay} returns the input \code{edges} appended with
+#'   a 'dir_corr_delay' column indicating the temporal delay (in units of
+#'   timegroups) at which ID1's direction of movement is most similar to
+#'   ID2's direction of movement, within the temporal window defined.
+#'
+#' @export
+#'
+#' @family Edge-list generation
+#'
+#' @examples
+#' # Load data.table
+#' library(data.table)
+#' \dontshow{data.table::setDTthreads(1)}
+#'
+#' # Read example data
+#' DT <- fread(system.file("extdata", "DT.csv", package = "spatsoc"))
+#'
+#' # Select only individuals A, B, C for this example
+#' DT <- DT[ID %in% c('A', 'B', 'C')]
+#'
+#' # Cast the character column to POSIXct
+#' DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+#'
+#' # Temporal grouping
+#' group_times(DT, datetime = 'datetime', threshold = '20 minutes')
+#'
+#' # Calculate direction
+#' direction_step(
+#'   DT = DT,
+#'   id = 'ID',
+#'   coords = c('X', 'Y'),
+#'   projection = 32736
+#' )
+#'
+#' # Distance based edge list generation
+#' edges <- edge_dist(
+#'   DT,
+#'   threshold = 100,
+#'   id = 'ID',
+#'   coords = c('X', 'Y'),
+#'   timegroup = 'timegroup',
+#'   returnDist = TRUE,
+#'   fillNA = FALSE
+#' )
+#'
+#' # Generate dyad id
+#' dyad_id(edges, id1 = 'ID1', id2 = 'ID2')
+#'
+#' # Generate fusion id
+#' fusion_id(edges, threshold = 100)
+#'
+#' # Directional correlation delay
+#' delay <- edge_delay(
+#'   edges = edges,
+#'   DT = DT,
+#'   window = 3,
+#'   id = 'ID'
+#' )
+#'
+#' print(delay)
+edge_delay <- function(
+    edges,
+    DT,
+    window = NULL,
+    id = NULL,
+    direction = 'direction') {
+
+  if (is.null(DT)) {
+    stop('input DT required')
+  }
+
+  if (is.null(edges)) {
+    stop('input edges required')
+  }
+
+  if (is.null(id)) {
+    stop('id column name required')
+  }
+
+  check_cols_edges <- c('ID1', 'ID2', 'timegroup')
+  if (any(!(check_cols_edges %in% colnames(edges)))) {
+    stop(paste0(
+      as.character(paste(setdiff(
+        check_cols_edges,
+        colnames(edges)
+      ), collapse = ', ')),
+      ' field(s) provided are not present in input DT'
+    ))
+  }
+
+  check_cols_DT <- c(id, 'timegroup', direction)
+  if (any(!(check_cols_DT %in% colnames(DT)
+  ))) {
+    stop(paste0(
+      as.character(paste(setdiff(
+        check_cols_DT,
+        colnames(DT)
+      ), collapse = ', ')),
+      ' field(s) provided are not present in input DT'
+    ))
+  }
+
+  if (is.null(window)) {
+    stop('window is required')
+  }
+
+  if (!is.numeric(window)) {
+    stop('window should be a numeric, in the units of timegroup')
+  }
+
+  if (!'fusionID' %in% colnames(edges)) {
+    stop('fusionID field not present in edges, did you run fusion_id?')
+  }
+
+  if (!'dyadID' %in% colnames(edges)) {
+    stop('dyadID field not present in edges, did you run dyad_id?')
+  }
+
+
+  if ('window' %in% colnames(DT)) {
+    setnames(DT, 'window', 'zzz_window')
+  }
+
+  if ('window' %in% colnames(edges)) {
+    setnames(edges, 'window', 'zzz_window')
+  }
+
+  data.table::setorderv(DT, 'timegroup')
+
+  id_tg <- edges[!is.na(fusionID), .(
+    tg = unique(timegroup),
+    dyadID = unique(dyadID),
+    ID1 = first(ID1),
+    ID2 = first(ID2)
+  ), by = c('fusionID')]
+
+  id_tg[, min_tg :=
+          data.table::fifelse(tg - window < min(tg), min(tg), tg - window),
+        by = c('fusionID')]
+
+  id_tg[, max_tg :=
+          data.table::fifelse(tg + window < min(tg), min(tg), tg + window),
+          # TODO: check max_tg
+          # data.table::fifelse(tg + window > max(tg), max(tg), tg + window),
+        by = c('fusionID')]
+
+  id_tg[, delay_tg := {
+    focal_direction <- DT[timegroup == .BY$tg &
+                            id == ID1, direction]
+    DT[between(timegroup, min_tg, max_tg) & id == ID2,
+       timegroup[which.min(delta_rad(focal_direction, direction))],
+       env = list(id = 'id')]
+  },
+  by = c('tg',  'dyadID')]
+
+  id_tg[, dir_corr_delay := tg - delay_tg]
+
+  data.table::setnames(id_tg,  'tg', 'timegroup')
+  data.table::set(id_tg, j = c('min_tg', 'max_tg','delay_tg'), value = NULL)
+
+  data.table::setorder(id_tg, timegroup, ID1, ID2, dir_corr_delay)
+
+  out <- data.table::rbindlist(list(
+    id_tg,
+    id_tg[, .(timegroup,  dyadID, fusionID,
+              ID1 = ID2, ID2 = ID1, dir_corr_delay = - dir_corr_delay)]
+  ), use.names = TRUE)
+
+  if ('zzz_window' %in% colnames(DT)) {
+    setnames(DT, 'zzz_window', 'window')
+  }
+
+  if ('zzz_window' %in% colnames(edges)) {
+    setnames(edges, 'zzz_window', 'window')
+  }
+
+  return(out)
+}
+

R/internal.R

@@ -0,0 +1,30 @@
+#' Difference of two angles measured in radians
+#'
+#' Internal function
+#'
+#' @param target angle in radians
+#' @param source angle in radians
+#' @param signed boolean if signed difference should be returned, default FALSE
+#'
+#' @return
+#' @references adapted from https://stackoverflow.com/a/7869457
+#'
+#' @examples
+delta_rad <- function(target, source,  signed = FALSE) {
+  if (!inherits(target, 'units') || units(target)$numerator != 'rad') {
+    stop('units(targets) is not radians')
+  }
+  if (!inherits(source, 'units') || units(source)$numerator != 'rad') {
+    stop('units(source) is not radians')
+  }
+
+  d <- source - target
+  pi_rad <- units::as_units(pi, 'rad')
+  d <- (d + pi_rad) %% (2 * pi_rad) - pi_rad
+
+  if (signed) {
+    return(d)
+  } else {
+    return(abs(d))
+  }
+}

tests/testthat/test-edge-delay.R

@@ -0,0 +1,99 @@
+# Test edge_delay
+context('test edge_delay')
+
+library(spatsoc)
+
+DT <- fread('../testdata/DT.csv')
+id <- 'ID'
+datetime <- 'datetime'
+timethreshold <- '20 minutes'
+threshold <- 50
+coords <- c('X', 'Y')
+timegroup <- 'timegroup'
+group <- 'group'
+projection <- 32736
+window <- 3
+
+
+DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+group_times(DT, datetime = datetime, threshold = timethreshold)
+direction_step(DT, id, coords, projection)
+edges <- edge_dist(DT, threshold = threshold, id = id,
+                   coords = coords, timegroup = timegroup,
+                   returnDist = TRUE, fillNA = FALSE)
+dyad_id(edges, id1 = 'ID1', id2 = 'ID2')
+fusion_id(edges, threshold = threshold)
+
+clean_DT <- copy(DT)
+clean_edges <- copy(edges)
+
+# edge_delay(DT = DT, edges = edges, id = id, window = window)
+
+test_that('edges, DT are required', {
+  expect_error(edge_delay(edges, DT = NULL))
+  expect_error(edge_delay(edges = NULL, DT))
+})
+
+test_that('arguments required, otherwise error detected', {
+  expect_error(edge_delay(edges, DT, id = NULL),
+               'id column name required')
+  expect_error(edge_delay(edges, DT, id = id, window = NULL),
+               'window is required')
+})
+
+test_that('window is numeric', {
+  expect_error(edge_delay(edges, DT, id = id, window = 'potato'),
+               'numeric')
+})
+
+test_that('column names must exist in DT', {
+  expect_error(edge_delay(edges, DT, id = 'potato'),
+               'potato field')
+
+  copy_edges <- copy(clean_edges)
+  copy_edges[, timegroup := NULL]
+  expect_error(edge_delay(copy_edges, DT, id = id, window = window),
+               'timegroup field')
+
+  copy_edges <- copy(clean_edges)
+  copy_edges[, fusionID := NULL]
+  expect_error(edge_delay(copy_edges, DT, id = id, window = window),
+               'fusionID field')
+
+  copy_edges <- copy(clean_edges)
+  copy_edges[, dyadID := NULL]
+  expect_error(edge_delay(copy_edges, DT, id = id, window = window),
+               'dyadID field')
+
+  expect_error(edge_delay(edges, DT, id = id, window = window,
+                          direction = 'potato'),
+               'potato field')
+
+  copy_DT <- copy(clean_DT)
+  copy_DT[, timegroup := NULL]
+  expect_error(edge_delay(edges, copy_DT, id = id, window = window),
+               'timegroup field')
+})
+
+test_that('no rows are added to the result edges', {
+  expect_equal(nrow(edges),
+               nrow(edge_delay(edges, DT, id = id, window = window)))
+})
+
+test_that('two columns added to the result DT', {
+  copyEdges <- copy(edges)
+
+  expect_equal(length(c('ID1', 'ID2', 'timegroup',
+                        'dyadID', 'fusionID', 'dir_corr_delay')),
+               ncol(edge_delay(edges, DT, id = id, window = window)))
+})
+
+test_that('column added to the result DT is integer', {
+  expect_type(edge_delay(edges, DT, id = id, window = window)$dir_corr_delay, 'integer')
+})
+
+test_that('returns a data.table', {
+  expect_s3_class(edge_delay(edges, DT, id = id, window = window), 'data.table')
+})
+
+# TODO: expected results tests