Ga naar hoofdinhoud

Graafoptimalisatie met NetworkX in Python

Deze NetworkX-tutorial laat je zien hoe je graafoptimalisatie in Python doet door het Chinese Postman-probleem op te lossen in Python.
· 15 min lezen

Verken met AI

datacamp-graphic

Met deze tutorial pak je een bekend probleem uit de grafentheorie aan: het Chinese postbodeprobleem. Sommige onderdelen van het algoritme zijn conceptueel eenvoudig, maar blijken rekentechnisch pittig. Voor deze tutorial heb je echter alleen wat basiskennis van Python nodig: geen zware wiskunde, informatica of grafentheorie-achtergrond vereist.

Deze tutorial behandelt eerst de basisbouwstenen van grafen (nodes, edges, paden, enz.) en lost het probleem op een echte graaf op (het trailnetwerk van een state park) met behulp van de NetworkX-bibliotheek in Python. Je focust op de kernconcepten en de implementatie. Voor de geïnteresseerde lezer is verdere lectuur over de optimalisatie onder de motorkap voorzien.

Motiverende graafoptimalisatie

Het probleem

Je hebt vast wel eens gehoord van het probleem van de rondreizende verkoper, dat neerkomt op het vinden van de kortste route (bijvoorbeeld wegen) die een reeks nodes (bijvoorbeeld steden) verbindt. Het minder bekende Chinese postbodeprobleem (CPP), ook wel het Route Inspection- of Arc Routing-probleem genoemd, lijkt daar sterk op. Het doel van het CPP is om het kortste pad te vinden dat alle verbindingen (wegen) op een graaf minstens één keer aflegt. Als dat kan zonder ergens twee keer over dezelfde weg terug te moeten, mooi; dat is het ideale scenario en dan is het probleem vrij simpel. Maar als sommige wegen meerdere keren genomen moeten worden, heb je wat wiskunde nodig om de kortste route te vinden die elke weg minstens één keer aandoet met de laagste totale afstand.

Persoonlijke motivatie

(Het volgende is een persoonlijke noot: zoetsappig, gevat en 100% niet nodig om graafoptimalisatie in Python te leren)

Ik had een toepassing in het echt voor dit probleem: de rang van Giantmaster Marathoner behalen.

Wat is een Giantmaster? Een Giantmaster is iemand (hond of mens) die in zijn of haar leven elk pad van Sleeping Giant State Park in Hamden, CT (buur van mijn geboorteplaats Wallingford) heeft gelopen. Een Giantmaster Marathoner is iemand die al deze paden op één dag heeft gelopen.

Dankzij de nauwgezette administratie van de Sleeping Giant Park Association is de volledige lijst van Giantmasters en hun niveau van Giantmastering hier te vinden. Ik moet toegeven dat dit me behoorlijk motiveerde om dit zijproject te starten en de paden op te gaan. Terwijl ikzelf Giantmaster-status behaalde in de winter van 2006, toen ik een beginnende jonge vrijwilliger was bij de Sleeping Giant Trail Crew (wat ik met plezier terugzag in het SG-archief), zijn er sindsdien nieuwe uitdagingen bijgekomen. Hoewel de 12-maanden- en 4-seizoenen-Giantmaster-categorieën indrukwekkend en aantrekkelijk zijn, zouden ze ook meer reizen vereisen van mijn huidige thuis (DC) naar mijn vormende thuis (CT) dan ik redelijkerwijs kon opbrengen... en ze zijn minder interessant voor graafoptimalisatie, dus dan wordt het de Giantmaster Marathon!

Ter referentie vind je hieronder de kaart van de paden van Sleeping Giant:

Introductie in grafen

Het fijne aan grafen is dat de concepten en terminologie over het algemeen intuïtief zijn. Toch hier wat basisjargon:

Grafen zijn structuren die relaties tussen objecten in kaart brengen. De objecten noem je in deze tutorial nodes en de verbindingen ertussen edges. Let op dat edges en nodes vaak andere namen hebben die doorgaans exact hetzelfde betekenen:

node == vertex == point
edge == arc == link

De uitgangsgraaf is ongericht. Dat wil zeggen: je edges hebben geen oriëntatie en zijn bidirectioneel. Bijvoorbeeld: A<--->B == B<--->A.
Daarentegen kan de graaf die je maakt om het kortste pad te bepalen om elk pad te lopen een gerichte graaf zijn, waar de volgorde en richting van edges ertoe doen. Bijvoorbeeld: A--->B != B--->A.

De graaf is ook een edge-gewogen graaf waarbij de afstand (in mijl) tussen elk paar aangrenzende nodes het gewicht van een edge weergeeft. Dit wordt beheerd als een edge-attribuut met de naam "distance".

Graad verwijst naar het aantal edges dat aan een node raakt. Nodes noem je oneven graad wanneer dit aantal oneven is en even graad wanneer het even is.

De oplossing van dit CPP wordt een Euleriaanse tour: een graaf waarin je een cyclus kunt maken die elk edge precies één keer doorloopt en weer terugkomt bij de startnode (zonder terug te lopen). Een Euler Tour staat ook bekend onder meerdere namen:

Eulerian tour == Eulerian circuit == Eulerian cycle

Een matching is een subset van edges waarin geen enkele node meer dan één keer voorkomt. Een minimumgewicht-matching vindt de matching met de laagst mogelijke som van edge-gewichten.

NetworkX: manipulatie en analyse van grafen

NetworkX is het populairste Python-pakket voor het manipuleren en analyseren van grafen. Meerdere pakketten bieden dezelfde basisfunctionaliteit, met name igraph, dat ook bindings heeft voor R en C++. Ik vond echter dat NetworkX de sterkste graafalgoritmen had die ik nodig had om het CPP op te lossen.

Pakketten installeren

Als je wel eens data-analyse in Python hebt gedaan of de Anaconda-distributie gebruikt, is de kans groot dat je pandas en matplotlib al hebt. Maar misschien heb je networkx nog niet. Dit zouden de enige afhankelijkheden buiten de Python Standard Library moeten zijn die je nodig hebt voor deze tutorial. Ze zijn eenvoudig te installeren met pip:

pip install pandas
pip install networkx
pip install matplotlib

Dit zijn alle pakketten die je voorlopig nodig hebt. imageio en numpy worden helemaal aan het eind geïmporteerd om de GIF-animatie van de CPP-oplossing te maken. De animatie is in dit bericht ingesloten, dus deze pakketten zijn optioneel.

import itertools
import copy
import networkx as nx
import pandas as pd
import matplotlib.pyplot as plt

Data laden

Edgelijst

De edgelijst is een eenvoudige datastructuur waarmee je de graaf maakt. Elke rij stelt één edge in de graaf voor met enkele edge-attributen.

  • node1 & node2: namen van de verbonden nodes.
  • trail: edge-attribuut dat de afgekorte naam van het pad voor elke edge aangeeft. Bijvoorbeeld: rs = red square
  • distance: edge-attribuut dat de padlengte in mijl aangeeft.
  • color: padkleur die voor het plotten wordt gebruikt.
  • estimate: edge-attribuut dat aangeeft of de edge-afstand is geschat door op de trailmap te turen (1=ja, 0=nee), omdat sommige afstanden niet zijn opgegeven. Dit is puur ter referentie; het wordt niet gebruikt voor de analyse.
# Grab edge list data hosted on Gist
edgelist = pd.read_csv('https://gist.githubusercontent.com/brooksandrew/e570c38bcc72a8d102422f2af836513b/raw/89c76b2563dbc0e88384719a35cba0dfc04cd522/edgelist_sleeping_giant.csv')

# Preview edgelist
edgelist.head(10)
 node1node2traildistancecolorestimate
