tests/testthat/test_fit_DoseResponseCurve.R

In Luminescence: Comprehensive Luminescence Dating Data Analysis

## load data
data(ExampleData.LxTxData, envir = environment())

## odd data set where the calculated De is negative for SSE
df_odd <- data.frame(
    dose = c(
      0, 0, 2.71828182845905, 2.74202785430992,
      2.76598131771852, 2.79014403079188, 2.814517821467, 2.83910453364916,
      2.86390602735137, 2.88892417883514, 2.91416088075237, 2.93961804228855,
      2.96529758930721, 2.99120146449558, 3.01733162751159, 3.04369005513209,
      3.07027874140241, 3.09709969778723, 3.1241549533227, 3.15144655477),
    LxTx = c(0.439, 0.851865307456887, 0.881530377359027, 0.881788046363334,
      0.882047940677405, 0.882310079298982, 0.882574481384203, 0.882841166248859,
      0.883110153369655, 0.883381462385488, 0.883655113098727, 0.883931125476509,
      0.884209519652028, 0.88449031592586, 0.884773534767266, 0.885059196815536,
      0.885347322881313, 0.885637933947948, 0.885931051172847, 0.886226695888845),
    LxTx.error = c(0.029, 1e-04, 0.0393190034426231, 0.0393304962836274,
      0.0393420883803856, 0.0393537805802564, 0.0393655737376631, 0.0393774687141504,
      0.03938946637844, 0.0394015676064878, 0.0394137732815415, 0.0394260842941973,
      0.0394385015424588, 0.0394510259317946, 0.0394636583751979, 0.0394763997932453,
      0.0394892511141562, 0.0395022132738537, 0.0395152872160239, 0.0395284738921779))

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

  ## object
  expect_error(
      fit_DoseResponseCurve("error"),
      "[fit_DoseResponseCurve()] 'object' should be of class 'data.frame'",
      fixed = TRUE)
  expect_error(fit_DoseResponseCurve(data.frame()),
               "'object' cannot be an empty data.frame")
  expect_error(fit_DoseResponseCurve(iris[, 1, drop = FALSE]),
               "'object' should have at least 2 columns")
  expect_error(fit_DoseResponseCurve(as.list(LxTxData)),
               "All elements of 'object' should be of class 'data.frame'")

  ## mode
  expect_error(
    fit_DoseResponseCurve(LxTxData, mode = "error"),
    "[fit_DoseResponseCurve()] 'mode' should be one of 'interpolation'",
    fixed = TRUE)
  expect_error(fit_DoseResponseCurve(LxTxData, mode = set_RLum("RLum.Results")),
               "'mode' should be one of 'interpolation', 'extrapolation' or")

  ## fit.method
  expect_error(
    fit_DoseResponseCurve(LxTxData, fit.method = "error"),
    "[fit_DoseResponseCurve()] 'fit.method' should be one of 'LIN', 'QDR'",
    fixed = TRUE)

  ## other arguments
  expect_error(fit_DoseResponseCurve(LxTxData,
                                     fit.force_through_origin = "error"),
               "'fit.force_through_origin' should be a single logical value")
  expect_error(fit_DoseResponseCurve(LxTxData, fit.weights = "error"),
               "'fit.weights' should be one of 'inverse_var', 'inverse_std', 'norm_inverse_std', a numeric vector or NULL")
  expect_error(fit_DoseResponseCurve(LxTxData, fit.weights = iris),
               "'fit.weights' should be of class 'character', 'numeric' or NULL")
  expect_error(fit_DoseResponseCurve(LxTxData, fit.weights = c(1, 2)),
               "'fit.weights' should have length 7")
  expect_error(fit_DoseResponseCurve(LxTxData,
                                     fit.includingRepeatedRegPoints = "error"),
               "'fit.includingRepeatedRegPoints' should be a single logical")
  expect_error(fit_DoseResponseCurve(LxTxData,
                                     fit.NumberRegPoints = "error"),
               "'fit.NumberRegPoints' should be a single positive integer value")
  expect_error(fit_DoseResponseCurve(LxTxData, fit.NumberRegPointsReal = "error"),
               "'fit.NumberRegPointsReal' should be a single positive integer value")
  expect_error(fit_DoseResponseCurve(LxTxData,
                                     fit.bounds = "error"),
               "'fit.bounds' should be a single logical value")
  expect_error(fit_DoseResponseCurve(LxTxData,
                                     n.MC = "error"),
               "'n.MC' should be a single positive integer value")
  expect_error(fit_DoseResponseCurve(LxTxData,
                                     verbose = "error"),
               "'verbose' should be a single logical value")

  ## shorten dataframe
  expect_warning(fit_DoseResponseCurve(LxTxData[1:2, ], verbose = FALSE),
                 "Fitting a non-linear least-squares model requires at least 3")
  expect_warning(fit_DoseResponseCurve(LxTxData[1:3, ], fit.method = "GOK",
                                       verbose = FALSE),
                 "requires at least 4 dose points besides the natural, 'fit.method'")
  expect_warning(fit_DoseResponseCurve(LxTxData[1:3, ], fit.method = "GOK",
                                       mode = "extrapolation",
                                       verbose = FALSE),
                 "requires at least 4 dose points, 'fit.method' changed to")

  ## wrong combination of fit.method and mode
  expect_error(
    fit_DoseResponseCurve(LxTxData, fit.method = "DSE",
                     mode = "extrapolation"),
    "Mode 'extrapolation' for fitting method 'DSE' not supported",
    fixed = TRUE)

  ## deprecated option
  SW({
  expect_warning(fit_DoseResponseCurve(LxTxData, fit.weights = FALSE),
                 "'fit.weight' no longer accepts a logical value, reset automatically to NULL")
  expect_warning(fit_DoseResponseCurve(LxTxData, fit.weights = TRUE),
                 "'fit.weight' no longer accepts a logical value, reset automatically to inverse_var")
  expect_warning(res <- fit_DoseResponseCurve(LxTxData, fit.method = "EXP"),
                 "'fit.method = \"EXP\"' was deprecated in v1.3.0, use 'fit.method = \"SSE\"' instead")
  expect_equal(res@data$De$Fit,
               "SSE")
  })
})

