#' @srrstats {G5.2} *Error and warning behaviour is explicitly demonstrated through tests.*
#' @srrstats {G5.2a} *Every error message is unique and tested.*
#' @srrstats {G5.2b} *Tests trigger every error message and compare with expected values.*
#' @srrstats {G5.3} *Return objects tested for absence of NA, NaN, Inf.*
#' @srrstats {G5.6} *Parameter recovery tests verify implementations produce expected results given data with known properties.*
#' @srrstats {G5.6a} *Parameter recovery tests succeed within defined tolerance rather than exact values.*
#' @srrstats {G5.6b} *Parameter recovery tests run with multiple random seeds when randomness is involved.*
#' @srrstats {G5.7} *Algorithm performance tests verify implementations perform correctly as data properties change.*
#' @srrstats {G5.8} *Edge condition tests verify appropriate behavior with extreme data properties.*
#' @srrstats {G5.8a} *Zero-length data tests trigger clear errors.*
#' @srrstats {G5.8b} *Unsupported data type tests trigger clear errors.*
#' @srrstats {G5.8c} *All-NA and all-identical data tests trigger clear errors or warnings.*
#' @srrstats {G5.8d} *Out-of-scope data tests verify appropriate behavior.*
#' @srrstats {G5.9} *Noise susceptibility tests verify stochastic behavior stability.*
#' @srrstats {G5.9a} *Trivial noise tests show results are stable at machine epsilon scale.*
#' @srrstats {G5.9b} *Random seed stability tests show consistent behavior across different seeds.*

# Sample data for testing
data <- data.frame(time = 1:10, completion = c(5, 15, 40, 60, 70, 75, 80, 85, 90, 95))

test_that("fit_sigmoidal returns a valid model for logistic type", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  expect_s3_class(fit, "nls")
})

test_that("fit_sigmoidal returns a valid model for pearl type", {
  fit <- fit_sigmoidal(data, "time", "completion", "pearl")
  expect_s3_class(fit, "nls")
})

test_that("fit_sigmoidal returns a valid model for gompertz type", {
  fit <- fit_sigmoidal(data, "time", "completion", "gompertz")
  expect_s3_class(fit, "nls")
})

test_that("fit_sigmoidal stops with an invalid model type", {
  expect_error(fit_sigmoidal(data, "time", "completion", "invalid_type"), "Invalid model type. Choose 'pearl', 'gompertz', or 'logistic'.")
})

test_that("predict_sigmoidal returns a data frame with predictions for logistic type", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  predictions <- predict_sigmoidal(fit, seq(min(data$time), max(data$time), length.out = 100), "logistic")
  expect_s3_class(predictions, "data.frame")
  expect_true(all(c("x", "pred") %in% names(predictions)))
})

test_that("predict_sigmoidal returns a data frame with predictions for pearl type", {
  fit <- fit_sigmoidal(data, "time", "completion", "pearl")
  predictions <- predict_sigmoidal(fit, seq(min(data$time), max(data$time), length.out = 100), "pearl")
  expect_s3_class(predictions, "data.frame")
  expect_true(all(c("x", "pred") %in% names(predictions)))
})

test_that("predict_sigmoidal returns a data frame with predictions for gompertz type", {
  fit <- fit_sigmoidal(data, "time", "completion", "gompertz")
  predictions <- predict_sigmoidal(fit, seq(min(data$time), max(data$time), length.out = 100), "gompertz")
  expect_s3_class(predictions, "data.frame")
  expect_true(all(c("x", "pred") %in% names(predictions)))
})

test_that("predict_sigmoidal stops with an invalid model type", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  expect_error(predict_sigmoidal(fit, seq(min(data$time), max(data$time), length.out = 100), "invalid_type"), "Invalid model type. Choose 'pearl', 'gompertz', or 'logistic'.")
})

# ============================================================================
# fit_sigmoidal Error Handling Tests
# ============================================================================
test_that("fit_sigmoidal validates column existence", {
  expect_error(
    fit_sigmoidal(data, "nonexistent", "completion", "logistic"),
    "Column nonexistent not found in data frame."
  )
  expect_error(
    fit_sigmoidal(data, "time", "nonexistent", "logistic"),
    "Column nonexistent not found in data frame."
  )
})

