Tutorial de optimización de grafos con Python NetworkX

En este tutorial vas a abordar un problema clásico de la teoría de grafos: el problema del cartero chino. Aunque algunos componentes del algoritmo son conceptualmente sencillos, su cálculo puede ser exigente. Aun así, para seguir este tutorial solo necesitas conocimientos previos básicos de Python: no hace falta dominar matemáticas avanzadas, informática ni teoría de grafos.

Primero repasaremos los elementos básicos de los grafos (nodos, aristas, caminos, etc.) y resolveremos el problema sobre un grafo real (la red de senderos de un parque estatal) usando la librería NetworkX en Python. Nos centraremos en los conceptos clave y en su implementación. Para quien quiera profundizar, incluimos lecturas adicionales sobre los entresijos de la optimización.

Por qué optimizar grafos

El problema

Probablemente te suene el problema del viajante, que consiste en encontrar la ruta más corta (por ejemplo, carreteras) que conecte un conjunto de nodos (por ejemplo, ciudades). Menos conocido pero muy parecido es el problema del cartero chino (CPP), también llamado problema de inspección de rutas o de enrutamiento por aristas. El objetivo del CPP es encontrar el recorrido más corto que cubra todos los enlaces (carreteras) de un grafo al menos una vez. Si se puede hacer sin pasar dos veces por la misma carretera, perfecto; ese es el escenario ideal y el problema es sencillo. Pero si es inevitable repetir algunos tramos, necesitarás algo de matemática para hallar la ruta más corta que recorra todas las carreteras al menos una vez con la menor distancia total.

Motivación personal

(Lo que sigue es una nota personal: cursi, gamberra y 100% innecesaria para aprender optimización de grafos en Python)

Tenía una aplicación real para resolver este problema: conseguir el rango de Giantmaster Marathoner.

¿Qué es un Giantmaster? Un Giantmaster es quien (canino o humano) ha recorrido todos los senderos del Sleeping Giant State Park en Hamden, CT (vecino de mi ciudad natal, Wallingford)... a lo largo de su vida. Un Giantmaster Marathoner es quien los recorre todos en un solo día.

Gracias al minucioso registro de la Sleeping Giant Park Association, puedes ver el listado completo de Giantmasters y su nivel aquí. Admito que esto me motivó bastante para arrancar este proyecto paralelo y salir a correr por los senderos. Yo mismo conseguí el estatus de Giantmaster en el invierno de 2006, cuando era un joven voluntario del Sleeping Giant Trail Crew (me alegró verlo registrado en el archivo de SG), pero desde entonces han surgido nuevos retos. Las categorías de Giantmaster de 12 meses y 4 estaciones son impresionantes y tentadoras, pero me exigirían más viajes desde mi casa actual (DC) a mi hogar de formación (CT) de los que puedo permitirme... y además son menos interesantes para la optimización de grafos, así que ¡Giantmaster Marathon al canto!

Como referencia, abajo tienes el mapa de senderos de Sleeping Giant:

Introducción a los grafos

Lo bueno de los grafos es que sus conceptos y su terminología suelen ser intuitivos. Aun así, aquí tienes algo de jerga básica:

Los grafos son estructuras que modelan relaciones entre objetos. En este tutorial llamaremos nodos a los objetos y aristas a las conexiones entre ellos. Ten en cuenta que aristas y nodos reciben distintos nombres que, en general, significan lo mismo:

node == vertex == point
edge == arc == link
Con la tecnología de

El grafo de partida es no dirigido. Es decir, las aristas no tienen orientación: son bidireccionales. Por ejemplo: A<--->B == B<--->A.
En cambio, el grafo que podrías crear para especificar el camino más corto para recorrer todos los senderos podría ser un grafo dirigido, donde el orden y la dirección de las aristas importan. Por ejemplo: A--->B != B--->A.

El grafo también es un grafo ponderado por aristas, donde la distancia (en millas) entre cada par de nodos adyacentes representa el peso de una arista. Esto se gestiona como un atributo de arista llamado "distance".

El grado es el número de aristas incidentes (que tocan) un nodo. Los nodos son de grado impar cuando este número es impar, y de grado par cuando es par.

La solución a este CPP será un recorrido euleriano: un grafo en el que se puede formar un ciclo que pasa por cada arista exactamente una vez y vuelve del nodo inicial a sí mismo (sin retroceder). Un recorrido euleriano también se conoce con varios nombres:

Eulerian tour == Eulerian circuit == Eulerian cycle
Con la tecnología de

Un emparejamiento es un subconjunto de aristas en el que ningún nodo aparece más de una vez. Un emparejamiento de peso mínimo encuentra el emparejamiento con la suma de pesos de arista más baja posible.

