Generating Steiner tree using Network X in Python?

I'm trying to generate a Steiner tree in Python using NetworkX. I have a range of geolocated nodes saved in the GeoJSON format, allowing me to visually validate the outputs in QGIS. A radially deployed straight-line network looks as follows:

This isn't very realistic however for a broadband network. A Steiner tree would reduce the cost for connecting all premises. Here is my artistic attempt for how I imagine it should look:

So far I've setup a function which loads in each GeoJSON and measures the distance from the central node for the edge weight:

``````import networkx as nx
import matplotlib.pyplot as plt
from shapely.geometry import Point

G = nx.Graph()

for dist_point in distribution_point_data:
dist_point_name = [dist_point['properties']['name']]
dist_point_coordinates = dist_point['geometry']['coordinates']

for premises in premises_data:
premises_name = [premises ['properties']['name']]
premises_coordinates = premises['geometry']['coordinates']

for dist_point in distribution_point_data:
for premises in premises_data:
dist_point_geom = Point(dist_point['geometry']['coordinates'])
premises_geom = Point(premises['geometry']['coordinates'])
distance = round(dist_point_geom.distance(premises_geom),2)

return G
``````

I then use the source code for the Network X Steiner tree (and metric_closure function) as follows:

``````from itertools import combinations, chain
from networkx.utils import pairwise, not_implemented_for

def metric_closure(G, weight='weight'):
"""  Return the metric closure of a graph.

The metric closure of a graph *G* is the complete graph in which each edge
is weighted by the shortest path distance between the nodes in *G* .

Parameters
----------
G : NetworkX graph

Returns
-------
NetworkX graph
Metric closure of the graph `G`.

"""
M = nx.Graph()

Gnodes = set(G)

# check for connected graph while processing first node
all_paths_iter = nx.all_pairs_dijkstra(G, weight=weight)
u, (distance, path) = next(all_paths_iter)
if Gnodes - set(distance):
msg = "G is not a connected graph. metric_closure is not defined."
raise nx.NetworkXError(msg)
Gnodes.remove(u)
for v in Gnodes:

# first node done -- now process the rest
for u, (distance, path) in all_paths_iter:
Gnodes.remove(u)
for v in Gnodes:

return M

def steiner_tree(G, terminal_nodes, weight='weight'):
""" Return an approximation to the minimum Steiner tree of a graph.

Parameters
----------
G : NetworkX graph

terminal_nodes : list
A list of terminal nodes for which minimum steiner tree is
to be found.

Returns
-------
NetworkX graph
Approximation to the minimum steiner tree of `G` induced by
`terminal_nodes` .

Notes
-----
Steiner tree can be approximated by computing the minimum spanning
tree of the subgraph of the metric closure of the graph induced by the
terminal nodes, where the metric closure of *G* is the complete graph in
which each edge is weighted by the shortest path distance between the
nodes in *G* .
This algorithm produces a tree whose weight is within a (2 - (2 / t))
factor of the weight of the optimal Steiner tree where *t* is number of
terminal nodes.

"""
# M is the subgraph of the metric closure induced by the terminal nodes of
# G.
M = metric_closure(G, weight=weight)
# Use the 'distance' attribute of each edge provided by the metric closure
# graph.
H = M.subgraph(terminal_nodes)
mst_edges = nx.minimum_spanning_edges(H, weight='distance', data=True)
# Create an iterator over each edge in each shortest path; repeats are okay
edges = chain.from_iterable(pairwise(d['path']) for u, v, d in mst_edges)
T = G.edge_subgraph(edges)
return T
``````

Here is the GeoJSON data for a reproducible example for the premises on the far right of the pictures above:

``````GEOJSON_PREMISES = [
{
'type': "Feature",
'geometry':{
"type": "Point",
"coordinates": [0.129960,52.220977]
},
'properties': {
'name': 'premises_1',
}
},
{
'type': "Feature",
'geometry':{
"type": "Point",
"coordinates": [0.130056,52.220895]
},
'properties': {
'name': 'premises_2',
}
},
{
'type': "Feature",
'geometry':{
"type": "Point",
"coordinates": [0.130318,52.220782]
},
'properties': {
'name': 'premises_3',
}
},
]
GEOJSON_DISTRIBUTION_POINT = [
{
'type': "Feature",
'geometry':{
"type": "Point",
"coordinates": [0.129829,52.220802]
},
'properties': {
'name': 'distribution_point'
}
}
]
``````

And I run the following to generate a graph, and then (supposedly) a Steiner tree:

``````graph = load_nodes(GEOJSON_PREMISES, GEOJSON_DISTRIBUTION_POINT)
graph2 = steiner_tree(graph, graph.nodes, weight = 'distance')

plt.figure()
nx.draw_networkx(graph2)
plt.show()
``````

This is the output I get, which was not what I was expecting:

Firstly, what am I doing wrong for generating a Steiner tree graph which allows the network to run between the premises points, like in my artistic impression?

Secondly, why does this output plot not use the geographic coordinates of the nodes? (I think I'm configuring this incorrectly because I think Network X can deal with node coordinates as well as weights).

Answers should be for Python and not QGIS/ArcMap, as this needs to eventually be run on a large-scale computer cluster.

• You are talking about Steiner tree, not mst Commented Jan 1, 2019 at 20:33
• @FelixIP Question now updated to reflect this. Commented Jan 2, 2019 at 11:35
• Number of nodes connected by Steiner tree << total nodes in network. Commented Jan 2, 2019 at 19:05
• Missing node at intersection of p2-p3 and p1-source. Commented Jan 2, 2019 at 19:20

Steiner tree connects some(!) of the network's nodes (terminals) shown as selected nodes:

However don't get over excited about this feature of networkX, there is a good reason they called it "approximation.steinertree.steiner_tree", e.g. tree length in below picture is 380 m (4%) less than in the first one:

Computation of Steiner tree is ginormous task, it involves search for so called Steiners' points. Result shown on second picture is my poor attempt to improve default outcome (1st picture) by considering 3 terminals at a time in a search for such points.

NetworkX algorithm not searching for that points. It assumes that Steiner tree is minimum spanning tree of graph represented by OD matrix for terminals. That is pragmatic approach working well enough with dense terminals, where chances of Steiner points sitting on a shortest path between 2 terminals are big, e.g. difference in tree length shown on 2 pictures below is only 2%:

To improve your results I second @Spacedman and suggest to add all possible edges between premises to your graph.

The minimum spanning tree is computed on the input graph. Your input graph is the star network from the distribution point to the three premises - it doesn't contain edges from the premises to the other premises, so the output MST can't have those links in it.

You want to create a full graph with all edges between all the pairs of locations - distribution point and premises points. Then the MST will span over that.

• I've taken on board the comment from FelixIP that this should be a Steiner Tree, and have edited the question to reflect this. In the code posted, the Steiner Tree function initially obtains the full metric closure for the graph, as you suggest. However, even after doing this I am still unable to generate the Steiner graph. Help would be much appreciated. Commented Jan 2, 2019 at 11:40