test_that("weird LxTx values", {
  testthat::skip_on_cran()

  LxTx <- structure(list(
    Dose = c(0, 250, 500, 750, 1000, 1500, 0, 500, 500),
    LxTx = c(1, Inf, 0, -Inf, Inf, 0, Inf, -0.25, 2),
    LxTx.Error = c(1.58133646008685, Inf, 0, Inf, Inf, 0, Inf, 1.41146256149428, 3.16267292017369)),
    class = "data.frame", row.names = c(NA, -9L))

  ##fit
  SW({
  expect_warning(
    fit_DoseResponseCurve(LxTx[, c("Dose", "LxTx", "LxTx.Error")]),
    "Inf values found, replaced by NA")
  })

  ##all points have the same dose ... error but NULL
  data(ExampleData.LxTxData, envir = environment())
  tmp_LxTx <- LxTxData
  tmp_LxTx$Dose <- 10

  expect_message(expect_null(
      fit_DoseResponseCurve(tmp_LxTx)),
      "Error: All points have the same dose, NULL returned")

  ## check input objects ... matrix
  SW({
  expect_s4_class(
    fit_DoseResponseCurve(as.matrix(LxTxData)),
    class = "RLum.Results")
  })

  ## shuffle column names
  SW({
    LxTxData_shuffle <- LxTxData[,c("LxTx.Error", "LxTx", "Dose")]
    expect_s4_class(
      fit_DoseResponseCurve(LxTxData_shuffle),
      class = "RLum.Results")
  })

  ## test case for only two columns
  expect_s4_class(
    suppressWarnings(fit_DoseResponseCurve(LxTxData[,1:2], verbose = FALSE)),
    class = "RLum.Results")

  ## test case with all NA
  tmp_LxTx <- LxTxData
  tmp_LxTx$LxTx <- NA
  expect_message(expect_null(
      suppressWarnings(fit_DoseResponseCurve(tmp_LxTx))),
      "Error: After NA removal, nothing is left from the data set")

  ## test case without TnTx column
  tmp_LxTx <- LxTxData
  tmp_LxTx$TnTx <- NULL
  SW({
  expect_s4_class(
    fit_DoseResponseCurve(tmp_LxTx, verbose = FALSE),
    "RLum.Results")
  })

  ## do not include reg point
  expect_s4_class(
    fit_DoseResponseCurve(
      LxTxData,
      verbose = FALSE,
      fit.includingRepeatedRegPoints = FALSE),
    class = "RLum.Results")
})

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

  ## see https://github.com/R-Lum/Luminescence/pull/308
  skip_on_os("windows")
  skip_on_os("mac")

  snapshot.tolerance <- 1.5e-6

  SW({
  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "SSE",
      verbose = FALSE,
      n.MC = 10
    ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "SSE",
      fit.weights = NULL,
      verbose = FALSE,
      n.MC = 10
  ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
    LxTxData,
    fit.method = "SSE",
    fit.weights = 1 / LxTxData[[3]]^2,
    verbose = FALSE,
    n.MC = 10
  ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "LIN",
      verbose = FALSE,
      n.MC = 10
  ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "LIN",
      mode = "extrapolation",
      fit.force_through_origin = TRUE,
      verbose = FALSE,
      n.MC = 10
  ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "SSE+LIN",
      fit.bounds = FALSE,
      fit.force_through_origin = TRUE,
      verbose = FALSE,
      n.MC = 10
  ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "DSE",
      verbose = TRUE,
      txtProgressBar = FALSE,
      n.MC = 10
  ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "QDR",
      verbose = TRUE,
      txtProgressBar = FALSE,
      n.MC = 10
  ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "GOK",
      verbose = FALSE,
      n.MC = 10
  ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "OTOR",
      verbose = FALSE,
      n.MC = 10
  ), tolerance = snapshot.tolerance)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "QDR",
      mode = "extrapolation",
      verbose = TRUE,
      txtProgressBar = FALSE,
      n.MC = 10
  ), tolerance = 2.0e-5)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "QDR",
      mode = "extrapolation",
      fit.force_through_origin = TRUE,
      verbose = TRUE,
      txtProgressBar = FALSE,
      n.MC = 10
  ), tolerance = 5.0e-5)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "LIN",
      mode = "extrapolation",
      verbose = FALSE,
      n.MC = 10
  ), tolerance = 5.0e-5)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      LxTxData,
      fit.method = "GOK",
      mode = "extrapolation",
      verbose = FALSE,
      n.MC = 10
  ), tolerance = 5e-5)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      cbind(LxTxData, Test_Dose = 17),
      fit.method = "OTORX",
      fit.force_through_origin = TRUE,
      mode = "interpolation",
      verbose = TRUE,
      n.MC = 10
  ), tolerance = 5.0e-5)

  ## negative De
  expect_snapshot_RLum(fit_DoseResponseCurve(
      df_odd,
      mode = "interpolation",
      verbose = FALSE,
      n.MC = 10
  ), tolerance = 5.0e-3)

  expect_snapshot_RLum(fit_DoseResponseCurve(
      df_odd,
      mode = "extrapolation",
      verbose = FALSE,
      n.MC = 10
  ), tolerance = 4.0e-2)

  ## weights
  expect_snapshot_RLum(fit_DoseResponseCurve(
    LxTxData,
    mode = "interpolation",
    fit.weights = NULL,
    verbose = FALSE,
    n.MC = 10
  ), tolerance = 4.0e-2)

  expect_snapshot_RLum(fit_DoseResponseCurve(
    LxTxData,
    mode = "interpolation",
    fit.weights = "inverse_var",
    verbose = FALSE,
    n.MC = 10
  ), tolerance = 4.0e-2)

  expect_snapshot_RLum(fit_DoseResponseCurve(
    LxTxData,
    mode = "interpolation",
    fit.weights = "norm_inverse_std",
    verbose = FALSE,
    n.MC = 10
  ), tolerance = 4.0e-2)

  expect_snapshot_RLum(fit_DoseResponseCurve(
    LxTxData,
    mode = "interpolation",
    fit.weights = "inverse_std",
    verbose = FALSE,
    n.MC = 10
  ), tolerance = 4.0e-2)
  })

})

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

  ## self-call
  res <- fit_DoseResponseCurve(
      list(LxTxData, LxTxData),
      fit.method = "LIN",
      verbose = FALSE,
      n.MC = 10)
  expect_type(res, "list")
  expect_length(res, 2)

  set.seed(1)
  SW({
  expect_output(temp_SSE <-
    fit_DoseResponseCurve(
      LxTxData,
      fit.method = "SSE",
      verbose = TRUE,
      n.MC = 10
    ), " | D01 = ", fixed = TRUE)
  temp_LIN <-
    fit_DoseResponseCurve(
      LxTxData,
      fit.method = "LIN",
      verbose = FALSE,
      n.MC = 10
    )
  temp_SSELIN <-
    fit_DoseResponseCurve(
      LxTxData,
      fit.method = "SSE+LIN",
      verbose = FALSE,
      n.MC = 10
    )
  temp_DSE <-
    fit_DoseResponseCurve(
      LxTxData,
      fit.method = "DSE",
      verbose = FALSE,
      n.MC = 10
    )
  temp_QDR <-
    fit_DoseResponseCurve(
      LxTxData,
      fit.method = "QDR",
      verbose = FALSE,
      n.MC = 10
    )
  temp_GOK <-
    fit_DoseResponseCurve(
      LxTxData,
      fit.method = "GOK",
      verbose = FALSE,
      n.MC = 10)

  ## force through the origin
  temp_LxTx <-LxTxData
  temp_LxTx$LxTx[[7]] <- 1
  expect_s4_class(fit_DoseResponseCurve(
    temp_LxTx,
    fit.method = "GOK",
    verbose = FALSE,
    n.MC = 10,
    fit.force_through_origin = TRUE
  ), "RLum.Results")
temp_OTOR <-
  fit_DoseResponseCurve(
    LxTxData,
    fit.method = "OTOR",
    verbose = FALSE,
    n.MC = 10
  )
temp_OTORX <-
  fit_DoseResponseCurve(
    cbind(LxTxData, Test_Dose = 17),## we have to set the TEST_DOSE
    fit.method = "OTORX",
    fit.force_through_origin = TRUE,
    verbose = FALSE,
    n.MC = 10
  )

  ## TnTx column containing only NAs
  LxTxData$TnTx <- NA
  expect_s4_class(fit_DoseResponseCurve(LxTxData),
                  "RLum.Results")

## test reference dataset from
## (https://raw.githubusercontent.com/jll2/LumDRC/refs/heads/main/otorx.py)
LxTxData_alt <- data.frame(
  Dose = c(0, 17, 47, 94, 201, 402, 804, 1.07e+03, 1.68e+03, 3.35e+03, 5.36e+03,
           6.7e+03, 8.04e+03, 1e+04, 47, 94, 3.35e+03, 5.36e+03),
  LxTx = c(10, 1.09, 2.3, 3.67, 5.63, 8.09, 10.3, 11.4, 13, 14.8, 15.8, 16.9,
           17, 17.1, 2.26, 3.65, 15.6, 15.9),
  LxTx.Error = 0.1,
  TEST_DOSE = 17)
temp_OTORX_alt <-
  fit_DoseResponseCurve(
    cbind(LxTxData_alt, Test_Dose = 17),## we have to set the TEST_DOSE
    fit.method = "OTORX",
    fit.force_through_origin = TRUE,
    verbose = FALSE,
    n.MC = 10)

  temp_OTORX_alt2 <-
    fit_DoseResponseCurve(
      cbind(LxTxData_alt, Test_Dose = 17),## we have to set the TEST_DOSE
      fit.method = "OTORX",
      fit.force_through_origin = FALSE,
      verbose = FALSE,
      n.MC = 10)
  })

  expect_s3_class(temp_SSE$Fit, class = "nls")
  expect_s3_class(temp_LIN$Fit, class = "lm")
  expect_s3_class(temp_SSELIN$Fit, class = "nls")
  expect_s3_class(temp_DSE$Fit, class = "nls")
  expect_s3_class(temp_QDR$Fit, class = "lm")
  expect_s3_class(temp_GOK$Fit, class = "nls")
  expect_s3_class(temp_OTOR$Fit, class = "nls")
  expect_s3_class(temp_OTORX$Fit, class = "nls")
  expect_s3_class(temp_OTORX_alt$Fit, class = "nls")
  expect_s3_class(temp_OTORX_alt2$Fit, class = "nls")

   expect_equal(round(temp_SSE$De[[1]], digits = 2), 1737.71)
   expect_equal(round(sum(temp_SSE$De.MC, na.rm = TRUE), digits = 0), 17563)
   expect_equal(round(temp_LIN$De[[1]], digits = 1), 1673)
   expect_equal(round(sum(temp_LIN$De.MC, na.rm = TRUE), digits = 0),16983)
   expect_equal(round(temp_SSELIN$De[[1]], digits = 1), 1793)
   expect_equal(round(sum(temp_SSELIN$De.MC, na.rm = TRUE), digits = 0), 18068)
   expect_equal(round(temp_DSE$De[[1]], digits = 2), 1786.98)
   expect_equal(round(sum(temp_DSE$De.MC, na.rm = TRUE), digits = 0), 7303,
                tolerance = 10)
   expect_equal(round(temp_QDR$De[[1]], digits = 1), 1646.8)
   expect_equal(round(sum(temp_QDR$De.MC, na.rm = TRUE), digits = 0), 16342)
   expect_equal(round(temp_GOK$De[[1]], digits = 0), 1786)
   ##fix for different R versions
   if (R.version$major > "3"){
     if(any(grepl("aarch64", sessionInfo()$platform))) {
       expect_equal(round(sum(temp_GOK$De.MC, na.rm = TRUE), digits = 1), 17796,
                    tolerance = 0.1)

     } else {
       expect_equal(round(sum(temp_GOK$De.MC, na.rm = TRUE), digits = 1), 17828.9,
                    tolerance = 0.1)
     }
   }

   expect_equal(round(temp_OTOR$De[[1]], digits = 1),  1784.4)
   expect_equal(round(temp_OTORX$De[[1]], digits = 1), 2469.8, tolerance = 0.1)
   expect_equal(round(temp_OTORX_alt$De[[1]], digits = 2),  758.280)
   expect_equal(round(temp_OTORX_alt2$De[[1]], digits = 2),  793.21, tolerance = 0.2)
   expect_equal(round(sum(temp_OTOR$De.MC, na.rm = TRUE), digits = 0), 17611)
   expect_equal(round(sum(temp_OTORX$De.MC, na.rm = TRUE), digits = 0), 24408, tolerance = 0.2)