NetworkX: manipulación y análisis de grafos

NetworkX es el paquete de Python más popular para manipular y analizar grafos. Hay otros con un nivel básico similar, como igraph (con bindings para R y C++). Sin embargo, NetworkX cuenta con los algoritmos de grafos más sólidos que necesitaba para resolver el CPP.

Instalar paquetes

Si ya has hecho algo de análisis de datos en Python o usas la distribución Anaconda, probablemente tengas pandas y matplotlib. Puede que te falte networkx. Estas deberían ser las únicas dependencias fuera de la librería estándar de Python que necesitas para seguir este tutorial. Se instalan fácilmente con pip:

pip install pandas
pip install networkx
pip install matplotlib
Con la tecnología de

Con esto tendrás lo necesario por ahora. imageio y numpy se importan al final para crear la animación GIF de la solución del CPP. La animación está incrustada en este artículo, así que estos paquetes son opcionales.

import itertools
import copy
import networkx as nx
import pandas as pd
import matplotlib.pyplot as pltCon la tecnología de

Cargar datos

Lista de aristas

La lista de aristas es una estructura sencilla que usarás para crear el grafo. Cada fila representa una arista con algunos atributos.

node1 y node2: nombres de los nodos conectados.
trail: atributo de arista con el nombre abreviado del sendero para cada arista. Por ejemplo: rs = red square
distance: atributo de arista con la longitud del tramo en millas.
color: color del sendero para el gráfico.
estimate: atributo que indica si la distancia de esa arista es una estimación a ojo del mapa (1=sí, 0=no) porque algunas distancias no se proporcionan. Es solo informativo; no se usa en el análisis.

# Cargar la lista de aristas alojada en Gist
edgelist = pd.read_csv('https://gist.githubusercontent.com/brooksandrew/e570c38bcc72a8d102422f2af836513b/raw/89c76b2563dbc0e88384719a35cba0dfc04cd522/edgelist_sleeping_giant.csv')
Con la tecnología de

# Vista previa
edgelist.head(10)
Con la tecnología de

	node1	node2	trail	distance	color
0	rs_end_north	v_rs	rs	0.30	red
1	v_rs	b_rs	rs	0.21	red
2	b_rs	g_rs	rs	0.11	red
3	g_rs	w_rs	rs	0.18	red
4	w_rs	o_rs	rs	0.21	red
5	o_rs	y_rs	rs	0.12	red
6	y_rs	rs_end_south	rs	0.39	red
7	rc_end_north	v_rc	rc	0.70	red
8	v_rc	b_rc	rc	0.04	red
9	b_rc	g_rc	rc	0.15	red

Lista de nodos

Las listas de nodos suelen ser opcionales en networkx y otras librerías cuando ya tienes la lista de aristas, porque los nombres de los nodos aparecen en las dos primeras columnas. En este caso, sin embargo, queremos añadir atributos de nodo: las coordenadas X, Y de los nodos (intersecciones de senderos) para poder trazar el grafo con el mismo diseño que el mapa.

Pasé una tarde anotándolas manualmente calcando la imagen con GIMP:

id: nombre del nodo que corresponde a node1 y node2 en la lista de aristas.
X: posición/coordenada horizontal del nodo respecto a la esquina superior izquierda.
Y posición/coordenada vertical del nodo respecto a la esquina superior izquierda.

Nota sobre la generación de las listas de nodos y aristas

Crear los nombres de los nodos también llevó trabajo manual. Cada nodo representa una intersección de dos o más senderos. Cuando fue posible, el nodo se nombró como trail1_trail2, poniendo trail1 antes que trail2 en orden alfabético.

Se complica cuando los mismos senderos se cruzan más de una vez. Por ejemplo, los senderos Orange y White. En esos casos añadí un sufijo _2 o _3 al nombre del nodo. Así, tienes dos nombres distintos para dos intersecciones distintas de Orange y White: o_w y o_w_2.

Esto requirió bastante prueba y error, comparando los gráficos generados con X,Y con el mapa real.

# Cargar la lista de nodos alojada en Gist
nodelist = pd.read_csv('https://gist.githubusercontent.com/brooksandrew/f989e10af17fb4c85b11409fea47895b/raw/a3a8da0fa5b094f1ca9d82e1642b384889ae16e8/nodelist_sleeping_giant.csv')
Con la tecnología de

# Vista previa
nodelist.head(5)
Con la tecnología de

	id	X	Y
0	b_bv	1486	732
1	b_bw	716	1357
2	b_end_east	3164	1111
3	b_end_west	141	1938
4	b_g	1725	771

Crear el grafo

Ahora usas la lista de aristas y la de nodos para crear un objeto grafo en networkx.

