test_that("Regularization works", {
  skip_on_cran()
  skip_on_ci()
  local_edition(2)
  d <- subset(pcod, year >= 2015)
  d$depth_scaled <- as.numeric(scale(d$depth_scaled))
  m1 <- sdmTMB(data = d,
    formula = density ~ 0 + depth_scaled + as.factor(year),
    mesh = make_mesh(d, c("X", "Y"), n_knots = 50, type = "kmeans"),
    family = tweedie(link = "log"))

  # Bypassing via NAs:
  m2 <- sdmTMB(data = d,
    formula = density ~ 0 + depth_scaled + as.factor(year),
    mesh = make_mesh(d, c("X", "Y"), n_knots = 50, type = "kmeans"),
    family = tweedie(link = "log"),
    priors = sdmTMBpriors(b = normal(c(NA, NA, NA), c(NA, NA, NA))))
  expect_equal(m1$sd_report, m2$sd_report)

  # Ridge regression on depth term:
  m2 <- sdmTMB(data = d,
    formula = density ~ 0 + depth_scaled + as.factor(year),
    mesh = make_mesh(d, c("X", "Y"), n_knots = 50, type = "kmeans"),
    family = tweedie(link = "log"),
    priors = sdmTMBpriors(b = normal(c(0, NA, NA), c(1, NA, NA))))
  b1 <- tidy(m1)
  b2 <- tidy(m2)
  expect_lt(b2$estimate[1], b1$estimate[2])
  expect_lt(b2$std.error[1], b1$std.error[2])
})

test_that("A time-varying model fits and predicts appropriately", {
  skip_on_cran()
  skip_on_ci()
  local_edition(2)
  # SEED <- 42
  # set.seed(SEED)
  # x <- stats::runif(60, -1, 1)
  # y <- stats::runif(60, -1, 1)
  # initial_betas <- 0.5
  # range <- 0.5
  # sigma_O <- 0
  # sigma_E <- 0.1
  # phi <- 0.1
  # sigma_V <- 0.3
  # loc <- data.frame(x = x, y = y)
  # spde <- make_mesh(loc, c("x", "y"), cutoff = 0.02)
  #
  # s <- sdmTMB_sim(
  #   x = x, y = y, mesh = spde, range = range,
  #   betas = initial_betas, time_steps = 12L, sigma_V = sigma_V,
  #   phi = phi, sigma_O = sigma_O, sigma_E = sigma_E,
  #   seed = SEED
  # )
  spde <- make_mesh(pcod, c("X", "Y"), cutoff = 20)
  m <- sdmTMB(data = pcod, formula = density ~ 0, spatial = TRUE,
    time_varying = ~ 1, time = "year", mesh = spde, family = tweedie(),
    spatiotemporal = "off")
  expect_equal(exp(m$model$par["ln_tau_V"])[[1]], 0.5971512, tolerance = 0.001)
  tidy(m, effects = "ran_par")
  # b_t <- dplyr::group_by(s, time) %>%
  #   dplyr::summarize(b_t = unique(b), .groups = "drop") %>%
  #   dplyr::pull(b_t)
  # r <- m$tmb_obj$report()
  # b_t_fit <- r$b_rw_t[,,1]
  # plot(b_t, b_t_fit, asp = 1);abline(a = 0, b = 1)
  # expect_equal(mean((b_t- b_t_fit)^2), 0, tolerance = 1e-4)
  p <- predict(m, newdata = NULL)
  # plot(p$est, s$observed, asp = 1);abline(a = 0, b = 1)
  # expect_equal(mean((p$est - s$observed)^2), 0, tolerance = 0.01)

  cols <- c("est", "est_non_rf", "est_rf", "omega_s")
  p_nd <- predict(m, newdata = pcod)
  expect_equal(p[,cols], p_nd[,cols], tolerance = 1e-4)
})

test_that("Year indexes get created correctly", {
  expect_identical(make_year_i(c(1, 2, 3)),    c(0L, 1L, 2L))
  expect_identical(make_year_i(c(1L, 2L, 3L)), c(0L, 1L, 2L))
  expect_identical(make_year_i(c(1L, 2L, 4L)), c(0L, 1L, 2L))
  expect_identical(make_year_i(c(1, 2, 4, 2)), c(0L, 1L, 2L, 1L))
  expect_identical(make_year_i(c(1, 4, 2)),    c(0L, 2L, 1L))
})