# Check extrapolation -----------------------------------------------------
  ## load data
  data(ExampleData.LxTxData, envir = environment())

  set.seed(1)
  LxTxData[1,2:3] <- c(0.5, 0.001)
  SW({
  LIN <- expect_s4_class(
    fit_DoseResponseCurve(LxTxData,mode = "extrapolation", fit.method = "LIN"),
    "RLum.Results")
  SSE <- expect_s4_class(
    fit_DoseResponseCurve(LxTxData,mode = "extrapolation", fit.method = "SSE"),
    "RLum.Results")
  SSELIN <- expect_s4_class(
    suppressWarnings(
      fit_DoseResponseCurve(
        LxTxData,mode = "extrapolation", fit.method = "SSE+LIN")),
    "RLum.Results")

  GOK <- expect_s4_class(
    fit_DoseResponseCurve(LxTxData,mode = "interpolation", fit.method = "GOK"),
    "RLum.Results")

  OTOR <- expect_s4_class(
    fit_DoseResponseCurve(LxTxData,mode = "extrapolation", fit.method = "OTOR"), "RLum.Results")

  ##OTORX
  OTORX <- expect_output(
    fit_DoseResponseCurve(
      object = cbind(LxTxData, Test_Dose = 17),
      mode = "extrapolation", fit.method = "OTORX"),
    "Fit failed for OTORX (extrapolation)", fixed = TRUE)

  ##OTORX ... trigger uniroot warning
  LxTxData[1,2:3] <- c(0.2, 0.001)
  expect_warning(
    OTORX <- fit_DoseResponseCurve(
      object = cbind(LxTxData, Test_Dose = 17),
      mode = "extrapolation", fit.method = "OTORX"),
    "Standard root estimation using stats::uniroot() failed", fixed = TRUE)
  })

  expect_equal(round(LIN$De$De,0), 184)
  expect_equal(round(SSE$De$De,0), 139)
  expect_equal(round(OTOR$De$De,0),  147)
  expect_equal(round(OTORX$De$De, 0), 1879)

  #it fails on some unix platforms for unknown reason.
  #expect_equivalent(round(SSELIN$De$De, 0), 110)

