library(distfreereg)
set.seed(20250509)
n <- 60
Sig <- rWishart(1, df = n, Sigma = diag(n))[,,1]
theta <- c(2,3)
true_func <- function(theta, x) theta[1]*x[1] + theta[2]*x[2]
test_func <- function(theta, x) theta[1]*x[1] + theta[2]*x[2]^2
true_data <- data.frame(a = rexp(n, rate = 1), b = rnorm(n))
true_X <- as.matrix(true_data)
X <- true_X + rexp(length(true_X))
m <- lm(b ~ a, data = true_data)

cdfr <- compare(true_mean = true_func,
                test_mean = test_func,
                true_X = true_X,
                X = X,
                true_covariance = list(Sigma = Sig),
                covariance = list(Sigma = Sig),
                theta_init = rep(1, length(theta)),
                theta = theta, prog = Inf, keep = 1,
                manual = function(x) sum(x$theta_hat),
                control = list(return_on_error = FALSE))

cdfr

edf <- function(n, reps) matrix(rnorm(n*reps), nrow = n)

cdfr_edf <- update(cdfr, err_dist_fun = "edf")
cdfr_edf

signif(rejection(cdfr, alpha = c(0.1, 0.05))[,2:3], digits = 3)

signif(ks.test(cdfr)[["p.value"]], digits = 3)

plot(cdfr, confband_args = NULL)
plot(cdfr, confband_args = list(N = 1e3))
plot(cdfr, which = "dens", confband_args = NULL)
plot(cdfr, which = "qq")
plot(cdfr, which = "qqp")



# No reps, Y warning, custom tolerances

cdfr_no_reps <- compare(true_mean = true_func,
                        test_mean = test_func,
                        true_X = true_X,
                        Y = rep(1, n),
                        X = X,
                        true_covariance = list(Sigma = Sig),
                        covariance = list(Sigma = Sig),
                        theta_init = rep(1, length(theta)),
                        theta = theta, prog = Inf,
                        manual = function(x) sum(x$theta_hat),
                        control = list(return_on_error = FALSE,
                                       matsqrt_tol = -1e-3,
                                       solve_tol = 1e-3,
                                       sym_tol = 1e-3,
                                       sym_tol1 = 1e-3),
                        reps = 0)



# Global override

smaller_n <- 1:30
cdfr_global_override <- compare(true_mean = true_func,
                                test_mean = test_func,
                                true_X = true_X[smaller_n,],
                                # Y = rep(1, smaller_n),
                                X = X[smaller_n,],
                                true_covariance = list(Sigma = Sig[smaller_n,smaller_n]),
                                covariance = list(Sigma = Sig[smaller_n,smaller_n]),
                                theta_init = rep(1, length(theta)),
                                theta = theta, prog = Inf,
                                manual = function(x) sum(x$theta_hat),
                                control = list(return_on_error = FALSE),
                                global_override = list(theta_hat = theta),
                                reps = 5e2)


plot(x = cdfr, y = cdfr_global_override, confband_args = NULL)
plot(x = cdfr, y = cdfr_global_override, confband_args = NULL, which = "dens")
plot(x = cdfr, y = cdfr_global_override, confband_args = NULL, which = "qq")



# Keep all

cdfr_keep_all <- compare(true_mean = true_func,
                         test_mean = test_func,
                         true_X = true_X,
                         Y = rep(1, n),
                         X = X, keep = "all",
                         true_covariance = list(Sigma = Sig),
                         covariance = list(Sigma = Sig),
                         theta_init = rep(1, length(theta)),
                         theta = theta, prog = Inf,
                         manual = function(x) sum(x$theta_hat),
                         control = list(return_on_error = FALSE),
                         global_override = list(theta_hat = theta),
                         reps = 10)

stopifnot(length(cdfr_keep_all[["dfrs"]]) == 10)


# n supplied

cdfr_n <- compare(true_mean = true_func,
                  test_mean = test_func,
                  n = n,
                  true_X = true_X,
                  X = X, keep = "all",
                  true_covariance = list(Sigma = Sig),
                  covariance = list(Sigma = Sig),
                  theta_init = rep(1, length(theta)),
                  theta = theta, prog = Inf,
                  manual = function(x) sum(x$theta_hat),
                  control = list(return_on_error = FALSE),
                  global_override = list(theta_hat = theta),
                  reps = 3)


# List supplied for Sigma

Sig_array <- rWishart(2, df = n/2, Sigma = diag(n/2))
Sig_list <- list(Sig_array[,,1], Sig_array[,,2])

