Approximative solution based on Hierholzer's algorithm

.gitignore

@@ -1,2 +1,3 @@
 __pycache__
-in/
+in/
+.pytest_cache/

graph-test.py

@@ -0,0 +1,36 @@
+from graph import Graph
+from input import parse_file
+
+
+# Refacto : have a fixture file
+def read_graph(in_file: str) -> Graph:
+    vertices, edges = parse_file(in_file)
+    print(f"#E={len(edges)}, #V={len(vertices)}")
+    return Graph(vertices, edges)
+
+
+triangle = read_graph("instances/triangle.txt")
+line = read_graph("instances/line.txt")
+paris = read_graph("instances/paris_map.txt")
+
+## End refacto needed
+
+
+def odd():
+    assert triangle.is_eulerian()
+    assert not line.is_eulerian()
+    assert triangle.odd_vertices() == []
+    assert line.odd_vertices() == [0, 2]
+    assert not paris.is_eulerian()
+
+
+def next():
+    assert triangle.next(0, 1) == [1, 2]
+    assert triangle.next(1, 2) == [2, 0]  # Next
+    assert triangle.next(0, 0) == [2, 2]  # Reverse
+    assert line.next(1, 2) == [-1, None]  # End of the line
+
+
+def tests():
+    odd()
+    next()

graph.py

@@ -1,10 +1,15 @@
 from __future__ import annotations
 
+from typing import Union
+
 import matplotlib.pyplot as plt
 from matplotlib.collections import LineCollection
+from scipy.sparse import csr_array
 
 Coordinates = tuple[float, float]
-Edge = tuple[int, int, int, Coordinates, Coordinates]
+True_Edge = tuple[int, int, int, Coordinates, Coordinates]
+Virtual_Edge = tuple[int, int, int, Coordinates, Coordinates, list[True_Edge]]
+Edge = Union[None, True_Edge, Virtual_Edge]
 
 
 class Graph:
@@ -25,19 +30,21 @@ def __init__(
         """
         self.vertices = vertices
         self.edges = edges
+        self.real_edges = edges
+        self.adjency_computed = False
 
     def plot(self):
         """
         Plot the graph.
         """
-        weights = list(set(edge[2] for edge in self.edges))
+        weights = list(set(edge[2] for edge in self.edges if edge is not None))
         colors = plt.cm.get_cmap("viridis", len(weights))
         _, ax = plt.subplots()
         for i, weight in enumerate(weights):
             lines = [
                 [edge[-2][::-1], edge[-1][::-1]]
                 for edge in self.edges
-                if edge[2] == weight
+                if edge is not None and edge[2] == weight
             ]
             ax.add_collection(
                 LineCollection(
@@ -60,10 +67,108 @@ def display_paths(*, paths: list[list[Edge]]):
         ax = figure.add_subplot()
         for path_index, path in enumerate(paths):
             for edge_index, edge in enumerate(path):
-                ax.annotate(
-                    str(edge_index),
-                    xytext=edge[3],
-                    xy=edge[4],
-                    arrowprops=dict(arrowstyle="->", color=colors(path_index)),
-                )
+                if edge is None:
+                    continue
+                if edge is True_Edge:
+                    ax.annotate(
+                        str(edge_index),
+                        xytext=edge[3],
+                        xy=edge[4],
+                        width=1,
+                        arrowprops=dict(arrowstyle="->", color=colors(path_index)),
+                    )
+                if edge is Virtual_Edge:
+                    # TODO : more elegant way to display edges travelled twice
+                    ax.annotate(
+                        f"virtual {edge_index}",
+                        xytext=edge[3],
+                        xy=edge[4],
+                        width=2,
+                        arrowprops=dict(arrowstyle="->", color=colors(path_index)),
+                    )
         plt.show()
+
+    def compute_adjency(self):
+        vertex_adjency = [0 for vertex in self.vertices]
+        for edge in self.edges:
+            if edge is not None:
+                [v1, v2] = edge[:2]
+                vertex_adjency[v1] += 1
+                vertex_adjency[v2] += 1
+        self.vertex_adjency = vertex_adjency
+        self.vertices_odd = [
+            i for i in range(len(self.vertices)) if vertex_adjency[i] % 2 == 1
+        ]
+        self.adjency_computed = True
+        self.is_eulerian_cache = len(self.vertices_odd) == 0
+
+    def is_eulerian(self):
+        if not self.adjency_computed:
+            self.compute_adjency()
+        return self.is_eulerian_cache
+
+    def odd_vertices(self):
+        if not self.adjency_computed:
+            self.compute_adjency()
+        return self.vertices_odd
+
+    def next(self, input_edge_index: int, vertex_index: int):
+        for i, edge in enumerate(self.edges):
+            if i == input_edge_index or edge is None:
+                continue  # Actually not necessary for Hierholzer since we remove the edge before computing the next
+            if edge[0] == vertex_index:
+                return [i, edge[1]]
+            if edge[1] == vertex_index:
+                return [i, edge[0]]
+        return [-1, None]
+
+    # TODO instead of None, add a boolean to mark the edge as visited
+    # no need to clone this way
+    # prevents all the "if edge is None" in the code
+    def remove_edge(self, index):
+        if not self.adjency_computed:
+            self.compute_adjency()
+        edge = self.edges[index]
+        self.vertex_adjency[edge[0]] -= 1
+        self.vertex_adjency[edge[1]] -= 1
+        self.edges[index] = None
+
+    def add_virtual_edge(
+        self, v1_index: int, v2_index: int, weight: int, path: list[True_Edge]
+    ):
+        [v1, v2] = [self.vertices[v1_index], self.vertices[v2_index]]
+        new_edge: Virtual_Edge = (v1_index, v2_index, weight, v1, v2, path)
+        self.edges.append(new_edge)
+        self.adjency_computed = False
+
+    def has_unvisited_edges(self):
+        if not self.adjency_computed:
+            return len(self.edges) > 0
+        for e in self.edges:
+            if e is not None:
+                return True
+        return False
+
+    def is_hub(self, vertex):
+        if not self.adjency_computed:
+            self.compute_adjency()
+        return self.vertex_adjency[vertex] > 0
+
+    # Probably a standard method somewhere
+    def matrix(self):
+        n = len(self.vertices)
+        matrix = [[0 for i in range(n)] for i in range(n)]
+        for edge in self.edges:
+            if edge is not None:
+                [i, j, weight] = [edge[0], edge[1], edge[2]]
+                matrix[i][j] = weight
+                matrix[j][i] = weight
+        return csr_array(matrix)
+
+    def clone(self):
+        g = Graph(
+            [vertex for vertex in self.vertices], [edge for edge in self.real_edges]
+        )
+        g.real_edges = g.edges
+        g.edges = [edge for edge in self.edges]
+        return g

instances/farfalle.txt

@@ -0,0 +1,12 @@
+5 6
+0 0
+0 1
+1 1
+1 2
+2 2
+0 1 2
+1 2 1
+2 0 1
+2 3 1
+3 4 1
+4 2 1

instances/line.txt

@@ -0,0 +1,6 @@
+3 2
+0.5 0
+0.5 1
+0.5 2
+0 1 1
+1 2 1

instances/triangle.txt

@@ -0,0 +1,7 @@
+3 3
+0 0
+0 1
+1 1
+0 1 2
+1 2 1
+2 0 1

requirements.txt

@@ -1 +1,3 @@
-matplotlib
+matplotlib
+pytest
+scipy

solver/chinesePostman-test.py

@@ -0,0 +1,51 @@
+import datetime
+
+from chinesePostman import chinese_postman_solver
+
+from graph import Graph
+from input import parse_file
+
+
+# Refacto : have a fixture file
+def read_graph(in_file: str) -> Graph:
+    vertices, edges = parse_file(in_file)
+    print(f"#E={len(edges)}, #V={len(vertices)}")
+    return Graph(vertices, edges)
+
+
+line = read_graph("instances/line.txt")
+triangle = read_graph("instances/triangle.txt")
+hard_to_choose = read_graph("instances/hard_to_choose.txt")
+islands = read_graph("instances/islands.txt")
+paris = read_graph("instances/paris_map.txt")
+
+# End Refacto
+
+
+def integration():
+    path = chinese_postman_solver(line)
+    assert path[0] == 4  # weight
+    assert len(path[1]) == 3  # two real one virtual
+    path = chinese_postman_solver(triangle)
+    assert path[0] == 4
+    assert len(path[1]) == 3
+
+
+def performance(graph):
+    begin = datetime.datetime.now()
+    path = chinese_postman_solver(graph)
+    duration = datetime.datetime.now() - begin
+    print(f"found path of weight {path[0]} in {duration}s")
+
+
+def performances():
+    performance(hard_to_choose)
+    """ both these instances fail with indexes bound when reconstructing Dijkstra paths
+    performance(islands)
+    performance(paris)
+    """
+
+
+def tests():
+    integration()
+    performances()

solver/chinesePostman.py

@@ -0,0 +1,29 @@
+from hierholzer import hierholzer
+from naiveMatching import naive_matching
+
+from graph import Graph
+
+
+def chinese_postman_solver_generic(g: Graph, eulerian_cycle, matching):
+    if not g.is_eulerian():
+        print("no Eulerian cycle in the graph ; adding missing edges", flush=True)
+        g = matching(g)
+        print("computing Eulerian cycle on enriched graph", flush=True)
+    else:
+        print("Graph is Eulerian, computing cycle directly")
+    cycle = eulerian_cycle(g)
+    return cycle_to_path(cycle)
+
+
+def cycle_to_path(cycle):
+    # TODO
+    # if there are virtual edges, remove the one of maximum weight
+    # present the path in a readable way
+    return cycle
+
+
+def chinese_postman_solver(g: Graph):
+    # TODO
+    # Refacto to inject smallest path strategy as well (here choice to use Dijkstra is hidden)
+    # Implement and use Edmond's blossom algorithm instead of naive matching : will increase computation time but will give optimal solution
+    return chinese_postman_solver_generic(g, hierholzer, naive_matching)

solver/dijkstra-test.py

@@ -0,0 +1,29 @@
+import datetime
+
+from dijkstra import compute_dijkstra
+
+from graph import Graph
+from input import parse_file
+
+
+# Refacto : have a fixture file
+def read_graph(in_file: str) -> Graph:
+    vertices, edges = parse_file(in_file)
+    print(f"#E={len(edges)}, #V={len(vertices)}")
+    return Graph(vertices, edges)
+
+
+paris = read_graph("instances/paris_map.txt")
+
+# End Refacto
+
+
+def perf():
+    begin = datetime.datetime.now()
+    compute_dijkstra(paris)
+    duration = datetime.datetime.now() - begin
+    print(f"computed dijkstra in {duration}")
+
+
+def tests():
+    perf()

solver/dijkstra.py

@@ -0,0 +1,15 @@
+from scipy.sparse.csgraph import dijkstra
+
+from graph import Graph
+
+
+# TODO make it a class so it can be a service called by every method used for perfect-matching
+def compute_dijkstra(g: Graph):
+    graph = g.matrix()
+    dist_matrix, predecessors = dijkstra(
+        csgraph=graph,
+        directed=False,
+        return_predecessors=True,
+        indices=g.odd_vertices(),
+    )
+    return [dist_matrix, predecessors]