# Crear grafo vacío
g = nx.Graph()
Con la tecnología de

Recorre las filas de la lista de aristas y añade cada arista y sus atributos correspondientes al grafo g.

# Añadir aristas y atributos
for i, elrow in edgelist.iterrows():
    g.add_edge(elrow[0], elrow[1], attr_dict=elrow[2:].to_dict())
Con la tecnología de

Para ilustrarlo, imprime los valores de la última fila añadida al grafo g:

# Ejemplo de lista de aristas
print(elrow[0]) # node1
print(elrow[1]) # node2
print(elrow[2:].to_dict()) # diccionario de atributos de arista
Con la tecnología de

o_gy2
y_gy2
{'color': 'yellowgreen', 'estimate': 0, 'trail': 'gy2', 'distance': 0.12}
Con la tecnología de

De forma similar, recorre las filas de la lista de nodos y añade estos atributos de nodo.

# Añadir atributos de nodo
for i, nlrow in nodelist.iterrows():
    g.node[nlrow['id']] = nlrow[1:].to_dict()
Con la tecnología de

Aquí tienes un ejemplo de la última fila de la lista de nodos:

# Ejemplo de lista de nodos
print(nlrow)
Con la tecnología de

id    y_rt
X      977
Y     1666
Name: 76, dtype: objectCon la tecnología de

Inspeccionar el grafo

Aristas

Las aristas se representan como una lista de tuplas de longitud 3. Los dos primeros elementos son los nombres de los nodos conectados por la arista. El tercero es el diccionario de atributos.

# Vista previa de las 5 primeras aristas
g.edges(data=True)[0:5]
Con la tecnología de

[('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'})]
Con la tecnología de

Nodos

Del mismo modo, los nodos se representan como una lista de tuplas de longitud 2. El primer elemento es el identificador del nodo, seguido del diccionario de atributos.

# Vista previa de los 10 primeros nodos
g.nodes(data=True)[0:10]
Con la tecnología de