0rs_end_northv_rsrs0.30red0
1v_rsb_rsrs0.21red0
2b_rsg_rsrs0.11red0
3g_rsw_rsrs0.18red0
4w_rso_rsrs0.21red0
5o_rsy_rsrs0.12red0
6y_rsrs_end_southrs0.39red0
7rc_end_northv_rcrc0.70red0
8v_rcb_rcrc0.04red0
9b_rcg_rcrc0.15red0

Nodelijst

Nodelijsten zijn meestal optioneel in networkx en andere graafbibliotheken wanneer edgelijsten worden aangeleverd, omdat de nodenamen in de eerste twee kolommen van de edgelijst staan. In dit geval willen we echter enkele node-attributen toevoegen: X- en Y-coördinaten van de nodes (kruisingen van paden), zodat je je graaf met dezelfde lay-out als de trailmap kunt plotten.

Ik heb een middag besteed aan het handmatig annoteren door over de afbeelding te tekenen met GIMP:

  • id: naam van de node die overeenkomt met node1 en node2 in de edgelijst.
  • X: horizontale positie/coördinaat van de node ten opzichte van de linkerbovenhoek.
  • Y verticale positie/coördinaat van de node ten opzichte van de linkerbovenhoek.

Opmerking over het genereren van de node- en edgelijsten

Het maken van de nodenamen vergde ook wat handwerk. Elke node vertegenwoordigt een kruising van twee of meer paden. Waar mogelijk is de node genoemd naar trail1_trail2, waarbij trail1 voorafgaat aan trail2 in alfabetische volgorde.

Het werd wat lastiger wanneer dezelfde paden elkaar meer dan eens kruisten. Bijvoorbeeld het Orange- en White-pad. In die gevallen heb ik _2 of _3 aan de nodenaam toegevoegd. Zo heb je twee verschillende nodenamen voor de twee verschillende kruisingen van Orange en White: o_w en o_w_2.

Dit vergde veel trial-and-error en het vergelijken van de met X,Y-coördinaten gegenereerde plots met de echte trailmap.

# Grab node list data hosted on Gist
nodelist = pd.read_csv('https://gist.githubusercontent.com/brooksandrew/f989e10af17fb4c85b11409fea47895b/raw/a3a8da0fa5b094f1ca9d82e1642b384889ae16e8/nodelist_sleeping_giant.csv')

# Preview nodelist
nodelist.head(5)
 idXY
0b_bv1486732
1b_bw7161357
2b_end_east31641111
3b_end_west1411938
4b_g1725771

Graaf maken

Nu gebruik je de edgelijst en de nodelijst om een graafobject in networkx te maken.

# Create empty graph
g = nx.Graph()

Loop door de rijen van de edgelijst en voeg elke edge en de bijbehorende attributen toe aan graaf g.

# Add edges and edge attributes
for i, elrow in edgelist.iterrows():
    g.add_edge(elrow[0], elrow[1], attr_dict=elrow[2:].to_dict())

Om te laten zien wat hier gebeurt, printen we de waarden van de laatste rij in de edgelijst die aan graaf g is toegevoegd:

# Edge list example
print(elrow[0]) # node1
print(elrow[1]) # node2
print(elrow[2:].to_dict()) # edge attribute dict

o_gy2
y_gy2
{'color': 'yellowgreen', 'estimate': 0, 'trail': 'gy2', 'distance': 0.12}

Op dezelfde manier loop je door de rijen in de nodelijst en voeg je deze node-attributen toe.

# Add node attributes
for i, nlrow in nodelist.iterrows():
    g.node[nlrow['id']] = nlrow[1:].to_dict()

Hier is een voorbeeld uit de laatste rij van de nodelijst:

# Node list example
print(nlrow)

id    y_rt
X      977
Y     1666
Name: 76, dtype: object

Graaf inspecteren

Edges

Je graafedges worden weergegeven als een lijst met tuples van lengte 3. De eerste twee elementen zijn de nodenamen die door de edge verbonden zijn. De derde is de dictionary met edge-attributen.

# Preview first 5 edges
g.edges(data=True)[0:5]

[('rs_end_south',
  'y_rs',
  {'color': 'red', 'distance': 0.39, 'estimate': 0, 'trail': 'rs'}),
 ('w_gy2',
  'park_east',
  {'color': 'gray', 'distance': 0.12, 'estimate': 0, 'trail': 'w'}),
 ('w_gy2',
  'g_gy2',
  {'color': 'yellowgreen', 'distance': 0.05, 'estimate': 0, 'trail': 'gy2'}),
 ('w_gy2',
  'b_w',
  {'color': 'gray', 'distance': 0.42, 'estimate': 0, 'trail': 'w'}),
 ('w_gy2',
  'b_gy2',
  {'color': 'yellowgreen', 'distance': 0.03, 'estimate': 1, 'trail': 'gy2'})]

Nodes

Evenzo worden je nodes weergegeven als een lijst met tuples van lengte 2. Het eerste element is de node-ID, gevolgd door de dictionary met node-attributen.

# Preview first 10 nodes
g.nodes(data=True)[0:10]

[('rs_end_south', {'X': 1865, 'Y': 1598}),
 ('w_gy2', {'X': 2000, 'Y': 954}),
 ('rd_end_south_dupe', {'X': 273, 'Y': 1869}),
 ('w_gy1', {'X': 1184, 'Y': 1445}),
 ('g_rt', {'X': 908, 'Y': 1378}),
 ('v_rd', {'X': 258, 'Y': 1684}),
 ('g_rs', {'X': 1676, 'Y': 775}),
 ('rc_end_north', {'X': 867, 'Y': 618}),
 ('v_end_east', {'X': 2131, 'Y': 921}),
 ('rh_end_south', {'X': 721, 'Y': 1925})]

Samenvattende statistieken

Print eerst wat samenvattende statistieken voordat je de graaf visualiseert.

print('# of edges: {}'.format(g.number_of_edges()))
print('# of nodes: {}'.format(g.number_of_nodes()))

# of edges: 123
# of nodes: 77

Graaf visualiseren

Kleuren en lay-out manipuleren

Posities: Eerst moet je de nodeposities uit de graaf omzetten in een dictionary. Zo kun je de graaf opnieuw maken met dezelfde lay-out als de echte trailmap. Y wordt negatief gemaakt om de oorsprong van de Y-as van linksboven naar linksonder te verplaatsen.

# Define node positions data structure (dict) for plotting
node_positions = {node[0]: (node[1]['X'], -node[1]['Y']) for node in g.nodes(data=True)}

# Preview of node_positions with a bit of hack (there is no head/slice method for dictionaries).
dict(list(node_positions.items())[0:5])

{'b_rd': (268, -1744),
 'g_rt': (908, -1378),
 'o_gy1': (1130, -1297),
 'rh_end_tt_2': (550, -1608),
 'rs_end_south': (1865, -1598)}

Kleuren: Nu zet je de edgekleuren uit de graaf om naar een eenvoudige lijst, zodat je de paden op kleur kunt visualiseren.

# Define data structure (list) of edge colors for plotting
edge_colors = [e[2]['color'] for e in g.edges(data=True)]

# Preview first 10
edge_colors[0:10]

['red',
 'gray',
 'yellowgreen',
 'gray',
 'yellowgreen',
 'blue',
 'black',
 'yellowgreen',
 'gray',
 'gray']

Plot

Nu kun je een nette plot maken die mooi uitlijnt met de trailmap van Sleeping Giant:

plt.figure(figsize=(8, 6))
nx.draw(g, pos=node_positions, edge_color=edge_colors, node_size=10, node_color='black')
plt.title('Graph Representation of Sleeping Giant Trail Map', size=15)
plt.show()

