R/plot.eyeris.R

In eyeris: Flexible, Extensible, & Reproducible Pupillometry Preprocessing

#' Plot pre-processed pupil data from `eyeris`
#'
#' S3 plotting method for objects of class `eyeris`. Plots a single-panel
#' timeseries for a subset of the pupil time series at each preprocessing step.
#' The intended use of this function is to provide a simple method for
#' qualitatively assessing the consequences of the preprocessing recipe and
#' parameters on the raw pupillary signal.
#'
#' @param x An object of class `eyeris` derived from [eyeris::load_asc()]
#' @param ... Additional arguments to be passed to `plot`
#' @param steps Which steps to plot; defaults to `all` (i.e., plot all steps).
#' Otherwise, pass in a vector containing the index of the step(s) you want to
#' plot, with index `1` being the original raw pupil time series
#' @param preview_n Number of random example "epochs" to generate for
#' previewing the effect of each preprocessing step on the pupil time series
#' @param preview_duration Time in seconds of each randomly selected preview
#' @param preview_window The start and stop raw timestamps used to subset the
#' preprocessed data from each step of the `eyeris` workflow for visualization
#' Defaults to NULL, meaning random epochs as defined by `preview_n` and
#' `preview_duration` will be plotted. To override the random epochs, set
#' `preview_window` here to a vector with relative start and stop times (in
#' seconds), for example -- `c(5,6)` -- to indicate the raw data from 5-6 secs
#' on data that were recorded at 1000 Hz). Note, the start/stop time values
#' indicated here are in seconds because `eyeris` automatically computes the
#' indices for the supplied range of seconds using the `$info$sample.rate`
#' metadata in the `eyeris` S3 class object
#' @param seed Random seed for current plotting session. Leave NULL to select
#' `preview_n` number of random preview "epochs" (of `preview_duration`) each
#' time. Otherwise, choose any seed-integer as you would normally select for
#' [base::set.seed()], and you will be able to continue re-plotting the same
#' random example pupil epochs each time -- which is helpful when adjusting
#' parameters within and across `eyeris` workflow steps
#' @param block For multi-block recordings, specifies which block to plot.
#' Defaults to 1. When a single `.asc` data file contains multiple
#' recording blocks, this parameter determines which block's time series to
#' visualize. Must be a positive integer not exceeding the total number of
#' blocks in the recording
#' @param plot_distributions Logical flag to indicate whether to plot both
#' diagnostic pupil time series *and* accompanying histograms of the pupil
#' samples at each processing step. Defaults to `FALSE`
#' @param suppress_prompt Logical flag to disable interactive confirmation
#' prompts during plotting. Defaults to `TRUE`, which avoids hanging behavior in
#' non-interactive or automated contexts (e.g., RMarkdown, scripts)
#' Set to `FALSE` only when running inside `glassbox()` with
#' `interactive_preview = TRUE`, where prompting after each step is desired, as
#' well as in the generation of interactive HTML reports with [eyeris::bidsify]
#' @param verbose A logical flag to indicate whether to print status messages to
#' the console. Defaults to `TRUE`. Set to `FALSE` to suppress messages about
#' the current processing step and run silently
#' @param add_progressive_summary Logical flag to indicate whether to add a
#' progressive summary plot after plotting. Defaults to `FALSE`. Set to `TRUE`
#' to enable the progressive summary plot (useful for interactive exploration).
#' Set to `FALSE` to disable the progressive summary plot (useful in automated
#' contexts like bidsify reports)
#' @param eye For binocular data, specifies which eye to plot: "left", "right",
#' or "both". Defaults to "left". For "both", currently plots left eye data
#' (use eye="right" for right eye data)
#' @param num_previews **(Deprecated)** Use `preview_n` instead
#'
#' @return No return value; iteratively plots a subset of the pupil time series
#' from each preprocessing step run
#'
#' @seealso [lifecycle::deprecate_warn()]
#'
#' @examples
#' # first, generate the preprocessed pupil data
#' my_eyeris_data <- system.file("extdata", "memory.asc", package = "eyeris") |>
#'   eyeris::load_asc() |>
#'   eyeris::deblink(extend = 50) |>
#'   eyeris::detransient() |>
#'   eyeris::interpolate() |>
#'   eyeris::lpfilt(plot_freqz = TRUE) |>
#'   eyeris::zscore()
#'
#' # controlling the time series range (i.e., preview window) in your plots:
#'
#' ## example 1: using the default 10000 to 20000 ms time subset
#' plot(my_eyeris_data, seed = 0, add_progressive_summary = TRUE)
#'
#' ## example 2: using a custom time subset (i.e., 1 to 500 ms)
#' plot(
#'   my_eyeris_data,
#'   preview_window = c(0.01, 0.5),
#'   seed = 0,
#'   add_progressive_summary = TRUE
#' )
#'
#' # controlling which block of data you would like to plot:
#'
#' ## example 1: plots first block (default)
#' plot(my_eyeris_data, seed = 0)
#'
#' ## example 2: plots a specific block
#' plot(my_eyeris_data, block = 1, seed = 0)
#'
#' ## example 3: plots a specific block along with a custom preview window
#' ##   (i.e., 1000 to 2000 ms)
#' plot(
#'   my_eyeris_data,
#'   block = 1,
#'   preview_window = c(1, 2),
#'   seed = 0
#' )
#'
#' @rdname plot.eyeris
#'
#' @export
plot.eyeris <- function(
  x,
  ...,
  steps = NULL,
  preview_n = NULL,
  preview_duration = NULL,
  preview_window = NULL,
  seed = NULL,
  block = 1,
  plot_distributions = FALSE,
  suppress_prompt = TRUE,
  verbose = TRUE,
  add_progressive_summary = FALSE,
  eye = c("left", "right", "both"),
  num_previews = deprecated()
) {
  # handle deprecated parameters
  if (is_present(num_previews)) {
    deprecate_warn("1.2.0", "plot(num_previews)", "plot(preview_n)")
    preview_n <- num_previews
  }

  eye_suffix <- NULL
  # handle binocular eyeris objects
  eye <- match.arg(eye)
  if (is_binocular_object(x)) {
    if (eye == "left") {
      x <- x$left
      eye_suffix <- "eye-L"
      log_info("Plotting left eye data", verbose = verbose)
    } else if (eye == "right") {
      x <- x$right
      eye_suffix <- "eye-R"
      log_info("Plotting right eye data", verbose = verbose)
    } else if (eye == "both") {
      x <- x$left
      log_info(
        "Plotting left eye data (use eye='right' for right eye)",
        verbose = verbose
      )
    }
  }

  # safely handle user's current options
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  # tests
  tryCatch(
    {
      check_data(x, "plot")
    },
    error = function(e) {
      error_handler(e, "input_data_type_error")
    }
  )

  tryCatch(
    {
      check_pupil_cols(x, "plot")
    },
    error = function(e) {
      error_handler(e, "missing_pupil_raw_error")
    }
  )

  params <- list(...)

  only_liner_trend <- if ("only_linear_trend" %in% names(params)) {
    params$only_linear_trend <- params$only_linear_trend
  } else {
    params$only_linear_trend <- FALSE
  }

  non_plot_params <- c(
    "preview_window",
    "seed",
    "steps",
    "num_previews",
    "preview_n",
    "preview_duration",
    "block",
    "suppress_prompt",
    "plot_distributions",
    "only_linear_trend",
    "next_step",
    "add_progressive_summary",
    "eye"
  )

  plot_params <- params[!(names(params) %in% non_plot_params)]

  # set param defaults outside of function declaration
  if (!is.null(preview_window)) {
    if (!is.null(preview_n) || !is.null(preview_duration)) {
      log_warn(
        "preview_n and/or preview_duration will be ignored,",
        "since preview_window was specified here."
      )
    }
  }

  if (is.null(steps)) {
    steps <- "all"
  }

  if (is.null(preview_n)) {
    preview_n <- 3
  }

  if (is.null(preview_duration)) {
    preview_duration <- 5 # seconds
  }

  hz <- if (!is.na(x$decimated.sample.rate)) {
    x$decimated.sample.rate
  } else {
    x$info$sample.rate
  }

  # handle random seed for this plotting session
  if (is.null(seed)) {
    seed <- sample.int(.Machine$integer.max, 1)
  }

  # blocks handler
  if (is.list(x$timeseries) && !is.data.frame(x$timeseries)) {
    available_blocks <- get_block_numbers(x)
    available_blocks_sorted <- sort(as.numeric(available_blocks), na.last = NA)

    if (block %in% available_blocks) {
      pupil_data <- x$timeseries[[paste0("block_", block)]]
      log_info(
        "Plotting block {block} with sampling rate {hz} Hz from possible blocks: {toString(available_blocks_sorted)}",
        verbose = verbose
      )
    } else {
      log_error(
        "Block {block} does not exist. Available blocks: {toString(available_blocks_sorted)}"
      )
    }
  } else {
    pupil_data <- x$timeseries$block_1
  }

  pupil_steps <- grep("^pupil_", names(pupil_data), value = TRUE)

  # when a downsample/bin step has run, the stored timeseries only retains the
  # decimated samples. plot earlier (pre-decimation) pipeline steps from the
  # preserved full-resolution data so they are not misleadingly shown at the
  # decimated sampling rate (issue #294)
  full_pupil_data <- get_pre_decimation_block(x, block)
  has_decimation <- !is.null(full_pupil_data)
  full_hz <- x$info$sample.rate

  colorpal <- eyeris_color_palette()
  colors <- c("black", colorpal)

  transparent_colors <- sapply(colors, function(x) {
    grDevices::adjustcolor(x, alpha.f = 0.5)
  })

  if (length(steps) == 1) {
    if (steps[1] == "all") {
      pupil_steps <- pupil_steps
      colors <- colors
    } else {
      pupil_steps <- pupil_steps[steps]
      colors <- colors[steps]
    }
  } else if (length(steps) > 1 && !is.null(preview_window)) {
    pupil_steps <- pupil_steps[steps]
    colors <- colors[steps]
  } else {
    pupil_steps <- pupil_steps
    colors <- colors
  }

  if (is.null(preview_window)) {
    withr::with_seed(seed, {
      random_epochs <- draw_random_epochs(
        pupil_data,
        preview_n,
        preview_duration,
        hz
      )
    })

    # build matching full-resolution epochs (over the same time windows) so
    # that pre-decimation steps can be plotted at their original sampling rate
    random_epochs_full <- NULL
    if (has_decimation) {
      epoch_starts <- vapply(
        random_epochs,
        function(e) min(e$time_secs, na.rm = TRUE),
        numeric(1)
      )
      random_epochs_full <- lapply(epoch_starts, function(s0) {
        slice_epoch_window(full_pupil_data, s0, preview_duration)
      })
    }

    par(mfrow = c(1, preview_n), oma = c(0, 0, 3, 0))
    detrend_plotted <- FALSE
    for (i in seq_along(pupil_steps)) {
      use_full <- has_decimation && !is_decimated_col(pupil_steps[i])
      for (n in 1:preview_n) {
        epoch_n <- if (use_full) random_epochs_full[[n]] else random_epochs[[n]]
        st <- min(epoch_n$time_orig, na.rm = TRUE)
        et <- max(epoch_n$time_orig, na.rm = TRUE)
        title <- paste0("\n[", st, " - ", et, "]")
        header <- paste0(
          gsub("_", " > ", gsub("pupil_", "", pupil_steps[i])),
          if (is.list(x$timeseries) && !is.data.frame(x$timeseries)) {
            paste(
              sprintf(" (Run %d)", block),
              if (!is.null(eye_suffix)) paste0(" (", eye_suffix, ")") else ""
            )
          } else {
            ""
          }
        )
        if (grepl("z", pupil_steps[i])) {
          y_units <- "(z)"
        } else {
          y_units <- "(a.u.)"
        }

        y_label <- paste("pupil size", y_units)

        # used when running `plot()` by itself (and thus plotting all steps)
        if (!only_liner_trend) {
          if (grepl("_detrend$", pupil_steps[i]) && !detrend_plotted) {
            # only attempt detrend overlay if detrend_fitted_values exists
            if ("detrend_fitted_values" %in% colnames(pupil_data)) {
              detrend_success <- plot_detrend_overlay(
                pupil_data,
                pupil_steps = pupil_steps,
                preview_n = preview_n,
                suppress_prompt = suppress_prompt
              )

              if (detrend_success) {
                detrend_plotted <- TRUE
              }
            } else {
              detrend_plotted <- TRUE
            }
          }
        } else {
          if (!detrend_plotted) {
            if ("detrend_fitted_values" %in% colnames(pupil_data)) {
              detrend_success <- plot_detrend_overlay(
                pupil_data,
                pupil_steps = pupil_steps,
                preview_n = preview_n,
                suppress_prompt = suppress_prompt
              )

              if (detrend_success) {
                detrend_plotted <- TRUE
              }
            } else {
              detrend_plotted <- TRUE
            }
          }
        }

        if (!is.null(params$next_step)) {
          plot_data <- epoch_n[[params$next_step[length(params$next_step)]]]
        } else {
          plot_data <- epoch_n[[pupil_steps[i]]]
        }

        is_placeholder <- "message" %in%
          colnames(epoch_n) &&
          any(epoch_n$message == "NO_VALID_SAMPLES")
        no_valid_data <- is.null(plot_data) || all(is.na(plot_data))

        if (is_placeholder || no_valid_data) {
          plot(
            NA,
            xlim = c(0, 1),
            ylim = c(0, 1),
            type = "n",
            xlab = "",
            ylab = "",
            main = title
          )
          text(
            0.5,
            0.5,
            "No valid samples\nin this segment.\n
            Please re-run with a different `report_seed`",
            cex = 0.8,
            col = "red"
          )
        } else {
          do.call(
            robust_plot,
            c(
              list(y = plot_data, x = epoch_n$time_scaled),
              plot_params,
              list(
                type = "l",
                col = colors[i],
                lwd = 2,
                main = title,
                xlab = "time (ms)",
                ylab = y_label
              )
            )
          )
        }
      }

      graphics::mtext(header, outer = TRUE, cex = 1.25, font = 2)

      if (plot_distributions) {
        plot_pupil_distribution(
          data = if (use_full) {
            full_pupil_data[[pupil_steps[i]]]
          } else {
            pupil_data[[pupil_steps[i]]]
          },
          color = colors[i],
          main = header,
          xlab = y_label,
          backuplab = "pupil size"
        )

        par(mfrow = c(1, preview_n), oma = c(0, 0, 3, 0))
      }
    }
    par(mfrow = c(1, preview_n), oma = c(0, 0, 3, 0))
  } else {
    for (i in seq_along(pupil_steps)) {
      # plot pre-decimation steps from the preserved full-resolution data at
      # the original sampling rate; decimated steps stay on the decimated data
      use_full <- has_decimation && !is_decimated_col(pupil_steps[i])
      this_data <- if (use_full) full_pupil_data else pupil_data
      this_hz <- if (use_full) full_hz else hz

      preview_window_indices <- round(preview_window * this_hz) + 1
      start_index <- preview_window_indices[1]
      end_index <- preview_window_indices[2]

      if (use_full) {
        # the preserved full-resolution frame can extend slightly past the
        # decimated frame the preview_window was derived from; clamp rather
        # than error so the same window can be drawn at full resolution
        start_index <- max(1, start_index)
        end_index <- min(nrow(this_data), end_index)
      }

      if (
        start_index < 1 ||
          start_index > nrow(this_data) ||
          end_index < 1 ||
          end_index > nrow(this_data) ||
          start_index >= end_index
      ) {
        log_error(
          "Invalid preview_window: start/end index out of range or invalid."
        )
      }

      sliced_pupil_data <- this_data[start_index:end_index, ]

      st <- this_data$time_orig[start_index]
      et <- this_data$time_orig[end_index]

      if (grepl("z", pupil_steps[i])) {
        y_units <- "(z)"
      } else {
        y_units <- "(a.u.)"
      }

      y_label <- paste("pupil size", y_units)

      do.call(
        robust_plot,
        c(
          list(
            y = sliced_pupil_data[[pupil_steps[i]]],
            x = sliced_pupil_data$time_scaled
          ),
          plot_params,
          list(
            type = "l",
            col = colors[i],
            lwd = 2,
            main = paste0(
              gsub("_", " > ", gsub("pupil_", "", pupil_steps[i])),
              if (is.list(x$timeseries) && !is.data.frame(x$timeseries)) {
                sprintf(" (Run %d)", block)
              } else {
                ""
              },
              "\n[",
              st,
              " - ",
              et,
              " ms] | ",
              "[index: ",
              preview_window_indices[1],
              " - ",
              preview_window_indices[2],
              "]"
            ),
            xlab = "time (secs)",
            ylab = y_label
          )
        )
      )

      if (plot_distributions) {
        plot_pupil_distribution(
          data = this_data[[pupil_steps[i]]],
          color = colors[i],
          main = paste(paste0(
            gsub("_", " > ", gsub("pupil_", "", pupil_steps[i])),
            if (is.list(x$timeseries) && !is.data.frame(x$timeseries)) {
              sprintf(" (Run %d)", block)
            } else {
              ""
            }
          )),
          xlab = y_label,
          backuplab = "pupil size"
        )
        par(mfrow = c(1, 1), oma = c(0, 0, 0, 0))
      }
    }

    par(mfrow = c(1, 1), oma = c(0, 0, 0, 0))
  }

  par(mfrow = c(1, 1), oma = c(0, 0, 0, 0))

  # add progressive summary plot at the end (if requested)
  if (add_progressive_summary) {
    log_info(
      "Creating progressive summary plot for block_{block}",
      verbose = verbose
    )

    tryCatch(
      {
        make_prog_summary_plot(
          pupil_data = pupil_data,
          pupil_steps = pupil_steps,
          preview_n = preview_n,
          plot_params = plot_params,
          run_id = if (is.list(x$timeseries) && !is.data.frame(x$timeseries)) {
            paste0("run-", sprintf("%02d", block))
          } else {
            "run-01"
          },
          cex = 1.15,
          full_pupil_data = full_pupil_data
        )

        log_success(
          "Progressive summary plot created successfully!",
          verbose = verbose
        )
      },
      error = function(e) {
        log_warn(
          "Could not create progressive summary plot: {e$message}",
          verbose = verbose
        )
      }
    )
  }

  # reset plotting parameters to prevent downstream issues
  par(mfrow = c(1, 1), oma = c(0, 0, 0, 0), mar = c(5, 4, 4, 2) + 0.1)
}

