Isotropic bond-percolation in G(n, p) networks#

The goal of this example is to visualise what happens to the largest component size when you randomly occupying edges in a static, undirected network.

A random G(n, p) network can generates a random network with each edge occupied with probability \(p\). So the exercise boils down to generating \(G(n, p)\) networks with various values of \(p\) and plotting the largest component sizes.

 1import numpy as np
 2import reticula as ret
 3import matplotlib.pyplot as plt
 5def lcc_size(n, p):
 6    state = ret.mersenne_twister()
 7    g = ret.random_gnp_graph[ret.int64](n, p, state)
 8    return len(ret.largest_connected_component(g))
10n = 1000000
11ps = np.linspace(0, 2./n)
13lccs = [lcc_size(n, p) for p in ps]
15plt.plot(ps, list(lccs))

While this script definitly works, it is only utilising one CPU core. We can use more CPU cores (and run the code faster) using multi-processing. This, however, would be wasteful, as it entails running multiple independent python processes, each loading the standard library and Reticula in mostly separate places in the memory. The separate processes also have to communicate with each other by pickling and unpickling Python objects, which can be inefficient.

More desirable would be to perform this in multi-threaded setting. This reduces the memory-use overhead significantly, and mitigates the inter-process communication issue mentioned above as threads can access shared memory. In general, this also increases the risk of race conditions. On Reticula, network and edge types are immutable, so the risk race conditions when sharing a network object across threads is minimised.

In this case, there is no need to share anything other than the parameter \(p\) across threads. In the next example, we run the function lcc_size on 8 threads. Note that each thread creates its own instance of pseudorandom number generator using the command state = ret.mersenne_twister().

 1import numpy as np
 2import reticula as ret
 3import matplotlib.pyplot as plt
 5from concurrent.futures import ThreadPoolExecutor
 6from functools import partial
 8def lcc_size(n, p):
 9    state = ret.mersenne_twister()
10    g = ret.random_gnp_graph[ret.int64](n, p, state)
11    return len(ret.largest_connected_component(g))
13n = 1000000
14ps = np.linspace(0, 2./n)
16with ThreadPoolExecutor(max_workers=8) as e:
17    lccs =, n), ps)
19plt.plot(ps, list(lccs))