context("Testing C++ headers from R") # assign model parameters and number of random effects and fixed effects q <- 4 p <- 5 beta <- sqrt((1:p) / sum(1:p)) Sigma <- diag(q) # # simulate data set # n_clusters <- 4L # set.seed(1) # # sim_dat <- replicate(n_clusters, { # n_members <- sample.int(20L, 1L) + 2L # X <- matrix(runif(p * n_members, -sqrt(6 / 2), sqrt(6 / 2)), # p) # u <- drop(rnorm(q) %*% chol(Sigma)) # Z <- matrix(runif(q * n_members, -sqrt(6 / 2 / q), sqrt(6 / 2 / q)), # q) # eta <- drop(beta %*% X + u %*% Z) # y <- as.numeric((1 + exp(-eta))^(-1) > runif(n_members)) # # list(X = X, Z = Z, y = y, u = u, Sigma_inv = solve(Sigma)) # }, simplify = FALSE) # saveRDS(sim_dat, file.path("tests", "testthat", "mixed-logit.RDS")) f <- file.path("tests", "testthat", "mixed-logit.RDS") if(!file.exists(f)) f <- "mixed-logit.RDS" sim_dat <- readRDS(f) test_that("mixed logit model gives the same", { skip_if_not_installed("RcppArmadillo") skip_if(!has_openmp()) skip_on_macOS() reset_info <- compile_cpp_file("mlogit-ex.cpp") on.exit(reset_compile_cpp_file(reset_info), add = TRUE) setwd(reset_info$old_wd) optimizer <- get_mlogit_optimizer(sim_dat, max_threads = 2L) val <- c(beta, sapply(sim_dat, function(x) x$u)) fv <- eval_mlogit(val = val, ptr = optimizer, n_threads = 1L) expect_equal(fv, 35.248256794375) fv <- eval_mlogit(val = val, ptr = optimizer, n_threads = 2L) expect_equal(fv, 35.248256794375) gr_res <- structure( c(-3.09880891801443, -3.11919456739973, -1.26969269606313, 5.50412085634227, -5.07396341960751, 0.407174373086653, 0.590879440022872, 0.122784612825495, 0.222039129397368, -0.335740954704027, 0.967036607234169, -2.05362226641706, -0.517311694420387, -2.33859175587387, 0.911175257815629, -1.80538928384093, -0.81631781729293, 3.14749812772018, -0.894143195465396, 0.647354085121652, 0.344172484580979), value = 35.248256794375) gr <- grad_mlogit(val = val, ptr = optimizer, n_threads = 1L) expect_equal(gr, gr_res) gr <- grad_mlogit(val = val, ptr = optimizer, n_threads = 2L) expect_equal(gr, gr_res) # hess_truth <- Matrix( # numDeriv::jacobian(grad_mlogit, val, ptr = optimizer, n_threads = 1L), # sparse = TRUE) # saveRDS(hess_truth, "test-cpp-api-true-hess.RDS") hess <- true_hess_sparse(optimizer, val) expect_true(isTRUE(all.equal(hess, readRDS("test-cpp-api-true-hess.RDS")))) rel_eps <- sqrt(.Machine$double.eps) opt <- optim_mlogit( val = val, ptr = optimizer, rel_eps = rel_eps, max_it = 100L, c1 = 1e-4, c2 = .9, n_threads = 1L) opt_res <- list(par = c(0.647310437469711, 0.773205671911178, 0.656983834118283, 0.119681528415682, 0.855849653168437, -0.484311298187747, 0.412484752266496, 0.766994317977443, 0.422874453513229, 0.605169744625452, -0.783701817869368, 0.096194899967712, -0.922934757392211, -0.483542761956257, 0.438193243670319, -0.493309514536335, 0.0274193762771719, 0.426816737611666, 0.0281704108821127, -0.180877909582702, -0.150512019701299), value = 25.1215015278071, info = 0L, counts = c(`function` = 12, gradient = 11, n_cg = 40 ), convergence = TRUE) tol <- sqrt(rel_eps) * 10 do_check <- !names(opt) %in% "counts" expect_equal(opt[do_check], opt_res[do_check], tolerance = tol) opt <- optim_mlogit( val = val, ptr = optimizer, rel_eps = rel_eps, max_it = 100L, c1 = 1e-4, c2 = .9, n_threads = 2L) expect_equal(opt[do_check], opt_res[do_check], tolerance = tol) # the Hessian yields the same result hess_ress_sparse <- get_sparse_Hess_approx_mlogit(optimizer) expect_known_value(hess_ress_sparse, "mlogit-ex-hess-res.RDS") hess_ress <- get_Hess_approx_mlogit(optimizer) expect_equal(as.matrix(hess_ress_sparse), hess_ress, check.attributes = FALSE) # works with other preconditioners for(i in 0:3){ opt_new <- optim_mlogit( val = val, ptr = optimizer, rel_eps = rel_eps, max_it = 100L, c1 = 1e-4, c2 = .9, n_threads = 2L, pre_method = i) expect_equal(opt$value, opt_new$value, info = i) expect_equal(opt$par, opt_new$par, info = i, tolerance = 4 * sqrt(rel_eps)) } # the gradient tolerance works start_priv <- opt$par start_priv[-seq_along(beta)] <- 0 gr_tol <- 1e-7 start <- optim_mlogit_private( start_priv, optimizer, rel_eps = 1, n_threads = 1, gr_tol = gr_tol, max_it = 1000, c1 = 1e-4, c2 = .9) gr_start <- grad_mlogit(val = start, ptr = optimizer, n_threads = 2L) gr_start <- gr_start[-seq_along(beta)] n_clusters <- length(sim_dat) gr_norm <- tapply(gr_start, gl(n_clusters, q), function(x) sqrt(sum(x^2))) expect_true(all(gr_norm < gr_tol)) opt <- optim_mlogit( val = val, ptr = optimizer, rel_eps = 1, max_it = 1000L, c1 = 1e-4, c2 = .9, n_threads = 2L, gr_tol = gr_tol) gr_opt <- grad_mlogit(val = opt$par, ptr = optimizer, n_threads = 2L) expect_lt(sqrt(sum(gr_opt^2)), gr_tol) # works with masking idx_mask <- c(7L, 0L, 20L, 15L) par_fix <- c(0.2582, 0.36515, 0.44721, 0.5164, 0.57735, -0.01619, 0.94384, 0.82122, 0.5939, 0.14377, -0.11775, -0.91207, -1.43759, -1.91436, 1.17658, -1.66497, -0.46353, 1.40856, -0.54176, 0.27866, -0.19397) set_masked(optimizer, idx_mask) opt_mask <- optim_mlogit( val = par_fix, ptr = optimizer, rel_eps = rel_eps, max_it = 100L, c1 = 1e-4, c2 = .9, n_threads = 2L) clear_masked(optimizer) opt_res <- list(par = c(0.2582, 0.830931465957626, 0.586329698291959, 0.136054541268922, 0.758404096966284, -0.508402916414152, 0.397840484268012, 0.82122, 0.385471793612842, 0.601355805758345, -0.772664836507305, 0.227470250425005, -0.990712066557003, -0.40427168450336, 0.473345494155202, -1.66497, 0.0985617406985485, 0.40708903513448, 0.105020277870488, -0.0635572013637352, -0.19397), value = 26.6443138686727, info = 0L, counts = c( `function` = 12, gradient = 10, n_cg = 16), convergence = TRUE) expect_equal(opt_mask[do_check], opt_res[do_check], tolerance = tol) expect_equal(opt_mask$par[idx_mask + 1L], par_fix[idx_mask + 1L]) # check the function to optimize the private parameters start_priv <- opt$par start_priv[-seq_along(beta)] <- 0 out <- optim_mlogit_private( val = start_priv, ptr = optimizer, rel_eps = rel_eps, max_it = 100L, c1 = 1e-4, c2 = .9, n_threads = 2L) expect_equal(opt$par, out, tolerance = tol, check.attributes = FALSE) expect_true(all(opt$par[seq_along(beta)] == out[seq_along(beta)])) expect_equal(opt$value, attr(out, "value"), tolerance = rel_eps) opt <- optim_mlogit( val = val, ptr = optimizer, rel_eps = rel_eps, max_it = 100L, c1 = 1e-4, c2 = .9, use_bfgs = FALSE, n_threads = 1L) opt_res <- list(par = c(0.647309984423475, 0.773099829472392, 0.657030091436329, 0.119542556826119, 0.856048420971233, -0.48432079434579, 0.412471124068345, 0.76694534092756, 0.422786612641562, 0.605124122567012, -0.783742636851702, 0.0960937729102007, -0.922965266682886, -0.483501533118553, 0.438327834070562, -0.493318173514831, 0.0274749226734802, 0.426944683060686, 0.0287992526677307, -0.180834889398274, -0.150315014163914), value = 25.121501682662, info = 0L, counts = c(`function` = 14, gradient = 11, n_cg = 22 ), convergence = TRUE) expect_equal(opt[do_check], opt_res[do_check], tolerance = tol) opt <- optim_mlogit( val = val, ptr = optimizer, rel_eps = rel_eps, max_it = 100L, c1 = 1e-4, c2 = .9, use_bfgs = FALSE, n_threads = 2L) expect_equal(opt[do_check], opt_res[do_check], tolerance = tol) opt <- optim_mlogit( val = val, ptr = optimizer, rel_eps = rel_eps, max_it = 100L, c1 = 1e-4, c2 = .9, n_threads = 1L, strong_wolfe = FALSE) opt_res <- list(par = c(0.647310437469711, 0.773205671911178, 0.656983834118283, 0.119681528415682, 0.855849653168437, -0.484311298187747, 0.412484752266496, 0.766994317977443, 0.422874453513229, 0.605169744625452, -0.783701817869368, 0.096194899967712, -0.922934757392211, -0.483542761956257, 0.438193243670319, -0.493309514536335, 0.0274193762771719, 0.426816737611666, 0.0281704108821127, -0.180877909582702, -0.150512019701299), value = 25.1215015278071, info = 0L, counts = c(`function` = 12, gradient = 11, n_cg = 40 ), convergence = TRUE) expect_equal(opt[do_check], opt_res[do_check], tolerance = tol) })