#' Identify decimated (downsample/bin) pupil columns
#'
#' A pupil column is considered "decimated" if it was produced by the
#' `downsample()`/`bin()` step or by any step that follows it. Such columns
#' are stored at the decimated sampling rate, whereas earlier columns should
#' be plotted from the preserved full-resolution data (see issue #294).
#'
#' @param col A character vector of pupil column names
#'
#' @return A logical vector that is `TRUE` for decimated columns
#'
#' @keywords internal
is_decimated_col <- function(col) {
  # match the `_downsample`/`_bin` operation suffix whether it is the final
  # column (e.g. `..._downsample`) or has later steps appended (e.g.
  # `..._downsample_detrend`, `..._bin_z`), while avoiding accidental matches
  # mid-token in unrelated names
  grepl("(_downsample|_bin)($|_)", col)
}

#' Retrieve the preserved full-resolution data for a block
#'
#' Returns the full-resolution (pre-decimation) time series for the requested
#' block, if it was preserved when a `downsample()`/`bin()` step ran. Used so
#' that diagnostic plots of earlier pipeline steps can be rendered at their
#' original sampling rate rather than the decimated rate (see issue #294).
#'
#' @param x An object of class `eyeris`
#' @param block The block number to retrieve
#'
#' @return A data frame of full-resolution time series data for the block, or
#' `NULL` if no pre-decimation data was preserved
#'
#' @keywords internal
get_pre_decimation_block <- function(x, block) {
  pre <- x$timeseries_pre_decimation
  if (is.null(pre)) {
    return(NULL)
  }
  if (is.data.frame(pre)) {
    return(pre)
  }
  pre[[paste0("block_", block)]]
}

