diff --git a/AGENTS.md b/AGENTS.md
new file mode 100644
index 0000000..e5d2a73
--- /dev/null
+++ b/AGENTS.md
@@ -0,0 +1,33 @@
+# Highlander package guide
+
+## Package structure
+
+- `R/` contains the package functions. `Highlander()` coordinates CMA-ES and
+  LaplacesDemon; `Lowlander()` performs an initial Latin-hypercube search;
+  `Lowlander_Z.R` contains helpers for reusing Lowlander results; and
+  `Specs_help()` builds algorithm-specific LaplacesDemon settings.
+- `man/` contains the package reference documentation, and `vignettes/` contains
+  longer usage examples.
+- `tests/testthat/` contains the package's unit tests. Testthat edition 3 is
+  configured in `DESCRIPTION`.
+
+## Tests and checks
+
+- Run the unit tests from the repository root with
+  `Rscript -e 'testthat::test_local()'`.
+- Run the full package check with `R CMD build .` followed by
+  `R CMD check --no-manual --as-cran Highlander_*.tar.gz`.
+- Keep unit tests deterministic and fast. Prefer supplied Latin-hypercube
+  matrices and small function inputs over random designs or full sampler runs.
+- The `R-CMD-check` GitHub Actions workflow runs `R CMD check`, which includes
+  the testthat suite. Add separate CI test steps only if they provide coverage
+  that the package check does not.
+
+## Change guidance
+
+- Keep user-facing function behavior and argument defaults aligned with the
+  reference documentation and README examples.
+- When changing exported functions, update their reference documentation and
+  add or adjust focused testthat coverage.
+- Avoid tests that require a long MCMC/CMA optimization run when the same
+  behavior can be checked with a small deterministic unit test.
diff --git a/R/Lowlander.R b/R/Lowlander.R
index 35320a1..dd26592 100644
--- a/R/Lowlander.R
+++ b/R/Lowlander.R
@@ -19,7 +19,8 @@ Lowlander = function(lower, upper, Data, likefunc, Nsamp = 'auto', pcut = 0.1,
     }
   }
 
