# Create examplary data:
ms_ex_data <- data.frame(time = seq(2, 10, by = 2), 
                     eye0_x = runif(5), 
                     eye0_y = runif(5), 
                     eye1_x = runif(5), 
                     eye1_y = runif(5))

kable(ms_ex_data, digits = 3) %>% 
  kable_styling(full_width = F, 
                position = "left") %>% 
  column_spec(1:5, width = "5em")

# Create examplary data:

data_gaze %>%  
  mutate(gaze_timestamp = gaze_timestamp - min(gaze_timestamp)) %>% 
  slice(50715:50725) %>% 
  select(gaze_timestamp, gaze_normal0_x, gaze_normal0_y, gaze_normal1_x, gaze_normal1_y) %>% 
  kable(digits = 8) %>% 
  kable_styling(full_width = F, 
                position = "left") %>% 
  column_spec(1:5, width = "5em")

data_gaze %>%  
  select(gaze_normal0_x, gaze_normal1_x) %>% 
  pivot_longer(everything()) %>% 
  mutate(name = recode(name, 
                       gaze_normal0_x = "eye0", 
                       gaze_normal1_x = "eye1")) %>% 
  group_by(name) %>% 
  summarise(n = n(), 
            NA_count = sum(is.na(value)), 
            NA_proportion = NA_count/n) %>% 
  kable() %>% 
  kable_styling()

data_stimon <- data_annot %>%  
  mutate(ttype = parse_number(label)) %>% 
  filter(ttype == 2) 

data_trials = list()

data_stimon <- data_annot %>%  
  mutate(ttype = parse_number(label)) %>% 
  filter(ttype == 2) 

data_trials = list()

for (i in 1:nrow(data_stimon)) {
  tts <- vector(length = 0)
  # epoch data:
  idx <- which((data_gaze$gaze_timestamp > (data_stimon$timestamp[i] - 0.8)) & data_gaze$gaze_timestamp < data_stimon$timestamp[i] + 2.2)
  tts <- append(tts, idx)
  data_fix <-  data_gaze[tts, ]
  data_trials[[i]] <- data_fix
}

data_clean <- bind_rows(data_trials, .id = "trial") %>% 
  mutate(trial = as.integer(as.character(trial)))

data_clean %>%  
  select(gaze_normal0_x, gaze_normal1_x) %>% 
  pivot_longer(c(gaze_normal0_x, gaze_normal1_x)) %>% 
  mutate(name = recode(name, 
                       gaze_normal0_x = "eye0", 
                       gaze_normal1_x = "eye1")) %>% 
  group_by(name) %>% 
  summarise(n = n(), 
            NA_count = sum(is.na(value)), 
            NA_proportion = NA_count/n) %>% 
  kable() %>% 
  kable_styling()

data_clean %>%  
  select(trial, gaze_normal0_x, gaze_normal1_x) %>% 
  pivot_longer(c(gaze_normal0_x, gaze_normal1_x)) %>% 
  mutate(name = recode(name, 
                       gaze_normal0_x = "eye0", 
                       gaze_normal1_x = "eye1"), 
         trial = as_factor(trial)) %>% 
  group_by(trial, name) %>% 
  summarise(n = n(), 
            NA_count = sum(is.na(value)), 
            NA_proportion = NA_count/n) %>% 
  kable() %>% 
  kable_styling() %>% 
  scroll_box(width = "100%", height = "200px")

Let’s check other approaches:

I will from now on work with the “clean” data (aka. the one cropped to the trial intervals).

data_gaze <- data_clean

Furthermore, problematic are samples whose gaze_timestamp is a duplicate, because these create ambiguous mappings.
Therefore, I suggest to kick out duplicate samples from gaze_positions.csv entirely:
[This bit could be moved up]

n_duplic_timestamps <- sum(duplicated(data_gaze$gaze_timestamp))

data_gaze <- data_gaze %>% 
  distinct_at(vars(gaze_timestamp), .keep_all = TRUE)

