tests/testthat/test-correct_spike.R

In OpenSpecy: Analyze, Process, Identify, and Share Raman and (FT)IR Spectra

make_spike_test_spec <- function(values, axis = seq_along(values),
                                 id = "sample") {
  spectra <- data.frame(values)
  names(spectra) <- id
  as_OpenSpecy(axis, spectra)
}

test_that("correct_spike() validates dispatch and method-specific inputs", {
  expect_error(correct_spike(1:20), "OpenSpecy")

  clean <- make_spike_test_spec(seq_len(51))
  expect_error(
    correct_spike(clean, method = "prominence_fwhm"),
    "prominence_threshold.*width_threshold"
  )
  expect_error(correct_spike(clean, residual_window = 0),
               "residual_window")
  expect_error(correct_spike(clean, noise_multiplier = 0),
               "noise_multiplier")
  expect_error(correct_spike(clean, rel_height = 1.1), "rel_height")
  expect_error(
    correct_spike(clean, threshold = 5),
    "unused argument.*'threshold'"
  )
})

test_that("residual correction handles both spike signs and is idempotent", {
  axis <- seq(400, 1800, length.out = 101)
  baseline <- sin(axis / 200)
  positive <- negative <- baseline
  positive[51] <- positive[51] + 20
  negative[61] <- negative[61] - 20
  original <- as_OpenSpecy(
    axis,
    data.frame(positive = positive, negative = negative)
  )
  attr(original, "source_tag") <- "residual fixture"

  corrected <- correct_spike(
    original, method = "residual", interpolation_points = 5L
  )
  diagnostic <- attr(corrected, "automatic_spike")

  expect_true(diagnostic$applied)
  expect_identical(diagnostic$method, "residual")
  expect_setequal(diagnostic$affected_spectra, c("positive", "negative"))
  expect_equal(diagnostic$before_count, 2L)
  expect_equal(diagnostic$after_count, 0L)
  expect_equal(unname(corrected$spectra[51, "positive"]), baseline[51],
               tolerance = 0.01)
  expect_equal(unname(corrected$spectra[61, "negative"]), baseline[61],
               tolerance = 0.01)
  expect_identical(
    correct_spike(corrected, method = "residual", interpolation_points = 5L),
    corrected
  )
})

test_that("descending axes preserve orientation for positive and negative residuals", {
  descending_axis <- rev(seq(400, 1800, length.out = 121))
  baseline <- 2 + 0.01 * descending_axis
  cases <- list(
    list(direction = "positive", delta = 30),
    list(direction = "negative", delta = -30)
  )

  for (case in cases) {
    values <- baseline
    values[61] <- values[61] + case$delta
    # as_OpenSpecy() canonicalizes new input; reverse a valid object explicitly
    # to exercise the orientation accepted from existing OpenSpecy workflows.
    original <- as_OpenSpecy(
      rev(descending_axis),
      data.frame(sample = rev(values))
    )
    original$wavenumber <- rev(original$wavenumber)
    original$spectra <- original$spectra[
      rev(seq_len(nrow(original$spectra))), , drop = FALSE
    ]
    attr(original, "axis_source") <- "descending fixture"

    corrected <- correct_spike(
      original,
      method = "residual",
      direction = case$direction,
      interpolation_points = 5L
    )
    diagnostic <- attr(corrected, "automatic_spike")

    expect_true(all(diff(corrected$wavenumber) < 0))
    expect_identical(corrected$wavenumber, original$wavenumber)
    expect_identical(corrected$metadata, original$metadata)
    expect_identical(dimnames(corrected$spectra),
                     dimnames(original$spectra))
    expect_identical(attr(corrected, "axis_source"), "descending fixture")
    expect_equal(as.numeric(corrected$spectra[, 1]), baseline,
                 tolerance = 1e-12)
    expect_true(diagnostic$applied)
    expect_true(all(diagnostic$corrected_regions$region_min <=
                      diagnostic$corrected_regions$region_max))
    expect_identical(
      correct_spike(
        corrected,
        method = "residual",
        direction = case$direction,
        interpolation_points = 5L
      ),
      corrected
    )
  }
})

test_that("residual correction preserves the OpenSpecy contract", {
  axis <- cumsum(seq(0.5, 1.5, length.out = 121))
  baseline <- 2 + 0.03 * axis
  values <- baseline
  values[70] <- values[70] + 25
  original <- make_spike_test_spec(values, axis, id = "irregular")
  original$metadata$sample_note <- "kept"
  attr(original, "source_tag") <- list(owner = "test")
  original_attributes <- attributes(original)

  detection <- OpenSpecy:::.detect_spikes(
    original,
    method = "residual",
    interpolation_points = 5L
  )
  corrected <- correct_spike(
    original, method = "residual", interpolation_points = 5L
  )
  changed <- corrected$spectra != original$spectra
  changed[is.na(changed)] <- FALSE

  expect_identical(corrected$wavenumber, original$wavenumber)
  expect_identical(dim(corrected$spectra), dim(original$spectra))
  expect_identical(dimnames(corrected$spectra), dimnames(original$spectra))
  expect_identical(corrected$metadata, original$metadata)
  expect_identical(attr(corrected, "source_tag"),
                   original_attributes$source_tag)
  expect_identical(class(corrected), class(original))
  expect_false(any(changed & !detection$flagged))
  expect_equal(unname(corrected$spectra[70, 1]), baseline[70],
               tolerance = 1e-12)
})