cdfr_list <- compare(true_mean = true_func,
                     test_mean = test_func,
                     true_X = true_X,
                     X = X,
                     true_covariance = list(Sigma = Sig_list),
                     covariance = list(Sigma = Sig),
                     theta_init = rep(1, length(theta)),
                     theta = theta, prog = Inf, reps = 3,
                     control = list(return_on_error = FALSE))



### Partial outputs

cdfr_partial <- compare(true_mean = true_func,
                        test_mean = test_func,
                        true_X = true_X,
                        Y = rep(1, n),
                        X = X, keep = "all",
                        true_covariance = list(Sigma = Sig),
                        covariance = list(Sigma = Sig),
                        theta_init = rep(1, length(theta)),
                        theta = theta, prog = Inf,
                        control = list(orth_tol = 1e-100),
                        reps = 10)
names(cdfr_partial)




### Failures

tryCatch(plot(cdfr, confband_args = 3), error = function(e) warning(e))
tryCatch(plot(cdfr, density_args = 3), error = function(e) warning(e))
tryCatch(plot(cdfr, hlines = 3), error = function(e) warning(e))
tryCatch(plot(cdfr, curve_args = 3), error = function(e) warning(e))
tryCatch(plot(cdfr, poly = 3), error = function(e) warning(e))
tryCatch(plot(cdfr, legend = 3), error = function(e) warning(e))
tryCatch(plot(cdfr, qqline = 3), error = function(e) warning(e))
tryCatch(plot(cdfr, stat = "hello"), error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE)),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = true_X,
                 X = 1:10,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE)),
         error = function(e) warning(e))

tryCatch(compare(true_mean = function(a) a^2,
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE)),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = NULL,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE)),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = true_X,
                 n = 0,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE)),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = NULL,
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE)),
         error = function(e) warning(e))

tryCatch(compare(true_mean = function(x) x,
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE)),
         error = function(e) warning(e))

err_dist_fun_vec <- function() 1:10
err_dist_fun_mat <- function() matrix(1:10)

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 control = list(return_on_error = FALSE),
                 err_dist_fun = c("err_dist_fun_vec", "err_dist_fun_mat")),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "pi"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 control = list(return_on_error = FALSE),
                 err_dist_args = "pi"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_vec"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 test_mean = test_func,
                 true_X = true_X,
                 true_data = as.data.frame(true_X),
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = b ~ a,
                 true_method = "lm",
                 true_method_args = "4",
                 test_mean = test_func,
                 true_X = true_X,
                 true_data = as.data.frame(true_X),
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = b ~ a,
                 true_method = "lm",
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = b ~ a,
                 true_method = "lm",
                 test_mean = test_func,
                 true_data = as.data.frame(true_X),
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 true_method = "lm",
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 true_method_args = "lm",
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 true_method_args = "lm",
                 test_mean = test_func,
                 true_X = true_X,
                 X = X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 theta = theta, prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = b ~ a,
                 true_method = "lm",
                 test_mean = test_func,
                 true_data = as.data.frame(true_X),
                 data = 1:n,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = b ~ a,
                 true_method = "lm",
                 test_mean = test_func,
                 true_data = as.data.frame(true_X),
                 data = data.frame(),
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = b ~ a,
                 true_method = "lm",
                 test_mean = test_func,
                 true_data = as.data.frame(true_X),
                 data = as.data.frame(true_X),
                 X = true_X,
                 true_covariance = list(Sigma = Sig),
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 prog = Inf, keep = 1,
                 manual = function(x) sum(x$theta_hat),
                 control = list(return_on_error = FALSE),
                 err_dist_fun = "err_dist_fun_mat"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = test_func,
                 true_data = true_data,
                 data = true_X,
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = test_func,
                 theta = 1,
                 data = true_X,
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = m,
                 data = true_X,
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = test_func,
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = test_func,
                 X = true_X,
                 data = true_X,
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = b ~ a,
                 X = true_X,
                 data = true_X,
                 covariance = list(Sigma = Sig),
                 theta_init = rep(1, length(theta)),
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = m,
                 true_covariance = list(Sigma = Sig),
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = m,
                 prog = Inf, keep = "hello"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = m,
                 prog = Inf, manual = "hello"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = m,
                 prog = Inf, err_dist_fun = "hello"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = m,
                 prog = Inf, err_dist_args = "hello"),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = m,
                 prog = Inf, keep = rep(1, 2e4)),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 method = "lm",
                 test_mean = e ~ a,
                 data = true_X,
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(plot(x = cdfr, y = cdfr), error = function(e) warning(e))
tryCatch(plot(x = cdfr, y = "cdfr"), error = function(e) warning(e))
tryCatch(plot(x = cdfr, y = cdfr_global_override, which = "qqp"),
         error = function(e) warning(e))