Graafoptimalisatie in Python

Deze graafweergave vangt natuurlijk niet alle bochten en kronkels van de paden, maar geen zorgen: die worden nauwkeurig vastgelegd in het edge-attribuut distance, dat voor de berekeningen wordt gebruikt. De visual laat wel de afstand tussen nodes (kruisingen) zien in vogelvlucht, wat een aardige benadering lijkt.

Overzicht van het CPP-algoritme

Oké, nu je de termen hebt gedefinieerd en de graaf hebt gemaakt: hoe vind je de kortste route erdoorheen?

Het Chinese postbodeprobleem oplossen is conceptueel vrij eenvoudig:

  1. Vind alle nodes met oneven graad (heel makkelijk).
    (Vind alle kruisingen waar het aantal paden dat de kruising raakt oneven is)

  2. Voeg edges toe aan de graaf zodat alle nodes met oneven graad even worden. Deze toegevoegde edges moeten duplicaten zijn van de oorspronkelijke graaf (we gaan voor dit probleem uit van geen bushwhacking). De set toegevoegde edges moet tot de minimale mogelijke afstand optellen (moeilijk... in feite np-hard).
    (Eenvoudiger gezegd: minimaliseer hoeveel je moet teruglopen op een route die elk pad aandoet)

  3. Vind, gegeven een startpunt, de Euleriaanse tour over de aangevulde dataset (redelijk makkelijk).
    (Zodra we weten over welke paden we gaan teruglopen, bereken dan de route van begin tot eind)

Aannames en vereenvoudigingen

Hoewel een korter en preciezer pad te genereren is door onderstaande aannames te versoepelen, voegt dat complexiteit toe buiten de scope van deze tutorial, die focust op de CPP.

Aannames 1: Alleen verplichte paden

Zoals je op de trailmap hierboven ziet, zijn er wegen langs de randen van het park die paden kunnen verbinden, met name de rode paden. Er zijn ook paden (Horseshoe en ongemarkeerde markeringen) die niet vereist zijn volgens het Giantmaster-log, maar die kunnen helpen om veel teruglopen te voorkomen. Het opnemen van optionele paden is een bekende variant van het CPP, het Rural Postman-probleem. We negeren optionele paden in deze tutorial en focussen alleen op verplichte paden.

Aannames 2: Bergop == bergaf

Het CPP gaat ervan uit dat de kosten van een pad gelijk zijn aan de afstand, ongeacht in welke richting je het loopt. Sommige paden zijn echter behoorlijk heuvelachtig en kosten meer energie omhoog dan omlaag. Een maat die zowel afstand als hoogteverschil combineert over een gerichte graaf, kan worden opgenomen in een uitbreiding van het CPP, het Windy Postman-probleem.

Aannames 3: Geen parallelle edges (paden)

Hoewel mogelijk, maakt het opnemen van parallelle edges (meerdere paden die dezelfde twee nodes verbinden) de berekening complexer. Gelukkig gebeurt dat hier maar twee keer (Blue <=> Red Diamond en Blue <=> Tower Trail). Dit is opgelost met een kleine hack in de edgelijst: dubbele nodes zijn opgenomen met een _dupe-suffix om elk pad te vangen en toch uniciteit in de edges te behouden. De CPP-implementatie in het postman_problems-pakket dat ik schreef, behandelt parallelle edges robuust op een elegantere manier als je het CPP op je eigen graaf met veel parallelle edges wilt oplossen.

CPP stap 1: nodes met oneven graad vinden

Dit is een vrij rechttoe rechtaan telberekening. Je ziet dat 36 van de 76 nodes een oneven graad hebben. Dit zijn vooral doodlopende paden (graad 1) en kruisingen van 3 paden. Er zijn een handvol nodes met graad 5.

# Calculate list of nodes with odd degree
nodes_odd_degree = [v for v, d in g.degree_iter() if d % 2 == 1]

# Preview
nodes_odd_degree[0:5]

['rs_end_south', 'rc_end_north', 'v_end_east', 'rh_end_south', 'b_end_east']

# Counts
print('Number of nodes of odd degree: {}'.format(len(nodes_odd_degree)))
print('Number of total nodes: {}'.format(len(g.nodes())))

Number of nodes of odd degree: 36
Number of total nodes: 77

CPP stap 2: paren met minimale afstand vinden

Dit is echt de kern van het probleem. Je splitst het op in 5 delen:

  1. Bereken alle mogelijke paren van nodes met oneven graad.
  2. Bereken het kortste pad tussen elk nodepaar uit 1..
  3. Maak een complete graaf die elk nodepaar uit 1. verbindt met kortste-pad-afstanden als attributen uit 2.
  4. Bereken een minimumgewicht-matching van de graaf uit 3.
    (Dit komt neer op bepalen hoe je de oneven nodes paart zodat de som van de afstanden tussen de paren zo klein mogelijk is).
  5. Vul de oorspronkelijke graaf aan met de kortste paden tussen de nodeparen uit 4.

Stap 2.1: nodeparen berekenen

Je gebruikt de functie itertools combination om alle mogelijke paren van nodes met oneven graad te berekenen. Je graaf is ongericht, dus de volgorde maakt niet uit: bijvoorbeeld (a,b) == (b,a).

# Compute all pairs of odd nodes. in a list of tuples
odd_node_pairs = list(itertools.combinations(nodes_odd_degree, 2))

# Preview pairs of odd degree nodes
odd_node_pairs[0:10]

[('rs_end_south', 'rc_end_north'),
 ('rs_end_south', 'v_end_east'),
 ('rs_end_south', 'rh_end_south'),
 ('rs_end_south', 'b_end_east'),
 ('rs_end_south', 'b_bv'),
 ('rs_end_south', 'rt_end_south'),
 ('rs_end_south', 'o_rt'),
 ('rs_end_south', 'y_rt'),
 ('rs_end_south', 'g_gy2'),
 ('rs_end_south', 'b_tt_3')]

# Counts
print('Number of pairs: {}'.format(len(odd_node_pairs)))

Number of pairs: 630

Laten we bevestigen dat dit aantal paren klopt met de combinatoriek hieronder. Gelukkig heb je hier maar 630 paren om je zorgen over te maken. Je rekentijd om dit CPP-voorbeeld op te lossen is verwaarloosbaar (een paar seconden).

Maar als je in plaats daarvan 3.600 oneven nodeparen had, zou je ~6,5 miljoen paren moeten optimaliseren. Dat is ~10.000x meer output bij een 100x grotere input.

\begin{equation*} \#\;of\;pairs = n\;choose\;r = {n \choose r} = \frac{n!}{r!(n-r)!} = \frac{36!}{2! (36-2)!} = 630 \end{equation*}

Stap 2.2: kortste paden tussen nodeparen berekenen

Dit is de eerste stap die wat echte berekening vergt. Gelukkig heeft networkx een handige implementatie van Dijkstra's algoritme om het kortste pad tussen twee nodes te berekenen. Je past deze functie toe op elk paar (alle 630) uit odd_node_pairs.

def get_shortest_paths_distances(graph, pairs, edge_weight_name):
    """Compute shortest distance between each pair of nodes in a graph.  Return a dictionary keyed on node pairs (tuples)."""
    distances = {}
    for pair in pairs:
        distances[pair] = nx.dijkstra_path_length(graph, pair[0], pair[1], weight=edge_weight_name)
    return distances

# Compute shortest paths.  Return a dictionary with node pairs keys and a single value equal to shortest path distance.
odd_node_pairs_shortest_paths = get_shortest_paths_distances(g, odd_node_pairs, 'distance')

# Preview with a bit of hack (there is no head/slice method for dictionaries).
dict(list(odd_node_pairs_shortest_paths.items())[0:10])

