test_that("A logistic threshold 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)
  m <- sdmTMB(density ~ 0 + as.factor(year) + logistic(depth_scaled), data = d,
    mesh = pcod_spde, family = tweedie(link = "log"),
    time = "year")
  expect_true(all(!is.na(summary(m$sd_report)[,"Std. Error"])))

  expect_true("depth_scaled-s50" %in% tidy(m)$term)
  expect_true("depth_scaled-s95" %in% tidy(m)$term)
  expect_true("depth_scaled-smax" %in% tidy(m)$term)
  expect_equal(tidy(m)[,"estimate",drop=TRUE], c(1.555 , 1.655 , 1.718 , 1.138, -0.979, -3.173 , 1.760), tolerance = 1e-3)
})

test_that("A linear threshold 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)
  m <- sdmTMB(density ~ 0 + as.factor(year) + breakpt(depth_scaled), data = d,
    mesh = pcod_spde, family = tweedie(link = "log"), spatial = "off")
  expect_true(all(!is.na(summary(m$sd_report)[,"Std. Error"])))

  expect_true("depth_scaled-slope" %in% tidy(m)$term)
  expect_true("depth_scaled-breakpt" %in% tidy(m)$term)
  expect_equal(tidy(m)[,"estimate",drop=TRUE], c(4.798 , 4.779 , 4.768 , 4.112 , 1.085 ,-1.328), tolerance = 1e-3)
})

test_that("A linear threshold *delta* model fits", {
  skip_on_cran()
  skip_on_ci()

  set.seed(1)
  predictor_dat <- data.frame(
    X = runif(1000), Y = runif(1000),
    a1 = rnorm(1000)
  )
  mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.2)
  s1 <- sdmTMB_simulate(
    formula = ~ 1 + breakpt(a1),
    data = predictor_dat,
    mesh = mesh,
    family = binomial(),
    range = 0.5,
    phi = 0.001,
    sigma_O = 0.1,
    seed = 42,
    B = 0,
    threshold_coefs = c(0.5, 0.3)
  )
  s2 <- sdmTMB_simulate(
    formula = ~ 1 + breakpt(a1),
    data = predictor_dat,
    mesh = mesh,
    family = Gamma(link = "log"),
    range = 0.5,
    phi = 1000,
    sigma_O = 0.1,
    seed = 42,
    B = 0,
    threshold_coefs = c(0.3, 0.3)
  )

  plot(predictor_dat$a1, s1$observed)
  plot(predictor_dat$a1, s2$observed)

  s <- s1
  s$observed <- s1$observed * s2$observed
  s$a1 <- predictor_dat$a1
  s1$a1 <- predictor_dat$a1
  s2$a1 <- predictor_dat$a1

  ctrl <- sdmTMBcontrol(newton_loops = 1L)

  # binomial works:
  fit1 <- sdmTMB(observed ~ breakpt(a1),
    data = s1,
    family = binomial(),
    # mesh = mesh,
    spatial = "off",
    control = ctrl
  )
  print(fit1)

  s2_pos <- subset(s2, s1$observed > 0)
  # mesh2 <- make_mesh(s2_pos, xy_cols = c("X", "Y"), mesh = mesh$mesh)
  # Gamma works:
  fit2 <- sdmTMB(observed ~ breakpt(a1),
    data = s2_pos,
    family = Gamma(link = "log"),
    spatial = "off",
    # mesh = mesh2,
    control = ctrl
  )
  print(fit2)

  fit <- sdmTMB(
    observed ~ breakpt(a1),
    data = s,
    family = delta_gamma(),
    # mesh = mesh,
    spatial = "off",
    control = ctrl
  )
  print(fit)
  sanity(fit)

  t1 <- tidy(fit1)
  t2 <- tidy(fit2)
  td1 <- tidy(fit, model = 1)
  td2 <- tidy(fit, model = 2)

  expect_equal(t1$estimate, td1$estimate, tolerance = 1e-5)
  expect_equal(t2$estimate, td2$estimate, tolerance = 1e-5)
  expect_equal(t1$std.error, td1$std.error, tolerance = 1e-5)
  expect_equal(t2$std.error, td2$std.error, tolerance = 1e-5)
})