test_that("A spatially varying coefficient model fits", {
  skip_on_cran()
  skip_on_ci()
  d <- subset(pcod, year >= 2011) # subset for speed
  pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 30)
  d$scaled_year <- (d$year - mean(d$year)) / sd(d$year)
  m <- sdmTMB(density ~ depth_scaled, data = d,
    mesh = pcod_spde, family = tweedie(link = "log"),
    spatial_varying = ~ 0 + scaled_year, time = "year")
  expect_true(all(!is.na(summary(m$sd_report)[,"Std. Error"])))
  nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
  nd$scaled_year <- scale(nd$year)
  p <- predict(m, newdata = subset(nd, year >= 2011))
})


test_that("SPDE as generated by make_mesh is consistent", {
  set.seed(42)
  skip_on_cran()
  skip_on_ci()
  set.seed(1)
  local_edition(2)
  x <- stats::runif(70, -1, 1)
  y <- stats::runif(70, -1, 1)
  loc <- data.frame(x = x, y = y)
  spde <- make_mesh(loc, c("x", "y"), cutoff = 0.02)

  # spot check first 10 locations
  loc_xy <- matrix(c(
    -0.468982674,-0.321854124,
    -0.255752201,0.678880700,
    0.145706727,-0.306633022,
    0.816415580,-0.332450138,
    -0.596636138,-0.047297510,
    0.796779370,0.784396672,
    0.889350537,0.728678941,
    0.321595585,-0.220020913,
    0.258228088,0.554641398,
    -0.876427459,0.921235994), ncol = 2, byrow = TRUE)

  expect_equal(c(spde$loc_xy[1:10,]), c(loc_xy[1:10,]), tolerance = 1e-8)

  # Depends on stable or testing INLA!?
  # test the n
  # expect_equal(spde$spde$param.inla$n, 141.0, tolerance = 1e-8)
  # expect_equal(spde$spde$param.inla$n, 150)

  # test M0 values - first 10
  # expect_equal(c(0.06889138, 0.04169878, 0.01547313, 0.05734409, 0.03271476, 0.03672571, 0.05200510, 0.01293629, 0.03631651, 0.04446844),
  #   as.numeric(spde$spde$param.inla$M0)[c(1,143,285,427,569,711,853,995,1137,1279)], tolerance=1e-4)
  # # test M1 values - first 10
  # expect_equal(c(3.602458, 5.105182, 4.571970, 4.353246, 4.496775, 4.802626, 4.505773, 4.935963, 4.419817, 5.169096),
  #   as.numeric(spde$spde$param.inla$M1)[c(1,143,285,427,569,711,853,995,1137,1279)], tolerance=1e-4)
  # # test M2 values - first 10
  # expect_equal(c(253.1738,814.3858,1754.2543,401.3752,865.0118,754.9644,543.8968,2425.8404,729.7511,1014.0049),
  #   as.numeric(spde$spde$param.inla$M2)[c(1,143,285,427,569,711,853,995,1137,1279)], tolerance=1e-4)

})

# test_that("The `map` argument works.", {
#   skip_on_cran()
#   skip_on_ci()
#   d <- subset(pcod, year >= 2013)
#   pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 50)
#
#   m <- sdmTMB(density ~ 0 + as.factor(year),
#     data = d, time = "year", mesh = pcod_spde,
#     family = tweedie(link = "log")
#   )
#   p <- predict(m)
#
#   .map <- m$tmb_map
#   .map$ln_tau_O <- as.factor(NA)
#   .map$ln_tau_E <- as.factor(NA)
#   .map$omega_s <- factor(rep(NA, length(m$tmb_params$omega_s)))
#   .map$epsilon_st <- factor(rep(NA, length(m$tmb_params$epsilon_st)))
#   .map$ln_kappa <- as.factor(NA)
#
#   m.map <- sdmTMB(density ~ 0 + as.factor(year),
#     data = d, time = "year", mesh = pcod_spde,
#     family = tweedie(link = "log"), map = .map
#   )
#   p.map <- predict(m.map)
#
#   expect_true(all(p.map$omega_s == 0))
#   expect_true(all(p.map$epsilon_st == 0))
#   expect_true(!identical(p$est, p.map$est))
#   })