if (n_duplic_timestamps > 0) {
  warning(sprintf("Removing %i samples with identical timestamps. Keeping first instance of each.", n_duplic_timestamps))
}

## Warning: Removing 22 samples with identical timestamps. Keeping first instance
## of each.

Regarding point “3. For each sample only one of the eyes gets updated while info for the other one remains identical.”:

The samples in the data in gaze_positions.csv are actually merged mappings from the data in pupil_positions.csv. (ideally) two of these samples are combined by the gaze mapper and merged into one sample in the gaze positions file. This merging is done for each incoming data message from one of the pupil cams (which both run at ~200 Hz). This leads to the observed behavior, namely more than 200 samples per second (almost double) and updates only for one eye. This is useful for binocular gaze mapping, but for saccade detection, we want to avoid this. To do so, we can use the information in the column base_data in pupil_positions.csv which identifies the pupil samples from which the gaze sample has been merged.
We can split up this column into two separate ones that list the respective timestamps of the original samples:

# get separate cols for eye samples: 
data_gaze <- data_gaze %>% 
  separate(base_data, c('timestamp_eye_0','eye0' ,'timestamp_eye_1', 'eye1'), '[- ]', convert = TRUE, remove = FALSE) %>% 
  mutate(timestamp_eye_1 = if_else(eye0 == 1, timestamp_eye_0, timestamp_eye_1), 
         timestamp_eye_0 = if_else(eye0 == 1, NA_real_ , timestamp_eye_0)) %>% 
  select(-c(eye0, eye1))

data_gaze %>%  
  mutate(gaze_timestamp = gaze_timestamp - min(gaze_timestamp), 
         timestamp_eye_0 = timestamp_eye_0 - min(timestamp_eye_0, na.rm = T), 
         timestamp_eye_1 = timestamp_eye_1 - min(timestamp_eye_1, na.rm = T), ) %>% 
  slice(50715:50725) %>% 
  select(gaze_timestamp, 
         gaze_normal0_x, 
         gaze_normal0_y, 
         gaze_normal1_x, 
         gaze_normal1_y, 
         timestamp_eye_0, 
         timestamp_eye_1, 
         trial) %>% 
  kable(digits = 8) %>% 
  kable_styling(full_width = T, 
                position = "left") %>% 
  column_spec(1:5, width = "5em")

gaze_timestamp	gaze_normal0_x	gaze_normal0_y	gaze_normal1_x	gaze_normal1_y	timestamp_eye_0	timestamp_eye_1	trial
287.9364	-0.02042761	-0.00006372	0.03030742	-0.00245592	287.9404	287.9324	43
287.9385	-0.02042761	-0.00006372	0.03013676	0.00029433	287.9404	287.9366	43
287.9404	-0.02042761	-0.00006372	0.03111813	-0.00025101	287.9404	287.9404	43
287.9444	-0.02042761	-0.00006372	0.03153463	0.00074236	287.9404	287.9484	43
287.9464	-0.02042761	-0.00006372	0.03153463	0.00074236	287.9444	287.9484	43
287.9484	-0.02042761	-0.00006372	0.03153463	0.00074236	287.9483	287.9484	43
287.9504	-0.02020681	-0.00022910	0.03153463	0.00074236	287.9523	287.9484	43
287.9524	-0.02020681	-0.00022910	0.03048902	0.00031344	287.9523	287.9524	43
287.9564	-0.02212029	0.00268128	0.03048902	0.00031344	287.9604	287.9524	43
287.9584	-0.02212029	0.00268128	0.03144657	-0.00055299	287.9604	287.9565	43
287.9604	-0.02212029	0.00268128	0.02966766	0.00044098	287.9604	287.9604	43

We can now get rid of the duplicates per eye and have a closer look:

In the following, I do:

• create a new df per eye with the relevant columns
• drop the rows with duplicated rows in the actually relevant columns (as there is no new information in these)
• drop the rows with NA timestamps (this should be only one per df, after the distinct command)
• sort the rows ascending in time
• calculate the distance between the samples (diff())
• add a column with the (5 point) running average (note that k in rollmean() here needs to be 1 step smaller as x in the wanted x-point average to account for the fact that it is run on differences)

timings_eye0 <- data_gaze %>% 
  select(timestamp_eye_0, trial, gaze_timestamp, gaze_normal0_x, gaze_normal0_y) %>% 
  distinct_at(vars(timestamp_eye_0, gaze_normal0_x, gaze_normal0_y), .keep_all = TRUE) %>% 
  drop_na(timestamp_eye_0, gaze_normal0_x, gaze_normal0_y) %>% 
  group_by(trial) %>% 
  arrange(timestamp_eye_0) %>% 
  mutate(sample_dist = c(NA_real_, diff(timestamp_eye_0)), 
         sample_dist_smooth = rollmean(sample_dist, 4, fill = NA))  # k = 

timings_eye1 <- data_gaze %>% 
  select(timestamp_eye_1, trial, gaze_timestamp, gaze_normal1_x, gaze_normal1_y) %>% 
  distinct_at(vars(timestamp_eye_1, gaze_normal1_x, gaze_normal1_y), .keep_all = TRUE) %>% 
  drop_na(timestamp_eye_1, gaze_normal1_x, gaze_normal1_y) %>% 
  group_by(trial) %>% 
  arrange(timestamp_eye_1) %>% 
  mutate(sample_dist = c(NA_real_, diff(timestamp_eye_1)), 
         sample_dist_smooth = rollmean(sample_dist, 4, fill = NA))

We dropped all relevant NAs by now. Is this problemeatic, i.e. did we kick out a lot of samples?
According to what we did above, we would kill samples where the relevant columns for both eyes would be NA. So these are no bad loss. (And actually, there shouldn’t be any: sum(is.na(data_gaze$timestamp_eye_0) & is.na(data_gaze$timestamp_eye_1)) = 0.)
All other samples that we loose must be such which are duplicates for both eyes, or NA for one eye and a duplicate for the other. Also these are no dramatic loss. We can check for their number/proportion:

remaining_samples <- unique(c(timings_eye0$gaze_timestamp, timings_eye1$gaze_timestamp))
n_dropped_samples <- length(setdiff(data_gaze$gaze_timestamp, remaining_samples))
prop_dropped_samples <- n_dropped_samples/nrow(data_gaze)

Number of discarded samples: 19
Which accounts for 0.022% of all samples.

So this is not so problematic.

Does it resolve the problem of the uneven SF?

mean_sample_dist_eye0 <- mean(timings_eye0$sample_dist, na.rm = T)
mean_sample_dist_eye1 <- mean(timings_eye1$sample_dist, na.rm = T)

# SF right eye:
mean_SF_eye0 <- 1/mean_sample_dist_eye0
# SF left eye:
mean_SF_eye1 <- 1/mean_sample_dist_eye1

Avg. SF right eye: 199.2590682
Avg. SF left eye: 199.1858299

These means are very close to the 200 Hz that we expected. However, the spacing unfortunately is not as uniform as one might hope:

timings_eye0 %>% 
  ggplot(aes(x = sample_dist)) + 
  geom_histogram(bins = 300) + 
  #xlim(c(0, 0.01)) + 
  labs(x = "Avg. time between samples (s)", 
       title = "Right eye") + 
  theme_minimal()

It seems like the inverse of the mean time betwwen two samples is not the best representative of the SF actually. For most of the samples the SF is actually a bit higher than expected (~ 250 Hz), but for around 25% of the samples, it looks like a camera frame/sample was lost and therefore the SF is about half (~125 Hz). Also, the automatic fit of the x axis shows that there seems to be at least 1 outlier with a delay of around 4x the average.

summary(timings_eye0$sample_dist)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## 0.00083 0.00396 0.00403 0.00502 0.00699 0.01974      72