tryCatch(plot(x = cdfr, y = cdfr_global_override, which = "hello"),
         error = function(e) warning(e))

tryCatch(rejection(cdfr, alpha = 1.1), error = function(e) warning(e))
tryCatch(rejection(cdfr, alpha = -1), error = function(e) warning(e))

tryCatch(plot(x = cdfr, which = "qqqq"), error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 method = "lm",
                 test_mean = b + e ~ a,
                 data = true_X,
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 method = "lm",
                 test_mean = b ~ a,
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = m,
                 test_mean = m,
                 true_X = true_X,
                 prog = Inf, keep = 1),
         error = function(e) warning(e))

tryCatch(compare(true_mean = true_func,
                 true_X = true_X,
                 true_covariance = list(Sigma = Sig),
                 theta = rep(1, length(theta) + 1),
                 test_mean = m,
                 prog = Inf, keep = 1),
         error = function(e) warning(e))






###########################

# These test various combinations of test mean specifications with compare(),
# specifically to test ordering options.

rm(list=ls())

all.equal.compare <- distfreereg:::all.equal.compare

set.seed(20250613)
n <- 20
reps <- 10
data <- data.frame(x1 = rnorm(n), x2 = rnorm(n), y = rnorm(n))

# lm

form_lm <- y ~ x1 + x2
m_lm_1 <- lm(y ~ x1 + x2, data = data)
m_lm_2 <- lm(form_lm, data = data)

set.seed(20250614)
cdfr_lm_1 <- compare(true_mean = m_lm_1, test_mean = m_lm_1, reps = 10, prog = Inf)
set.seed(20250614)
cdfr_lm_2 <- compare(true_mean = m_lm_1, test_mean = m_lm_2, reps = 10, prog = Inf)
signif(cdfr_lm_1[["observed_stats"]], 3)
signif(cdfr_lm_1[["p"]][["KS"]], 3)
stopifnot(all.equal(cdfr_lm_1, cdfr_lm_2))

set.seed(20250614)
cdfr_lm_1_optimal <- update(cdfr_lm_1, ordering = "optimal")
set.seed(20250614)
cdfr_lm_2_optimal <- update(cdfr_lm_2, ordering = "optimal")
signif(cdfr_lm_1_optimal[["observed_stats"]], 3)
signif(cdfr_lm_1_optimal[["p"]][["KS"]], 3)
stopifnot(all.equal(cdfr_lm_1_optimal, cdfr_lm_2_optimal))

set.seed(20250614)
cdfr_lm_1_natural <- update(cdfr_lm_1, ordering = "natural")
set.seed(20250614)
cdfr_lm_2_natural <- update(cdfr_lm_2, ordering = "natural")
signif(cdfr_lm_1_natural[["observed_stats"]], 3)
signif(cdfr_lm_1_natural[["p"]][["KS"]], 3)
stopifnot(all.equal(cdfr_lm_1_natural, cdfr_lm_2_natural))


# glm

m_glm_1 <- glm(y ~ x1 + x2, data = data)
m_glm_2 <- glm(form_lm, data = data)

set.seed(20250614)
cdfr_glm_1 <- compare(true_mean = m_glm_1, test_mean = m_glm_1, reps = 10, prog = Inf)
set.seed(20250614)
cdfr_glm_2 <- compare(true_mean = m_glm_1, test_mean = m_glm_2, reps = 10, prog = Inf)
signif(cdfr_glm_1[["observed_stats"]], 3)
signif(cdfr_glm_1[["p"]][["KS"]], 3)
stopifnot(all.equal(cdfr_glm_1, cdfr_glm_2))

set.seed(20250614)
cdfr_glm_1_optimal <- update(cdfr_glm_1, ordering = "optimal")
set.seed(20250614)
cdfr_glm_2_optimal <- update(cdfr_glm_2, ordering = "optimal")
signif(cdfr_glm_1_optimal[["observed_stats"]], 3)
signif(cdfr_glm_1_optimal[["p"]][["KS"]], 3)
stopifnot(all.equal(cdfr_glm_1_optimal, cdfr_glm_2_optimal))