{('b_bv', 'y_gy1'): 1.22,
 ('b_bw', 'rc_end_south'): 1.35,
 ('b_end_east', 'b_bw'): 3.0400000000000005,
 ('b_end_east', 'rd_end_north'): 3.83,
 ('g_gy1', 'nature_end_west'): 0.9900000000000001,
 ('o_rt', 'y_gy1'): 0.53,
 ('rc_end_north', 'rd_end_south'): 2.21,
 ('rc_end_north', 'rs_end_north'): 1.79,
 ('rs_end_north', 'o_tt'): 2.0999999999999996,
 ('w_bw', 'rd_end_north'): 1.02}

Stap 2.3: complete graaf maken

Een complete graaf is simpelweg een graaf waarin elke node met elke andere node is verbonden door een unieke edge.

Hier is een basisvoorbeeld van Wikipedia van een complete graaf met 7 nodes en 21 (7 over 2) edges:

networkx-tutorial python

De graaf die je hieronder maakt heeft 36 nodes en 630 edges met hun bijbehorende edgegewicht (afstand).

create_complete_graph is gedefinieerd om deze te berekenen. De parameter flip_weights wordt gebruikt om distance te transformeren naar het attribuut weight, waarbij kleinere getallen grote afstanden weergeven en grote getallen korte afstanden. Dit klinkt wat contra-intuïtief, maar is nodig voor stap 2.4, waar je de minimumgewicht-matching op de complete graaf berekent.

Idealiter zou je de minimumgewicht-matching direct berekenen, maar NetworkX implementeert alleen een functie max_weight_matching die het edgegewicht maximaliseert in plaats van minimaliseert. We hacken dit door het attribuut distance te negateren (vermenigvuldigen met -1) om weight te krijgen. Zo blijven volgorde en schaal naar afstand behouden, maar dan omgekeerd.

def create_complete_graph(pair_weights, flip_weights=True):
    """
    Create a completely connected graph using a list of vertex pairs and the shortest path distances between them
    Parameters:
        pair_weights: list[tuple] from the output of get_shortest_paths_distances
        flip_weights: Boolean. Should we negate the edge attribute in pair_weights?
    """
    g = nx.Graph()
    for k, v in pair_weights.items():
        wt_i = - v if flip_weights else v
        g.add_edge(k[0], k[1], attr_dict={'distance': v, 'weight': wt_i})
    return g

# Generate the complete graph
g_odd_complete = create_complete_graph(odd_node_pairs_shortest_paths, flip_weights=True)

# Counts
print('Number of nodes: {}'.format(len(g_odd_complete.nodes())))
print('Number of edges: {}'.format(len(g_odd_complete.edges())))

Number of nodes: 36
Number of edges: 630

Als visuele ondersteuning is de volledig verbonden graaf van nodeparen met oneven graad hieronder geplot. Let op dat je de X- en Y-coördinaten van elke node behoudt, maar de edges vertegenwoordigen niet per se echte paden. Bijvoorbeeld: twee nodes kunnen in deze graaf door één edge verbonden zijn, maar het kortste pad ertussen kan 5 sprongen via nodes met even graad bevatten (hier niet getoond).

# Plot the complete graph of odd-degree nodes
plt.figure(figsize=(8, 6))
pos_random = nx.random_layout(g_odd_complete)
nx.draw_networkx_nodes(g_odd_complete, node_positions, node_size=20, node_color="red")
nx.draw_networkx_edges(g_odd_complete, node_positions, alpha=0.1)
plt.axis('off')
plt.title('Complete Graph of Odd-degree Nodes')
plt.show()

graaf van nodes met oneven graad

Stap 2.4: minimumgewicht-matching berekenen

Dit is de meest complexe stap in de CPP. Je moet de nodeparen met oneven graad vinden waarvan de gecombineerde som (van de afstand ertussen) zo klein mogelijk is. Voor jouw probleem komt dit neer op het selecteren van de optimale 18 edges (36 nodes met oneven graad / 2) uit de wirwar van een graaf die in 2.3 is gegenereerd.

Zowel de implementatie als de intuïtie achter deze optimalisatie vallen buiten de scope van deze tutorial... denk aan 800+ regels code en een corpus aan academische literatuur dat verder gaat dan hier.

Toch een korte zijstap voor de geïnteresseerde lezer:

Groot dank aan Joris van Rantwijk voor het schrijven van de oorspronkelijke implementatie op zijn blog in 2008. Ik belandde op een vergelijkbare manier bij het probleem met dezelfde intentie als Joris. Uit Joris’ post uit 2008:

Since I did not find any Perl implementations of maximum weighted matching, I lightly decided to write some code myself. It turned out that I had underestimated the problem, but by the time I realized my mistake, I was so obsessed with the problem that I refused to give up.

Ik heb echter wél opgegeven. Gelukkig deed Joris dat niet.

Deze Maximum Weight Matching is sindsdien opgenomen in en onderhouden binnen het NetworkX-pakket. Nogmaals dank aan de 10+ bijdragers op GitHub die deze forse codebase hebben onderhouden.

Dit is een zware en intensieve berekening. De eerste doorbraak in 1965 bewees dat het Maximum Matching-probleem in polynomiale tijd kon worden opgelost. Het werd gepubliceerd door Jack Edmonds met misschien wel een van de mooiste titels voor een academisch artikel ooit: "Paths, trees, and flowers" [1]. Een corpus aan literatuur heeft hierop voortgebouwd en de optimalisatieprocedure verbeterd. De code die is geïmplementeerd in de NetworkX-functie max_weight_matching is gebaseerd op Galil, Zvi (1986) [2] en gebruikt een algoritme met O(n3)-tijd.

# Compute min weight matching.
# Note: max_weight_matching uses the 'weight' attribute by default as the attribute to maximize.
odd_matching_dupes = nx.algorithms.max_weight_matching(g_odd_complete, True)

print('Number of edges in matching: {}'.format(len(odd_matching_dupes)))

Number of edges in matching: 36

De matching-output (odd_matching_dupes) is een dictionary. Hoewel er 36 edges in deze matching zitten, wil je er maar 18. Elk edge-paar komt twee keer voor (één keer met node 1 als sleutel en een tweede keer met node 2 als sleutel van de dictionary).

# Preview of matching with dupes
odd_matching_dupes

{'b_bv': 'v_bv',
 'b_bw': 'rh_end_tt_1',
 'b_end_east': 'g_gy2',
 'b_end_west': 'rd_end_south',
 'b_tt_3': 'rt_end_north',
 'b_v': 'v_end_west',
 'g_gy1': 'rc_end_north',
 'g_gy2': 'b_end_east',
 'g_w': 'w_bw',
 'nature_end_west': 'o_y_tt_end_west',
 'o_rt': 'o_w_1',
 'o_tt': 'rh_end_tt_2',
 'o_w_1': 'o_rt',
 'o_y_tt_end_west': 'nature_end_west',
 'rc_end_north': 'g_gy1',
 'rc_end_south': 'y_gy1',
 'rd_end_north': 'rh_end_north',
 'rd_end_south': 'b_end_west',
 'rh_end_north': 'rd_end_north',
 'rh_end_south': 'y_rh',
 'rh_end_tt_1': 'b_bw',
 'rh_end_tt_2': 'o_tt',
 'rh_end_tt_3': 'rh_end_tt_4',
 'rh_end_tt_4': 'rh_end_tt_3',
 'rs_end_north': 'v_end_east',
 'rs_end_south': 'y_gy2',
 'rt_end_north': 'b_tt_3',
 'rt_end_south': 'y_rt',
 'v_bv': 'b_bv',
 'v_end_east': 'rs_end_north',
 'v_end_west': 'b_v',
 'w_bw': 'g_w',
 'y_gy1': 'rc_end_south',
 'y_gy2': 'rs_end_south',
 'y_rh': 'rh_end_south',
 'y_rt': 'rt_end_south'}

