Generate graphs from a stochastic block model (#2).

Conflicts:
	src/Makefile.am

include/igraph_games.h

@@ -184,6 +184,11 @@ int igraph_k_regular_game(igraph_t *graph,
     igraph_integer_t no_of_nodes, igraph_integer_t k,
     igraph_bool_t directed, igraph_bool_t multiple);
 
+int igraph_sbm_game(igraph_t *graph, igraph_integer_t n, 
+		    const igraph_matrix_t *pref_matrix,
+		    const igraph_vector_int_t *block_sizes,
+		    igraph_bool_t directed, igraph_bool_t loops);
+
 __END_DECLS
 
 #endif

interfaces/R/igraph/inst/tests/test_sbm.game.R

@@ -0,0 +1,29 @@
+
+context("Stochastic block models")
+
+test_that("Generating stochastic block models works", {
+
+  library(igraph)
+  pm <- matrix(1, nrow=2, ncol=2)
+  bs <- c(4,6)
+  g1 <- sbm.game(10, pref.matrix=pm, block.sizes=bs,
+                 directed=FALSE, loops=FALSE)
+  expect_that(graph.isomorphic(g1, graph.full(10, directed=FALSE, loops=FALSE)),
+              is_true())
+
+  g2 <- sbm.game(10, pref.matrix=pm, block.sizes=bs,
+                directed=FALSE, loops=TRUE)
+  g2x <- graph.full(10, directed=FALSE, loops=TRUE)
+  expect_that(g2[sparse=FALSE], equals(g2x[sparse=FALSE]))
+
+  g3 <- sbm.game(10, pref.matrix=pm, block.sizes=bs,
+                 directed=TRUE, loops=FALSE)
+  g3x <- graph.full(10, directed=TRUE, loops=FALSE)
+  expect_that(g3[sparse=FALSE], equals(g3x[sparse=FALSE]))
+
+  g4 <- sbm.game(10, pref.matrix=pm, block.sizes=bs,
+                 directed=TRUE, loops=TRUE)
+  g4x <- graph.full(10, directed=TRUE, loops=TRUE)
+  expect_that(g4[sparse=FALSE], equals(g4x[sparse=FALSE]))
+
+})

interfaces/R/types-C.def

@@ -243,6 +243,20 @@ VECTOR_LONG:
          igraph_vector_long_destroy(&%C%); \
          IGRAPH_FINALLY_CLEAN(1);
 
+VECTOR_INT:
+  CTYPE: igraph_vector_int_t
+  CALL: &%C%
+  INCONV:
+    IN: R_SEXP_to_vector_int(%I%, &%C%);
+    OUT: if (0 != igraph_vector_int_init(&%C%, 0)) { \
+             igraph_error("", __FILE__, __LINE__, IGRAPH_ENOMEM); \
+         } \
+         IGRAPH_FINALLY(igraph_vector_int_destroy, &%C%);
+  OUTCONV:
+    OUT: PROTECT(%I%=R_igraph_vector_int_to_SEXP(&%C%)); \
+         igraph_vector_int_destroy(&%C%); \
+         IGRAPH_FINALLY_CLEAN(1);
+
 VECTOR_LONG_M1:
   CTYPE: igraph_vector_long_t
   CALL: &%C%

interfaces/R/types-R.def

@@ -37,6 +37,9 @@ VERTEXINDEX_OR_0:
 VECTOR_LONG:
     INCONV: %I% <- as.numeric(%I%)
 
+VECTOR_INT:
+    INCONV: %I% <- as.integer(%I%)
+
 VECTORM1:
     INCONV: %I% <- as.numeric(%I%)-1
 

interfaces/functions.def

@@ -311,6 +311,14 @@ igraph_k_regular_game:
         GATTR-R: name IS k-regular graph
         GATTR-PARAM-R: k
 
+igraph_sbm_game:
+        PARAMS: OUT GRAPH graph, INTEGER n, MATRIX pref_matrix, \
+                VECTOR_INT block_sizes, BOOLEAN directed=False, \
+                BOOLEAN loops=False
+        NAME-R: sbm.game
+        GATTR-R: name IS Stochastic block-model
+        GATTR-PARAM-R: loops
+
 #######################################
 # Basic query functions
 #######################################

src/Makefile.am

@@ -262,7 +262,8 @@ SOURCES = 		     basic_query.c games.c cocitation.c iterators.c \
 			     scg_kmeans.c scg_utils.c scg_optimal_method.c \
 			     qsort.c qsort_r.c types.c lad.c hacks.c \
 			     embedding.c scan.c qsort.c qsort_r.c types.c \
-			     lad.c hacks.c triangles.c glet.c maximal_cliques.c
+			     lad.c hacks.c triangles.c glet.c \
+			     maximal_cliques.c sbm.c
 
 libigraph_la_SOURCES       = $(SOURCES) $(HEADERS_PRIVATE)
 libigraph_la_CFLAGS        = -I$(top_srcdir)/include -I$(top_builddir)/include $(WARNING_CFLAGS)

src/sbm.c

@@ -0,0 +1,182 @@
+/* -*- mode: C -*-  */
+/* vim:set ts=8 sw=2 sts=2 et: */
+/* 
+   IGraph R library.
+   Copyright (C) 2003-2013  Gabor Csardi <csardi.gabor@gmail.com>
+   334 Harvard street, Cambridge, MA 02139 USA
+   
+   This program is free software; you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation; either version 2 of the License, or
+   (at your option) any later version.
+   
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+   
+   You should have received a copy of the GNU General Public License
+   along with this program; if not, write to the Free Software
+   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
+   02110-1301 USA
+
+*/
+
+#include "igraph_interface.h"
+#include "igraph_vector.h"
+#include "igraph_matrix.h"
+#include "igraph_random.h"
+#include "igraph_constructors.h"
+#include "igraph_games.h"
+
+/**
+ * \function igraph_sbm_game
+ */
+
+int igraph_sbm_game(igraph_t *graph, igraph_integer_t n, 
+		    const igraph_matrix_t *pref_matrix,
+		    const igraph_vector_int_t *block_sizes,
+		    igraph_bool_t directed, igraph_bool_t loops) {
+
+  int no_blocks=igraph_matrix_nrow(pref_matrix);
+  int from, to, fromoff=0;
+  igraph_real_t minp, maxp;
+  igraph_vector_t edges;
+
+  /* ------------------------------------------------------------ */
+  /* Check arguments                                              */
+  /* ------------------------------------------------------------ */
+
+  if (igraph_matrix_ncol(pref_matrix) != no_blocks) {
+    IGRAPH_ERROR("Preference matrix is not square", 
+		 IGRAPH_NONSQUARE);
+  }
+
+  igraph_matrix_minmax(pref_matrix, &minp, &maxp);
+  if (minp < 0 || maxp > 1) { 
+    IGRAPH_ERROR("Connection probabilities must in [0,1]", IGRAPH_EINVAL);
+  }
+
+  if (n < 0) { 
+    IGRAPH_ERROR("Number of vertices must be non-negative", IGRAPH_EINVAL);
+  }
+
+  if (!directed && !igraph_matrix_is_symmetric(pref_matrix)) {
+    IGRAPH_ERROR("Preference matrix must be symmetric for undirected graphs",
+		 IGRAPH_EINVAL);
+  }
+
+  if (igraph_vector_int_size(block_sizes) != no_blocks) {
+    IGRAPH_ERROR("Invalid block size vector length", IGRAPH_EINVAL);
+  }
+
+  if (igraph_vector_int_min(block_sizes) < 0) {
+    IGRAPH_ERROR("Block size must be non-negative", IGRAPH_EINVAL);
+  }
+
+  if (igraph_vector_int_sum(block_sizes) != n) {
+    IGRAPH_ERROR("Block sizes must sum up to number of vertices", 
+		 IGRAPH_EINVAL);
+  }
+
+  IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);  
+
+  RNG_BEGIN();
+
+  for (from = 0; from < no_blocks; from++) {
+    int fromsize = VECTOR(*block_sizes)[from];
+    int start = directed ? 0 : from;
+    int i, tooff=0;
+    for (i=0; i<start; i++) {
+      tooff += VECTOR(*block_sizes)[i];
+    }
+    for (to = start; to < no_blocks; to++) {
+      int tosize = VECTOR(*block_sizes)[to];
+      igraph_real_t prob=MATRIX(*pref_matrix, from, to);
+      double maxedges, last=RNG_GEOM(prob);
+      if (directed && loops) {
+	maxedges = fromsize * tosize;
+	while (last < maxedges) {
+	  int vto=floor(last/fromsize);
+	  int vfrom=last - (igraph_real_t)vto * fromsize;
+	  igraph_vector_push_back(&edges, fromoff + vfrom);
+	  igraph_vector_push_back(&edges, tooff + vto);
+	  last += RNG_GEOM(prob);
+	  last += 1;
+	}
+      } else if (directed && !loops && from!=to) {
+	maxedges = fromsize * tosize;
+	while (last < maxedges) {
+	  int vto=floor(last/fromsize);
+	  int vfrom=last - (igraph_real_t)vto * fromsize;
+	  igraph_vector_push_back(&edges, fromoff + vfrom);
+	  igraph_vector_push_back(&edges, tooff + vto);
+	  last += RNG_GEOM(prob);
+	  last += 1;
+	}	
+      } else if (directed && !loops && from==to) {
+	maxedges = fromsize * (fromsize-1);
+	while (last < maxedges) {
+	  int vto=floor(last/fromsize);
+	  int vfrom=last - (igraph_real_t)vto * fromsize;
+	  if (vfrom == vto) { vto=fromsize-1; }
+	  igraph_vector_push_back(&edges, fromoff + vfrom);
+	  igraph_vector_push_back(&edges, tooff + vto);
+	  last += RNG_GEOM(prob);
+	  last += 1;
+	}
+      } else if (!directed && loops && from!=to) {
+	maxedges = fromsize * tosize;	
+	while (last < maxedges) {
+	  int vto=floor(last/fromsize);
+	  int vfrom=last - (igraph_real_t)vto * fromsize;
+	  igraph_vector_push_back(&edges, fromoff + vfrom);
+	  igraph_vector_push_back(&edges, tooff + vto);
+	  last += RNG_GEOM(prob);
+	  last += 1;
+	}
+      } else if (!directed && loops && from==to) {
+	maxedges = fromsize * (fromsize+1) / 2.0;
+	while (last < maxedges) {
+	  long int vto=floor((sqrt(8*last+1)-1)/2);
+	  long int vfrom=last-(((igraph_real_t)vto)*(vto+1))/2;
+	  igraph_vector_push_back(&edges, fromoff + vfrom);
+	  igraph_vector_push_back(&edges, tooff + vto);
+	  last += RNG_GEOM(prob);
+	  last += 1;
+	}	
+      } else if (!directed && !loops && from!=to) {
+	maxedges = fromsize * tosize;	
+	while (last < maxedges) {
+	  int vto=floor(last/fromsize);
+	  int vfrom=last - (igraph_real_t)vto * fromsize;
+	  igraph_vector_push_back(&edges, fromoff + vfrom);
+	  igraph_vector_push_back(&edges, tooff + vto);
+	  last += RNG_GEOM(prob);
+	  last += 1;
+	}
+      } else /*!directed && !loops && from==to */ {
+	maxedges = fromsize * (fromsize-1) / 2.0;
+	while (last < maxedges) {
+	  int vto=floor((sqrt(8*last+1)+1)/2);
+	  int vfrom=last-(((igraph_real_t)vto)*(vto-1))/2;
+	  igraph_vector_push_back(&edges, fromoff + vfrom);
+	  igraph_vector_push_back(&edges, tooff + vto);
+	  last += RNG_GEOM(prob);
+	  last += 1;
+	}
+      }
+
+      tooff += tosize;
+    }
+    fromoff += fromsize;
+  }
+
+  RNG_END();
+
+  igraph_create(graph, &edges, n, directed);
+  igraph_vector_destroy(&edges);
+  IGRAPH_FINALLY_CLEAN(1);
+
+  return 0;
+}