test_that("a clean spectrum is an exact no-op", {
  axis <- cumsum(seq(0.8, 1.2, length.out = 101))
  clean <- make_spike_test_spec(3 + 0.01 * axis, axis)
  attr(clean, "custom") <- "unchanged"

  expect_identical(correct_spike(clean, interpolation_points = 5L), clean)
  expect_null(attr(clean, "automatic_spike"))
})

test_that("the detector exposes a stable reusable result structure", {
  values <- rep(0, 101)
  values[51] <- 30
  detected <- OpenSpecy:::.detect_spikes(
    make_spike_test_spec(values),
    interpolation_points = 5L
  )

  expect_identical(
    names(detected),
    c("method", "parameters", "candidates", "flagged",
      "candidate_count", "correctable_count", "reason")
  )
  expect_identical(
    names(detected$candidates),
    c("spectrum_index", "spectrum_id", "direction", "peak_index",
      "peak_wavenumber", "start_index", "end_index", "region_min",
      "region_max", "residual", "score", "prominence", "width",
      "prominence_width_ratio", "correctable", "reason")
  )
  expect_identical(dim(detected$flagged), c(length(values), 1L))
  expect_equal(detected$correctable_count, 1L)
  expect_true(detected$flagged[51, 1])
})

test_that("the default MAD prominence method matches the supplied workflow", {
  baseline <- rep(5, 101)
  positive <- negative <- baseline
  positive[51] <- 105
  negative[61] <- -95
  original <- as_OpenSpecy(
    seq_len(101), data.frame(positive = positive, negative = negative)
  )
  attr(original, "source_tag") <- "NCL fixture"

  corrected <- correct_spike(original)
  diagnostic <- attr(corrected, "automatic_spike")

  expect_true(diagnostic$applied)
  expect_identical(diagnostic$method, "mad_prominence_width")
  expect_identical(diagnostic$parameters$noise_multiplier, 10)
  expect_identical(diagnostic$parameters$width_threshold, 2)
  expect_identical(diagnostic$parameters$interpolation_points, 5L)
  expect_setequal(diagnostic$affected_spectra, c("positive", "negative"))
  expect_equal(corrected$spectra[49:53, "positive"], rep(5, 5))
  expect_equal(corrected$spectra[59:63, "negative"], rep(5, 5))
  expect_identical(corrected$metadata, original$metadata)
  expect_identical(attr(corrected, "source_tag"), "NCL fixture")
  expect_identical(correct_spike(corrected), corrected)
})

test_that("base peak intervals reproduce frozen pracma findpeaks results", {
  signal <- c(0, 1, 3, 1, 0, 0, -1, -3, -1, 0)
  positive <- OpenSpecy:::.spike_findpeak_intervals(signal, threshold = 2)
  negative <- OpenSpecy:::.spike_findpeak_intervals(-signal, threshold = 2)

  expect_equal(positive, data.frame(
    peak = 3L, start = 1L, end = 5L, prominence = 3, width = 4
  ))
  expect_equal(negative, data.frame(
    peak = 8L, start = 6L, end = 10L, prominence = 3, width = 4
  ))
  expect_equal(OpenSpecy:::.spike_raw_mad_diff(c(0, 1, 2, 20, 3)), 8.5)
})

test_that("default MAD correction uses conservative one-sided edge values", {
  values <- rep(0, 31)
  values[2] <- 50
  corrected <- correct_spike(make_spike_test_spec(values))
  diagnostic <- attr(corrected, "automatic_spike")

  expect_true(diagnostic$applied)
  expect_equal(corrected$spectra[1:4, 1], rep(0, 4))
  expect_false(any(!is.finite(corrected$spectra)))
})

