tests/testthat/test_calc_Huntley2006.R

In Luminescence: Comprehensive Luminescence Dating Data Analysis

## load data
set.seed(1)
data("ExampleData.Fading", envir = environment())
fading_data <- ExampleData.Fading$fading.data$IR50
data <- ExampleData.Fading$equivalentDose.data$IR50
ddot <- c(7.00, 0.004)
readerDdot <- c(0.134, 0.0067)
rhop <-
  analyse_FadingMeasurement(fading_data,
                            plot = FALSE,
                            verbose = FALSE,
                            n.MC = 10)

test_that("input validation", {
  testthat::skip_on_cran()

  rhop.test <- rhop
  rhop.test@originator <- "unexpected"

  expect_error(calc_Huntley2006(),
               "'data' should be of class 'data.frame'")
  expect_error(calc_Huntley2006("test"),
               "'data' should be of class 'data.frame'")
  expect_error(calc_Huntley2006(data.frame(a = 1, b = 2, c = 3)),
               "'data' should contain at least 4 dose points")
  expect_error(calc_Huntley2006(iris[, 1, drop = FALSE]),
               "'data' should have at least 2 columns")

  expect_error(calc_Huntley2006(data, fit.method = "test"),
               "'fit.method' should be one of 'SSE' or 'GOK'")
  expect_error(calc_Huntley2006(data, fit.method = "GOK", lower.bounds = 0),
               "'lower.bounds' should be of class 'numeric' and have length 4")

  expect_error(calc_Huntley2006(data, LnTn = list()),
               "'LnTn' should be of class 'data.frame'")
  expect_error(calc_Huntley2006(data, LnTn = data),
               "'LnTn' should be a data frame with 2 columns")
  expect_error(calc_Huntley2006(cbind(data, data), LnTn = data[, 1:2]),
               "When 'LnTn' is specified, 'data' should have only 2 or 3")
  expect_error(calc_Huntley2006(cbind(data, data[, 1])),
               "The number of columns in 'data' must be a multiple of 3")

  expect_error(calc_Huntley2006(data, rhop = 1),
               "'rhop' should have length 2")
  expect_error(calc_Huntley2006(data, rhop = "test"),
               "'rhop' should be of class 'numeric' or 'RLum.Results'")
  expect_error(calc_Huntley2006(data, rhop = rhop.test),
               "'rhop' has an unsupported originator")
  expect_error(calc_Huntley2006(data, rhop = c(0, 4e-7)),
               "'rhop' must be a positive number, the provided value was 0 \u00B1 4e-07")

  expect_error(calc_Huntley2006(data, rhop = rhop),
               "'ddot' should be of class 'numeric' and have length 2")
  expect_error(calc_Huntley2006(data, rhop = rhop,
                                ddot = "test", readerDdot = readerDdot),
               "'ddot' should be of class 'numeric' and have length 2")
  expect_error(calc_Huntley2006(data, rhop = rhop,
                                ddot = 0.4, readerDdot = readerDdot),
               "'ddot' should be of class 'numeric' and have length 2")
  expect_error(calc_Huntley2006(data, rhop = rhop, ddot = ddot),
               "'readerDdot' should be of class 'numeric' and have length 2")
  expect_error(calc_Huntley2006(data, rhop = rhop,
                                ddot = ddot, readerDdot = "test"),
               "'readerDdot' should be of class 'numeric' and have length 2")
  expect_error(calc_Huntley2006(data, rhop = rhop,
                                ddot = ddot, readerDdot = 0.13),
               "'readerDdot' should be of class 'numeric' and have length 2")
  expect_error(calc_Huntley2006(data, rhop = rhop,
                                ddot = ddot, readerDdot = readerDdot,
                                n.MC = 0),
               "'n.MC' should be a single positive integer value")
  expect_error(calc_Huntley2006(data, rhop = rhop,
                                ddot = ddot, readerDdot = readerDdot,
                                rprime = list()),
               "'rprime' should be of class 'numeric'")

  SW({
  expect_warning(calc_Huntley2006(data[, 1:2], rhop = rhop, n.MC = 2,
                                  ddot = ddot, readerDdot = readerDdot,
                                  fit.method = "GOK", plot = FALSE),
                 "'data' has only two columns: we assume that the errors")
  })

  set.seed(1)
  expect_warning(
      calc_Huntley2006(data = data[1:20, ], rhop = rhop, ddot = ddot,
                       readerDdot = c(0.002, 0.003), n.MC = 2,
                       plot = FALSE, verbose = FALSE),
      "Ln/Tn is smaller than the minimum computed LxTx value")

  expect_error(calc_Huntley2006(data = data[1:20, ], LnTn = data[1, c(2, 3)],
                                rhop = c(2, 2), ddot = c(7.00, 0.004),
                                readerDdot = c(0.01, 0.02), n.MC = 2,
                                verbose = FALSE),
               "Could not fit simulated curve, check suitability of model")
})