test_that("profile option works", {
  skip_on_cran()
  skip_on_ci()
  d <- subset(pcod, year == 2013)
  pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 10)
  mp <- sdmTMB(density ~ depth_scaled + depth_scaled2,
    data = d, mesh = pcod_spde, control = sdmTMBcontrol(profile = TRUE),
    family = tweedie(link = "log")
  )
  m <- sdmTMB(density ~ depth_scaled + depth_scaled2,
    data = d, mesh = pcod_spde, control = sdmTMBcontrol(profile = FALSE),
    family = tweedie(link = "log")
  )
  expect_lt(length(mp$model$par), length(m$model$par))
  bp <- tidy(mp)
  b <- tidy(m)
  expect_equal(bp$estimate, b$estimate, tolerance = 0.05)
})

test_that("The mapping off spatial and spatiotemporal fields works.", {
  skip_on_cran()
  skip_on_ci()
  d <- subset(pcod, year >= 2013)
  pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 50)
  m <- sdmTMB(density ~ 0 + as.factor(year),
    data = d, time = "year", mesh = pcod_spde,
    family = tweedie(link = "log")
  )
  p <- predict(m)
  m.map <- sdmTMB(density ~ 0 + as.factor(year),
    data = d, time = "year", mesh = pcod_spde,
    family = tweedie(link = "log"), spatial = "off", spatiotemporal = "off"
  )
  p.map <- predict(m.map)

  expect_true(is.na(m.map$tmb_map$ln_tau_E))
  expect_true(is.na(m.map$tmb_map$ln_tau_O))
  expect_true(!identical(p$est, p.map$est))
  expect_true(length(unique(p.map$est)) == length(unique(d$year)))

  # basic lm test:
  dpos <- d[d$density > 0, ]
  pcod_spde <- make_mesh(dpos, c("X", "Y"), cutoff = 50)
  m.sdmTMB.map <- sdmTMB(log(density) ~ depth, data = dpos,
    family = gaussian(), spatial = "off", spatiotemporal = "off",
    mesh = pcod_spde)
  m.stats.glm <- stats::lm(log(density) ~ depth, data = dpos)
  m.glmmTMB <- glmmTMB::glmmTMB(log(density) ~ depth, data = dpos)
  .t <- tidy(m.sdmTMB.map)
  expect_equal(.t$estimate, as.numeric(coef(m.stats.glm)), tolerance = 1e-5)
  expect_equal(exp(m.sdmTMB.map$model$par[["ln_phi"]]),
    glmmTMB::sigma(m.glmmTMB), tolerance = 1e-5)

  # Bernoulli:
  pcod_binom <- pcod_2011
  pcod_binom$present <- ifelse(pcod_binom$density > 0, 1L, 0L)
  pcod_spde <- make_mesh(pcod_binom, c("X", "Y"), cutoff = 50)
  m.sdmTMB.map <- sdmTMB(present ~ depth, data = pcod_binom,
    family = binomial(link = "logit"), spatial = "off", spatiotemporal = "off", mesh = pcod_spde)
  m.stats.glm <- stats::glm(present ~ depth, data = pcod_binom,
    family = binomial(link = "logit"))
  m.glmmTMB <- glmmTMB::glmmTMB(present ~ depth, data = pcod_binom,
    family = binomial(link = "logit"))
  b1 <- as.numeric(unclass(glmmTMB::fixef(m.glmmTMB))$cond)
  b2 <- tidy(m.sdmTMB.map)$estimate
  b3 <- as.numeric(coef(m.stats.glm))
  expect_equal(b1, b2, tolerance = 1e-6)
  expect_equal(b3, b2, tolerance = 1e-6)
})

test_that("Large coordinates cause a warning.", {
  skip_on_cran()
  d <- subset(pcod, year == 2017)
  d$X <- d$X * 1000
  d$Y <- d$Y * 1000
  expect_warning(pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 70))
})

