In JimMcL/trajr: Animal Trajectory Analysis

knitr::opts_chunk$set(echo = TRUE)
library(trajr)
library(microbenchmark)

Some miscellaneous examples

Mean vector characterising changes in direction:

# Using vector arithmetic. This is roughly twice as fast as using calculus
angularMeanVectorArith <- function(angles) mean(complex(modulus = 1, argument = angles))
# Using calculus
angularMeanCalculus <- function(angles) {
  phi <- atan2(sum(sin(angles)), sum(cos(angles)))
  r <- sqrt(sum(cos(angles)) ^ 2 + sum(sin(angles)) ^ 2) / length(angles)
  complex(modulus = r, argument = phi)
}

angles <- runif(10000, min = -pi, max = pi)

microbenchmark(angularMeanVectorArith(angles), angularMeanCalculus(angles))

Emax

Empirical testing of Emax with generated trajectories.

set.seed(1)
n <- 5
angularErrorSd <- runif(n, 0, 2)

# Assess Emax for some different types of trajectories, all with the same angular error sd
emaxr <- sapply(1:n, function(i) { TrajEmax(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i])) })
emaxd <- sapply(1:n, function(i) { TrajEmax(TrajGenerate(500, random = FALSE, angularErrorSd = angularErrorSd[i]), compass.direction = 0) })
emaxSm <- sapply(1:n, function(i) { TrajEmax(TrajRediscretize(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i]), R = 1)) })
emaxBig <- sapply(1:n, function(i) { TrajEmax(TrajRediscretize(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i]), R = 10)) })
emaxBig[emaxBig < 0] <- NA # can't take log of negative numbers

par(mar = c(4, 4, 1, 1))
plot(rep(angularErrorSd, 4), c(emaxr, emaxd, emaxSm, emaxBig),
     log = 'xy', pch = '.', 
     col = c(rep('red', n), rep('blue', n), rep('green', n), rep('black', n)),
     xlab = expression(sigma["angular error"]), ylab = expression("E"["max"]))
       labels <- c("Random walk", 
                   "Directed walk", 
                   "Random walk, step length 1", 
                   "Random walk, step length 10")
legend("bottomleft", labels, pch = 16, col = c("red", "blue", "green", "black"), inset = 0.01)

Sinuosity

Empirical testing of Sinuosity with generated trajectories.

set.seed(1)
n <- 1000
angularErrorSd <- runif(n, 0, 2)

# Assess sinuosity for some different types of trajectories, all with the same angular error sd
sinur <- sapply(1:n, function(i) { 
  TrajSinuosity(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i])) 
})
sinud <- sapply(1:n, function(i) { 
  TrajSinuosity(TrajGenerate(500, random = FALSE, angularErrorSd = angularErrorSd[i]),
                compass.direction = 0) 
})
sinuSm <- sapply(1:n, function(i) { 
  TrajSinuosity(TrajRediscretize(TrajGenerate(500, random = TRUE, 
                                              angularErrorSd = angularErrorSd[i]), 
                                 R = 1)) 
})
sinuBig <- sapply(1:n, function(i) { 
  TrajSinuosity(TrajRediscretize(TrajGenerate(500, random = TRUE, 
                                              angularErrorSd = angularErrorSd[i]), 
                                 R = 10)) 
})

par(mar = c(4, 4, 1, 1))
plot(rep(angularErrorSd, 4), c(sinur, sinud, sinuSm, sinuBig),
     pch = '.', 
     col = c(rep('red', n), rep('blue', n), rep('green', n), rep('black', n)),
     xlab = expression(sigma["angular error"]), ylab = "Sinuosity")
legend("bottomright", c("Random walk", "Directed walk", "Random walk, step length 1", "Random walk, step length 10"), pch = 16, col = c("red", "blue", "green", "black"), inset = 0.01)

Relationship between Sinuosity and Sinuosity 2

set.seed(1)
n <- 1000

angularErrorSd <- runif(n, 0, pi)

# Assess sinuosities for some different types of trajectories, all with the same angular error sd
sins <- sapply(1:n, function(i) { 
  c(TrajSinuosity(TrajGenerate(500, angularErrorSd = angularErrorSd[i])),
    TrajSinuosity2(TrajGenerate(500, angularErrorSd = angularErrorSd[i])) )
})
plot(sins[2,] ~ sins[1,], pch = 16, cex = .2, xlab = 'TrajSinuosity', ylab = 'TrajSinuosity2')
abline(a = 0, b = 1, col = "red")

