NOTE: Run all the code included on the “Base Code” page prior to running the following.

Setting Up

Load eigenworms and worm data of choice.

# working directory
setwd("")  # enter appropriate working directory

# load eigenworms [Broekmans et al. 2016]
ew = read.csv("eigenworms.csv", header=F, sep=",") 

# load coefficients from worm of choice
data1 = read.csv("91 Escaping Worms/50.txt", header=F, sep="")[,1:5]
data2 = read.csv("91 Escaping Worms/69.txt", header=F, sep="")[,1:5]
data = rbind(data1, data2)

Variables

These are the optimal parameters for the data produced by Broekmans et al. 2016.

E = 5  # embedding dimension
tp = 1
theta = 2  # linearity
tau = 1

# library and prediction selection
lib = c(1,600)
pred = c(601,1200)

# convert frames to seconds (i.e. indicate frames per second as fps)
fps = 20  # escaping worms

Run Functions

Run the functions created on the “Base Code” page.

# run the embedding function
matrix1 <- make_embed(data1, E, tau, tp)
matrix2 <- make_embed(data2, E, tau, tp)
matrix = rbind(matrix1, matrix2)

# run the prediction function
new_pred <- make_pred(E, matrix, theta, lib, pred)
observations_total = new_pred$obs
predictions_total = new_pred$pred

# run the EDM error function
errors <- edm_error(pred, observations_total, ew, predictions_total)

# run the cp error function
cp_mean <- cp_error(pred, data, ew)

# run the data classification function
new_class <- data_class(data, pred)
forward = new_class$forw
backward = new_class$back
turns = new_class$turn
change = new_class$chan

Error Plot

Create an RMS error vs time plot with colored classifications.

# setup the RMS error vs time plot
par(mai=c(0.9,0.95,0.15,0.2))
plot(c(1:length(errors))/fps, errors, col="white", xlab="time (s)", ylab="prediction error (rad)",
     cex.lab=1.5, cex.axis=1.5, xaxs="i", yaxs="i", ylim=c(0, 0.067))

# add colored rectangles for forwards, backwards, and turns
for (i in 1:length(forward)){
  rect(c(1:length(errors))[forward[i]]/fps, 0, c(1:length(errors))[forward[i]+2]/fps, 0.067,
       col="skyblue", border=NA)
}
for (i in 1:length(backward)){
  rect(c(1:length(errors))[backward[i]]/fps, 0, c(1:length(errors))[backward[i]+2]/fps, 0.067,
       col="yellow2", border=NA)
}
for (i in 1:length(turns)){
  rect(c(1:length(errors))[turns[i]]/fps, 0, c(1:length(errors))[turns[i]+2]/fps, 0.067,
       col="salmon", border=NA)
}

# add the error line on top
lines(c(1:length(errors))/fps, errors, lwd=3)
abline(h=cp_mean, lty="dashed")  # average cp error
box(col="black")