test_that("Random walk fields work", {
  skip_on_cran()
  skip_on_ci()
  d <- subset(pcod, year >= 2013)
  pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 15)
  m_rw <- sdmTMB(density ~ 1,
    data = d, time = "year", mesh = pcod_spde,
    control = sdmTMBcontrol(newton_loops = 1),
    family = tweedie(link = "log"), spatiotemporal = "RW"
  )
  expect_identical(m_rw$tmb_data$rw_fields, TRUE)
  p_rw <- predict(m_rw)

  # close to AR1 with high rho:
  m_ar1 <- sdmTMB(density ~ 1,
    data = d, time = "year", mesh = pcod_spde,
    family = tweedie(link = "log"), spatiotemporal = "AR1",
    control = sdmTMBcontrol(
      newton_loops = 1,
      start = list(ar1_phi = qlogis((0.98 + 1) / 2)),
      map = list(ar1_phi = factor(NA)))
  )
  expect_identical(m_ar1$tmb_data$ar1_fields, TRUE)
  expect_identical(m_ar1$tmb_data$rw_fields, FALSE)
  p_ar1 <- predict(m_ar1)

  expect_gt(stats::cor(p_ar1$est, p_rw$est), 0.9)
})


test_that("start works", {
  skip_on_cran()
  skip_on_ci()

    expect_error({
      m2 <- sdmTMB(density ~ poly(depth, 2),
        data = pcod_2011,
        mesh = pcod_mesh_2011, family = tweedie(),
        control = sdmTMBcontrol(start = list(ln_kappa = -1.78)))
    }, regexp = "kappa")

  expect_message({
    m2 <- sdmTMB(density ~ poly(depth, 2),
      data = pcod_2011,
      mesh = pcod_mesh_2011, family = tweedie(),
      control = sdmTMBcontrol(start = list(ln_kappa = matrix(c(-1.78, -1.78), ncol = 1L))))
  }, regexp = "ln_kappa")
})

test_that("Multiple SVC works", {
  skip_on_cran()
  skip_on_ci()

  mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 15)

  pcod$syear <- as.numeric(scale(pcod$year))
  fit <- sdmTMB(
    density ~ 1,
    spatial_varying = ~ 0 + syear + depth_scaled,
    mesh = mesh,
    family = delta_gamma(),
    data = pcod
  )
  fit$sd_report
  s <- as.list(fit$sd_report, "Std. Error")
  expect_true(sum(is.na(s$b_j)) == 0L)
  fit
  nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
  nd$syear <- as.numeric(scale(nd$year))
  p <- predict(fit, newdata = nd)
  # p <- predict(fit, newdata = NULL)
})

test_that("More esoteric prediction options work", {
  skip_on_cran()
  skip_on_ci()
  fit <- sdmTMB(
    density ~ depth_scaled,
    data = pcod_2011, mesh = pcod_mesh_2011,
    family = delta_gamma()
  )
  expect_error(p <- predict(fit, nsim = 5, model = "aaa"),
    regexp = "model"
  )
  expect_error(p <- predict(fit, nsim = 5, model = 99),
    regexp = "model"
  )

  p <- predict(fit, nsim = 5, model = NA)
  head(p)
  p <- predict(fit, nsim = 5, model = 1)
  head(p)
  # p <- predict(fit, nsim = 5, model = 1, type = "response")
  # head(p)
  # expect_true(all(p <= 1 & p >= 0))
  p <- predict(fit, nsim = 5, model = 2)
  head(p)
  # p <- predict(fit, nsim = 5, model = 2, type = "response")
  # head(p)
  # expect_true(all(p > 0))

  p <- predict(fit,
    nsim = 5, model = 2,
    return_tmb_report = TRUE
  )
  expect_length(p, 5L)

  fit <- sdmTMB(
    density ~ depth_scaled,
    data = pcod_2011, mesh = pcod_mesh_2011,
    family = tweedie()
  )
  p <- predict(fit, nsim = 5)
  # p <- predict(fit, nsim = 5, type = "response")
  # head(p)
  # expect_true(all(p > 0))
})

test_that("update() works", {
  skip_on_cran()
  skip_on_ci()
  fit <- sdmTMB(
    density ~ depth_scaled,
    data = pcod_2011, mesh = pcod_mesh_2011,
    family = tweedie()
  )
  fit2 <- update(fit)
  expect_equal(fit$model, fit2$model)
})