Relationship between Emax and Sinuosity

set.seed(1)
n <- 1000
angularErrorSd <- runif(n, 0, 2)

xy <- sapply(1:n, function(i) {
  trj <- TrajGenerate(500, angularErrorSd = angularErrorSd[i])
  c(Emax = TrajEmax(trj), sinuosity = TrajSinuosity(trj))
})
xy <- xy[, which(xy[1,] > 0)]
plot(log(xy[1,]), xy[2,], pch = '.', xlab = expression(log(E[max])), ylab = "Sinuosity")

Fractal dimension

set.seed(1)
n <- 20
angularErrorSd <- runif(n, 0, 2)
stepCount <- round(runif(n, 20, 100))

# Use the same step sizes for all trajectories
#trj <- TrajGenerate((max(stepCount) + min(stepCount)) / 2)
stepSizes <- exp(log(10) * seq(log10(.1), log10(7), length.out = 20))


dim <- sapply(1:n, function(i) { 
  trj <- TrajGenerate(stepCount[i], random = TRUE, angularErrorSd = angularErrorSd[i])
  c(TrajFractalDimension(trj, stepSizes), 
    TrajFractalDimension(trj, stepSizes, FALSE, FALSE))
})
plot(angularErrorSd, dim[1,], 
     pch = 16, 
     col = rgb((stepCount - min(stepCount)) / (max(stepCount) - min(stepCount)), 0.2, 0.1),
     ylim = range(dim),
     xlab = expression(sigma["angular error"]), ylab = "Fractal dimension")
col <- rgb(0.2, (stepCount - min(stepCount)) / (max(stepCount) - min(stepCount)), .01)
points(angularErrorSd, dim[2,], pch = 3, col = col)

Plot of path length to divider size for a range of divider sizes, Levy flight and correlated random walk.

set.seed(1)

.plotFractalValues <- function(trj, ...) {
  muL <- mean(TrajStepLengths(trj))
  stepSizes <- TrajLogSequence(0.5 * muL, 5 * muL, 20)
  v <- TrajFractalDimensionValues(trj, stepSizes)
  d <- TrajFractalDimension(trj, stepSizes)

  l <- stats::lm(pathlength ~ stepsize, data = v)

  plot(v, log = "xy", pch = 16, ...)
  #abline(l, untf = TRUE)
  mtext(sprintf("Fractal dimension = %g", d), 1, -1.1, adj = 0.05)
}

lf <- TrajGenerate(linearErrorDist = rcauchy)
crw <- TrajGenerate()

par(mfrow = c(2, 1))
.plotFractalValues(lf, main = "Levy flight")
.plotFractalValues(crw, main = "Correlated random walk", col = "red")

Example theoretical anlysis

Expect travel distance

Following is an implementation of a simulation study from (Cheung at el., 2007). It plots the expected endpoints of multiple random trajectories after different numbers of steps, illustrating the difference in expected distance travelled between directed and random walks. It demonstrates that, for travel without a compass (i.e. an external means to determine direction), the distance travelled along an intended direction does not increase in proportion to the number of steps, which implies an upper limit to the distance that can be travelled without a compass.

set.seed(1)
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
ntraj <- 1000 # Cheung at el. used 20000, but it's slow

# Define a function to simulate journeys of different lengths
simulateJourneys <- function(random, nTrajectories = ntraj) {

  generate <- function(i) {
    # These parameters replicate the simulations in (Cheung et al., 2007)
    TrajGenerate(n = 500, random = random, stepLength = 2, 
                 linearErrorSd = 0.2, angularErrorSd = 0.5)
  }

  trjs <- lapply(1:nTrajectories, generate)

  # Use same plot axes as in the original article
  if (random) {
    xlims <- list(c(-20, 20), c(-20, 20), c(-20, 20), c(-40, 40), c(-170, 170), c(-400, 400))
    ylims <- xlims
  } else {
    xlims <- list(c(-20, 20), c(-20, 20), c(-20, 20), c(20, 50), c(140, 200), c(810, 950))
    ylims <- list(c(-20, 20), c(-20, 20), c(-20, 20), c(-20, 20), c(-30, 30), c(-70, 70))
  }

  # 6 plots in 2 rows
  par(mfrow = c(2, 3))
  steps <- c(1, 2, 5, 20, 100, 500)
  for(i in 1:length(steps)) {
    step <- steps[i]
    title <- sprintf("%d %s", step, ifelse(step == 1, "step", "steps"))
    # Plot the end-points of each journey after "step" steps
    endPoints <- t(sapply(trjs, function(t) t[step + 1, c('x', 'y')]))
    plot(endPoints, pch = '.', asp = 1, main = title, xlim = xlims[[i]], ylim = ylims[[i]])
  }
}