test_that("fit_sigmoidal returns correct coefficients for logistic", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  coefs <- coef(fit)
  expect_true("K" %in% names(coefs))
  expect_true("r" %in% names(coefs))
  expect_true("t0" %in% names(coefs))
})

test_that("fit_sigmoidal returns correct coefficients for gompertz", {
  fit <- fit_sigmoidal(data, "time", "completion", "gompertz")
  coefs <- coef(fit)
  expect_true("A" %in% names(coefs))
  expect_true("b" %in% names(coefs))
  expect_true("c" %in% names(coefs))
})

test_that("fit_sigmoidal returns correct coefficients for pearl", {
  fit <- fit_sigmoidal(data, "time", "completion", "pearl")
  coefs <- coef(fit)
  expect_true("K" %in% names(coefs))
  expect_true("r" %in% names(coefs))
  expect_true("t0" %in% names(coefs))
})

# ============================================================================
# predict_sigmoidal Error Handling Tests
# ============================================================================
test_that("predict_sigmoidal validates NULL fit", {
  expect_error(predict_sigmoidal(NULL, 1:10, "logistic"), "Fitted model is NULL.")
})

test_that("predict_sigmoidal validates numeric x_range", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  expect_error(predict_sigmoidal(fit, "not numeric", "logistic"), "x_range must be a numeric vector.")
})

test_that("predict_sigmoidal returns predictions with correct dimensions", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  x_range <- seq(1, 10, length.out = 50)
  predictions <- predict_sigmoidal(fit, x_range, "logistic")

  expect_equal(nrow(predictions), 50)
  expect_true("x" %in% names(predictions))
  expect_true("pred" %in% names(predictions))
})

# ============================================================================
# predict_sigmoidal Confidence Bounds Tests
# ============================================================================
test_that("predict_sigmoidal with confidence bounds returns lwr and upr columns", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  x_range <- seq(1, 10, length.out = 50)
  predictions <- predict_sigmoidal(fit, x_range, "logistic", conf_level = 0.95)

  expect_true("lwr" %in% names(predictions))
  expect_true("upr" %in% names(predictions))
  expect_true(all(predictions$lwr <= predictions$pred))
  expect_true(all(predictions$upr >= predictions$pred))
})

test_that("predict_sigmoidal validates conf_level range", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  x_range <- seq(1, 10, length.out = 50)

  expect_error(
    predict_sigmoidal(fit, x_range, "logistic", conf_level = 0),
    "conf_level must be between 0 and 1"
  )
  expect_error(
    predict_sigmoidal(fit, x_range, "logistic", conf_level = 1),
    "conf_level must be between 0 and 1"
  )
  expect_error(
    predict_sigmoidal(fit, x_range, "logistic", conf_level = 1.5),
    "conf_level must be between 0 and 1"
  )
  expect_error(
    predict_sigmoidal(fit, x_range, "logistic", conf_level = -0.5),
    "conf_level must be between 0 and 1"
  )
})

test_that("predict_sigmoidal confidence bounds work for pearl model", {
  fit <- fit_sigmoidal(data, "time", "completion", "pearl")
  x_range <- seq(1, 10, length.out = 50)
  predictions <- predict_sigmoidal(fit, x_range, "pearl", conf_level = 0.95)

  expect_true("lwr" %in% names(predictions))
  expect_true("upr" %in% names(predictions))
})

test_that("predict_sigmoidal confidence bounds work for gompertz model", {
  fit <- fit_sigmoidal(data, "time", "completion", "gompertz")
  x_range <- seq(1, 10, length.out = 50)
  predictions <- predict_sigmoidal(fit, x_range, "gompertz", conf_level = 0.95)

  expect_true("lwr" %in% names(predictions))
  expect_true("upr" %in% names(predictions))
})

test_that("predict_sigmoidal different confidence levels produce different bounds", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  x_range <- seq(1, 10, length.out = 50)

  pred_90 <- predict_sigmoidal(fit, x_range, "logistic", conf_level = 0.90)
  pred_95 <- predict_sigmoidal(fit, x_range, "logistic", conf_level = 0.95)
  pred_99 <- predict_sigmoidal(fit, x_range, "logistic", conf_level = 0.99)

  # Higher confidence level should produce wider bounds
  expect_true(all(pred_95$upr - pred_95$lwr >= pred_90$upr - pred_90$lwr - 1e-10))
  expect_true(all(pred_99$upr - pred_99$lwr >= pred_95$upr - pred_95$lwr - 1e-10))
})