Je zet deze dictionary om naar een lijst met tuples, omdat je een ongerichte graaf hebt en volgorde niet uitmaakt. Duplicaten verwijderen levert de unieke 18 edge-paren op die samen tot de kleinst mogelijke afstand optellen.

# Convert matching to list of deduped tuples
odd_matching = list(pd.unique([tuple(sorted([k, v])) for k, v in odd_matching_dupes.items()]))

# Counts
print('Number of edges in matching (deduped): {}'.format(len(odd_matching)))

Number of edges in matching (deduped): 18

# Preview of deduped matching
odd_matching

[('rs_end_south', 'y_gy2'),
 ('b_end_west', 'rd_end_south'),
 ('b_bv', 'v_bv'),
 ('rh_end_tt_3', 'rh_end_tt_4'),
 ('b_bw', 'rh_end_tt_1'),
 ('o_tt', 'rh_end_tt_2'),
 ('g_w', 'w_bw'),
 ('b_end_east', 'g_gy2'),
 ('nature_end_west', 'o_y_tt_end_west'),
 ('g_gy1', 'rc_end_north'),
 ('o_rt', 'o_w_1'),
 ('rs_end_north', 'v_end_east'),
 ('rc_end_south', 'y_gy1'),
 ('rh_end_south', 'y_rh'),
 ('rt_end_south', 'y_rt'),
 ('b_tt_3', 'rt_end_north'),
 ('rd_end_north', 'rh_end_north'),
 ('b_v', 'v_end_west')]

Laten we deze paren visualiseren op de complete graaf die eerder in stap 2.3 is geplot. Zoals eerder geven de nodeposities hier de echte graaf (trailmap) weer, maar de getoonde edgedistances (blauwe lijnen) zijn in vogelvlucht. Het daadwerkelijke kortste pad van de ene naar de andere node kan uit meerdere edges bestaan die kronkelen en aanzienlijk langer zijn.

plt.figure(figsize=(8, 6))

# Plot the complete graph of odd-degree nodes
nx.draw(g_odd_complete, pos=node_positions, node_size=20, alpha=0.05)

# Create a new graph to overlay on g_odd_complete with just the edges from the min weight matching
g_odd_complete_min_edges = nx.Graph(odd_matching)
nx.draw(g_odd_complete_min_edges, pos=node_positions, node_size=20, edge_color='blue', node_color='red')

plt.title('Min Weight Matching on Complete Graph')
plt.show()

min weight matching op complete graaf

Om te laten zien hoe dit past bij de oorspronkelijke graaf, plot je dezelfde minimumgewicht-paren (blauwe lijnen), maar dan over de trailmap (vervaagd) in plaats van over de complete graaf. Merk opnieuw op dat de blauwe lijnen de bushwhacking-route zijn (edges als de vogel vliegt, geen echte paden). Je hebt nog wat werk te doen om de edges te vinden die het kortste pad tussen elk paar vormen in stap 3.

plt.figure(figsize=(8, 6))

# Plot the original trail map graph
nx.draw(g, pos=node_positions, node_size=20, alpha=0.1, node_color='black')

# Plot graph to overlay with just the edges from the min weight matching
nx.draw(g_odd_complete_min_edges, pos=node_positions, node_size=20, alpha=1, node_color='red', edge_color='blue')

plt.title('Min Weight Matching on Orginal Graph')
plt.show()

min weight matching op oorspronkelijke graaf

Stap 2.5: de oorspronkelijke graaf aanvullen

Nu vul je de oorspronkelijke graaf aan met de edges uit de matching van 2.4. Hieronder staat een eenvoudige functie om dit te doen, die ook noteert dat deze nieuwe edges uit de aangevulde graaf komen. Dat heb je nodig in 3. wanneer je de Euleriaanse cyclus daadwerkelijk door de graaf maakt.

def add_augmenting_path_to_graph(graph, min_weight_pairs):
    """
    Add the min weight matching edges to the original graph
    Parameters:
        graph: NetworkX graph (original graph from trailmap)
        min_weight_pairs: list[tuples] of node pairs from min weight matching
    Returns:
        augmented NetworkX graph
    """

    # We need to make the augmented graph a MultiGraph so we can add parallel edges
    graph_aug = nx.MultiGraph(graph.copy())
    for pair in min_weight_pairs:
        graph_aug.add_edge(pair[0],
                           pair[1],
                           attr_dict={'distance': nx.dijkstra_path_length(graph, pair[0], pair[1]),
                                      'trail': 'augmented'}
                          )
    return graph_aug

Bevestigen we even dat je aangevulde graaf het verwachte aantal (18) edges toevoegt:

# Create augmented graph: add the min weight matching edges to g
g_aug = add_augmenting_path_to_graph(g, odd_matching)

# Counts
print('Number of edges in original graph: {}'.format(len(g.edges())))
print('Number of edges in augmented graph: {}'.format(len(g_aug.edges())))

Number of edges in original graph: 123
Number of edges in augmented graph: 141

Laten we ook bevestigen dat nu elke node een even graad heeft:

pd.value_counts(g_aug.degree())

4    54
2    18
6     5
dtype: int64

CPP stap 3: Euleriaanse cyclus berekenen

Nu je een graaf met even graad hebt, is het zware optimalisatiewerk voorbij. Zoals Euler in 1736 beroemd stelde met het probleem van de zeven bruggen van Königsberg, bestaat er een pad dat elke edge precies één keer bezoekt als alle nodes even graad hebben. Carl Hierholzer bewees dit formeel in de jaren 1870.

Er zijn veel Euleriaanse cycli met dezelfde afstand die je kunt construeren. Je komt 90% in de buurt met de NetworkX-functie eulerian_circuit. Er zijn echter enkele beperkingen.

Beperkingen die je gaat fixen:

  1. De aangevulde graaf kan (en zal waarschijnlijk) edges bevatten die niet in de oorspronkelijke graaf voorkwamen. Om de cyclus te krijgen (zonder bushwhacking), moet je deze aangevulde edges opsplitsen in het kortste pad via de edges die wel bestaan.

  2. eulerian_circuit geeft alleen de volgorde terug waarin we elke node aandoen. Het geeft niet de attributes terug van de edges die nodig zijn om de cyclus te voltooien. Dat is nodig omdat je moet bijhouden welke edges al bewandeld zijn wanneer er meerdere edges tussen twee nodes bestaan.

Beperkingen die je niet gaat fixen:

  1. Om je benen wat werk te besparen, zou je de aanname van de Euleriaanse cyclus kunnen versoepelen dat je op dezelfde node start en eindigt. Een Euleriaans pad (de algemene vorm van de Euleriaanse cyclus) kan ook worden gevonden als er precies twee nodes met oneven graad zijn. Dat zou je wat teruglopen schelen... ervan uitgaande dat je een lift terug kunt regelen vanaf de andere kant van het park. Op het moment van schrijven biedt NetworkX echter geen Euler Path-algoritme. De eulerian_circuit-code is niet al te ingewikkeld en zou hiervoor aangepast kunnen worden, maar we houden het hier simpel.

Naïeve cyclus

Toch beginnen we met de simpele maar onvolledige oplossing:

naive_euler_circuit = list(nx.eulerian_circuit(g_aug, source='b_end_east'))

Zoals verwacht is de lengte van de naïeve Euleriaanse cyclus gelijk aan het aantal edges in de aangevulde graaf.