#' Slice a fixed time window out of a time series for previewing
#'
#' Extracts the rows of `df` whose `time_secs` fall within
#' `[start_secs, start_secs + d)` and (re)computes a `time_scaled` column (in
#' milliseconds, relative to the window start). Used to draw a full-resolution
#' preview epoch over the same time window selected on the decimated data.
#'
#' @param df A data frame containing a `time_secs` column
#' @param start_secs The start of the window in seconds
#' @param d The window duration in seconds
#'
#' @return A data frame containing the rows that fall within the window
#'
#' @keywords internal
slice_epoch_window <- function(df, start_secs, d) {
  in_window <- df$time_secs >= start_secs & df$time_secs < (start_secs + d)
  epoch <- df[in_window, , drop = FALSE]
  epoch$time_scaled <- (epoch$time_secs - start_secs) * 1000
  epoch
}

#' Draw random epochs for plotting
#'
#' Generates random time segments from the time series data for preview plotting.
#'
#' @param x A data frame containing time series data
#' @param n Number of random epochs to draw
#' @param d Duration of each epoch in seconds
#' @param hz Sampling rate in Hz
#'
#' @return A list of data frames, each containing a random epoch segment
#'
#' @keywords internal
draw_random_epochs <- function(x, n, d, hz) {
  # get number of samples needed for specified duration
  n_samples <- ceiling(d * hz)

  min_time_secs <- min(x$time_secs, na.rm = TRUE)
  max_time_secs <- max(x$time_secs, na.rm = TRUE)

  if ((max_time_secs - min_time_secs) < d) {
    log_error("Example duration is longer than the duration of data.")
  }

  # get step size and ensure it's valid for the time range
  step_size <- 1 / hz
  time_range <- max_time_secs - d - min_time_secs

  # case: if step size is larger than available time range, adjust it
  if (step_size > time_range) {
    step_size <- time_range / 10 # use 10 steps as a reasonable minimum?
    if (step_size <= 0) {
      step_size <- 0.001 # fallback to 1ms if still invalid?
    }
  }

  drawn_epochs <- list()
  max_attempts <- 100 # prevent looping forever

  for (i in 1:n) {
    attempts <- 0
    valid_epoch_found <- FALSE

    while (attempts < max_attempts && !valid_epoch_found) {
      rand_start_secs <- sample(
        seq(min_time_secs, max_time_secs - d, by = step_size),
        1
      )
      rand_end_secs <- rand_start_secs + d

      epoch_data <- x |>
        dplyr::filter(time_secs >= rand_start_secs & time_secs < rand_end_secs)

      # ensure proper x-axis scaling with the time_orig column in ms
      epoch_data$time_scaled <- (epoch_data$time_secs - rand_start_secs) * 1000

      pupil_cols <- grep("^pupil_", colnames(epoch_data), value = TRUE)
      if (length(pupil_cols) > 0) {
        has_valid_data <- any(sapply(pupil_cols, function(col) {
          any(is.finite(epoch_data[[col]]))
        }))

        if (has_valid_data) {
          drawn_epochs[[i]] <- epoch_data
          valid_epoch_found <- TRUE
        } else {
          attempts <- attempts + 1
        }
      } else {
        drawn_epochs[[i]] <- epoch_data
        valid_epoch_found <- TRUE
      }
    }

    if (!valid_epoch_found) {
      placeholder_data <- data.frame(
        time_secs = c(rand_start_secs, rand_end_secs),
        time_orig = c(rand_start_secs, rand_end_secs),
        time_scaled = c(0, d * 1000),
        message = c("NO_VALID_SAMPLES", "NO_VALID_SAMPLES")
      )
      drawn_epochs[[i]] <- placeholder_data

      log_warn(
        "Randomly selected plot segment {i} had no valid samples. Please re-run with a different `report_seed`."
      )
    }
  }

  drawn_epochs
}