[('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})]
Con la tecnología de

Estadísticas resumidas

Imprime algunas estadísticas antes de visualizar el grafo.

print('# of edges: {}'.format(g.number_of_edges()))
print('# of nodes: {}'.format(g.number_of_nodes()))
Con la tecnología de

# of edges: 123
# of nodes: 77Con la tecnología de

Visualizar el grafo

Ajustar colores y distribución

Posiciones: Primero necesitas transformar las posiciones de los nodos del grafo en un diccionario. Esto te permitirá recrear el grafo con el mismo diseño que el mapa de senderos. Se niega Y para transformar el origen del eje Y de la esquina superior izquierda a la inferior izquierda.

# Definir el diccionario de posiciones de nodos para el gráfico
node_positions = {node[0]: (node[1]['X'], -node[1]['Y']) for node in g.nodes(data=True)}

# Vista previa de node_positions con un pequeño truco (no hay head/slice para dicts)
dict(list(node_positions.items())[0:5])
Con la tecnología de

{'b_rd': (268, -1744),
 'g_rt': (908, -1378),
 'o_gy1': (1130, -1297),
 'rh_end_tt_2': (550, -1608),
 'rs_end_south': (1865, -1598)}
Con la tecnología de

Colores: Ahora extraes los colores de las aristas a una lista simple para visualizar los senderos por su color.

# Lista de colores de aristas para el gráfico
edge_colors = [e[2]['color'] for e in g.edges(data=True)]

# Vista previa de las 10 primeras
edge_colors[0:10]
Con la tecnología de

['red',
 'gray',
 'yellowgreen',
 'gray',
 'yellowgreen',
 'blue',
 'black',
 'yellowgreen',
 'gray',
 'gray']
Con la tecnología de

Gráfico

Ya puedes crear un gráfico que se alinea bien con el mapa de senderos de 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('Representación en grafo del mapa de senderos de Sleeping Giant', size=15)
plt.show()
Con la tecnología de

Este grafo obviamente no refleja todas las curvas y meandros de los senderos, pero no pasa nada: eso queda recogido en el atributo distance, que es el que usamos para calcular. La visualización sí captura la distancia entre nodos (intersecciones) en línea recta, que parece una aproximación razonable.

Resumen del algoritmo del CPP

Bien, ahora que ya definiste conceptos y creaste el grafo, ¿cómo encuentras el camino más corto?

Resolver el problema del cartero chino es conceptualmente sencillo:

Encuentra todos los nodos de grado impar (muy fácil).
(Todas las intersecciones donde el número de senderos que llegan es impar)
Añade aristas al grafo de forma que todos los nodos de grado impar pasen a ser pares. Estas aristas añadidas deben duplicar aristas del grafo original (supondremos que aquí no hay campo a través). El conjunto de aristas añadidas debe sumar la mínima distancia posible (difícil... de hecho, NP-difícil).
(En pocas palabras: minimiza cuánto tienes que desandar en una ruta que recorre todos los senderos)
Dado un punto de inicio, calcula el recorrido euleriano sobre el conjunto de datos aumentado (moderadamente fácil).
(Una vez sabemos por qué senderos habrá que desandar, calculamos la ruta de principio a fin)

Supuestos y simplificaciones

Aunque relajando los supuestos de abajo podrías obtener una ruta algo más corta y precisa, añadiría complejidad fuera del objetivo de este tutorial, centrado en el CPP.

Supuesto 1: solo senderos obligatorios

Como ves en el mapa, hay carreteras en los bordes del parque que podrían conectar senderos, especialmente los rojos. También hay senderos (Horseshoe y no señalizados) que no son obligatorios según el registro Giantmaster, pero podrían ayudar a evitar largos retrocesos. Incluir senderos opcionales es, de hecho, una variante establecida del CPP llamada Rural Postman Problem. En este tutorial ignoramos los opcionales y nos ceñimos a los obligatorios.

Supuesto 2: subir == bajar

El CPP asume que el coste de caminar un sendero equivale a su distancia, independientemente de la dirección. Sin embargo, algunos senderos son bastante empinados y subir puede requerir más esfuerzo que bajar. Podrías incorporar una métrica que combine distancia y desnivel sobre un grafo dirigido en una extensión del CPP llamada Windy Postman Problem.

Supuesto 3: sin aristas paralelas (senderos)

Aunque es posible, incluir aristas paralelas (varios senderos entre los mismos dos nodos) complica el cálculo. Por suerte aquí solo pasa dos veces (Blue <=> Red Diamond y Blue <=> Tower Trail). Se resuelve con un pequeño truco en la lista de aristas: se incluyen nodos duplicados con sufijo _dupe para capturar todos los senderos manteniendo la unicidad en las aristas. La implementación del CPP en el paquete postman_problems que escribí gestiona bien las aristas paralelas si quieres resolver el CPP en tu propio grafo con muchas aristas paralelas.

Paso 1 del CPP: encontrar los nodos de grado impar

Es básicamente un recuento. Verás que 36 de los 76 nodos tienen grado impar. La mayoría son senderos sin salida (grado 1) e intersecciones de 3 senderos. Hay unos pocos de grado 5.

# Calcular la lista de nodos de grado impar
nodes_odd_degree = [v for v, d in g.degree_iter() if d % 2 == 1]

# Vista previa
nodes_odd_degree[0:5]
Con la tecnología de

['rs_end_south', 'rc_end_north', 'v_end_east', 'rh_end_south', 'b_end_east']
Con la tecnología de

# Recuentos
print('Number of nodes of odd degree: {}'.format(len(nodes_odd_degree)))
print('Number of total nodes: {}'.format(len(g.nodes())))
Con la tecnología de

Number of nodes of odd degree: 36
Number of total nodes: 77
Con la tecnología de

Paso 2 del CPP: encontrar los pares a mínima distancia

Aquí está el meollo del problema. Lo dividiremos en 5 partes:

Calcula todos los pares posibles de nodos de grado impar.
Calcula el camino más corto entre cada par del punto 1.
Crea un grafo completo conectando cada par del 1. con el atributo de distancia de camino más corto del 2.
Calcula un emparejamiento de peso mínimo del grafo del 3.
(Se reduce a emparejar los nodos impares para que la suma de distancias entre pares sea lo más pequeña posible).
Aumenta el grafo original con los caminos más cortos entre los pares del 4.

Paso 2.1: calcular los pares de nodos

Usas la función itertools combination para obtener todos los pares posibles de nodos de grado impar. El grafo es no dirigido, así que el orden no importa: por ejemplo, (a,b) == (b,a).

# Calcular todos los pares de nodos impares (lista de tuplas)
odd_node_pairs = list(itertools.combinations(nodes_odd_degree, 2))

# Vista previa de pares de nodos impares
odd_node_pairs[0:10]
Con la tecnología de

[('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')]
Con la tecnología de

# Recuento
print('Number of pairs: {}'.format(len(odd_node_pairs)))
Con la tecnología de

Number of pairs: 630
Con la tecnología de

Confirmemos que este número de pares es correcto con la combinatoria de abajo. Por suerte, solo tienes 630 pares que optimizar. El tiempo de cómputo para este ejemplo es trivial (un par de segundos).

Sin embargo, si tuvieras 3.600 nodos impares, saldrían ~6,5 millones de pares a optimizar. Es decir, ~10.000× más salida ante una entrada 100× mayor.

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

Paso 2.2: calcular los caminos más cortos entre pares

Este es el primer paso con cálculo real. Por suerte, networkx implementa de forma cómoda el algoritmo de Dijkstra para calcular el camino más corto entre dos nodos. Lo aplicas a cada par (los 630) calculado en 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
Con la tecnología de

# Calcular caminos más cortos. Devuelve un diccionario con pares de nodos como clave y la distancia mínima como valor.
odd_node_pairs_shortest_paths = get_shortest_paths_distances(g, odd_node_pairs, 'distance')

# Vista previa con un pequeño truco (no hay head/slice para dicts)
dict(list(odd_node_pairs_shortest_paths.items())[0:10])
Con la tecnología de

{('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}
Con la tecnología de

Paso 2.3: crear el grafo completo

Un grafo completo es aquel en el que todos los nodos están conectados entre sí por una arista única.

Aquí tienes un ejemplo básico de Wikipedia con 7 nodos y 21 aristas (7 sobre 2):

El grafo que creas a continuación tiene 36 nodos y 630 aristas con su peso (distancia) correspondiente.

Definimos create_complete_graph para calcularlo. El parámetro flip_weights transforma el atributo distance en el atributo weight negándolo, de forma que números pequeños (distancias grandes) se convierten en pesos bajos y números grandes (distancias pequeñas) en pesos altos. Suena contraintuitivo, pero es necesario para el paso 2.4, donde calculamos el emparejamiento de peso mínimo en el grafo completo.

Idealmente calcularíamos directamente el emparejamiento de peso mínimo, pero NetworkX solo implementa max_weight_matching, que maximiza el peso. Lo "hackeamos" negando (multiplicando por -1) distance para obtener weight. Así preservamos el orden y la escala por distancia, pero invertidos.

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
Con la tecnología de

# Generar el grafo completo
g_odd_complete = create_complete_graph(odd_node_pairs_shortest_paths, flip_weights=True)

# Recuentos
print('Number of nodes: {}'.format(len(g_odd_complete.nodes())))
print('Number of edges: {}'.format(len(g_odd_complete.edges())))
Con la tecnología de

Number of nodes: 36
Number of edges: 630
Con la tecnología de

Como apoyo visual, graficamos el grafo totalmente conectado de los nodos de grado impar. Observa que preservamos las coordenadas X, Y de cada nodo, pero las aristas no representan necesariamente senderos reales. Por ejemplo, dos nodos pueden estar conectados por una única arista aquí, mientras que el camino más corto entre ellos puede requerir 5 saltos por nodos de grado par (no mostrados).

# Graficar el grafo completo de nodos de grado impar
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('Grafo completo de nodos de grado impar')
plt.show()
Con la tecnología de

Paso 2.4: calcular el emparejamiento de peso mínimo

Este es el paso más complejo del CPP. Necesitas encontrar los pares de nodos de grado impar cuya suma (de distancias entre ellos) sea lo más pequeña posible. En tu caso, se trata de seleccionar las 18 aristas óptimas (36 nodos impares / 2) del grafo enmarañado generado en el 2.3.

Tanto la implementación como la intuición de esta optimización quedan fuera del alcance del tutorial... hablamos de 800+ líneas de código y mucha literatura académica.

Aun así, un apunte para quien tenga curiosidad:

Mil gracias a Joris van Rantwijk por escribir la implementación original en su blog allá por 2008. Me topé con el problema de forma similar y con la misma intención que Joris. De su post de 2008:

Como no encontré implementaciones en Perl de emparejamiento con peso máximo, me decidí alegremente a escribir algo yo mismo. Resultó que había infravalorado el problema, pero cuando me di cuenta del error ya estaba tan obsesionado que me negué a rendirme.

Yo sí me rendí. Por suerte, Joris no.

Desde entonces, este algoritmo de emparejamiento con peso máximo forma parte de NetworkX y se mantiene allí. Gracias también a los 10+ contribuyentes en GitHub que lo han mantenido.

Es un cálculo difícil e intensivo. El primer gran avance, en 1965, probó que el problema del emparejamiento máximo podía resolverse en tiempo polinómico. Lo publicó Jack Edmonds con quizá uno de los títulos de paper más bonitos: "Paths, trees, and flowers" [1]. Desde entonces se ha construido mucha literatura mejorando el procedimiento. El código de max_weight_matching en NetworkX se basa en Galil, Zvi (1986) [2] con un algoritmo de tiempo O(n³).

# Calcular el emparejamiento de peso mínimo.
# Nota: max_weight_matching usa por defecto el atributo 'weight' como peso a maximizar.
odd_matching_dupes = nx.algorithms.max_weight_matching(g_odd_complete, True)

print('Number of edges in matching: {}'.format(len(odd_matching_dupes)))
Con la tecnología de

Number of edges in matching: 36
Con la tecnología de

La salida (odd_matching_dupes) es un diccionario. Aunque hay 36 aristas en este emparejamiento, solo quieres 18. Cada pareja aparece dos veces (una con el nodo 1 como clave y otra con el nodo 2 como clave).

# Vista previa del emparejamiento con duplicados
odd_matching_dupes
Con la tecnología de

{'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'}
Con la tecnología de

Conviertes este diccionario en una lista de tuplas, ya que el grafo es no dirigido y el orden no importa. Al eliminar duplicados obtienes las 18 parejas únicas cuya suma de distancias es mínima.

# Convertir el emparejamiento a lista de tuplas deduplicadas
odd_matching = list(pd.unique([tuple(sorted([k, v])) for k, v in odd_matching_dupes.items()]))

# Recuento
print('Number of edges in matching (deduped): {}'.format(len(odd_matching)))
Con la tecnología de

Number of edges in matching (deduped): 18
Con la tecnología de

# Vista previa del emparejamiento deduplicado
odd_matching
Con la tecnología de

[('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')]
Con la tecnología de

Visualicemos estas parejas en el grafo completo del paso 2.3. Como antes, aunque las posiciones de nodos reflejan el grafo real (mapa), las distancias de las aristas (líneas azules) son en línea recta. La ruta real más corta entre dos nodos puede implicar varias aristas que giran y se alargan.

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

# Graficar el grafo completo de nodos de grado impar
nx.draw(g_odd_complete, pos=node_positions, node_size=20, alpha=0.05)

# Grafo superpuesto solo con las aristas del emparejamiento de peso mínimo
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('Emparejamiento de peso mínimo en el grafo completo')
plt.show()
Con la tecnología de

Para ver cómo encaja con el grafo original, trazamos esas mismas parejas mínimas (líneas azules) sobre el mapa de senderos (atenuado) en lugar del grafo completo. De nuevo, las líneas azules son en línea recta (no senderos reales). Aún queda trabajo para encontrar las aristas que componen el camino más corto entre cada pareja en el paso 3.

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

# Graficar el grafo original del mapa de senderos
nx.draw(g, pos=node_positions, node_size=20, alpha=0.1, node_color='black')

# Superponer las aristas del emparejamiento de peso mínimo
nx.draw(g_odd_complete_min_edges, pos=node_positions, node_size=20, alpha=1, node_color='red', edge_color='blue')

plt.title('Emparejamiento de peso mínimo sobre el grafo original')
plt.show()
Con la tecnología de

Paso 2.5: aumentar el grafo original

Ahora aumentas el grafo original con las aristas del emparejamiento del 2.4. Definimos una función sencilla para hacerlo, marcando además que estas nuevas aristas provienen del grafo aumentado. Necesitarás saberlo en el 3. al crear el circuito euleriano.

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
    """

    # Hacemos el grafo aumentado un MultiGraph para permitir aristas paralelas
    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
Con la tecnología de

Confirmemos que el grafo aumentado añade el número esperado (18) de aristas:

# Crear el grafo aumentado: añadir las aristas del emparejamiento mínimo a g
g_aug = add_augmenting_path_to_graph(g, odd_matching)

# Recuentos
print('Number of edges in original graph: {}'.format(len(g.edges())))
print('Number of edges in augmented graph: {}'.format(len(g_aug.edges())))
Con la tecnología de

Number of edges in original graph: 123
Number of edges in augmented graph: 141
Con la tecnología de

Y que ahora todos los nodos tienen grado par:

pd.value_counts(g_aug.degree())
Con la tecnología de

4    54
2    18
6     5
dtype: int64
Con la tecnología de

Paso 3 del CPP: calcular el circuito euleriano

Con un grafo de grado par, lo más duro ya está hecho. Como postuló Euler en 1736 con el problema de los Siete puentes de Königsberg, existe un camino que visita cada arista exactamente una vez si todos los nodos tienen grado par. Carl Hierholzer lo probó formalmente en la década de 1870.

Puede haber muchos circuitos eulerianos con la misma distancia. Con la función eulerian_circuit de NetworkX puedes cubrir el 90% del trabajo, pero hay limitaciones.

Limitaciones que vamos a resolver:

El grafo aumentado puede (y probablemente lo hará) contener aristas que no existen en el grafo original. Para obtener el circuito (sin campo a través), debes descomponer esas aristas aumentadas en el camino más corto a través de aristas existentes.
eulerian_circuit solo devuelve el orden en que visitamos cada nodo. No devuelve los atributos de las aristas necesarios para completar el circuito. Es importante porque necesitas llevar la cuenta de qué aristas ya has recorrido cuando exista más de una arista entre dos nodos.

Limitaciones que no resolveremos:

Para ahorrarte piernas, podrías relajar el supuesto del circuito euleriano de empezar y terminar en el mismo nodo. Un camino euleriano (caso general del circuito) también existe si hay exactamente dos nodos de grado impar. Ahorraría algo de doble recorrido... si te pudieran recoger en el otro extremo del parque. A día de hoy, NetworkX no ofrece un algoritmo de camino euleriano. El código de eulerian_circuit no es muy farragoso y podría adaptarse, pero aquí lo mantenemos simple.

Circuito ingenuo

Empecemos con la solución simple pero incompleta:

naive_euler_circuit = list(nx.eulerian_circuit(g_aug, source='b_end_east'))
Con la tecnología de

Como era de esperar, la longitud del circuito euleriano ingenuo es igual al número de aristas del grafo aumentado.

print('Length of eulerian circuit: {}'.format(len(naive_euler_circuit)))
Con la tecnología de

Length of eulerian circuit: 141
Con la tecnología de

La salida es una lista de tuplas que representan pares de nodos. Observa que el primer nodo de cada par coincide con el segundo del par anterior.

# Vista previa del circuito euleriano ingenuo
naive_euler_circuit[0:10]
Con la tecnología de

[('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')]
Con la tecnología de

Circuito correcto

Ahora definimos una función que usa el grafo original para indicarte qué senderos usar de A a B. Aunque el código sea algo verboso, la lógica es simple. Transformamos el circuito ingenuo —que incluía aristas inexistentes en el grafo original— en un circuito euleriano que usa solo aristas existentes.

Recorres cada arista del circuito euleriano ingenuo (naive_euler_circuit). Siempre que encuentres una arista que no exista en el grafo original, la sustituyes por la secuencia de aristas que forman el camino más corto entre sus nodos en el grafo original.

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':
            # Si `edge` existe en el grafo original, toma sus atributos y añade al circuito
            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))

            # Si `edge` no existe en el grafo original, toma el camino más corto y añade sus aristas
            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
Con la tecnología de

"Hackeamos" un poco la limitación 3 empezando el circuito euleriano en el extremo este del parque, en el sendero azul (nodo "b_end_east"). A la hora de correrlo, podrías saltarte la última indicación, que desanda ese tramo.

Añadimos impresiones verbosas para mostrar qué ocurre al sustituir aristas inexistentes del grafo aumentado por el camino más corto entre nodos con aristas reales.

# Crear el circuito euleriano
euler_circuit = create_eulerian_circuit(g_aug, g, 'b_end_east')
Con la tecnología de

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')]
Con la tecnología de

Verás que la longitud del circuito euleriano es mayor que la del circuito ingenuo, lo cual tiene sentido.

print('Length of Eulerian circuit: {}'.format(len(euler_circuit)))
Con la tecnología de

Length of Eulerian circuit: 158
Con la tecnología de

Calcular la solución del CPP

Texto

Aquí tienes un volcado en texto de la solución:

# Vista previa de las primeras 20 indicaciones de la solución del CPP
for i, edge in enumerate(euler_circuit[0:20]):
    print(i, edge)
Con la tecnología de

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})
Con la tecnología de

Se aprecia rápido que el algoritmo no es muy fiel a un sendero concreto: salta entre senderos con frecuencia. Una extensión interesante sería introducir cierta «fidelidad al sendero» en la función objetivo para hacer la ruta más manejable al correrla.

Estadísticas

Echemos un vistazo a la solución para ver si es razonable.
(No hace falta detenerse en el código verboso, solo en la salida impresa)

# Cálculo de algunas estadísticas
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})

# Impresión de estadísticas
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))
Con la tecnología de

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
Con la tecnología de

Visualizar la solución del CPP

Aunque NetworkX también permite visualizar grafos, son bastante humildes en este apartado:

NetworkX ofrece funcionalidad básica para visualizar grafos, pero su objetivo principal es el análisis. En el futuro, la visualización podría eliminarse o quedar como paquete adicional.

Una buena visualización de grafos es difícil; recomendamos usar herramientas dedicadas como Cytoscape, Gephi, Graphviz o, para LaTeX, PGF/TikZ.

Dicho esto, el dibujo integrado de NetworkX con matplotlib es suficiente para explorar visualmente grafos básicos, así que aquí seguimos con draw de NetworkX.

Usé graphviz y el lenguaje dot para visualizar la solución en mi paquete de Python postman_problems. Aunque costó convertir la estructura de NetworkX a dot, te da más calidad y control.

Crear el grafo del CPP

El primer paso es convertir la lista de aristas a recorrer en el circuito euleriano en una lista de aristas con atributos amigables para el gráfico.

create_cpp_edgelist crea una lista de aristas con atributos adicionales que usarás al trazar:

sequence: registra el orden en que caminamos cada arista.
visits: número de veces que caminamos una arista.

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())
Con la tecnología de

Creemos la lista de aristas del CPP:

cpp_edgelist = create_cpp_edgelist(euler_circuit)
Con la tecnología de

Como era de esperar, tu lista de aristas tiene el mismo número de aristas que el grafo original.

print('Number of edges in CPP edge list: {}'.format(len(cpp_edgelist)))
Con la tecnología de

Number of edges in CPP edge list: 123
Con la tecnología de

La lista del CPP se parece a euler_circuit, con algunos atributos extra.

# Vista previa de la lista de aristas para el gráfico
cpp_edgelist[0:3]
Con la tecnología de

[('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})]
Con la tecnología de

Ahora creamos el grafo:

# Crear el grafo de la solución del CPP
g_cpp = nx.Graph(cpp_edgelist)
Con la tecnología de

Visualización 1: pasos repetidos

Aquí ilustras qué aristas se recorren una vez (gris) y cuáles más de una (azul). Es la versión "correcta" de la visualización del 2.4, que mostraba las conexiones ingenuas (en línea recta) entre pares de nodos impares (rojo). Ahora lo corregimos siguiendo el camino más corto por aristas reales para cada par de nodos impares.

Si la optimización es buena, estas líneas azules deberían representar la menor distancia posible. En concreto, la distancia mínima necesaria para lograr un emparejamiento de los nodos impares.

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()
Con la tecnología de

Visualización 2: secuencia de la solución del CPP

Ahora trazas el grafo original (mapa de senderos) anotado con los números de secuencia en los que caminamos los tramos según la solución del CPP. Varios números indican tramos que hay que desandar.

Empiezas en el sendero azul abajo a la derecha (direcciones 0 y 157).

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]}  # truco para etiquetar aristas encima de la línea
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()
Con la tecnología de

Visualización 3: animación

La animación que recorre el circuito euleriano de principio a fin está incrustada abajo. Las aristas se colorean en negro la primera vez que se caminan y en rojo la segunda.

Ten en cuenta que este gif no hace plena justicia visual a aristas superpuestas o demasiado pequeñas. Una librería más robusta como graphviz podría resolverlo trazando splines en vez de líneas rectas.

Incluimos el código que la crea como referencia.

Primero se genera una imagen PNG por cada dirección (arista recorrida) de la solución del CPP.

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

    # Grafo completo (atenuado de fondo)
    nx.draw_networkx(g_cpp, pos=node_positions, node_size=6, node_color='gray', with_labels=False, alpha=0.07)

    # Aristas recorridas hasta la iteración 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()
Con la tecnología de

Después, las PNG se unen para crear el GIF.

Primero se ordenan las PNG de 0 a 157. Luego se ensamblan con imageio a 3 fotogramas por segundo para crear el gif.

import glob
import numpy as np
import imageio
import os

def make_circuit_video(image_path, movie_filename, fps=5):
    # ordenar nombres de archivo
    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]

    # crear vídeo
    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)Con la tecnología de

Siguientes pasos

Enhorabuena, has terminado este tutorial resolviendo el problema del cartero chino en Python. Hemos cubierto mucho terreno (33,6 millas de senderos, para ser exactos). Si quieres profundizar en los fundamentos de redes, puede interesarte el curso Network Analysis in Python de Datacamp, que trata los conceptos básicos con más detalle.

No dudes en consultar la documentación de NetworkX para aprender a crear, manipular y recorrer estas redes complejas. La documentación es completa y tiene bastantes ejemplos y una serie de tutoriales.

Si te apetece resolver el CPP en tu propio grafo, he empaquetado la funcionalidad de este tutorial en el paquete de Python postman_problems en GitHub. También puedes recomponer los bloques de código con otras listas de aristas y nodos, pero probablemente el paquete te lleve más rápido y con menos fricción.

Algún día me gustaría implementar aquí las extensiones del CPP (Rural y Windy Postman Problem). También tengo la ambición de escribir sobre estas extensiones y de probar las rutas en los senderos en mi blog aquí. Otra aplicación que quiero explorar y documentar es incorporar coordenadas lat/long para desarrollar (o usar) un mecanismo que envíe indicaciones giro a giro a mi reloj Garmin.

Y, por supuesto, un último siguiente paso: ¡salir ahí fuera a correr la ruta!

Si te gustaría aprender más sobre redes en Python, echa un vistazo a estos cursos de DataCamp:

Introduction to Network Analysis in Python

Intermediate Network Analysis in Python

Referencias

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.

Temas

Python Ciencia de datos Visualización de datos