# Check alternate ---------------------------------------------------------
  ## load data
  data(ExampleData.LxTxData, envir = environment())

  set.seed(1)
  LxTxData[1,2:3] <- c(0.5, 0.001)

  SW({
  ##LIN
  expect_s4_class(
    fit_DoseResponseCurve(LxTxData, mode = "alternate", fit.method = "LIN"),
    "RLum.Results")

  ## SSE
  SSE <- expect_s4_class(
    fit_DoseResponseCurve(LxTxData, mode = "alternate", fit.method = "SSE"),
    "RLum.Results")

  ## SSE+LIN
  SSELIN <- expect_s4_class(
    suppressWarnings(
      fit_DoseResponseCurve(LxTxData, mode = "alternate", fit.method = "SSE+LIN")),
    "RLum.Results")
  })

  ## GOK
  expect_s4_class(
    fit_DoseResponseCurve(
      LxTxData,
      mode = "alternate",
      fit.method = "GOK",
      verbose = FALSE
    ),
    "RLum.Results"
  )

  ## OTOR
  expect_s4_class(
    fit_DoseResponseCurve(
      LxTxData,
      mode = "alternate",
      fit.method = "OTOR",
      verbose = FALSE
    ),
    "RLum.Results"
  )

  ## OTORX
  expect_s4_class(
    fit_DoseResponseCurve(
      cbind(LxTxData, Test_Dose = 17),
      mode = "alternate",
      fit.method = "OTORX",
      verbose = FALSE
    ),
    "RLum.Results"
  )

  ### check fit.weight option ------------
  expect_s4_class(
    fit_DoseResponseCurve(
      LxTxData,
      mode = "alternate",
      fit.method = "SSE",
      fit.weights = 1,
      verbose = FALSE
    ),
    "RLum.Results"
  )

  ## trigger OTOR related warning for
  ## extrapolation mode
  tmp <- structure(list(
    dose = c(
      0,
      1388.88888888889,
      2777.77777777778,
      4166.66666666667,
      8333.33333333333,
      16666.6666666667,
      33333.3333333333,
      0,
      1388.88888888889,
      2777.77777777778,
      4166.66666666667,
      8333.33333333333,
      16666.6666666667,
      33333.3333333333,
      0,
      1388.88888888889,
      2777.77777777778,
      4166.66666666667,
      8333.33333333333,
      16666.6666666667,
      33333.3333333333
    ),
    LxTx = c(
      1.54252220145258,
      4.43951568403849,
      6.23268064543138,
      7.84372723139206,
      12.1816246695694,
      16.220421545207,
      19.9805214420208,
      1.5693958789807,
      4.01446969642433,
      6.50442121919275,
      8.13912565845306,
      11.2791435536017,
      14.2739718127927,
      17.7646886436743,
      1.55083317135234,
      4.10327222363961,
      6.1705969614814,
      8.30005789933367,
      12.7612004529065,
      14.807776070804,
      17.1563663039162
    ),
    LxTx_X = c(
      0.130074482379272,
      2.59694106608334,
      4.46970034588506,
      3.0630786645803,
      0.744512263874143,
      6.0383153231303,
      0.785060450424326,
      3.16210365279,
      0.0425273193228004,
      2.9667194222907,
      0.187174353876429,
      4.29989597009486,
      4.19802308979151,
      2.77791088935002,
      0.248412040945932,
      0.626745230335262,
      3.80396486752602,
      16.1846310553925,
      4.14921514089229,
      1.40190110413806,
      7.74406545663656
    )
  ),
  class = "data.frame",
  row.names = c(NA, -21L))

  expect_warning(fit_DoseResponseCurve(
    tmp,
    mode = "extrapolation",
    fit.method = "OTOR",
    verbose = FALSE),
    "[fit_DoseResponseCurve()] Standard root estimation using stats::uniroot()",
    fixed = TRUE)

  ## only two valid points provided: this generates two warnings, hence
  ## we cannot use expect_warning(), which can only capture one at a time
  SW({
  warnings <- capture_warnings(expect_message(fit_DoseResponseCurve(
    data.frame(
        dose = c(0, 1388.88888888889, NA),
        LxTx = c(1.54252220145258, 4.43951568403849, NA),
        LxTx_X = c(0.130074482379272, 2.59694106608, NA)),
    verbose = TRUE),
    "'fit.method' changed to 'LIN'"))
  })
  expect_match(warnings, "1 NA values removed",
               all = FALSE, fixed = TRUE)
  expect_match(warnings, "Fitting a non-linear least-squares model requires",
               all = FALSE, fixed = TRUE)

  ## more coverage
  tmp$dose <- c(0:6, 100 + 0:6, 200000 + 0:6)
  expect_output(fit_DoseResponseCurve(
      tmp[4:8, ],
      fit.method = "GOK",
      fit.weights = "norm_inverse_std",
      verbose = TRUE,
      n.MC = 10),
      "Fit failed for GOK")
  expect_output(fit_DoseResponseCurve(
      cbind(tmp[4:8, ], Test_Dose = 17),
      fit.method = "OTORX",
      verbose = TRUE,
      n.MC = 10),
      "Fit failed for OTORX")
  expect_error(fit_DoseResponseCurve(
      LxTxData[-7, ],
      fit.method = "DSE",
      mode = "extrapolation"),
      "Mode 'extrapolation' for fitting method 'DSE' not supported",
      fixed = TRUE)

  set.seed(1)
  df <- data.frame(DOSE = c(0, 5, 10, 20, 25),
                   LxTx = c(40, -10, 30, -5, -20),
                   LxTx_X = c(2, 1, 1, 10, 1))
  expect_output(fit_DoseResponseCurve(df, fit.method = "SSE"),
                "De = NaN")
  expect_output(fit_DoseResponseCurve(rbind(df, c(3, 2, 3)), fit.method = "DSE"),
                "Fit failed for DSE")
})

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

  ## issue 374 --------------------------------------------------------------

  ## odd data that cause NaN but must not fail
  df <- data.frame(DOSE = c(0,5,10,20,30), LxTx = c(10,5,-20,-30,-40), LxTx_X = c(1, 1,1,1,1))
  SW({
  expect_s4_class(
    fit_DoseResponseCurve(df, fit.method = "SSE"), "RLum.Results")
  expect_s4_class(
    fit_DoseResponseCurve(df, fit.method = "SSE+LIN"), "RLum.Results")
  expect_s4_class(
    fit_DoseResponseCurve(df, fit.method = "DSE"), "RLum.Results")
  expect_s4_class(
    fit_DoseResponseCurve(df, fit.method = "OTOR"), "RLum.Results")
  })

  ## issue 723
  set.seed(1)
  df <- data.frame(DOSE = c(0, 5, 10, 20, 25),
                   LxTx = c(40, -10, 30, -5, -20),
                   LxTx_X = c(2, 2, 1, 2, 1))
  SW({
  expect_s4_class(fit_DoseResponseCurve(df, fit.method = "SSE"),
                  "RLum.Results")
  })

  ## issue 1539: test Berger's reference data -------------------------------

  QNL84_2_bleached <-
    read.table(system.file("extdata/QNL84_2_bleached.txt", package = "Luminescence"))
  QNL84_2_unbleached <-
    read.table(system.file("extdata/QNL84_2_unbleached.txt", package = "Luminescence"))
  STRB87_1_bleached <-
    read.table(system.file("extdata/STRB87_1_bleached.txt", package = "Luminescence"))
  STRB87_1_unbleached <-
    read.table(system.file("extdata/STRB87_1_unbleached.txt", package = "Luminescence"))

  ## add uncertainties of 2% to counts
  QNL84_2_bleached <- cbind(QNL84_2_bleached, QNL84_2_bleached[[2]] * 0.02)
  QNL84_2_unbleached <- cbind(QNL84_2_unbleached, QNL84_2_unbleached[[2]] * 0.02)
  STRB87_1_bleached <- cbind(STRB87_1_bleached, STRB87_1_bleached[[2]] * 0.02)
  STRB87_1_unbleached <- cbind(STRB87_1_unbleached, STRB87_1_unbleached[[2]] * 0.02)

  set.seed(1234)
  t_QNL84_2_bleached <- suppressWarnings(fit_DoseResponseCurve(
      QNL84_2_bleached,
      mode = "extrapolation",
      verbose = FALSE))
  t_QNL84_2_unbleached <- suppressWarnings(fit_DoseResponseCurve(
      QNL84_2_unbleached,
      mode = "extrapolation",
      verbose = FALSE))
  t_STRB87_1_bleached <- suppressWarnings(fit_DoseResponseCurve(
      STRB87_1_bleached,
      mode = "extrapolation",
      verbose = FALSE))
  t_STRB87_1_unbleached <- suppressWarnings(fit_DoseResponseCurve(
      STRB87_1_unbleached,
      mode = "extrapolation",
      verbose = FALSE))

  ## values are double-checked with
  ## Hayes, R.B., Haskell, E.H., Kenner, G.H., 1998. An assessment
  ## of the Levenberg-Marquardt fitting algorithm on saturating exponential
  ## data sets. Ancient TL 16, 57–62. https://doi.org/10.26034/la.atl.1998.294
  expect_equal(sum(t_QNL84_2_bleached$De[,c(1:2)]), expected = 204, tolerance = 0.01)
  expect_equal(sum(t_QNL84_2_unbleached$De[,c(1:2)]), expected = 126, tolerance = 0.01)
  expect_equal(sum(t_STRB87_1_bleached$De[,c(1:2)]), expected = 0.7, tolerance = 0.01)
  expect_equal(sum(t_STRB87_1_unbleached$De[,c(1:2)]), expected = 0.6, tolerance = 0.01)

  ## issue 1541
  expect_output(fit_DoseResponseCurve(df_odd, fit.method = "QDR"),
                "Fit:    QDR (interpolation) | De = 35.08", fixed = TRUE)

  ## issue 1543
  expect_output(fit_DoseResponseCurve(LxTxData, fit.method = "LIN", n.MC = 1),
                "Fit:    LIN (interpolation) | De = 1673.02", fixed = TRUE)
  expect_output(fit_DoseResponseCurve(LxTxData, fit.method = "QDR", n.MC = 1),
                "Fit:    QDR (interpolation) | De = 1646.83", fixed = TRUE)

  ## issue 1570
  df <- data.frame(Dose = c(0, 940.4, 2821.2, 4702, 0, 940.4),
                   LxTx = c(1.69, 91.72, 13.71, 16.92, -0.50, 4.34),
                   LxTx.Error = c(2.26, 1306.72, 9.12, 5.89, 0.37, 1.87),
                   TnTx = c(7.02, 0.45, 15.56, 29.38, 25.03, 19.07))
  expect_output(fit_DoseResponseCurve(df),
                "Fit:    SSE (interpolation) | De = 268.26 | D01 = 2612.50",
                fixed = TRUE)

  ## issue 1591
  data(ExampleData.BINfileData, envir = environment())
  object <- Risoe.BINfileData2RLum.Analysis(CWOSL.SAR.Data, pos=1)
  results <- analyse_SAR.CWOSL(
    object = object,
    signal_integral = 1:2,
    background_integral = 900:1000,
    log = "x",
    fit.method = "SSE",
    plot = FALSE,
    verbose = FALSE
  )

  t <- expect_s4_class(
    object = fit_DoseResponseCurve(results$LnLxTnTx.table, verbose = FALSE),
    class = "RLum.Results")
  expect_equal(results$data$De, t$De$De)

  ## issue 1634
  expect_warning(fit_DoseResponseCurve(LxTxData[1:5, ], fit.method = "DSE",
                                       verbose = FALSE),
                 "requires at least 5 dose points besides the natural, 'fit.method'")

  ## issue 1636
  expect_warning(fit_DoseResponseCurve(LxTxData[1:3, ], fit.method = "QDR",
                                       verbose = FALSE),
                 "requires at least 3 dose points besides the natural, 'fit.method'")

  ## issue 1730
  set.seed(2)
  fit <- fit_DoseResponseCurve(df_odd, mode = "interpolation",
                               verbose = FALSE, n.MC = 10)
  expect_equal(sum(is.na(fit$De.MC)), 6)
  expect_equal(fit$De$D01.ERROR, 1.369969, tolerance = 5.0e-3)
})

test_that("test internal functions", {
  testthat::skip_on_cran()

  ## reference tests from
  ## https://github.com/jll2/LumDRC/blob/main/otorx.py
  expect_equal(sum(Luminescence:::.nN2D(
    nN = 1-exp(-1),
    Q = c(0.5,0.6,0.7,0.8),
    D63 = c(1,10,100,1000))), expected = 1111)

  expect_equal(sum(Luminescence:::.nN2D(
    nN = 0.5,
    Q = c(0.5,0.6,0.7,0.8),
    D63 = c(1,10,100,1000))),
    expected = 661, tolerance = 2)

  expect_equal(sum(Luminescence:::.D2nN(
    D = 1,
    Q = c(-10,-3,0.1,1),
    D63 = 1)), expected = 2.5, tolerance = 1)

  expect_error(Luminescence:::.D2nN(D = 1, Q = c(-10, 0.1, 0, 0), D63 = 1),
               "Unsupported zero and non-zero Q")
})