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.

import numpy as np
import reticula as ret
import matplotlib.pyplot as plt

def lcc_size(n, p):
    state = ret.mersenne_twister()
    g = ret.random_gnp_graph[ret.int64](n, p, state)
    return len(ret.largest_connected_component(g))

n = 1000000
ps = np.linspace(0, 2./n)

lccs = [lcc_size(n, p) for p in ps]

plt.plot(ps, list(lccs))
plt.show()

This script, after doing a lot of computation, produces a figure that looks like this, minus the axis labels and other niceties:

Results of the isotropic bond percolation in random \(G(n, p)\) graphs. The evolution of the largest connected component size shows that the connectivity of the system changes rapidly around \(p = 10^{-6}\).

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().

import numpy as np
import reticula as ret
import matplotlib.pyplot as plt

from concurrent.futures import ThreadPoolExecutor
from functools import partial

def lcc_size(n, p):
    state = ret.mersenne_twister()
    g = ret.random_gnp_graph[ret.int64](n, p, state)
    return len(ret.largest_connected_component(g))

n = 1000000
ps = np.linspace(0, 2./n)

with ThreadPoolExecutor(max_workers=8) as e:
    lccs = e.map(partial(lcc_size, n), ps)

plt.plot(ps, list(lccs))
plt.show()