#' Robust plotting function with error handling
#'
#' A wrapper around base plotting functions that handles errors and missing
#' data gracefully.
#'
#' @param y The y-axis data to plot
#' @param x The x-axis data (optional, defaults to sequence)
#' @param ... Additional arguments passed to plot()
#'
#' @return No return value; creates a plot or displays warning messages
#'
#' @keywords internal
robust_plot <- function(y, x = NULL, ...) {
  tryCatch(
    {
      if (length(y) == 0 || all(is.na(y))) {
        log_warn("No finite data to plot.")
        return(invisible(NULL))
      }

      dots <- list(...)
      col_user <- if ("col" %in% names(dots)) dots$col else "blue"

      # store original y for getting NA positions
      y_orig <- y

      # if x is NULL, use 1:length(y)
      if (is.null(x)) {
        x_seq <- seq_along(y_orig)
      } else {
        x_seq <- x
      }

      # init placeholder line
      # handle case where x_seq has no finite values
      x_range <- range(x_seq, na.rm = TRUE, finite = TRUE)
      if (any(!is.finite(x_range))) {
        # fallback to default range if no finite values
        x_range <- c(0, length(x_seq))
      }

      plot(x_seq, ifelse(is.na(y_orig), NA, y_orig), xlim = x_range, ...)

      # add vertical lines where there are NAs (using x values if available)
      na_idx <- which(is.na(y_orig))
      if (length(na_idx) > 0) {
        abline(
          v = if (!is.null(x)) x_seq[na_idx] else na_idx,
          col = "black",
          lty = 2
        )
      }

      # replace NA with -1 after drawing NA lines for continuity
      y_clean <- y_orig
      y_clean[is.na(y_clean)] <- -1
      lines(x_seq, y_clean, col = col_user)
    },
    error = function(e) {
      log_warn("An error occurred during plotting: {e$message}")
    },
    warning = function(w) {
      log_warn("A warning occurred during plotting: {w$message}")
    }
  )
}