-  if (length(Nsamp) != 1L || is.na(Nsamp) || !is.finite(Nsamp)) {
+  if (length(Nsamp) != 1L || is.na(Nsamp) || !is.finite(Nsamp) ||
+      Nsamp < 1 || Nsamp != floor(Nsamp) || Nsamp > .Machine$integer.max) {
     stop("Nsamp must be a positive integer (>= 1) or 'auto'.")
   }
   Nsamp = as.integer(Nsamp)
diff --git a/tests/testthat/helper-lowlander-fixture.R b/tests/testthat/helper-lowlander-fixture.R
new file mode 100644
index 0000000..36bf926
--- /dev/null
+++ b/tests/testthat/helper-lowlander-fixture.R
@@ -0,0 +1,11 @@
+lowlander_fixture <- function() {
+  list(
+    latin_mod = matrix(c(0, 10, 1, 11, 2, 12, 3, 13), ncol = 2, byrow = TRUE),
+    keep = 1:4,
+    output = c(4, 1, 3, 2),
+    best = c(first = 1, second = 11),
+    lower = c(0, 10),
+    upper = c(3, 13),
+    liketype = "min"
+  )
+}
diff --git a/tests/testthat/test-likelihood-converters.R b/tests/testthat/test-likelihood-converters.R
new file mode 100644
index 0000000..717c6ec
--- /dev/null
+++ b/tests/testthat/test-likelihood-converters.R
@@ -0,0 +1,35 @@
+test_that("CMA converters apply optimization direction and prior consistently", {
+  scalar_like <- function(parm, Data) parm[1]^2
+  prior <- function(parm, Data) -0.5
+
+  expect_equal(
+    Highlander:::.convert_CMA2CMA(
+      parm = 2, Data = list(), likefunc = scalar_like,
+      liketype = "max", prior = prior
+    ),
+    -3.5
+  )
+
+  converted <- Highlander:::.convert_CMA2LD(
+    parm = 2,
+    Data = list(applyconstraints = FALSE, applyintervals = FALSE),
+    likefunc = scalar_like,
+    liketype = "min",
+    prior = prior
+  )
+  expect_equal(converted$LP, -4.5)
+  expect_equal(converted$Dev, 8)
+  expect_equal(converted$parm, 2)
+})
+
+test_that("LD likelihoods are converted to the CMA objective", {
+  ld_like <- function(parm, Data) list(LP = -4)
+
+  expect_equal(
+    Highlander:::.convert_LD2CMA(
+      parm = 1, Data = list(), likefunc = ld_like,
+      liketype = "max", prior = function(parm, Data) -0.5
+    ),
+    4.5
+  )
+})
diff --git a/tests/testthat/test-lowlander-helpers.R b/tests/testthat/test-lowlander-helpers.R
index 50f4be7..fe10e8a 100644
--- a/tests/testthat/test-lowlander-helpers.R
+++ b/tests/testthat/test-lowlander-helpers.R
@@ -1,14 +1,58 @@
-test_that("Lowlander_kept ranks finite samples in the recorded direction", {
-  result <- list(
-    latin_mod = matrix(c(0, 10, 1, 11, 2, 12, 3, 13), ncol = 2, byrow = TRUE),
-    keep = 1:4,
-    output = c(4, 1, 3, 2),
-    liketype = "min"
-  )
+test_that("Lowlander_kept filters and ranks samples in the recorded direction", {
+  result <- lowlander_fixture()
+  result$keep <- c(1, 2, 3, 4)
   result$output[2] <- Inf
 
-  expect_equal(
-    Lowlander_kept(result),
-    matrix(c(3, 13, 2, 12, 0, 10), ncol = 2, byrow = TRUE)
+  kept <- Lowlander_kept(result, parm.names = c("first", "second"))
+
+  expect_equal(kept, matrix(c(3, 13, 2, 12, 0, 10), ncol = 2, byrow = TRUE,
+                           dimnames = list(NULL, c("first", "second"))))
+
+  result$output <- c(4, 1, 3, 2)
+  result$liketype <- "max"
+  expect_equal(Lowlander_kept(result)[1, ], c(0, 10))
+  expect_null(Lowlander_kept(NULL))
+})
+
+test_that("Lowlander_Z returns spread, bounded starting points", {
+  result <- lowlander_fixture()
+
+  z <- Lowlander_Z(result, N = 2, parm.names = c("first", "second"))
+  expect_equal(dim(z), c(2, 2))
+  expect_equal(colnames(z), c("first", "second"))
+  expect_equal(z, matrix(c(1, 11, 2, 12), ncol = 2, byrow = TRUE,
+                         dimnames = list(NULL, c("first", "second"))))
+
+  set.seed(10)
+  larger_z <- Lowlander_Z(result, N = 6)
+  expect_equal(dim(larger_z), c(6, 2))
+  expect_true(all(larger_z[, 1] >= 0 & larger_z[, 1] <= 3))
+  expect_true(all(larger_z[, 2] >= 10 & larger_z[, 2] <= 13))
+  expect_equal(anyDuplicated(as.data.frame(larger_z)), 0L)
+})
+
+test_that("Lowlander helpers handle unusable inputs and name the best point", {
+  expect_warning(
+    expect_null(Lowlander_kept(list(keep = 1))),
+    "missing latin_mod/keep/output"
+  )
+
+  invalid_keep <- lowlander_fixture()
+  invalid_keep$keep <- 5L
+  expect_warning(
+    expect_null(Lowlander_kept(invalid_keep)),
+    "no valid kept samples"
   )
+
+  nonfinite <- lowlander_fixture()
+  nonfinite$output[nonfinite$keep] <- Inf
+  expect_warning(
+    expect_null(Lowlander_kept(nonfinite)),
+    "All kept Lowlander samples have non-finite likelihood"
+  )
+
+  best <- Lowlander_best(lowlander_fixture(), parm.names = c("first", "second"))
+  expect_equal(best, c(first = 1, second = 11))
+  expect_null(Lowlander_best(NULL))
+  expect_warning(expect_null(Lowlander_best(list())), "no .*best element")
 })
diff --git a/tests/testthat/test-lowlander.R b/tests/testthat/test-lowlander.R
new file mode 100644
index 0000000..9f9e419
--- /dev/null
+++ b/tests/testthat/test-lowlander.R
@@ -0,0 +1,92 @@
+test_that("Lowlander maps design points to bounds and selects the best sample", {
+  latin <- cbind(seq(0, 1, length.out = 6), seq(0, 1, length.out = 6))
+  colnames(latin) <- c("x", "y")
+  result <- Lowlander(
+    lower = c(0, 0),
+    upper = c(1, 1),
+    Data = list(),
+    likefunc = function(parm, Data) sum((parm - c(0.2, 0.2))^2),
+    Nsamp = nrow(latin),
+    pcut = 0.75,
+    latin = latin,
+    parm.names = c("x", "y")
+  )
+
+  expect_equal(result$latin_mod, latin)
+  expect_equal(result$best, c(x = 0.2, y = 0.2))
+  expect_equal(result$keep, 1:4)
+  expect_equal(result$lower, c(x = 0, y = 0))
+  expect_equal(result$upper, c(x = 0.6, y = 0.6))
+  expect_equal(result$frac_keep, 4 / 6)
+  expect_equal(result$bad, integer())
+})
+
+test_that("Lowlander supports maximization and excludes non-finite samples", {
+  latin <- cbind(seq(0, 1, length.out = 6), seq(0, 1, length.out = 6))
+  result <- Lowlander(
+    lower = c(0, 0),
+    upper = c(1, 1),
+    Data = list(),
+    likefunc = function(parm, Data) {
+      if (all(parm == 0.4)) return(NA_real_)
+      sum(parm)
+    },
+    Nsamp = nrow(latin),
+    pcut = 0.8,
+    liketype = "max",
+    latin = latin
+  )
+
+  expect_equal(result$liketype, "max")
+  expect_equal(result$bad, 3L)
+  expect_equal(result$best, c(1, 1))
+  expect_true(all(result$keep %in% setdiff(seq_len(nrow(latin)), result$bad)))
+})
+
+test_that("Lowlander validates its design and scalar likelihood output", {
+  latin <- matrix(seq(0, 1, length.out = 8), ncol = 2)
+  args <- list(
+    lower = c(0, 0),
+    upper = c(1, 1),
+    Data = list(),
+    likefunc = function(parm, Data) sum(parm),
+    Nsamp = nrow(latin),
+    pcut = 0.5,
+    latin = latin
+  )
+
+  expect_error(
+    do.call(Lowlander, utils::modifyList(args, list(upper = 1))),
+    "'lower' and 'upper'"
+  )
+  expect_error(
+    do.call(Lowlander, utils::modifyList(args, list(pcut = 1))),
+    "'pcut'"
+  )
+  expect_error(
+    do.call(Lowlander, utils::modifyList(args, list(Nsamp = 0))),
+    "Nsamp must be a positive integer"
+  )
+  expect_error(
+    do.call(Lowlander, utils::modifyList(args, list(Nsamp = 1.5))),
+    "Nsamp must be a positive integer"
+  )
+  expect_error(
+    do.call(Lowlander, utils::modifyList(
+      args, list(latin = matrix(2, nrow = nrow(latin), ncol = 2))
+    )),
+    "values in 'latin'"
+  )
+  expect_error(
+    do.call(Lowlander, utils::modifyList(
+      args, list(likefunc = function(parm, Data) list(LP = 1, Dev = 2))
+    )),
+    "single numeric value"
+  )
+  expect_error(
+    do.call(Lowlander, utils::modifyList(
+      args, list(likefunc = function(parm, Data) "not numeric")
+    )),
+    "must return numeric values"
+  )
+})
diff --git a/tests/testthat/test-specs-help.R b/tests/testthat/test-specs-help.R
new file mode 100644
index 0000000..6644202
--- /dev/null
+++ b/tests/testthat/test-specs-help.R
@@ -0,0 +1,28 @@
+test_that("Specs_help provides algorithm defaults and seeds ensemble samplers", {
+  data <- list(parm.names = c("first", "second"))
+  lowlander <- lowlander_fixture()
+
+  expect_null(Specs_help("CHARM", data))
+
+  aies <- Specs_help("AIES", data, Lowlander = lowlander, Nc = 5, CPUs = 2)
+  expect_equal(aies$Nc, 6L)
+  expect_equal(dim(aies$Z), c(6, 2))
+  expect_true(all(aies$Z[, 1] >= lowlander$lower[1] &
+                  aies$Z[, 1] <= lowlander$upper[1]))
+
+  demc <- Specs_help("DEMC", data, Lowlander = lowlander,
+                     Iterations = 12, Thinning = 2)
+  expect_equal(dim(demc$Z), c(7, 2, 4))
+  expect_true(all(demc$Z[, 1, ] >= lowlander$lower[1] &
+                  demc$Z[, 1, ] <= lowlander$upper[1]))
+})
+
+test_that("Specs_help scales Griddy-Gibbs grids to plausible bounds", {
+  data <- list(parm.names = c("first", "second"))
+  specs <- Specs_help("AGG", data, Lowlander = lowlander_fixture())
+
+  expect_length(specs$Grid, 2)
+  expect_equal(specs$Grid[[1]][1], -0.75)
+  expect_equal(specs$Grid[[1]][length(specs$Grid[[1]])], 0.75)
+  expect_null(specs$dparm)
+})