Added new classes to return the edges of polyhedra

[janbridley]

Description

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.

Types of Changes

• Documentation update
• Bug fix
• New feature

Checklist:

• I am familiar with the Contributing Guidelines.
• I agree with the terms of the Contributor Agreement.
• The changes introduced by this pull request are covered by existing or newly introduced tests.
• I have updated the changelog and the credits.

[codecov]

Codecov Report

Merging #171 (341cd53) into master (e1253af) will decrease coverage by 0.06%.
The diff coverage is 47.05%.

@@            Coverage Diff             @@
##           master     #171      +/-   ##
==========================================
- Coverage   54.83%   54.78%   -0.06%     
==========================================
  Files          26       26              
  Lines        2604     2621      +17     
==========================================
+ Hits         1428     1436       +8     
- Misses       1176     1185       +9     

Files Changed	Coverage Δ
coxeter/shapes/convex_polyhedron.py	48.21% <33.33%> (-0.85%)	⬇️
coxeter/shapes/polyhedron.py	63.83% <50.00%> (-0.64%)	⬇️

[vyasr]

Nice work! I have some suggestions for improvement here.

[vyasr]

This is an odd construct. You're using a list comprehension to mimic a for loop, but you're populating it with a bunch of None. Instead, you can do a nested list comprehension where you construct edges in the comprehension rather than appending to it, something like return [(i, j) for ... for face in self.faces.

[janbridley]

Based on your above feedback, I have come up with these options. The first option is an improved version of the append structure in the previous commit, which implements your suggestion of a nested list comprehension. It takes approximately 86us per iteration, as tested with timeit and 10 million iterations.

The second uses two list comprehensions, one to generate the edges and their reverses, and one to create the final list without reverse elements. It takes approximately 96 microseconds per iteration, tested in the same manner as above.

Note that both methods were formatted with black, so the images above are the formatted structure of the code as written.

Although I tried to come up with a cleaner algorithm for this, I could not make one that (1) returns the edges in an order corresponding to the input faces and (2) was faster than the appending method. I would be happy to hear your input on this if you have any, and if you like either of the above methods I will commit them.

[bdice]

Let's avoid using a list comprehension if the result is thrown out -- all those edges.append calls are returning a list of [None, None, None, ...] which is unused. Just use a for loop instead of a comprehension in this case, because this logic depends on the current value of the list (if (j, i) not in edges). It's a bit faster (no allocation for the throwaway result) and more readable.

Secondly, we can make this faster if we avoid np.roll and "roll our own." 😉 For educational purposes, I wrote up a microbenchmark linked below. Try it out locally! For me, "option 3" is around 3x faster than the other implementations. That said, I would not choose the final implementation based solely on speed. Readability is important.

https://github.com/bdice/python_microbenchmarks/blob/main/microbenchmarks/faces_to_edges.py

(Another small request -- please write formatted snippets with ```python code blocks instead of screenshots. It makes it easier to read and propose alternatives.)

[bdice]

@vyasr I think the reason a plain list comprehension doesn't work is mentioned above -- the edges are being de-duplicated in a way that depends on the current state of the edge list, and a comprehension can't reference its own state.

[vyasr]

Yup you're right I didn't account for the deduplication process needing to access the previously existing ones. A simple for loop is probably easiest alternative. I would definitely go for higher readability over reimplementing something like roll. coxeter is pretty suboptimal performance-wise anyway.

[janbridley]

Looks good - The roll1 function is impressive but I agree it is probably not the cleanest option.

(Another small request -- please write formatted snippets with ```python code blocks instead of screenshots. It makes it easier to read and propose alternatives.)

Noted - thank you for the information. The following code has been formatted as requested.

Option 1: most readable

def edges():
    return {
        (i, j) for face in self.faces for i, j in zip(face, np.roll(face, -1)) if i < j
    }

Option 2: fastest (by about an order of magnitude)

def edges():
  return{
      (face[i], face[(i + 1) % len(face)])
      for face in self.faces
      for i in range(len(face))
      if face[i] < face[(i + 1) % len(face)]
  }

As above, I think the more readable option 1 is the better choice for this package.

[janbridley]

Looks like CircleCI is failing because of inaccuracies in the dodecahedron's stored data. I feel that we should wait for #177 to be merged before merging this branch (#171), so that the small fixes currently in test_edges can be reverted before this branch goes live. This will result in a cleaner merge.

[bdice]

The docs should mention that edges are de-duplicated and that pairs are unordered. However, I'm not really sure if that's what we want. I would think that users would prefer to have control of the result, and de-duplicate and/or determine symmetry (i, j) == (j, i) for themselves. Should edges just be return [(i, j) for f in self.faces for i, j in zip(f, np.roll(f, -1))]?

[bdice]

I can see it either way. Both properties would be useful for different things. If you wanted to draw wireframes or make 3D models, you might need the full set of directed face edges (including "duplicates" belonging to different faces), but you can get that easily from the faces property. Perhaps de-duplicating and accounting for symmetry is the right move here.

[janbridley]

My personal feeling is that the 'cleanest' output should be accessible directly from the package, and the (slightly) more complicated case can be done by the user. That being said, I agree with your sentiment that neither method is incorrect. @vyasr I would appreciate your input here as well!

[bdice]

@janbridley Can you also give context on your motivation for this PR? What does your use case require?

[janbridley]

Also note that edges generated here are ordered: edges move clockwise around the surface of a face, and faces are ordered counterclockwise with respect to an outward normal. The is actually not the case, you are correct. The roll operation means that each face has at least one "jump" as edges are drawn. I will look at a method that better orders the edges before I update the docstring

Note that using j,i pairs reverses the direction of each vector, but keeps their clockwise ordering. I will update the docstring to clarify!

[bdice]

Sets of sets would be the most Python-friendly way to indicate this that I can think of. However, you'll run into a problem: sets are mutable, and thus not hashable. You'd have to have a set of frozenset objects (immutable sets, which are hashable). We could return a NumPy array, sorted and with only pairs where i < j.

[bdice]

I'd advise against multiple properties unless we have a really clear use case. I'm still not sure if I fully understand what the intended use case is for this property. How you do want to use the edges property?

[janbridley]

Edges allows for the calculation of each component of the asphericity of a polyhedra, allowing for a determination of the relative contribution of each feature. Edge_vectors allows for the calculation of angles and distances between specific edges on a polyhedron, which can be useful for simulations of highly anisotropic particles. @Tobias-Dwyer if you have any contributions I'd love to hear them!

[vyasr]

My personal feeling is that the 'cleanest' output should be accessible directly from the package, and the (slightly) more complicated case can be done by the user.

This is the right way to think about it. The package shouldn't be overloaded with a dozen variations on essentially the same property. From that perspective, given that as you've pointed out there are multiple possibilities for exactly what is required in different situations, would we be better served by including more documentation of how to compute edges in the documentation of the faces property? Otherwise we could easily proliferate multiple options based on considerations of ordering and uniqueness.

[janbridley]

Overall, I feel that some variant of edges should be included with the package - along with vertices and faces, edges are a defining element of polyhedra, and the number of either (i,j) or (j,i) edges allows for the correct calculation of the Euler characteristic, an intrinsic property of polyhedra in general. I agree that cluttering the package is undesirable, but I feel that including documentation on how to calculate edges is no cleaner than just efficiently implementing that method.

Looking at other widely used geometry packages, it seems like selecting edges with (i<j) is probably the best way forward for Coxeter. Mathematica's PolyhedronData[] is ordered this way, and it is probably the least confusing approach. CGAL provides the Halfedgetrait for polyhedra, which is equivalent to our (i,j) or (j,i) pairs, and allow users to select the direction and ordering of edges. However, CGAL's Polyhedron_3 does not calculate edges for an entire polyhedra at once: rather, there are independent member functions to find local half-edges from a vertex, edge, or face on the surface. This approach is probably too verbose for Coxeter, but I am not opposed to renaming our method halfedges if that would avoid confusion.

Based on this, I think the best approach is probably to return a set of edges with (i<j) and each edge included once. The documentation will be updated to better describe the method, and I have included two options for the new code above.

[bdice]

To test the edges property, I would consider some of its invariants / properties:

• assert that the first column of i indices is sorted
• assert that the the indices are ordered i < j
• assert that the array has an integer dtype
• etc.

[janbridley]

See 986fa8d

[DomFijan]

@bdice, any chance you can take a look at this before the end of Wednesday next week if you have the time? If not, I can jump in so we can finish this PR. There are a couple more tasks that depend on this that we want to start moving along. :)