print('Length of eulerian circuit: {}'.format(len(naive_euler_circuit)))

Length of eulerian circuit: 141

De output is gewoon een lijst met tuples die nodeparen voorstellen. Merk op dat de eerste node van elk paar dezelfde is als de tweede node van het voorgaande paar.

# Preview naive Euler circuit
naive_euler_circuit[0:10]

[('b_end_east', 'g_gy2'),
 ('g_gy2', 'b_g'),
 ('b_g', 'b_w'),
 ('b_w', 'b_gy2'),
 ('b_gy2', 'w_gy2'),
 ('w_gy2', 'b_w'),
 ('b_w', 'w_rs'),
 ('w_rs', 'g_rs'),
 ('g_rs', 'b_g'),
 ('b_g', 'b_rs')]

Correcte cyclus

Nu definiëren we een functie die de oorspronkelijke graaf gebruikt om je te vertellen welke paden je moet nemen van node A naar node B. Hoewel de code uitvoerig is, is de logica eigenlijk heel simpel. Je transformeert de naïeve cyclus, die edges bevatte die niet in de oorspronkelijke graaf bestonden, naar een Euleriaanse cyclus die alleen edges gebruikt die in de oorspronkelijke graaf bestaan.

Je loopt elk edge in de naïeve Euleriaanse cyclus (naive_euler_circuit) langs. Waar je een edge tegenkomt die niet in de oorspronkelijke graaf bestaat, vervang je die door de reeks edges die samen het kortste pad vormen tussen de nodes met behulp van de oorspronkelijke graaf.

def create_eulerian_circuit(graph_augmented, graph_original, starting_node=None):
    """Create the eulerian path using only edges from the original graph."""
    euler_circuit = []
    naive_circuit = list(nx.eulerian_circuit(graph_augmented, source=starting_node))

    for edge in naive_circuit:
        edge_data = graph_augmented.get_edge_data(edge[0], edge[1])    

        if edge_data[0]['trail'] != 'augmented':
            # If `edge` exists in original graph, grab the edge attributes and add to eulerian circuit.
            edge_att = graph_original[edge[0]][edge[1]]
            euler_circuit.append((edge[0], edge[1], edge_att))
        else:
            aug_path = nx.shortest_path(graph_original, edge[0], edge[1], weight='distance')
            aug_path_pairs = list(zip(aug_path[:-1], aug_path[1:]))

            print('Filling in edges for augmented edge: {}'.format(edge))
            print('Augmenting path: {}'.format(' => '.join(aug_path)))
            print('Augmenting path pairs: {}\n'.format(aug_path_pairs))

            # If `edge` does not exist in original graph, find the shortest path between its nodes and
            #  add the edge attributes for each link in the shortest path.
            for edge_aug in aug_path_pairs:
                edge_aug_att = graph_original[edge_aug[0]][edge_aug[1]]
                euler_circuit.append((edge_aug[0], edge_aug[1], edge_aug_att))

    return euler_circuit

Je omzeilt beperking 3 een beetje door de Euleriaanse cyclus te starten aan de oostkant van het park op het Blue-pad (node "b_end_east"). In de praktijk kun je simpelweg de laatste aanwijzing overslaan die hierop terugloopt.

Er zijn uitgebreide printstatements toegevoegd om te laten zien wat er gebeurt wanneer je niet-bestaande edges uit de aangevulde graaf vervangt door het kortste pad met edges die wél bestaan.

# Create the Eulerian circuit
euler_circuit = create_eulerian_circuit(g_aug, g, 'b_end_east')

Filling in edges for augmented edge: ('b_end_east', 'g_gy2')
Augmenting path: b_end_east => b_y => b_o => b_gy2 => w_gy2 => g_gy2
Augmenting path pairs: [('b_end_east', 'b_y'), ('b_y', 'b_o'), ('b_o', 'b_gy2'), ('b_gy2', 'w_gy2'), ('w_gy2', 'g_gy2')]

Filling in edges for augmented edge: ('b_bw', 'rh_end_tt_1')
Augmenting path: b_bw => b_tt_1 => rh_end_tt_1
Augmenting path pairs: [('b_bw', 'b_tt_1'), ('b_tt_1', 'rh_end_tt_1')]

Filling in edges for augmented edge: ('b_tt_3', 'rt_end_north')
Augmenting path: b_tt_3 => b_tt_2 => tt_rt => v_rt => rt_end_north
Augmenting path pairs: [('b_tt_3', 'b_tt_2'), ('b_tt_2', 'tt_rt'), ('tt_rt', 'v_rt'), ('v_rt', 'rt_end_north')]

Filling in edges for augmented edge: ('rc_end_north', 'g_gy1')
Augmenting path: rc_end_north => v_rc => b_rc => g_rc => g_gy1
Augmenting path pairs: [('rc_end_north', 'v_rc'), ('v_rc', 'b_rc'), ('b_rc', 'g_rc'), ('g_rc', 'g_gy1')]

Filling in edges for augmented edge: ('y_gy1', 'rc_end_south')
Augmenting path: y_gy1 => y_rc => rc_end_south
Augmenting path pairs: [('y_gy1', 'y_rc'), ('y_rc', 'rc_end_south')]

Filling in edges for augmented edge: ('b_end_west', 'rd_end_south')
Augmenting path: b_end_west => b_v => rd_end_south
Augmenting path pairs: [('b_end_west', 'b_v'), ('b_v', 'rd_end_south')]

Filling in edges for augmented edge: ('rh_end_north', 'rd_end_north')
Augmenting path: rh_end_north => v_rh => v_rd => rd_end_north
Augmenting path pairs: [('rh_end_north', 'v_rh'), ('v_rh', 'v_rd'), ('v_rd', 'rd_end_north')]

Filling in edges for augmented edge: ('v_end_east', 'rs_end_north')
Augmenting path: v_end_east => v_rs => rs_end_north
Augmenting path pairs: [('v_end_east', 'v_rs'), ('v_rs', 'rs_end_north')]

Filling in edges for augmented edge: ('y_gy2', 'rs_end_south')
Augmenting path: y_gy2 => y_rs => rs_end_south
Augmenting path pairs: [('y_gy2', 'y_rs'), ('y_rs', 'rs_end_south')]

Je ziet dat de lengte van de Euleriaanse cyclus langer is dan die van de naïeve cyclus, wat logisch is.

print('Length of Eulerian circuit: {}'.format(len(euler_circuit)))

Length of Eulerian circuit: 158

CPP-oplossing berekenen

Tekst

Hier is een tekstuele weergave van de oplossing:

# Preview first 20 directions of CPP solution
for i, edge in enumerate(euler_circuit[0:20]):
    print(i, edge)

