NOTE: Run all the following code to create the functions prior to creating the figures.
Install and import the package “rEDM” from the CRAN Repository. CRAN - Package rEDM: Empirical Dynamic Modeling (‘EDM’)
# install this specific version using devtools
require(devtools)
install_version("rEDM", version="1.2.3", repos="http://cran.us.r-project.org")
library(rEDM)Generate an embedding matrix of data based on E, tau, and tp.
# create a function to make an embedding matrix
make_embed <- function(data, E, tau, tp){
# data rows are observations and columns are variables
data = as.matrix(data)
# calculate sizes
numSamples = dim(data)[1]
inputDim = dim(data)[2]
finalDim = inputDim*(E+1)
output = array(NA, dim=c(numSamples,finalDim))
# insert the values to be predicted in the first "inputDim" columns
output[1:numSamples,1:inputDim] <- data[1:numSamples,1:inputDim]
for (i in c(0:(E-1))){
# shift and insert the input data to create the lags
output[((i*tau)+1+tp):numSamples,((i+1)*inputDim+1):((i+2)*inputDim)] <-
data[1:(numSamples-(i*tau)-tp),1:inputDim]
}
return(output)
}Use the “s-map” method in the rEDM package to make predictions on the data.
# create a function to make predictions
make_pred <- function(E, matrix, theta, lib, pred){
# create empty vectors for observations/predictions output
observations_total = {}
predictions_total = {}
# loop through the embedding matrix and generate predictions using block_lnlp
for(target in c(1:5)){
# the optimal exclusion radius is 10
# tp is 0 because it is already accounted for in the embedding
output <- block_lnlp(matrix, target_column=target, columns=c(6:((E+1)*5)),
exclusion_radius=10, stats_only=FALSE, method="s-map",
tp=0, theta=theta, lib=lib, pred=pred)
outputs <- output$model_output[[1]]
predictions <- outputs$pred
observations <- outputs$obs
# add observations/predictions to the vectors created above
observations_total = cbind(observations_total, observations)
predictions_total = cbind(predictions_total, predictions)
}
# combine and return observations and predictions
new_list = list("obs"=observations_total, "pred"=predictions_total)
return(new_list)
}Calculate RMS error based on the observations and the predictions.
# create a function to calculate error
edm_error <- function(pred, observations_total, ew, predictions_total){
# define the length of the predictions and an empty array for error
range_loop = c(1:(pred[2]-pred[1]+1))
errors = array(NA, dim=c(length(range_loop),1))
for(t in range_loop){
# a linear combination of eigenworms will generate body angles
curr = observations_total[t,]
pos_obs = array(0, dim=c(100,1))
for(i in c(1:5)){
# multiply each coefficient by its corresponding eigenworm
pos_obs = pos_obs + ew[,i]*as.numeric(curr[i])
}
curr = predictions_total[t,]
pos_pred = array(0, dim=c(100,1))
for(i in c(1:5)){
# multiply each coefficient by its corresponding eigenworm
pos_pred = pos_pred + ew[,i]*as.numeric(curr[i])
}
# calculate the RMS error
errors[t] = (mean((pos_obs - pos_pred)^2))^(1/2)
}
return(errors)
}Calculate RMS error based on the observations and the predictions for each section of the worm.
# create a function to calculate error
edm_error_section <- function(pred, observations_total, ew, predictions_total){
# define the length of the predictions and an empty array for error
range_loop = c(1:(pred[2] - pred[1] + 1))
errors = array(NA, dim=c(length(range_loop),5))
# break the worm body angles into five equal segments
angles = c(1,20,21,40,41,60,61,80,81,100)
for(ang in 1:5){
for(t in range_loop){
# a linear combination of eigenworms will generate body angles
curr = observations_total[t,]
pos_obs = array(0, dim=c(100,1))
for(i in c(1:5)){
# multiply each coefficient by its corresponding eigenworm
pos_obs = pos_obs + ew[,i]*as.numeric(curr[i])
}
pos_obs = pos_obs[angles[(ang*2-1)]:angles[(ang*2)]]
curr = predictions_total[t,]
pos_pred = array(0, dim=c(100,1))
for(i in c(1:5)){
# multiply each coefficient by its corresponding eigenworm
pos_pred = pos_pred + ew[,i]*as.numeric(curr[i])
}
pos_pred = pos_pred[angles[(ang*2-1)]:angles[(ang*2)]]
# calculate the RMS error
errors[t,ang] = (mean((pos_obs - pos_pred)^2))^(1/2)
}
}
return(errors)
}Calculate constant predictor error based on one time step earlier.
# create a function to calculate error
cp_error <- function(pred, data, ew){
# define the length of the predictions and an empty array for error
range_loop = c(1:(pred[2] - pred[1]))
cp_errors = array(NA, dim=c(length(range_loop),1))
for(t in range_loop){
# a linear combination of eigenworms will generate body angles
curr = data[t + pred[1],]
pos_obs = array(0, dim=c(100,1))
for(i in c(1:5)){
# multiply each coefficient by its corresponding eigenworm
pos_obs = pos_obs + ew[,i]*as.numeric(curr[i])
}
curr = data[t + pred[1] - 1,]
pos_pred = array(0, dim=c(100,1))
for(i in c(1:5)){
# multiply each coefficient by its corresponding eigenworm
pos_pred = pos_pred + ew[,i]*as.numeric(curr[i])
}
# calculate the RMS error
cp_errors[t] = (mean((pos_obs - pos_pred)^2))^(1/2)
}
# find the average value of the constant predictor error time series
cp_mean = mean(cp_errors, na.rm=TRUE)
return(cp_mean)
}Use thresholds to filter data and classify sections of the time series.
# create a function to classify data
data_class <- function(data, pred){
# thresholds
delta_thresh = 3
epsilon_thresh = 0.01
turns_thresh = 10
# define eigenworms a1, a2, and a3
a1_before = data[(pred[1]:(pred[2]-2)),1]
a1 = data[((pred[1]+1):(pred[2]-1)),1]
a1_after = data[((pred[1]+2):pred[2]),1]
a2_before = data[(pred[1]:(pred[2]-2)),2]
a2 = data[((pred[1]+1):(pred[2]-1)),2]
a2_after = data[((pred[1]+2):pred[2]),2]
a3 = data[((pred[1]+1):(pred[2]-1)),3]
# estimate phase velocity in the a1,a2 plane
ang_change = atan(a2_after/a1_after) - atan(a2_before/a1_before)
change = pmin(abs(ang_change), (pi - abs(ang_change)))
thresh = which(((a1^2 + a2^2)^(1/2)) > delta_thresh)
# if phase velocity is smaller than pi/2
place_a = which(abs(ang_change) < (pi - abs(ang_change)))
new_a = c()
for (i in 1:length(place_a)){
# only use these values if they exceed the threshold
if (place_a[i] %in% thresh){new_a = append(new_a, place_a[i])}
}
change[new_a] = change[new_a]*sign(ang_change)[new_a] # return the same sign
# if phase velocity is larger than pi/2
place_b = which(abs(ang_change) > (pi - abs(ang_change)))
new_b = c()
for (i in 1:length(place_b)){
# only use these values if they exceed the threshold
if (place_b[i] %in% thresh){new_b = append(new_b, place_b[i])}
}
change[new_b] = change[new_b]*-sign(ang_change)[new_b] # reverse the sign
# find values that do not meet thresholds
not_c = c()
for (i in 1:length(change)){
if (!(i %in% append(new_a, new_b))){not_c = append(not_c, i)}
}
change[not_c] = NA # set them to NA
# classify forwards, backwards, and turns
forward = which(change < -epsilon_thresh)
backward = which(change > epsilon_thresh)
turns = which(abs(a3) > turns_thresh)
change[turns] = NA # remove phase velocity at turns
# combine and return forward, backward, turns, and change
new_list = list("forw"=forward, "back"=backward, "turn"=turns, "chan"=change)
return(new_list)
}