set.seed(20250614)
cdfr_glm_1_natural <- update(cdfr_glm_1, ordering = "natural")
set.seed(20250614)
cdfr_glm_2_natural <- update(cdfr_glm_2, ordering = "natural")
signif(cdfr_glm_1_natural[["observed_stats"]], 3)
signif(cdfr_glm_1_natural[["p"]][["KS"]], 3)
stopifnot(all.equal(cdfr_glm_1_natural, cdfr_glm_2_natural))



# nls

form_nls <- y ~ a*x1 + b*x2
m_nls_1 <- nls(y ~ a*x1 + b*x2, data = data, start = c(a = 1, b = 1))
m_nls_2 <- nls(form_nls, data = data, start = c(a = 1, b = 1))

set.seed(20250614)
cdfr_nls_1 <- compare(true_mean = m_nls_1, test_mean = m_nls_1,
                      true_covariance = list(Sigma = 1), reps = 10, prog = Inf)
set.seed(20250614)
cdfr_nls_2 <- compare(true_mean = m_nls_1, test_mean = m_nls_2,
                      true_covariance = list(Sigma = 1), reps = 10, prog = Inf)
signif(cdfr_nls_1[["observed_stats"]], 3)
signif(cdfr_nls_1[["p"]][["KS"]], 3)
stopifnot(all.equal(cdfr_nls_1, cdfr_nls_2))

set.seed(20250614)
cdfr_nls_1_optimal <- update(cdfr_nls_1, ordering = "optimal")
set.seed(20250614)
cdfr_nls_2_optimal <- update(cdfr_nls_2, ordering = "optimal")
signif(cdfr_nls_1_optimal[["observed_stats"]], 3)
signif(cdfr_nls_1_optimal[["p"]][["KS"]], 3)
stopifnot(all.equal(cdfr_nls_1_optimal, cdfr_nls_2_optimal))

set.seed(20250614)
cdfr_nls_1_natural <- update(cdfr_nls_1, ordering = "natural")
set.seed(20250614)
cdfr_nls_2_natural <- update(cdfr_nls_2, ordering = "natural")
signif(cdfr_nls_1_natural[["observed_stats"]], 3)
signif(cdfr_nls_1_natural[["p"]][["KS"]], 3)
stopifnot(all.equal(cdfr_nls_1_natural, cdfr_nls_2_natural))





# lmer

if(require(lme4)){
  set.seed(20250614)
  data$g <- factor(sample(1:2, size = n, replace = TRUE))
  
  control_list <- list(action = "ignore", tol = 1)
  control_spec <- lmerControl(check.conv.singular = control_list,
                              check.conv.grad = control_list,
                              check.conv.hess = control_list)
  form_lmer <- y ~ x1 + x2 + (1 | g)
  
  suppressWarnings({
    m_lmer_1 <- lmer(y ~ x1 + x2 + (1 | g), data = data, control = control_spec)
    m_lmer_2 <- lmer(form_lmer, data = data, control = control_spec)
    
    set.seed(20250614)
    cdfr_lmer_1 <- compare(true_mean = m_lmer_1, test_mean = m_lmer_1,
                           reps = 10, prog = Inf)
    set.seed(20250614)
    cdfr_lmer_2 <- compare(true_mean = m_lmer_1, test_mean = m_lmer_2,
                           reps = 10, prog = Inf)
    print(signif(cdfr_lmer_1[["observed_stats"]], 3))
    print(signif(cdfr_lmer_1[["p"]][["KS"]], 3))
    stopifnot(all.equal(cdfr_lmer_1, cdfr_lmer_2))
    
    set.seed(20250614)
    cdfr_lmer_1_optimal <- update(cdfr_lmer_1, ordering = "optimal")
    set.seed(20250614)
    cdfr_lmer_2_optimal <- update(cdfr_lmer_2, ordering = "optimal")
    print(signif(cdfr_lmer_1_optimal[["observed_stats"]], 3))
    print(signif(cdfr_lmer_1_optimal[["p"]][["KS"]], 3))
    stopifnot(all.equal(cdfr_lmer_1_optimal, cdfr_lmer_2_optimal))
    
    set.seed(20250614)
    cdfr_lmer_1_natural <- update(cdfr_lmer_1, ordering = "natural")
    set.seed(20250614)
    cdfr_lmer_2_natural <- update(cdfr_lmer_2, ordering = "natural")
    print(signif(cdfr_lmer_1_natural[["observed_stats"]], 3))
    print(signif(cdfr_lmer_1_natural[["p"]][["KS"]], 3))
    stopifnot(all.equal(cdfr_lmer_1_natural, cdfr_lmer_2_natural))
  })
}