test_that("snapshot tests", {
  testthat::skip_on_cran()

  os <- tolower(Sys.info()[["sysname"]])
  snapshot.tolerance <- switch(os,
                               "linux" = 1.5e-5,
                               "darwin" = 8.0e-2,
                               "windows" = 8.0e-2)

  ## check extrapolation
  expect_snapshot_RLum(
      calc_Huntley2006(
        data = data,
        rhop = rhop,
        ddot = ddot,
        readerDdot = readerDdot,
        n.MC = 100,
        fit.method = "GOK",
        mode = "extrapolation",
        plot = FALSE),
      expect_snapshot_output = TRUE,
      tolerance = snapshot.tolerance)

  ## check force through origin SSE with wrong mode settings
  expect_snapshot_RLum(
      calc_Huntley2006(
        data = data,
        rhop = rhop,
        ddot = ddot,
        readerDdot = readerDdot,
        n.MC = 100,
        fit.method = "SSE",
        fit.force_through_origin = TRUE,
        mode = "extrapolation",
        plot = FALSE),
      expect_snapshot_output = TRUE,
      tolerance = max(snapshot.tolerance, 1.0e-2))

  ## SSE ... normal
  expect_snapshot_RLum(
      calc_Huntley2006(
        data = data,
        rhop = rhop,
        ddot = ddot,
        readerDdot = readerDdot,
        n.MC = 100,
        fit.method = "SSE",
        fit.force_through_origin = TRUE,
        mode = "interpolation",
        plot = FALSE),
      expect_snapshot_output = TRUE,
      tolerance = snapshot.tolerance)

  ## GOK normal
  expect_snapshot_RLum(
      calc_Huntley2006(
        data = data,
        rhop = rhop,
        ddot = ddot,
        readerDdot = readerDdot,
        n.MC = 100,
        fit.method = "GOK",
        fit.force_through_origin = TRUE,
        mode = "interpolation",
        plot = FALSE),
      expect_snapshot_output = TRUE,
      tolerance = if (os == "darwin") 1.7e-1 else max(snapshot.tolerance, 1.5e-3))
})

test_that("graphical snapshot tests", {
  testthat::skip_on_cran()
  testthat::skip_if_not_installed("vdiffr")

  set.seed(1)

  SW({
  vdiffr::expect_doppelganger("interpolation SSE",
                              calc_Huntley2006(data = data,
                                               rhop = rhop,
                                               ddot = ddot,
                                               readerDdot = readerDdot,
                                               fit.method = "SSE",
                                               mode = "interpolation",
                                               plot_all_DRC = FALSE,
                                               n.MC = 50))
  vdiffr::expect_doppelganger("extrapolation GOK",
                              calc_Huntley2006(data = data,
                                               rhop = rhop,
                                               ddot = ddot,
                                               readerDdot = readerDdot,
                                               fit.method = "GOK",
                                               mode = "extrapolation",
                                               plot_all_DRC = FALSE,
                                               n.MC = 50))
  })
})

