diff --git a/src/pgraph/PGraph.py b/src/pgraph/PGraph.py index b73a9926..b409b7a0 100644 --- a/src/pgraph/PGraph.py +++ b/src/pgraph/PGraph.py @@ -1379,6 +1379,11 @@ def iscyclic(self) -> bool: vertices' outgoing edges via :meth:`BaseVertex.neighbours` the same way :meth:`path_BFS` and friends do. + .. note:: A self-loop is a cycle of length one, so it makes either kind + of graph cyclic. Two parallel edges between the same pair of + vertices are a cycle in an undirected graph, but a repeated + arrow in a directed graph is not. + .. note:: A matrix-based test also exists in theory -- a digraph is acyclic iff its adjacency matrix is nilpotent (all eigenvalues zero) -- but that needs an O(N^3) eigendecomposition plus a @@ -1433,6 +1438,10 @@ def average_degree(self) -> float: :math:`E / N` for a directed graph where :math:`E` is the total number of edges and :math:`N` is the number of vertices. + Parallel edges and self-loops are each counted as one edge in + :math:`E`, which is consistent with :meth:`degree`: a self-loop + contributes 2 to the degree sum of an undirected graph. + .. runblock:: pycon >>> from pgraph import UGraph @@ -1498,9 +1507,16 @@ def Laplacian(self) -> NDArray: hold. For example a simple weakly-connected out-tree (one component) already has two zero eigenvalues, not one. + .. rubric:: For graphs with parallel edges or loops + + Both :meth:`degree` and :meth:`adjacency` count edges, so the + identity ``L = degree() - adjacency()`` and the zero row sums hold + unchanged. A self-loop adds the same amount to both and cancels out + of the Laplacian entirely. + :seealso: :meth:`adjacency` :meth:`incidence` :meth:`degree` """ - return self.degree() - (self.adjacency() > 0) + return self.degree() - self.adjacency() def connectivity(self, vertices: Iterable[BaseVertex] | None = None) -> list[int]: """ @@ -1535,6 +1551,10 @@ def connectivity(self, vertices: Iterable[BaseVertex] | None = None) -> list[int >>> print(c) + .. note:: This is the diagonal of :meth:`degree`, so parallel edges are + counted individually and a self-loop counts twice for a + ``UGraph`` and once for a ``DGraph``. See the discussion there. + :seealso: :meth:`degree` """ @@ -1572,6 +1592,13 @@ def degree(self) -> NDArray: >>> d = g.degree() >>> print(d) + .. rubric:: For graphs with parallel edges or loops + + Every edge is counted, so ``k`` parallel edges between two vertices + add ``k`` to the degree of each end. A self-loop adds 2 to the degree + of a ``UGraph`` vertex, since it has two ends, and 1 to the + (out-)degree of a ``DGraph`` vertex. + :seealso: :meth:`adjacency` :meth:`incidence` :meth:`laplacian` """ @@ -1607,11 +1634,20 @@ def adjacency(self) -> NDArray: can be resolved to a vertex reference by ``graph[i]``. - for an undirected graph the matrix is symmetric - Eigenvalues of ``A`` are real and are known as the spectrum of the graph. - - The element ``A[i,j]`` can be considered the number of walks of length one - edge from vertex ``i`` to vertex ``j`` (either zero or one). + - The element ``A[i,j]`` is the number of walks of length one + edge from vertex ``i`` to vertex ``j``, which is the number of + edges between them (zero or one for a simple graph). - If ``Ak = A ** k`` the element ``Ak[i,j]`` is the number of walks of length ``k`` from vertex ``i`` to vertex ``j``. + .. rubric:: For graphs with parallel edges or loops + + The matrix holds edge *counts*, so ``k`` parallel edges from ``i`` to + ``j`` give ``A[i,j] = k`` and ``Ak = A ** k`` still counts walks + correctly. A self-loop adds to the diagonal: 2 for a ``UGraph`` + (consistent with :meth:`degree`) and 1 for a ``DGraph``. In every case + each row of ``A`` sums to the degree of that vertex. + :seealso: :meth:`Laplacian` :meth:`incidence` :meth:`degree` """ # create a dict mapping vertex to an id @@ -1622,7 +1658,7 @@ def adjacency(self) -> NDArray: A = np.zeros((self.n, self.n)) for vertex in self: for n in vertex.neighbours(): - A[vdict[vertex], vdict[n]] = 1 + A[vdict[vertex], vdict[n]] += 1 return A def incidence(self) -> NDArray: @@ -1661,6 +1697,13 @@ def incidence(self) -> NDArray: :meth:`BaseVertex.edges` instead would silently drop such vertices for a directed graph. + .. rubric:: For graphs with parallel edges or loops + + Parallel edges each get their own column. For a ``UGraph`` the column + of a self-loop holds a single 2 (it touches its vertex at both ends, + so row sums still equal :meth:`degree`); for a ``DGraph`` it holds a + single 1. + :seealso: :meth:`Laplacian` :meth:`adjacency` :meth:`degree` """ edges = self.edges() @@ -1674,6 +1717,8 @@ def incidence(self) -> NDArray: assert e.v1 is not None and e.v2 is not None I[vdict[e.v1], j] = 1 I[vdict[e.v2], j] = 1 + if e.v1 is e.v2 and isinstance(self, UGraph): + I[vdict[e.v1], j] = 2 # self-loop, counted at both ends return I @@ -1702,6 +1747,12 @@ def distance(self) -> NDArray: >>> d = g.distance() >>> print(d) + .. rubric:: For graphs with parallel edges or loops + + If there are parallel edges the element is the cost of the *cheapest* + of them. A self-loop gives a non-zero diagonal element equal to its + cost. + :seealso: :meth:`BaseVertex.distance` """ # create a dict mapping vertex to an id @@ -1709,10 +1760,12 @@ def distance(self) -> NDArray: for i, vert in enumerate(self): vdict[vert] = i - A = np.zeros((self.n, self.n)) + A = np.full((self.n, self.n), np.inf) for v1 in self: for v2, edge in v1.incidences(): - A[vdict[v1], vdict[v2]] = self._require_cost(edge) + i, j = vdict[v1], vdict[v2] + A[i, j] = min(A[i, j], self._require_cost(edge)) + A[np.isinf(A)] = 0 # not connected return A # GRAPH COMPONENTS @@ -1789,6 +1842,9 @@ def path_BFS( :type G: BaseVertex subclass :param verbose: print search progress, defaults to False :param summary: print a one-line search summary, defaults to False + :raises TypeError: ``S`` or ``G`` is neither a vertex nor a string + :raises KeyError: ``S`` or ``G`` is a name not in the graph + :raises ValueError: an edge on the path has no cost :return: list of vertices from S to G inclusive, path length :rtype: list of BaseVertex subclass, float @@ -1808,6 +1864,9 @@ def path_BFS( >>> print(path) >>> print(length) + .. note:: Self-loops are ignored. If there are parallel edges the + cheapest is used for the path length. + :seealso: :meth:`path_UCS` :meth:`path_Astar` """ S = self._resolve_vertex(S, "start") @@ -1832,6 +1891,8 @@ def path_BFS( # expand the vertex for n in x.neighbours(): + if n is x: + continue # self-loop, never part of a shortest path if n is G: if verbose: print(" goal", n.name, "reached") @@ -1860,7 +1921,7 @@ def path_BFS( while x is not S: p = parent[x] - length += self._require_cost(x.edgeto(p)) + length += self._require_cost(p.edgeto(x)) path.insert(0, p) x = p @@ -1883,6 +1944,9 @@ def path_UCS( :type G: BaseVertex subclass :param verbose: print search progress, defaults to False :param summary: print a one-line search summary, defaults to False + :raises TypeError: ``S`` or ``G`` is neither a vertex nor a string + :raises KeyError: ``S`` or ``G`` is a name not in the graph + :raises ValueError: an edge examined during the search has no cost :return: list of vertices from S to G inclusive, path length, tree :rtype: list of BaseVertex subclass, float, dict @@ -1907,6 +1971,9 @@ def path_UCS( >>> print(path) >>> print(length) + .. note:: Self-loops are ignored. If there are parallel edges the + cheapest is used for the path length. + :seealso: :meth:`path_BFS` :meth:`path_Astar` """ S = self._resolve_vertex(S, "start") @@ -1933,6 +2000,8 @@ def path_UCS( break # expand the vertex for n, e in x.incidences(): + if n is x: + continue # self-loop, never part of a shortest path fnew = f[x] + self._require_cost(e) if n not in frontier and n not in explored: # add it to the frontier @@ -1995,6 +2064,9 @@ def path_Astar( :type G: BaseVertex subclass :param verbose: print search progress, defaults to False :param summary: print a one-line search summary, defaults to False + :raises TypeError: ``S`` or ``G`` is neither a vertex nor a string + :raises KeyError: ``S`` or ``G`` is a name not in the graph + :raises ValueError: an edge examined during the search has no cost :return: list of vertices from S to G inclusive, path length, tree :rtype: list of BaseVertex subclass, float, dict @@ -2019,6 +2091,9 @@ def path_Astar( >>> print(path) >>> print(length) + .. note:: Self-loops are ignored. If there are parallel edges the + cheapest is used for the path length. + :seealso: :meth:`heuristic` :meth:`path_BFS` :meth:`path_UCS` """ S = self._resolve_vertex(S, "start") @@ -2046,6 +2121,8 @@ def path_Astar( break # expand the vertex for n, e in x.incidences(): + if n is x: + continue # self-loop, never part of a shortest path if n not in frontier and n not in explored: # add it to the frontier frontier.append(n) @@ -2646,6 +2723,10 @@ def neighbours(self) -> list[BaseVertex]: .. note:: For a directed graph the neighbours are those on edges leaving this vertex + .. note:: There is one entry per edge, so a neighbour joined by ``k`` + parallel edges appears ``k`` times. In an undirected graph a + self-loop lists the vertex itself twice, once for each end. + .. runblock:: pycon >>> from pgraph import UGraph @@ -2864,6 +2945,8 @@ def edgeto(self, dest: BaseVertex) -> Edge: .. note:: - For a directed graph ``dest`` must be at the arrow end of the edge + - If there are parallel edges to ``dest`` the cheapest is returned + (the first, if none of them has a cost) .. runblock:: pycon @@ -2874,10 +2957,13 @@ def edgeto(self, dest: BaseVertex) -> Edge: >>> g.add_edge(v1, v2, cost=1.414) >>> print(v1.edgeto(v2)) """ - for n, e in self.incidences(): - if n is dest: - return e - raise ValueError("dest is not a neighbour") + edges = [e for n, e in self.incidences() if n is dest] + if not edges: + raise ValueError("dest is not a neighbour") + costed = [e for e in edges if e.cost is not None] + if costed: + return min(costed, key=lambda e: e.cost) # type: ignore[arg-type,return-value] + return edges[0] def edges(self) -> list[Edge]: """ @@ -2985,6 +3071,12 @@ def degree(self) -> int: >>> g.add_edge(v1, v3) >>> print(v1.degree) + .. rubric:: For graphs with parallel edges or loops + + Every edge is counted, so ``k`` parallel edges to the same neighbour + add ``k``. A self-loop adds 2 for a ``UGraph`` vertex, since it has two + ends, and 1 for a ``DGraph`` vertex. + :seealso: :meth:`edges` """ return len(self.edges()) diff --git a/tests/test_multigraph.py b/tests/test_multigraph.py new file mode 100644 index 00000000..a309b6de --- /dev/null +++ b/tests/test_multigraph.py @@ -0,0 +1,296 @@ +""" +Self-loops and parallel edges. + +pgraph does not forbid either, so every quantity derived from the edge set +(degree, the matrix representations, cycle detection, path planning) has to +give a defined, consistent answer when they are present. + +The same narrative is run against ``UGraph`` and ``DGraph``: + +1. a simple acyclic network, the baseline +2. add a self-loop on a vertex that lies on the planned route +3. remove it, which must restore the baseline exactly +4. add parallel edges, planning must take the cheapest +5. add the self-loop back in on top of the parallel edges + +Conventions being tested (see the matrix methods' docstrings): + +- a self-loop counts twice towards the degree of a ``UGraph`` vertex, once + towards the (out-)degree of a ``DGraph`` vertex +- ``adjacency()`` counts edges, so ``A[i, j]`` is the number of edges + from ``i`` to ``j`` and each row of ``A`` sums to the degree of the vertex +- ``Laplacian()`` is therefore ``degree() - adjacency()``, rows sum to zero +- ``distance()`` and path lengths use the cheapest of any parallel edges +""" + +import contextlib +import signal +import unittest + +import numpy as np + +from pgraph import DGraph, UGraph + + +@contextlib.contextmanager +def within(seconds=5): + """Fail, rather than hang the whole suite, if the body does not return.""" + if not hasattr(signal, "SIGALRM"): # pragma: no cover - not on Windows + yield + return + + def handler(signum, frame): + raise AssertionError(f"did not finish within {seconds}s (infinite loop?)") + + old = signal.signal(signal.SIGALRM, handler) + signal.alarm(seconds) + try: + yield + finally: + signal.alarm(0) + signal.signal(signal.SIGALRM, old) + + +class _Multigraph: + """ + Shared narrative; subclasses set ``Graph`` and the expected values that + differ between undirected and directed graphs. + + The network is the line A--B--C--D with unit costs. Vertices are 0.5 + apart so the default Euclidean heuristic stays admissible for A*. + """ + + Graph: type + # connectivity() / degree-matrix diagonal, in vertex creation order + deg_base: list + deg_loop: list + deg_parallel: list + deg_both: list + # incidence column sum of a self-loop (2 undirected, 1 directed) + loop_colsum: int + # adjacency A[B, B] for a single self-loop on B + loop_adjacency: int + directed: bool + + # ------------------------------------------------------------------ # + + def setUp(self): + g = self.Graph() + self.g = g + self.A = g.add_vertex([0.0, 0], name="A") + self.B = g.add_vertex([0.5, 0], name="B") + self.C = g.add_vertex([1.0, 0], name="C") + self.D = g.add_vertex([1.5, 0], name="D") + self.AB = self.A.connect(self.B, cost=1.0) + self.BC = self.B.connect(self.C, cost=1.0) + self.CD = self.C.connect(self.D, cost=1.0) + + def add_loop(self, at=None): + v = self.B if at is None else at + return v.connect(v, cost=2.0) + + def add_parallels(self): + # AB already exists with cost 1.0; the cheap one is added last on + # purpose, so "first found" and "last written" both give a wrong answer + self.A.connect(self.B, cost=3.0) + self.A.connect(self.B, cost=0.5) + + # ------------------------------------------------------------------ # + # generic checks + + def check_degrees(self, expected): + g = self.g + self.assertEqual(g.connectivity(), expected) + np.testing.assert_array_equal(np.diag(g.degree()), expected) + self.assertAlmostEqual(g.average_degree(), sum(expected) / g.n) + + def check_matrices_consistent(self): + """Relationships that must hold for *any* edge multiset.""" + g = self.g + deg = np.array(g.connectivity()) + + A = g.adjacency() + np.testing.assert_array_equal(A.sum(axis=1), deg) + if not self.directed: + np.testing.assert_array_equal(A, A.T) + + L = g.Laplacian() + np.testing.assert_array_equal(L, g.degree() - A) + np.testing.assert_allclose(L.sum(axis=1), 0) + + I = g.incidence() + self.assertEqual(I.shape, (g.n, g.ne)) + if not self.directed: + # an undirected incidence row sums to the vertex degree + np.testing.assert_array_equal(I.sum(axis=1), deg) + + def check_baseline(self): + g = self.g + self.assertEqual((g.n, g.ne, g.nc), (4, 3, 1)) + self.assertFalse(g.iscyclic()) + self.check_degrees(self.deg_base) + self.check_matrices_consistent() + + A = g.adjacency() + self.assertEqual(A.sum(), 6 if not self.directed else 3) + self.assertEqual(A.max(), 1) + self.assertEqual(np.trace(A), 0) + np.testing.assert_array_equal(g.incidence().sum(axis=0), [2, 2, 2]) + + D = g.distance() + self.assertEqual(D[0, 1], 1.0) + self.assertEqual(D[1, 2], 1.0) + self.assertEqual(D[2, 3], 1.0) + + self.check_planning(length=2.0) + + def check_planning(self, length): + """Every planner must reach C from A via B, with this length.""" + g, A, C = self.g, self.A, self.C + names = lambda path: [v.name for v in path] + + with within(): + path, cost, _ = g.path_UCS(A, C) + self.assertEqual(names(path), ["A", "B", "C"]) + self.assertAlmostEqual(cost, length) + + with within(): + path, cost, _ = g.path_Astar(A, C) + self.assertEqual(names(path), ["A", "B", "C"]) + self.assertAlmostEqual(cost, length) + + with within(): + path, cost = g.path_BFS(A, C) + self.assertEqual(names(path), ["A", "B", "C"]) + self.assertAlmostEqual(cost, length) + + # ------------------------------------------------------------------ # + # the narrative + + def test_baseline(self): + self.check_baseline() + + def test_self_loop(self): + g = self.g + loop = self.add_loop() + + self.assertEqual((g.n, g.ne, g.nc), (4, 4, 1)) + self.assertTrue(g.iscyclic()) + self.check_degrees(self.deg_loop) + self.check_matrices_consistent() + + A = g.adjacency() + self.assertEqual(A[1, 1], self.loop_adjacency) + self.assertEqual(np.trace(A), self.loop_adjacency) + + I = g.incidence() + self.assertEqual(sorted(I.sum(axis=0)), sorted([2, 2, 2, self.loop_colsum])) + + D = g.distance() + self.assertEqual(D[1, 1], 2.0) # cost of the loop + self.assertEqual(D[0, 1], 1.0) # others unchanged + + # the loop lies on the route A -> B -> C, planning must ignore it + self.check_planning(length=2.0) + + # the loop is a real edge of the graph + self.assertIn(loop, g.edges()) + + def test_self_loop_removed_restores_baseline(self): + loop = self.add_loop() + self.g.remove_edge(loop) + self.check_baseline() + + def test_self_loop_off_route(self): + # loops at the start and beyond the goal never mattered, guard them too + self.add_loop(self.A) + self.add_loop(self.D) + self.assertTrue(self.g.iscyclic()) + self.check_matrices_consistent() + self.check_planning(length=2.0) + + def test_parallel_edges(self): + g = self.g + self.add_parallels() + + self.assertEqual((g.n, g.ne, g.nc), (4, 5, 1)) + # parallel edges are a cycle in an undirected graph, but a repeated + # arrow in a directed graph is not + self.assertEqual(g.iscyclic(), not self.directed) + self.check_degrees(self.deg_parallel) + self.check_matrices_consistent() + + A = g.adjacency() + self.assertEqual(A[0, 1], 3) # three A->B edges, not "connected: 1" + self.assertEqual(A[1, 2], 1) + + D = g.distance() + self.assertEqual(D[0, 1], 0.5) # the cheapest, whatever the order + + # edgeto() is what the planners use to turn a path into a length + self.assertEqual(self.A.edgeto(self.B).cost, 0.5) + + # route A -> B -> C must use the 0.5 edge: 0.5 + 1.0 + self.check_planning(length=1.5) + + def test_parallel_edges_removed_restores_baseline(self): + before = set(self.g.edges()) + self.add_parallels() + for e in set(self.g.edges()) - before: + self.g.remove_edge(e) + self.check_baseline() + + def test_parallel_edges_and_self_loop(self): + g = self.g + self.add_parallels() + self.add_loop() + + self.assertEqual((g.n, g.ne, g.nc), (4, 6, 1)) + self.assertTrue(g.iscyclic()) + self.check_degrees(self.deg_both) + self.check_matrices_consistent() + + A = g.adjacency() + self.assertEqual(A[0, 1], 3) + self.assertEqual(A[1, 1], self.loop_adjacency) + + D = g.distance() + self.assertEqual(D[0, 1], 0.5) + self.assertEqual(D[1, 1], 2.0) + + self.check_planning(length=1.5) + + +class TestUGraphMultigraph(_Multigraph, unittest.TestCase): + Graph = UGraph + directed = False + deg_base = [1, 2, 2, 1] + deg_loop = [1, 4, 2, 1] # a loop contributes 2 + deg_parallel = [3, 4, 2, 1] + deg_both = [3, 6, 2, 1] + loop_colsum = 2 + loop_adjacency = 2 + + def test_parallel_edges_bfs(self): + # BFS minimises hops, but must still report the cheap edge's cost + self.add_parallels() + self.add_loop() + with within(): + path, cost = self.g.path_BFS(self.A, self.C) + self.assertEqual([v.name for v in path], ["A", "B", "C"]) + self.assertAlmostEqual(cost, 1.5) + + +class TestDGraphMultigraph(_Multigraph, unittest.TestCase): + Graph = DGraph + directed = True + deg_base = [1, 1, 1, 0] # out-degree only + deg_loop = [1, 2, 1, 0] # a loop contributes 1 + deg_parallel = [3, 1, 1, 0] + deg_both = [3, 2, 1, 0] + loop_colsum = 1 + loop_adjacency = 1 + + +if __name__ == "__main__": + unittest.main()