diff --git a/NAMESPACE b/NAMESPACE
index d1c40956c1..62529efb22 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -186,6 +186,8 @@ export(bipartite.mapping)
 export(bipartite.projection)
 export(bipartite.projection.size)
 export(bipartite.random.game)
+export(bipartite_gnm)
+export(bipartite_gnp)
 export(bipartite_graph)
 export(bipartite_mapping)
 export(bipartite_projection)
@@ -746,6 +748,8 @@ export(running_mean)
 export(sample_)
 export(sample_asym_pref)
 export(sample_bipartite)
+export(sample_bipartite_gnm)
+export(sample_bipartite_gnp)
 export(sample_chung_lu)
 export(sample_cit_cit_types)
 export(sample_cit_types)
diff --git a/R/games.R b/R/games.R
index 27ddea4eb8..f62e3d5dd9 100644
--- a/R/games.R
+++ b/R/games.R
@@ -265,7 +265,10 @@ callaway.traits.game <- function(nodes, types, edge.per.step = 1, type.dist = re
 #' @keywords internal
 #' @export
 bipartite.random.game <- function(n1, n2, type = c("gnp", "gnm"), p, m, directed = FALSE, mode = c("out", "in", "all")) { # nocov start
-  lifecycle::deprecate_soft("2.0.0", "bipartite.random.game()", "sample_bipartite()")
+  lifecycle::deprecate_warn(
+    "2.0.0", "bipartite.random.game()",
+    details = "Use sample_bipartite_gnp() or sample_bipartite_gnm()"
+  )
   sample_bipartite(n1 = n1, n2 = n2, type = type, p = p, m = m, directed = directed, mode = mode)
 } # nocov end
 
@@ -1695,13 +1698,8 @@ cit_cit_types <- function(...) constructor_spec(sample_cit_cit_types, ...)
 
 #' Bipartite random graphs
 #'
-#' Generate bipartite graphs using the Erdős-Rényi model
-#'
-#' Similarly to unipartite (one-mode) networks, we can define the \eqn{G(n,p)}, and
-#' \eqn{G(n,m)} graph classes for bipartite graphs, via their generating process.
-#' In \eqn{G(n,p)} every possible edge between top and bottom vertices is realized
-#' with probability \eqn{p}, independently of the rest of the edges. In \eqn{G(n,m)}, we
-#' uniformly choose \eqn{m} edges to realize.
+#' `r lifecycle::badge("deprecated")` Generate bipartite graphs using the Erdős-Rényi model.
+#' Use [`sample_bipartite_gnm()`] and [`sample_bipartite_gnp()`] instead.
 #'
 #' @param n1 Integer scalar, the number of bottom vertices.
 #' @param n2 Integer scalar, the number of top vertices.
@@ -1725,72 +1723,161 @@ cit_cit_types <- function(...) constructor_spec(sample_cit_cit_types, ...)
 #' @family games
 #' @export
 #' @keywords graphs
+sample_bipartite <- function(n1, n2, type = c("gnp", "gnm"), p, m,
+                             directed = FALSE, mode = c("out", "in", "all")) {
+
+  type <- igraph.match.arg(type)
+
+  if (type == "gnp") {
+    lifecycle::deprecate_soft(
+      "2.1.3",
+      "sample_bipartite()",
+      "sample_bipartite_gnp()"
+    )
+    sample_bipartite_gnp(n1, n2, p, directed = directed, mode = mode)
+  } else if (type == "gnm") {
+    lifecycle::deprecate_soft(
+      "2.1.3",
+      "sample_bipartite()",
+      "sample_bipartite_gnm()"
+    )
+    sample_bipartite_gnm(n1, n2, m, directed = directed, mode = mode)
+  }
+}
+
+#' @rdname sample_bipartite
+#' @param ... Passed to `sample_bipartite()`.
+#' @export
+bipartite <- function(...) {
+  if (type == "gnp") {
+    lifecycle::deprecate_soft(
+      "2.1.3",
+      "bipartite()",
+      "bipartite_gnp()"
+    )
+    bipartite_gnp(...)
+  } else if (type == "gnm") {
+    lifecycle::deprecate_soft(
+      "2.1.3",
+      "bipartite()",
+      "bipartite_gnm()"
+    )
+    bipartite_gnm(...)
+  }
+
+}
+
+#' @rdname sample_bipartite_gnm
+#' @param ... Passed to `sample_bipartite_gnm()`.
+#' @export
+bipartite_gnm <- function(...) constructor_spec(sample_bipartite_gnm, ...)
+
+#' @rdname sample_bipartite_gnm
+#' @param ... Passed to `sample_bipartite_gnp()`.
+#' @export
+bipartite_gnp <- function(...) constructor_spec(sample_bipartite_gnp, ...)
+
+#' Bipartite random graphs
+#'
+#' Generate bipartite graphs using the Erdős-Rényi model
+#'
+#' Similarly to unipartite (one-mode) networks, we can define the \eqn{G(n,p)}, and
+#' \eqn{G(n,m)} graph classes for bipartite graphs, via their generating process.
+#' In \eqn{G(n,p)} every possible edge between top and bottom vertices is realized
+#' with probability \eqn{p}, independently of the rest of the edges. In \eqn{G(n,m)}, we
+#' uniformly choose \eqn{m} edges to realize.
+#'
+#'
+#' @param n1 Integer scalar, the number of bottom vertices.
+#' @param n2 Integer scalar, the number of top vertices.
+#' @param p Real scalar, connection probability for \eqn{G(n,p)} graphs.
+#' @param m Integer scalar, the number of edges for \eqn{G(n,m)} graphs.
+#' @param directed Logical scalar, whether to create a directed graph. See also
+#'   the `mode` argument.
+#' @param mode Character scalar, specifies how to direct the edges in directed
+#'   graphs. If it is \sQuote{out}, then directed edges point from bottom
+#'   vertices to top vertices. If it is \sQuote{in}, edges point from top
+#'   vertices to bottom vertices. \sQuote{out} and \sQuote{in} do not generate
+#'   mutual edges. If this argument is \sQuote{all}, then each edge direction is
+#'   considered independently and mutual edges might be generated. This argument
+#'   is ignored for undirected graphs.
+#' @inheritParams rlang::args_dots_empty
 #' @examples
 #'
 #' ## empty graph
-#' sample_bipartite(10, 5, p = 0)
+#' sample_bipartite_gnp(10, 5, p = 0)
 #'
 #' ## full graph
-#' sample_bipartite(10, 5, p = 1)
+#' sample_bipartite_gnp(10, 5, p = 1)
 #'
 #' ## random bipartite graph
-#' sample_bipartite(10, 5, p = .1)
+#' sample_bipartite_gnp(10, 5, p = .1)
 #'
 #' ## directed bipartite graph, G(n,m)
-#' sample_bipartite(10, 5, type = "Gnm", m = 20, directed = TRUE, mode = "all")
+#' sample_bipartite_gnm(10, 5, m = 20, directed = TRUE, mode = "all")
 #'
-sample_bipartite <- function(n1, n2, type = c("gnp", "gnm"), p, m,
-                             directed = FALSE, mode = c("out", "in", "all")) {
+#' @family games
+#' @export
+sample_bipartite_gnm <- function(n1, n2, m,
+                                ...,
+                                directed = FALSE,
+                                mode = c("out", "in", "all")) {
+  check_dots_empty()
+
   n1 <- as.numeric(n1)
   n2 <- as.numeric(n2)
-  type <- igraph.match.arg(type)
-  if (!missing(p)) {
-    p <- as.numeric(p)
-  }
-  if (!missing(m)) {
-    m <- as.numeric(m)
-  }
+
+  m <- as.numeric(m)
+
   directed <- as.logical(directed)
+
   mode <- switch(igraph.match.arg(mode),
     "out" = 1,
     "in" = 2,
     "all" = 3
   )
 
-  if (type == "gnp" && missing(p)) {
-    stop("Connection probability `p' is not 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")
-  }
-  if (type == "gnm" && !missing(p)) {
-    cli::cli_warn("Connection probability {.arg p} is ignored for Gnp graph.")
-  }
-
   on.exit(.Call(R_igraph_finalizer))
-  if (type == "gnp") {
-    res <- .Call(R_igraph_bipartite_game_gnp, n1, n2, p, directed, mode)
-    res <- set_vertex_attr(res$graph, "type", value = res$types)
-    res$name <- "Bipartite Gnp random graph"
-    res$p <- p
-  } else if (type == "gnm") {
-    res <- .Call(R_igraph_bipartite_game_gnm, n1, n2, m, directed, mode)
-    res <- set_vertex_attr(res$graph, "type", value = res$types)
-    res$name <- "Bipartite Gnm random graph"
-    res$m <- m
-  }
+
+  res <- .Call(R_igraph_bipartite_game_gnm, n1, n2, m, directed, mode)
+  res <- set_vertex_attr(res$graph, "type", value = res$types)
+  res$name <- "Bipartite Gnm random graph"
+  res$m <- m
 
   res
-}
 
-#' @rdname sample_bipartite
-#' @param ... Passed to `sample_bipartite()`.
+}
+#' @rdname sample_bipartite_gnm
 #' @export
-bipartite <- function(...) constructor_spec(sample_bipartite, ...)
+sample_bipartite_gnp <- function(n1, n2, p,
+                                ...,
+                                directed = FALSE,
+                                mode = c("out", "in", "all")) {
+  check_dots_empty()
+
+  n1 <- as.numeric(n1)
+  n2 <- as.numeric(n2)
+
+  p <- as.numeric(p)
+
+  directed <- as.logical(directed)
 
+  mode <- switch(igraph.match.arg(mode),
+    "out" = 1,
+    "in" = 2,
+    "all" = 3
+  )
+
+  on.exit(.Call(R_igraph_finalizer))
+
+  res <- .Call(R_igraph_bipartite_game_gnp, n1, n2, p, directed, mode)
+  res <- set_vertex_attr(res$graph, "type", value = res$types)
+  res$name <- "Bipartite Gnp random graph"
+  res$p <- p
+
+  res
+
+}
 
 #' Sample stochastic block model
 #'
diff --git a/man/closeness.Rd b/man/closeness.Rd
index 0502a7f28f..dc860d003e 100644
--- a/man/closeness.Rd
+++ b/man/closeness.Rd
@@ -43,7 +43,7 @@ Numeric vector with the closeness values of all the vertices in
 \code{v}.
 }
 \description{
-Closeness centrality measures how many steps is required to access every other
+Closeness centrality measures how many steps are required to access every other
 vertex from a given vertex.
 }
 \details{
diff --git a/man/erdos.renyi.game.Rd b/man/erdos.renyi.game.Rd
index 3390cd3c4d..caa0e0f924 100644
--- a/man/erdos.renyi.game.Rd
+++ b/man/erdos.renyi.game.Rd
@@ -52,6 +52,7 @@ Mathematicae} 6, 290--297 (1959).
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
 \code{\link{sample_chung_lu}()},