# ============================================================================
# plot_sigmoidal Tests
# ============================================================================
test_that("plot_sigmoidal validates NULL fit", {
  expect_error(
    plot_sigmoidal(NULL, data, "time", "completion", "logistic"),
    "Fitted model is NULL."
  )
})

test_that("plot_sigmoidal validates column existence", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  expect_error(
    plot_sigmoidal(fit, data, "nonexistent", "completion", "logistic"),
    "Column nonexistent not found in data frame."
  )
  expect_error(
    plot_sigmoidal(fit, data, "time", "nonexistent", "logistic"),
    "Column nonexistent not found in data frame."
  )
})

test_that("plot_sigmoidal returns predictions invisibly", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  result <- plot_sigmoidal(fit, data, "time", "completion", "logistic")

  expect_s3_class(result, "data.frame")
  expect_true("x" %in% names(result))
  expect_true("pred" %in% names(result))
})

test_that("plot_sigmoidal returns confidence bounds when requested", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  result <- plot_sigmoidal(fit, data, "time", "completion", "logistic", conf_level = 0.95)

  expect_true("lwr" %in% names(result))
  expect_true("upr" %in% names(result))
})

test_that("plot_sigmoidal works with custom parameters", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")

  # Should not error with custom parameters
  result <- plot_sigmoidal(fit, data, "time", "completion", "logistic",
    n_points = 50,
    main = "Test Plot",
    xlab = "X Label",
    ylab = "Y Label",
    line_col = "blue",
    ci_col = "gray"
  )

  expect_s3_class(result, "data.frame")
  expect_equal(nrow(result), 50)
})

test_that("plot_sigmoidal works with all model types", {
  # Logistic
  fit_log <- fit_sigmoidal(data, "time", "completion", "logistic")
  result_log <- plot_sigmoidal(fit_log, data, "time", "completion", "logistic")
  expect_s3_class(result_log, "data.frame")

  # Pearl
  fit_pearl <- fit_sigmoidal(data, "time", "completion", "pearl")
  result_pearl <- plot_sigmoidal(fit_pearl, data, "time", "completion", "pearl")
  expect_s3_class(result_pearl, "data.frame")

  # Gompertz
  fit_gomp <- fit_sigmoidal(data, "time", "completion", "gompertz")
  result_gomp <- plot_sigmoidal(fit_gomp, data, "time", "completion", "gompertz")
  expect_s3_class(result_gomp, "data.frame")
})

# ============================================================================
# Edge Cases and Integration Tests
# ============================================================================
test_that("sigmoidal functions work with minimal data", {
  minimal_data <- data.frame(time = 1:5, completion = c(10, 30, 50, 70, 90))
  fit <- fit_sigmoidal(minimal_data, "time", "completion", "logistic")
  expect_s3_class(fit, "nls")

  predictions <- predict_sigmoidal(fit, 1:5, "logistic")
  expect_equal(nrow(predictions), 5)
})

test_that("sigmoidal predictions are monotonically reasonable", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  x_range <- seq(1, 10, length.out = 100)
  predictions <- predict_sigmoidal(fit, x_range, "logistic")

  # For logistic/sigmoidal curves, predictions should generally increase
  # (allowing for some numerical tolerance)
  diffs <- diff(predictions$pred)
  expect_true(sum(diffs >= -0.1) > length(diffs) * 0.9)
})

test_that("fitted values are close to original data", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  predictions <- predict_sigmoidal(fit, data$time, "logistic")

  # Predictions should be reasonably close to actual values
  residuals <- abs(predictions$pred - data$completion)
  expect_true(mean(residuals) < 10) # Average error < 10
})

