# examples already tested # Setup the S matrix for every test perm_size <- 6 mu <- numeric(perm_size) # sigma_matrix is a matrix invariant under permutation (1,2,3,4,5,6) sigma_matrix <- matrix( data = c( 1.1, 0.8, 0.6, 0.4, 0.6, 0.8, 0.8, 1.1, 0.8, 0.6, 0.4, 0.6, 0.6, 0.8, 1.1, 0.8, 0.6, 0.4, 0.4, 0.6, 0.8, 1.1, 0.8, 0.6, 0.6, 0.4, 0.6, 0.8, 1.1, 0.8, 0.8, 0.6, 0.4, 0.6, 0.8, 1.1 ), nrow = perm_size, byrow = TRUE ) number_of_observations <- 13 Z <- MASS::mvrnorm(number_of_observations, mu = mu, Sigma = sigma_matrix) # this place is random S <- (t(Z) %*% Z) / number_of_observations # the theoretical mean is 0 test_that("Setting custom permutation in gips constructor works", { custom_perm1 <- gips_perm("(1,2)(3,4,5,6)", 6) g1 <- gips( S, number_of_observations, was_mean_estimated = FALSE, perm = custom_perm1 ) custom_perm2 <- permutations::permutation("(1,2)(3,4,5)(6)") g2 <- gips( S, number_of_observations, was_mean_estimated = FALSE, perm = custom_perm2 ) expect_identical(custom_perm1, g1[[1]]) expect_identical(gips_perm(custom_perm2, ncol(S)), g2[[1]]) # gips as the `perm` parameter g3 <- gips( S, number_of_observations, was_mean_estimated = FALSE, perm = g2 ) expect_equal(g2, g3) }) test_that("new_gips() works or throws an erron on wrong arguments", { expect_silent(new_gips( list(gips_perm("(1,2)(3,4,5,6)", 6)), S, number_of_observations, 3, diag(nrow = ncol(S)), FALSE, NULL )) expect_error(new_gips( gips_perm("(1,2)(3,4,5,6)", 6), S, number_of_observations, 3, diag(nrow = ncol(S)), FALSE, NULL )) expect_error(new_gips( list(gips_perm("(1,2)(3,4,5,6)", 6)), 7, number_of_observations, 3, diag(nrow = ncol(S)), FALSE, NULL )) expect_error(new_gips( list(gips_perm("(1,2)(3,4,5,6)", 6)), S, number_of_observations + 0.1, 3, diag(nrow = ncol(S)), FALSE, NULL )) expect_error(new_gips( list(gips_perm("(1,2)(3,4,5,6)", 6)), S, "r", 3, diag(nrow = ncol(S)), FALSE, NULL )) expect_error(new_gips( list(gips_perm("(1,2)(3,4,5,6)", 6)), S, number_of_observations, "r", diag(nrow = ncol(S)), NULL )) expect_error(new_gips( list(gips_perm("(1,2)(3,4,5,6)", 6)), S, number_of_observations, 3, 7, FALSE, NULL )) expect_error(new_gips( list(gips_perm("(1,2)(3,4,5,6)", 6)), S, number_of_observations, 3, "diag(nrow = ncol(S))", FALSE, NULL )) expect_error(new_gips( list(gips_perm("(1,2)(3,4,5,6)", 6)), S, number_of_observations, 3, diag(nrow = ncol(S)), "FALSE", NULL )) expect_error(new_gips( list(gips_perm("(1,2)(3,4,5,6)", 6)), S, number_of_observations, 3, diag(nrow = ncol(S)), FALSE, "NULL" )) }) test_that("Properly validate the gips class with no optimization or after a single optimization", { custom_perm1 <- gips_perm("(1,2)(3,4,5,6)", 6) g1 <- gips( S, number_of_observations, was_mean_estimated = FALSE, perm = custom_perm1 ) g2 <- find_MAP(g1, max_iter = 3, show_progress_bar = FALSE, optimizer = "MH", return_probabilities = TRUE, save_all_perms = TRUE ) while (attr(g2, "optimization_info")[["acceptance_rate"]] == 0) { # Around 4% of time, in the optimization all permutations were rejected. Is such a case, try again. We want g2 to have at least 1 accepted transposition. g2 <- find_MAP(g1, max_iter = 3, show_progress_bar = FALSE, optimizer = "MH", return_probabilities = TRUE, save_all_perms = TRUE ) } g3 <- find_MAP(g1, max_iter = 3, show_progress_bar = FALSE, optimizer = "HC", return_probabilities = FALSE ) g_BF <- find_MAP( gips( matrix(c( 1, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 1 ), nrow = 3), number_of_observations, was_mean_estimated = FALSE ), optimizer = "BF", show_progress_bar = FALSE ) g_BF_prob <- find_MAP( gips( matrix(c( 1, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 1 ), nrow = 3), number_of_observations, was_mean_estimated = FALSE ), optimizer = "BF", show_progress_bar = FALSE, return_probabilities = TRUE, save_all_perms = TRUE ) # tests: expect_silent(validate_gips(g1)) expect_silent(validate_gips(g2)) expect_silent(validate_gips(g3)) expect_silent(validate_gips(g_BF)) expect_silent(validate_gips(g_BF_prob)) # NaNs or Infs in D_matrix expect_error(gips( S[1:4, 1:4], number_of_observations, was_mean_estimated = FALSE, D_matrix = diag(4) * 1e350 )) expect_error(gips( S[1:4, 1:4], number_of_observations, was_mean_estimated = FALSE, D_matrix = diag(Inf, 4) )) # Other tests g_err <- g2 class(g_err[[1]]) <- "test" expect_error( validate_gips(g_err), "must be an object of a `gips_perm` class." ) g_err <- g2 g_err[[2]] <- g_err[[1]] expect_error( validate_gips(g_err), "The `length\\(g\\)` must be `1`." ) g_err <- "(1,2)(3,4,5)" attributes(g_err) <- attributes(g2) expect_error( validate_gips(g_err), "The `g` must be a list." ) g_err <- list("(1,2)(3,4,5)") class(g_err[[1]]) <- "gips_perm" attributes(g_err) <- attributes(g2) expect_error( validate_gips(g_err), "You provided `g\\[\\[1\\]\\]` with `class\\(g\\[\\[1\\]\\]\\) == 'gips_perm'`, but your g\\[\\[1\\]\\] does not pass `validate_gips_perm\\(g\\[\\[1\\]\\]\\)`." ) # test of "optimization_info" validation: g_err <- g2 attr(g_err, "optimization_info") <- "test" expect_error( validate_gips(g_err), "The `optimization_info` value must be either a `NULL`, or a list." ) g_err <- g2 attr(g_err, "optimization_info")[["non_existing"]] <- "test" expect_error( validate_gips(g_err), "You have a list of 16 elements." ) g_err <- g2 attr(g_err, "optimization_info")[["acceptance_rate"]] <- NULL expect_error( validate_gips(g_err), "You have a list of 14 elements." ) g_err <- g2 attr(g_err, "optimization_info")[["non_existing"]] <- "test" attr(g_err, "optimization_info")[["acceptance_rate"]] <- NULL expect_error( validate_gips(g_err), "You have a list of 15 elements." ) # this one showed an error that one have the list of proper number of elements, which is actually expected, but the names of the fields are not expected. g_err <- g2 attr(g_err, "optimization_info")[["acceptance_rate"]] <- -0.1 expect_error( validate_gips(g_err), "must be a number in range \\[0, 1\\]." ) g_err <- g2 attr(g_err, "optimization_info")[["optimization_algorithm_used"]] <- "brute_force" # in general, brute_force is all right, but then the acceptance_rate has to be NULL expect_error( validate_gips(g_err), "When brute force algorithm was used for optimization, `attr\\(g, 'optimization_info'\\)\\[\\['acceptance_rate'\\]\\]` must be a `NULL`" ) g_err <- g2 attr(g_err, "optimization_info")[["log_posteriori_values"]] <- as.character(attr(g_err, "optimization_info")[["log_posteriori_values"]]) expect_error( validate_gips(g_err), "must be a vector of numbers." ) g_err <- g2 attr(g_err, "optimization_info")[["visited_perms"]] <- "text" expect_error( validate_gips(g_err), "must be a list or `NA`." ) g_err <- g2 class(attr(g_err, "optimization_info")[["visited_perms"]][[1]]) <- "test" expect_error( validate_gips(g_err), "must be of a `gips_perm` class." ) g_err <- g2 attr(g_err, "optimization_info")[["visited_perms"]] <- list() expect_error( validate_gips(g_err), "is a list, but of a length 0." ) g_err <- g_BF attr(g_err, "optimization_info")[["visited_perms"]] <- list() expect_error( validate_gips(g_err), "is a list, but of a length 0" ) g_err <- g_BF attr(g_err, "optimization_info")[["visited_perms"]] <- "()" expect_error( validate_gips(g_err), "You have `attr\\(g, 'optimization_info'\\)\\[\\['visited_perms'\\]\\]` of type character." ) g_err <- g_BF attr(g_err, "optimization_info")[["last_perm"]] <- "()" expect_error( validate_gips(g_err), "You have `attr\\(g, 'optimization_info'\\)\\[\\['last_perm'\\]\\]` of type character." ) g_err <- g2 attr(g_err, "optimization_info")[["last_perm"]] <- 7 expect_error( validate_gips(g_err), "must be the last element of" ) g_err <- g2 attr(g_err, "optimization_info")[["last_perm"]] <- gips_perm("", 6) expect_error( validate_gips(g_err), "must be the log_posteriori of" ) g_err <- g2 fake_gips_perm <- "()" class(fake_gips_perm) <- "gips_perm" attr(g_err, "optimization_info")[["last_perm"]] <- fake_gips_perm expect_error( validate_gips(g_err), ", but your attr\\(g, 'optimization_info'\\)\\[\\['last_perm'\\]\\] does not pass `validate_gips_perm" ) g_err <- g2 attr(g_err, "optimization_info")[["last_perm_log_posteriori"]] <- 7 expect_error( validate_gips(g_err), "must be the log_posteriori of" ) g_err <- g2 attr(g_err, "optimization_info")[["last_perm_log_posteriori"]] <- "7" expect_error( validate_gips(g_err), "must be an object of a `gips_perm` class." ) g_err <- g2 attr(g_err, "optimization_info")[["iterations_performed"]] <- 4 # there were 2 iterations + init (+ additional 2 iterations if repeated); so it cannot be 4 expect_error( validate_gips(g_err), "In every iteration at least one value of log_posteriori is calculated." ) g_err <- g2 attr(g_err, "optimization_info")[["iterations_performed"]] <- 4.5 expect_error( validate_gips(g_err), "attr\\(g, 'optimization_info'\\)\\[\\['iterations_performed'\\]\\]` must be a vector of whole numbers" ) g_err <- g2 attr(g_err, "optimization_info")[["optimization_algorithm_used"]] <- "MH" # Even if MH was used, it would produce the text "Metropolis_Hastings" expect_error( validate_gips(g_err), "The available optimization algorithms are 'Metropolis_Hastings', 'hill_climbing' and 'brute_force'." ) g_err <- g2 attr(g_err, "optimization_info")[["optimization_algorithm_used"]] <- "hill_climbing" # hill_climbing is legal, but the post_probabilities are not with this optimization_algorithm_used expect_error( validate_gips(g_err), "`post_probabilities` can only be obtained with 'Metropolis_Hastings' or 'brute_force' optimization method." ) g_err <- g2 attr(g_err, "optimization_info")[["post_probabilities"]] <- attr(g_err, "optimization_info")[["post_probabilities"]] + 0.1 expect_error( validate_gips(g_err), "must have properties of probability. All elements in range" ) g_err <- g2 len_post_prob <- length(attr(g_err, "optimization_info")[["post_probabilities"]]) attr(g_err, "optimization_info")[["post_probabilities"]] <- attr(g_err, "optimization_info")[["post_probabilities"]] + c(0, rep(1 / (len_post_prob - 1), len_post_prob - 1)) expect_error( validate_gips(g_err), "must have properties of probability. All elements in range" ) g_good <- g2 attr(g_good, "optimization_info")[["post_probabilities"]] <- c(1, rep(0, length(attr(g_good, "optimization_info")[["post_probabilities"]]) - 1)) # this could underflow to 0 names(attr(g_good, "optimization_info")[["post_probabilities"]]) <- names(attr(g2, "optimization_info")[["post_probabilities"]]) expect_silent( validate_gips(g_good) ) g_err <- g2 attr(g_err, "optimization_info")[["did_converge"]] <- TRUE expect_error( validate_gips(g_err), "`did_converge` can only be obtained with 'hill_climbing' or 'brute_force' optimization method." ) g_err <- g3 attr(g_err, "optimization_info")[["did_converge"]] <- "test" expect_error( validate_gips(g_err), "When 'hill_climbing' optimization method, the `did_converge` must be `TRUE` or `FALSE`." ) g_err <- g3 attr(g_err, "optimization_info")[["did_converge"]] <- NA expect_error( validate_gips(g_err), "When 'hill_climbing' optimization method, the `did_converge` must be `TRUE` or `FALSE`." ) g_err <- g2 attr(g_err, "optimization_info")[["best_perm_log_posteriori"]] <- 7 expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["optimization_time"]] <- I(NA) expect_error( validate_gips(g_err), "is initially set to `NA`, but that state of the gips object should not be available to the user." ) g_err <- g2 time_now <- Sys.time() attr(g_err, "optimization_info")[["optimization_time"]] <- time_now expect_error( validate_gips(g_err), "has to be of a class 'difftime'." ) g_err <- g2 time_now <- Sys.time() attr(g_err, "optimization_info")[["optimization_time"]] <- 1 expect_error( validate_gips(g_err), " has to be of a class 'difftime'." ) g_err <- g2 time_now2 <- Sys.time() attr(g_err, "optimization_info")[["optimization_time"]] <- (time_now - time_now2) expect_error( validate_gips(g_err), "has to be a non negative time difference." ) g_err <- g2 attr(g_err, "optimization_info")[["whole_optimization_time"]] <- NA expect_error( validate_gips(g_err), "You have `is.na\\(attr\\(g, 'optimization_info'\\)\\[\\['whole_optimization_time'\\]\\]\\)" ) g_err <- g2 attr(g_err, "optimization_info")[["whole_optimization_time"]] <- 7 expect_error( validate_gips(g_err), "You have `attr\\(g, 'optimization_info'\\)\\[\\['whole_optimization_time'\\]\\]` of a class \\(numeric\\)" ) g_err <- g2 time_now1 <- Sys.time() attr(g_err, "optimization_info")[["whole_optimization_time"]] <- time_now2 - time_now1 expect_error( validate_gips(g_err), "`attr\\(g, 'optimization_info'\\)\\[\\['whole_optimization_time'\\]\\]` has to be a non negative time difference" ) # more than 5 problems at the time: g_err <- g2 attr(g_err, "optimization_info")[["acceptance_rate"]] <- -0.1 attr(g_err, "optimization_info")[["log_posteriori_values"]] <- as.character(attr(g_err, "optimization_info")[["log_posteriori_values"]]) attr(g_err, "optimization_info")[["last_perm_log_posteriori"]] <- 7 attr(g_err, "optimization_info")[["iterations_performed"]] <- 40 attr(g_err, "optimization_info")[["post_probabilities"]] <- attr(g_err, "optimization_info")[["post_probabilities"]] + c(0, rep(1 / (len_post_prob - 1), len_post_prob - 1)) attr(g_err, "optimization_info")[["best_perm_log_posteriori"]] <- 7 attr(g_err, "optimization_info")[["optimization_time"]] <- I(NA) expect_error( validate_gips(g_err), "and 3 more problems" ) }) test_that("Properly validate the gips class after multiple optimizations", { custom_perm1 <- gips_perm("(1,2)(3,4,5,6)", 6) g1 <- gips(S, number_of_observations, was_mean_estimated = FALSE, perm = custom_perm1 ) g2 <- find_MAP(g1, max_iter = 3, show_progress_bar = FALSE, optimizer = "MH", return_probabilities = TRUE, save_all_perms = TRUE) g2 <- find_MAP(g2, max_iter = 3, show_progress_bar = FALSE, optimizer = "MH", return_probabilities = TRUE, save_all_perms = TRUE) g2 <- find_MAP(g2, max_iter = 3, show_progress_bar = FALSE, optimizer = "MH", return_probabilities = TRUE, save_all_perms = TRUE) while (attr(g2, "optimization_info")[["acceptance_rate"]] == 0) { # Around 4% of time, in the optimization all permutations were rejected. Is such a case, try again. g2 <- find_MAP(g2, max_iter = 3, show_progress_bar = FALSE, optimizer = "MH", return_probabilities = TRUE, save_all_perms = TRUE) g2 <- find_MAP(g2, max_iter = 3, show_progress_bar = FALSE, optimizer = "MH", return_probabilities = TRUE, save_all_perms = TRUE) g2 <- find_MAP(g2, max_iter = 3, show_progress_bar = FALSE, optimizer = "MH", return_probabilities = TRUE, save_all_perms = TRUE) } g3 <- find_MAP(g1, max_iter = 3, show_progress_bar = FALSE, optimizer = "HC", return_probabilities = FALSE) g3 <- find_MAP(g3, max_iter = 3, show_progress_bar = FALSE, optimizer = "HC", return_probabilities = FALSE) g3 <- find_MAP(g3, max_iter = 3, show_progress_bar = FALSE, optimizer = "HC", return_probabilities = FALSE) expect_warning(expect_message(expect_message( g_MH_MH <- find_MAP( find_MAP( gips( matrix(c(1, 0.5, 0.5, 5), nrow = 2), number_of_observations, was_mean_estimated = FALSE ), optimizer = "MH", show_progress_bar = FALSE, return_probabilities = FALSE, max_iter = 3 ), optimizer = "MH", show_progress_bar = FALSE, return_probabilities = TRUE, max_iter = 3, save_all_perms = TRUE ) ))) # 2 messages and a warning # tests: expect_silent(validate_gips(g1)) expect_silent(validate_gips(g2)) expect_silent(validate_gips(g3)) expect_silent(validate_gips(g_MH_MH)) expect_true(length(attr(g2, "opt")[["visited_perms"]]) >= 9) expect_false(is.null(attr(g2, "optimization_info")[["post_probabilities"]])) expect_true(is.null(attr(g_MH_MH, "optimization_info")[["post_probabilities"]])) # first one was optimized without saving, so for the second one it was forgotten g_err <- g2 class(g_err[[1]]) <- "test" expect_error(validate_gips(g_err)) g_err <- g2 g_err[[2]] <- g_err[[1]] expect_error(validate_gips(g_err)) g_err <- "(1,2)(3,4,5)" attributes(g_err) <- attributes(g2) expect_error(validate_gips(g_err)) # test of "optimization_info" validation: g_err <- g2 attr(g_err, "optimization_info") <- "test" expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["non_existing"]] <- "test" expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["acceptance_rate"]] <- NULL expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["non_existing"]] <- "test" attr(g_err, "optimization_info")[["acceptance_rate"]] <- NULL expect_error(validate_gips(g_err)) # this one shows an error that one have the list of 10 elements, which is actually expected, but the names of the fields are not expected. g_err <- g2 attr(g_err, "optimization_info")[["acceptance_rate"]] <- -0.1 expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["log_posteriori_values"]] <- as.character(attr(g_err, "optimization_info")[["log_posteriori_values"]]) expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["visited_perms"]] <- "text" expect_error(validate_gips(g_err)) g_err <- g2 class(attr(g_err, "optimization_info")[["visited_perms"]][[1]]) <- "test" expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["last_perm"]] <- 7 expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["last_perm"]] <- gips_perm("", 6) expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["last_perm_log_posteriori"]] <- 7 expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["last_perm_log_posteriori"]] <- "7" expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["iterations_performed"]] <- 40 expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["optimization_algorithm_used"]] <- "MH" # Even if MH was used, it would produce the text "Metropolis_Hastings" expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["optimization_algorithm_used"]] <- "hill_climbing" # hill_climbing is legal, but the post_probabilities are not with this optimization_algorithm_used expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["post_probabilities"]] <- attr(g_err, "optimization_info")[["post_probabilities"]] + 0.1 expect_error(validate_gips(g_err)) g_err <- g2 len_post_prob <- length(attr(g_err, "optimization_info")[["post_probabilities"]]) attr(g_err, "optimization_info")[["post_probabilities"]] <- attr(g_err, "optimization_info")[["post_probabilities"]] + c(0, rep(1 / (len_post_prob - 1), len_post_prob - 1)) expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["did_converge"]] <- TRUE expect_error(validate_gips(g_err)) g_err <- g3 attr(g_err, "optimization_info")[["did_converge"]] <- "test" expect_error(validate_gips(g_err)) g_err <- g3 attr(g_err, "optimization_info")[["did_converge"]] <- NA expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["best_perm_log_posteriori"]] <- 7 expect_error(validate_gips(g_err)) g_err <- g2 attr(g_err, "optimization_info")[["optimization_time"]] <- I(NA) expect_error(validate_gips(g_err)) g_err <- g2 time_now <- Sys.time() attr(g_err, "optimization_info")[["optimization_time"]] <- time_now expect_error(validate_gips(g_err)) g_err <- g2 time_now <- Sys.time() attr(g_err, "optimization_info")[["optimization_time"]] <- 1 expect_error(validate_gips(g_err)) g_err <- g2 time_now2 <- Sys.time() attr(g_err, "optimization_info")[["optimization_time"]] <- (time_now - time_now2) expect_error(validate_gips(g_err)) # more than 5 problems at the time: g_err <- g2 attr(g_err, "optimization_info")[["acceptance_rate"]] <- -0.1 attr(g_err, "optimization_info")[["log_posteriori_values"]] <- as.character(attr(g_err, "optimization_info")[["log_posteriori_values"]]) attr(g_err, "optimization_info")[["last_perm_log_posteriori"]] <- 7 attr(g_err, "optimization_info")[["iterations_performed"]] <- c(2, 2) attr(g_err, "optimization_info")[["post_probabilities"]] <- attr(g_err, "optimization_info")[["post_probabilities"]] + c(0, rep(1 / (len_post_prob - 1), len_post_prob - 1)) attr(g_err, "optimization_info")[["best_perm_log_posteriori"]] <- 7 attr(g_err, "optimization_info")[["optimization_time"]] <- I(NA) expect_error(validate_gips(g_err)) }) test_that("Process proper parameters", { expect_silent(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations = number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_silent(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations = number_of_observations, max_iter = 10, start_perm = gips_perm(example_perm, 6), delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) }) test_that("check_correctness_of_arguments() properly validates arguments", { expect_silent(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( 6, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( matrix(1:30, ncol = 5), number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( matrix(c(LETTERS, LETTERS)[1:36], ncol = 6), number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) S_nonsymetric <- S S_nonsymetric[1, 2] <- S[1, 2] - 1 expect_error(check_correctness_of_arguments( S_nonsymetric, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) S_non_positive_semi_definite <- S - diag(ncol(S)) * eigen(S, symmetric = TRUE, only.values = TRUE)[["values"]][2] expect_error(check_correctness_of_arguments( S_non_positive_semi_definite, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, NULL, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, 0, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations + 0.1, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30.1, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 1, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, "(1,3)(2,4)(5,6)", 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, gips_perm("(1,3)(2,4)(5,6)", 7), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), NULL, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_silent(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 1.1, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 0.9, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 1, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S) + 1), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, 7, FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, matrix(1:30, nrow = ncol(S)), FALSE, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), "FALSE", FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, "FALSE", FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, "FALSE", FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, FALSE, FALSE, "FALSE" )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, diag(nrow = ncol(S)), FALSE, TRUE, FALSE, FALSE )) # return_probabilities can be TRUE only when save_all_perms is also TRUE expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, matrix(1:30, nrow = ncol(S)), NA, FALSE, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, matrix(1:30, nrow = ncol(S)), FALSE, NA, FALSE, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, matrix(1:30, nrow = ncol(S)), FALSE, FALSE, NA, FALSE )) expect_error(check_correctness_of_arguments( S, number_of_observations, 30, permutations::permutation("(1,3)(2,4)(5,6)"), 3, matrix(1:30, nrow = ncol(S)), FALSE, FALSE, FALSE, NA )) # A lot of problems at once expect_error(check_correctness_of_arguments( S, number_of_observations + 0.1, 1, "(1,3)(2,4)(5,6)", 1, diag(nrow = ncol(S)), "FALSE", "FALSE", "FALSE", "FALSE" ), "8 problems identified with the provided arguments") # old tests: # A single problem at the same time: expect_error(check_correctness_of_arguments( S = NULL, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments( S = matrix(c(1:30), nrow = 6), number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations = NULL, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations = 0, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations = 15.5, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10.5, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 0, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 1, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = NULL, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = NULL, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 1, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = 7, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = matrix(c(1:30), nrow = 6), was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = "FALSE", return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = NA, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = 7, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = NULL, save_all_perms = FALSE, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = "FALSE", show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = NA, show_progress_bar = FALSE )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = 7 )) expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations, max_iter = 10, start_perm = example_perm, delta = 3, D_matrix = NULL, was_mean_estimated = FALSE, return_probabilities = FALSE, save_all_perms = FALSE, show_progress_bar = NULL )) # A number of problems at the same time. Not all are printed: expect_error(check_correctness_of_arguments(matrix_invariant_by_example_perm, number_of_observations = -1, max_iter = 1, start_perm = example_perm, delta = 1, D_matrix = 7, was_mean_estimated = NA, return_probabilities = "FALSE", save_all_perms = 7, show_progress_bar = NULL ), "\\.\\.\\. and 3 more problems") }) test_that("print.gips() works", { expect_identical( convert_log_diff_to_str(1009.5, 3), "2.632e+438" ) expect_identical( convert_log_diff_to_str(16.1, 3), "9820670.922" ) expect_identical( convert_log_diff_to_str(16.2, 3), "1.085e+7" ) expect_identical( convert_log_diff_to_str(-7.677, 3), "4.634e-4" ) expect_identical( convert_log_diff_to_str(Inf, 3), "Inf" ) expect_identical( convert_log_diff_to_str(-Inf, 3), "-Inf" ) expect_identical( convert_log_diff_to_str(0, 3), "1" ) g <- gips(S, number_of_observations, was_mean_estimated = FALSE) g_map <- find_MAP(g, 10, show_progress_bar = FALSE, optimizer = "MH") expect_output( print(g), "The permutation \\(\\):\n - is 1 times more likely than the \\(\\) permutation" ) expect_output( print(g, log_value = TRUE), "The permutation \\(\\):\n - is 1 times more likely than the \\(\\) permutation;\n - has log posteriori" ) expect_output( print(g_map), "\n - was found after 10 posteriori calculations;\n - is" ) # oneline: expect_output( print(g, oneline = TRUE), "The permutation \\(\\): is 1 times more likely" ) expect_output( print(g, oneline = TRUE, log_value = TRUE), "permutation; has log posteriori " ) expect_output( print(g_map, oneline = TRUE), ": was found after 10 posteriori calculations; is" ) }) test_that("plot.gips() works or abords for wrong arguments", { custom_perm1 <- gips_perm("(1,2)(3,4,5,6)", 6) g1 <- gips(S, number_of_observations, was_mean_estimated = FALSE, perm = custom_perm1 ) expect_error(plot.gips(custom_perm1)) expect_error(plot(g1, type = "both")) expect_error(plot(g1, type = c("du", "pa"))) expect_message( plot(g1), "`type = NA` was automatically changed to `type = 'heatmap'`" ) expect_silent(my_ggplot1 <- plot(g1, type = "heatmap")) expect_silent(my_ggplot2 <- plot(g1, type = "MLE")) expect_true(all.equal(my_ggplot1, my_ggplot2)) # cannot use expect_equal(), because it checks for equal enviroments, but in R all enviroments are different even if have the same elements in itself g1_found <- find_MAP(g1, 3, show_progress_bar = FALSE, optimizer = "MH") expect_message( plot(g1_found), "`type = NA` was automatically changed to `type = 'both'`" ) expect_silent(plot(g1_found, type = "both")) expect_silent(plot(g1_found, type = "all", logarithmic_y = FALSE)) expect_silent(plot(g1_found, type = "best")) expect_silent(plot(g1_found, type = "best", logarithmic_x = TRUE, ylim = range(attr(g1_found, "optimization_info")["log_posteriori_values"]) + c(-10, 10) )) expect_error( plot.gips(g1_found, type = "non_existing"), "Did You misspell the 'type' argument?" ) }) test_that("plot.gips() works for books example", { Z <- DAAG::oddbooks[, c(1, 2, 3)] Z$height <- Z$height / sqrt(2) g <- gips(cov(Z), 7, D_matrix = 1 * diag(3)) expect_silent(plot(g, type = "heatmap")) }) test_that("get_diagonalized_matrix_for_heatmap() works", { custom_perm1 <- gips_perm("(1,2)(3,4,5)(6)", 6) g1 <- gips(S, number_of_observations, was_mean_estimated = FALSE, perm = custom_perm1 ) actual <- get_diagonalized_matrix_for_heatmap(g1) # block_ends: 3,5,6 actual[1:3, 1:3] <- NA actual[4:5, 4:5] <- NA actual[6, 6] <- NA expect_true(all(is.na(actual))) }) test_that("summary.gips() works", { p <- 6 custom_perm1 <- gips_perm("(1,2)(3,4,5,6)", p) g1 <- gips(S, number_of_observations, was_mean_estimated = FALSE, perm = custom_perm1, D_matrix = diag(1, p) ) start_permutation_log_posteriori <- log_posteriori_of_gips(g1) log_posteriori_id <- log_posteriori_of_perm( perm_proposal = "", S = S, number_of_observations = number_of_observations, delta = attr(g1, "delta"), D_matrix = attr(g1, "D_matrix") ) likelihood_ratio_test_statistics <- 13*(determinant(project_matrix(S, custom_perm1))$modulus - determinant(S)$modulus) attributes(likelihood_ratio_test_statistics) <- NULL df_chisq <- p*(p+1)/2 - sum(get_structure_constants(custom_perm1)[["dim_omega"]]) likelihood_ratio_test_p_value <- 1 - pchisq(likelihood_ratio_test_statistics, df_chisq) my_sum <- structure(list( optimized = FALSE, start_permutation = structure(list( c(1, 2), c(3, 4, 5, 6) ), size = 6, class = "gips_perm"), start_permutation_log_posteriori = start_permutation_log_posteriori, times_more_likely_than_id = exp(start_permutation_log_posteriori - log_posteriori_id), log_times_more_likely_than_id = start_permutation_log_posteriori - log_posteriori_id, likelihood_ratio_test_statistics = likelihood_ratio_test_statistics, likelihood_ratio_test_p_value = likelihood_ratio_test_p_value, n0 = 2, S_matrix = S, number_of_observations = 13, was_mean_estimated = FALSE, delta = 3, D_matrix = structure(c( 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1 ), .Dim = c(6L, 6L)), n_parameters = 7, AIC = AIC(g1), BIC = BIC(g1) ), class = "summary.gips") expect_identical(summary(g1), my_sum) # all in my_sum is calculated by the same functions the summary should use expect_output( print(summary(g1)), "The number of observations is bigger than n0 for this permutation,\nso " ) expect_output( print(summary(g1)), "free parameters" ) expect_output( print(summary(g1)), "AIC" ) expect_output( print(summary(g1)), "BIC" ) g2 <- find_MAP(g1, max_iter = 10, optimizer = "MH", show_progress_bar = FALSE) expect_output( print(summary(g2)), "Optimization time:" ) g3 <- find_MAP(g2, max_iter = 10, optimizer = "continue", show_progress_bar = FALSE ) expect_output( print(summary(g3)), "Optimization algorithms:" ) # Optimized with BF; those 3 are computed differently for optimized with BF: g4 <- find_MAP(gips(S[1:4, 1:4], number_of_observations), optimizer = "BF", show_progress_bar = FALSE) expect_true(is.null(summary(g4)[["when_was_best"]])) expect_true(is.null(summary(g4)[["log_posteriori_calls_after_best"]])) expect_s3_class(summary(g4)[["start_permutation"]], "gips_perm") # proper log_posteriori_calls_after_best: # g_fake started in (1,2,3) and then comed to (1,3,2), which had bigger log_posteriori value, but only a little bigger g_fake <- structure(list(structure(list(c(1, 2, 3)), size = 3L, class = "gips_perm")), S = structure(c( 1, 0, 0, 0, 1, 0, 0, 0, 1 ), dim = c(3L, 3L)), number_of_observations = 10, delta = 3, D_matrix = structure(c( 1, 0, 0, 0, 1, 0, 0, 0, 1 ), dim = c(3L, 3L)), was_mean_estimated = TRUE, optimization_info = list( original_perm = gips_perm("", 3), acceptance_rate = 0.333333333333333, log_posteriori_values = c( -16.01209771488621000000, -18.9125198086359, -16.01209771488621421085 ), visited_perms = list( structure(list(c(1, 2, 3)), size = 3L, class = "gips_perm"), structure(list(1, c(2, 3)), size = 3L, class = "gips_perm"), structure(list(c(1, 3, 2)), size = 3L, class = "gips_perm") ), start_perm = structure(list(c(1, 2, 3)), size = 3L, class = "gips_perm"), last_perm = structure(list(c(1, 3, 2)), size = 3L, class = "gips_perm"), last_perm_log_posteriori = -16.0120977148862, iterations_performed = 2L, optimization_algorithm_used = "Metropolis_Hastings", post_probabilities = NULL, did_converge = NULL, best_perm_log_posteriori = -16.0120977148862, optimization_time = structure(0.00564193725585938, class = "difftime", units = "secs"), whole_optimization_time = structure(0.00564193725585938, class = "difftime", units = "secs"), all_n0 = c(2, 3, 2) ), class = "gips") expect_equal( summary(g_fake)[["log_posteriori_calls_after_best"]], 2 ) }) test_that("start_permutation_log_posteriori was calculated correctly", { g <- gips(S, number_of_observations, was_mean_estimated = TRUE) g_map <- find_MAP(g, max_iter = 10, show_progress_bar = FALSE, optimizer = "MH" ) optimization_info <- attr(g_map, "optimization_info") expect_equal( optimization_info[["log_posteriori_values"]][1], log_posteriori_of_perm( "", attr(g_map, "S"), attr(g_map, "number_of_observations") - attr(g_map, "was_mean_estimated"), attr(g_map, "delta"), attr(g_map, "D_matrix") ) ) }) test_that("summary.gips() returns proper n0 for estimated and unestimated mean", { g_no_em <- gips(S, number_of_observations, was_mean_estimated = FALSE) g_em <- gips(S, number_of_observations, was_mean_estimated = TRUE) expect_equal(summary.gips(g_no_em)[["n0"]], ncol(S)) # for known mean and perm id, one needs n >= p expect_equal(summary.gips(g_em)[["n0"]], ncol(S) + 1) # for estimated mean and perm id, one needs n >= (p + 1) }) test_that("summary.gips() returns proper Times more likely than identity permutation", { g_no_em <- gips(S, number_of_observations, was_mean_estimated = FALSE) g_em <- gips(S, number_of_observations, was_mean_estimated = TRUE) expect_equal( summary.gips(g_no_em)[["times_more_likely_than_id"]], 1 # for known mean ) expect_equal( summary.gips(g_em)[["times_more_likely_than_id"]], 1 # for estimated mean ) }) test_that("print.summary.gips() will properly print BIC and AIC when n < n0", { p <- 6 number_of_observations <- 2 Z <- MASS::mvrnorm(number_of_observations, mu = mu, Sigma = sigma_matrix) # this place is random S <- (t(Z) %*% Z) / number_of_observations # the theoretical mean is 0 g <- gips(S, number_of_observations, perm = "()") expect_silent(my_sum <- summary.gips(g)) expect_null(my_sum$AIC) expect_null(my_sum$BIC) expect_output(print(my_sum), "BIC:\n The number of observations is smaller than n0 for this permutation,\n so the gips model based on the found permutation does not exist.") expect_output(print(my_sum), "AIC:\n The number of observations is smaller than n0 for this permutation,\n so the gips model based on the found permutation does not exist.") }) test_that("get_probabilities_from_gips() works", { g <- gips(matrix(c(1, 0.5, 0.5, 1.3), nrow = 2), 13, was_mean_estimated = FALSE) g_map <- find_MAP(g, optimizer = "BF", show_progress_bar = FALSE, return_probabilities = TRUE, save_all_perms = TRUE ) expect_silent(probs <- get_probabilities_from_gips(g_map)) expect_type(probs, "double") expect_named(probs) expect_error( get_probabilities_from_gips(g), "Did You forget to optimize `g`?" ) g_map_no_prob <- find_MAP(g, optimizer = "BF", show_progress_bar = FALSE, return_probabilities = FALSE, save_all_perms = TRUE ) expect_message( out <- get_probabilities_from_gips(g_map_no_prob), "on the `gips` object that does not have saved probabilities." ) expect_null(out) g_map_no_prob_no_save <- find_MAP(g, optimizer = "BF", show_progress_bar = FALSE, return_probabilities = FALSE, save_all_perms = FALSE ) expect_message( out <- get_probabilities_from_gips(g_map_no_prob_no_save), "on the `gips` object that does not have saved probabilities." ) expect_null(out) # sorted probs <- get_probabilities_from_gips(g_map) expect_equal(order(probs, decreasing = TRUE), c(1, 2)) }) test_that("forget_perms() works properly", { g <- gips(S, number_of_observations, was_mean_estimated = FALSE) g_map <- find_MAP(g, 10, show_progress_bar = FALSE, optimizer = "MH", save_all_perms = TRUE ) expect_silent(validate_gips(g_map)) expect_message( forget_perms(g), "Provided \\`g\\` is a \\`gips\\` object, but it was not optimized yet\\." ) expect_false(all(is.na(attr(g_map, "optimization_info")[["visited_perms"]]))) expect_silent(g_map_no_perms <- forget_perms(g_map)) expect_true(all(is.na(attr(g_map_no_perms, "optimization_info")[["visited_perms"]]))) expect_silent(validate_gips(g_map_no_perms)) expect_message( g_map_no_perms_again <- forget_perms(g_map_no_perms), "Provided \\`g\\` is an optimized \\`gips\\` object that already has forgotten all permutations\\." ) }) test_that("logLik.gips() works", { # logLik calculated by hand: Z <- matrix(c( -1, 2, -3, 4, 1, -1, 1, 1, -1, 2, 2, 3, 2, -3, 2, -3, 3, -2, 2, -3 ), nrow = 5, byrow = TRUE) p <- ncol(Z) # 4 n <- nrow(Z) # 5 # ================== # Mean is (0,0,0,0) U <- t(Z) %*% Z perm <- gips_perm("(12)(34)", 4) S <- project_matrix(U, perm) / n # logLik from definition: loglikelihoods <- mvtnorm::dmvnorm( Z, rep(0, 4), S, log = TRUE ) logLik_definition <- sum(loglikelihoods) expect_equal(logLik_definition, -35.2883973048347) attr(logLik_definition, "df") <- 6 # (choose(p, 2) + p) - 4 parameters; choose(p, 2) + p is standard CoV; 4 is how much equalities is with (12)(34) attr(logLik_definition, "nobs") <- n class(logLik_definition) <- "logLik" # logLik.gips: expect_equal( logLik(gips(U / n, n, perm = perm, was_mean_estimated = FALSE)), logLik_definition ) # ================== # mean was estimated U <- cov(Z) * (n - 1) perm <- gips_perm("(12)(34)", 4) S <- project_matrix(U, perm) / (n - 1) logLik_expected <- structure(-28.8015774226105, df = 6, nobs = 5L, class = "logLik") # logLik.gips: expect_equal( logLik(gips(U / (n - 1), n, perm = perm, was_mean_estimated = TRUE)), logLik_expected ) # ================== # NUll: U <- t(Z) %*% Z expect_warning( expect_equal( logLik(gips(U / n, 2, perm = perm)), NULL ), class = "likelihood_does_not_exists" ) # ================== # -Inf: g <- gips(diag(0, 4), n) expect_warning( expect_equal(logLik(g), -Inf), class = "singular_matrix" ) }) test_that("AIC.gips() works", { Z <- matrix(c( -1, 2, -3, 4, 1, -1, 1, 1, -1, 2, 2, 3, 2, -3, 2, -3, 3, -2, 2, -3 ), nrow = 5, byrow = TRUE) p <- ncol(Z) # 4 n <- nrow(Z) # 5 g <- gips(t(Z) %*% Z / n, n, perm = "(12)(34)", was_mean_estimated = FALSE) expect_equal(AIC(g), 82.5767946096694) expect_equal(BIC(g), 80.233422084274) # ================== # NUll S <- t(Z) %*% Z / n g <- gips(S, 2) expect_warning( expect_equal( AIC(g), NULL ), class = "likelihood_does_not_exists" ) expect_warning( expect_equal( BIC(g), NULL ), class = "likelihood_does_not_exists" ) # ================== # Inf g <- gips(diag(0, 4), n) expect_warning(expect_equal(AIC(g), Inf), class = "singular_matrix") expect_warning(expect_equal(BIC(g), Inf), class = "singular_matrix") }) test_that("as.character.gips() work", { A <- matrix(rnorm(4 * 4), nrow = 4) S <- t(A) %*% A g <- gips(S, 14, perm = "(123)") expect_equal(as.character(g), "(1,2,3)") }) test_that("print.gips() will show the original perm for BF", { g <- gips(S[1:4,1:4], number_of_observations, perm = "(123)") g_map <- find_MAP(g, optimizer = "BF", show_progress_bar = FALSE) expect_output(print(g_map), "than the \\(1,2,3\\)") }) test_that("print.summary.gips() will not compare with original the unoptimized gips that is in id", { g <- gips(S[1:4,1:4], number_of_observations, perm = "(123)") expect_output(print(summary(g)), "Times more likely than identity permutation", fixed = TRUE) g <- gips(S[1:4,1:4], number_of_observations, perm = "()") g_map <- find_MAP(g, optimizer = "BF", show_progress_bar = FALSE) expect_output(print(summary(g_map)), "Times more likely than starting permutation", fixed = TRUE) pattern <- "Log_posteriori:\\s*(-?\\d+\\.\\d+)\\s\\sThe current permutation is id" expect_output(print(summary(g)), pattern) # The "Times more likely than starting permutation:" is skipped })