Move poisson_edges and allow_self_loops logic into model object const…

…ructors (#29)

NAMESPACE

@@ -19,23 +19,14 @@ S3method(print,undirected_overlapping_sbm)
 S3method(print,undirected_planted_partition)
 S3method(print,undirected_sbm)
 S3method(sample_edgelist,Matrix)
-S3method(sample_edgelist,directed_erdos_renyi)
 S3method(sample_edgelist,directed_factor_model)
 S3method(sample_edgelist,matrix)
-S3method(sample_edgelist,undirected_erdos_renyi)
 S3method(sample_edgelist,undirected_factor_model)
-S3method(sample_edgelist,undirected_sbm)
-S3method(sample_igraph,directed_erdos_renyi)
 S3method(sample_igraph,directed_factor_model)
-S3method(sample_igraph,undirected_erdos_renyi)
 S3method(sample_igraph,undirected_factor_model)
-S3method(sample_sparse,directed_erdos_renyi)
 S3method(sample_sparse,directed_factor_model)
-S3method(sample_sparse,undirected_erdos_renyi)
 S3method(sample_sparse,undirected_factor_model)
-S3method(sample_tidygraph,directed_erdos_renyi)
 S3method(sample_tidygraph,directed_factor_model)
-S3method(sample_tidygraph,undirected_erdos_renyi)
 S3method(sample_tidygraph,undirected_factor_model)
 S3method(svds,directed_factor_model)
 S3method(svds,undirected_factor_model)

NEWS.md

@@ -1,5 +1,11 @@
 # fastRG (development version)
 
+## Breaking changes
+
+- Users must now pass `poisson_edges` and `allow_self_loops` arguments to model object constructors (i.e. `sbm()`) rather than `sample_*()` methods. Additionally, when `poisson_edges = FALSE`, the mixing matrix `S` is taken (after degree-scaling and possible symmetrization for undirected models) to represent desired inter-factor connection probabilities, and thus should be between zero and one. This Bernoulli-parameterized `S` is then transformed into the equivalent (or approximately equivalent) Poisson `S`. See Section 2.3 of Rohe et al. (2017) for additional details about this conversion and approximation of Bernoulli graphs by Poisson graphs (#29).
+
+## Other news
+
 * Add overlapping stochastic blockmodel (#7, #25)
 * Add directed degree-corrected stochastic blockmodels (#18)
 * Allow rank 1 undirected stochastic block models

R/directed_dcsbm.R

@@ -213,6 +213,7 @@ validate_directed_dcsbm <- function(x) {
 #'   block. Defaults to `TRUE`.
 #'
 #' @inheritDotParams directed_factor_model expected_in_degree expected_density expected_out_degree
+#' @inheritParams directed_factor_model
 #'
 #' @return A `directed_dcsbm` S3 object, a subclass of the
 #'   [directed_factor_model()] with the following additional
@@ -292,15 +293,15 @@ validate_directed_dcsbm <- function(x) {
 #' ## Edge formulation
 #'
 #' Once we know the block memberships \eqn{x} and \eqn{y}
-#' and the degree  heterogeneity parameters \eqn{\theta_in} and
-#' \eqn{\theta_out}, we need one more
+#' and the degree  heterogeneity parameters \eqn{\theta_{in}} and
+#' \eqn{\theta_{out}}, we need one more
 #' ingredient, which is the baseline intensity of connections
 #' between nodes in block `i` and block `j`. Then each edge forms
 #' independently according to a Poisson distribution with
 #' parameters
 #'
 #' \deqn{
-#'   \lambda = \theta_in * B[x, y] * \theta_out.
+#'   \lambda = \theta_{in} * B_{x, y} * \theta_{out}.
 #' }
 #'
 #' @examples
@@ -330,7 +331,9 @@ directed_dcsbm <- function(
   pi_in = rep(1 / k_in, k_in),
   pi_out = rep(1 / k_out, k_out),
   sort_nodes = TRUE,
-  force_identifiability = TRUE) {
+  force_identifiability = TRUE,
+  poisson_edges = TRUE,
+  allow_self_loops = TRUE) {
 
   ### heterogeneity parameters
 
@@ -500,6 +503,8 @@ directed_dcsbm <- function(
     pi_in = pi_in,
     pi_out = pi_out,
     sorted = sort_nodes,
+    poisson_edges = poisson_edges,
+    allow_self_loops = allow_self_loops,
     ...
   )
 
@@ -532,6 +537,9 @@ print.directed_dcsbm <- function(x, ...) {
   cat("S:", dim_and_class(x$S), "\n")
   cat("Y:", dim_and_class(x$Y), "\n\n")
 
+  cat("Poisson edges:", as.character(x$poisson_edges), "\n")
+  cat("Allow self loops:", as.character(x$allow_self_loops), "\n\n")
+
   cat(glue("Expected edges: {round(expected_edges(x))}\n", .trim = FALSE))
   cat(glue("Expected in degree: {round(expected_in_degree(x), 1)}\n", .trim = FALSE))
   cat(glue("Expected out degree: {round(expected_out_degree(x), 1)}\n", .trim = FALSE))

R/directed_erdos_renyi.R

@@ -1,8 +1,10 @@
-new_directed_erdos_renyi <- function(X, S, Y, p, ...) {
+new_directed_erdos_renyi <- function(X, S, Y, p, poisson_edges, allow_self_loops,...) {
 
   er <- directed_factor_model(
     X, S, Y, ...,
-    subclass = "directed_erdos_renyi"
+    subclass = "directed_erdos_renyi",
+    poisson_edges = poisson_edges,
+    allow_self_loops = allow_self_loops
   )
 
   er$p <- p
@@ -21,8 +23,8 @@ validate_directed_erdos_renyi <- function(x) {
     stop("`Y` must have a single column.", call. = FALSE)
   }
 
-  if (values$p <= 0 || 1 <= values$p) {
-    stop("`p` must be strictly between zero and one.", call. = FALSE)
+  if (values$p < 0) {
+    stop("`p` must be strictly non-negative.", call. = FALSE)
   }
 
   x
@@ -36,12 +38,10 @@ validate_directed_erdos_renyi <- function(x) {
 #'   either `p`, `expected_in_degree`, or `expected_out_degree`.
 #'
 #' @inheritDotParams directed_factor_model expected_in_degree expected_out_degree
-#'
-#' @return Never returns Poisson edges.
+#' @inherit directed_factor_model params return
 #'
 #' @export
 #'
-#' @family bernoulli graphs
 #' @family erdos renyi
 #' @family directed graphs
 #'
@@ -59,7 +59,9 @@ validate_directed_erdos_renyi <- function(x) {
 #' A
 #'
 directed_erdos_renyi <- function(
-  n, ..., p = NULL) {
+  n, ..., p = NULL,
+  poisson_edges = TRUE,
+  allow_self_loops = TRUE) {
 
   X <- Matrix(1, nrow = n, ncol = 1)
   Y <- Matrix(1, nrow = n, ncol = 1)
@@ -77,51 +79,14 @@ directed_erdos_renyi <- function(
     p <- 0.5  # doesn't matter, will get rescaled anyway
   }
 
-  poisson_p <- -log(1 - p)
-  S <- matrix(poisson_p, nrow = 1, ncol = 1)
+  S <- matrix(p, nrow = 1, ncol = 1)
 
-  er <- new_directed_erdos_renyi(X, S, Y, p = p, ...)
-  validate_directed_erdos_renyi(er)
-}
-
-# dispatch hacks to always avoid Poisson edges ---------------------------------
-
-#' @rdname sample_edgelist
-#' @export
-sample_edgelist.directed_erdos_renyi <- function(
-  factor_model,
-  ...,
-  allow_self_loops = TRUE) {
-
-  NextMethod("sample_edgelist", factor_model, ..., poisson_edges = FALSE)
-}
-
-#' @rdname sample_sparse
-#' @export
-sample_sparse.directed_erdos_renyi <- function(
-  factor_model,
-  ...,
-  allow_self_loops = TRUE) {
-
-  NextMethod("sample_sparse", factor_model, ..., poisson_edges = FALSE)
-}
-
-#' @rdname sample_igraph
-#' @export
-sample_igraph.directed_erdos_renyi <- function(
-  factor_model,
-  ...,
-  allow_self_loops = TRUE) {
-
-  NextMethod("sample_igraph", factor_model, ..., poisson_edges = FALSE)
-}
-
-#' @rdname sample_tidygraph
-#' @export
-sample_tidygraph.directed_erdos_renyi <- function(
-  factor_model,
-  ...,
-  allow_self_loops = TRUE) {
+  er <- new_directed_erdos_renyi(
+    X, S, Y, p = p,
+    poisson_edges = poisson_edges,
+    allow_self_loops = allow_self_loops,
+    ...
+  )
 
-  NextMethod("sample_tidygraph", factor_model, ..., poisson_edges = FALSE)
+  validate_directed_erdos_renyi(er)
 }

R/directed_factor_model.R

@@ -1,5 +1,7 @@
 new_directed_factor_model <- function(
   X, S, Y,
+  poisson_edges,
+  allow_self_loops,
   ...,
   subclass = character()) {
 
@@ -17,7 +19,9 @@ new_directed_factor_model <- function(
     n = n,
     k1 = k1,
     d = d,
-    k2 = k2
+    k2 = k2,
+    poisson_edges = poisson_edges,
+    allow_self_loops = allow_self_loops
   )
 
   class(model) <- c(subclass, "directed_factor_model")
@@ -43,6 +47,14 @@ validate_directed_factor_model <- function(x) {
     stop("`k2` must equal the number of columns in `S`", call. = FALSE)
   }
 
+  if (!is.logical(values$poisson_edges)) {
+    stop("`poisson_edges` must be a logical(1) vector.", call. = FALSE)
+  }
+
+  if (!is.logical(values$allow_self_loops)) {
+    stop("`allow_self_loops` must be a logical(1) vector.", call. = FALSE)
+  }
+
   x
 }
 
@@ -86,6 +98,21 @@ validate_directed_factor_model <- function(x) {
 #'   to achieve this. Defaults to `NULL`. Specify only one of
 #'   `expected_in_degree`, `expected_out_degree`, and `expected_density`.
 #'
+#' @param poisson_edges Logical indicating whether or not
+#'   multiple edges are allowed to form between a pair of
+#'   nodes. Defaults to `TRUE`. When `FALSE`, sampling proceeds
+#'   as usual, and duplicate edges are removed afterwards. Further,
+#'   when `FALSE`, we assume that `S` specifies a desired between-factor
+#'   connection probability, and back-transform this `S` to the
+#'   appropriate Poisson intensity parameter to approximate Bernoulli
+#'   factor connection probabilities. See Section 2.3 of Rohe et al. (2017)
+#'   for some additional details.
+#'
+#' @param allow_self_loops Logical indicating whether or not
+#'   nodes should be allowed to form edges with themselves.
+#'   Defaults to `TRUE`. When `FALSE`, sampling proceeds allowing
+#'   self-loops, and these are then removed after the fact.
+#'
 #' @return A `directed_factor_model` S3 class based on a list
 #'   with the following elements:
 #'
@@ -107,6 +134,12 @@ validate_directed_factor_model <- function(x) {
 #'   - `k2`: The dimension of the latent node position vectors
 #'     encoding outgoing latent communities (i.e. in `Y`).
 #'
+#'   - `poisson_edges`: Whether or not the graph is taken to be have
+#'     Poisson or Bernoulli edges, as indicated by a logical vector
+#'     of length 1.
+#'
+#'   - `allow_self_loops`: Whether or not self loops are allowed.
+#'
 #' @export
 #'
 #' @examples
@@ -133,7 +166,9 @@ directed_factor_model <- function(
   ...,
   expected_in_degree = NULL,
   expected_out_degree = NULL,
-  expected_density = NULL) {
+  expected_density = NULL,
+  poisson_edges = TRUE,
+  allow_self_loops = TRUE) {
 
   X <- Matrix(X)
   S <- Matrix(S)
@@ -153,7 +188,12 @@ directed_factor_model <- function(
     )
   }
 
-  fm <- new_directed_factor_model(X, S, Y, ...)
+  fm <- new_directed_factor_model(
+    X, S, Y,
+    poisson_edges = poisson_edges,
+    allow_self_loops = allow_self_loops,
+    ...
+  )
 
   if (!is.null(expected_in_degree)) {
 
@@ -193,6 +233,24 @@ directed_factor_model <- function(
 
   fm$S <- S
 
+  if (!poisson_edges) {
+
+    # when poisson_edges = FALSE, S is the desired Bernoulli edge probability.
+    # we must
+    # back-transform it to a Poisson parameterization of S. see section 2.3
+    # of the paper and issue #20 for details.
+
+    if (max(fm$S) > 1) {
+      stop(
+        "Elements of `S` (after symmetrizing and scaling to achieve expected ",
+        "degree) must not exceed 1 for Bernoulli graphs.",
+        call. = FALSE
+      )
+    }
+
+    fm$S <- -log(1 - fm$S)
+  }
+
   validate_directed_factor_model(fm)
 }
 
@@ -220,6 +278,9 @@ print.directed_factor_model <- function(x, ...) {
   cat("S:", dim_and_class(x$S), "\n")
   cat("Y:", dim_and_class(x$Y), "\n\n")
 
+  cat("Poisson edges:", as.character(x$poisson_edges), "\n")
+  cat("Allow self loops:", as.character(x$allow_self_loops), "\n\n")
+
   cat(
     glue("Expected edges: {round(expected_edges(x))}"),
     glue("Expected density: {round(expected_density(x), 5)}"),

R/expected-degrees.R

@@ -25,7 +25,7 @@
 #'
 #' B <- matrix(c(a,b,b,a), nrow = 2)
 #'
-#' b_model <- fastRG::sbm(n = n, k = 2, B = B,  edge_distribution = "bernoulli")
+#' b_model <- fastRG::sbm(n = n, k = 2, B = B, poisson_edges = FALSE)
 #'
 #' b_model
 #'