test_that("the Coca-Lopez manual thresholds use sample-unit width", {
  axis <- seq_len(201)
  baseline <- sin(axis / 30)
  values <- baseline
  values[101] <- values[101] + 50
  original <- make_spike_test_spec(values, axis)

  detection <- OpenSpecy:::.detect_spikes(
    original,
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 40,
    width_threshold = 4,
    rel_height = 0.8,
    interpolation_points = 5L
  )
  corrected <- correct_spike(
    original,
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 40,
    width_threshold = 4,
    rel_height = 0.8,
    interpolation_points = 5L
  )

  expect_equal(detection$candidates$peak_index, 101L)
  expect_lt(detection$candidates$width, 4)
  expect_gt(detection$candidates$prominence, 40)
  expect_equal(unname(corrected$spectra[101, 1]), baseline[101],
               tolerance = 0.01)
  expect_identical(attr(corrected, "automatic_spike")$parameters$rel_height,
                   0.8)
})

test_that("prominence/FWHM ratio mode applies the paper's upper Z rule", {
  axis <- seq_len(801)
  baseline <- sin(2 * pi * axis / 25) + 0.15 * sin(2 * pi * axis / 7)
  values <- baseline
  values[401] <- values[401] + 30
  original <- make_spike_test_spec(values, axis)

  detection <- OpenSpecy:::.detect_spikes(
    original,
    method = "prominence_fwhm_ratio",
    direction = "positive",
    min_peaks = 20L,
    interpolation_points = 5L
  )
  corrected <- correct_spike(
    original,
    method = "prominence_fwhm_ratio",
    direction = "positive",
    min_peaks = 20L,
    interpolation_points = 5L
  )

  expect_identical(detection$reason, "detected")
  expect_equal(detection$candidates$peak_index, 401L)
  expect_gt(detection$candidates$score, 3.5)
  expect_lt(abs(unname(corrected$spectra[401, 1]) - baseline[401]), 0.1)

  too_few <- OpenSpecy:::.detect_spikes(
    make_spike_test_spec(c(0, 1, 0, 1, 0, 20, 0, 1, 0)),
    method = "prominence_fwhm_ratio",
    direction = "positive",
    min_peaks = 20L,
    interpolation_points = 1L
  )
  expect_identical(too_few$reason, "insufficient_peaks")
  expect_equal(too_few$candidate_count, 0L)
})

test_that("adjacent paper spikes share clean interpolation neighbors", {
  values <- rep(5, 201)
  values[100:101] <- 105
  original <- make_spike_test_spec(values)

  corrected <- correct_spike(
    original,
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 40,
    width_threshold = 4,
    interpolation_points = 5L
  )
  diagnostic <- attr(corrected, "automatic_spike")

  expect_true(diagnostic$applied)
  expect_equal(corrected$spectra[100:101, 1], c(5, 5))
  expect_equal(diagnostic$corrected_regions$start_index, 100L)
  expect_equal(diagnostic$corrected_regions$end_index, 101L)
  expect_false(any(!is.finite(corrected$spectra)))
})

test_that("boundary intervals are rejected without wrapping", {
  values <- rep(0, 101)
  values[2] <- 50
  original <- make_spike_test_spec(values)

  corrected <- correct_spike(
    original,
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 10,
    width_threshold = 4,
    interpolation_points = 5L
  )
  diagnostic <- attr(corrected, "automatic_spike")

  expect_identical(corrected$spectra, original$spectra)
  expect_false(diagnostic$applied)
  expect_identical(diagnostic$reason, "no_correctable_regions")
  expect_true("boundary_interval" %in% diagnostic$rejected_regions$reason)
})

test_that("narrow real bands are rejected conservatively", {
  axis <- seq_len(201)
  narrow_band <- 100 * exp(-0.5 * ((axis - 101) / 1.5)^2)
  original <- make_spike_test_spec(narrow_band, axis)

  residual <- correct_spike(
    original, method = "residual", interpolation_points = 5L
  )
  residual_diagnostic <- attr(residual, "automatic_spike")
  expect_identical(residual$spectra, original$spectra)
  expect_false(residual_diagnostic$applied)
  expect_setequal(
    residual_diagnostic$rejected_regions$reason,
    c("candidate_too_wide", "spectral_band_shoulder")
  )

  calibrated <- correct_spike(
    original,
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 10,
    width_threshold = 1,
    interpolation_points = 5L
  )
  expect_identical(calibrated, original)
})

test_that("the Fig. 6-style broad band is not truncated", {
  axis <- seq_len(201)
  truth <- 100 * exp(-((axis - 101) / 18)^2)
  values <- truth
  values[110:112] <- values[110:112] + c(200, 400, 200)
  original <- make_spike_test_spec(values, axis)

  corrected <- correct_spike(
    original,
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 40,
    width_threshold = 4,
    rel_height = 0.8,
    interpolation_points = 10L
  )

  expect_true(attr(corrected, "automatic_spike")$applied)
  expect_equal(corrected$spectra[110:112, 1], truth[110:112],
               tolerance = 1)
  expect_identical(corrected$spectra[-(110:112), 1],
                   original$spectra[-(110:112), 1])
})