# Call the function
simulateJourneys(random = FALSE)

set.seed(1)
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
simulateJourneys(random = TRUE)

Random or directed?

Benhamou (2006) defines a method, called backward pattern analysis, to determine whether a trajectory is a random or directed walk. Here we test it on a number of randomly generated trajectories, and report the failure rate.

# Returns TRUE if trj is a random trajectory according to backward pattern analysis (Benhamou, 2006)
isRandom <- function(trj) {
  # For each pair of points along the trajectory,
  # calculate distance between the points
  D <- sapply(1:(nrow(trj) - 1), function(j) TrajDistance(trj, j, nrow(trj)))
  # calculate path length between points
  L <- sapply(1:(nrow(trj) - 1), function(j) TrajLength(trj, j, nrow(trj)))
  # Calculate linear regression with 0 intercept: D vs L
  dir <- lm(D ~ L - 1)
  # D vs sqrt(L)
  sL <- sqrt(L)
  rnd <- lm(D ~ sL - 1)

  # Pick the linear regression with the smallest sum of squared residuals
  sum(rnd$residuals ^ 2) < sum(dir$residuals ^ 2)
}

# Generate some random trajectories, randomly choosing 
# between random search path and directed walks
set.seed(1)
search <- sample(c(TRUE,FALSE), 1000, TRUE) # search/directed flag
results <- sapply(1:length(search), function(i) {
  c(search[i], isRandom(TrajGenerate(100, random = search[i]))) 
  })
# Flags whether the algorithm correctly identified the trajectory type
correct <- results[1,] == results[2,]
t <- table(data.frame(correct, results[1,]))

cat(sprintf("Backward pattern analysis misclassified %g%% of random walks, and %g%% of directed walks\n", 
            100 * t[1,2] / sum(t), 100 * t[1,1] / sum(t)))

An actual fractal - the Hilbert curve

Following is a plot of the values used to calculate the fractal dimension of an approximation of an actual fractal curve. The points should lie on a straight line if the curve is a fractal. Also shows an inflexion point (i.e. horizontal line, then sloped line) since the value of n used to create the curve defines the scale down to which it is fractal.

# Adapted from http://www.soc.napier.ac.uk/~andrew/hilbert.html
Hilbert <- function(x0, y0, xis, xjs, yis, yjs, n) {
  # Preallocate a matrix
  result <- matrix(NA, 4 ^ n, 2)
  row <- 1

  .h <- function(x0, y0, xis, xjs, yis, yjs, n) {
    # x0 and y0 are the coordinates of the bottom left corner
    # xis & xjs are the i & j components of the unit x vector this frame
    # similarly yis and yjs
    if (n <= 0) {
      result[row, 1] <<- x0+(xis+yis)/2
      result[row, 2] <<- y0+(xjs+yjs)/2
      row <<- row + 1
    } else {
      .h(x0, y0, yis/2, yjs/2, xis/2, xjs/2, n-1)
      .h(x0+xis/2, y0+xjs/2 ,xis/2, xjs/2, yis/2, yjs/2, n-1)
      .h(x0+xis/2+yis/2, y0+xjs/2+yjs/2, xis/2, xjs/2, yis/2, yjs/2,n-1)
      .h(x0+xis/2+yis, y0+xjs/2+yjs, -yis/2,-yjs/2, -xis/2, -xjs/2,n-1)
    }
  }

  .h(x0, y0, xis, xjs, yis, yjs, n)
  result
}

trj <- TrajFromCoords(as.data.frame(Hilbert(0, 0, 300, 0, 0, 300, 8)))
plot(TrajFractalDimensionValues(trj, TrajLogSequence(.1, 400, 20)), log = 'xy', pch = 16)
trj <- TrajFromCoords(as.data.frame(Hilbert(0, 0, 300, 0, 0, 300, 4)))
points(TrajFractalDimensionValues(trj, TrajLogSequence(.1, 300, 40)), pch = 16, col = "red")