## Tests for GALAHAD 2.0.0
## testthat edition 3

test_that("GALAHAD solves a simple quadratic (all euclidean)", {
  set.seed(1)
  p          <- 5
  theta_true <- c(1, -2, 3, -1, 0.5)
  V    <- function(th) sum((th - theta_true)^2)
  gV   <- function(th) 2 * (th - theta_true)
  parts <- list(positive = integer(0), euclidean = 1:5)
  fit  <- GALAHAD(V, gV, rep(0, p), parts,
                  control = list(max_iter = 500, tol_g = 1e-8))
  expect_true(fit$converged)
  expect_lt(max(abs(fit$theta - theta_true)), 1e-5)
})

test_that("GALAHAD enforces positivity via softplus", {
  ## Objective: negative log-likelihood of exponential distribution
  ## theta = rate > 0
  set.seed(2)
  x    <- rexp(30, rate = 2)
  V    <- function(th) -sum(dexp(x, rate = th[1], log = TRUE))
  gV   <- function(th) galahad_numgrad(V, th)
  parts <- list(positive = 1L, euclidean = integer(0))

  fit <- GALAHAD(V, gV, c(rate = 1), parts,
                 control = list(max_iter = 1000, tol_g = 1e-6))
  expect_true(fit$converged)
  expect_gt(fit$theta[1], 0)
  expect_lt(abs(fit$theta[1] - 2), 0.3)   # within sampling noise
})

test_that("GALAHAD fits exponential decay: A > 0, k > 0", {
  set.seed(3)
  t   <- seq(0, 4, by = 0.5)
  y   <- 3 * exp(-0.8 * t) + rnorm(length(t), sd = 0.1)
  obj <- function(theta) sum((y - theta[1] * exp(-theta[2] * t))^2)
  grd <- function(theta) {
    r  <- y - theta[1] * exp(-theta[2] * t)
    c(-2 * sum(r * exp(-theta[2] * t)),
      -2 * sum(r * (-t) * theta[1] * exp(-theta[2] * t)))
  }
  parts <- list(positive = c(1L, 2L), euclidean = integer(0))
  fit <- GALAHAD(obj, grd, c(2, 0.5), parts,
                 control = list(max_iter = 2000, tol_g = 1e-6))
  expect_true(fit$converged)
  expect_gt(fit$theta[1], 0)
  expect_gt(fit$theta[2], 0)
  expect_lt(abs(fit$theta[1] - 3), 0.5)
  expect_lt(abs(fit$theta[2] - 0.8), 0.3)
})

test_that("legacy T/P/E partition is accepted and gives same result", {
  set.seed(4)
  p          <- 4
  theta_true <- c(1, 2, -1, 0)
  V  <- function(th) sum((th - theta_true)^2)
  gV <- function(th) 2 * (th - theta_true)

  parts_new <- list(positive = c(1L, 2L), euclidean = c(3L, 4L))
  parts_old <- list(T = 1L, P = 2L, E = c(3L, 4L))

  fit_new <- GALAHAD(V, gV, c(0.5, 0.5, 0, 0), parts_new,
                     control = list(max_iter = 500, tol_g = 1e-8))
  fit_old <- GALAHAD(V, gV, c(0.5, 0.5, 0, 0), parts_old,
                     control = list(max_iter = 500, tol_g = 1e-8))

  expect_true(fit_new$converged)
  expect_true(fit_old$converged)
  expect_lt(max(abs(fit_new$theta - fit_old$theta)), 1e-6)
})

test_that("history data.frame has correct columns and length", {
  V  <- function(th) sum(th^2)
  gV <- function(th) 2 * th
  fit <- GALAHAD(V, gV, c(1, 1), list(positive = integer(0), euclidean = 1:2),
                 control = list(max_iter = 50))
  hist <- fit$history
  expect_s3_class(hist, "data.frame")
  expected_cols <- c("iter", "f", "g_inf", "step_norm", "df",
                     "eta", "method", "armijo_iters", "pred_red", "rho")
  expect_true(all(expected_cols %in% names(hist)))
  expect_equal(nrow(hist), fit$iterations)
})

test_that("convergence reason is one of the documented strings", {
  V  <- function(th) sum(th^2)
  gV <- function(th) 2 * th
  fit <- GALAHAD(V, gV, c(1, 1), list(positive = integer(0), euclidean = 1:2))
  expect_true(fit$reason %in%
    c("GRAD_TOL", "FUNC_STALL_ABS", "FUNC_STALL_REL", "MAX_ITER"))
})

test_that("diagnostics includes parameterization = 'softplus'", {
  V  <- function(th) sum(th^2)
  gV <- function(th) 2 * th
  fit <- GALAHAD(V, gV, c(1, 1), list(positive = integer(0), euclidean = 1:2))
  expect_equal(fit$diagnostics$parameterization, "softplus")
})