test_that("further tests", {
  testthat::skip_on_cran()

  ## check warning for failed fits
  ## dataset provided by Christine Neudorf
  df <- structure(list(V1 = c(0L, 0L, 200L, 800L, 1500L, 3000L, 8000L
  ), V2 = c(0.439, -0.046720922, 1.988131642, 7.577744961, 12.87699795,
            18.50187874, 32.72443771),
  V3 = c(0.029, 0.01269548, 0.090232208,
         0.322242141, 0.546595156, 0.739308178, 1.285979033)),
  class = "data.frame", row.names = c(NA, -7L))

  set.seed(1)
  SW({
  expect_warning(calc_Huntley2006(
    data = df,
    LnTn = NULL,
    rhop = c(0.0000121549740899913, 4.91596040125088E-07),
    ddot = c(6.96, 0.29),
    readerDdot = c(0.094, 0.01),
    normalise = FALSE,
    fit.method = "SSE",
    summary = TRUE,
    plot = TRUE,
    n.MC = 100),
    regexp = "\\[calc\\_Huntley2006\\(\\)\\] Ln\\/Tn is smaller than the minimum computed LxTx value.")
  })

  ## failing to fit unfaded SSE model
  ## test derived from data provided by Salome Oehler
  input <- data.frame(dose = c(0, 1550, 9300, 37200, 71500,
                               0, 1550, 9300, 37200, 71500),
                      LxTx = c(0.79, 4.67, 14.41, 26.59, 24.88,
                               0.95, 4.87, 14.17, 21.98, 25.12),
                      LxTx_error = c(0.002, 0.006, 0.01, 0.04, 0.02,
                                     0.002, 0.011, 0.02, 0.03, 0.03))
  set.seed(1)
  expect_error(calc_Huntley2006(input,
                                rhop = c(6.5e-06, 2.0e-08),
                                ddot = c(8.5, 1.5),
                                fit.weights = "norm_inverse_std",
                                fit.method = "SSE",
                                mode = "extrapolation",
                                readerDdot = c(0.154, 0.1),
                                n.MC = 2),
               "Could not fit unfaded curve, check suitability of model and")

  ## more coverage
  expect_output(
    calc_Huntley2006(
      data = data[1:10, ],
      LnTn = data[1:10, c(2, 3)],
      rhop = rhop, ddot = ddot, readerDdot = readerDdot,
      rprime = c(0.01, 2.2, length.out = 500),
      n.MC = 2, plot = FALSE, verbose = FALSE, maxiter = 50, trace = TRUE))

  expect_error(
    calc_Huntley2006(
      data = iris[, 2:4],
      rhop = rhop, ddot = ddot, readerDdot = readerDdot,
      n.MC = 2, plot = FALSE, verbose = FALSE),
    "Unable to fit growth curve to measured data, try setting 'fit.bounds = FALSE'")
  expect_error(
    calc_Huntley2006(
      data = iris[, 2:4],
      rhop = rhop, ddot = ddot, readerDdot = readerDdot,
      n.MC = 2, plot = FALSE, verbose = FALSE, fit.bounds = FALSE),
    "Unable to fit growth curve to measured data$")
  expect_warning(expect_error(
    calc_Huntley2006(
      data = data,
      rhop = c(4e-5, 5e-7), ddot = c(8, 0.04), readerDdot = c(0.1, 0.006),
      n.MC = 2, fit.method = "GOK", plot = FALSE, verbose = FALSE),
    "Could not fit unfaded curve, check suitability of model and parameters"),
    "Ln is >10 % larger than the maximum computed LxTx value")
  expect_warning(expect_error(
    calc_Huntley2006(
      data = data[8:11, ],
      rhop = c(4e-5, 5e-7), ddot = c(8, 0.04), readerDdot = c(0.0001, 0.006),
      n.MC = 2, mode = "extrapolation", plot = FALSE, verbose = FALSE),
    "Simulated D0 is NA: either your input values are unsuitable"),
    "Ln is >10 % larger than the maximum computed LxTx value")
  expect_error(
    calc_Huntley2006(
      data = data,
      rhop = c(1e-3, 1e-7), ddot = ddot, readerDdot = readerDdot,
      n.MC = 2, plot = FALSE, verbose = FALSE),
    "All simulated Lx/Tx values are identical and approximately zero")
})

test_that("regression tests", {
  testthat::skip_on_cran()

  ## issue 660
  DRC <- data.frame(
      Dose = c(rep(0, 7), rep(500, 12), rep(1000, 6), rep(2000, 6)),
      Lx.Tx = c(1.79, rep(0.06, 6), rep(0.97, 12), rep(1.80, 6), rep(3.20, 6)),
      Lx.Tx.Err = c(0.038, rep(0.001, 6), rep(0.021, 12), rep(0.038, 6), rep(0.068, 6))
  )
  expect_s4_class(
      calc_Huntley2006(
          data = DRC,
          rhop = c(7.67e-7, 2.84e-7),
          fit.method = "GOK",
          readerDdot = c(0.0868, 0.005),
          ddot = c(2.372, 0.199),
      n.MC = 2, plot = FALSE, verbose = FALSE),
      "RLum.Results")

  ## issue 733
  expect_s4_class(
      calc_Huntley2006(data, rhop = c(4e-6, 5e-7), ddot = c(7, 0.004),
                       readerDdot = c(0.134, 0.0067), n.MC = 1,
                       mode = "extrapolation", plot = FALSE, verbose = FALSE),
      "RLum.Results")

  ## issue 1048
  set.seed(1)
  expect_silent(
      calc_Huntley2006(data, rhop = c(1e-7, 5e-7), ddot = c(7, 0.004),
                       readerDdot = c(0.134, 0.0067), n.MC = 1,
                       mode = "extrapolation", verbose = FALSE,
                       plot_all_DRC = FALSE, plot = TRUE))

  ## issue 1672
  df <- data.frame(dose = c(151.4, 152.0, 165.8, 136.1, 144.4, 123.4, 123.6, 127.0,
                          124.4, 118.6, 128.0, 110.7, 121.0, 124.0, 124.7, 123.6,
                          128.5, 131.4, 127.7, 131.0, 126.3, 115.4, 119.5, 331.6),
                 LxTx = c(5.334, 5.144, 6.805, 4.608, 4.642, 4.471, 4.227, 4.396,
                          4.256, 4.049, 4.408, 3.701, 4.187, 4.129, 4.043, 4.262,
                          4.254, 4.448, 4.330, 5.023, 4.317, 3.479, 3.815, 13.61),
                 LxTx.error = 0.5)
  set.seed(1)
  expect_s4_class(suppressWarnings(
      calc_Huntley2006(df, rhop = c(4e-5, 5e-7),
                       ddot = c(8, 0.04), readerDdot = c(0.01, 0.06),
                       mode = "extrapolation", n.MC = 1, verbose = FALSE)),
      "RLum.Results")
})