[bdice]

A few small suggestions / verifications but otherwise LGTM. Feel free to merge once you are happy with it.

[bdice]

Suggested change

	* Added edges, edge_vectors, and num_edges methods
	* Added edges, edge_vectors, and num_edges methods.

[bdice]

Just making sure — the face indices are ordered consistently, right? In other words, this code relies on face indices being sorted in such a way that any shared edge is in the order (i, j) for one face and (j, i) for the other face sharing that edge. Is that a safe assumption from the rest of the code (it has been a while since I’ve read it). If this isn’t handled properly by the construction of the shape, then you’ll get missing edges. Are there any degenerate cases where this doesn’t happen automatically (can a polyhedron have one face in coxeter)?

[janbridley]

I am unsure exactly what the limits of Polyhedron are: for example, if self-intersecting shapes, shapes with holes, etc work properly. That being said, faces are ordered and edges are built using another method during initialization so it must work for the vast majority of use cases! I can take a deeper look in the next few days.

[janbridley]

Complex polyhedra (within limits) seem to work - self intersecting polyhedra, polyhedra with holes, and simple non convex (e.g. stellated) polyhedra all compute edges as expected. However, there are cases where volumes, moments of inertia, and other properties are incorrect. There are some cases where edges does NOT work - but in these cases, many other methods break down. Polygons (one-face polyhedron), dihedra, and nonclosed solids tend to break this method but they also return null values for any method that includes volume. Some also have erroneous behavior in the centroid or inertia tensor calculation. I will create an issue for some of these edge cases in general, but overall the new edge methods are at least as robust as the other core methods.

[bdice]

Good. It’s good to know the limitations — and which cases you are willing to treat as undefined behavior vs. raising an error or “fix” in some way.

[bdice]

This can be simpler if you want. I would just use pytest’s “assert raises” context manager, wrapping the assignment.

[DomFijan]

@bdice thanks for the prompt reaction :)