From f0b8598e445f50fdbf9bbdfc887b9b7c6cb4a08d Mon Sep 17 00:00:00 2001 From: Brandon Butler Date: Wed, 25 Jan 2023 18:30:33 -0500 Subject: [PATCH 01/73] test: Fix assumption on quaternions. Only test quaternions that are not effectively zero. --- tests/test_polygon.py | 4 ++-- tests/test_spheropolygon.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_polygon.py b/tests/test_polygon.py index 4c6d78be..9973f785 100644 --- a/tests/test_polygon.py +++ b/tests/test_polygon.py @@ -189,7 +189,7 @@ def test_convex_area(points): @given(random_quat=arrays(np.float64, (4,), elements=floats(-1, 1, width=64))) def test_rotation_signed_area(random_quat): """Ensure that rotating does not change the signed area.""" - assume(not np.all(random_quat == 0)) + assume(not np.allclose(random_quat, 0)) random_quat = rowan.normalize(random_quat) rotated_points = rowan.rotate(random_quat, get_square_points()) poly = Polygon(rotated_points) @@ -253,7 +253,7 @@ def test_bounding_circle_radius_random_hull(points): ) def test_bounding_circle_radius_random_hull_rotation(points, rotation): """Test that rotating vertices does not change the bounding radius.""" - assume(not np.all(rotation == 0)) + assume(not np.allclose(rotation, 0)) hull = ConvexHull(points) poly = Polygon(points[hull.vertices]) diff --git a/tests/test_spheropolygon.py b/tests/test_spheropolygon.py index a9df4a0e..a1d34b6e 100644 --- a/tests/test_spheropolygon.py +++ b/tests/test_spheropolygon.py @@ -152,7 +152,7 @@ def test_convex_signed_area(square_points): ) ) def testfun(random_quat): - assume(not np.all(random_quat == 0)) + assume(not np.allclose(random_quat, 0)) random_quat = rowan.normalize(random_quat) rotated_points = rowan.rotate(random_quat, square_points) r = 1 From ab212791e2f016a71a4b6e582813326527d28d10 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 6 Feb 2023 11:01:17 -0500 Subject: [PATCH 02/73] 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. --- .gitignore | 3 ++- coxeter/shapes/polyhedron.py | 28 +++++++++++++++++++++++++++- tests/test_polyhedron.py | 18 ++++++++++++++++++ 3 files changed, 47 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index 028cfa56..22b6b6b8 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,6 @@ *~ *.swp - +.DS_Store # Packages *.egg *.egg-info* @@ -17,6 +17,7 @@ lib lib64 __pycache__ .hypothesis +.coverage *.pyc doc/source/bibtex.json .ipynb_checkpoints diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 074db4b6..13481c42 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -356,9 +356,35 @@ 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 + @property + def edges(self): + """list(:class:`numpy.ndarray`): Get the polyhedron's edges. + + Results returned as vertex index pairs. + """ + edges = [] + for face in self.faces: + [ + edges.append((i, j)) + for i, j in zip(np.append(face, face[0]), np.append(face[1:], face[0])) + if (j, i) not in edges + ] + return edges + + @property + def get_edge_vectors(self): + """list(:class:`numpy.ndarray`): Get the polyhedron's edges as vectors.""" + vecs = [] + vertices = self.vertices + for edge in self.edges: + vecs.append(vertices[edge[1]] - vertices[edge[0]]) + @property def volume(self): """float: Get or set the polyhedron's volume.""" diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 0f3286c7..017a6208 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -224,6 +224,24 @@ def test_dihedrals(): assert np.isclose(poly.get_dihedral(i, j), dihedral) +def test_edges(): + # The shapes in the PlatonicFamily are normalized to unit volume + 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) + if name == "Dodecahedron": + poly.merge_faces(rtol=1) + for edge in poly.get_edge_vectors: + # rtol must be lowered to accomodate small inaccuracies in vertex locations + assert np.isclose(np.linalg.norm(edge), edgelength, rtol=1e-4) + + def test_curvature(): """Regression test against values computed with older method.""" # The shapes in the PlatonicFamily are normalized to unit volume. From 22399d8cfc3912dc9ad6fdcda4e2f3588f43293c Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 7 Feb 2023 14:11:05 -0500 Subject: [PATCH 03/73] Fixed return for get_edge_vectors method. --- coxeter/shapes/polyhedron.py | 1 + 1 file changed, 1 insertion(+) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 13481c42..f8dcda73 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -384,6 +384,7 @@ def get_edge_vectors(self): vertices = self.vertices for edge in self.edges: vecs.append(vertices[edge[1]] - vertices[edge[0]]) + return vecs @property def volume(self): From 7c941004342f69b80d53fe4514fb80177899ee09 Mon Sep 17 00:00:00 2001 From: Jen Bradley <55467578+janbridley@users.noreply.github.com> Date: Fri, 24 Feb 2023 13:50:58 -0500 Subject: [PATCH 04/73] Renamed get_edge_vectors property to edge_vectors Co-authored-by: Vyas Ramasubramani --- coxeter/shapes/polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index f8dcda73..42dd8c80 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -378,7 +378,7 @@ def edges(self): return edges @property - def get_edge_vectors(self): + def edge_vectors(self): """list(:class:`numpy.ndarray`): Get the polyhedron's edges as vectors.""" vecs = [] vertices = self.vertices From 0e6501a19c23980dfb98b911b3a1116d42b9ec3c Mon Sep 17 00:00:00 2001 From: Jen Bradley <55467578+janbridley@users.noreply.github.com> Date: Mon, 27 Feb 2023 11:03:58 -0500 Subject: [PATCH 05/73] Vectorized edge_vectors return with numpy Co-authored-by: Vyas Ramasubramani --- coxeter/shapes/polyhedron.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 42dd8c80..87e7146a 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -380,11 +380,8 @@ def edges(self): @property def edge_vectors(self): """list(:class:`numpy.ndarray`): Get the polyhedron's edges as vectors.""" - vecs = [] - vertices = self.vertices - for edge in self.edges: - vecs.append(vertices[edge[1]] - vertices[edge[0]]) - return vecs + edges = np.asarray(self.edges) + return list(self.vertices[edges[:, 1]] - self.vertices[edges[:, 0]]) @property def volume(self): From fdec57c7e76d028296487a7a28382ea702c80681 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 27 Feb 2023 11:12:38 -0500 Subject: [PATCH 06/73] Updated test_polyhedron to be compatible with the renamed edge_vectors property --- tests/test_polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 621f841f..dddda5b5 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -238,7 +238,7 @@ def test_edges(): poly = PlatonicFamily.get_shape(name) if name == "Dodecahedron": poly.merge_faces(rtol=1) - for edge in poly.get_edge_vectors: + for edge in poly.edge_vectors: # rtol must be lowered to accomodate small inaccuracies in vertex locations assert np.isclose(np.linalg.norm(edge), edgelength, rtol=1e-4) From dd47ed2cb8c5416e8e5281a91dc5eeaa96ba0cff Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 27 Feb 2023 11:24:49 -0500 Subject: [PATCH 07/73] Updated test_polyhedron to explicitly cover the edges property --- tests/test_polyhedron.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index dddda5b5..750c5811 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -238,10 +238,17 @@ def test_edges(): poly = PlatonicFamily.get_shape(name) if name == "Dodecahedron": poly.merge_faces(rtol=1) + # Test edge_vectors property for edge in poly.edge_vectors: # rtol must be lowered to accomodate small inaccuracies in vertex locations assert np.isclose(np.linalg.norm(edge), edgelength, rtol=1e-4) + # Test edges property + edges = np.asarray(poly.edges) + vertices = poly.vertices + veclens = np.linalg.norm(vertices[edges[:, 1]] - vertices[edges[:, 0]], axis=1) + assert np.allclose(veclens, edgelength) + def test_curvature(): """Regression test against values computed with older method.""" From 210f0099e04ce5b69ad77d53c458ea18b6333abf Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 27 Feb 2023 13:05:44 -0500 Subject: [PATCH 08/73] 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 --- tests/test_polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 750c5811..1c9c0469 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -247,7 +247,7 @@ def test_edges(): edges = np.asarray(poly.edges) vertices = poly.vertices veclens = np.linalg.norm(vertices[edges[:, 1]] - vertices[edges[:, 0]], axis=1) - assert np.allclose(veclens, edgelength) + assert np.allclose(veclens, edgelength, rtol=1e-4) def test_curvature(): From a2cfd3ca839aa6527695fe98ddc4f32631f2fa20 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 27 Feb 2023 13:14:48 -0500 Subject: [PATCH 09/73] 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. --- tests/test_polyhedron.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 1c9c0469..71540de1 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -236,18 +236,16 @@ def test_edges(): } for name, edgelength in known_shapes.items(): poly = PlatonicFamily.get_shape(name) - if name == "Dodecahedron": - poly.merge_faces(rtol=1) # Test edge_vectors property for edge in poly.edge_vectors: # rtol must be lowered to accomodate small inaccuracies in vertex locations - assert np.isclose(np.linalg.norm(edge), edgelength, rtol=1e-4) + assert np.isclose(np.linalg.norm(edge), edgelength) # Test edges property edges = np.asarray(poly.edges) vertices = poly.vertices veclens = np.linalg.norm(vertices[edges[:, 1]] - vertices[edges[:, 0]], axis=1) - assert np.allclose(veclens, edgelength, rtol=1e-4) + assert np.allclose(veclens, edgelength) def test_curvature(): From 1f4a10673d0928aa1ab8e7126f37855ccfafe3c7 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 27 Feb 2023 22:42:10 -0500 Subject: [PATCH 10/73] Removed unnecessary comment from ``test_edges`` --- tests/test_polyhedron.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 71540de1..375e492b 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -238,7 +238,6 @@ def test_edges(): poly = PlatonicFamily.get_shape(name) # Test edge_vectors property for edge in poly.edge_vectors: - # rtol must be lowered to accomodate small inaccuracies in vertex locations assert np.isclose(np.linalg.norm(edge), edgelength) # Test edges property From d387fca343149839412d6994ff1c832c95b31d72 Mon Sep 17 00:00:00 2001 From: Jen Bradley <55467578+janbridley@users.noreply.github.com> Date: Mon, 27 Feb 2023 22:50:50 -0500 Subject: [PATCH 11/73] Changed ``edge_vectors`` property to return a numpy array Co-authored-by: Bradley Dice --- coxeter/shapes/polyhedron.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 87e7146a..b19e6513 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -380,8 +380,8 @@ def edges(self): @property def edge_vectors(self): """list(:class:`numpy.ndarray`): Get the polyhedron's edges as vectors.""" - edges = np.asarray(self.edges) - return list(self.vertices[edges[:, 1]] - self.vertices[edges[:, 0]]) + edges = self.edges + return self.vertices[edges[:, 1]] - self.vertices[edges[:, 0]] @property def volume(self): From f473e27c380adb96da4e1f6bdb33cd6d7112d32e Mon Sep 17 00:00:00 2001 From: janbridley Date: Wed, 22 Mar 2023 14:39:48 -0400 Subject: [PATCH 12/73] 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 Date: Wed, 22 Mar 2023 15:48:02 -0400 Subject: [PATCH 13/73] Updated ``test_polyhedron`` to work with set edges --- tests/test_polyhedron.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 375e492b..03dc8a06 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -234,6 +234,7 @@ def test_edges(): "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) # Test edge_vectors property @@ -241,7 +242,11 @@ def test_edges(): assert np.isclose(np.linalg.norm(edge), edgelength) # Test edges property - edges = np.asarray(poly.edges) + edges = np.asarray( + [ + *poly.edges, + ] + ) vertices = poly.vertices veclens = np.linalg.norm(vertices[edges[:, 1]] - vertices[edges[:, 0]], axis=1) assert np.allclose(veclens, edgelength) From 01077ccd3d29c08a21ee198cfa0a674beb7ff268 Mon Sep 17 00:00:00 2001 From: janbridley Date: Wed, 22 Mar 2023 17:47:06 -0400 Subject: [PATCH 14/73] Fixed docstring formatting for ``edges`` and ``edge_vectors`` --- coxeter/shapes/polyhedron.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 6e0f34d4..7eefdf3e 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -364,15 +364,14 @@ def faces(self): @property def edges(self): - """set(:class:`tuple`): Get the polyhedron's edges. + """set(tuple(int,int)): 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 Date: Thu, 25 May 2023 10:45:21 -0400 Subject: [PATCH 15/73] 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 Date: Thu, 25 May 2023 11:12:29 -0400 Subject: [PATCH 16/73] test_polyhedron::test_edges now checks edge count --- tests/test_polyhedron.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 03dc8a06..4c5eccfe 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -234,9 +234,20 @@ def test_edges(): "Dodecahedron": np.power(2, 2 / 3) * np.cbrt(1 / (15 + np.sqrt(245))), "Icosahedron": np.cbrt(9 / 5 - 3 / 5 * np.sqrt(5)), } + number_of_edges = { + "Tetrahedron": 6, + "Cube": 12, + "Octahedron": 12, + "Dodecahedron": 30, + "Icosahedron": 30, + } for name, edgelength in known_shapes.items(): poly = PlatonicFamily.get_shape(name) + # Test that the correct number of edges has been found + assert(number_of_edges[name]==len(poly.edges)) + assert(number_of_edges[name]==len(poly.edge_vectors)) + # Test edge_vectors property for edge in poly.edge_vectors: assert np.isclose(np.linalg.norm(edge), edgelength) From 770199b5cabd2355bcede1574158a05800963bd9 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 25 May 2023 15:27:21 +0000 Subject: [PATCH 17/73] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- coxeter/shapes/polyhedron.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 5efd66f5..f93e1b49 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -381,19 +381,21 @@ def edges(self): for i, j in zip(face, np.roll(face, -1)) if i < j ], - key=lambda x: (x[0],x[1]), + key=lambda x: (x[0], x[1]), ) ) @property def edge_vectors(self): """list(tuple(float,float,float)): Get the polyhedron's edges as vectors.""" - return np.array([ - np.subtract( - *self.vertices[[j, i]], - ) - for (i, j) in self.edges - ]) + return np.array( + [ + np.subtract( + *self.vertices[[j, i]], + ) + for (i, j) in self.edges + ] + ) @property def volume(self): From b4f60f9e9c5f2b5c64ccd455916241f70f772a6c Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 4 Aug 2023 12:21:18 -0400 Subject: [PATCH 18/73] Updated edges property to cache result --- coxeter/shapes/polyhedron.py | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index f93e1b49..530e370e 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -368,22 +368,12 @@ def edges(self): 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 Date: Mon, 7 Aug 2023 10:40:53 -0400 Subject: [PATCH 19/73] Updated edges and edge_vectors to make better use of numpy functions and functools caching --- coxeter/shapes/polyhedron.py | 45 +++++++++++++++--------------------- 1 file changed, 18 insertions(+), 27 deletions(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 530e370e..05ab3d34 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -4,6 +4,7 @@ """Defines a polyhedron.""" import warnings +from functools import cached_property import numpy as np import rowan @@ -362,43 +363,33 @@ def faces(self): """ return self._faces - @property + @cached_property def edges(self): """set(tuple(int,int)): 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 Date: Mon, 7 Aug 2023 12:47:22 -0400 Subject: [PATCH 20/73] Updated documentation for ConvexPolyhedron rewrite --- coxeter/shapes/convex_polyhedron.py | 55 ++++++++++++++++------------- 1 file changed, 31 insertions(+), 24 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index bb3472e0..4a7b45d7 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -16,35 +16,42 @@ class ConvexPolyhedron(Polyhedron): """A convex polyhedron. A convex polyhedron is defined as the convex hull of its vertices. The - class is a simple extension of :class:`~.Polyhedron` that builds the - faces from the simplices of the convex hull. This class also includes - various additional properties that can be used to characterize the - geometric features of the polyhedron. + class is an extension of :class:`~.Polyhedron` that builds the + faces from the simplices of the convex hull. Simplices are stored and class methods + are optimized to make use of the triangulation, as well as special properties of + convex solids in three dimensions. This class also includes various additional + properties that can be used to characterize geometric features of the polyhedron. Args: vertices (:math:`(N, 3)` :class:`numpy.ndarray`): The vertices of the polyhedron. + fast (bool, optional): + Creation mode for the polyhedron. fast=False (default) will perform all + checks and precalculate most properties. fast=True will precalculate some + properties, but will not find face neighbors or sort face indices. These + calculations will instead be performed when required by another method. + (Default value: False). Example: >>> cube = coxeter.shapes.ConvexPolyhedron( - ... [[1, 1, 1], [1, -1, 1], [1, 1, -1], [1, -1, -1], - ... [-1, 1, 1], [-1, -1, 1], [-1, 1, -1], [-1, -1, -1]]) + ... [[-1, -1, -1], [-1, -1, 1], [-1, 1, -1], [-1, 1, 1], + ... [1, -1, -1], [1, -1, 1], [1, 1, -1], [1, 1, 1]]) >>> import numpy as np >>> assert np.isclose(cube.asphericity, 1.5) >>> bounding_sphere = cube.minimal_bounding_sphere >>> assert np.isclose(bounding_sphere.radius, np.sqrt(3)) - >>> cube.center + >>> cube.centroid array([0., 0., 0.]) >>> circumsphere = cube.circumsphere >>> assert np.isclose(circumsphere.radius, np.sqrt(3)) >>> cube.faces - [array([4, 5, 1, 0], dtype=int32), array([0, 2, 6, 4], dtype=int32), - array([6, 7, 5, 4], dtype=int32), array([0, 1, 3, 2], dtype=int32), - array([5, 7, 3, 1], dtype=int32), array([2, 3, 7, 6], dtype=int32)] + [array([0, 2, 6, 4], dtype=int32), array([0, 4, 5, 1], dtype=int32), + array([4, 6, 7, 5], dtype=int32), array([0, 1, 3, 2], dtype=int32), + array([2, 3, 7, 6], dtype=int32), array([1, 5, 7, 3], dtype=int32)] >>> cube.gsd_shape_spec - {'type': 'ConvexPolyhedron', 'vertices': [[1.0, 1.0, 1.0], [1.0, -1.0, 1.0], - [1.0, 1.0, -1.0], [1.0, -1.0, -1.0], [-1.0, 1.0, 1.0], [-1.0, -1.0, 1.0], - [-1.0, 1.0, -1.0], [-1.0, -1.0, -1.0]]} + {'type': 'ConvexPolyhedron', 'vertices': [[-1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], + [-1.0, 1.0, -1.0], [-1.0, 1.0, 1.0], [1.0, -1.0, -1.0], [1.0, -1.0, 1.0], + [1.0, 1.0, -1.0], [1.0, 1.0, 1.0]]} >>> assert np.allclose( ... cube.inertia_tensor, ... np.diag([16. / 3., 16. / 3., 16. / 3.])) @@ -57,12 +64,12 @@ class is a simple extension of :class:`~.Polyhedron` that builds the [array([1, 2, 3, 4]), array([0, 2, 3, 5]), array([0, 1, 4, 5]), array([0, 1, 4, 5]), array([0, 2, 3, 5]), array([1, 2, 3, 4])] >>> cube.normals - array([[-0., -0., 1.], - [-0., 1., -0.], - [-1., 0., -0.], + array([[-0., -0., -1.], + [ 0., -1., 0.], [ 1., -0., -0.], - [-0., -1., 0.], - [-0., -0., -1.]]) + [-1., -0., -0.], + [ 0., 1., -0.], + [-0., -0., 1.]]) >>> cube.num_faces 6 >>> cube.num_vertices @@ -71,14 +78,14 @@ class is a simple extension of :class:`~.Polyhedron` that builds the 24.0 >>> assert np.isclose(cube.tau, 3. / 8. * np.pi) >>> cube.vertices - array([[ 1., 1., 1.], - [ 1., -1., 1.], - [ 1., 1., -1.], - [ 1., -1., -1.], - [-1., 1., 1.], + array([[-1., -1., -1.], [-1., -1., 1.], [-1., 1., -1.], - [-1., -1., -1.]]) + [-1., 1., 1.], + [ 1., -1., -1.], + [ 1., -1., 1.], + [ 1., 1., -1.], + [ 1., 1., 1.]]) >>> assert np.isclose(cube.volume, 8.) """ From 2f22ab2159bdc90b18102c03187be928a56071d3 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:22:46 -0400 Subject: [PATCH 21/73] Updated __init__ and shape creation helpers --- coxeter/shapes/convex_polyhedron.py | 219 +++++++++++++++++++++++++++- 1 file changed, 215 insertions(+), 4 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 4a7b45d7..f957a7fe 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -4,6 +4,7 @@ """Defines a convex polyhedron.""" import warnings +from collections import defaultdict import numpy as np from scipy.spatial import ConvexHull @@ -90,10 +91,220 @@ class is an extension of :class:`~.Polyhedron` that builds the """ - def __init__(self, vertices): - hull = ConvexHull(vertices) - super().__init__(vertices, hull.simplices, True) - self.merge_faces() + def __init__(self, vertices, fast=True): + self._vertices = np.array(vertices, dtype=np.float64) + self._ndim = self._vertices.shape[1] + self._convex_hull = ConvexHull(self._vertices) + self._faces_are_convex = True + + if not len(self._convex_hull.vertices) == len(self._vertices): + warnings.warn("Not all vertices used for hull.", RuntimeWarning) + + # Remove internal points + self._vertices = self._vertices[self._convex_hull.vertices] + + # Transfer data in from convex hull, then clean up the results. + self._consume_hull() + self._combine_simplices() + + # Sort simplices. This method also calculates simplex equations and the centroid + self._sort_simplices() + + # If mode is NOT fast, perform additional initializiation steps + if fast is False: + self._find_coplanar_simplices() + self.sort_faces() + else: + self._faces_are_sorted = False + + def _consume_hull(self): + """Extract data from ConvexHull. + + Data is moved from convex hull into private variables. This method deletes the + original hull in order to avoid double storage, and chec. + """ + assert ( + self._ndim == self._convex_hull.ndim + ), "Input points are coplanar or close to coplanar." + + self._simplices = self._convex_hull.simplices[:] + self._simplex_equations = self._convex_hull.equations[:] + self._simplex_neighbors = self._convex_hull.neighbors[:] + self._volume = self._convex_hull.volume + self._area = self._convex_hull.area + self._maximal_extents = np.array( + [self._convex_hull.min_bound, self._convex_hull.max_bound] + ) + + # Clean up the result. + del self._convex_hull + + def _combine_simplices(self, rounding: int = 15): + """Combine simplices into faces, merging based on simplex equations. + + Coplanar faces will have identical equations (within rounding tolerance). Values + should be slightly larger than machine epsilon (e.g. rounding=15 for ~1e-15) + + Args: + rounding (int, optional): + Integer number of decimal places to round equations to. + (Default value: 15). + + """ + equation_groups = defaultdict(list) + + # Iterate over all simplex equations + for i, equation in enumerate(self._simplex_equations): + # Convert to hashable key + equation_key = tuple(equation.round(rounding)) + + # Store vertex indices from the new simplex under the correct equation key + equation_groups[equation_key].extend(self._simplices[i]) + + # Combine elements with the same plan equation and remove duplicate indices + ragged_faces = [ + np.fromiter(set(group), np.int32) for group in equation_groups.values() + ] + self._faces = ragged_faces + self._equations = np.array(list(equation_groups.keys())) + + def _find_coplanar_simplices(self, rounding: int = 15): + """ + Get lists of simplex indices for coplanar simplices. + + Args: + rounding (int, optional): + Integer number of decimal places to round equations to. + (Default value: 15). + + + """ + # Combine simplex indices + equation_groups = defaultdict(list) + + # Iterate over all simplex equations + for i, equation in enumerate(self._simplex_equations): + # Convert equation into hashable tuple + equation_key = tuple(equation.round(rounding)) + equation_groups[equation_key].append(i) + ragged_coplanar_indices = [ + np.fromiter(set(group), np.int32) for group in equation_groups.values() + ] + + self._coplanar_simplices = ragged_coplanar_indices + + def _sort_simplices(self): + """Reorder simplices counterclockwise relatative to the plane they lie on. + + This does NOT change the *order* of simplices in the list. + """ + # Get correct-quadrant angles about the simplex normal + vertices = self._vertices[self._simplices] + + # Get the absolute angles of each vertex and fit to unit circle + angles = np.arctan2(vertices[:, 1], vertices[:, 0]) + angles = np.mod(angles - angles[0], 2 * np.pi) + + # Calculate distances + distances = np.linalg.norm(vertices, axis=1) + + # Create a tuple of distances and angles to use for lexicographical sorting + vert_order = np.lexsort((distances, angles)) + + # Apply orientation reordering to every simplex + self._simplices = np.take_along_axis( + self._simplices, np.argsort(vert_order), axis=1 + ) + + # Compute N,3,2 array of simplex edges from an N,3 array of simplices + def _tri_to_edge_tuples(simplices, shift=-3): + # edges = np.column_stack((simplices, np.roll(simplices, shift, axis=1))) + edges = np.roll( + np.repeat(simplices, repeats=2, axis=1), shift=shift, axis=1 + ) + edges = edges.reshape(-1, 3, 2) + # Convert to tuples for fast comparison + return [[tuple(edge) for edge in triedge] for triedge in edges] + + simplex_edges = _tri_to_edge_tuples(self._simplices, shift=-1) + + # Now, reorient simplices to match the orientation of simplex 0 + visited_simplices = [] + remaining_simplices = [0] + while len(remaining_simplices): + current_simplex = remaining_simplices[-1] + visited_simplices.append(current_simplex) + remaining_simplices.pop() + + # Search for common edges between pairs of simplices, then check the + # ordering of the edge to determine relative face orientation. + current_edges = simplex_edges[current_simplex] + for neighbor in self._simplex_neighbors[current_simplex]: + if neighbor in visited_simplices: + continue + remaining_simplices.append(neighbor) + + # Two faces can only share a single edge (otherwise they would + # be coplanar), so we can break as soon as we find the + # neighbor. Flip the neighbor if the edges are identical. + for edge in simplex_edges[neighbor]: + if edge in current_edges: + self._simplices[neighbor] = self._simplices[neighbor][::-1] + simplex_edges[neighbor] = [ + (j, i) for i, j in simplex_edges[neighbor] + ][::-1] + break + elif edge[::-1] in current_edges: + break + visited_simplices.append(neighbor) + + # Flip if calcualted volume is negative + if self._calculate_signed_volume() < 0: + self._simplices = self._simplices[:, ::-1, ...] + + # Recompute simplex equations and centroid from the new triangulation. + self._find_simplex_equations() + self._centroid_from_triangulated_surface() + + def _calculate_signed_volume(self): + """Internal method to calculate the signed volume of the polyhedron. + + This class splits the shape into tetrahedra, then sums their contributing + volumes. The external volume property will always be a positive value, but + accessing the signed volume can be useful for some mathematical operations. + + Returns: + float: Signed volume of the polyhedron. + """ + signed_volume = np.sum(np.linalg.det(self._vertices[self._simplices]) / 6) + self._volume = abs(signed_volume) + return signed_volume + + def _find_simplex_equations(self): + """Find the plane equations of the polyhedron simplices.""" + abc = self._vertices[self._simplices] + a = abc[:, 0] + b = abc[:, 1] + c = abc[:, 2] + n = np.cross((b - a), (c - a)) + n /= np.linalg.norm(n, axis=1)[:, None] + + _equations = np.empty((len(self._simplices), 4)) + _equations[:, :3] = n + _equations[:, 3] = -np.einsum("ij,ij->i", n, a) + self._simplex_equations = _equations + + def _centroid_from_triangulated_surface(self): + abc = self._vertices[self._simplices] + a = abc[:, 0] + b = abc[:, 1] + c = abc[:, 2] + n = np.cross((b - a), (c - a)) + self._centroid = ( + 1 + / (48 * self._volume) + * np.sum(n * ((a + b) ** 2 + (b + c) ** 2 + (a + c) ** 2), axis=0) + ) @property def mean_curvature(self): From f598272ee63839b3014fdaed57953af4b811c6eb Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:26:26 -0400 Subject: [PATCH 22/73] Updated centroid getter and setter --- coxeter/shapes/convex_polyhedron.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index f957a7fe..5b5c4ac3 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -294,6 +294,23 @@ def _find_simplex_equations(self): _equations[:, 3] = -np.einsum("ij,ij->i", n, a) self._simplex_equations = _equations + @property + def centroid(self): + """:math:`(3, )` :class:`numpy.ndarray` of float: Get or set the center of mass. + + The centroid is calculated using the curl theorem over the surface simplices. + """ + return self._centroid + + @centroid.setter + def centroid(self, value): + assert len(value) == 3, "Centroid must be a point in 3-space." + self._vertices += np.asarray(value) - self.centroid + self._find_equations() + self._find_simplex_equations() + self._centroid_from_triangulated_surface() + self._calculate_signed_volume() + def _centroid_from_triangulated_surface(self): abc = self._vertices[self._simplices] a = abc[:, 0] From 5dfef18f43020afb503801c99465ffef9a3a9ad8 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:28:58 -0400 Subject: [PATCH 23/73] Updated volume getter, setter, and _rescale function --- coxeter/shapes/convex_polyhedron.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 5b5c4ac3..e1ceab83 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -5,6 +5,7 @@ import warnings from collections import defaultdict +from numbers import Number import numpy as np from scipy.spatial import ConvexHull @@ -266,6 +267,33 @@ def _tri_to_edge_tuples(simplices, shift=-3): self._find_simplex_equations() self._centroid_from_triangulated_surface() + @property + def volume(self): + """float: Get or set the polyhedron's volume.""" + return self._volume + + @volume.setter + def volume(self, value: Number): + scale_factor = np.cbrt(value / self._volume) + self._rescale(scale_factor) + + def _rescale(self, scale_factor: Number): + """Scale polytope by changing the length of the edges. + + Args: + scale_factor (int or float): + Multiplier to scale edges by. Volume and surface area setters preconvert + the scale_factor to the correct value for the desired property. + """ + self._vertices *= scale_factor + self._equations[:, 3] *= scale_factor + self._simplex_equations[:, 3] *= scale_factor + self._volume = self._volume * scale_factor**3 + self._area = self._area * scale_factor**2 + + # Recalculate centroid of shape + self._centroid_from_triangulated_surface() + def _calculate_signed_volume(self): """Internal method to calculate the signed volume of the polyhedron. From 0d21d9825cb082bc6f82de01731d8da8792843e1 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:34:12 -0400 Subject: [PATCH 24/73] Improved _find_equations method --- coxeter/shapes/convex_polyhedron.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index e1ceab83..d14b1cbe 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -308,6 +308,31 @@ def _calculate_signed_volume(self): self._volume = abs(signed_volume) return signed_volume + def _find_equations(self): + """Find the plane equations of the polyhedron faces.""" + # This method only takes the first three items from each face, so it is + # unaffected by face structure or ordering. + point_on_face_indices = [] + for face in self.faces: + point_on_face_indices.append(face[0:3]) + vertices = self._vertices[point_on_face_indices] + + # Calculate the directions of the normals + v1 = vertices[:, 2] - vertices[:, 1] + v2 = vertices[:, 0] - vertices[:, 1] + normals = np.cross(v1, v2) + + # Normalize the normals + norms = np.linalg.norm(normals, axis=1) + normals /= norms[:, None] + + _equations = np.empty((len(self._faces), 4)) + _equations[:, :3] = normals + _equations[:, 3] = -np.einsum( + "ij,ij->i", normals, vertices[:, 0] + ) # dot product + self._equations = _equations + def _find_simplex_equations(self): """Find the plane equations of the polyhedron simplices.""" abc = self._vertices[self._simplices] From d5b202cc2c5da32548d94a8cd5efc149ef934253 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:36:15 -0400 Subject: [PATCH 25/73] Added equations property --- coxeter/shapes/convex_polyhedron.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index d14b1cbe..e2bf62b4 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -333,6 +333,14 @@ def _find_equations(self): ) # dot product self._equations = _equations + @property + def equations(self): + """:math:`(N, 4)` :class:`numpy.ndarray`: Get plane equations for each face. + + Sign convention matches Scipy Convex Hull (ax + by + cz + d = 0). + """ + return self._equations + def _find_simplex_equations(self): """Find the plane equations of the polyhedron simplices.""" abc = self._vertices[self._simplices] From 2e3914c6ee763e102261e7750d60ae0784a37a9e Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:37:21 -0400 Subject: [PATCH 26/73] Updated faces property to work with optional sorting --- coxeter/shapes/convex_polyhedron.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index e2bf62b4..297bde8f 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -384,6 +384,13 @@ def _centroid_from_triangulated_surface(self): * np.sum(n * ((a + b) ** 2 + (b + c) ** 2 + (a + c) ** 2), axis=0) ) + @property + def faces(self): + """list(:class:`numpy.ndarray`): Get the polyhedron's faces.""" + if self._faces_are_sorted is False: + self.sort_faces() + return self._faces + @property def mean_curvature(self): r"""float: The integrated, normalized mean curvature. From ed33ca62a9798b7c94676d20369534cce8a306a4 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:38:08 -0400 Subject: [PATCH 27/73] Added normals property --- coxeter/shapes/convex_polyhedron.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 297bde8f..32c531be 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -341,6 +341,11 @@ def equations(self): """ return self._equations + @property + def normals(self): + """:math:`(N, 3)` :class:`numpy.ndarray`: Get normal vectors for each face.""" + return self._equations[:, :3] + def _find_simplex_equations(self): """Find the plane equations of the polyhedron simplices.""" abc = self._vertices[self._simplices] From 0dbc2681cf07b47d49ac685cfeb705c143a21783 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:40:19 -0400 Subject: [PATCH 28/73] Updated sort_faces --- coxeter/shapes/convex_polyhedron.py | 37 +++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 32c531be..07dc73a3 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -8,6 +8,7 @@ from numbers import Number import numpy as np +import rowan from scipy.spatial import ConvexHull from .polyhedron import Polyhedron @@ -396,6 +397,42 @@ def faces(self): self.sort_faces() return self._faces + def sort_faces(self): + """Reorder faces counterclockwise relatative to the plane they lie on. + + This does NOT change the *order* of faces in the list. + """ + self._faces_are_sorted = True + + # Get correct-quadrant angles about the face normal + sorted_faces = [] + + for i, face in enumerate(self._faces): + vertices = self._vertices[face] + + # Rotate the face's points into the XY plane + normal = self._equations[i][:3] + rotation, _ = rowan.mapping.kabsch( + [normal, -normal], [[0, 0, 1], [0, 0, -1]] + ) + vertices = np.dot(vertices - np.mean(vertices, axis=0), rotation.T) + + # Get the absolute angles of each vertex and fit to unit circle + angles = np.arctan2(vertices[:, 1], vertices[:, 0]) + angles = np.mod(angles - angles[0], 2 * np.pi) + + # Calculate distances + distances = np.linalg.norm(vertices, axis=1) + + # Create a tuple of distances and angles to use for lexicographical sorting + vert_order = np.lexsort((distances, angles)) + + # Apply reordering to every simplex + sorted_faces.append(face[vert_order]) + + self._faces = sorted_faces + self._find_neighbors() + @property def mean_curvature(self): r"""float: The integrated, normalized mean curvature. From e1a700d13c99c79d8e4c32e19c3321a4b5135e03 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:46:27 -0400 Subject: [PATCH 29/73] Updated neighbors property --- coxeter/shapes/convex_polyhedron.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 07dc73a3..c059f8b7 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -397,6 +397,18 @@ def faces(self): self.sort_faces() return self._faces + @property + def neighbors(self): + r"""list(:class:`numpy.ndarray`): Get neighboring pairs of faces. + + The neighbors are provided as a list where the :math:`i^{\text{th}}` + element is an array of indices of faces that are neighbors of face + :math:`i`. + """ + if self._faces_are_sorted is False: + self.sort_faces() + return self._neighbors + def sort_faces(self): """Reorder faces counterclockwise relatative to the plane they lie on. From d5534bec8b423f3f7a7943cd51cf6c537e6e15aa Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:47:33 -0400 Subject: [PATCH 30/73] Added simplices property and updated _surface_triangulation --- coxeter/shapes/convex_polyhedron.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index c059f8b7..9f94b472 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -445,6 +445,25 @@ def sort_faces(self): self._faces = sorted_faces self._find_neighbors() + def _surface_triangulation(self): + """Output the vertices of simplices composing the polyhedron's surface. + + Returns: + :math:`(N,3,3)` :class:`numpy.ndarray`: + Array of vertices for simplices composing the polyhedron's surface. + """ + return self.vertices[self._simplices] + + @property + def simplices(self): + """Output the vertex indices of simplices composing the polyhedron's surface. + + Returns: + :math:`(N,3)` :class:`numpy.ndarray`: + Array of vertex indices of simplices making up the polyhedron's surface. + """ + return self._simplices + @property def mean_curvature(self): r"""float: The integrated, normalized mean curvature. From a728486b3f74969df69d039f3ee3f9cb5ab35636 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:51:48 -0400 Subject: [PATCH 31/73] Added _find_triangle_array_area helper --- coxeter/shapes/convex_polyhedron.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 9f94b472..1210f365 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -464,6 +464,32 @@ def simplices(self): """ return self._simplices + def _find_triangle_array_area(self, triangle_vertices, sum_result=True): + """ + Get the areas of each triangle in an input array. + + Args: + angles (:math:`(N, 3, 3)` :class:`numpy.ndarray`): + Array of vertices for the triangles that will have their area computed. + sum_result (bool): + Whether the output should be summed. + (Default value: True). + + Returns: + :math:`(N, )` :class:`numpy.ndarray` or float: + Boolean array indicating which points are contained in thepolyhedron. + If sum_result is True, a single value is returned. + + """ + v1 = triangle_vertices[:, 2] - triangle_vertices[:, 1] + v2 = triangle_vertices[:, 0] - triangle_vertices[:, 1] + unnormalized_normals = np.cross(v1,v2,axis=1) + + if sum_result: + return np.sum(np.linalg.norm(unnormalized_normals, axis=1)) / 2 + else: + return np.linalg.norm(unnormalized_normals, axis=1) / 2 + @property def mean_curvature(self): r"""float: The integrated, normalized mean curvature. From 135e0a12c3663da7295bc768a599e63c6cb77312 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:55:43 -0400 Subject: [PATCH 32/73] Add faster surface area getter, setter, and calculation --- coxeter/shapes/convex_polyhedron.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 1210f365..d4684daf 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -278,6 +278,21 @@ def volume(self, value: Number): scale_factor = np.cbrt(value / self._volume) self._rescale(scale_factor) + @property + def surface_area(self): + """float: Get or set the surface area.""" + return self._area + + @surface_area.setter + def surface_area(self, value: Number): + scale_factor = np.sqrt(value / self._area) + self._rescale(scale_factor) + + def _calculate_surface_area(self): + new_area = self._find_triangle_array_area(self._vertices[self._simplices]) + self._area = new_area + return self._area + def _rescale(self, scale_factor: Number): """Scale polytope by changing the length of the edges. @@ -483,8 +498,8 @@ def _find_triangle_array_area(self, triangle_vertices, sum_result=True): """ v1 = triangle_vertices[:, 2] - triangle_vertices[:, 1] v2 = triangle_vertices[:, 0] - triangle_vertices[:, 1] - unnormalized_normals = np.cross(v1,v2,axis=1) - + unnormalized_normals = np.cross(v1, v2, axis=1) + if sum_result: return np.sum(np.linalg.norm(unnormalized_normals, axis=1)) / 2 else: From f2c55e1ddc1442d3b5ae54b2fce11c25646abfcc Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 13:57:02 -0400 Subject: [PATCH 33/73] Added get_face_area method --- coxeter/shapes/convex_polyhedron.py | 45 +++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index d4684daf..7baf907d 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -324,6 +324,51 @@ def _calculate_signed_volume(self): self._volume = abs(signed_volume) return signed_volume + def get_face_area(self, face=None): + """Get the total surface area of a set of faces. + + Args: + faces (int, sequence, or None): + The index of a face or a set of face indices for which to + find the area. If None, finds the area of all faces. + (Default value: None). + + Returns: + :class:`numpy.ndarray`: The area of each face. + + Example: + >>> cube = coxeter.shapes.ConvexPolyhedron( + ... [[1, 1, 1], [1, -1, 1], [1, 1, -1], [1, -1, -1], + ... [-1, 1, 1], [-1, -1, 1], [-1, 1, -1], [-1, -1, -1]]) + >>> cube = coxeter.shapes.Polyhedron( + ... vertices=cube.vertices,faces=cube.faces) + >>> import numpy as np + >>> assert np.allclose( + ... cube.get_face_area([1, 2, 3]), + ... [4., 4., 4.]) + + """ + self._simplex_areas = self._find_triangle_array_area( + self._vertices[self._simplices], sum_result=False + ) + if face is None: + return [ + np.sum(self._simplex_areas[self._coplanar_simplices[fac]]) + for fac in range(self.num_faces) + ] + elif face == "total": # return total surface area + return np.sum(self._simplex_areas) + elif hasattr(face, "__len__"): + # For list of input vertices + # Combine coplanar simplices + return [ + np.sum(self._simplex_areas[self._coplanar_simplices[fac]]) + for fac in face + ] + else: + # For integer input (single face) + return np.sum(self._simplex_areas[self._coplanar_simplices[face]]) + def _find_equations(self): """Find the plane equations of the polyhedron faces.""" # This method only takes the first three items from each face, so it is From 1721a9e44f6b8b9dff06e20bede720f37ee0889d Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 14:00:45 -0400 Subject: [PATCH 34/73] Added face_centroids property and _get_face_centroids --- coxeter/shapes/convex_polyhedron.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 7baf907d..79343558 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -457,6 +457,33 @@ def faces(self): self.sort_faces() return self._faces + def _find_face_centroids(self): + simplex_centroids = np.mean(self._vertices[self._simplices], axis=1) # (N,3) + self._simplex_areas = self._find_triangle_array_area( + self._vertices[self._simplices], sum_result=False + ) # (N,) + self._face_centroids = [] + for face in self._coplanar_simplices: + self._face_centroids.append( + np.sum( + simplex_centroids[face] * self._simplex_areas[face][:, None], + axis=0, + ) + / np.sum(self._simplex_areas[face]), # Rescale by area of face + ) + self._face_centroids = np.array(self._face_centroids) + + @property + def face_centroids(self): + """Calculate the centroid (center of mass) of each polygonal face. + + Returns: + :math:`(N,3)` :class:`numpy.ndarray`: + Array of centroids for each face. + """ + self._find_face_centroids() + return self._face_centroids + @property def neighbors(self): r"""list(:class:`numpy.ndarray`): Get neighboring pairs of faces. From b62a3aa3eb2f0669bd032aa99e1cd9a312568499 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 14:37:15 -0400 Subject: [PATCH 35/73] Updated inertia_tensor property + _compute and diagonalize methods Compute_inertia_tensor needed to be updated because the previous algorithm does not work if the tetrahedral decomposition of the shape includes tetrahedra with negative volumes. This new method, implemented based on an algorithm published by A. M. Messner in 1980 is significantly faster and provides slightly more accurate results. Because of this, the relevant pytest has been updated to better check the method. --- coxeter/shapes/convex_polyhedron.py | 118 ++++++++++++++++++++++++++++ doc/source/coxeter.bib | 12 ++- tests/test_polyhedron.py | 26 +++--- 3 files changed, 145 insertions(+), 11 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 79343558..6de19712 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -13,6 +13,7 @@ from .polyhedron import Polyhedron from .sphere import Sphere +from .utils import translate_inertia_tensor class ConvexPolyhedron(Polyhedron): @@ -577,6 +578,123 @@ def _find_triangle_array_area(self, triangle_vertices, sum_result=True): else: return np.linalg.norm(unnormalized_normals, axis=1) / 2 + @property + def inertia_tensor(self): + """:math:`(3, 3)` :class:`numpy.ndarray`: Get the inertia tensor. + + The inertia tensor for convex shapes is computed using the algorithm described + in :cite:`Messner1980`. + + Note: + For improved stability, the inertia tensor is computed about the + center of mass and then shifted rather than directly computed in + the global frame. + """ + it = self._compute_inertia_tensor() + return translate_inertia_tensor(self.center, it, self.volume) + + def _compute_inertia_tensor(self, centered=True): + """Compute the inertia tensor. + + Internal function for computing the inertia tensor that supports both + centered and uncentered calculations. + """ + + def _quadrature_points(abc): + """Compute quadrature points at which to evaluate the integral.""" + scalars = [ + [5 / 3, 5 / 3, 5 / 3], + [[1], [1], [3]], + [[3], [1], [1]], + [[1], [3], [1]], + ] + + q = np.zeros((abc.shape[0], 3, 4)) + + for i in range(4): + q[:, :, i] = np.sum(abc * scalars[i], axis=1) + q /= 5 + return q + + # Define triangle array + abc = self.vertices[self.simplices] + if centered: + abc -= self.centroid + n = self._simplex_equations[:, :3] + nt = n.T + q = _quadrature_points(abc) + q2 = q**2 + q3 = q**3 + + # Define quadrature weights (from integrating a third order polynomial over a + # triangular domain). The weights have been normalized such that their sum is 1. + w = np.array([[-9 / 16, 25 / 48, 25 / 48, 25 / 48]]).T + + # Face simplex areas + at = self._find_triangle_array_area(abc, sum_result=False) * 2 + + # Calculate and sum on-diagonal contributions to the moment of inertia. + def i_nn(nt, q3, w, at, sub): + return np.einsum("ij, j, jik, kl ->", nt[sub, :], at, q3[:, sub, :], w) / 6 + + # Calculate and sum off-diagonal contributions to the moment of inertia. + def i_nm(n, q, q2, w, at, sub): + return ( + -( + np.einsum( + "ij,ki,k,k->", + w, + q2[:, sub[0], :] * q[:, sub[1], :], + n[:, sub[0]], + at, + ) + + np.einsum( + "ij,ki,k,k->", + w, + q[:, sub[0], :] * q2[:, sub[1], :], + n[:, sub[1]], + at, + ) + ) + / 8 + ) + + i_xx = i_nn(nt, q3, w, at, sub=[1, 2]) + i_xy = i_nm(n, q, q2, w, at, sub=[0, 1]) + i_xz = i_nm(n, q, q2, w, at, sub=[0, 2]) + i_yy = i_nn(nt, q3, w, at, sub=[0, 2]) + i_yz = i_nm(n, q, q2, w, at, sub=[1, 2]) + i_zz = i_nn(nt, q3, w, at, sub=[0, 1]) + + return np.array([[i_xx, i_xy, i_xz], [i_xy, i_yy, i_yz], [i_xz, i_yz, i_zz]]) + + def diagonalize_inertia(self): + """Orient the shape along its principal axes. + + The principal axes of a shape are defined by the eigenvectors of the inertia + tensor. This method computes the inertia tensor of the shape, diagonalizes it, + and then rotates the shape by the corresponding orthogonal transformation. + + Example: + >>> cube = coxeter.shapes.ConvexPolyhedron( + ... [[1, 1, 1], [1, -1, 1], [1, 1, -1], [1, -1, -1], + ... [-1, 1, 1], [-1, -1, 1], [-1, 1, -1], [-1, -1, -1]]) + >>> cube.diagonalize_inertia() + >>> cube.vertices + array([[ 1., 1., 1.], + [ 1., -1., 1.], + [ 1., 1., -1.], + [ 1., -1., -1.], + [-1., 1., 1.], + [-1., -1., 1.], + [-1., 1., -1.], + [-1., -1., -1.]]) + + """ + principal_moments, principal_axes = np.linalg.eigh(self.inertia_tensor) + self._vertices = np.dot(self._vertices, principal_axes) + self._sort_simplices() + @property def mean_curvature(self): r"""float: The integrated, normalized mean curvature. diff --git a/doc/source/coxeter.bib b/doc/source/coxeter.bib index 041d6fcf..6a69af07 100644 --- a/doc/source/coxeter.bib +++ b/doc/source/coxeter.bib @@ -1,3 +1,14 @@ +@article{Messner1980, +author={Messner, A. M. and Taylor, G. Q.}, +title={Algorithm 550: Solid polyhedron measures [z]}, +journal={ACM Transactions on Mathematical Software}, +volume={6}, +number={1}, +pages={121–130}, +year={1980}, +doi={10.1145/355873.355885}, +} + @article{Irrgang2017, author = {Irrgang, M. Eric and Engel, Michael and Schultz, Andrew J. and Kofke, David A. and Glotzer, Sharon C.}, title = {Virial Coefficients and Equations of State for Hard Polyhedron Fluids}, @@ -35,7 +46,6 @@ @article{Kallay2006 doi = {10.1080/2151237X.2006.10129220}, } - @article{Chen2014, title = {Complexity in Surfaces of Densest Packings for Families of Polyhedra}, author = {Chen, Elizabeth R. and Klotsa, Daphne and Engel, Michael and Damasceno, Pablo F. and Glotzer, Sharon C.}, diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 34aab06f..16397221 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -193,27 +193,33 @@ def test_surface_area_shapes(poly): assert np.isclose(poly.surface_area, hull.area) -# This test is a bit slow (a couple of minutes), so skip running it locally. -@pytest.mark.skipif( - os.getenv("CI", "false") != "true" and os.getenv("CIRCLECI", "false") != "true", - reason="Test is too slow to run during rapid development", +# This test is slow at high precisions. Run fast locally, and test in detail on CircleCI +@pytest.mark.parametrize( + "atol", + [ + 5e-2 + if os.getenv("CI", "false") == "true" + or os.getenv("CIRCLECI", "false") == "true" + else 1e-1 + ], ) @named_damasceno_shapes_mark -def test_moment_inertia_damasceno_shapes(shape): +def test_moment_inertia_damasceno_shapes(shape, atol): # These shapes pass the test for a sufficiently high number of samples, but # the number is too high to be worth running them regularly. bad_shapes = [ "Augmented Truncated Dodecahedron", "Deltoidal Hexecontahedron", "Disdyakis Triacontahedron", - "Truncated Dodecahedron", - "Truncated Icosidodecahedron", "Metabiaugmented Truncated Dodecahedron", - "Pentagonal Hexecontahedron", + "Parabiaugmented Truncated Dodecahedron", "Paragyrate Diminished Rhombicosidodecahedron", + "Pentagonal Hexecontahedron", + "Rhombic Enneacontahedron", "Square Cupola", "Triaugmented Truncated Dodecahedron", - "Parabiaugmented Truncated Dodecahedron", + "Truncated Dodecahedron", + "Truncated Icosidodecahedron", ] if shape["name"] in ["RESERVED", "Sphere"] + bad_shapes: return @@ -228,7 +234,7 @@ def test_moment_inertia_damasceno_shapes(shape): while num_samples < 1e8: try: mc_result = compute_inertia_mc(poly.vertices, volume, num_samples) - assert np.allclose(coxeter_result, mc_result, atol=1e-1) + assert np.allclose(coxeter_result, mc_result, atol=atol) accept = True break except AssertionError: From 2fe4d05e19b197ca7d9a1bd42da916dce8f4405d Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 14:40:04 -0400 Subject: [PATCH 36/73] Added get_dihedral --- coxeter/shapes/convex_polyhedron.py | 31 +++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 6de19712..389ac93d 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -798,6 +798,37 @@ def maximal_centered_bounded_sphere(self): min_distance = -np.max(distances) return Sphere(min_distance, center) + def get_dihedral(self, a, b): + """Get the dihedral angle between a pair of faces. + + The dihedral is computed from the dot product of the face normals. + + Args: + a (int): + The index of the first face. + b (int): + The index of the second face. + + Returns: + float: The dihedral angle in radians. + + Example: + >>> cube = coxeter.shapes.ConvexPolyhedron( + ... [[1, 1, 1], [1, -1, 1], [1, 1, -1], [1, -1, -1], + ... [-1, 1, 1], [-1, -1, 1], [-1, 1, -1], [-1, -1, -1]]) + >>> cube = coxeter.shapes.Polyhedron( + ... vertices=cube.vertices, faces=cube.faces) + >>> import numpy as np + >>> assert np.isclose(cube.get_dihedral(1, 2), np.pi / 2.) + + """ + if self._faces_are_sorted is False: + self.sort_faces() + if b not in self.neighbors[a]: + raise ValueError("The two faces are not neighbors.") + n1, n2 = self._equations[[a, b], :3] + return np.arccos(np.dot(-n1, n2)) + def _plato_primitive(self, backend): return backend.ConvexPolyhedra( positions=np.array([self.center]), From 158ac57e7bf848abea05f8dfbae43b4f42d86728 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 15:12:04 -0400 Subject: [PATCH 37/73] Added more stringent pytests for centroid and centroid setter The centroid methods are core functionality to this class, and prior to this commit was not being tested. This fixes that, with a two-part test that checks for initial and repeated setting behavior. The monte carlo centroid compute has also been updated to test more samples on CircleCI. Test sampling has been tuned to account for the tighter tolerances --- tests/test_polyhedron.py | 34 ++++++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 16397221..e0787dcd 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -799,16 +799,25 @@ def test_repr_convex(convex_cube): assert str(convex_cube), str(eval(repr(convex_cube))) +@pytest.mark.parametrize( + "atol", + [ + 5e-3 + if os.getenv("CI", "false") == "true" + or os.getenv("CIRCLECI", "false") == "true" + else 1e-2 + ], +) @named_damasceno_shapes_mark -def test_center(shape): +def test_center(shape, atol): poly = ConvexPolyhedron(shape["vertices"]) coxeter_result = poly.center - num_samples = 1000 + num_samples = 5000 accept = False - while num_samples < 1e8: + while num_samples < 5e7: try: mc_result = compute_centroid_mc(shape["vertices"], num_samples) - assert np.allclose(coxeter_result, mc_result, atol=1e-1) + assert np.allclose(coxeter_result, mc_result, atol=atol) accept = True break except AssertionError: @@ -821,3 +830,20 @@ def test_center(shape): shape["name"], mc_result, coxeter_result ) ) + + +@combine_marks( + named_platonic_mark, + named_archimedean_mark, + named_catalan_mark, + named_johnson_mark, + named_prismantiprism_mark, + named_pyramiddipyramid_mark, +) +@given(arrays(np.float64, (3,), elements=floats(-10, 10, width=64), unique=True)) +def test_set_centroid(poly, centroid_vector): + poly.centroid = centroid_vector + coxeter_result = poly.centroid + assert np.allclose(coxeter_result, centroid_vector, atol=1e-12) + poly.centroid = [0, 0, 0] + assert np.allclose(poly.centroid, [0, 0, 0], atol=1e-12) From ffb09d8934f190e69d7a806fe2aa7bc2ddec6515 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 15:43:37 -0400 Subject: [PATCH 38/73] Added in-depth tests for surface area and volume setters --- tests/test_polyhedron.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index e0787dcd..e34afcd6 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -155,6 +155,22 @@ def test_volume_damasceno_shapes(shape): assert np.isclose(poly.volume, hull.volume) +@settings(deadline=1000) +@named_damasceno_shapes_mark +@given(v_test=floats(0, 10)) +def test_set_volume_damasceno_shapes(shape, v_test): + if shape["name"] in ("RESERVED", "Sphere"): + return + if v_test == 0: + return + vertices = shape["vertices"] + poly = ConvexPolyhedron(vertices) + poly.volume = v_test + # Recalculate volume from simplices + calculated_volume = poly._calculate_signed_volume() + assert np.isclose(calculated_volume, v_test) + + @named_damasceno_shapes_mark def test_surface_area_damasceno_shapes(shape): if shape["name"] in ("RESERVED", "Sphere"): @@ -165,6 +181,21 @@ def test_surface_area_damasceno_shapes(shape): assert np.isclose(poly.surface_area, hull.area) +@settings(deadline=1000) +@named_damasceno_shapes_mark +@given(a_test=floats(0, 10)) +def test_set_surface_area_damasceno_shapes(shape, a_test): + if shape["name"] in ("RESERVED", "Sphere"): + return + if a_test == 0: + return + vertices = shape["vertices"] + poly = ConvexPolyhedron(vertices) + poly.surface_area = a_test + calculated_area = poly._calculate_surface_area() + assert np.isclose(calculated_area, a_test) + + @combine_marks( named_platonic_mark, named_archimedean_mark, From 11f5bd3108d5bd3001bf1bc3df6793fda5271bd6 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 16:14:28 -0400 Subject: [PATCH 39/73] Added circleci checker function to speed up local pytest runs --- tests/conftest.py | 10 ++++++++++ tests/test_polyhedron.py | 23 ++++++++++------------- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index 453ef69d..ad482d02 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,6 +1,8 @@ # Copyright (c) 2021 The Regents of the University of Michigan # All rights reserved. # This software is licensed under the BSD 3-Clause License. +import os + import numpy as np import pytest import rowan @@ -37,6 +39,14 @@ def combine_marks(*marks): return combined_mark +# Define a function that checks if the script is running on CircleCI +def is_not_ci_circleci(): + if os.getenv("CI", "false") == "true" or os.getenv("CIRCLECI", "false") == "true": + return False + else: + return True + + # Need to declare this outside the fixture so that it can be used in multiple # fixtures (pytest does not allow fixtures to be called). def get_cube_points(): diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index e34afcd6..e15c9945 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -1,7 +1,6 @@ # Copyright (c) 2021 The Regents of the University of Michigan # All rights reserved. # This software is licensed under the BSD 3-Clause License. -import os import numpy as np import pytest @@ -19,6 +18,7 @@ combine_marks, get_oriented_cube_faces, get_oriented_cube_normals, + is_not_ci_circleci, named_archimedean_mark, named_catalan_mark, named_damasceno_shapes_mark, @@ -155,6 +155,9 @@ def test_volume_damasceno_shapes(shape): assert np.isclose(poly.volume, hull.volume) +@pytest.mark.skipif( + is_not_ci_circleci(), reason="Test is too slow to run during rapid development" +) @settings(deadline=1000) @named_damasceno_shapes_mark @given(v_test=floats(0, 10)) @@ -181,6 +184,9 @@ def test_surface_area_damasceno_shapes(shape): assert np.isclose(poly.surface_area, hull.area) +@pytest.mark.skipif( + is_not_ci_circleci(), reason="Test is too slow to run during rapid development" +) @settings(deadline=1000) @named_damasceno_shapes_mark @given(a_test=floats(0, 10)) @@ -227,12 +233,7 @@ def test_surface_area_shapes(poly): # This test is slow at high precisions. Run fast locally, and test in detail on CircleCI @pytest.mark.parametrize( "atol", - [ - 5e-2 - if os.getenv("CI", "false") == "true" - or os.getenv("CIRCLECI", "false") == "true" - else 1e-1 - ], + [1e-1 if is_not_ci_circleci() else 5e-2], ) @named_damasceno_shapes_mark def test_moment_inertia_damasceno_shapes(shape, atol): @@ -830,14 +831,10 @@ def test_repr_convex(convex_cube): assert str(convex_cube), str(eval(repr(convex_cube))) +# Test fast locally, and in more depth on CircleCI @pytest.mark.parametrize( "atol", - [ - 5e-3 - if os.getenv("CI", "false") == "true" - or os.getenv("CIRCLECI", "false") == "true" - else 1e-2 - ], + [1e-2 if is_not_ci_circleci() else 5e-3], ) @named_damasceno_shapes_mark def test_center(shape, atol): From ac2b10b7c050933f9352aa22d883fca3cc9ff207 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 16:15:34 -0400 Subject: [PATCH 40/73] Added Truncated Icosahedron to bad_shapes for test MOI This shape is individually responsible for ~25% of the total runtime of the test. The test passes (even at high precisions) but it has been marked bad due to the poor performance --- tests/test_polyhedron.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index e15c9945..dc4aed74 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -251,6 +251,7 @@ def test_moment_inertia_damasceno_shapes(shape, atol): "Square Cupola", "Triaugmented Truncated Dodecahedron", "Truncated Dodecahedron", + "Truncated Icosahedron", "Truncated Icosidodecahedron", ] if shape["name"] in ["RESERVED", "Sphere"] + bad_shapes: From 8d22516b2bf178f2e3171724e1e14738d67d8ae5 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 16:29:53 -0400 Subject: [PATCH 41/73] Added pytest for test_face_centroids and fixes for shape creation with Fast=True --- coxeter/shapes/convex_polyhedron.py | 4 ++++ tests/test_polyhedron.py | 11 +++++++++++ 2 files changed, 15 insertions(+) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 389ac93d..8e1c9d3f 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -352,6 +352,8 @@ def get_face_area(self, face=None): self._simplex_areas = self._find_triangle_array_area( self._vertices[self._simplices], sum_result=False ) + if not hasattr(self, "_coplanar_simplices"): + self._find_coplanar_simplices() if face is None: return [ np.sum(self._simplex_areas[self._coplanar_simplices[fac]]) @@ -459,6 +461,8 @@ def faces(self): return self._faces def _find_face_centroids(self): + if not hasattr(self, "_coplanar_simplices"): + self._find_coplanar_simplices() simplex_centroids = np.mean(self._vertices[self._simplices], axis=1) # (N,3) self._simplex_areas = self._find_triangle_array_area( self._vertices[self._simplices], sum_result=False diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index dc4aed74..37e1cd9d 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -876,3 +876,14 @@ def test_set_centroid(poly, centroid_vector): assert np.allclose(coxeter_result, centroid_vector, atol=1e-12) poly.centroid = [0, 0, 0] assert np.allclose(poly.centroid, [0, 0, 0], atol=1e-12) + + +@named_platonic_mark +def test_face_centroids(poly): + # For platonic solids, the centroid of a face is equal to the mean of its vertices + poly.centroid = [0, 0, 0] + coxeter_result = poly.face_centroids + for i, face in enumerate(poly.faces): + face_vertices = poly.vertices[face] + vertex_mean = np.mean(face_vertices, axis=0) + assert np.allclose(vertex_mean, coxeter_result[i]) From eee3cf584442be2a467dccc257aca30ee908d467 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 16:30:49 -0400 Subject: [PATCH 42/73] Set fast flag to properly default to False --- coxeter/shapes/convex_polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 8e1c9d3f..49bfa249 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -94,7 +94,7 @@ class is an extension of :class:`~.Polyhedron` that builds the """ - def __init__(self, vertices, fast=True): + def __init__(self, vertices, fast=False): self._vertices = np.array(vertices, dtype=np.float64) self._ndim = self._vertices.shape[1] self._convex_hull = ConvexHull(self._vertices) From 184ec366eddf3028f1b10c04419ca2063210c900 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 7 Aug 2023 16:33:30 -0400 Subject: [PATCH 43/73] Fix typo in _consume_hull --- coxeter/shapes/convex_polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 49bfa249..f32ea576 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -124,7 +124,7 @@ def _consume_hull(self): """Extract data from ConvexHull. Data is moved from convex hull into private variables. This method deletes the - original hull in order to avoid double storage, and chec. + original hull in order to avoid double storage. """ assert ( self._ndim == self._convex_hull.ndim From f783f67ac2c7ba9c33a0d0d5b447deb8c61a742f Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 8 Aug 2023 10:34:19 -0400 Subject: [PATCH 44/73] Fixed documentation for Polyhedron class The Polyhedron example bases the input off of a ConvexPolyhedron object. Because of this, the ordering of faces can be different. Documentation has been updated to match --- coxeter/shapes/polyhedron.py | 38 ++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 05ab3d34..5548bb2e 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -71,8 +71,8 @@ class Polyhedron(Shape3D): Example: >>> cube = coxeter.shapes.ConvexPolyhedron( - ... [[1, 1, 1], [1, -1, 1], [1, 1, -1], [1, -1, -1], - ... [-1, 1, 1], [-1, -1, 1], [-1, 1, -1], [-1, -1, -1]]) + ... [[-1, -1, -1], [-1, -1, 1], [-1, 1, -1], [-1, 1, 1], + ... [1, -1, -1], [1, -1, 1], [1, 1, -1], [1, 1, 1]]) >>> cube = coxeter.shapes.Polyhedron( ... vertices=cube.vertices, faces=cube.faces) >>> bounding_sphere = cube.minimal_bounding_sphere @@ -81,16 +81,16 @@ class Polyhedron(Shape3D): >>> cube.center array([0., 0., 0.]) >>> cube.faces - [array([4, 5, 1, 0], dtype=int32), array([0, 2, 6, 4], dtype=int32), - array([6, 7, 5, 4], dtype=int32), array([0, 1, 3, 2], dtype=int32), - array([5, 7, 3, 1], dtype=int32), array([2, 3, 7, 6], dtype=int32)] + [array([0, 2, 6, 4], dtype=int32), array([0, 4, 5, 1], dtype=int32), + array([4, 6, 7, 5], dtype=int32), array([0, 1, 3, 2], dtype=int32), + array([2, 3, 7, 6], dtype=int32), array([1, 5, 7, 3], dtype=int32)] >>> cube.gsd_shape_spec {'type': 'Mesh', 'vertices': [[1.0, 1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, -1.0], [1.0, -1.0, -1.0], [-1.0, 1.0, 1.0], [-1.0, -1.0, 1.0], [-1.0, 1.0, -1.0], [-1.0, -1.0, -1.0]], 'indices': - [array([4, 5, 1, 0], dtype=int32), array([0, 2, 6, 4], dtype=int32), - array([6, 7, 5, 4], dtype=int32), array([0, 1, 3, 2], dtype=int32), - array([5, 7, 3, 1], dtype=int32), array([2, 3, 7, 6], dtype=int32)]} + [array([0, 2, 6, 4], dtype=int32), array([0, 4, 5, 1], dtype=int32), + array([4, 6, 7, 5], dtype=int32), array([0, 1, 3, 2], dtype=int32), + array([2, 3, 7, 6], dtype=int32), array([1, 5, 7, 3], dtype=int32)]} >>> assert np.allclose( ... cube.inertia_tensor, ... np.diag([16. / 3., 16. / 3., 16. / 3.])) @@ -99,26 +99,26 @@ class Polyhedron(Shape3D): [array([1, 2, 3, 4]), array([0, 2, 3, 5]), array([0, 1, 4, 5]), array([0, 1, 4, 5]), array([0, 2, 3, 5]), array([1, 2, 3, 4])] >>> cube.normals - array([[ 0., 0., 1.], - [ 0., 1., -0.], - [-1., 0., 0.], - [ 1., -0., 0.], + array([[ 0., 0., -1.], [ 0., -1., 0.], - [ 0., 0., -1.]]) + [ 1., 0., -0.], + [-1., 0., 0.], + [-0., 1., 0.], + [ 0., -0., 1.]]) >>> cube.num_faces 6 >>> cube.num_vertices 8 >>> assert np.isclose(cube.surface_area, 24.0) >>> cube.vertices - array([[ 1., 1., 1.], - [ 1., -1., 1.], - [ 1., 1., -1.], - [ 1., -1., -1.], - [-1., 1., 1.], + array([[-1., -1., -1.], [-1., -1., 1.], [-1., 1., -1.], - [-1., -1., -1.]]) + [-1., 1., 1.], + [ 1., -1., -1.], + [ 1., -1., 1.], + [ 1., 1., -1.], + [ 1., 1., 1.]]) >>> assert np.isclose(cube.volume, 8.0) """ From cfa29566d8a44c90d5e1ec9e863b0221e5555a13 Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 8 Aug 2023 10:37:55 -0400 Subject: [PATCH 45/73] Docstring formatting --- coxeter/shapes/polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 5548bb2e..73e2d43f 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -71,7 +71,7 @@ class Polyhedron(Shape3D): Example: >>> cube = coxeter.shapes.ConvexPolyhedron( - ... [[-1, -1, -1], [-1, -1, 1], [-1, 1, -1], [-1, 1, 1], + ... [[-1, -1, -1], [-1, -1, 1], [-1, 1, -1], [-1, 1, 1], ... [1, -1, -1], [1, -1, 1], [1, 1, -1], [1, 1, 1]]) >>> cube = coxeter.shapes.Polyhedron( ... vertices=cube.vertices, faces=cube.faces) From 83408cfd5d7f6c2c868c450c5fbac5cec543ff4a Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 8 Aug 2023 11:04:58 -0400 Subject: [PATCH 46/73] Docstring formatting fix --- coxeter/shapes/polyhedron.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 73e2d43f..63032302 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -85,9 +85,9 @@ class Polyhedron(Shape3D): array([4, 6, 7, 5], dtype=int32), array([0, 1, 3, 2], dtype=int32), array([2, 3, 7, 6], dtype=int32), array([1, 5, 7, 3], dtype=int32)] >>> cube.gsd_shape_spec - {'type': 'Mesh', 'vertices': [[1.0, 1.0, 1.0], [1.0, -1.0, 1.0], - [1.0, 1.0, -1.0], [1.0, -1.0, -1.0], [-1.0, 1.0, 1.0], - [-1.0, -1.0, 1.0], [-1.0, 1.0, -1.0], [-1.0, -1.0, -1.0]], 'indices': + {'type': 'Mesh', 'vertices': [[-1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], + [-1.0, 1.0, -1.0], [-1.0, 1.0, 1.0], [1.0, -1.0, -1.0], + [1.0, -1.0, 1.0], [1.0, 1.0, -1.0], [1.0, 1.0, 1.0]], 'indices':, 'indices': [array([0, 2, 6, 4], dtype=int32), array([0, 4, 5, 1], dtype=int32), array([4, 6, 7, 5], dtype=int32), array([0, 1, 3, 2], dtype=int32), array([2, 3, 7, 6], dtype=int32), array([1, 5, 7, 3], dtype=int32)]} From 8a61dd5975d159a8524c29f4758a7a2fe889503a Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 8 Aug 2023 11:36:44 -0400 Subject: [PATCH 47/73] Another docstring fix --- coxeter/shapes/polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 63032302..158df5ee 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -87,7 +87,7 @@ class Polyhedron(Shape3D): >>> cube.gsd_shape_spec {'type': 'Mesh', 'vertices': [[-1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], [-1.0, 1.0, -1.0], [-1.0, 1.0, 1.0], [1.0, -1.0, -1.0], - [1.0, -1.0, 1.0], [1.0, 1.0, -1.0], [1.0, 1.0, 1.0]], 'indices':, 'indices': + [1.0, -1.0, 1.0], [1.0, 1.0, -1.0], [1.0, 1.0, 1.0]], 'indices':, [array([0, 2, 6, 4], dtype=int32), array([0, 4, 5, 1], dtype=int32), array([4, 6, 7, 5], dtype=int32), array([0, 1, 3, 2], dtype=int32), array([2, 3, 7, 6], dtype=int32), array([1, 5, 7, 3], dtype=int32)]} From 8df746ae67dd71337917a8e26e6bf341ced1f6a2 Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 8 Aug 2023 12:56:26 -0400 Subject: [PATCH 48/73] One more docstring fix --- coxeter/shapes/polyhedron.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 158df5ee..db732982 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -86,8 +86,8 @@ class Polyhedron(Shape3D): array([2, 3, 7, 6], dtype=int32), array([1, 5, 7, 3], dtype=int32)] >>> cube.gsd_shape_spec {'type': 'Mesh', 'vertices': [[-1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], - [-1.0, 1.0, -1.0], [-1.0, 1.0, 1.0], [1.0, -1.0, -1.0], - [1.0, -1.0, 1.0], [1.0, 1.0, -1.0], [1.0, 1.0, 1.0]], 'indices':, + [-1.0, 1.0, -1.0], [-1.0, 1.0, 1.0], [1.0, -1.0, -1.0], [1.0, -1.0, 1.0], + [1.0, 1.0, -1.0], [1.0, 1.0, 1.0]], 'indices': [array([0, 2, 6, 4], dtype=int32), array([0, 4, 5, 1], dtype=int32), array([4, 6, 7, 5], dtype=int32), array([0, 1, 3, 2], dtype=int32), array([2, 3, 7, 6], dtype=int32), array([1, 5, 7, 3], dtype=int32)]} From f0355e081b925dc70bd50378d9e274a6c08727d6 Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 8 Aug 2023 13:13:41 -0400 Subject: [PATCH 49/73] Added tests for simplex equations, face equations, and normals --- tests/test_polyhedron.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 37e1cd9d..1d8d4f5d 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -887,3 +887,35 @@ def test_face_centroids(poly): face_vertices = poly.vertices[face] vertex_mean = np.mean(face_vertices, axis=0) assert np.allclose(vertex_mean, coxeter_result[i]) + +@given(EllipsoidSurfaceStrategy) +def test_find_simplex_equations(points): + hull = ConvexHull(points) + poly = ConvexPolyhedron(points[hull.vertices]) + # Check simplex equations are stored properly + assert np.allclose(hull.equations, poly._simplex_equations) + + # Now recalculate and check answers are still correct + poly._find_simplex_equations() + assert np.allclose(hull.equations, poly._simplex_equations) + + +@combine_marks( + named_platonic_mark, + named_archimedean_mark, + named_catalan_mark, + named_johnson_mark, + named_prismantiprism_mark, + named_pyramiddipyramid_mark, +) +def test_find_equations_and_normals(poly): + ppoly = Polyhedron(poly.vertices, poly.faces) + # Check face equations are stored properly + assert np.allclose(poly.equations, ppoly._equations) + assert np.allclose(poly.normals, ppoly.normals) + + # Now recalculate and check answers are still correct + poly._find_equations() + ppoly._find_equations() + assert np.allclose(poly.equations, ppoly._equations) + assert np.allclose(poly.normals, ppoly.normals) \ No newline at end of file From c676f18966e0444025192be16b2d09df48c9db4b Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 8 Aug 2023 13:20:20 -0400 Subject: [PATCH 50/73] Combined pytest marks to clean up testing code --- tests/conftest.py | 9 ++++++ tests/test_polyhedron.py | 61 ++++++---------------------------------- 2 files changed, 18 insertions(+), 52 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index ad482d02..3c0ed318 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -337,6 +337,15 @@ def _test_get_set_minimal_bounding_sphere_radius(shape, centered=False): ], ) +named_solids_mark = combine_marks( + named_platonic_mark, + named_archimedean_mark, + named_catalan_mark, + named_johnson_mark, + named_prismantiprism_mark, + named_pyramiddipyramid_mark, +) + data_filenames_mark = pytest.mark.parametrize( argnames="family", argvalues=[ diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 1d8d4f5d..41ba245c 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -22,10 +22,8 @@ named_archimedean_mark, named_catalan_mark, named_damasceno_shapes_mark, - named_johnson_mark, named_platonic_mark, - named_prismantiprism_mark, - named_pyramiddipyramid_mark, + named_solids_mark, sphere_isclose, ) from coxeter.families import DOI_SHAPE_REPOSITORIES, PlatonicFamily @@ -202,28 +200,14 @@ def test_set_surface_area_damasceno_shapes(shape, a_test): assert np.isclose(calculated_area, a_test) -@combine_marks( - named_platonic_mark, - named_archimedean_mark, - named_catalan_mark, - named_johnson_mark, - named_prismantiprism_mark, - named_pyramiddipyramid_mark, -) +@named_solids_mark def test_volume_shapes(poly): vertices = poly.vertices hull = ConvexHull(vertices) assert np.isclose(poly.volume, hull.volume) -@combine_marks( - named_platonic_mark, - named_archimedean_mark, - named_catalan_mark, - named_johnson_mark, - named_prismantiprism_mark, - named_pyramiddipyramid_mark, -) +@named_solids_mark def test_surface_area_shapes(poly): vertices = poly.vertices hull = ConvexHull(vertices) @@ -781,14 +765,7 @@ def test_form_factor(cube): ) -@combine_marks( - named_platonic_mark, - named_archimedean_mark, - named_catalan_mark, - named_johnson_mark, - named_prismantiprism_mark, - named_pyramiddipyramid_mark, -) +@named_solids_mark @pytest.mark.xfail( reason=( "Numerical precision problems with miniball. " @@ -799,14 +776,7 @@ def test_get_set_minimal_bounding_sphere_radius(poly): _test_get_set_minimal_bounding_sphere_radius(poly) -@combine_marks( - named_platonic_mark, - named_archimedean_mark, - named_catalan_mark, - named_johnson_mark, - named_prismantiprism_mark, - named_pyramiddipyramid_mark, -) +@named_solids_mark def test_get_set_minimal_centered_bounding_sphere_radius(poly): _test_get_set_minimal_bounding_sphere_radius(poly, True) @@ -861,14 +831,7 @@ def test_center(shape, atol): ) -@combine_marks( - named_platonic_mark, - named_archimedean_mark, - named_catalan_mark, - named_johnson_mark, - named_prismantiprism_mark, - named_pyramiddipyramid_mark, -) +@named_solids_mark @given(arrays(np.float64, (3,), elements=floats(-10, 10, width=64), unique=True)) def test_set_centroid(poly, centroid_vector): poly.centroid = centroid_vector @@ -888,6 +851,7 @@ def test_face_centroids(poly): vertex_mean = np.mean(face_vertices, axis=0) assert np.allclose(vertex_mean, coxeter_result[i]) + @given(EllipsoidSurfaceStrategy) def test_find_simplex_equations(points): hull = ConvexHull(points) @@ -900,14 +864,7 @@ def test_find_simplex_equations(points): assert np.allclose(hull.equations, poly._simplex_equations) -@combine_marks( - named_platonic_mark, - named_archimedean_mark, - named_catalan_mark, - named_johnson_mark, - named_prismantiprism_mark, - named_pyramiddipyramid_mark, -) +@named_solids_mark def test_find_equations_and_normals(poly): ppoly = Polyhedron(poly.vertices, poly.faces) # Check face equations are stored properly @@ -918,4 +875,4 @@ def test_find_equations_and_normals(poly): poly._find_equations() ppoly._find_equations() assert np.allclose(poly.equations, ppoly._equations) - assert np.allclose(poly.normals, ppoly.normals) \ No newline at end of file + assert np.allclose(poly.normals, ppoly.normals) From 9bb5a10b6c40e55069430163aa8ad9cfbf9a9d50 Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 8 Aug 2023 13:32:16 -0400 Subject: [PATCH 51/73] Update Credits.rst --- Credits.rst | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/Credits.rst b/Credits.rst index a2376e64..fa342638 100644 --- a/Credits.rst +++ b/Credits.rst @@ -112,6 +112,16 @@ 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 optional ``fast`` flag for ConvexPolyhedron generation. +* Reimplemented ``find_equations``, ``_volume``, ``surface_area``, ``centroid``, + ``_compute_inertia_tensor``, ``rescale``, and ``get_face_area`` methods for convex + polyhedra using NumPy vectorized operations and polyhedron simplices. +* Added the ``combine_simplices``, ``find_simplex_equations``, ``_find_face_centroids``, + ``find_coplanar_simplices``, ``_find_face_centroids``, and ``calculate_signed_volume`` + methods to the ConvexPolyhedron class. +* Added ``simplices``, ``equations``, and ``face_centroids`` properties to the + ConvexPolyhedron class. +* Optimized pytest configurations for more efficient use of local and remote resources. Domagoj Fijan * Rewrote point in polygon check to use NumPy vectorized operations. From 0cf79cb9b78ad40b6eaa4fe2cb02183fff490c40 Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 8 Aug 2023 17:01:44 -0400 Subject: [PATCH 52/73] Added test for heavily degenerate input vertices Also checks appropriate warning behavior for such a case --- tests/test_polyhedron.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 41ba245c..7ebe3f8b 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -2,6 +2,8 @@ # All rights reserved. # This software is licensed under the BSD 3-Clause License. +import warnings + import numpy as np import pytest import rowan @@ -876,3 +878,37 @@ def test_find_equations_and_normals(poly): ppoly._find_equations() assert np.allclose(poly.equations, ppoly._equations) assert np.allclose(poly.normals, ppoly.normals) + + +def test_degenerate_cube(convex_cube): + np.random.seed(0) + convex_cube.center = [0, 0, 0] + + # Get a point on the plane for each face + xyz_planes = convex_cube.equations[:, -1][:, None] * convex_cube.equations[:, 0:-1] + + # Extend so there are 10_000 copies of each plane point + xyz_planes = np.repeat(xyz_planes[:, None], 10_000, axis=1) + + # Generate a (10_000,3) array of random points with values ranging from -0.5 to 0.5 + random_points = np.random.rand(*xyz_planes.shape) - 0.5 + + # Move points onto the faces + nonzero_mask = np.where(xyz_planes == 0, 1, 0) + random_points *= nonzero_mask + plane_points = xyz_planes + random_points + + # Add actual cube vertices to the list of degenerate points + degenerate_points = convex_cube.vertices + degenerate_points = np.vstack([degenerate_points, plane_points.reshape(-1, 3)]) + + with warnings.catch_warnings(record=True) as caught: + degenerate_cube = ConvexPolyhedron(degenerate_points) + + # Should catch one warning - the one raised when degenerate points are pruned + assert len(caught) == 1 + + assert degenerate_cube.num_faces == 6 + assert degenerate_cube.num_edges == 12 + assert degenerate_cube.num_vertices == 8 + assert len(degenerate_cube.simplices) == 12 From 83e664214ae8b1d46d50e1e94d60e844ba39b8a7 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 18 Sep 2023 14:37:34 -0400 Subject: [PATCH 53/73] Changed [test] to [tests] --- .circleci/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index 85c55cac..410b853b 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -64,7 +64,7 @@ jobs: - run: name: Install wheel command: - python -m pip install coxeter[test] --progress-bar off --user -U --force-reinstall -f dist/ + python -m pip install coxeter[tests] --progress-bar off --user -U --force-reinstall -f dist/ - run: <<: *run-tests From 39490bbcf6fede4c8be38328ebe76e319c46fea0 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 23 Oct 2023 11:40:11 -0400 Subject: [PATCH 54/73] Updated ConvexPolyhedron to properly handle nonconvex inputs. --- coxeter/shapes/convex_polyhedron.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index d3b5f1c7..4827a9a6 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -101,10 +101,15 @@ def __init__(self, vertices, fast=False): self._faces_are_convex = True if not len(self._convex_hull.vertices) == len(self._vertices): - warnings.warn("Not all vertices used for hull.", RuntimeWarning) + # Identify vertices that do not contribute to the convex hull + nonconvex_vertices = sorted( + set(range(len(self._vertices))) - set(self._convex_hull.vertices) + ) - # Remove internal points - self._vertices = self._vertices[self._convex_hull.vertices] + raise ValueError( + "Input vertices must be a convex set. " + + f"Vertices {nonconvex_vertices} are inside the shape (or coplanar)." + ) # Transfer data in from convex hull, then clean up the results. self._consume_hull() From 2fbe7bf05a676de3b53415b814e88323f28b8921 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 23 Oct 2023 11:43:52 -0400 Subject: [PATCH 55/73] Removed use of Fast flag from class --- coxeter/shapes/convex_polyhedron.py | 24 +++--------------------- 1 file changed, 3 insertions(+), 21 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 4827a9a6..5ebb75c8 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -29,12 +29,6 @@ class is an extension of :class:`~.Polyhedron` that builds the Args: vertices (:math:`(N, 3)` :class:`numpy.ndarray`): The vertices of the polyhedron. - fast (bool, optional): - Creation mode for the polyhedron. fast=False (default) will perform all - checks and precalculate most properties. fast=True will precalculate some - properties, but will not find face neighbors or sort face indices. These - calculations will instead be performed when required by another method. - (Default value: False). Example: >>> cube = coxeter.shapes.ConvexPolyhedron( @@ -94,7 +88,7 @@ class is an extension of :class:`~.Polyhedron` that builds the """ - def __init__(self, vertices, fast=False): + def __init__(self, vertices): self._vertices = np.array(vertices, dtype=np.float64) self._ndim = self._vertices.shape[1] self._convex_hull = ConvexHull(self._vertices) @@ -117,13 +111,8 @@ def __init__(self, vertices, fast=False): # Sort simplices. This method also calculates simplex equations and the centroid self._sort_simplices() - - # If mode is NOT fast, perform additional initializiation steps - if fast is False: - self._find_coplanar_simplices() - self.sort_faces() - else: - self._faces_are_sorted = False + self._find_coplanar_simplices() + self.sort_faces() def _consume_hull(self): """Extract data from ConvexHull. @@ -461,8 +450,6 @@ def _centroid_from_triangulated_surface(self): @property def faces(self): """list(:class:`numpy.ndarray`): Get the polyhedron's faces.""" - if self._faces_are_sorted is False: - self.sort_faces() return self._faces def _find_face_centroids(self): @@ -502,8 +489,6 @@ def neighbors(self): element is an array of indices of faces that are neighbors of face :math:`i`. """ - if self._faces_are_sorted is False: - self.sort_faces() return self._neighbors def sort_faces(self): @@ -511,7 +496,6 @@ def sort_faces(self): This does NOT change the *order* of faces in the list. """ - self._faces_are_sorted = True # Get correct-quadrant angles about the face normal sorted_faces = [] @@ -837,8 +821,6 @@ def get_dihedral(self, a, b): >>> assert np.isclose(cube.get_dihedral(1, 2), np.pi / 2.) """ - if self._faces_are_sorted is False: - self.sort_faces() if b not in self.neighbors[a]: raise ValueError("The two faces are not neighbors.") n1, n2 = self._equations[[a, b], :3] From f3561a58f48398083667602035cfadbc2a2eace5 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 23 Oct 2023 11:44:36 -0400 Subject: [PATCH 56/73] Removed fast flag from credits --- Credits.rst | 1 - 1 file changed, 1 deletion(-) diff --git a/Credits.rst b/Credits.rst index 1cadc7df..955cee02 100644 --- a/Credits.rst +++ b/Credits.rst @@ -113,7 +113,6 @@ Jen Bradley * Added shape family for prisms and antiprisms. * Added shape family for equilateral pyramids and dipyramids. * Added edges, edge_vectors, and num_edges methods. -* Added optional ``fast`` flag for ConvexPolyhedron generation. * Reimplemented ``find_equations``, ``_volume``, ``surface_area``, ``centroid``, ``_compute_inertia_tensor``, ``rescale``, and ``get_face_area`` methods for convex polyhedra using NumPy vectorized operations and polyhedron simplices. From 266a15f9753a6527fc95916bd2d30d7c2b7a0b91 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 23 Oct 2023 11:57:07 -0400 Subject: [PATCH 57/73] Updated ci flags for pytest --- tests/conftest.py | 9 ++++++--- tests/test_polyhedron.py | 18 ++++++++++-------- 2 files changed, 16 insertions(+), 11 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index 3c0ed318..abad3ab6 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -39,9 +39,12 @@ def combine_marks(*marks): return combined_mark -# Define a function that checks if the script is running on CircleCI -def is_not_ci_circleci(): - if os.getenv("CI", "false") == "true" or os.getenv("CIRCLECI", "false") == "true": +# Define a function that checks if the script is running on Github Actions +def is_not_ci(): + if ( + os.getenv("CI", "false") == "true" + or os.getenv("GITHUB_ACTIONS", "false") == "true" + ): return False else: return True diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index a6da0e73..21f77da7 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -20,7 +20,7 @@ combine_marks, get_oriented_cube_faces, get_oriented_cube_normals, - is_not_ci_circleci, + is_not_ci, named_archimedean_mark, named_catalan_mark, named_damasceno_shapes_mark, @@ -156,7 +156,7 @@ def test_volume_damasceno_shapes(shape): @pytest.mark.skipif( - is_not_ci_circleci(), reason="Test is too slow to run during rapid development" + is_not_ci(), reason="Test is too slow to run during rapid development" ) @settings(deadline=1000) @named_damasceno_shapes_mark @@ -185,7 +185,7 @@ def test_surface_area_damasceno_shapes(shape): @pytest.mark.skipif( - is_not_ci_circleci(), reason="Test is too slow to run during rapid development" + is_not_ci(), reason="Test is too slow to run during rapid development" ) @settings(deadline=1000) @named_damasceno_shapes_mark @@ -216,10 +216,11 @@ def test_surface_area_shapes(poly): assert np.isclose(poly.surface_area, hull.area) -# This test is slow at high precisions. Run fast locally, and test in detail on CircleCI -@pytest.mark.parametrize( - "atol", - [1e-1 if is_not_ci_circleci() else 5e-2], +# This test is a bit slow (a couple of minutes), so skip running it locally. +@pytest.mark.skipif( + os.getenv("CI", "false") != "true" + and os.getenv("GITHUB_ACTIONS", "false") != "true", + reason="Test is too slow to run during rapid development", ) @named_damasceno_shapes_mark def test_moment_inertia_damasceno_shapes(shape, atol): @@ -412,6 +413,7 @@ def test_edge_lengths(): poly.vertices[poly.edges[:, 1]] - poly.vertices[poly.edges[:, 0]], axis=1 ) assert np.allclose(veclens, edgelength) + assert np.allclose(poly.edge_lengths, edgelength) assert np.allclose(veclens, np.linalg.norm(poly.edge_vectors, axis=1)) @@ -860,7 +862,7 @@ def test_repr_convex(convex_cube): # Test fast locally, and in more depth on CircleCI @pytest.mark.parametrize( "atol", - [1e-2 if is_not_ci_circleci() else 5e-3], + [1e-2 if is_not_ci() else 5e-3], ) @named_damasceno_shapes_mark def test_center(shape, atol): From 23244fd808849c0da0c68ceb7be9192b8db2bf14 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 23 Oct 2023 11:59:06 -0400 Subject: [PATCH 58/73] Fixed D202 in convex_polyhedron --- coxeter/shapes/convex_polyhedron.py | 1 - 1 file changed, 1 deletion(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 5ebb75c8..69056186 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -496,7 +496,6 @@ def sort_faces(self): This does NOT change the *order* of faces in the list. """ - # Get correct-quadrant angles about the face normal sorted_faces = [] From 3e416404aee03d0202151aff1fb52dc516213c73 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 23 Oct 2023 12:04:27 -0400 Subject: [PATCH 59/73] Added ci check to MI test --- tests/test_polyhedron.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 21f77da7..275d1a6a 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -218,9 +218,7 @@ def test_surface_area_shapes(poly): # This test is a bit slow (a couple of minutes), so skip running it locally. @pytest.mark.skipif( - os.getenv("CI", "false") != "true" - and os.getenv("GITHUB_ACTIONS", "false") != "true", - reason="Test is too slow to run during rapid development", + is_not_ci(), reason="Test is too slow to run during rapid development" ) @named_damasceno_shapes_mark def test_moment_inertia_damasceno_shapes(shape, atol): From 7cdcdf78ed6924a87d43eae0b6681d202f2e6bf4 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 23 Oct 2023 16:11:22 +0000 Subject: [PATCH 60/73] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- coxeter/shapes/convex_polyhedron.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 69056186..4f50fe7e 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -312,7 +312,8 @@ def _calculate_signed_volume(self): volumes. The external volume property will always be a positive value, but accessing the signed volume can be useful for some mathematical operations. - Returns: + Returns + ------- float: Signed volume of the polyhedron. """ signed_volume = np.sum(np.linalg.det(self._vertices[self._simplices]) / 6) @@ -328,7 +329,8 @@ def get_face_area(self, face=None): find the area. If None, finds the area of all faces. (Default value: None). - Returns: + Returns + ------- :class:`numpy.ndarray`: The area of each face. Example: @@ -474,7 +476,8 @@ def _find_face_centroids(self): def face_centroids(self): """Calculate the centroid (center of mass) of each polygonal face. - Returns: + Returns + ------- :math:`(N,3)` :class:`numpy.ndarray`: Array of centroids for each face. """ @@ -528,7 +531,8 @@ def sort_faces(self): def _surface_triangulation(self): """Output the vertices of simplices composing the polyhedron's surface. - Returns: + Returns + ------- :math:`(N,3,3)` :class:`numpy.ndarray`: Array of vertices for simplices composing the polyhedron's surface. """ @@ -538,7 +542,8 @@ def _surface_triangulation(self): def simplices(self): """Output the vertex indices of simplices composing the polyhedron's surface. - Returns: + Returns + ------- :math:`(N,3)` :class:`numpy.ndarray`: Array of vertex indices of simplices making up the polyhedron's surface. """ @@ -555,7 +560,8 @@ def _find_triangle_array_area(self, triangle_vertices, sum_result=True): Whether the output should be summed. (Default value: True). - Returns: + Returns + ------- :math:`(N, )` :class:`numpy.ndarray` or float: Boolean array indicating which points are contained in thepolyhedron. If sum_result is True, a single value is returned. @@ -807,7 +813,8 @@ def get_dihedral(self, a, b): b (int): The index of the second face. - Returns: + Returns + ------- float: The dihedral angle in radians. Example: From ad30552fc4c279ee142c4e6221a9c0d4e404eefe Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 23 Oct 2023 16:35:14 -0400 Subject: [PATCH 61/73] Updated tests --- tests/test_polyhedron.py | 25 ++++++------------------- 1 file changed, 6 insertions(+), 19 deletions(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 275d1a6a..5c2946dc 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -155,17 +155,12 @@ def test_volume_damasceno_shapes(shape): assert np.isclose(poly.volume, hull.volume) -@pytest.mark.skipif( - is_not_ci(), reason="Test is too slow to run during rapid development" -) -@settings(deadline=1000) +@settings(max_examples=10 if is_not_ci() else 50) @named_damasceno_shapes_mark -@given(v_test=floats(0, 10)) +@given(v_test=floats(0, 10, exclude_min=True)) def test_set_volume_damasceno_shapes(shape, v_test): if shape["name"] in ("RESERVED", "Sphere"): return - if v_test == 0: - return vertices = shape["vertices"] poly = ConvexPolyhedron(vertices) poly.volume = v_test @@ -184,17 +179,12 @@ def test_surface_area_damasceno_shapes(shape): assert np.isclose(poly.surface_area, hull.area) -@pytest.mark.skipif( - is_not_ci(), reason="Test is too slow to run during rapid development" -) -@settings(deadline=1000) +@settings(max_examples=10 if is_not_ci() else 50) @named_damasceno_shapes_mark -@given(a_test=floats(0, 10)) +@given(a_test=floats(0, 10, exclude_min=True)) def test_set_surface_area_damasceno_shapes(shape, a_test): if shape["name"] in ("RESERVED", "Sphere"): return - if a_test == 0: - return vertices = shape["vertices"] poly = ConvexPolyhedron(vertices) poly.surface_area = a_test @@ -216,12 +206,9 @@ def test_surface_area_shapes(poly): assert np.isclose(poly.surface_area, hull.area) -# This test is a bit slow (a couple of minutes), so skip running it locally. -@pytest.mark.skipif( - is_not_ci(), reason="Test is too slow to run during rapid development" -) @named_damasceno_shapes_mark -def test_moment_inertia_damasceno_shapes(shape, atol): +def test_moment_inertia_damasceno_shapes(shape, atol=1e-1): + # Values of atol up to 5e-2 work as expected, but take much longer to run. # These shapes pass the test for a sufficiently high number of samples, but # the number is too high to be worth running them regularly. bad_shapes = [ From 815e7313bcea0dbee4ec2c0fab1de223b338d348 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 23 Oct 2023 16:36:38 -0400 Subject: [PATCH 62/73] Updated docstring for _calculate_signed_volume --- coxeter/shapes/convex_polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 3965e5e6..d642fa45 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -306,7 +306,7 @@ def _rescale(self, scale_factor: Number): self._centroid_from_triangulated_surface() def _calculate_signed_volume(self): - """Internal method to calculate the signed volume of the polyhedron. + """Calculate the signed volume of the polyhedron. This class splits the shape into tetrahedra, then sums their contributing volumes. The external volume property will always be a positive value, but From 23de270a38f6fde724904808bcf356c69bf5f44c Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 23 Oct 2023 17:14:18 -0400 Subject: [PATCH 63/73] Removed degenerate cube test --- tests/test_polyhedron.py | 35 ----------------------------------- 1 file changed, 35 deletions(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 5c2946dc..c9cfac48 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -2,7 +2,6 @@ # All rights reserved. # This software is licensed under the BSD 3-Clause License. -import warnings import numpy as np import pytest @@ -918,37 +917,3 @@ def test_find_equations_and_normals(poly): ppoly._find_equations() assert np.allclose(poly.equations, ppoly._equations) assert np.allclose(poly.normals, ppoly.normals) - - -def test_degenerate_cube(convex_cube): - np.random.seed(0) - convex_cube.center = [0, 0, 0] - - # Get a point on the plane for each face - xyz_planes = convex_cube.equations[:, -1][:, None] * convex_cube.equations[:, 0:-1] - - # Extend so there are 10_000 copies of each plane point - xyz_planes = np.repeat(xyz_planes[:, None], 10_000, axis=1) - - # Generate a (10_000,3) array of random points with values ranging from -0.5 to 0.5 - random_points = np.random.rand(*xyz_planes.shape) - 0.5 - - # Move points onto the faces - nonzero_mask = np.where(xyz_planes == 0, 1, 0) - random_points *= nonzero_mask - plane_points = xyz_planes + random_points - - # Add actual cube vertices to the list of degenerate points - degenerate_points = convex_cube.vertices - degenerate_points = np.vstack([degenerate_points, plane_points.reshape(-1, 3)]) - - with warnings.catch_warnings(record=True) as caught: - degenerate_cube = ConvexPolyhedron(degenerate_points) - - # Should catch one warning - the one raised when degenerate points are pruned - assert len(caught) == 1 - - assert degenerate_cube.num_faces == 6 - assert degenerate_cube.num_edges == 12 - assert degenerate_cube.num_vertices == 8 - assert len(degenerate_cube.simplices) == 12 From 5a13055bfb572e05709eb2b5381862aa4a68a934 Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 24 Oct 2023 15:47:39 -0400 Subject: [PATCH 64/73] Removed remaining circleci files --- .circleci/config.yml | 93 -------------------------------------------- 1 file changed, 93 deletions(-) delete mode 100644 .circleci/config.yml diff --git a/.circleci/config.yml b/.circleci/config.yml deleted file mode 100644 index 410b853b..00000000 --- a/.circleci/config.yml +++ /dev/null @@ -1,93 +0,0 @@ -# Python CircleCI 2.0 configuration file -# -# Check https://circleci.com/docs/2.0/language-python/ for more details -# -version: 2.1 - -orbs: - codecov: codecov/codecov@3.2.2 - -jobs: - test-310: &test-template - docker: - - image: cimg/python:3.10 - - working_directory: ~/repo - - steps: - - - checkout - - - run: &install - name: Install package - command: | - python -m pip install .[tests] --progress-bar off --user - - - run: &run-tests - name: Run tests - command: | - python --version - python -c "import numpy; print('numpy', numpy.__version__)" - python -c "import scipy; print('scipy', scipy.__version__)" - python -m pytest -v --cov=coxeter/ --cov-report=xml - - test-39: - <<: *test-template - docker: - - image: cimg/python:3.9 - - test-38: - <<: *test-template - docker: - - image: cimg/python:3.8 - - pypi_wheels: &pypi_wheels - docker: - - image: cimg/python:3.10 - - working_directory: ~/repo - - steps: - - checkout - - - run: - name: Create source distribution and wheels - command: | - python --version - python -m pip --version - python -m pip install --progress-bar off --user -U twine wheel setuptools - python -m twine --version - python -m wheel version - python setup.py sdist - python setup.py bdist_wheel - - - run: - name: Install wheel - command: - python -m pip install coxeter[tests] --progress-bar off --user -U --force-reinstall -f dist/ - - - run: - <<: *run-tests - - - run: - name: Upload source distribution and wheels - command: | - python -m twine upload --username ${PYPI_USERNAME} --password ${PYPI_PASSWORD} dist/* - -workflows: - version: 2 - test: - jobs: - - test-310: - post-steps: - - codecov/upload - - test-39 - - test-38 - deploy: - jobs: - - pypi_wheels: - filters: - tags: - only: /^v.*/ - branches: - ignore: /.*/ From 261e3da256b8952fa73eec29ec651150245a322f Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 24 Oct 2023 16:02:08 -0400 Subject: [PATCH 65/73] Began updating changelog --- ChangeLog.rst | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/ChangeLog.rst b/ChangeLog.rst index a4d03401..d330b9d9 100644 --- a/ChangeLog.rst +++ b/ChangeLog.rst @@ -7,6 +7,11 @@ Added ~~~~~ - New `edge_lengths` method. +- ``combine_simplices``, ``find_simplex_equations``, ``_find_face_centroids``, + ``find_coplanar_simplices``, ``_find_face_centroids``, and ``calculate_signed_volume`` + methods for the ConvexPolyhedron class. +- ``simplices``, ``equations``, and ``face_centroids`` properties for the + ConvexPolyhedron class. Changed ~~~~~~~ @@ -14,6 +19,15 @@ Changed - Pre-commit now uses ruff instead of flake8, pydocstyle, pyupgrade and isort. - CI now uses GitHub Actions. - Docs ported to furo theme. +- Reimplemented ``find_equations``, ``_volume``, ``surface_area``, ``centroid``, + ``_compute_inertia_tensor``, ``rescale``, and ``get_face_area`` methods for convex + polyhedra using NumPy vectorized operations and polyhedron simplices. +- [breaking] ``ConvexPolyhedron._surface_triangulation`` now returns sorted simplices, + rather than running polytri. This can change the order of vertices and/or triangles. +- [breaking] ``sort_faces`` may return faces in a different order than previously. Sorted faces will still be sorted such that curl and divergence theorems work properly. +- ``volume``, ``surface_area``, and ``centroid`` properties now return stored values, rather than computing the quantity at each call. +- ``rescale`` now computes the centroid to ensure the correct value is available when ``centroid`` is called. +- Optimized pytest configurations for more efficient use of local and remote resources. v0.7.0 - 2023-09-18 ------------------- From 16505be286bc85aafbcd4654c3a991ea50978285 Mon Sep 17 00:00:00 2001 From: janbridley Date: Tue, 24 Oct 2023 16:04:01 -0400 Subject: [PATCH 66/73] Updating changelog --- ChangeLog.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/ChangeLog.rst b/ChangeLog.rst index d330b9d9..133066d3 100644 --- a/ChangeLog.rst +++ b/ChangeLog.rst @@ -12,6 +12,7 @@ Added methods for the ConvexPolyhedron class. - ``simplices``, ``equations``, and ``face_centroids`` properties for the ConvexPolyhedron class. +- Additional pytests for surface area, volume, centroid, moment of inertia, and equations properties. Changed ~~~~~~~ From d62154afb432cd0066e9d7aa607740b60e00f4e2 Mon Sep 17 00:00:00 2001 From: janbridley Date: Wed, 8 Nov 2023 10:03:59 -0500 Subject: [PATCH 67/73] Clarified language on faces breaking change --- ChangeLog.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ChangeLog.rst b/ChangeLog.rst index 133066d3..d2e70122 100644 --- a/ChangeLog.rst +++ b/ChangeLog.rst @@ -25,7 +25,7 @@ Changed polyhedra using NumPy vectorized operations and polyhedron simplices. - [breaking] ``ConvexPolyhedron._surface_triangulation`` now returns sorted simplices, rather than running polytri. This can change the order of vertices and/or triangles. -- [breaking] ``sort_faces`` may return faces in a different order than previously. Sorted faces will still be sorted such that curl and divergence theorems work properly. +- [breaking] ``faces`` may return faces in a different order than previously. Faces are still sorted with ``sort_faces``, and will still be ordered such that curl and divergence theorems work properly. - ``volume``, ``surface_area``, and ``centroid`` properties now return stored values, rather than computing the quantity at each call. - ``rescale`` now computes the centroid to ensure the correct value is available when ``centroid`` is called. - Optimized pytest configurations for more efficient use of local and remote resources. From 571d5b51f37844eb5def545aa71285f91b24bf88 Mon Sep 17 00:00:00 2001 From: janbridley Date: Wed, 8 Nov 2023 14:07:09 -0500 Subject: [PATCH 68/73] Cleaned up convex hull handling --- coxeter/shapes/convex_polyhedron.py | 32 ++++++++++++----------------- 1 file changed, 13 insertions(+), 19 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index d642fa45..4ff66cbb 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -91,13 +91,13 @@ class is an extension of :class:`~.Polyhedron` that builds the def __init__(self, vertices): self._vertices = np.array(vertices, dtype=np.float64) self._ndim = self._vertices.shape[1] - self._convex_hull = ConvexHull(self._vertices) + hull = ConvexHull(self._vertices) self._faces_are_convex = True - if not len(self._convex_hull.vertices) == len(self._vertices): + if not len(hull.vertices) == len(self._vertices): # Identify vertices that do not contribute to the convex hull nonconvex_vertices = sorted( - set(range(len(self._vertices))) - set(self._convex_hull.vertices) + set(range(len(self._vertices))) - set(hull.vertices) ) raise ValueError( @@ -106,7 +106,7 @@ def __init__(self, vertices): ) # Transfer data in from convex hull, then clean up the results. - self._consume_hull() + self._consume_hull(hull) self._combine_simplices() # Sort simplices. This method also calculates simplex equations and the centroid @@ -114,27 +114,21 @@ def __init__(self, vertices): self._find_coplanar_simplices() self.sort_faces() - def _consume_hull(self): + def _consume_hull(self, hull): """Extract data from ConvexHull. - Data is moved from convex hull into private variables. This method deletes the - original hull in order to avoid double storage. + Data is moved from convex hull into private variables. """ assert ( - self._ndim == self._convex_hull.ndim + self._ndim == hull.ndim ), "Input points are coplanar or close to coplanar." - self._simplices = self._convex_hull.simplices[:] - self._simplex_equations = self._convex_hull.equations[:] - self._simplex_neighbors = self._convex_hull.neighbors[:] - self._volume = self._convex_hull.volume - self._area = self._convex_hull.area - self._maximal_extents = np.array( - [self._convex_hull.min_bound, self._convex_hull.max_bound] - ) - - # Clean up the result. - del self._convex_hull + self._simplices = hull.simplices[:] + self._simplex_equations = hull.equations[:] + self._simplex_neighbors = hull.neighbors[:] + self._volume = hull.volume + self._area = hull.area + self._maximal_extents = np.array([hull.min_bound, hull.max_bound]) def _combine_simplices(self, rounding: int = 15): """Combine simplices into faces, merging based on simplex equations. From 615231362939740535f472a32a1069bb3777aec5 Mon Sep 17 00:00:00 2001 From: Jen Bradley <55467578+janbridley@users.noreply.github.com> Date: Wed, 8 Nov 2023 14:34:52 -0500 Subject: [PATCH 69/73] Update coxeter/shapes/convex_polyhedron.py Co-authored-by: Domagoj Fijan <50439291+DomFijan@users.noreply.github.com> --- coxeter/shapes/convex_polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 4ff66cbb..cef7081b 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -185,7 +185,7 @@ def _find_coplanar_simplices(self, rounding: int = 15): self._coplanar_simplices = ragged_coplanar_indices def _sort_simplices(self): - """Reorder simplices counterclockwise relatative to the plane they lie on. + """Reorder simplices counterclockwise relative to the plane they lie on. This does NOT change the *order* of simplices in the list. """ From 9aa4a59d2422f11edd8392d3ca036be33820a1bb Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 10 Nov 2023 10:35:12 -0500 Subject: [PATCH 70/73] Removed unnecessary call to principal_moments --- coxeter/shapes/convex_polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 4ff66cbb..7565447d 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -683,7 +683,7 @@ def diagonalize_inertia(self): [-1., -1., -1.]]) """ - principal_moments, principal_axes = np.linalg.eigh(self.inertia_tensor) + _, principal_axes = np.linalg.eigh(self.inertia_tensor) self._vertices = np.dot(self._vertices, principal_axes) self._sort_simplices() From 1ed3735ed15f092c7000da27a5f93045411ba59d Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 10 Nov 2023 11:00:03 -0500 Subject: [PATCH 71/73] Fixed pytest fixture call --- tests/test_polyhedron.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index c9cfac48..42374cae 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -349,9 +349,7 @@ def test___repr__(): repr(icosidodecahedron) -@combine_marks( - named_solids_mark, -) +@named_solids_mark def test_edges(poly): # Check that the first column is in ascending order. assert np.all(np.diff(poly.edges[:, 0]) >= 0) From 4fe9d8846b26b755bc8adbdbb66828c49513404e Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 4 Dec 2023 12:26:08 -0500 Subject: [PATCH 72/73] Fixed bug in _find_coplanar_simplices --- coxeter/shapes/convex_polyhedron.py | 79 ++++++++++------------------- tests/test_polyhedron.py | 19 +++++++ 2 files changed, 47 insertions(+), 51 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index ffc0f511..d2ab4333 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -4,7 +4,6 @@ """Defines a convex polyhedron.""" import warnings -from collections import defaultdict from numbers import Number import numpy as np @@ -111,7 +110,6 @@ def __init__(self, vertices): # Sort simplices. This method also calculates simplex equations and the centroid self._sort_simplices() - self._find_coplanar_simplices() self.sort_faces() def _consume_hull(self, hull): @@ -130,59 +128,41 @@ def _consume_hull(self, hull): self._area = hull.area self._maximal_extents = np.array([hull.min_bound, hull.max_bound]) - def _combine_simplices(self, rounding: int = 15): + def _combine_simplices(self, tol: float = 2e-15): """Combine simplices into faces, merging based on simplex equations. Coplanar faces will have identical equations (within rounding tolerance). Values - should be slightly larger than machine epsilon (e.g. rounding=15 for ~1e-15) + should be about an order of magnitude greater than machine epsilon. Args: - rounding (int, optional): - Integer number of decimal places to round equations to. - (Default value: 15). + tol (float, optional): + Floating point tolerance within which values are considered identical. + (Default value: 2e-15). """ - equation_groups = defaultdict(list) + is_coplanar = np.all( + np.abs(self._simplex_equations[:, None] - self._simplex_equations) < tol, + axis=2, + ) + coplanar_indices = [[] for _ in range(len(self._simplices))] - # Iterate over all simplex equations - for i, equation in enumerate(self._simplex_equations): - # Convert to hashable key - equation_key = tuple(equation.round(rounding)) + # Iterate over coplanar indices to build face index lists + for face, index in zip(*is_coplanar.nonzero()): + coplanar_indices[face].append(index) - # Store vertex indices from the new simplex under the correct equation key - equation_groups[equation_key].extend(self._simplices[i]) + # Remove duplicate faces, then sort the face indices by their minimum value + coplanar_indices = sorted(set(map(tuple, coplanar_indices)), key=lambda x: x[0]) - # Combine elements with the same plan equation and remove duplicate indices - ragged_faces = [ - np.fromiter(set(group), np.int32) for group in equation_groups.values() - ] - self._faces = ragged_faces - self._equations = np.array(list(equation_groups.keys())) + # Extract vertex indices from simplex indices and remove duplicates + faces = [np.unique(self.simplices[[ind]]) for ind in coplanar_indices] + self._faces = faces - def _find_coplanar_simplices(self, rounding: int = 15): - """ - Get lists of simplex indices for coplanar simplices. - - Args: - rounding (int, optional): - Integer number of decimal places to round equations to. - (Default value: 15). - - - """ - # Combine simplex indices - equation_groups = defaultdict(list) - - # Iterate over all simplex equations - for i, equation in enumerate(self._simplex_equations): - # Convert equation into hashable tuple - equation_key = tuple(equation.round(rounding)) - equation_groups[equation_key].append(i) - ragged_coplanar_indices = [ - np.fromiter(set(group), np.int32) for group in equation_groups.values() + # Copy the simplex equation for one of the simplices on each face + self._equations = self._simplex_equations[ + [equation_index[0] for equation_index in coplanar_indices] ] - - self._coplanar_simplices = ragged_coplanar_indices + # Convert the simplex indices to numpy arrays and save + self._coplanar_simplices = list(map(np.array, coplanar_indices)) def _sort_simplices(self): """Reorder simplices counterclockwise relative to the plane they lie on. @@ -342,24 +322,23 @@ def get_face_area(self, face=None): self._simplex_areas = self._find_triangle_array_area( self._vertices[self._simplices], sum_result=False ) - if not hasattr(self, "_coplanar_simplices"): - self._find_coplanar_simplices() if face is None: + # Return face area for every face. return [ np.sum(self._simplex_areas[self._coplanar_simplices[fac]]) for fac in range(self.num_faces) ] - elif face == "total": # return total surface area + elif face == "total": + # Return total surface area return np.sum(self._simplex_areas) elif hasattr(face, "__len__"): - # For list of input vertices - # Combine coplanar simplices + # Return face areas for a list of faces return [ np.sum(self._simplex_areas[self._coplanar_simplices[fac]]) for fac in face ] else: - # For integer input (single face) + # Return face area of a single face return np.sum(self._simplex_areas[self._coplanar_simplices[face]]) def _find_equations(self): @@ -449,8 +428,6 @@ def faces(self): return self._faces def _find_face_centroids(self): - if not hasattr(self, "_coplanar_simplices"): - self._find_coplanar_simplices() simplex_centroids = np.mean(self._vertices[self._simplices], axis=1) # (N,3) self._simplex_areas = self._find_triangle_array_area( self._vertices[self._simplices], sum_result=False diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index 42374cae..25cd2a8b 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -205,6 +205,25 @@ def test_surface_area_shapes(poly): assert np.isclose(poly.surface_area, hull.area) +@named_solids_mark +def test_surface_area_per_face(poly): + # Sum over all simplices + total_area = poly.get_face_area(face="total") + assert np.isclose(poly.surface_area, total_area) + + # Compute per-face areas and check + face_areas = poly.get_face_area(face=None) + total_face_area = np.sum(face_areas) + assert np.isclose(poly.surface_area, total_face_area) + + # Compute areas of each face and check that ordering is correct + for face_number in range(poly.num_faces): + current_face_area = poly.get_face_area(face=face_number) + assert current_face_area == face_areas[face_number] + total_face_area -= current_face_area + assert np.isclose(total_face_area, 0) + + @named_damasceno_shapes_mark def test_moment_inertia_damasceno_shapes(shape, atol=1e-1): # Values of atol up to 5e-2 work as expected, but take much longer to run. From 778f69f50b566876928d48f5efaf4361a21694f8 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 4 Dec 2023 12:36:26 -0500 Subject: [PATCH 73/73] Removed copy of inherited method get_dihedral --- coxeter/shapes/convex_polyhedron.py | 30 ----------------------------- 1 file changed, 30 deletions(-) diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index d2ab4333..bbfcd6bc 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -776,36 +776,6 @@ def maximal_centered_bounded_sphere(self): min_distance = -np.max(distances) return Sphere(min_distance, center) - def get_dihedral(self, a, b): - """Get the dihedral angle between a pair of faces. - - The dihedral is computed from the dot product of the face normals. - - Args: - a (int): - The index of the first face. - b (int): - The index of the second face. - - Returns - ------- - float: The dihedral angle in radians. - - Example: - >>> cube = coxeter.shapes.ConvexPolyhedron( - ... [[1, 1, 1], [1, -1, 1], [1, 1, -1], [1, -1, -1], - ... [-1, 1, 1], [-1, -1, 1], [-1, 1, -1], [-1, -1, -1]]) - >>> cube = coxeter.shapes.Polyhedron( - ... vertices=cube.vertices, faces=cube.faces) - >>> import numpy as np - >>> assert np.isclose(cube.get_dihedral(1, 2), np.pi / 2.) - - """ - if b not in self.neighbors[a]: - raise ValueError("The two faces are not neighbors.") - n1, n2 = self._equations[[a, b], :3] - return np.arccos(np.dot(-n1, n2)) - def _plato_primitive(self, backend): return backend.ConvexPolyhedra( positions=np.array([self.center]),