#' Plot pupil distribution histogram
#'
#' Creates a histogram of pupil size distribution with customizable parameters.
#'
#' @param data The pupil data to plot
#' @param color The color for the histogram bars
#' @param main The main title for the plot
#' @param xlab The x-axis label
#' @param backuplab A backup label if xlab is NULL
#'
#' @return No return value; creates a histogram plot
#'
#' @keywords internal
plot_pupil_distribution <- function(data, color, main, xlab, backuplab = NULL) {
  # safely handle user's current options
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  par(mfrow = c(1, 1), oma = c(0, 0, 0, 0))

  new_xlab <- if (!is.null(xlab)) {
    xlab
  } else if (!is.null(backuplab)) {
    backuplab
  } else {
    "pupil size"
  }

  # keep only finite samples so the Freedman-Diaconis rule (and the plot) do not
  # choke on the NA/NaN/Inf values that are present in the raw pupil signal
  finite_data <- data[is.finite(data)]

  if (length(finite_data) < 2) {
    # not enough data to build a histogram -- draw an informative empty panel
    # rather than letting hist() error out
    plot(
      NA,
      xlim = c(0, 1),
      ylim = c(0, 1),
      type = "n",
      xlab = new_xlab,
      ylab = "frequency (count)",
      main = main
    )
    text(
      0.5,
      0.5,
      "Not enough data\nto plot distribution",
      cex = 0.9,
      col = "red"
    )
    return(invisible(NULL))
  }

  # fall back to the default (Sturges) breaks if Freedman-Diaconis fails, e.g.
  # for a near-constant signal where the IQR is 0
  h <- tryCatch(
    hist(finite_data, breaks = "FD", plot = FALSE),
    error = function(e) hist(finite_data, plot = FALSE)
  )

  # a white bar outline gives nice separation when there are only a handful of
  # bars, but once the bars get thin the outline completely covers the fill and
  # the histogram renders blank. the raw step has the widest spread (and thus
  # the most Freedman-Diaconis bins), so its histogram was the one disappearing
  # in multi-run reports. drop the outline once there are too many bars so the
  # distribution always stays visible.
  bar_border <- if (length(h$counts) <= 100) "white" else NA

  plot(
    h,
    main = main,
    xlab = new_xlab,
    ylab = "frequency (count)",
    col = color,
    border = bar_border,
    freq = TRUE
  )
}