test_that("Irregular time gets detected", {
  skip_on_cran()

  mesh <- make_mesh(pcod, c("X", "Y"), n_knots = 50)
  expect_message({
    fit <- sdmTMB(density ~ 1, time = "year",
      data = pcod, mesh = mesh,
      family = tweedie(), spatiotemporal = "ar1",
      do_fit = FALSE
    )
  }, regexp = "irregular time")
  expect_message({
    fit <- sdmTMB(density ~ 1, time = "year",
      data = pcod, mesh = mesh,
      family = tweedie(), spatiotemporal = "ar1",
      do_fit = FALSE
    )
  }, regexp = "c\\(2006, 2008, 2010, 2012, 2014, 2016\\)")

  expect_silent({
    fit <- sdmTMB(density ~ 1, time = "year",
      data = pcod, mesh = mesh,
      family = tweedie(), spatiotemporal = "ar1",
      extra_time = c(2006, 2008, 2010, 2012, 2014, 2016),
      do_fit = FALSE
    )
  })
})

test_that("find_missing_time works", {
  expect_identical(find_missing_time(c(1, 3, 4)), 2)
  expect_identical(find_missing_time(c(1, 3, 4, 5, 7)), c(2, 6))
  expect_identical(find_missing_time(c(1, 2, 3)), numeric(0))
  expect_identical(find_missing_time(c(1, 4, 3)), 2)
  expect_identical(find_missing_time(c(1L, 4L, 3L)), 2L)
})

test_that("offset() throws an error", {
  expect_error({
    fit <- sdmTMB(
      density ~ 1 + offset(depth),
      data = pcod_2011, mesh = pcod_mesh_2011,
      family = tweedie(link = "log")
    )
  }, regexp = "offset")
  expect_error({
    fit <- sdmTMB(
      density ~ 1 + offset(log(depth)),
      data = pcod_2011, mesh = pcod_mesh_2011,
      family = tweedie(link = "log")
    )
  }, regexp = "offset")
})

test_that("sdmTMB throws error on AR1/RW with non-numeric time", {
  skip_on_cran()

  pcod_2011$fyear <- as.factor(pcod_2011$year)
  fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
    time = "fyear", spatiotemporal = "iid", data = pcod_2011,
    mesh = pcod_mesh_2011)
  expect_true(inherits(fit, "sdmTMB"))

  expect_error(
    fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
      time = "fyear", spatiotemporal = "ar1", data = pcod_2011,
      mesh = pcod_mesh_2011), regexp = "Time"
  )

  expect_error(
    fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
      time = "fyear", spatiotemporal = "rw", data = pcod_2011,
      mesh = pcod_mesh_2011), regexp = "Time"
  )

  pcod_2011$cyear <- as.character(pcod_2011$year)
  fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
    time = "cyear", spatiotemporal = "iid", data = pcod_2011,
    mesh = pcod_mesh_2011)
  expect_true(inherits(fit, "sdmTMB"))

  expect_error(
    fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
      time = "cyear", spatiotemporal = "ar1", data = pcod_2011,
      mesh = pcod_mesh_2011), regexp = "Time"
  )

  expect_error(
    fit <- sdmTMB(present ~ 1, family = binomial(), spatial = 'off',
      time = "cyear", spatiotemporal = "rw", data = pcod_2011,
      mesh = pcod_mesh_2011), regexp = "Time"
  )

  fit <- sdmTMB(
    present ~ 0 + depth,
    family = binomial(),
    spatial = 'off',
    time_varying = ~ 1,
    time = "year",
    data = pcod_2011,
    spatiotemporal = "off",
    mesh = pcod_mesh_2011
  )

  expect_error(
    fit <- sdmTMB(
      present ~ 0,
      family = binomial(),
      spatial = 'off',
      time_varying = ~ 1,
      time = "fyear",
      data = pcod_2011,
      spatiotemporal = "off",
      mesh = pcod_mesh_2011
    ))

  expect_error(
    fit <- sdmTMB(
      present ~ 0,
      family = binomial(),
      spatial = 'off',
      time_varying = ~ 1,
      time = "cyear",
      data = pcod_2011,
      spatiotemporal = "off",
      mesh = pcod_mesh_2011
    ))
})