# ============================================================================
# NaN/NA/Inf Error Tests (G5.2, G5.2b)
# ============================================================================
test_that("fit_sigmoidal rejects NaN in data columns", {
  bad_data <- data.frame(time = c(1, NaN, 3), completion = c(5, 15, 40))
  expect_error(
    fit_sigmoidal(bad_data, "time", "completion", "logistic"),
    "Data columns must not contain NaN values."
  )
})

test_that("fit_sigmoidal rejects NA in data columns", {
  bad_data <- data.frame(time = c(1, NA, 3), completion = c(5, 15, 40))
  expect_error(
    fit_sigmoidal(bad_data, "time", "completion", "logistic"),
    "Data columns must not contain NA values."
  )
})

test_that("fit_sigmoidal rejects Inf in data columns", {
  bad_data <- data.frame(time = c(1, Inf, 3), completion = c(5, 15, 40))
  expect_error(
    fit_sigmoidal(bad_data, "time", "completion", "logistic"),
    "Data columns must not contain infinite values."
  )
})

test_that("predict_sigmoidal rejects NaN in x_range", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  expect_error(
    predict_sigmoidal(fit, c(1, NaN, 3), "logistic"),
    "x_range must not contain NaN values."
  )
})

test_that("predict_sigmoidal rejects NA in x_range", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  expect_error(
    predict_sigmoidal(fit, c(1, NA, 3), "logistic"),
    "x_range must not contain NA values."
  )
})

test_that("predict_sigmoidal rejects Inf in x_range", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  expect_error(
    predict_sigmoidal(fit, c(1, Inf, 3), "logistic"),
    "x_range must not contain infinite values."
  )
})

test_that("plot_sigmoidal rejects NaN in data columns for plotting", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  bad_data <- data.frame(time = c(1, NaN, 3), completion = c(5, 15, 40))
  expect_error(
    plot_sigmoidal(fit, bad_data, "time", "completion", "logistic"),
    "Data columns must not contain NaN values for plotting."
  )
})

test_that("plot_sigmoidal rejects NA in data columns for plotting", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  bad_data <- data.frame(time = c(1, NA, 3), completion = c(5, 15, 40))
  expect_error(
    plot_sigmoidal(fit, bad_data, "time", "completion", "logistic"),
    "Data columns must not contain NA values for plotting."
  )
})

test_that("plot_sigmoidal rejects Inf in data columns for plotting", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  bad_data <- data.frame(time = c(1, Inf, 3), completion = c(5, 15, 40))
  expect_error(
    plot_sigmoidal(fit, bad_data, "time", "completion", "logistic"),
    "Data columns must not contain infinite values for plotting."
  )
})

# ============================================================================
# G5.3: Return value tests
# ============================================================================
test_that("predict_sigmoidal result contains no NA, NaN, or Inf", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  predictions <- predict_sigmoidal(fit, seq(1, 10, length.out = 50), "logistic")
  expect_false(anyNA(predictions$x))
  expect_false(any(is.nan(predictions$x)))
  expect_false(any(is.infinite(predictions$x)))
  expect_false(anyNA(predictions$pred))
  expect_false(any(is.nan(predictions$pred)))
  expect_false(any(is.infinite(predictions$pred)))
})

# ============================================================================
# Parameter Recovery Tests (G5.6, G5.6a, G5.6b)
# ============================================================================

test_that("fit_sigmoidal recovers logistic parameters from synthetic data (seed 123)", {
  set.seed(123)

  # Generate synthetic data from known logistic curve
  K_true <- 100
  r_true <- 0.5
  t0_true <- 10

  x <- seq(1, 20, length.out = 50)
  y_true <- K_true / (1 + exp(-r_true * (x - t0_true)))
  # Add small noise
  y <- y_true + rnorm(50, mean = 0, sd = 1)

  synthetic_data <- data.frame(time = x, completion = y)
  fit <- fit_sigmoidal(synthetic_data, "time", "completion", "logistic")
  coefs <- coef(fit)

  # Recovered parameters should be close to true values
  expect_equal(as.numeric(coefs["K"]), K_true, tolerance = 5)
  expect_equal(as.numeric(coefs["r"]), r_true, tolerance = 0.1)
  expect_equal(as.numeric(coefs["t0"]), t0_true, tolerance = 2)
})

