The reliability polynomial (specifically, the all-terminal reliability polynomial) is, given a value \$0 \le p \le 1\$, outputs the probability that a given graph \$G\$ is connected.
Note: at the link given, the probability is given to an edge being deleted, and not to an edge existing (which is the definition we will use here).
Of course, if the graph is already not connected, then it has probability 0 that it is connected. If the graph has exactly 1 vertex and no edges, it always will be connected (so has reliability polynomial 1).
The usual "brute-force" calculation of the polynomial is by the "Factoring Theorem": for any edge \$e\$ in the graph \$G\$:
$$Rel(G; p) = p*Rel(G*e; p) + (1-p)*Rel(G/e; p)$$
where \$G*e\$ is the resulting graph of contracting the edge \$e\$, and \$G/e\$ is that of deleting the edge \$e\$.
I implemented this calculation using Sympy (for multiplying polynomials) and networkx (for using graphs). I would like to get some feedback on the structure of the code as well as readability.
Note: I actually do need MultiGraph
here because contracting an edge will introduce parallel edges, which do make a difference in the calculation.
import networkx as nx
import random
import sympy
p = sympy.symbols('p')
def relpoly(G):
H = nx.MultiGraph(G)
rel = _recursive(H)
return sympy.simplify(rel)
def _recursive(G):
# If the graph is not connected, then it has a rel poly of 0
if not nx.is_connected(G):
return sympy.Poly(0, p)
# if # edges > 0, then we perform the two subcases of the
# Factoring theorem.
if len(G.edges()) > 0:
e = random.choice(G.edges())
contracted = nx.contracted_edge(G, e, self_loops=False)
G.remove_edge(*e)
rec_deleted = _recursive(G)
rec_contracted = _recursive(contracted)
s = sympy.Poly(p)*(rec_contracted) + sympy.Poly(1-p)*(rec_deleted)
return s
# Otherwise, we only have 0 edges and 1 vertex, which is connected,
# so we return 1.
return sympy.Poly(1, p)
print(relpoly(nx.petersen_graph()))