#' Draw vertical lines at NA positions
#'
#' Adds vertical dashed lines at positions where y values are NA.
#'
#' @param x The x-axis values
#' @param y The y-axis values
#' @param ... Additional arguments passed to abline()
#'
#' @return No return value; adds lines to the current plot
#'
#' @keywords internal
draw_na_lines <- function(x, y, ...) {
  na_idx <- which(is.na(y))
  abline(v = x[na_idx], col = "black", lty = 2, ...)
}

#' Internal helper to plot detrending overlay
#'
#' This function replicates the exact detrending visualization from the
#' `glassbox()` interactive preview mode. It uses `robust_plot()` to show the
#' most recent detrended pupil signal overlaid with the fitted linear trend.
#'
#' @param pupil_data A single block of pupil time series data
#' (e.g. `eyeris$timeseries$block_1`)
#' @param preview_n Number of columns for `par(mfrow)`. Default = 3.
#' @param plot_params A named list of additional parameters to forward to
#' `robust_plot()`
#' @param suppress_prompt Logical. Whether to skip prompting. Default = TRUE.
#'
#' @return Logical indicating whether detrend overlay was plotted successfully
#'
#' @keywords internal
plot_detrend_overlay <- function(
  pupil_data,
  pupil_steps,
  preview_n = preview_n,
  plot_params = list(),
  suppress_prompt = TRUE
) {
  # store current par settings to restore them in case func returns early
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par), add = TRUE)

  par(mfrow = c(1, 1), oma = c(0, 0, 0, 0))

  detrend_step <- grep("_detrend$", pupil_steps, value = TRUE)

  all_cols <- colnames(pupil_data)
  detrend_fitted_index <- which(all_cols == "detrend_fitted_values")

  # guard if detrend_fitted_values exists and has a valid previous column
  if (length(detrend_fitted_index) == 0) {
    log_warn("detrend_fitted_values not found in eyeris S3 object.")
    par(mfrow = c(1, preview_n), oma = c(0, 0, 3, 0))
    return(FALSE)
  }

  if (detrend_fitted_index <= 1) {
    log_warn(
      "No previous pupil column found to plot detrend overlay against.",
      "This can happen when detrend is the only preprocessing step enabled."
    )
    # restore main plotting func layout
    par(mfrow = c(1, preview_n), oma = c(0, 0, 3, 0))
    return(FALSE)
  }

  prev_col <- all_cols[detrend_fitted_index - 1]

  # ensure prev col is a pupil col
  if (!grepl("^pupil_", prev_col)) {
    log_warn(
      "Previous column is not a pupil column. Cannot plot detrend overlay."
    )
    # restore main plotting func layout
    par(mfrow = c(1, preview_n), oma = c(0, 0, 3, 0))
    return(FALSE)
  }

  ydat <- pupil_data[[prev_col]]
  xdat <- pupil_data$time_secs

  do.call(
    robust_plot,
    c(
      list(y = ydat, x = xdat),
      plot_params,
      list(
        type = "l",
        col = "black",
        lwd = 2,
        main = paste0(
          "detrend:\n",
          gsub("_", " > ", gsub("pupil_", "", detrend_step))
        ),
        xlab = "tracker time (s)",
        ylab = "pupil size (a.u.)"
      )
    )
  )

  lines(
    pupil_data$time_secs,
    pupil_data$detrend_fitted_values,
    type = "l",
    col = "blue",
    lwd = 2,
    lty = 1
  )

  legend(
    "topleft",
    legend = c("pupil time series", "linear trend"),
    col = c("black", "blue"),
    lwd = 2,
    lty = c(1, 1)
  )

  par(mfrow = c(1, preview_n), oma = c(0, 0, 3, 0))
  if (!suppress_prompt) {
    prompt_user()
  }

  return(TRUE)
}