test_that("fit_sigmoidal recovers logistic parameters from synthetic data (seed 42)", {
  set.seed(42)

  # Same test with different seed to verify robustness (G5.6b)
  K_true <- 100
  r_true <- 0.5
  t0_true <- 10

  x <- seq(1, 20, length.out = 50)
  y_true <- K_true / (1 + exp(-r_true * (x - t0_true)))
  y <- y_true + rnorm(50, mean = 0, sd = 1)

  synthetic_data <- data.frame(time = x, completion = y)
  fit <- fit_sigmoidal(synthetic_data, "time", "completion", "logistic")
  coefs <- coef(fit)

  expect_equal(as.numeric(coefs["K"]), K_true, tolerance = 5)
  expect_equal(as.numeric(coefs["r"]), r_true, tolerance = 0.1)
  expect_equal(as.numeric(coefs["t0"]), t0_true, tolerance = 2)
})

test_that("fit_sigmoidal recovers gompertz parameters from synthetic data (seed 123)", {
  set.seed(123)

  # Generate synthetic Gompertz data
  A_true <- 100
  b_true <- 2
  c_true <- 0.3

  x <- seq(1, 20, length.out = 50)
  y_true <- A_true * exp(-b_true * exp(-c_true * x))
  y <- y_true + rnorm(50, mean = 0, sd = 1)

  synthetic_data <- data.frame(time = x, completion = y)
  fit <- fit_sigmoidal(synthetic_data, "time", "completion", "gompertz")
  coefs <- coef(fit)

  expect_equal(as.numeric(coefs["A"]), A_true, tolerance = 5)
  expect_equal(as.numeric(coefs["b"]), b_true, tolerance = 0.5)
  expect_equal(as.numeric(coefs["c"]), c_true, tolerance = 0.1)
})

test_that("fit_sigmoidal recovers gompertz parameters from synthetic data (seed 42)", {
  set.seed(42)

  # Same test with different seed (G5.6b)
  A_true <- 100
  b_true <- 2
  c_true <- 0.3

  x <- seq(1, 20, length.out = 50)
  y_true <- A_true * exp(-b_true * exp(-c_true * x))
  y <- y_true + rnorm(50, mean = 0, sd = 1)

  synthetic_data <- data.frame(time = x, completion = y)
  fit <- fit_sigmoidal(synthetic_data, "time", "completion", "gompertz")
  coefs <- coef(fit)

  expect_equal(as.numeric(coefs["A"]), A_true, tolerance = 5)
  expect_equal(as.numeric(coefs["b"]), b_true, tolerance = 0.5)
  expect_equal(as.numeric(coefs["c"]), c_true, tolerance = 0.1)
})

# ============================================================================
# Algorithm Performance Tests (G5.7)
# ============================================================================

test_that("fit_sigmoidal improves with more data points", {
  set.seed(123)

  K_true <- 100
  r_true <- 0.5
  t0_true <- 10

  # Test with different numbers of data points
  for (n_points in c(10, 30, 50)) {
    x <- seq(1, 20, length.out = n_points)
    y_true <- K_true / (1 + exp(-r_true * (x - t0_true)))
    y <- y_true + rnorm(n_points, mean = 0, sd = 1)

    data_test <- data.frame(time = x, completion = y)
    fit <- fit_sigmoidal(data_test, "time", "completion", "logistic")
    coefs <- coef(fit)

    # Store error for comparison
    error_K <- abs(as.numeric(coefs["K"]) - K_true)

    # More data should generally give better fits (with some randomness)
    expect_true(is.numeric(error_K))
    expect_true(error_K < 20) # Should be reasonably close
  }
})

test_that("fit_sigmoidal convergence improves with less noisy data", {
  set.seed(42)

  K_true <- 100
  r_true <- 0.5
  t0_true <- 10

  x <- seq(1, 20, length.out = 50)
  y_true <- K_true / (1 + exp(-r_true * (x - t0_true)))

  # Test with different noise levels
  noise_levels <- c(5, 2, 0.5)
  errors_K <- numeric(length(noise_levels))

  for (i in seq_along(noise_levels)) {
    y <- y_true + rnorm(50, mean = 0, sd = noise_levels[i])
    data_test <- data.frame(time = x, completion = y)
    fit <- fit_sigmoidal(data_test, "time", "completion", "logistic")
    coefs <- coef(fit)
    errors_K[i] <- abs(as.numeric(coefs["K"]) - K_true)
  }

  # Lower noise should give better parameter estimates
  expect_true(errors_K[3] < errors_K[1])
})

