Added new classes to return the edges of polyhedra (#171)

* test: Fix assumption on quaternions.

Only test quaternions that are not effectively zero.

* Added new classes to return the edges of polyhedra

The new class ``polyhedron.edges`` returns the a list of the edges of a polyhedron as vertex-index pairs, similar to the current ``polyhedron.faces`` class. The class ``polyhedron.get_edge_vectors`` returns a list of edges as vectors in 3D space.

* Fixed return for get_edge_vectors method.

* Renamed get_edge_vectors property to edge_vectors

Co-authored-by: Vyas Ramasubramani <[email protected]>

* Vectorized edge_vectors return with numpy

Co-authored-by: Vyas Ramasubramani <[email protected]>

* Updated test_polyhedron to be compatible with the renamed edge_vectors property

* Updated test_polyhedron to explicitly cover the edges property

* Updated rtol for edge property test

Until the precision improvements in #177 are added, a decrease in rtol and the merge_faces call are required for pytest to succeed on the dodecahedron's edges. Once the precision is increased, these temporary solutions can be reverted

* Reverted small fixes that account for inaccuracies in polyhedron stored data

Currently, the stored data for polyhedra is not accurate enough for the ``edges`` and ``edge_vectors`` properties to function as expected. Once #177 is merged, the accuracy increase should fix the issue and assertion tolerances for testing can be tightened. In addition, ``merge_faces`` should no longer be required for any of the polyhedra included with this package.

* Removed unnecessary comment from ``test_edges``

* Changed ``edge_vectors`` property to return a numpy array

Co-authored-by: Bradley Dice <[email protected]>

* Rewrote ``edges`` method and updated ``edge_vectors`` for compatibility

After discussing possible options for returning polyhedral edges, it was determined that a set of i<j pairs of vertices was the best option for output in the ``edges`` method. A description of how to generate (j,i) pairs was added to the docstring, and the edge_vectors method was updated to work with the new set/tuple structure.

* Updated ``test_polyhedron`` to work with set edges

* Fixed docstring formatting for ``edges`` and ``edge_vectors``

* Updated edges method to return Numpy array

It was determined that edges should be returned as an ordered Numpy array. This commit ensures the final output is a numpy array, sorted first by the (i) element and then by the (j) element where i<j. This mimics Mathematica's edges output format and should aid in usability. Documentation was updated accordingly.

* test_polyhedron::test_edges now checks edge count

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Updated edges property to cache result

* Updated edges and edge_vectors to make better use of numpy functions and functools caching

* Added num_edges method

* Updated credits

* Updated test_polyhedron to test num_edges property

* Updated edges documentation

Co-authored-by: Bradley Dice <[email protected]>

* Updated edge_vectors documentation

Co-authored-by: Bradley Dice <[email protected]>

* Refactored test_edges to be more comprehensive

Added explicit tests for sorting, double checks for edge length, and an additional test for the number of edges based on the euler characteristic

* Added fast num_edges calculation for convex polyhedron and improved pytests

* test_num_edges now covers nonconvex Polyhedron class

* Update Credits.rst

Co-authored-by: Bradley Dice <[email protected]>

* Test i-i edges

Co-authored-by: Bradley Dice <[email protected]>

* Removed changes to .gitignore

---------

Co-authored-by: Brandon Butler <[email protected]>
Co-authored-by: Vyas Ramasubramani <[email protected]>
Co-authored-by: Bradley Dice <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Credits.rst

@@ -112,6 +112,7 @@ Jen Bradley
 * Added shape families for Archimedean, Catalan, and Johnson solids.
 * Added shape family for prisms and antiprisms.
 * Added shape family for equilateral pyramids and dipyramids.
+* Added edges, edge_vectors, and num_edges methods.
 
 Domagoj Fijan
 * Rewrote point in polygon check to use NumPy vectorized operations.

coxeter/shapes/convex_polyhedron.py

@@ -126,6 +126,12 @@ def asphericity(self):
         """float: Get the asphericity as defined in :cite:`Irrgang2017`."""
         return self.mean_curvature * self.surface_area / (3 * self.volume)
 
+    @property
+    def num_edges(self):
+        """int: Get the number of edges."""
+        # Calculate number of edges from Euler Characteristic
+        return self.num_vertices + self.num_faces - 2
+
     def is_inside(self, points):
         """Determine whether points are contained in this polyhedron.
 

coxeter/shapes/polyhedron.py

@@ -4,6 +4,7 @@
 """Defines a polyhedron."""
 
 import warnings
+from functools import cached_property
 
 import numpy as np
 import rowan
@@ -356,9 +357,45 @@ def vertices(self):
 
     @property
     def faces(self):
-        """list(:class:`numpy.ndarray`): Get the polyhedron's faces."""
+        """list(:class:`numpy.ndarray`): Get the polyhedron's faces.
+
+        Results returned as vertex index lists.
+        """
         return self._faces
 
+    @cached_property
+    def edges(self):
+        """:class:`numpy.ndarray`: Get the polyhedron's edges.
+
+        Results returned as vertex index pairs,  with each edge of the polyhedron
+        included exactly once.  Edge (i,j) pairs are ordered by vertex index with i<j.
+        """
+        ij_pairs = np.array(
+            [
+                [i, j]
+                for face in self.faces
+                for i, j in zip(face, np.roll(face, -1))
+                if i < j
+            ]
+        )
+        sorted_indices = np.lexsort(ij_pairs.T[::-1])
+        sorted_ij_pairs = ij_pairs[sorted_indices]
+        # Make edge data read-only so that the cached property of this instance
+        # cannot be edited
+        sorted_ij_pairs.flags.writeable = False
+
+        return sorted_ij_pairs
+
+    @property
+    def edge_vectors(self):
+        """:class:`numpy.ndarray`: Get the polyhedron's edges as vectors."""
+        return self.vertices[self.edges[:, 1]] - self.vertices[self.edges[:, 0]]
+
+    @property
+    def num_edges(self):
+        """int: Get the number of edges."""
+        return len(self.edges)
+
     @property
     def volume(self):
         """float: Get or set the polyhedron's volume."""

tests/test_polyhedron.py

@@ -28,7 +28,7 @@
     named_pyramiddipyramid_mark,
     sphere_isclose,
 )