#' Create gaze heatmap of eye coordinates
#'
#' Creates a heatmap showing the distribution of eye_x and eye_y coordinates
#' across the entire screen area. The heatmap shows where the participant
#' looked most frequently during the recording period.
#'
#' @param eyeris An object of class `eyeris` derived from [eyeris::load_asc()]
#' @param block Block number to plot (default: 1)
#' @param screen_width Screen width in pixels from `eyeris$info$screen.x`
#' @param screen_height Screen height in pixels from `eyeris$info$screen.y`
#' @param n_bins Number of bins for the heatmap grid (default: 50)
#' @param col_palette Color palette for the heatmap (default: "viridis")
#' @param main Title for the plot (default: "Fixation Heatmap")
#' @param xlab X-axis label (default: "Screen X (pixels)")
#' @param ylab Y-axis label (default: "Screen Y (pixels)")
#' @param sample_rate Sample rate in Hz (optional)
#' @param eye_suffix Eye suffix for binocular data (default: NULL)
#'
#' @return No return value; creates a heatmap plot
#'
#' @examples
#' demo_data <- eyelink_asc_demo_dataset()
#' eyeris_preproc <- glassbox(demo_data)
#' plot_gaze_heatmap(eyeris = eyeris_preproc, block = 1)
#'
#' @export
plot_gaze_heatmap <- function(
  eyeris,
  block = 1,
  screen_width = NULL,
  screen_height = NULL,
  n_bins = 50,
  col_palette = "viridis",
  main = "Gaze Heatmap",
  xlab = "Screen X (pixels)",
  ylab = "Screen Y (pixels)",
  sample_rate = NULL,
  eye_suffix = NULL
) {
  if (inherits(eyeris, "eyeris")) {
    block_str <- paste0("block_", block)
    if (is.null(screen_width)) {
      screen_width <- eyeris$info$screen.x
    }
    if (is.null(screen_height)) {
      screen_height <- eyeris$info$screen.y
    }

    df <- eyeris$timeseries[[block_str]]
    if (!is.data.frame(df)) {
      log_warn("Block not found in eyeris object.")
      return(invisible(NULL))
    }
  } else {
    df <- eyeris
    if (is.null(screen_width) || is.null(screen_height)) {
      log_error(
        "Screen width and height must be provided with data frame inputs."
      )
    }
  }

  if (!all(c("eye_x", "eye_y") %in% colnames(df))) {
    log_warn("eye_x and/or eye_y columns not found in input data.")
    return(invisible(NULL))
  }

  valid_coords <- !is.na(df$eye_x) & !is.na(df$eye_y)
  if (sum(valid_coords) == 0) {
    log_warn("No valid eye coordinates found")
    return(invisible(NULL))
  }

  x_coords <- df$eye_x[valid_coords]
  y_coords <- df$eye_y[valid_coords]

  if (!is.null(eye_suffix)) {
    main <- paste0(main, " - ", eye_suffix)
  }

  tryCatch(
    {
      dens <- MASS::kde2d(
        x_coords,
        y_coords,
        n = n_bins,
        lims = c(0, screen_width, 0, screen_height)
      )
      norm_density <- dens$z / max(dens$z, na.rm = TRUE)

      if (col_palette == "viridis") {
        colors <- viridis::viridis(100)
      } else if (col_palette == "plasma") {
        colors <- viridis::plasma(100)
      } else if (col_palette == "inferno") {
        colors <- viridis::inferno(100)
      } else if (col_palette == "magma") {
        colors <- viridis::magma(100)
      } else {
        colors <- grDevices::heat.colors(100)
      }

      fields::image.plot(
        x = dens$x,
        y = dens$y,
        z = t(norm_density)[, rev(seq_len(nrow(norm_density)))],
        col = colors,
        main = main,
        xlab = xlab,
        ylab = ylab,
        xlim = c(0, screen_width),
        ylim = c(screen_height, 0),
        legend.lab = "Normalized density",
        legend.line = 2.5,
        zlim = c(0, 1)
      )
      rect(0, 0, screen_width, screen_height, border = "black", lwd = 2)
      points(
        screen_width / 2,
        screen_height / 2,
        pch = 3,
        col = "red",
        cex = 1.5
      )
    },
    error = function(e) {
      plot(
        x_coords,
        y_coords,
        pch = 16,
        cex = 0.5,
        col = grDevices::adjustcolor("blue", alpha.f = 0.6),
        main = main,
        xlab = xlab,
        ylab = ylab,
        xlim = c(0, screen_width),
        ylim = c(screen_height, 0)
      )
      rect(0, 0, screen_width, screen_height, border = "black", lwd = 2)
      points(
        screen_width / 2,
        screen_height / 2,
        pch = 3,
        col = "red",
        cex = 1.5
      )
    }
  )
}