# ============================================================================
# Edge Condition Tests (G5.8c) - All-Identical Data
# ============================================================================

test_that("fit_sigmoidal fails gracefully with all identical y values", {
  data_constant <- data.frame(time = 1:10, completion = rep(50, 10))

  # Should error or fail to converge since there's no curve to fit
  # Platform-dependent: nlsLM may error or succeed with degenerate fit
  result <- tryCatch(
    {
      fit <- fit_sigmoidal(data_constant, "time", "completion", "logistic")
      # If it doesn't error, check that it's a degenerate fit
      # (constant data means fit should have issues)
      list(success = TRUE, fit = fit)
    },
    error = function(e) {
      list(success = FALSE, error = e)
    }
  )

  # Test passes if either:
  # 1. An error was thrown (success = FALSE), OR
  # 2. The fit succeeded but is degenerate (all predictions ≈ constant)
  if (result$success) {
    # Check that the fit is degenerate (predictions are essentially constant)
    preds <- predict(result$fit, newdata = data_constant)
    # All predictions should be very close to each other (constant)
    expect_true(sd(preds) < 1e-6)
  } else {
    # An error was thrown - this is expected behavior
    expect_true(!result$success)
  }
})

test_that("fit_sigmoidal handles all NA in completion column", {
  data_na <- data.frame(time = 1:10, completion = rep(NA_real_, 10))

  expect_error(
    fit_sigmoidal(data_na, "time", "completion", "logistic"),
    "Data columns must not contain NA values"
  )
})

# ============================================================================
# Noise Susceptibility Tests (G5.9b) - Random Seed Stability
# ============================================================================

test_that("fit_sigmoidal produces consistent results across seeds", {
  K_true <- 100
  r_true <- 0.5
  t0_true <- 10

  x <- seq(1, 20, length.out = 50)
  y_true <- K_true / (1 + exp(-r_true * (x - t0_true)))

  # Fit with different noise realizations
  set.seed(111)
  y1 <- y_true + rnorm(50, 0, 1)
  data1 <- data.frame(time = x, completion = y1)
  fit1 <- fit_sigmoidal(data1, "time", "completion", "logistic")

  set.seed(222)
  y2 <- y_true + rnorm(50, 0, 1)
  data2 <- data.frame(time = x, completion = y2)
  fit2 <- fit_sigmoidal(data2, "time", "completion", "logistic")

  set.seed(333)
  y3 <- y_true + rnorm(50, 0, 1)
  data3 <- data.frame(time = x, completion = y3)
  fit3 <- fit_sigmoidal(data3, "time", "completion", "logistic")

  # Parameter estimates should be reasonably consistent
  coef1 <- coef(fit1)
  coef2 <- coef(fit2)
  coef3 <- coef(fit3)

  # All K estimates should be within reasonable range
  expect_true(abs(as.numeric(coef1["K"]) - as.numeric(coef2["K"])) < 10)
  expect_true(abs(as.numeric(coef2["K"]) - as.numeric(coef3["K"])) < 10)
})

# ============================================================================
# print.pra_sigmoidal_fit
# ============================================================================

test_that("print.pra_sigmoidal_fit produces output and returns invisibly", {
  fit <- fit_sigmoidal(data, "time", "completion", "logistic")
  expect_s3_class(fit, "pra_sigmoidal_fit")
  out <- capture.output(result <- print(fit))
  expect_true(any(grepl("Sigmoidal Model Fit Summary", out)))
  # Should include coefficient information from summary
  expect_true(any(grepl("K|r|t0", out)))
})

test_that("fit_sigmoidal assigns pra_sigmoidal_fit class", {
  fit <- fit_sigmoidal(data, "time", "completion", "pearl")
  expect_s3_class(fit, "pra_sigmoidal_fit")
  expect_s3_class(fit, "nls")
  # Class should have pra_sigmoidal_fit as first element
  expect_equal(class(fit)[1], "pra_sigmoidal_fit")
})