diff --git a/man/sample_.Rd b/man/sample_.Rd
index ebea0a7d2d..212a9fe71b 100644
--- a/man/sample_.Rd
+++ b/man/sample_.Rd
@@ -46,6 +46,7 @@ blocky3 <- pref_matrix \%>\%
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_bipartite}()},
 \code{\link{sample_chung_lu}()},
diff --git a/man/sample_bipartite.Rd b/man/sample_bipartite.Rd
index 313bc7b16a..767047d151 100644
--- a/man/sample_bipartite.Rd
+++ b/man/sample_bipartite.Rd
@@ -48,32 +48,12 @@ is ignored for undirected graphs.}
 A bipartite igraph graph.
 }
 \description{
-Generate bipartite graphs using the Erdős-Rényi model
-}
-\details{
-Similarly to unipartite (one-mode) networks, we can define the \eqn{G(n,p)}, and
-\eqn{G(n,m)} graph classes for bipartite graphs, via their generating process.
-In \eqn{G(n,p)} every possible edge between top and bottom vertices is realized
-with probability \eqn{p}, independently of the rest of the edges. In \eqn{G(n,m)}, we
-uniformly choose \eqn{m} edges to realize.
-}
-\examples{
-
-## empty graph
-sample_bipartite(10, 5, p = 0)
-
-## full graph
-sample_bipartite(10, 5, p = 1)
-
-## random bipartite graph
-sample_bipartite(10, 5, p = .1)
-
-## directed bipartite graph, G(n,m)
-sample_bipartite(10, 5, type = "Gnm", m = 20, directed = TRUE, mode = "all")
-
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} Generate bipartite graphs using the Erdős-Rényi model.
+Use \code{\link[=sample_bipartite_gnm]{sample_bipartite_gnm()}} and \code{\link[=sample_bipartite_gnp]{sample_bipartite_gnp()}} instead.
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_chung_lu}()},
diff --git a/man/sample_bipartite_gnm.Rd b/man/sample_bipartite_gnm.Rd
new file mode 100644
index 0000000000..a6ec63a4d2
--- /dev/null
+++ b/man/sample_bipartite_gnm.Rd
@@ -0,0 +1,108 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/games.R
+\name{bipartite_gnm}
+\alias{bipartite_gnm}
+\alias{bipartite_gnp}
+\alias{sample_bipartite_gnm}
+\alias{sample_bipartite_gnp}
+\title{Bipartite random graphs}
+\usage{
+bipartite_gnm(...)
+
+bipartite_gnp(...)
+
+sample_bipartite_gnm(
+  n1,
+  n2,
+  m,
+  ...,
+  directed = FALSE,
+  mode = c("out", "in", "all")
+)
+
+sample_bipartite_gnp(
+  n1,
+  n2,
+  p,
+  ...,
+  directed = FALSE,
+  mode = c("out", "in", "all")
+)
+}
+\arguments{
+\item{...}{Passed to \code{sample_bipartite_gnp()}.}
+
+\item{n1}{Integer scalar, the number of bottom vertices.}
+
+\item{n2}{Integer scalar, the number of top vertices.}
+
+\item{m}{Integer scalar, the number of edges for \eqn{G(n,m)} graphs.}
+
+\item{directed}{Logical scalar, whether to create a directed graph. See also
+the \code{mode} argument.}
+
+\item{mode}{Character scalar, specifies how to direct the edges in directed
+graphs. If it is \sQuote{out}, then directed edges point from bottom
+vertices to top vertices. If it is \sQuote{in}, edges point from top
+vertices to bottom vertices. \sQuote{out} and \sQuote{in} do not generate
+mutual edges. If this argument is \sQuote{all}, then each edge direction is
+considered independently and mutual edges might be generated. This argument
+is ignored for undirected graphs.}
+
+\item{p}{Real scalar, connection probability for \eqn{G(n,p)} graphs.}
+}
+\description{
+Generate bipartite graphs using the Erdős-Rényi model
+}
+\details{
+Similarly to unipartite (one-mode) networks, we can define the \eqn{G(n,p)}, and
+\eqn{G(n,m)} graph classes for bipartite graphs, via their generating process.
+In \eqn{G(n,p)} every possible edge between top and bottom vertices is realized
+with probability \eqn{p}, independently of the rest of the edges. In \eqn{G(n,m)}, we
+uniformly choose \eqn{m} edges to realize.
+}
+\examples{
+
+## empty graph
+sample_bipartite_gnp(10, 5, p = 0)
+
+## full graph
+sample_bipartite_gnp(10, 5, p = 1)
+
+## random bipartite graph
+sample_bipartite_gnp(10, 5, p = .1)
+
+## directed bipartite graph, G(n,m)
+sample_bipartite_gnm(10, 5, m = 20, directed = TRUE, mode = "all")
+
+}
+\seealso{
+Random graph models (games)
+\code{\link{erdos.renyi.game}()},
+\code{\link{sample_}()},
+\code{\link{sample_bipartite}()},
+\code{\link{sample_chung_lu}()},
+\code{\link{sample_correlated_gnp}()},
+\code{\link{sample_correlated_gnp_pair}()},
+\code{\link{sample_degseq}()},
+\code{\link{sample_dot_product}()},
+\code{\link{sample_fitness}()},
+\code{\link{sample_fitness_pl}()},
+\code{\link{sample_forestfire}()},
+\code{\link{sample_gnm}()},
+\code{\link{sample_gnp}()},
+\code{\link{sample_grg}()},
+\code{\link{sample_growing}()},
+\code{\link{sample_hierarchical_sbm}()},
+\code{\link{sample_islands}()},
+\code{\link{sample_k_regular}()},
+\code{\link{sample_last_cit}()},
+\code{\link{sample_pa}()},
+\code{\link{sample_pa_age}()},
+\code{\link{sample_pref}()},
+\code{\link{sample_sbm}()},
+\code{\link{sample_smallworld}()},
+\code{\link{sample_traits_callaway}()},
+\code{\link{sample_tree}()}
+}
+\concept{games}
diff --git a/man/sample_chung_lu.Rd b/man/sample_chung_lu.Rd
index 529bd7fa22..d636643229 100644
--- a/man/sample_chung_lu.Rd
+++ b/man/sample_chung_lu.Rd
@@ -168,6 +168,7 @@ with sharply specified degrees. \code{\link[=sample_gnp]{sample_gnp()}} creates
 fixed connection probability \eqn{p} between all vertex pairs.
 
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_correlated_gnp.Rd b/man/sample_correlated_gnp.Rd
index c180df2459..563730a349 100644
--- a/man/sample_correlated_gnp.Rd
+++ b/man/sample_correlated_gnp.Rd
@@ -56,6 +56,7 @@ graph matching for correlated Erdős-Rényi graphs.
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_correlated_gnp_pair.Rd b/man/sample_correlated_gnp_pair.Rd
index 473733336a..f6ba5cfd95 100644
--- a/man/sample_correlated_gnp_pair.Rd
+++ b/man/sample_correlated_gnp_pair.Rd
@@ -49,6 +49,7 @@ graph matching for correlated Erdős-Rényi graphs.
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_degseq.Rd b/man/sample_degseq.Rd
index 8000211b9d..0f381f59d6 100644
--- a/man/sample_degseq.Rd
+++ b/man/sample_degseq.Rd
@@ -195,6 +195,7 @@ all(degree(powerlaw_vl_graph) == powerlaw_degrees)
 \code{\link[=realize_degseq]{realize_degseq()}} for a deterministic variant.
 
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_dot_product.Rd b/man/sample_dot_product.Rd
index e84b56a390..f658f4e08b 100644
--- a/man/sample_dot_product.Rd
+++ b/man/sample_dot_product.Rd
@@ -57,6 +57,7 @@ for social networks. Dissertation, Johns Hopkins University, Maryland, USA,
 and \code{\link[=sample_sphere_volume]{sample_sphere_volume()}} for sampling position vectors.
 
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_fitness.Rd b/man/sample_fitness.Rd
index f6ab1840d9..6b8c34bf65 100644
--- a/man/sample_fitness.Rd
+++ b/man/sample_fitness.Rd
@@ -79,6 +79,7 @@ distribution in scale-free networks. \emph{Phys Rev Lett} 87(27):278701,
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_fitness_pl.Rd b/man/sample_fitness_pl.Rd
index 3e092097d3..18e547e3f7 100644
--- a/man/sample_fitness_pl.Rd
+++ b/man/sample_fitness_pl.Rd
@@ -88,6 +88,7 @@ scale-free networks under the Achlioptas process. \emph{Phys Rev Lett}
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_forestfire.Rd b/man/sample_forestfire.Rd
index 6d44e99b1e..1d6c6f9765 100644
--- a/man/sample_forestfire.Rd
+++ b/man/sample_forestfire.Rd
@@ -77,6 +77,7 @@ conference on Knowledge discovery in data mining}, 177--187, 2005.
 model.
 
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_gnm.Rd b/man/sample_gnm.Rd
index 5f71695a37..aaf6353a42 100644
--- a/man/sample_gnm.Rd
+++ b/man/sample_gnm.Rd
@@ -43,6 +43,7 @@ Mathematicae} 6, 290--297 (1959).
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_gnp.Rd b/man/sample_gnp.Rd
index 698f159aa8..77d1347526 100644
--- a/man/sample_gnp.Rd
+++ b/man/sample_gnp.Rd
@@ -54,6 +54,7 @@ Mathematicae} 6, 290--297 (1959).
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_grg.Rd b/man/sample_grg.Rd
index e891532cc9..335d91b29f 100644
--- a/man/sample_grg.Rd
+++ b/man/sample_grg.Rd
@@ -45,6 +45,7 @@ g2 <- sample_grg(1000, 0.05, torus = TRUE)
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_growing.Rd b/man/sample_growing.Rd
index 7621ede2b1..ab7826d489 100644
--- a/man/sample_growing.Rd
+++ b/man/sample_growing.Rd
@@ -42,6 +42,7 @@ g2 <- sample_growing(500, citation = TRUE)
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_hierarchical_sbm.Rd b/man/sample_hierarchical_sbm.Rd
index fe5aebc0bb..48c0ab53d1 100644
--- a/man/sample_hierarchical_sbm.Rd
+++ b/man/sample_hierarchical_sbm.Rd
@@ -57,6 +57,7 @@ if (require(Matrix)) {
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_islands.Rd b/man/sample_islands.Rd
index 913e15b40f..be031b3f02 100644
--- a/man/sample_islands.Rd
+++ b/man/sample_islands.Rd
@@ -36,6 +36,7 @@ oc
 \code{\link[=sample_gnp]{sample_gnp()}}
 
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_k_regular.Rd b/man/sample_k_regular.Rd
index 66f2df1876..8258329274 100644
--- a/man/sample_k_regular.Rd
+++ b/man/sample_k_regular.Rd
@@ -48,6 +48,7 @@ sapply(k10, plot, vertex.label = NA)
 sequence.
 
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_last_cit.Rd b/man/sample_last_cit.Rd
index ac9b9dc09d..a5a55001a9 100644
--- a/man/sample_last_cit.Rd
+++ b/man/sample_last_cit.Rd
@@ -80,6 +80,7 @@ vertex only.
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_pa.Rd b/man/sample_pa.Rd
index d50f175ff6..ec8c4f8ada 100644
--- a/man/sample_pa.Rd
+++ b/man/sample_pa.Rd
@@ -127,6 +127,7 @@ de Solla Price, D. J. 1965. Networks of Scientific Papers \emph{Science},
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_pa_age.Rd b/man/sample_pa_age.Rd
index dc7a0b45ae..598cb93af4 100644
--- a/man/sample_pa_age.Rd
+++ b/man/sample_pa_age.Rd
@@ -136,6 +136,7 @@ max(degree(g3))
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_pref.Rd b/man/sample_pref.Rd
index 2b1140ec00..43aadaa2fa 100644
--- a/man/sample_pref.Rd
+++ b/man/sample_pref.Rd
@@ -99,6 +99,7 @@ tkplot(g, layout = layout_in_circle)
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_sbm.Rd b/man/sample_sbm.Rd
index 9985a7f18a..c771995ca8 100644
--- a/man/sample_sbm.Rd
+++ b/man/sample_sbm.Rd
@@ -53,6 +53,7 @@ and evaluation. \emph{Social Networks}, 14, 5--61.
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_smallworld.Rd b/man/sample_smallworld.Rd
index 8bd2e7a815..28f9adc0b4 100644
--- a/man/sample_smallworld.Rd
+++ b/man/sample_smallworld.Rd
@@ -65,6 +65,7 @@ Duncan J Watts and Steven H Strogatz: Collective dynamics of
 \code{\link[=make_lattice]{make_lattice()}}, \code{\link[=rewire]{rewire()}}
 
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_traits_callaway.Rd b/man/sample_traits_callaway.Rd
index 4a5e65ffa0..7d7886783e 100644
--- a/man/sample_traits_callaway.Rd
+++ b/man/sample_traits_callaway.Rd
@@ -79,6 +79,7 @@ g2 <- sample_traits(1000, 2, k = 2, pref.matrix = matrix(c(1, 0, 0, 1), ncol = 2
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/man/sample_tree.Rd b/man/sample_tree.Rd
index 046ca29e35..a5c600ba8f 100644
--- a/man/sample_tree.Rd
+++ b/man/sample_tree.Rd
@@ -37,6 +37,7 @@ g <- sample_tree(100, method = "lerw")
 }
 \seealso{
 Random graph models (games)
+\code{\link{bipartite_gnm}()},
 \code{\link{erdos.renyi.game}()},
 \code{\link{sample_}()},
 \code{\link{sample_bipartite}()},
diff --git a/tests/testthat/_snaps/games.md b/tests/testthat/_snaps/games.md
index 922ef9a505..90167ce0cc 100644
--- a/tests/testthat/_snaps/games.md
+++ b/tests/testthat/_snaps/games.md
@@ -14,3 +14,21 @@
       Error in `sample_degseq()`:
       ! At vendor/cigraph/src/games/degree_sequence_vl/gengraph_mr-connected.cpp:<linenumber> : Cannot realize the given degree sequence as an undirected, simple graph. Invalid value
 
+# sample_bipartite() deprecation
+
+    Code
+      s <- sample_bipartite(10, 5, type = "gnp", p = 0)
+    Condition
+      Warning:
+      `sample_bipartite()` was deprecated in igraph 2.1.3.
+      i Please use `sample_bipartite_gnp()` instead.
+
+---
+
+    Code
+      s <- sample_bipartite(10, 5, type = "gnm", m = 0)
+    Condition
+      Warning:
+      `sample_bipartite()` was deprecated in igraph 2.1.3.
+      i Please use `sample_bipartite_gnm()` instead.
+
diff --git a/tests/testthat/test-ba.game.R b/tests/testthat/test-ba.game.R
deleted file mode 100644
index e66d45a3e0..0000000000
--- 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.projection.R b/tests/testthat/test-bipartite.projection.R
index 24a82b879d..fed5648794 100644
--- a/tests/testthat/test-bipartite.projection.R
+++ b/tests/testthat/test-bipartite.projection.R
@@ -40,7 +40,7 @@ test_that("bipartite_projection works", {
 test_that("bipartite_projection can calculate only one projection", {
   withr::local_seed(42)
 
-  g <- sample_bipartite(5, 10, p = .3)
+  g <- sample_bipartite_gnp(5, 10, p = .3)
   proj <- bipartite_projection(g)
   proj1 <- bipartite_projection(g, which = "false")
   proj2 <- bipartite_projection(g, which = "true")
diff --git a/tests/testthat/test-bipartite.random.game.R b/tests/testthat/test-bipartite.random.game.R
deleted file mode 100644
index 9d2df975e3..0000000000
--- 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 3070de7635..0000000000
--- 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-games.R b/tests/testthat/test-games.R
index 9b0c9991ca..76b7350b14 100644
--- a/tests/testthat/test-games.R
+++ b/tests/testthat/test-games.R
@@ -158,3 +158,330 @@ test_that("sample_degseq works() -- old method names", {
     "must be"
   )
 })
+
+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))
+})
+
+test_that("sample_bipartite() works", {
+  rlang::local_options(lifecycle_verbosity = "quiet")
+  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)
+})
+
+test_that("sample_bipartite() deprecation", {
+  expect_snapshot(s <- sample_bipartite(10, 5, type = "gnp", p = 0))
+  expect_snapshot(s <- sample_bipartite(10, 5, type = "gnm", m = 0))
+})
+
+test_that("sample_bipartite_gnp() works", {
+  withr::local_seed(42)
+  g1 <- sample_bipartite_gnp(10, 5, 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_gnp(10, 5, 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_gnp(10, 5, p = .1, directed = TRUE, mode = "in")
+  expect_output(print_all(g3), "11->3")
+
+  g7 <- sample_bipartite_gnp(10, 5,
+   p = 0.9999, directed = TRUE,
+    mode = "all"
+  )
+  expect_equal(ecount(g7), 100)
+
+})
+test_that("sample_bipartite_gnm() works", {
+  withr::local_seed(42)
+
+  g4 <- sample_bipartite_gnm(10, 5, 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_gnm(10, 5, 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->11 7->11 8->11 8->12 4->13 5->13 6->13 9->13")
+
+  g6 <- sample_bipartite_gnm(10, 5, 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), "11-> 4 11-> 5 12-> 7 12-> 8 12-> 9 14-> 6 14->10 15-> 4")
+
+  g8 <- sample_bipartite_gnm(10, 5,
+    m = 99, directed = TRUE,
+    mode = "all"
+  )
+  expect_equal(ecount(g8), 99)
+})
+
+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[])
+})
+
+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-hsbm.R b/tests/testthat/test-hsbm.R
deleted file mode 100644
index 6eb9732a9e..0000000000
--- 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[])
-})