-from coxeter.families import DOI_SHAPE_REPOSITORIES, PlatonicFamily
+from coxeter.families import DOI_SHAPE_REPOSITORIES, ArchimedeanFamily, PlatonicFamily
 from coxeter.shapes import ConvexPolyhedron, Polyhedron
 from coxeter.shapes.utils import rotate_order2_tensor, translate_inertia_tensor
 from utils import compute_centroid_mc, compute_inertia_mc
@@ -339,6 +339,96 @@ def test___repr__():
     repr(icosidodecahedron)
 
 
+@combine_marks(
+    named_platonic_mark,
+    named_archimedean_mark,
+    named_catalan_mark,
+    named_johnson_mark,
+    named_prismantiprism_mark,
+    named_pyramiddipyramid_mark,
+)
+def test_edges(poly):
+    # Check that the first column is in ascending order.
+    assert np.all(np.diff(poly.edges[:, 0]) >= 0)
+
+    # Check that all items in the first column are greater than those in the second.
+    assert np.all(np.diff(poly.edges, axis=1) > 0)
+
+    # Check the second column is in ascending order for each unique item in the first.
+    # For example, [[0,1],[0,3],[1,2]] is permitted but [[0,1],[0,3],[0,2]] is not.
+    edges = poly.edges
+    unique_values = unique_values = np.unique(edges[:, 0])
+    assert all(
+        [
+            np.all(np.diff(edges[edges[:, 0] == value, 1]) >= 0)
+            for value in unique_values
+        ]
+    )
+
+    # Check that there are no duplicate edges. This also double-checks the sorting
+    assert np.all(np.unique(poly.edges, axis=1) == poly.edges)
+
+    # Check that the edges are immutable
+    try:
+        poly.edges[1] = [99, 99]
+        # If the assignment works, catch that:
+        assert poly.edges[1] != [99, 99]
+    except ValueError as ve:
+        assert "read-only" in str(ve)
+
+
+def test_edge_lengths():
+    known_shapes = {
+        "Tetrahedron": np.sqrt(2) * np.cbrt(3),
+        "Cube": 1,
+        "Octahedron": np.power(2, 5 / 6) * np.cbrt(3 / 8),
+        "Dodecahedron": np.power(2, 2 / 3) * np.cbrt(1 / (15 + np.sqrt(245))),
+        "Icosahedron": np.cbrt(9 / 5 - 3 / 5 * np.sqrt(5)),
+    }
+    for name, edgelength in known_shapes.items():
+        poly = PlatonicFamily.get_shape(name)
+        # Check that edge lengths are correct
+        veclens = np.linalg.norm(
+            poly.vertices[poly.edges[:, 1]] - poly.vertices[poly.edges[:, 0]], axis=1
+        )
+        assert np.allclose(veclens, edgelength)
+        assert np.allclose(veclens, np.linalg.norm(poly.edge_vectors, axis=1))
+
+
+def test_num_edges_archimedean():
+    known_shapes = {
+        "Cuboctahedron": 24,
+        "Icosidodecahedron": 60,
+        "Truncated Tetrahedron": 18,
+        "Truncated Octahedron": 36,
+        "Truncated Cube": 36,
+        "Truncated Icosahedron": 90,
+        "Truncated Dodecahedron": 90,
+        "Rhombicuboctahedron": 48,
+        "Rhombicosidodecahedron": 120,
+        "Truncated Cuboctahedron": 72,
+        "Truncated Icosidodecahedron": 180,
+        "Snub Cuboctahedron": 60,
+        "Snub Icosidodecahedron": 150,
+    }
+    for name, num_edges in known_shapes.items():
+        poly = ArchimedeanFamily.get_shape(name)
+        assert poly.num_edges == num_edges
+
+
+@given(
+    EllipsoidSurfaceStrategy,
+)
+def test_num_edges_polyhedron(points):
+    hull = ConvexHull(points)
+    poly = ConvexPolyhedron(points[hull.vertices])
+    ppoly = Polyhedron(poly.vertices, poly.faces)
+
+    # Calculate correct number of edges from euler characteristic
+    euler_characteristic_edge_count = ppoly.num_vertices + ppoly.num_faces - 2
+    assert ppoly.num_edges == euler_characteristic_edge_count
+
+
 def test_curvature():
     """Regression test against values computed with older method."""
     # The shapes in the PlatonicFamily are normalized to unit volume.