0 ('b_end_east', 'b_y', {'color': 'blue', 'estimate': 0, 'trail': 'b', 'distance': 1.32})
1 ('b_y', 'b_o', {'color': 'blue', 'estimate': 0, 'trail': 'b', 'distance': 0.08})
2 ('b_o', 'b_gy2', {'color': 'blue', 'estimate': 1, 'trail': 'b', 'distance': 0.05})
3 ('b_gy2', 'w_gy2', {'color': 'yellowgreen', 'estimate': 1, 'trail': 'gy2', 'distance': 0.03})
4 ('w_gy2', 'g_gy2', {'color': 'yellowgreen', 'estimate': 0, 'trail': 'gy2', 'distance': 0.05})
5 ('g_gy2', 'b_g', {'color': 'green', 'estimate': 0, 'trail': 'g', 'distance': 0.45})
6 ('b_g', 'b_w', {'color': 'blue', 'estimate': 0, 'trail': 'b', 'distance': 0.16})
7 ('b_w', 'b_gy2', {'color': 'blue', 'estimate': 0, 'trail': 'b', 'distance': 0.41})
8 ('b_gy2', 'w_gy2', {'color': 'yellowgreen', 'estimate': 1, 'trail': 'gy2', 'distance': 0.03})
9 ('w_gy2', 'b_w', {'color': 'gray', 'estimate': 0, 'trail': 'w', 'distance': 0.42})
10 ('b_w', 'w_rs', {'color': 'gray', 'estimate': 1, 'trail': 'w', 'distance': 0.06})
11 ('w_rs', 'g_rs', {'color': 'red', 'estimate': 0, 'trail': 'rs', 'distance': 0.18})
12 ('g_rs', 'b_g', {'color': 'green', 'estimate': 1, 'trail': 'g', 'distance': 0.05})
13 ('b_g', 'b_rs', {'color': 'blue', 'estimate': 0, 'trail': 'b', 'distance': 0.07})
14 ('b_rs', 'g_rs', {'color': 'red', 'estimate': 0, 'trail': 'rs', 'distance': 0.11})
15 ('g_rs', 'g_rc', {'color': 'green', 'estimate': 0, 'trail': 'g', 'distance': 0.45})
16 ('g_rc', 'g_gy1', {'color': 'green', 'estimate': 0, 'trail': 'g', 'distance': 0.37})
17 ('g_gy1', 'g_rt', {'color': 'green', 'estimate': 0, 'trail': 'g', 'distance': 0.26})
18 ('g_rt', 'g_w', {'color': 'green', 'estimate': 0, 'trail': 'g', 'distance': 0.31})
19 ('g_w', 'o_w_1', {'color': 'gray', 'estimate': 0, 'trail': 'w', 'distance': 0.18})

Je ziet al snel dat het algoritme niet erg trouw is aan een specifieke route: het springt vrij vlot van de ene naar de andere. Een uitbreiding van deze aanpak zou wat slimmer kunnen zijn door een vorm van routetrouw in de doelfunctie te verwerken, zodat het daadwerkelijk lopen van deze route beter behapbaar wordt.

Statistieken

Laten we even in je oplossing kijken om te zien hoe redelijk die eruitziet.
(Niet nodig om bij deze breedsprakige code stil te staan, het gaat om de geprinte output)

# Computing some stats
total_mileage_of_circuit = sum([edge[2]['distance'] for edge in euler_circuit])
total_mileage_on_orig_trail_map = sum(nx.get_edge_attributes(g, 'distance').values())
_vcn = pd.value_counts(pd.value_counts([(e[0]) for e in euler_circuit]), sort=False)
node_visits = pd.DataFrame({'n_visits': _vcn.index, 'n_nodes': _vcn.values})
_vce = pd.value_counts(pd.value_counts([sorted(e)[0] + sorted(e)[1] for e in nx.MultiDiGraph(euler_circuit).edges()]))
edge_visits = pd.DataFrame({'n_visits': _vce.index, 'n_edges': _vce.values})

# Printing stats
print('Mileage of circuit: {0:.2f}'.format(total_mileage_of_circuit))
print('Mileage on original trail map: {0:.2f}'.format(total_mileage_on_orig_trail_map))
print('Mileage retracing edges: {0:.2f}'.format(total_mileage_of_circuit-total_mileage_on_orig_trail_map))
print('Percent of mileage retraced: {0:.2f}%\n'.format((1-total_mileage_of_circuit/total_mileage_on_orig_trail_map)*-100))

print('Number of edges in circuit: {}'.format(len(euler_circuit)))
print('Number of edges in original graph: {}'.format(len(g.edges())))
print('Number of nodes in original graph: {}\n'.format(len(g.nodes())))

print('Number of edges traversed more than once: {}\n'.format(len(euler_circuit)-len(g.edges())))  

print('Number of times visiting each node:')
print(node_visits.to_string(index=False))

print('\nNumber of times visiting each edge:')
print(edge_visits.to_string(index=False))

Mileage of circuit: 33.59
Mileage on original trail map: 25.76
Mileage retracing edges: 7.83
Percent of mileage retraced: 30.40%

Number of edges in circuit: 158
Number of edges in original graph: 123
Number of nodes in original graph: 77

Number of edges traversed more than once: 35

Number of times visiting each node:
n_nodes  n_visits
     18         1
     38         2
     20         3
      1         4

Number of times visiting each edge:
n_edges  n_visits
     88         1
     35         2

CPP-oplossing visualiseren

Hoewel NetworkX ook functionaliteit biedt om grafen te visualiseren, zijn ze daarin opvallend bescheiden:

NetworkX biedt basisfunctionaliteit voor het visualiseren van grafen, maar het belangrijkste doel is graafanalyse, niet visualisatie. In de toekomst kan visualisatiefunctionaliteit uit NetworkX worden verwijderd of alleen beschikbaar zijn als add-onpakket.

Goede graafvisualisatie is lastig, en we raden sterk aan om grafen te visualiseren met tools die daarvoor zijn gemaakt. Bekende voorbeelden van volwaardige visualisatietools zijn Cytoscape, Gephi, Graphviz en, voor LaTeX-zetwerk, PGF/TikZ.

Dat gezegd hebbende: de ingebouwde NetworkX-tekenfunctionaliteit met matplotlib is krachtig genoeg om globaal te verkennen en te verifiëren, dus in deze tutorial blijf je bij NetworkX draw.

Ik heb graphviz en de dot-beschrijvingstaal gebruikt om de oplossing te visualiseren in mijn Python-pakket postman_problems. Het kostte wat werk om de NetworkX-graafstructuur naar een dot-graaf om te zetten, maar het levert wel betere kwaliteit en meer controle over visualisaties op.

CPP-graaf maken

Je eerste stap is om de lijst met te lopen edges in de Euler-cyclus om te zetten naar een edge-lijst met plotvriendelijke attributen.

create_cpp_edgelist maakt een edge-lijst met extra attributen die je gebruikt voor het plotten:

  • sequence: registreert in welke volgorde we elke edge lopen.
  • visits: is simpelweg het aantal keer dat we een bepaalde edge lopen.
def create_cpp_edgelist(euler_circuit):
    """
    Create the edgelist without parallel edge for the visualization
    Combine duplicate edges and keep track of their sequence and # of walks
    Parameters:
        euler_circuit: list[tuple] from create_eulerian_circuit
    """
    cpp_edgelist = {}

    for i, e in enumerate(euler_circuit):
        edge = frozenset([e[0], e[1]])

        if edge in cpp_edgelist:
            cpp_edgelist[edge][2]['sequence'] += ', ' + str(i)
            cpp_edgelist[edge][2]['visits'] += 1

        else:
            cpp_edgelist[edge] = e
            cpp_edgelist[edge][2]['sequence'] = str(i)
            cpp_edgelist[edge][2]['visits'] = 1

    return list(cpp_edgelist.values())

Laten we de CPP-edge-lijst maken:

cpp_edgelist = create_cpp_edgelist(euler_circuit)

Zoals verwacht heeft je edge-lijst hetzelfde aantal edges als de oorspronkelijke graaf.

print('Number of edges in CPP edge list: {}'.format(len(cpp_edgelist)))

Number of edges in CPP edge list: 123

De CPP-edge-lijst lijkt op euler_circuit, maar dan met een paar extra attributen.

# Preview CPP plot-friendly edge list
cpp_edgelist[0:3]

[('rh_end_tt_4',
  'nature_end_west',
  {'color': 'black',
   'distance': 0.2,
   'estimate': 0,
   'sequence': '73',
   'trail': 'tt',
   'visits': 1}),
 ('rd_end_south',
  'b_rd',
  {'color': 'red',
   'distance': 0.13,
   'estimate': 0,
   'sequence': '95',
   'trail': 'rd',
   'visits': 1}),
 ('w_gy1',
  'w_rc',
  {'color': 'gray',
   'distance': 0.33,
   'estimate': 0,
   'sequence': '151',
   'trail': 'w',
   'visits': 1})]