prop_outliers_larger_10ms <- sum(timings_eye0$sample_dist > 0.01, na.rm = T)/nrow(timings_eye1)
idx_outliers_larger_10ms <- which(timings_eye0$sample_dist > 0.01)

plot(1:89, as_vector(timings_eye0[(idx_outliers_larger_10ms[15]-44):(idx_outliers_larger_10ms[15]+44),6]))

plot(1:89, as_vector(timings_eye0[(idx_outliers_larger_10ms[15]-44):(idx_outliers_larger_10ms[15]+44),4]))
lines(c(45, 45), c(-100, 100))
lines(c(25, 25), c(-100, 100))

So roughly 0.9% of samples have a delay larger than 10ms (ideal would be 5ms). Knowing that MSTB applies a 5-point running average to smoothen the data, I applied the same procedure above to the time delays between samples. Here’s the result:

timings_eye0 %>% 
  ggplot(aes(x = sample_dist_smooth)) + 
  geom_histogram(bins = 3000) + 
  #xlim(c(0, 0.01)) + 
  labs(x = "Avg. time between samples (s)", 
       title = "Right eye") + 
  theme_minimal()

summary(timings_eye0$sample_dist_smooth)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## 0.00262 0.00499 0.00500 0.00502 0.00502 0.00894     288

prop_outliers_larger_10ms <- sum(timings_eye0$sample_dist_smooth > 0.01, na.rm = T)/nrow(timings_eye1)

The smoothing removes the outlier problem (at least for this data set) and most widely gets us rid of the 2 peaks in the distribution. But note that rollmean() (standard 5 point avg over the sample distances) does not give the same result as the “running average” implemented in MSTB (see vecvel.R). vecvel.R with smoothing TYPE==2 is not a plain implementation of a 5p moving average over velocities:

$\text{v}\text{→}{}_{\mathrm{n}}\text{=}\text{x}\text{→}{}_{\mathrm{n+2}}\text{+}\text{x}\text{→}{}_{\mathrm{n+1}}\text{−}\text{x}\text{→}{}_{\mathrm{n-1}}\text{−}\text{x}\text{→}{}_{\mathrm{n-2}}\text{6}\text{Δ}\text{t}\text{,}$

(from: Engbert & Kliegl, 2003)

Here’s a first draft for a pipeable version of vecvel.R. For now it only calculates the smoothened version of the timestamps. (It’s important that it is pipeable, so that it can take into account the group nature of the data frame.)

vectime_custom <- function(.data, ...) {

  
  if (dplyr::is_grouped_df(.data)) {
    return(dplyr::do(.data, vectime_custom(., ...)))
  }

  
  args <- list(...)
  
  d <- dim(.data)
  N <- d[1]
  v <- matrix(rep(0,N),ncol=1)
  
  if (args$type == 2) {
    slice_ <- list()
    for (i in c(1,2,4,5)) {
      end <- N-(5-i)
      slice_[[i]] <- slice(.data, i:end)
    }
    
    v[3:(N-2), ] <- 1/6* (slice_[[5]]$timestamp_eye_0 
                          + slice_[[4]]$timestamp_eye_0 
                          - slice_[[2]]$timestamp_eye_0
                          - slice_[[1]]$timestamp_eye_0)
    
    
    #v[3:(N-2),] <- 1/6*(x[5:N,] + x[4:(N-1),] - x[2:(N-3),] - x[1:(N-4),])
    
  
    v[2,] = 1/2*(.data$timestamp_eye_0[3] - .data$timestamp_eye_0[1])
    v[(N-1),] = 1/2*(.data$timestamp_eye_0[N] - .data$timestamp_eye_0[(N-2)])  
  } else {
    
    if ("winsize" %in% names(args)) {
      if (!(args$winsize %% 2)) {
        stop("Window size must be an odd number.")
      }
      winsize_ <- args$winsize
    } else {
      winsize_ <- 3
    }
    
    slice_ <- list()
    slice_[['lower']] <- slice(.data, 1:(N-(winsize_ - 1)))
    slice_[['upper']] <- slice(.data, winsize_:(N))
    
   #v <- rollmean(c(NA_real_, diff(.data$timestamp_eye_0)), winsize_ - 1, na.pad = T)   # this is for demo purposes; yields the same result
    
   v[ceiling(winsize_/2):(N-floor(winsize_/2)), ] <- 1/(winsize_ - 1) * (slice_[['upper']]$timestamp_eye_0 - slice_[['lower']]$timestamp_eye_0)
    
  }
  
  
  .data$timestamp_eye_0 <- v
  return(.data)
}