test_that("L2 regularization shifts optimum toward zero", {
  V  <- function(th) (th[1] - 5)^2 + (th[2] - 5)^2
  gV <- function(th) c(2 * (th[1] - 5), 2 * (th[2] - 5))

  fit_no_reg <- GALAHAD(V, gV, c(0, 0),
                         list(positive = integer(0), euclidean = 1:2),
                         control = list(lambda = 0))
  fit_reg    <- GALAHAD(V, gV, c(0, 0),
                         list(positive = integer(0), euclidean = 1:2),
                         control = list(lambda = 1))

  expect_lt(max(abs(fit_no_reg$theta - 5)), 0.01)     # unregularized -> 5
  expect_lt(max(abs(fit_reg$theta)),         4)        # regularized -> < 5
})

test_that("callback is called and receives correct fields", {
  iters_seen <- integer(0)
  cb <- function(info) {
    iters_seen <<- c(iters_seen, info$iter)
    expect_true(all(c("iter", "theta", "value", "grad_norm") %in% names(info)))
  }
  V  <- function(th) sum(th^2)
  gV <- function(th) 2 * th
  fit <- GALAHAD(V, gV, c(3, 3), list(positive = integer(0), euclidean = 1:2),
                 control = list(max_iter = 20), callback = cb)
  expect_equal(length(iters_seen), fit$iterations)
  expect_equal(iters_seen, seq_len(fit$iterations))
})

test_that("bad parts raises informative error", {
  V  <- function(th) sum(th^2)
  gV <- function(th) 2 * th
  expect_error(
    GALAHAD(V, gV, c(1, 1), list(positive = 1L, euclidean = 1L)),
    regexp = "partition"
  )
  expect_error(
    GALAHAD(V, gV, c(1, 1), list(positive = 1L, euclidean = integer(0))),
    regexp = "partition"
  )
})

test_that("non-positive theta0 for positive partition raises error", {
  V  <- function(th) (th[1] - 2)^2
  gV <- function(th) c(2 * (th[1] - 2))
  expect_error(
    GALAHAD(V, gV, c(-1), list(positive = 1L, euclidean = integer(0))),
    regexp = "> 0"
  )
})

test_that("galahad_numgrad approximates analytical gradient", {
  f   <- function(th) th[1]^2 + 3 * th[2]^3 + th[1] * th[2]
  gf  <- function(th) c(2 * th[1] + th[2], 9 * th[2]^2 + th[1])
  th  <- c(2.1, -1.3)
  ng  <- galahad_numgrad(f, th)
  ag  <- gf(th)
  expect_lt(max(abs(ng - ag)), 1e-5)
})

test_that("galahad_parts validates and returns correct structure", {
  p <- galahad_parts(positive = c(1L, 2L), euclidean = 3L, p = 3)
  expect_equal(p$positive,  c(1L, 2L))
  expect_equal(p$euclidean, 3L)

  expect_error(galahad_parts(positive = 1L, euclidean = 1L),
               regexp = "overlap")
  expect_error(galahad_parts(positive = 1L, euclidean = 3L, p = 3),
               regexp = "cover")
})

test_that("Polyak step used when V_star supplied", {
  V  <- function(th) sum(th^2)
  gV <- function(th) 2 * th
  fit <- GALAHAD(V, gV, c(3, 3),
                 list(positive = integer(0), euclidean = 1:2),
                 control = list(V_star = 0, max_iter = 100))
  expect_true("POLYAK" %in% fit$history$method)
})

test_that("rho is finite after accepted steps", {
  V  <- function(th) sum(th^2)
  gV <- function(th) 2 * th
  fit <- GALAHAD(V, gV, c(5, 5),
                 list(positive = integer(0), euclidean = 1:2),
                 control = list(max_iter = 30))
  rhos <- fit$history$rho
  finite_rhos <- rhos[is.finite(rhos)]
  expect_true(length(finite_rhos) > 0)
})

test_that("max_iter respected when no convergence", {
  ## Rosenbrock from hard start — won't converge in 3 iters with tight tols
  V  <- function(th) 100 * (th[2] - th[1]^2)^2 + (1 - th[1])^2
  gV <- function(th) c(
    -400 * th[1] * (th[2] - th[1]^2) - 2 * (1 - th[1]),
     200 * (th[2] - th[1]^2)
  )
  fit <- GALAHAD(V, gV, c(-1.2, 1),
                 list(positive = integer(0), euclidean = 1:2),
                 control = list(max_iter  = 3,
                                tol_g     = 1e-16,
                                tol_x     = 1e-16,
                                tol_f     = 1e-16,
                                tol_f_rel = 1e-16))
  expect_equal(fit$iterations, 3)
  expect_false(fit$converged)
  expect_equal(fit$reason, "MAX_ITER")
})