diff --git a/R/games.R b/R/games.R
index 27ddea4eb8b..29ba7631a31 100644
--- a/R/games.R
+++ b/R/games.R
@@ -1,4 +1,3 @@
-
 #' The Watts-Strogatz small-world model
 #'
 #' @description
@@ -480,7 +479,7 @@ sample_pa <- function(n, power = 1, m = NULL, out.dist = NULL, out.seq = NULL,
                       ),
                       start.graph = NULL) {
   if (!is.null(start.graph) && !is_igraph(start.graph)) {
-    stop("`start.graph' not an `igraph' object")
+    cli::cli_abort("{.arg start.graph} must be an {.cls igraph} object, not {.obj_type_friendly {start.graph}}.")
   }
 
   # Checks
@@ -727,7 +726,7 @@ erdos.renyi.game <- function(n, p.or.m, type = c("gnp", "gnm"),
 #' @family games
 #' @export
 random.graph.game <- function(n, p.or.m, type = c("gnp", "gnm"),
-                             directed = FALSE, loops = FALSE) {
+                              directed = FALSE, loops = FALSE) {
   type <- igraph.match.arg(type)
 
   if (type == "gnp") {
@@ -1124,16 +1123,16 @@ sample_pa_age <- function(n, pa.exp, aging.exp, m = NULL, aging.bin = 300,
     m <- NULL
   }
   if (!is.null(out.seq) && length(out.seq) != n) {
-    stop("`out.seq' should be of length `n'")
+    cli::cli_abort("{.arg out.seq} must have length {.val n}, not {.val {length( out.seq)}}.'")
   }
   if (!is.null(out.seq) && min(out.seq) < 0) {
-    stop("negative elements in `out.seq'")
+    cli::cli_abort("{.arg out.seq} must not contain negative elements.")
   }
   if (!is.null(m) && m < 0) {
-    stop("`m' is negative")
+    cli::cli_abort("{.arg m} must be positive or 0.")
   }
   if (!is.null(time.window) && time.window <= 0) {
-    stop("time window size should be positive")
+    cli::cli_abort("{.arg time.window} must be positive.")
   }
   if (!is.null(m) && m == 0) {
     cli::cli_warn("{.arg m} is zero, graph will be empty.")
@@ -1425,7 +1424,10 @@ sample_pref <- function(nodes, types, type.dist = rep(1, types),
                         pref.matrix = matrix(1, types, types),
                         directed = FALSE, loops = FALSE) {
   if (nrow(pref.matrix) != types || ncol(pref.matrix) != types) {
-    stop("Invalid size for preference matrix")
+    cli::cli_abort(c(
+      "{.arg pref.matrix} must have {.arg types} rows and columns.",
+      i = "See {.fun igraph::sample_pref}'s manual."
+    ))
   }
 
   on.exit(.Call(R_igraph_finalizer))
@@ -1460,10 +1462,16 @@ sample_asym_pref <- function(nodes, types,
                              pref.matrix = matrix(1, types, types),
                              loops = FALSE) {
   if (nrow(pref.matrix) != types || ncol(pref.matrix) != types) {
-    stop("Invalid size for preference matrix")
+    cli::cli_abort(c(
+      "{.arg pref.matrix} must have {.arg types} rows and columns.",
+      i = "See {.fun igraph::sample_asym_pref}'s manual."
+    ))
   }
   if (nrow(type.dist.matrix) != types || ncol(type.dist.matrix) != types) {
-    stop("Invalid size for type distribution matrix")
+    cli::cli_abort(c(
+      "{.arg type.dist.matrix} must have {.arg types} rows and columns.",
+      i = "See {.fun igraph::sample_asym_pref}'s manual."
+    ))
   }
 
   on.exit(.Call(R_igraph_finalizer))
@@ -1758,16 +1766,16 @@ sample_bipartite <- function(n1, n2, type = c("gnp", "gnm"), p, m,
   )
 
   if (type == "gnp" && missing(p)) {
-    stop("Connection probability `p' is not given for Gnp graph")
+    cli::cli_abort("Connection probability {.arg p} must be given for Gnp graph")
   }
   if (type == "gnp" && !missing(m)) {
     cli::cli_warn("Number of edges {.arg m} is ignored for Gnp graph.")
   }
   if (type == "gnm" && missing(m)) {
-    stop("Number of edges `m' is not given for Gnm graph")
+    cli::cli_abort("Number of edges {.arg m} must be given for Gnm graph")
   }
   if (type == "gnm" && !missing(p)) {
-    cli::cli_warn("Connection probability {.arg p} is ignored for Gnp graph.")
+    cli::cli_warn("Connection probability {.arg p} is ignored for Gnm graph.")
   }
 
   on.exit(.Call(R_igraph_finalizer))
@@ -1888,7 +1896,7 @@ sample_hierarchical_sbm <- function(n, m, rho, C, p) {
   } else {
     commonlen <- setdiff(commonlen, 1)
     if (length(commonlen) != 1) {
-      stop("Lengths of `m', `rho' and `C' must match")
+      cli::cli_abort("Lengths of {.arg m}, {.arg rho} and {.arg C} must match.")
     }
     m <- rep(m, length.out = commonlen)
     rho <- if (is.list(rho)) {
@@ -2165,7 +2173,7 @@ sample_k_regular <- k_regular_game_impl
 #'
 #' rowMeans(replicate(
 #'   100,
-#'   degree(sample_chung_lu(c(1, 3, 2, 1), c(2, 1, 2, 2), variant = "maxent"), mode='out')
+#'   degree(sample_chung_lu(c(1, 3, 2, 1), c(2, 1, 2, 2), variant = "maxent"), mode = "out")
 #' ))
 #' @export
 #' @cdocs igraph_chung_lu_game
@@ -2178,8 +2186,7 @@ chung_lu <- function(
     in.weights = NULL,
     ...,
     loops = TRUE,
-    variant = c("original", "maxent", "nr")
-) {
+    variant = c("original", "maxent", "nr")) {
   variant <- rlang::arg_match(variant)
   constructor_spec(
     sample_chung_lu,
diff --git a/tests/testthat/test-ba.game.R b/tests/testthat/test-ba.game.R
deleted file mode 100644
index e66d45a3e00..00000000000
--- a/tests/testthat/test-ba.game.R
+++ /dev/null
@@ -1,76 +0,0 @@
-test_that("sample_pa() works", {
-  withr::local_seed(20240209)
-
-  g <- sample_pa(100, m = 2)
-  expect_equal(ecount(g), 197)
-  expect_equal(vcount(g), 100)
-  expect_true(is_simple(g))
-
-  g2 <- sample_pa(100, m = 2, algorithm = "psumtree-multiple")
-  expect_equal(ecount(g2), 198)
-  expect_equal(vcount(g2), 100)
-  expect_false(is_simple(g2))
-
-  g3 <- sample_pa(100, m = 2, algorithm = "bag")
-  expect_equal(ecount(g3), 198)
-  expect_equal(vcount(g3), 100)
-  expect_false(is_simple(g3))
-
-  g4 <- sample_pa(3, out.seq = 0:2, directed = FALSE)
-  expect_equal(degree(g4), rep(2, 3))
-
-  g5 <- sample_pa(3, out.dist = rep(2, 1000), directed = FALSE)
-  expect_equal(degree(g5), rep(2, 3))
-})
-
-test_that("sample_pa can start from a graph", {
-  withr::local_seed(20231029)
-
-  g4 <- sample_pa(10, m = 1, algorithm = "bag", start.graph = make_empty_graph(5))
-  expect_equal(ecount(g4), 5)
-  expect_equal(vcount(g4), 10)
-  #  0    1    2    3    4
-  # 24  809 3904 4240 1023
-  is_degree_zero <- (degree(g4) == 0)
-  expect_true(sum(is_degree_zero) %in% 0:4)
-  #  2    3    4    5    6    7    8   10
-  # 25  302 1820 2563 3350 1093  816   31
-  is_degree_one <- (degree(g4) == 1)
-  expect_true(sum(is_degree_one) %in% c(2:8, 10L))
-  #   0    1    2    3    4
-  # 879 2271 5289 1532   29
-  is_degree_two_or_three <- (degree(g4) %in% 2:3)
-  expect_true(sum(is_degree_two_or_three) %in% 0:4)
-
-  g6 <- sample_pa(10, m = 1, algorithm = "bag", start.graph = make_star(10))
-  expect_isomorphic(g6, make_star(10))
-
-  g7 <- sample_pa(10,
-    m = 3, algorithm = "psumtree-multiple",
-    start.graph = make_empty_graph(5)
-  )
-  expect_equal(degree(g7, mode = "out"), c(0, 0, 0, 0, 0, 3, 3, 3, 3, 3))
-
-  g8 <- sample_pa(10,
-    m = 3, algorithm = "psumtree-multiple",
-    start.graph = make_star(5)
-  )
-  expect_equal(degree(g8, mode = "out"), c(0, 1, 1, 1, 1, 3, 3, 3, 3, 3))
-  expect_isomorphic(induced_subgraph(g8, 1:5), make_star(5))
-
-  g9 <- sample_pa(10,
-    m = 3, algorithm = "psumtree-multiple",
-    start.graph = make_star(10)
-  )
-  expect_isomorphic(g9, make_star(10))
-
-  g10 <- sample_pa(10, m = 3, start.graph = make_empty_graph(5))
-  expect_equal(degree(g10, mode = "out"), c(0, 0, 0, 0, 0, 3, 3, 3, 3, 3))
-
-  g11 <- sample_pa(10, m = 3, start.graph = make_star(5))
-  expect_equal(degree(g11, mode = "out"), c(0, 1, 1, 1, 1, 3, 3, 3, 3, 3))
-  expect_isomorphic(induced_subgraph(g11, 1:5), make_star(5))
-
-  g12 <- sample_pa(10, m = 3, start.graph = make_star(10))
-  expect_isomorphic(g12, make_star(10))
-})
diff --git a/tests/testthat/test-bipartite.random.game.R b/tests/testthat/test-bipartite.random.game.R
deleted file mode 100644
index 9d2df975e32..00000000000
--- a/tests/testthat/test-bipartite.random.game.R
+++ /dev/null
@@ -1,53 +0,0 @@
-test_that("sample_bipartite works", {
-  withr::local_seed(42)
-  g1 <- sample_bipartite(10, 5, type = "gnp", p = .1)
-  expect_equal(g1$name, "Bipartite Gnp random graph")
-  expect_equal(vcount(g1), 15)
-  expect_equal(ecount(g1), 7)
-  expect_true(bipartite_mapping(g1)$res)
-  expect_false(is_directed(g1))
-
-  g2 <- sample_bipartite(10, 5, type = "gnp", p = .1, directed = TRUE)
-  expect_equal(vcount(g2), 15)
-  expect_equal(ecount(g2), 6)
-  expect_true(bipartite_mapping(g2)$res)
-  expect_true(is_directed(g2))
-  expect_output(print_all(g2), "5->11")
-
-  g3 <- sample_bipartite(10, 5, type = "gnp", p = .1, directed = TRUE, mode = "in")
-  expect_output(print_all(g3), "11->3")
-
-  g4 <- sample_bipartite(10, 5, type = "gnm", m = 8)
-  expect_equal(vcount(g4), 15)
-  expect_equal(ecount(g4), 8)
-  expect_true(bipartite_mapping(g4)$res)
-  expect_false(is_directed(g4))
-
-  g5 <- sample_bipartite(10, 5, type = "gnm", m = 8, directed = TRUE)
-  expect_equal(vcount(g5), 15)
-  expect_equal(ecount(g5), 8)
-  expect_true(bipartite_mapping(g5)$res)
-  expect_true(is_directed(g5))
-  expect_output(print_all(g5), "5->12")
-
-  g6 <- sample_bipartite(10, 5, type = "gnm", m = 8, directed = TRUE, mode = "in")
-  expect_equal(vcount(g6), 15)
-  expect_equal(ecount(g6), 8)
-  expect_true(bipartite_mapping(g6)$res)
-  expect_true(is_directed(g6))
-  expect_output(print_all(g6), "12->10")
-
-  #####
-
-  g7 <- sample_bipartite(10, 5,
-    type = "gnp", p = 0.9999, directed = TRUE,
-    mode = "all"
-  )
-  expect_equal(ecount(g7), 100)
-
-  g8 <- sample_bipartite(10, 5,
-    type = "gnm", m = 99, directed = TRUE,
-    mode = "all"
-  )
-  expect_equal(ecount(g8), 99)
-})
diff --git a/tests/testthat/test-chung_lu.R b/tests/testthat/test-chung_lu.R
deleted file mode 100644
index 3070de76353..00000000000
--- a/tests/testthat/test-chung_lu.R
+++ /dev/null
@@ -1,13 +0,0 @@
-test_that("sample_chung_lu works", {
-  g <- sample_chung_lu(c(3, 3, 2, 2, 1, 1))
-  expect_false(any_multiple(g))
-
-  g <- sample_chung_lu(c(3, 3, 2, 2, 1, 1), loops = FALSE, variant = 'original')
-  expect_true(is_simple(g))
-
-  g <- sample_chung_lu(c(3, 3, 2, 2, 1, 1), loops = FALSE, variant = 'maxent')
-  expect_true(is_simple(g))
-
-  g <- sample_chung_lu(c(3, 3, 2, 2, 1, 1), loops = FALSE, variant = 'nr')
-  expect_true(is_simple(g))
-})
diff --git a/tests/testthat/test-correlated.R b/tests/testthat/test-correlated.R
deleted file mode 100644
index aa579b0cbf4..00000000000
--- a/tests/testthat/test-correlated.R
+++ /dev/null
@@ -1,88 +0,0 @@
-## Not very meaningful tests. They good for testing that the
-## functions run, but not much more
-
-test_that("sample_correlated_gnp works", {
-  withr::local_seed(42)
-
-  g <- sample_gnp(10, .1)
-  g2 <- sample_correlated_gnp(g, corr = 1, p = g$p, permutation = NULL)
-  expect_equal(g[], g2[])
-
-  g3 <- sample_correlated_gnp(g, corr = 0, p = g$p, permutation = NULL)
-  c3 <- cor(as.vector(g[]), as.vector(g3[]))
-  expect_true(abs(c3) < .3)
-})
-
-test_that("sample_correlated_gnp works when p is not given", {
-  withr::local_seed(42)
-
-  g <- sample_gnp(10, .1)
-  g2 <- sample_correlated_gnp(g, corr = 1)
-  expect_equal(g[], g2[])
-
-  g3 <- sample_correlated_gnp(g, corr = 0)
-  c3 <- cor(as.vector(g[]), as.vector(g3[]))
-  expect_true(abs(c3) < .3)
-})
-
-test_that("sample_correlated_gnp works even for non-ER graphs", {
-  withr::local_seed(42)
-
-  g <- sample_grg(100, 0.2)
-  g2 <- sample_correlated_gnp(g, corr = 1)
-  expect_equal(g[], g2[])
-
-  g3 <- sample_correlated_gnp(g, corr = 0)
-  c3 <- cor(as.vector(g[]), as.vector(g3[]))
-  expect_true(abs(c3) < .3)
-})
-
-test_that("sample_correlated_gnp_pair works", {
-  withr::local_seed(42)
-
-  gp <- sample_correlated_gnp_pair(10, corr = .95, p = .1, permutation = NULL)
-  expect_true(abs(ecount(gp[[1]]) - ecount(gp[[2]])) < 3)
-})
-
-## Some corner cases
-
-test_that("sample_correlated_gnp corner cases work", {
-  withr::local_seed(42)
-
-  is.full <- function(g) {
-    g2 <- make_full_graph(vcount(g), directed = is_directed(g))
-    isomorphic(g, g2)
-  }
-
-  g <- sample_gnp(10, .3)
-  g2 <- sample_correlated_gnp(g, corr = 0.000001, p = .99999999)
-  expect_true(is.full(g2))
-
-  g3 <- sample_correlated_gnp(g, corr = 0.000001, p = 0.0000001)
-  expect_equal(ecount(g3), 0)
-  expect_equal(vcount(g3), 10)
-
-  gg <- sample_gnp(10, .3, directed = TRUE)
-  gg2 <- sample_correlated_gnp(gg, corr = 0.000001, p = .99999999)
-  expect_true(is.full(gg2))
-
-  gg3 <- sample_correlated_gnp(gg, corr = 0.000001, p = 0.0000001)
-  expect_equal(ecount(gg3), 0)
-  expect_equal(vcount(gg3), 10)
-})
-
-test_that("permutation works for sample_correlated_gnp", {
-  withr::local_seed(42)
-
-  g <- sample_gnp(10, .3)
-  perm <- sample(vcount(g))
-  g2 <- sample_correlated_gnp(g, corr = .99999, p = .3, permutation = perm)
-  g <- permute(g, perm)
-  expect_equal(g[], g2[])
-
-  g <- sample_gnp(10, .3)
-  perm <- sample(vcount(g))
-  g2 <- sample_correlated_gnp(g, corr = 1, p = .3, permutation = perm)
-  g <- permute(g, perm)
-  expect_equal(g[], g2[])
-})
diff --git a/tests/testthat/test-forestfire.R b/tests/testthat/test-forestfire.R
deleted file mode 100644
index cea841c71d3..00000000000
--- a/tests/testthat/test-forestfire.R
+++ /dev/null
@@ -1,34 +0,0 @@
-test_that("sample_forestfire() works -- sparse", {
-  withr::local_seed(20231029)
-  N <- 5000
-  xv <- log(2:N)
-
-  g1 <- sample_forestfire(N, fw.prob = 0.35, bw.factor = 0.2 / 0.35)
-  yv1 <- log(cumsum(degree(g1, mode = "out"))[-1])
-
-  expect_equal(coef(lm(yv1 ~ xv))[[2]], 1.04, tolerance = 0.05)
-})
-
-test_that("sample_forestfire() works -- densifying", {
-  withr::local_seed(20231029)
-
-  N <- 5000
-  xv <- log(2:N)
-
-  g2 <- sample_forestfire(N, fw.prob = 0.37, bw.factor = 0.32 / 0.37)
-  yv2 <- log(cumsum(degree(g2, mode = "out"))[-1])
-
-  expect_equal(coef(lm(yv2 ~ xv))[[2]], 1.21, tolerance = 0.05)
-})
-
-test_that("sample_forestfire() works -- dense", {
-  withr::local_seed(20231029)
-
-  N <- 5000
-  xv <- log(2:N)
-
-  g3 <- sample_forestfire(N, fw.prob = 0.38, bw.factor = 0.38 / 0.37)
-  yv3 <- log(cumsum(degree(g3, mode = "out"))[-1])
-
-  expect_equal(coef(lm(yv3 ~ xv))[[2]], 1.9, tolerance = 0.05)
-})
diff --git a/tests/testthat/test-games.R b/tests/testthat/test-games.R
index 9b0c9991ca9..fe602cbdbfe 100644
--- a/tests/testthat/test-games.R
+++ b/tests/testthat/test-games.R
@@ -140,7 +140,6 @@ test_that("sample_degseq supports the sample_(...) syntax", {
 })
 
 test_that("sample_degseq works() -- old method names", {
-
   withr::local_options("lifecycle_verbosity" = "warning")
 
   expect_warning(
@@ -158,3 +157,449 @@ test_that("sample_degseq works() -- old method names", {
     "must be"
   )
 })
+
+test_that("sample_chung_lu works", {
+  chung_lu_small <- sample_chung_lu(c(3, 3, 2, 2, 1, 1))
+  expect_false(any_multiple(chung_lu_small))
+
+  chung_lu_no_loop_1 <- sample_chung_lu(c(3, 3, 2, 2, 1, 1), loops = FALSE, variant = "original")
+  expect_true(is_simple(chung_lu_no_loop_1))
+
+  chung_lu_no_loop_2 <- sample_chung_lu(c(3, 3, 2, 2, 1, 1), loops = FALSE, variant = "maxent")
+  expect_true(is_simple(chung_lu_no_loop_2))
+
+  chung_lu_no_loop_3 <- sample_chung_lu(c(3, 3, 2, 2, 1, 1), loops = FALSE, variant = "nr")
+  expect_true(is_simple(chung_lu_no_loop_3))
+})
+
+test_that("sample_forestfire() works -- sparse", {
+  withr::local_seed(20231029)
+  N <- 5000
+  xv <- log(2:N)
+
+  forest_fire <- sample_forestfire(N, fw.prob = 0.35, bw.factor = 0.2 / 0.35)
+  yv1 <- log(cumsum(degree(forest_fire, mode = "out"))[-1])
+
+  expect_equal(coef(lm(yv1 ~ xv))[[2]], 1.04, tolerance = 0.05)
+})
+
+test_that("sample_forestfire() works -- densifying", {
+  withr::local_seed(20231029)
+
+  N <- 5000
+  xv <- log(2:N)
+
+  forest_fire <- sample_forestfire(N, fw.prob = 0.37, bw.factor = 0.32 / 0.37)
+  yv2 <- log(cumsum(degree(forest_fire, mode = "out"))[-1])
+
+  expect_equal(coef(lm(yv2 ~ xv))[[2]], 1.21, tolerance = 0.05)
+})
+
+test_that("sample_forestfire() works -- dense", {
+  withr::local_seed(20231029)
+
+  N <- 5000
+  xv <- log(2:N)
+
+  forest_fire <- sample_forestfire(N, fw.prob = 0.38, bw.factor = 0.38 / 0.37)
+  yv3 <- log(cumsum(degree(forest_fire, mode = "out"))[-1])
+
+  expect_equal(coef(lm(yv3 ~ xv))[[2]], 1.9, tolerance = 0.05)
+})
+
+test_that("Generating stochastic block models works", {
+  pm <- matrix(1, nrow = 2, ncol = 2)
+  bs <- c(4, 6)
+  sbm_small <- sample_sbm(10,
+    pref.matrix = pm, block.sizes = bs,
+    directed = FALSE, loops = FALSE
+  )
+  expect_isomorphic(sbm_small, make_full_graph(10, directed = FALSE, loops = FALSE))
+
+  sbm_small_loops <- sample_sbm(10,
+    pref.matrix = pm, block.sizes = bs,
+    directed = FALSE, loops = TRUE
+  )
+  full_graph_loops <- make_full_graph(10, directed = FALSE, loops = TRUE)
+  expect_equal(sbm_small_loops[sparse = FALSE], full_graph_loops[sparse = FALSE])
+
+  sbm_small_directed <- sample_sbm(10,
+    pref.matrix = pm, block.sizes = bs,
+    directed = TRUE, loops = FALSE
+  )
+  full_graph_directed <- make_full_graph(10, directed = TRUE, loops = FALSE)
+  expect_equal(sbm_small_directed[sparse = FALSE], full_graph_directed[sparse = FALSE])
+
+  sbm_small_all <- sample_sbm(10,
+    pref.matrix = pm, block.sizes = bs,
+    directed = TRUE, loops = TRUE
+  )
+  full_graph_all <- make_full_graph(10, directed = TRUE, loops = TRUE)
+  expect_equal(sbm_small_all[sparse = FALSE], full_graph_all[sparse = FALSE])
+})
+
+test_that("sample_smallworld works", {
+  for (i in 1:50) {
+    p <- runif(1)
+    d <- sample(1:3, 1)
+    nei <- sample(2:5, 1)
+    g <- sample_smallworld(d, 10, nei, p, loops = FALSE)
+    expect_false(any(which_loop(g)))
+  }
+})
+
+test_that("sample_pa() works", {
+  withr::local_seed(20240209)
+
+  g_pa <- sample_pa(100, m = 2)
+  expect_equal(ecount(g_pa), 197)
+  expect_equal(vcount(g_pa), 100)
+  expect_true(is_simple(g_pa))
+
+  g_pa2 <- sample_pa(100, m = 2, algorithm = "psumtree-multiple")
+  expect_equal(ecount(g_pa2), 198)
+  expect_equal(vcount(g_pa2), 100)
+  expect_false(is_simple(g_pa2))
+
+  g_pa3 <- sample_pa(100, m = 2, algorithm = "bag")
+  expect_equal(ecount(g_pa3), 198)
+  expect_equal(vcount(g_pa3), 100)
+  expect_false(is_simple(g_pa3))
+
+  g_pa4 <- sample_pa(3, out.seq = 0:2, directed = FALSE)
+  expect_equal(degree(g_pa4), rep(2, 3))
+
+  g_pa5 <- sample_pa(3, out.dist = rep(2, 1000), directed = FALSE)
+  expect_equal(degree(g_pa5), rep(2, 3))
+})
+
+test_that("sample_pa can start from a graph", {
+  withr::local_seed(20231029)
+
+  g_pa1 <- sample_pa(10, m = 1, algorithm = "bag", start.graph = make_empty_graph(5))
+  expect_equal(ecount(g_pa1), 5)
+  expect_equal(vcount(g_pa1), 10)
+
+  is_degree_zero <- (degree(g_pa1) == 0)
+  expect_true(sum(is_degree_zero) %in% 0:4)
+  #  2    3    4    5    6    7    8   10
+  # 25  302 1820 2563 3350 1093  816   31
+  is_degree_one <- (degree(g_pa1) == 1)
+  expect_true(sum(is_degree_one) %in% c(2:8, 10L))
+  #   0    1    2    3    4
+  # 879 2271 5289 1532   29
+  is_degree_two_or_three <- (degree(g_pa1) %in% 2:3)
+  expect_true(sum(is_degree_two_or_three) %in% 0:4)
+
+  g_pa2 <- sample_pa(10, m = 1, algorithm = "bag", start.graph = make_star(10))
+  expect_isomorphic(g_pa2, make_star(10))
+
+  g_pa3 <- sample_pa(10,
+    m = 3, algorithm = "psumtree-multiple",
+    start.graph = make_empty_graph(5)
+  )
+  expect_equal(degree(g_pa3, mode = "out"), c(0, 0, 0, 0, 0, 3, 3, 3, 3, 3))
+
+  g_pa4 <- sample_pa(10,
+    m = 3, algorithm = "psumtree-multiple",
+    start.graph = make_star(5)
+  )
+  expect_equal(degree(g_pa4, mode = "out"), c(0, 1, 1, 1, 1, 3, 3, 3, 3, 3))
+  expect_isomorphic(induced_subgraph(g_pa4, 1:5), make_star(5))
+
+  g_pa5 <- sample_pa(10,
+    m = 3, algorithm = "psumtree-multiple",
+    start.graph = make_star(10)
+  )
+  expect_isomorphic(g_pa5, make_star(10))
+
+  g_pa6 <- sample_pa(10, m = 3, start.graph = make_empty_graph(5))
+  expect_equal(degree(g_pa6, mode = "out"), c(0, 0, 0, 0, 0, 3, 3, 3, 3, 3))
+
+  g_pa7 <- sample_pa(10, m = 3, start.graph = make_star(5))
+  expect_equal(degree(g_pa7, mode = "out"), c(0, 1, 1, 1, 1, 3, 3, 3, 3, 3))
+  expect_isomorphic(induced_subgraph(g_pa7, 1:5), make_star(5))
+
+  g_pa8 <- sample_pa(10, m = 3, start.graph = make_star(10))
+  expect_isomorphic(g_pa8, make_star(10))
+})
+
+test_that("sample_bipartite works -- undirected gnp", {
+  withr::local_seed(42)
+
+  g_rand_bip <- sample_bipartite(10, 5, type = "gnp", p = .1)
+  expect_equal(g_rand_bip$name, "Bipartite Gnp random graph")
+  expect_equal(vcount(g_rand_bip), 15)
+  expect_equal(ecount(g_rand_bip), 7)
+  expect_true(bipartite_mapping(g_rand_bip)$res)
+  expect_false(is_directed(g_rand_bip))
+})
+
+test_that("sample_bipartite works -- directed gnp", {
+  g_rand_bip_dir <- sample_bipartite(10, 5, type = "gnp", p = .1, directed = TRUE)
+  expect_equal(vcount(g_rand_bip_dir), 15)
+  expect_equal(ecount(g_rand_bip_dir), 6)
+  expect_true(bipartite_mapping(g_rand_bip_dir)$res)
+  expect_true(is_directed(g_rand_bip_dir))
+  expect_output(print_all(g_rand_bip_dir), "5->11")
+
+  g_rand_bip_in <- sample_bipartite(10, 5, type = "gnp", p = .1, directed = TRUE, mode = "in")
+  expect_output(print_all(g_rand_bip_in), "11->3")
+})
+
+test_that("sample_bipartite works -- undirected gnm", {
+  g_rand_bip_gnm <- sample_bipartite(10, 5, type = "gnm", m = 8)
+  expect_equal(vcount(g_rand_bip_gnm), 15)
+  expect_equal(ecount(g_rand_bip_gnm), 8)
+  expect_true(bipartite_mapping(g_rand_bip_gnm)$res)
+  expect_false(is_directed(g_rand_bip_gnm))
+})
+test_that("sample_bipartite works -- directed gnm", {
+  g_rand_bip_gnm_dir <- sample_bipartite(10, 5, type = "gnm", m = 8, directed = TRUE)
+  expect_equal(vcount(g_rand_bip_gnm_dir), 15)
+  expect_equal(ecount(g_rand_bip_gnm_dir), 8)
+  expect_true(bipartite_mapping(g_rand_bip_gnm_dir)$res)
+  expect_true(is_directed(g_rand_bip_gnm_dir))
+  expect_output(print_all(g_rand_bip_gnm_dir), "5->12")
+
+  g_rand_bip_gnm_in <- sample_bipartite(10, 5, type = "gnm", m = 8, directed = TRUE, mode = "in")
+  expect_equal(vcount(g_rand_bip_gnm_in), 15)
+  expect_equal(ecount(g_rand_bip_gnm_in), 8)
+  expect_true(bipartite_mapping(g_rand_bip_gnm_in)$res)
+  expect_true(is_directed(g_rand_bip_gnm_in))
+  expect_output(print_all(g_rand_bip_gnm_in), "12->10")
+
+  g_rand_bip_full <- sample_bipartite(10, 5,
+    type = "gnp", p = 0.9999, directed = TRUE,
+    mode = "all"
+  )
+  expect_equal(ecount(g_rand_bip_full), 100)
+
+  g_rand_bip_edges <- sample_bipartite(10, 5,
+    type = "gnm", m = 99, directed = TRUE,
+    mode = "all"
+  )
+  expect_equal(ecount(g_rand_bip_edges), 99)
+})
+
+
+test_that("sample_correlated_gnp works", {
+  withr::local_seed(42)
+
+  gnp_graph <- sample_gnp(10, .1)
+  cor_gnp_graph_1 <- sample_correlated_gnp(gnp_graph, corr = 1, p = gnp_graph$p, permutation = NULL)
+  expect_equal(gnp_graph[], cor_gnp_graph_1[])
+
+  cor_gnp_graph_0 <- sample_correlated_gnp(gnp_graph, corr = 0, p = gnp_graph$p, permutation = NULL)
+  graph_cor_1 <- cor(as.vector(gnp_graph[]), as.vector(cor_gnp_graph_0[]))
+  expect_true(abs(graph_cor_1) < .3)
+
+  cor_gnp_no_p_1 <- sample_correlated_gnp(gnp_graph, corr = 1)
+  expect_equal(gnp_graph[], cor_gnp_no_p_1[])
+
+  cor_gnp_no_p_0 <- sample_correlated_gnp(gnp_graph, corr = 0)
+  graph_cor_2 <- cor(as.vector(gnp_graph[]), as.vector(cor_gnp_no_p_0[]))
+  expect_true(abs(graph_cor_2) < .3)
+})
+
+
+test_that("sample_correlated_gnp works even for non-ER graphs", {
+  withr::local_seed(42)
+
+  grg_graph <- sample_grg(100, 0.2)
+  cor_gnp_graph_1 <- sample_correlated_gnp(grg_graph, corr = 1)
+  expect_equal(grg_graph[], cor_gnp_graph_1[])
+
+  cor_gnp_graph_0 <- sample_correlated_gnp(grg_graph, corr = 0)
+  graph_cor <- cor(as.vector(grg_graph[]), as.vector(cor_gnp_graph_0[]))
+  expect_true(abs(graph_cor) < .3)
+})
+
+test_that("sample_correlated_gnp_pair works", {
+  withr::local_seed(42)
+
+  cor_gnp_pair <- sample_correlated_gnp_pair(10, corr = .95, p = .1, permutation = NULL)
+  expect_true(abs(ecount(cor_gnp_pair[[1]]) - ecount(cor_gnp_pair[[2]])) < 3)
+})
+
+## Some corner cases
+
+test_that("sample_correlated_gnp corner cases work", {
+  withr::local_seed(42)
+
+  is_full <- function(g) {
+    g_full <- make_full_graph(vcount(g), directed = is_directed(g))
+    isomorphic(g, g_full)
+  }
+
+  gnp_graph <- sample_gnp(10, .3)
+  cor_gnp_full <- sample_correlated_gnp(gnp_graph, corr = 0.000001, p = .99999999)
+  expect_true(is_full(cor_gnp_full))
+
+  cor_gnp_empty <- sample_correlated_gnp(gnp_graph, corr = 0.000001, p = 0.0000001)
+  expect_equal(ecount(cor_gnp_empty), 0)
+  expect_equal(vcount(cor_gnp_empty), 10)
+
+  gnp_graph_directed <- sample_gnp(10, .3, directed = TRUE)
+  cor_gnp_directed <- sample_correlated_gnp(gnp_graph_directed, corr = 0.000001, p = .99999999)
+  expect_true(is_full(cor_gnp_directed))
+
+  cor_gnp_directed_empty <- sample_correlated_gnp(gnp_graph_directed, corr = 0.000001, p = 0.0000001)
+  expect_equal(ecount(cor_gnp_directed_empty), 0)
+  expect_equal(vcount(cor_gnp_directed_empty), 10)
+})
+
+test_that("permutation works for sample_correlated_gnp", {
+  withr::local_seed(42)
+
+  gnp_graph <- sample_gnp(10, .3)
+  perm <- sample(vcount(gnp_graph))
+  cor_gnp_graph <- sample_correlated_gnp(gnp_graph, corr = .99999, p = .3, permutation = perm)
+  gnp_graph <- permute(gnp_graph, perm)
+  expect_equal(gnp_graph[], cor_gnp_graph[])
+
+  perm <- sample(vcount(gnp_graph))
+  cor_gnp_graph <- sample_correlated_gnp(gnp_graph, corr = 1, p = .3, permutation = perm)
+  gnp_graph <- permute(gnp_graph, perm)
+  expect_equal(gnp_graph[], cor_gnp_graph[])
+})
+
+test_that("HSBM works", {
+  withr::local_seed(42)
+
+  C <- matrix(c(
+    1, 1 / 2, 0,
+    1 / 2, 0, 1 / 2,
+    0, 1 / 2, 1 / 2
+  ), nrow = 3)
+
+  g_hsbm1 <- sample_hierarchical_sbm(100, 10, rho = c(3, 3, 4) / 10, C = C, p = 0)
+  expect_equal(ecount(g_hsbm1), 172)
+  expect_equal(vcount(g_hsbm1), 100)
+  expect_false(is_directed(g_hsbm1))
+
+  withr::local_seed(42)
+
+  g_hsbm2 <- sample_hierarchical_sbm(100, 10, rho = c(3, 3, 4) / 10, C = C, p = 1)
+  expect_equal(ecount(g_hsbm2), ecount(g_hsbm1) + 10 * 9 * (90 + 10) / 2)
+  expect_equal(vcount(g_hsbm2), 100)
+  expect_true(is_simple(g_hsbm2))
+
+  withr::local_seed(42)
+
+  g_hsbm3 <- sample_hierarchical_sbm(100, 10, rho = c(3, 3, 4) / 10, C = C, p = 1e-15)
+  expect_equal(ecount(g_hsbm3), ecount(g_hsbm1))
+  expect_equal(vcount(g_hsbm3), 100)
+  expect_true(is_simple(g_hsbm3))
+
+  withr::local_seed(42)
+
+  g_hsbm4 <- sample_hierarchical_sbm(100, 10, rho = c(3, 3, 4) / 10, C = C, p = 1 - 1e-15)
+  expect_equal(ecount(g_hsbm4), ecount(g_hsbm2))
+  expect_equal(vcount(g_hsbm4), 100)
+  expect_true(is_simple(g_hsbm4))
+})
+
+test_that("HSBM with 1 cluster per block works", {
+  res <- Matrix::Matrix(0, nrow = 10, ncol = 10, doDiag = FALSE)
+  res[6:10, 1:5] <- res[1:5, 6:10] <- 1
+  g_hsbm <- sample_hierarchical_sbm(10, 5, rho = 1, C = matrix(0), p = 1)
+  expect_equal(g_hsbm[], res)
+})
+
+test_that("HSBM with list arguments works", {
+  blocks <- 5
+  C <- matrix(c(
+    1, 1 / 2, 0,
+    1 / 2, 0, 1 / 2,
+    0, 1 / 2, 1 / 2
+  ), nrow = 3)
+  vertices_per_block <- 10
+  rho <- c(3, 3, 4) / 10
+
+  withr::local_seed(42)
+  g_hsbm1 <- sample_hierarchical_sbm(
+    blocks * vertices_per_block,
+    vertices_per_block,
+    rho = rho, C = C, p = 0
+  )
+
+  withr::local_seed(42)
+  g_hsbm2 <- sample_hierarchical_sbm(
+    blocks * vertices_per_block,
+    rep(vertices_per_block, blocks),
+    rho = rho, C = C, p = 0
+  )
+  expect_equal(g_hsbm1[], g_hsbm2[])
+
+  withr::local_seed(42)
+  g_hsbm3 <- sample_hierarchical_sbm(
+    blocks * vertices_per_block,
+    vertices_per_block,
+    rho = replicate(blocks, rho, simplify = FALSE), C = C, p = 0
+  )
+  expect_equal(g_hsbm1[], g_hsbm3[])
+
+  withr::local_seed(42)
+  g_hsbm4 <- sample_hierarchical_sbm(
+    blocks * vertices_per_block,
+    vertices_per_block,
+    rho = rho, C = replicate(blocks, C, simplify = FALSE), p = 0
+  )
+
+  expect_equal(g_hsbm1[], g_hsbm4[])
+
+  expect_error(
+    sample_hierarchical_sbm(
+      blocks * vertices_per_block,
+      rep(vertices_per_block, blocks),
+      rho = list(rho, rho), C = C, p = 0
+    )
+  )
+
+  ###
+
+  n <- function(x) x / sum(x)
+
+  rho1 <- n(c(1, 2))
+  C1 <- matrix(0, nrow = 2, ncol = 2)
+  rho2 <- n(c(3, 3, 4))
+  C2 <- matrix(0, nrow = 3, ncol = 3)
+  rho3 <- 1
+  C3 <- matrix(0)
+  rho4 <- n(c(2, 1))
+  C4 <- matrix(0, nrow = 2, ncol = 2)
+
+  g_hsbm5 <- sample_hierarchical_sbm(21,
+    m = c(3, 10, 5, 3), rho = list(rho1, rho2, rho3, rho4),
+    C = list(C1, C2, C3, C4), p = 1
+  )
+  expect_true(is_simple(g_hsbm5))
+
+  withr::local_seed(42)
+  g_hsbm6 <- sample_hierarchical_sbm(21,
+    m = c(3, 10, 5, 3), rho = list(rho1, rho2, rho3, rho4),
+    C = list(C1, C2, C3, C4), p = 1 - 1e-10
+  )
+  expect_equal(g_hsbm5[], g_hsbm6[])
+
+  rho1 <- n(c(1, 2))
+  C1 <- matrix(1, nrow = 2, ncol = 2)
+  rho2 <- n(c(3, 3, 4))
+  C2 <- matrix(1, nrow = 3, ncol = 3)
+  rho3 <- 1
+  C3 <- matrix(1)
+  rho4 <- n(c(2, 1))
+  C4 <- matrix(1, nrow = 2, ncol = 2)
+  g_hsbm7 <- sample_hierarchical_sbm(21,
+    m = c(3, 10, 5, 3), rho = list(rho1, rho2, rho3, rho4),
+    C = list(C1, C2, C3, C4), p = 0
+  )
+  expect_true(is_simple(g_hsbm7))
+
+  g_hsbm8 <- sample_hierarchical_sbm(21,
+    m = c(3, 10, 5, 3), rho = list(rho1, rho2, rho3, rho4),
+    C = list(C1, C2, C3, C4), p = 1
+  )
+  expect_equal(g_hsbm5[] + g_hsbm7[], g_hsbm8[])
+})
diff --git a/tests/testthat/test-hsbm.R b/tests/testthat/test-hsbm.R
deleted file mode 100644
index 6eb9732a9e4..00000000000
--- a/tests/testthat/test-hsbm.R
+++ /dev/null
@@ -1,119 +0,0 @@
-test_that("HSBM works", {
-  withr::local_seed(42)
-
-  C <- matrix(c(
-    1, 1 / 2, 0,
-    1 / 2, 0, 1 / 2,
-    0, 1 / 2, 1 / 2
-  ), nrow = 3)
-
-  g <- sample_hierarchical_sbm(100, 10, rho = c(3, 3, 4) / 10, C = C, p = 0)
-  expect_equal(ecount(g), 172)
-  expect_equal(vcount(g), 100)
-  expect_false(is_directed(g))
-
-  withr::local_seed(42)
-
-  g2 <- sample_hierarchical_sbm(100, 10, rho = c(3, 3, 4) / 10, C = C, p = 1)
-  expect_equal(ecount(g2), ecount(g) + 10 * 9 * (90 + 10) / 2)
-  expect_equal(vcount(g2), 100)
-  expect_true(is_simple(g2))
-
-  withr::local_seed(42)
-
-  g3 <- sample_hierarchical_sbm(100, 10, rho = c(3, 3, 4) / 10, C = C, p = 1e-15)
-  expect_equal(ecount(g3), ecount(g))
-  expect_equal(vcount(g3), 100)
-  expect_true(is_simple(g3))
-
-  withr::local_seed(42)
-
-  g4 <- sample_hierarchical_sbm(100, 10, rho = c(3, 3, 4) / 10, C = C, p = 1 - 1e-15)
-  expect_equal(ecount(g4), ecount(g2))
-  expect_equal(vcount(g4), 100)
-  expect_true(is_simple(g4))
-})
-
-test_that("HSBM with 1 cluster per block works", {
-  res <- Matrix::Matrix(0, nrow = 10, ncol = 10, doDiag = FALSE)
-  res[6:10, 1:5] <- res[1:5, 6:10] <- 1
-  g <- sample_hierarchical_sbm(10, 5, rho = 1, C = matrix(0), p = 1)
-  expect_equal(g[], res)
-})
-
-test_that("HSBM with list arguments works", {
-  b <- 5
-  C <- matrix(c(
-    1, 1 / 2, 0,
-    1 / 2, 0, 1 / 2,
-    0, 1 / 2, 1 / 2
-  ), nrow = 3)
-  m <- 10
-  rho <- c(3, 3, 4) / 10
-
-  withr::local_seed(42)
-  g <- sample_hierarchical_sbm(b * m, m, rho = rho, C = C, p = 0)
-
-  withr::local_seed(42)
-  g2 <- sample_hierarchical_sbm(b * m, rep(m, b), rho = rho, C = C, p = 0)
-  expect_equal(g[], g2[])
-
-  withr::local_seed(42)
-  g3 <- sample_hierarchical_sbm(b * m, m, rho = replicate(b, rho, simplify = FALSE), C = C, p = 0)
-  expect_equal(g[], g3[])
-
-  withr::local_seed(42)
-  g4 <- sample_hierarchical_sbm(b * m, m, rho = rho, C = replicate(b, C, simplify = FALSE), p = 0)
-  expect_equal(g[], g4[])
-
-  expect_error(
-    sample_hierarchical_sbm(b * m, rep(m, b), rho = list(rho, rho), C = C, p = 0),
-    "Lengths of `m', `rho' and `C' must match"
-  )
-
-  ###
-
-  n <- function(x) x / sum(x)
-
-  rho1 <- n(c(1, 2))
-  C1 <- matrix(0, nrow = 2, ncol = 2)
-  rho2 <- n(c(3, 3, 4))
-  C2 <- matrix(0, nrow = 3, ncol = 3)
-  rho3 <- 1
-  C3 <- matrix(0)
-  rho4 <- n(c(2, 1))
-  C4 <- matrix(0, nrow = 2, ncol = 2)
-
-  gg1 <- sample_hierarchical_sbm(21,
-    m = c(3, 10, 5, 3), rho = list(rho1, rho2, rho3, rho4),
-    C = list(C1, C2, C3, C4), p = 1
-  )
-  expect_true(is_simple(gg1))
-
-  withr::local_seed(42)
-  gg11 <- sample_hierarchical_sbm(21,
-    m = c(3, 10, 5, 3), rho = list(rho1, rho2, rho3, rho4),
-    C = list(C1, C2, C3, C4), p = 1 - 1e-10
-  )
-  expect_equal(gg1[], gg11[])
-
-  rho1 <- n(c(1, 2))
-  C1 <- matrix(1, nrow = 2, ncol = 2)
-  rho2 <- n(c(3, 3, 4))
-  C2 <- matrix(1, nrow = 3, ncol = 3)
-  rho3 <- 1
-  C3 <- matrix(1)
-  rho4 <- n(c(2, 1))
-  C4 <- matrix(1, nrow = 2, ncol = 2)
-  gg2 <- sample_hierarchical_sbm(21,
-    m = c(3, 10, 5, 3), rho = list(rho1, rho2, rho3, rho4),
-    C = list(C1, C2, C3, C4), p = 0
-  )
-  expect_true(is_simple(gg2))
-
-  gg22 <- sample_hierarchical_sbm(21,
-    m = c(3, 10, 5, 3), rho = list(rho1, rho2, rho3, rho4),
-    C = list(C1, C2, C3, C4), p = 1
-  )
-  expect_equal(gg1[] + gg2[], gg22[])
-})
diff --git a/tests/testthat/test-sbm.game.R b/tests/testthat/test-sbm.game.R
deleted file mode 100644
index b2b43bbb544..00000000000
--- a/tests/testthat/test-sbm.game.R
+++ /dev/null
@@ -1,30 +0,0 @@
-test_that("Generating stochastic block models works", {
-  pm <- matrix(1, nrow = 2, ncol = 2)
-  bs <- c(4, 6)
-  g1 <- sample_sbm(10,
-    pref.matrix = pm, block.sizes = bs,
-    directed = FALSE, loops = FALSE
-  )
-  expect_isomorphic(g1, make_full_graph(10, directed = FALSE, loops = FALSE))
-
-  g2 <- sample_sbm(10,
-    pref.matrix = pm, block.sizes = bs,
-    directed = FALSE, loops = TRUE
-  )
-  g2x <- make_full_graph(10, directed = FALSE, loops = TRUE)
-  expect_equal(g2[sparse = FALSE], g2x[sparse = FALSE])
-
-  g3 <- sample_sbm(10,
-    pref.matrix = pm, block.sizes = bs,
-    directed = TRUE, loops = FALSE
-  )
-  g3x <- make_full_graph(10, directed = TRUE, loops = FALSE)
-  expect_equal(g3[sparse = FALSE], g3x[sparse = FALSE])
-
-  g4 <- sample_sbm(10,
-    pref.matrix = pm, block.sizes = bs,
-    directed = TRUE, loops = TRUE
-  )
-  g4x <- make_full_graph(10, directed = TRUE, loops = TRUE)
-  expect_equal(g4[sparse = FALSE], g4x[sparse = FALSE])
-})
diff --git a/tests/testthat/test-watts.strogatz.game.R b/tests/testthat/test-watts.strogatz.game.R
deleted file mode 100644
index c8dcce7008c..00000000000
--- a/tests/testthat/test-watts.strogatz.game.R
+++ /dev/null
@@ -1,9 +0,0 @@
-test_that("sample_smallworld works", {
-  for (i in 1:50) {
-    p <- runif(1)
-    d <- sample(1:3, 1)
-    nei <- sample(2:5, 1)
-    g <- sample_smallworld(d, 10, nei, p, loops = FALSE)
-    expect_false(any(which_loop(g)))
-  }
-})