Engineering a High-Performance, Opinionated Pupillometry Preprocessing Framework in R

eyeris_object <- list(
  timeseries = list(
    block_1 = data.frame(...),   # Columnar timeseries with transformations
    block_2 = data.frame(...)
  ),
  events = data.frame(...),      # Event markers from experiment
  blinks = data.frame(...),      # Detected blink intervals
  params = list(...),            # Full parameter history
  latest = "pupil_raw_deblink_detransient_...",  # Current column name
  metadata = list(...)           # Sampling rate, eye tracked, etc

pupil_raw
pupil_raw_deblink
pupil_raw_deblink_detransient
pupil_raw_deblink_detransient_interpolate

deblink <- function(eyeris, extend = 50, ...) {
  pipeline_handler(
    eyeris,
    operation = deblink_pupil,
    new_suffix = "deblink",
    extend = extend,
    ...
  )
}

#' @param wp The end of passband frequency in Hz (desired lowpass cutoff).
#' @param ws The start of stopband frequency in Hz (required lowpass cutoff).
#' @param rp Required maximal ripple within passband in dB.
#' @param rs Required minimal attenuation within stopband in dB.
lpfilt <- function(eyeris, wp = 4, ws = 8, rp = 1, rs = 35, plot_freqz = TRUE, ...) {
  # design a Butterworth filter with minimum order to meet requirements
  fs_nq <- fs / 2
  foo <- gsignal::buttord(wp / fs_nq, ws / fs_nq, rp, rs)  # gsignal package
  filt <- gsignal::butter(foo, output = "Sos")
  
  # users can inspect filter characteristics before applying
  if (plot_freqz) {
    # generate Bode plot
  }

  # filter twice (forward and backward) to preserve phase information
  gsignal::filtfilt(filt, prev_pupil)
}

# load your EyeLink eye-tracking data ASC file
eyeris_data <- load_asc("path/to/your/data.asc")

# preprocess and epoch your data with eyeris glassbox functions
processed_data <- eyeris_data %>%
  glassbox() %>%
  epoch(
    events = "TRIAL_START_{trial_type}_{trial_number}",
    limits = c(-0.5, 2.0),
    label = "trial_epochs"
  )

# enable eyeris database alongside CSV files
bidsify(
  processed_data,
  bids_dir = "~/my_eyetracking_study",
  participant_id = "001",
  session_num = "01",
  task_name = "memory_task",
  csv_enabled = FALSE,   # disable generation of CSV files
  db_enabled = TRUE,     # while also creating your eyeris project database
  db_path = "study_database"  # creates study_database.eyerisdb
)

eyeris_db_collect <- function(
  bids_dir,
  db_path = "my-project",
  subjects = NULL,
  data_types = NULL,
  sessions = NULL,
  tasks = NULL,
  epoch_labels = NULL,
  eye_suffixes = NULL,
  verbose = TRUE
) {
  # connect to database
  log_info("Connecting to eyeris database...", verbose = verbose)

  # first check if duckdb is installed
  if (!check_duckdb()) {
    log_error(
      "DuckDB is required for this feature. See installation instructions above.",
      verbose = TRUE
    )
  }

  con <- tryCatch(
    {
      eyeris_db_connect(bids_dir, db_path)
    },
    error = function(e) {
      log_error("Failed to connect to database: {e$message}")
    }
  )

  # ensure disconnection on exit
  on.exit(eyeris_db_disconnect(con))

  # get available tables
  all_tables <- eyeris_db_list_tables(con)

  if (length(all_tables) == 0) {
    log_warn("No tables found in database", verbose = TRUE)
    return(list())
  }

  log_info("Found {length(all_tables)} tables in database", verbose = verbose)

  # define all possible data types
  all_data_types <- c(
    "blinks",
    "events",
    "timeseries",
    "epochs",
    "epoch_summary",
    "run_confounds",
    "confounds_events",
    "confounds_summary"
  )

  # use all data types if none specified
  if (is.null(data_types)) {
    data_types <- all_data_types
  } else {
    # validate specified data types
    invalid_types <- setdiff(data_types, all_data_types)
    if (length(invalid_types) > 0) {
      log_warn(
        "Invalid data types ignored: {paste(invalid_types, collapse = ', ')}",
        verbose = TRUE
      )
      data_types <- intersect(data_types, all_data_types)
    }
  }

  log_info(
    "Extracting data types: {paste(data_types, collapse = ', ')}",
    verbose = verbose
  )

  # extract data for each type
  result_list <- list()

  for (data_type in data_types) {
    log_info("Processing {data_type}...", verbose = verbose)

    tryCatch(
      {
        # handle epoch-specific data types
        if (
          data_type %in% c("epochs", "confounds_events", "confounds_summary")
        ) {
          if (is.null(epoch_labels)) {
            # get all available epoch labels for this data type
            type_tables <- all_tables[grepl(
              paste0("^", data_type, "_"),
              all_tables
            )]
            if (length(type_tables) > 0) {
              # extract unique epoch labels from table names
              # handles both eye-L/eye-R and eyeL/eyeR formats
              epoch_pattern <- paste0(
                data_type,
                "_[^_]+_[^_]+_[^_]+_[^_]+_(.+?)(?:_eye[LR]|_eye-[LR])?$"
              )
              extracted_labels <- unique(gsub(
                epoch_pattern,
                "\\1",
                type_tables
              ))
              # remove failed matches (when pattern doesn't match, return original string)
              extracted_labels <- extracted_labels[
                extracted_labels != type_tables &
                  !is.na(extracted_labels) &
                  extracted_labels != ""
              ]
              epoch_labels_to_use <- if (length(extracted_labels) > 0) {
                extracted_labels
              } else {
                NULL
              }
            } else {
              epoch_labels_to_use <- NULL
            }
          } else {
            epoch_labels_to_use <- epoch_labels
          }

          # extract data for each epoch label
          epoch_data_list <- list()
          if (!is.null(epoch_labels_to_use)) {
            for (epoch_label in epoch_labels_to_use) {
              epoch_data <- eyeris_db_read(
                con = con,
                data_type = data_type,
                subject = subjects,
                session = sessions,
                task = tasks,
                epoch_label = epoch_label,
                eye_suffix = eye_suffixes
              )
              if (!is.null(epoch_data) && nrow(epoch_data) > 0) {
                epoch_data_list[[epoch_label]] <- epoch_data
              }
            }
          }

          # combine all epoch data
          if (length(epoch_data_list) > 0) {
            combined_data <- do.call(rbind, epoch_data_list)
            result_list[[data_type]] <- combined_data
          }
          # check if there are any tables for this data type at all
          type_tables <- all_tables[grepl(
            paste0("^", data_type, "_"),
            all_tables
          )]
          if (length(type_tables) == 0) {
            log_warn(
              "No tables found for data type: {data_type}",
              verbose = verbose
            )
          }
        } else {
          # handle non-epoch data types
          data <- eyeris_db_read(
            con = con,
            data_type = data_type,
            subject = subjects,
            session = sessions,
            task = tasks,
            eye_suffix = eye_suffixes
          )

          if (!is.null(data) && nrow(data) > 0) {
            result_list[[data_type]] <- data
          }
          # check if there are any tables for this data type at all
          type_tables <- all_tables[grepl(
            paste0("^", data_type, "_"),
            all_tables
          )]
          if (length(type_tables) == 0) {
            log_warn(
              "No tables found for data type: {data_type}",
              verbose = verbose
            )
          }
        }
      },
      error = function(e) {
        log_warn(
          "Failed to extract {data_type}: {e$message}",
          verbose = verbose
        )
      }
    )
  }

  # filter out empty results
  result_list <- result_list[lengths(result_list) > 0]

  log_success(
    "Successfully extracted {length(result_list)} data types",
    verbose = verbose
  )
  for (dtype in names(result_list)) {
    n_rows <- nrow(result_list[[dtype]])
    n_subjects <- length(unique(result_list[[dtype]]$subject_id))
    log_info(
      "  {dtype}: {n_rows} rows across {n_subjects} subjects",
      verbose = verbose
    )
  }

  return(result_list)
}

plot_gaze_heatmap <- function(eyeris, block = 1, bins = 100, ...) {
  # Extract x,y gaze coordinates
  # Compute 2D KDE using MASS::kde2d (MASS package)
  # Overlay on screen dimensions
  # Render with viridis (or another requested) color scale
}

pipeline_handler <- function(
  eyeris,
  operation,
  new_suffix,
  ...
) {
  # 1. Validate eyeris object
  # 2. Get previous operation name
  # 3. Apply core operation (func) to each block
  # 4. Create new column with updated name
  # 5. Record parameters
  # 6. Update latest pointer
  # 7. Generate plots if plot_fun provided
  # 8. Return modified object
}

#' Winsorize pupil values
#'
#' Applies winsorization to extreme pupil values within each block.
#'
#' @param eyeris An `eyeris` object created by [load_asc()].
#' @param lower Lower quantile threshold. Default is 0.01.
#' @param upper Upper quantile threshold. Default is 0.99.
#' @param call_info A list of call information and parameters. If not provided,
#'   it will be generated from the function call.
#'
#' @return Updated `eyeris` object with new winsorized pupil column.
winsorize <- function(eyeris, lower = 0.01, upper = 0.99, call_info = NULL) {
  # create call_info if not provided
  call_info <- if (is.null(call_info)) {
    list(
      call_stack = match.call(),
      parameters = list(lower = lower, upper = upper)
    )
  } else {
    call_info
  }

  # handle binocular objects
  if (is_binocular_object(eyeris)) {
    # process left and right eyes independently
    left_result <- eyeris$left |>
      pipeline_handler(
        winsorize_pupil,
        "winsorize",
        lower = lower,
        upper = upper,
        call_info = call_info
      )
    
    right_result <- eyeris$right |>
      pipeline_handler(
        winsorize_pupil,
        "winsorize",
        lower = lower,
        upper = upper,
        call_info = call_info
      )
    
    # return combined structure
    list_out <- list(
      left = left_result,
      right = right_result,
      original_file = eyeris$original_file,
      raw_binocular_object = eyeris$raw_binocular_object
    )

    class(list_out) <- "eyeris"

    return(list_out)
  } else {
    # regular eyeris object, process normally
    eyeris |>
      pipeline_handler(
        winsorize_pupil,
        "winsorize",
        lower = lower,
        upper = upper,
        call_info = call_info
      )
  }
}

test_that("bin() validates time series monotonicity", {
  # create test data with valid monotonically increasing time series
  valid_data <- data.frame(
    time_orig = seq(0, 10000, by = 100), # original time in ms
    time_secs = seq(0, 10, by = 0.1), # converted time in seconds
    time_scaled = seq(0, 10, by = 0.1), # scaled time
    block = rep(1, 101), # block identifier
    pupil_raw = rnorm(101),
    stringsAsFactors = FALSE
  )

  # create mock eyeris object for testing
  mock_eyeris <- list(
    timeseries = list(block_1 = valid_data),
    latest = list(block_1 = "pupil_raw"),
    info = list(sample.rate = 10)
  )
  class(mock_eyeris) <- "eyeris"

  # test that valid data passes
  expect_no_error({
    result <- eyeris::bin(mock_eyeris, bins_per_second = 5, method = "mean")
  })

  # test with non-monotonic time series (should fail)
  non_monotonic_data <- data.frame(
    time_orig = c(0, 1000, 500, 2000, 3000), # decreasing at index 3
    time_secs = c(0, 1, 0.5, 2, 3), # decreasing at index 3
    time_scaled = c(0, 1, 0.5, 2, 3), # decreasing at index 3
    block = rep(1, 5),
    pupil_raw = rnorm(5),
    stringsAsFactors = FALSE
  )

  mock_eyeris_bad <- list(
    timeseries = list(block_1 = non_monotonic_data),
    latest = list(block_1 = "pupil_raw"),
    info = list(sample.rate = 10)
  )
  class(mock_eyeris_bad) <- "eyeris"

  expect_error(
    {
      eyeris::bin(mock_eyeris_bad, bins_per_second = 5, method = "mean")
    },
    "Time series is not monotonically increasing"
  )

  # test the check_time_monotonic function directly
  # test with empty time vector (should fail)
  expect_error(
    {
      eyeris:::check_time_monotonic(numeric(0))
    },
    "Time vector is NULL or empty"
  )

  # test with single time point (should fail)
  expect_error(
    {
      eyeris:::check_time_monotonic(1)
    },
    "Insufficient non-NA time points to validate monotonicity"
  )

  # test with all NA time values (should fail)
  expect_error(
    {
      eyeris:::check_time_monotonic(rep(NA_real_, 5))
    },
    "Insufficient non-NA time points to validate monotonicity"
  )

  # test with some NA values but valid monotonic sequence
  expect_no_error({
    eyeris:::check_time_monotonic(c(0, NA, 1, 2, NA, 3))
  })

  # test with valid monotonically increasing sequence
  expect_no_error({
    eyeris:::check_time_monotonic(c(0, 1, 2, 3, 4, 5))
  })

  # test with non-monotonic sequence
  expect_error(
    {
      eyeris:::check_time_monotonic(c(0, 1, 0.5, 2, 3))
    },
    "Time series is not monotonically increasing"
  )
})