Nu maken we de graaf:

# Create CPP solution graph
g_cpp = nx.Graph(cpp_edgelist)

Visualisatie 1: Terugvolgen

Hier laat je zien welke edges één keer worden gelopen (grijs) en meer dan één keer (blauw). Dit is de "juiste" versie van de visualisatie uit 2.4 die de naïeve (hemelsbreed) verbindingen tussen de oneven knooppaaren liet zien (rood). Dat is hier gecorrigeerd door voor elk paar oneven knopen het kortste pad te volgen via edges die daadwerkelijk bestaan.

Als de optimalisatie goed is, zouden deze blauwe lijnen de kortst mogelijke afstand moeten vertegenwoordigen. Concreet: de minimale afstand die nodig is om een matching van de oneven knopen te genereren.

plt.figure(figsize=(14, 10))

visit_colors = {1:'lightgray', 2:'blue'}
edge_colors = [visit_colors[e[2]['visits']] for e in g_cpp.edges(data=True)]
node_colors = ['red'  if node in nodes_odd_degree else 'lightgray' for node in g_cpp.nodes()]

nx.draw_networkx(g_cpp, pos=node_positions, node_size=20, node_color=node_colors, edge_color=edge_colors, with_labels=False)
plt.axis('off')
plt.show()

graaf

Visualisatie 2: Volgorde van de CPP-oplossing

Hier plot je de oorspronkelijke graaf (het routenetwerk) met annotaties van de volgordenummers waarin we de paden lopen volgens de CPP-oplossing. Meerdere nummers geven paden aan waarop we moeten terugkeren.

Je begint op de blauwe route rechtsonder (de 0e en 157e richting).

plt.figure(figsize=(14, 10))

edge_colors = [e[2]['color'] for e in g_cpp.edges(data=True)]
nx.draw_networkx(g_cpp, pos=node_positions, node_size=10, node_color='black', edge_color=edge_colors, with_labels=False, alpha=0.5)

bbox = {'ec':[1,1,1,0], 'fc':[1,1,1,0]}  # hack to label edges over line (rather than breaking up line)
edge_labels = nx.get_edge_attributes(g_cpp, 'sequence')
nx.draw_networkx_edge_labels(g_cpp, pos=node_positions, edge_labels=edge_labels, bbox=bbox, font_size=6)

plt.axis('off')
plt.show()

graaf

Visualisatie 3: Animatie

De animatie hieronder volgt de Euler-cyclus van begin tot eind. Edges kleuren zwart de eerste keer dat ze gelopen worden en rood de tweede keer.

Let op dat deze gif overlappende of zeer korte edges niet altijd goed laat zien. Een robuustere visualisatiebibliotheek zoals graphviz kan dit verbeteren door splines te tekenen in plaats van rechte lijnen tussen knopen.

De code die dit maakt staat hieronder als referentie.

graaf

Eerst wordt voor elke richting (gelopen edge) uit de CPP-oplossing een PNG-afbeelding gemaakt.

visit_colors = {1:'black', 2:'red'}
edge_cnter = {}
g_i_edge_colors = []
for i, e in enumerate(euler_circuit, start=1):

    edge = frozenset([e[0], e[1]])
    if edge in edge_cnter:
        edge_cnter[edge] += 1
    else:
        edge_cnter[edge] = 1

    # Full graph (faded in background)
    nx.draw_networkx(g_cpp, pos=node_positions, node_size=6, node_color='gray', with_labels=False, alpha=0.07)

    # Edges walked as of iteration i
    euler_circuit_i = copy.deepcopy(euler_circuit[0:i])
    for i in range(len(euler_circuit_i)):
        edge_i = frozenset([euler_circuit_i[i][0], euler_circuit_i[i][1]])
        euler_circuit_i[i][2]['visits_i'] = edge_cnter[edge_i]
    g_i = nx.Graph(euler_circuit_i)
    g_i_edge_colors = [visit_colors[e[2]['visits_i']] for e in g_i.edges(data=True)]

    nx.draw_networkx_nodes(g_i, pos=node_positions, node_size=6, alpha=0.6, node_color='lightgray', with_labels=False, linewidths=0.1)
    nx.draw_networkx_edges(g_i, pos=node_positions, edge_color=g_i_edge_colors, alpha=0.8)

    plt.axis('off')
    plt.savefig('fig/png/img{}.png'.format(i), dpi=120, bbox_inches='tight')
    plt.close()

Daarna worden de PNG-afbeeldingen aan elkaar gezet tot de gif hierboven.

Eerst sorteer je de PNG's in de volgorde van 0 tot en met 157. Vervolgens zet je ze met imageio aan 3 frames per seconde aan elkaar tot een gif.

import glob
import numpy as np
import imageio
import os

def make_circuit_video(image_path, movie_filename, fps=5):
    # sorting filenames in order
    filenames = glob.glob(image_path + 'img*.png')
    filenames_sort_indices = np.argsort([int(os.path.basename(filename).split('.')[0][3:]) for filename in filenames])
    filenames = [filenames[i] for i in filenames_sort_indices]

    # make movie
    with imageio.get_writer(movie_filename, mode='I', fps=fps) as writer:
        for filename in filenames:
            image = imageio.imread(filename)
            writer.append_data(image)

make_circuit_video('fig/png/', 'fig/gif/cpp_route_animation.gif', fps=3)

Volgende stappen

Gefeliciteerd, je hebt deze tutorial voltooid en het Chinese Postman-probleem in Python opgelost. Je hebt veel terrein bestreken (33,6 mijl aan paden om precies te zijn). Wil je dieper in de netwerkfundamentals duiken, bekijk dan de Datacamp-cursus Network Analysis in Python, die de kernconcepten grondiger behandelt.

Aarzel niet om de NetworkX-documentatie te raadplegen voor meer over het creëren, manipuleren en doorlopen van deze complexe netwerken. De documentatie is uitgebreid, met veel voorbeelden en een reeks tutorials.

Als je de CPP op je eigen graaf wilt oplossen, heb ik de functionaliteit uit deze tutorial verpakt in het Python-pakket postman_problems op Github. Je kunt ook de codeblokken uit deze tutorial combineren met een andere edge- en knopenlijst, maar met het postman_problems-pakket kom je er waarschijnlijk sneller en netter.

Op een dag wil ik hier ook de uitbreidingen van de CPP (Rural en Windy Postman Problem) implementeren. Ik heb ook grootse plannen om over deze uitbreidingen en mijn ervaringen met het testen van de routes op de paden te schrijven op mijn blog hier. Een andere toepassing die ik wil verkennen en beschrijven is het gebruiken van lat/long-coördinaten om een mechanisme te ontwikkelen (of gebruiken) dat turn-by-turn-navigatie naar mijn Garmin-horloge stuurt.

En natuurlijk nog één laatste volgende stap: naar buiten gaan en de route trailrunnen!

Wil je meer leren over netwerken in Python, bekijk dan deze DataCamp-cursussen:

Introduction to Network Analysis in Python

Intermediate Network Analysis in Python

Referenties

1: Edmonds, Jack (1965). "Paths, trees, and flowers". Canad. J. Math. 17: 449–467.
2: Galil, Z. (1986). "Efficient algorithms for finding maximum matching in graphs". ACM Computing Surveys. Vol. 18, No. 1: 23-38.

Onderwerpen

PythonData ScienceDatavisualisatie

Gerelateerd

Meer bekijken →