fix: synchronize R RNG state in hand-written C glue (#2211)

src/rinterface_extra.c

@@ -132,6 +132,15 @@ void igraph_vector_int_list_destroy_pv(void *pv_ptr)
   igraph_vector_int_list_destroy((igraph_vector_int_list_t*) pv_ptr);
 }
 
+/* Restores R's RNG state from the igraph finally stack, so it runs even when an
+ * igraph call longjmps out via IGRAPH_R_CHECK. Pair every GetRNGstate() call
+ * with IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL), then on the success path call
+ * PutRNGstate() and IGRAPH_FINALLY_CLEAN(1). */
+static void Rx_PutRNGstate_pv(void *pv_ptr) {
+  (void) pv_ptr;
+  PutRNGstate();
+}
+
 igraph_error_t Rw_get_int_scalar(SEXP sexp, R_xlen_t index, igraph_integer_t *res)
 {
   if (Rf_xlength(sexp) <= index)
@@ -1672,7 +1681,7 @@ SEXP Rx_igraph_ac_random_numeric(SEXP attr,
   PROTECT(attr2=AS_NUMERIC(attr));
   PROTECT(res=NEW_NUMERIC(len));
 
-  RNG_BEGIN();
+  GetRNGstate();
 
   for (igraph_integer_t i=0; i<len; i++) {
     igraph_vector_int_t *v=igraph_vector_int_list_get_ptr(merges, i);
@@ -1687,7 +1696,7 @@ SEXP Rx_igraph_ac_random_numeric(SEXP attr,
     }
   }
 
-  RNG_END();
+  PutRNGstate();
 
   UNPROTECT(2);
   return res;
@@ -4304,7 +4313,11 @@ SEXP Rx_igraph_layout_fruchterman_reingold(SEXP graph, SEXP coords,
   if (!Rf_isNull(miny)) { Rz_SEXP_to_vector(miny, &c_miny); }
   if (!Rf_isNull(maxy)) { Rz_SEXP_to_vector(maxy, &c_maxy); }
   /* Call igraph */
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
   IGRAPH_R_CHECK(igraph_layout_fruchterman_reingold(&c_graph, &c_coords, !Rf_isNull(coords), c_niter, c_start_temp, c_grid, (Rf_isNull(weights) ? 0 : &c_weights), (Rf_isNull(minx) ? 0 : &c_minx), (Rf_isNull(maxx) ? 0 : &c_maxx), (Rf_isNull(miny) ? 0 : &c_miny), (Rf_isNull(maxy) ? 0 : &c_maxy)));
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
 
   /* Convert output */
   PROTECT(coords=Ry_igraph_matrix_to_SEXP(&c_coords));
@@ -4356,7 +4369,11 @@ SEXP Rx_igraph_layout_fruchterman_reingold_3d(SEXP graph, SEXP coords,
   if (!Rf_isNull(minz)) { Rz_SEXP_to_vector(minz, &c_minz); }
   if (!Rf_isNull(maxz)) { Rz_SEXP_to_vector(maxz, &c_maxz); }
   /* Call igraph */
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
   IGRAPH_R_CHECK(igraph_layout_fruchterman_reingold_3d(&c_graph, &c_coords, !Rf_isNull(coords), c_niter, c_start_temp, (Rf_isNull(weights) ? 0 : &c_weights), (Rf_isNull(minx) ? 0 : &c_minx), (Rf_isNull(maxx) ? 0 : &c_maxx), (Rf_isNull(miny) ? 0 : &c_miny), (Rf_isNull(maxy) ? 0 : &c_maxy), (Rf_isNull(minz) ? 0 : &c_minz), (Rf_isNull(maxz) ? 0 : &c_maxz)));
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
 
   /* Convert output */
   PROTECT(coords=Ry_igraph_matrix_to_SEXP(&c_coords));
@@ -4405,7 +4422,11 @@ SEXP Rx_igraph_layout_kamada_kawai(SEXP graph, SEXP coords, SEXP maxiter,
   if (!Rf_isNull(miny)) { Rz_SEXP_to_vector(miny, &c_miny); }
   if (!Rf_isNull(maxy)) { Rz_SEXP_to_vector(maxy, &c_maxy); }
   /* Call igraph */
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
   IGRAPH_R_CHECK(igraph_layout_kamada_kawai(&c_graph, &c_coords, !Rf_isNull(coords), c_maxiter, c_epsilon, c_kkconst, (Rf_isNull(weights) ? 0 : &c_weights), (Rf_isNull(minx) ? 0 : &c_minx), (Rf_isNull(maxx) ? 0 : &c_maxx), (Rf_isNull(miny) ? 0 : &c_miny), (Rf_isNull(maxy) ? 0 : &c_maxy)));
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
 
                                         /* Convert output */
   PROTECT(coords=Ry_igraph_matrix_to_SEXP(&c_coords));
@@ -4460,7 +4481,11 @@ SEXP Rx_igraph_layout_kamada_kawai_3d(SEXP graph, SEXP coords, SEXP maxiter,
   if (!Rf_isNull(minz)) { Rz_SEXP_to_vector(minz, &c_minz); }
   if (!Rf_isNull(maxz)) { Rz_SEXP_to_vector(maxz, &c_maxz); }
   /* Call igraph */
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
   IGRAPH_R_CHECK(igraph_layout_kamada_kawai_3d(&c_graph, &c_coords, !Rf_isNull(coords), c_maxiter, c_epsilon, c_kkconst, (Rf_isNull(weights) ? 0 : &c_weights), (Rf_isNull(minx) ? 0 : &c_minx), (Rf_isNull(maxx) ? 0 : &c_maxx), (Rf_isNull(miny) ? 0 : &c_miny), (Rf_isNull(maxy) ? 0 : &c_maxy), (Rf_isNull(minz) ? 0 : &c_minz), (Rf_isNull(maxz) ? 0 : &c_maxz)));
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
 
                                         /* Convert output */
   PROTECT(coords=Ry_igraph_matrix_to_SEXP(&c_coords));
@@ -4492,7 +4517,11 @@ SEXP Rx_igraph_layout_graphopt(SEXP graph, SEXP pniter, SEXP pcharge,
   } else {
     Rz_SEXP_to_igraph_matrix_copy(start, &res);
   }
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
   IGRAPH_R_CHECK(igraph_layout_graphopt(&g, &res, niter, charge, mass, spring_length, spring_constant, max_sa_movement, !Rf_isNull(start)));
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
   PROTECT(result=Ry_igraph_matrix_to_SEXP(&res));
   igraph_matrix_destroy(&res);
 
@@ -4517,7 +4546,11 @@ SEXP Rx_igraph_layout_lgl(SEXP graph, SEXP pmaxiter, SEXP pmaxdelta,
 
   Rz_SEXP_to_igraph(graph, &g);
   igraph_matrix_init(&res, 0, 0);
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
   IGRAPH_R_CHECK(igraph_layout_lgl(&g, &res, maxiter, maxdelta, area, coolexp, repulserad, cellsize, root));
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
   PROTECT(result=Ry_igraph_matrix_to_SEXP(&res));
   igraph_matrix_destroy(&res);
 
@@ -5049,7 +5082,11 @@ SEXP Rx_igraph_layout_merge_dla(SEXP graphs, SEXP layouts) {
     igraph_matrix_update(dest, &source);
   }
   igraph_matrix_init(&res, 0, 0);
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
   IGRAPH_R_CHECK(igraph_layout_merge_dla(&graphvec, &matrixlist, &res));
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
   igraph_vector_ptr_destroy(&graphvec);
   igraph_matrix_list_destroy(&matrixlist);
   PROTECT(result=Ry_igraph_matrix_to_SEXP(&res));
@@ -6132,7 +6169,11 @@ SEXP Rx_igraph_walktrap_community(SEXP graph, SEXP pweights,
   igraph_vector_init(&modularity, 0);
   igraph_vector_int_init(&membership, 0);
 
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
   IGRAPH_R_CHECK(igraph_community_walktrap(&g, ppweights, steps, &merges, &modularity, &membership));
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
 
   PROTECT(result=NEW_LIST(3));
   if (LOGICAL(pmerges)[0]) {
@@ -7841,9 +7882,21 @@ SEXP Rx_igraph_vcount(SEXP graph) {
 }
 
 SEXP Rx_igraph_layout_drl(SEXP graph, SEXP res, SEXP use_seed, SEXP options, SEXP weights) {
-  return R_igraph_layout_drl(graph, res, use_seed, options, weights);
+  SEXP result;
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
+  result = R_igraph_layout_drl(graph, res, use_seed, options, weights);
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
+  return result;
 }
 
 SEXP Rx_igraph_layout_drl_3d(SEXP graph, SEXP res, SEXP use_seed, SEXP options, SEXP weights) {
-  return R_igraph_layout_drl_3d(graph, res, use_seed, options, weights);
+  SEXP result;
+  GetRNGstate();
+  IGRAPH_FINALLY(Rx_PutRNGstate_pv, NULL);
+  result = R_igraph_layout_drl_3d(graph, res, use_seed, options, weights);
+  PutRNGstate();
+  IGRAPH_FINALLY_CLEAN(1);
+  return result;
 }

tests/testthat/test-layout.R

@@ -1,17 +1,70 @@
+check_matrix <- function(mat, nrow = 1, ncol = 2) {
+  expect_equal(dim(mat), c(nrow, ncol))
+  if (nrow > 0) {
+    expect_true(all(is.finite(mat)))
+  }
+}
+
 test_that("layout_with_fr() works", {
   skip_on_os("solaris")
 
+  g <- make_ring(10)
   withr::with_seed(42, {
-    g <- make_ring(10)
     l <- layout_with_fr(g, niter = 50, start.temp = sqrt(10) / 10)
   })
-  expect_equal(sum(l), 4.57228, tolerance = 0.1)
+  check_matrix(l, nrow = 10, ncol = 2)
 
+  g <- make_star(30)
   withr::with_seed(42, {
-    g <- make_star(30)
     l <- layout_with_fr(g, niter = 500, dim = 3, start.temp = 20)
   })
-  expect_equal(sum(l), -170.9312, tolerance = 0.1)
+  check_matrix(l, nrow = 30, ncol = 3)
+})
+
+test_that("stochastic layouts are reproducible with set.seed()", {
+  g <- make_ring(20) + make_star(10, mode = "undirected")
+
+  expect_identical(
+    withr::with_seed(42, layout_with_fr(g, niter = 50)),
+    withr::with_seed(42, layout_with_fr(g, niter = 50))
+  )
+  expect_identical(
+    withr::with_seed(42, layout_with_fr(g, niter = 50, dim = 3)),
+    withr::with_seed(42, layout_with_fr(g, niter = 50, dim = 3))
+  )
+  expect_identical(
+    withr::with_seed(42, layout_with_kk(g, maxiter = 50)),
+    withr::with_seed(42, layout_with_kk(g, maxiter = 50))
+  )
+  expect_identical(
+    withr::with_seed(42, layout_with_kk(g, maxiter = 50, dim = 3)),
+    withr::with_seed(42, layout_with_kk(g, maxiter = 50, dim = 3))
+  )
+  expect_identical(
+    withr::with_seed(42, layout_with_dh(g, maxiter = 5)),
+    withr::with_seed(42, layout_with_dh(g, maxiter = 5))
+  )
+  expect_identical(
+    withr::with_seed(42, layout_with_gem(g, maxiter = 50)),
+    withr::with_seed(42, layout_with_gem(g, maxiter = 50))
+  )
+  expect_identical(
+    withr::with_seed(42, layout_with_graphopt(g, niter = 50)),
+    withr::with_seed(42, layout_with_graphopt(g, niter = 50))
+  )
+  g_connected <- make_ring(20)
+  expect_identical(
+    withr::with_seed(42, layout_with_lgl(g_connected, maxiter = 50)),
+    withr::with_seed(42, layout_with_lgl(g_connected, maxiter = 50))
+  )
+  expect_identical(
+    withr::with_seed(42, layout_with_drl(g)),
+    withr::with_seed(42, layout_with_drl(g))
+  )
+  expect_identical(
+    withr::with_seed(42, layout_with_drl(g, dim = 3)),
+    withr::with_seed(42, layout_with_drl(g, dim = 3))
+  )
 })
 
 test_that("layout_with_fr() deprecated argument", {
@@ -103,13 +156,6 @@ test_that("layout algorithms work for null graphs", {
 test_that("layout algorithms work for singleton graphs", {
   g <- make_empty_graph(1)
 
-  check_matrix <- function(mat, nrow = 1, ncol = 2) {
-    expect_equal(dim(mat), c(nrow, ncol))
-    if (nrow > 0) {
-      expect_true(all(is.finite(mat)))
-    }
-  }
-
   expect_silent(layout_as_tree(g))
   check_matrix(layout_as_tree(g))
 

tests/testthat/test-rng-seeding.R

@@ -0,0 +1,51 @@
+# Compare two graphs by their edge structure (ignoring the auto-generated
+# `myid` UUID and `igraph` pointer that differ per object).
+expect_same_graph <- function(g1, g2) {
+  expect_identical(as_edgelist(g1), as_edgelist(g2))
+  expect_identical(vcount(g1), vcount(g2))
+  expect_identical(is_directed(g1), is_directed(g2))
+}
+
+test_that("random graph generators are reproducible with set.seed()", {
+  expect_same_graph(
+    withr::with_seed(42, sample_gnp(50, 0.1)),
+    withr::with_seed(42, sample_gnp(50, 0.1))
+  )
+  expect_same_graph(
+    withr::with_seed(42, sample_gnm(50, 100)),
+    withr::with_seed(42, sample_gnm(50, 100))
+  )
+  expect_same_graph(
+    withr::with_seed(42, sample_pa(50, m = 2)),
+    withr::with_seed(42, sample_pa(50, m = 2))
+  )
+  expect_same_graph(
+    withr::with_seed(42, sample_smallworld(1, 20, 5, 0.05)),
+    withr::with_seed(42, sample_smallworld(1, 20, 5, 0.05))
+  )
+})
+
+test_that("random_walk() is reproducible with set.seed()", {
+  g <- make_ring(20)
+  expect_identical(
+    withr::with_seed(42, as.integer(random_walk(g, start = 1, steps = 10))),
+    withr::with_seed(42, as.integer(random_walk(g, start = 1, steps = 10)))
+  )
+})
+
+test_that("community detection is reproducible with set.seed()", {
+  g <- make_graph("Zachary")
+
+  expect_identical(
+    withr::with_seed(42, membership(cluster_walktrap(g))),
+    withr::with_seed(42, membership(cluster_walktrap(g)))
+  )
+  expect_identical(
+    withr::with_seed(42, membership(cluster_louvain(g))),
+    withr::with_seed(42, membership(cluster_louvain(g)))
+  )
+  expect_identical(
+    withr::with_seed(42, membership(cluster_leiden(g))),
+    withr::with_seed(42, membership(cluster_leiden(g)))
+  )
+})