test_that("correction preserves existing non-finite values and adds none", {
  axis <- seq_len(151)
  values <- sin(axis / 20)
  values[20] <- NA_real_
  values[90] <- values[90] + 30
  original <- make_spike_test_spec(values, axis)

  corrected <- correct_spike(
    original, method = "residual", interpolation_points = 5L
  )

  expect_identical(which(!is.finite(corrected$spectra)),
                   which(!is.finite(original$spectra)))
  expect_true(attr(corrected, "automatic_spike")$applied)
})

test_that("iterative correction retains safe test-map passes", {
  map_path <- system.file(
    "extdata", "CA_tiny_map.zip", package = "OpenSpecy"
  )
  expect_true(nzchar(map_path) && file.exists(map_path))
  original <- read_any(map_path) |>
    c_spec(range = "common", res = 6) |>
    manage_na(ig = c(NA, 0), type = "remove")

  corrected <- correct_spike(
    original,
    method = "residual", direction = "both",
    residual_threshold = 8, residual_window = 5
  )
  diagnostic <- attr(corrected, "automatic_spike")

  expect_true(diagnostic$applied)
  expect_identical(diagnostic$reason, "corrected_with_safeguards")
  expect_identical(diagnostic$before_count, 43L)
  expect_identical(diagnostic$after_count, 1L)
  expect_identical(diagnostic$pass_count, 3L)
  expect_equal(nrow(diagnostic$corrected_regions), 49L)
  expect_true("interpolation_no_change" %in%
                diagnostic$rejected_regions$reason)
  expect_false(identical(diagnostic$reason, "correctable_spikes_remain"))
  expect_equal(
    sum(as.matrix(original$spectra) != as.matrix(corrected$spectra)), 49L
  )
  expect_identical(corrected$wavenumber, original$wavenumber)
  expect_identical(dimnames(corrected$spectra), dimnames(original$spectra))
  expect_identical(corrected$metadata, original$metadata)
  expect_identical(
    correct_spike(
      corrected,
      method = "residual", direction = "both",
      residual_threshold = 8, residual_window = 5
    ),
    corrected
  )
})

test_that("idempotency preserves applied diagnostics after metadata changes", {
  values <- rep(0, 101)
  values[c(2, 51)] <- 50
  corrected <- correct_spike(
    make_spike_test_spec(values),
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 10,
    width_threshold = 4,
    interpolation_points = 5L
  )
  original_diagnostic <- attr(corrected, "automatic_spike")
  expect_true(original_diagnostic$applied)
  expect_true("boundary_interval" %in%
                original_diagnostic$rejected_regions$reason)

  annotated <- corrected
  annotated$metadata$note <- "metadata does not affect spike detection"
  repeated <- correct_spike(
    annotated,
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 10,
    width_threshold = 4,
    interpolation_points = 5L
  )

  expect_identical(repeated, annotated)
  expect_identical(attr(repeated, "automatic_spike"), original_diagnostic)
})

test_that("idempotency does not preserve stale diagnostics after data changes", {
  values <- rep(0, 101)
  values[51] <- 50
  corrected <- correct_spike(
    make_spike_test_spec(values),
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 10,
    width_threshold = 4,
    interpolation_points = 5L
  )
  expect_true(attr(corrected, "automatic_spike")$applied)

  changed <- corrected
  changed$spectra[2, 1] <- 50
  attempted <- correct_spike(
    changed,
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 10,
    width_threshold = 4,
    interpolation_points = 5L
  )
  diagnostic <- attr(attempted, "automatic_spike")

  expect_false(diagnostic$applied)
  expect_identical(diagnostic$method, "prominence_fwhm")
  expect_identical(diagnostic$reason, "no_correctable_regions")
  expect_true("boundary_interval" %in% diagnostic$rejected_regions$reason)
  expect_false(identical(diagnostic$result_signature,
                         attr(corrected, "automatic_spike")$result_signature))
})

test_that("idempotency does not preserve stale diagnostics after method changes", {
  values <- rep(0, 101)
  values[c(2, 51)] <- 50
  corrected <- correct_spike(
    make_spike_test_spec(values),
    method = "residual",
    residual_window = 5L,
    interpolation_points = 5L
  )
  expect_true(attr(corrected, "automatic_spike")$applied)

  attempted <- correct_spike(
    corrected,
    method = "prominence_fwhm",
    direction = "positive",
    prominence_threshold = 10,
    width_threshold = 4,
    interpolation_points = 5L
  )
  diagnostic <- attr(attempted, "automatic_spike")

  expect_false(diagnostic$applied)
  expect_identical(diagnostic$method, "prominence_fwhm")
  expect_identical(diagnostic$reason, "no_correctable_regions")
  expect_true("boundary_interval" %in% diagnostic$rejected_regions$reason)
  expect_false(identical(diagnostic$result_signature,
                         attr(corrected, "automatic_spike")$result_signature))
})