#' Plot binocular correlation between left and right eye data
#'
#' Creates correlation plots showing the relationship between left and right eye
#' measurements for pupil size, x-coordinates, and y-coordinates. This function
#' is useful for validating binocular data quality and assessing the correlation
#' between the two eyes.
#'
#' @param eyeris An object of class `eyeris` derived from [eyeris::load_asc()]
#'   with binocular data, or a list containing `left` and `right` eyeris objects
#'   (from `binocular_mode = "both"`)
#' @param block Block number to plot (default: 1)
#' @param variables Variables to plot correlations for. Defaults to
#'   `c("pupil", "x", "y")` for pupil size, x-coordinates, and y-coordinates
#' @param main Title for the overall plot (default: "Binocular Correlation")
#' @param col_palette Color palette for the plots (default: "viridis")
#' @param sample_rate Sample rate in Hz (optional, for time-based sampling)
#' @param verbose Logical flag to indicate whether to print status messages
#'   (default: TRUE)
#'
#' @return No return value; creates correlation plots
#'
#' @examples
#' # For binocular data loaded with binocular_mode = "both"
#' binocular_data <- load_asc(eyelink_asc_binocular_demo_dataset(), binocular_mode = "both")
#' plot_binocular_correlation(binocular_data)
#'
#' # For binocular data loaded with binocular_mode = "average"
#' # (correlation plot will show original left vs right before averaging)
#' avg_data <- load_asc(eyelink_asc_binocular_demo_dataset(), binocular_mode = "average")
#' plot_binocular_correlation(avg_data$raw_binocular_object)
#'
#' @export
plot_binocular_correlation <- function(
  eyeris,
  block = 1,
  variables = c("pupil", "x", "y"),
  main = "",
  col_palette = "viridis",
  sample_rate = NULL,
  verbose = TRUE
) {
  # check if a binocular object (from binocular_mode = "both")
  if (is_binocular_object(eyeris)) {
    left_data <- eyeris$left
    right_data <- eyeris$right
    has_binocular <- TRUE
  } else {
    # check if a regular eyeris object with binocular columns
    if (all(c("left", "right") %in% names(eyeris))) {
      left_data <- eyeris$left
      right_data <- eyeris$right
      has_binocular <- TRUE
    } else {
      left_data <- eyeris
      right_data <- eyeris
      has_binocular <- isTRUE(eyeris$binocular)
    }
  }

  block_str <- paste0("block_", block)

  if (has_binocular) {
    if (!block_str %in% names(left_data$timeseries)) {
      log_warn("Block {block} not found in left eye data")
    }
    if (!block_str %in% names(right_data$timeseries)) {
      log_warn("Block {block} not found in right eye data")
    }

    left_df <- left_data$timeseries[[block_str]]
    right_df <- right_data$timeseries[[block_str]]

    # require exact match for pupil_raw; if not present, skip plot with message
    if (!"pupil_raw" %in% colnames(left_df)) {
      plot.new()
      title(main = "Skipped: No pupil_raw column found in left eye data")
      log_warn("Skipped: No pupil_raw column found in left eye data")
      return(invisible(NULL))
    }
    pupil_col <- grep("^pupil_", colnames(left_df), value = TRUE)
    if (length(pupil_col) == 0) {
      log_warn("No pupil columns found in left eye data")
    }
    pupil_col <- pupil_col[1] # use the first pupil column

    pupil_col <- "pupil_raw"
    left_pupil <- left_df[[pupil_col]]
    right_pupil <- right_df[[pupil_col]]
    left_x <- left_df$eye_x
    left_y <- left_df$eye_y
    right_x <- right_df$eye_x
    right_y <- right_df$eye_y
  } else {
    # for regular eyeris objects, check for binocular columns
    if (!block_str %in% names(left_data$timeseries)) {
      log_warn("Block {block} not found in eyeris data")
    }

    df <- left_data$timeseries[[block_str]]

    if (!has_binocular) {
      log_warn(
        "No binocular columns (psl, psr, xpl, xpr, ypl, ypr) found in data.",
        "Use binocular_mode = 'both' when loading data to enable this function."
      )
    }

    left_pupil <- df$psl
    right_pupil <- df$psr
    left_x <- df$xpl
    left_y <- df$ypl
    right_x <- df$xpr
    right_y <- df$ypr
  }

  n_vars <- length(variables)
  if (n_vars == 1) {
    par(mfrow = c(1, 1))
  } else if (n_vars == 2) {
    par(mfrow = c(1, 2))
  } else {
    par(mfrow = c(1, 3))
  }

  if (col_palette == "viridis") {
    colors <- viridis::viridis(100)
  } else if (col_palette == "plasma") {
    colors <- viridis::plasma(100)
  } else if (col_palette == "inferno") {
    colors <- viridis::inferno(100)
  } else if (col_palette == "magma") {
    colors <- viridis::magma(100)
  } else {
    colors <- grDevices::heat.colors(100)
  }

  # create correlation plots for each variable
  for (var in variables) {
    if (var == "pupil") {
      left_var <- left_pupil
      right_var <- right_pupil
      xlab <- "Left Eye Pupil Size\n"
      ylab <- "Right Eye Pupil Size"
      title <- ""
    } else if (var == "x") {
      left_var <- left_x
      right_var <- right_x
      xlab <- "Left Eye X-Coordinate\n"
      ylab <- "Right Eye X-Coordinate"
      title <- ""
    } else if (var == "y") {
      left_var <- left_y
      right_var <- right_y
      xlab <- "Left Eye Y-Coordinate\n"
      ylab <- "Right Eye Y-Coordinate"
      title <- ""
    } else {
      log_warn("Unknown variable '{var}', skipping")
      next
    }

    # remove NA values for correlation calculation
    valid_data <- !is.na(left_var) & !is.na(right_var)
    if (sum(valid_data) == 0) {
      log_warn("No valid data for {var} correlation")
      next
    }

    left_clean <- left_var[valid_data]
    right_clean <- right_var[valid_data]

    cor_value <- cor(left_clean, right_clean, use = "complete.obs")

    tryCatch(
      {
        plot(
          left_clean,
          right_clean,
          pch = 16,
          cex = 0.5,
          main = sprintf("%s\nr = %.3f", title, cor_value),
          xlab = xlab,
          ylab = ylab,
          xlim = c(
            min(min(left_clean), min(right_clean)),
            max(max(left_clean), max(right_clean))
          ),
          ylim = c(
            min(min(left_clean), min(right_clean)),
            max(max(left_clean), max(right_clean))
          ),
          col = grDevices::adjustcolor("blue", alpha.f = 0.6)
        )
        abline(0, 1, col = "red", lwd = 2, lty = 2)
      },
      error = function(e) {
        log_warn("Error creating correlation plot for {var}: {e$message}")
      }
    )
  }

  mtext(main, outer = TRUE, cex = 1.25, font = 2, line = -1)

  # reset plotting parameters
  par(mfrow = c(1, 1))

  log_success(
    "Created binocular correlation plots for block {block}",
    verbose = verbose
  )
}