Reliability polynomial calculation

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) = pRel(Ge; 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,

Some suggestions:

-
_recursive isn't a particularly descriptive name. A lot of things can be recursive. _relpoly_recursive may be a better name.

-
You can factor out the Poly call

s = sympy.Poly(p*rec_contracted + (1-p)*rec_deleted, p)

also note that I've removed the redundant parentheses.

-
Calling simplify on the result is completely useless, since rel is a Poly and Poly instances are always stored in canonical (expanded) form.

-
What is the value of going through the edges randomly? Don't you end up going through all of them.

-
As I noted in a comment, modifying the graph in-place seems problematic. I can't tell for sure just by looking at the code if it does the right thing or not. It seems that you should rather move the

H = nx.MultiGraph(G)

line to the top of the recursive function, and operate on H instead of G in the function.

StackExchange Code Review Q#131709, answer score: 3