Here’s the original version, to compare:

vecvel <- function(x,SAMPLING=500,TYPE=2) {
  d <- dim(x)
  N <- d[1]
  v <- matrix(rep(0,2*N),ncol=2)
  
  if ( TYPE==2 ) {
    v[3:(N-2),] <- SAMPLING/6*(x[5:N,] + x[4:(N-1),] - x[2:(N-3),] - x[1:(N-4),])
    v[2,] = SAMPLING/2*(x[3,] - x[1,])
    v[(N-1),] = SAMPLING/2*(x[N,] - x[(N-2),])   
  }  else  {
    v[2:(N-1),] <- SAMPLING/2*(x[3:N,] - x[1:(N-2),])
  }
  return(v)
}

here is the result of applying vectime_custom.R to our gaze data frame:

timings_eye0 %>% 
  select(timestamp_eye_0, trial) %>% 
  vectime_custom(., type = 2) %>% 
  ggplot(aes(x = timestamp_eye_0)) +
  geom_histogram(bins = 3000) +
  labs(x = "Avg. time between samples (s)",
       title = "Right eye") +
  theme_minimal()

The two peaks are now closer together but still very prominent.

My custom adaptation of vecvel allows for a similar usage of the typeargument as the original. So if we use type == 1 (actually just any number but 2), we get the result for running a plain 3-point average:

timings_eye0 %>% 
  select(timestamp_eye_0, trial) %>% 
  vectime_custom(., type = 1) %>% 
  ggplot(aes(x = timestamp_eye_0)) +
  geom_histogram(bins = 3000) +
  labs(x = "Avg. time between samples (s)",
       title = "Right eye") +
  theme_minimal()

This does not make a big difference. Interesting though is the difference to the plain 5-point average that I showed above. Bug in the code?
I don’t think so. My vectime_custom now also allows to adaptively set a window size with the parameter winsize (must be an odd number). So here’s the result for a 5-point window:

timings_eye0 %>% 
  select(timestamp_eye_0, trial) %>% 
  vectime_custom(., type = 1, winsize = 5) %>% 
  ggplot(aes(x = timestamp_eye_0)) +
  geom_histogram(bins = 3000) +
  labs(x = "Avg. time between samples (s)",
       title = "Right eye") +
  theme_minimal()

That’s the same result as shown above with rollmean() over diff(timestamp). So the difference between vecvel.R’s 5-point running average and the classical 5-point average seems to play a role here.

Preliminary summary:
The choice of the windowing/averaging parameter seems to play an essential role. I see three options:

1. Use the classical vecvel.R algorithm and pretend the SF was constantly 200 Hz, which is pretty much in the center of the two peaks and they are not too wide apart.
2. Use “normal” 5-point averaging (so set TYPEto 1 in vecvel.R). This leads to a centered distribution in the time domain (which does not necessarily say something about what it does to the velocities) and therefore hurt the assumption of a stable SF less. This would require tweaking the implementations of the MSTB just a bit because only a 3-point averaging is implemented (for TYPE != 2).
3. Implement a custom version of vecvel.R and all other scripts in the MSTB that expect a stable SF so that they respect the timestamps and calculate the velocity in respect of the actual temporal distance between two samples. This would probably be the cleanest option but requires digging into MSTB further, finding out which scripts depend on the SF and fix them accordingly.