Add random kernel graph generator

Generator of random graphs with a specified kernel.
Includes a simple test and example to generate a G(n,p) graph using a kernel.
The method is detailed in a (currently unpublished) paper.

3 files changed, 93 insertions, 4 deletions

diff --git a/doc/source/reference/generators.rst b/doc/source/reference/generators.rst
index 14453b10..b2a3d39e 100644
--- a/doc/source/reference/generators.rst
+++ b/doc/source/reference/generators.rst
@@ -102,6 +102,7 @@ Random Graphs
    barabasi_albert_graph
    powerlaw_cluster_graph
    duplication_divergence_graph
+   random_kernel_graph
    random_lobster
    random_shell_graph
    random_powerlaw_tree
diff --git a/networkx/generators/random_graphs.py b/networkx/generators/random_graphs.py
index f4f95716..a00ee3db 100644
--- a/networkx/generators/random_graphs.py
+++ b/networkx/generators/random_graphs.py
@@ -9,7 +9,7 @@ Generators for random graphs.
 #    Pieter Swart <swart@lanl.gov>
 #    All rights reserved.
 #    BSD license.
-from collections import defaultdict
+from __future__ import division
 import itertools
 import math
 import random
@@ -17,8 +17,7 @@ import random
 import networkx as nx
 from .classic import empty_graph, path_graph, complete_graph
 from .degree_seq import degree_sequence_tree
-
-
+from collections import defaultdict
 __author__ = "\n".join(['Aric Hagberg (hagberg@lanl.gov)',
                         'Pieter Swart (swart@lanl.gov)',
                         'Dan Schult (dschult@colgate.edu)'])
@@ -39,7 +38,8 @@ __all__ = ['fast_gnp_random_graph',
            'random_lobster',
            'random_shell_graph',
            'random_powerlaw_tree',
-           'random_powerlaw_tree_sequence']
+           'random_powerlaw_tree_sequence',
+           'random_kernel_graph']
 
 
 #-------------------------------------------------------------------------
@@ -998,3 +998,82 @@ def random_powerlaw_tree_sequence(n, gamma=3, seed=None, tries=100):
 
     raise nx.NetworkXError('Exceeded max (%d) attempts for a valid tree'
                            ' sequence.' % tries)
+
+
+def random_kernel_graph(n, kernel_integral, kernel_root=None, seed=None):
+    """Return an random graph based on the specified kernel.
+
+    The algorithm chooses each of the `[n(n-1)]/2` possible edges with
+    probability specified by a kernel `\kappa(x,y)` [1]_.  The kernel 
+    `\kappa(x,y)` must be a symmetric (in `x,y`), non-negative, 
+    bounded function.
+
+    Parameters
+    ----------
+    n : int
+        The number of nodes
+    kernal_integral : function
+        Function that returns the definite integral of the kernel `\kappa(x,y)`,
+        `F(y,a,b) := \int_a^b \kappa(x,y)dx`
+    kernel_root: function (optional)
+        Function that returns the root `b` of the equation `F(y,a,b) = r`.
+        If ``None``, the root is found using :func:`scipy.optimize.brentq` 
+        (this requires SciPy).
+    seed : int, optional
+        Seed for random number generator (default=None)
+
+    Notes
+    -----
+    The kernel is specified through its definite integral which must be 
+    provided as one of the arguments. If the integral and root of the 
+    kernel integral can be found in `O(1)` time then this algorithm runs in 
+    time `O(n+m)` where m is the expected number of edges [2]_.
+
+    The nodes are set to integers from 0 to n-1.
+
+    Examples
+    --------
+    Generate an Erdős–Rényi random graph `G(n,c/n)`, with kernel
+    `\kappa(x,y)=c` where `c` is the mean expected degree.
+
+    >>> def integral(u, w, z):
+    ...     return c*(z-w)
+    >>> def root(u, w, r):
+    ...     return r/c+w
+    >>> c = 1
+    >>> graph = random_kernel_graph(1000, integral, root)
+
+    See Also
+    --------
+    gnp_random_graph
+    expected_degree_graph
+
+    References
+    ----------
+    .. [1] Bollobás, Béla,  Janson, S. and Riordan, O.
+       "The phase transition in inhomogeneous random graphs",
+       *Random Structures Algorithms*, 31, 3--122, 2007.
+
+    .. [2] Hagberg A, Lemons N (2015),
+       "Fast Generation of Sparse Random Kernel Graphs".
+       PLoS ONE 10(9): e0135177, 2015. doi:10.1371/journal.pone.0135177
+    """
+    if not seed is None:
+        random.seed(seed)
+    if kernel_root is None:
+        import scipy.optimize as optimize
+        def kernel_root(y, a, r):
+            def my_function(b):
+                return kernel_integral(y, a, b) - r
+            return optimize.brentq(my_function, a, 1)
+    graph = nx.Graph()
+    graph.add_nodes_from(range(n))
+    (i, j) = (1, 1)
+    while i < n:
+        r = -math.log(1 - random.random()) # (1-random.random()) in (0, 1]
+        if kernel_integral(i/n, j/n, 1) <= r:
+            i, j = i+1, i+1
+        else:
+            j = int(math.ceil(n*kernel_root(i/n, j/n, r)))
+            graph.add_edge(i-1, j-1)
+    return graph
diff --git a/networkx/generators/tests/test_random_graphs.py b/networkx/generators/tests/test_random_graphs.py
index 8b17b2de..631b65f0 100644
--- a/networkx/generators/tests/test_random_graphs.py
+++ b/networkx/generators/tests/test_random_graphs.py
@@ -135,3 +135,12 @@ class TestGeneratorsRandom():
         # infinite loop used to occur when a node has degree n-1 and needs to rewire
         watts_strogatz_graph(10, 9, 0.25, seed=0)
         newman_watts_strogatz_graph(10, 9, 0.5, seed=0)
+
+    def test_random_kernel_graph(self):
+        def integral(u, w, z):
+            return c*(z-w)
+        def root(u, w, r):
+            return r/c+w
+        c = 1
+        graph = random_kernel_graph(1000, integral, root)
+        assert_equal(len(graph), 1000)