"look up" a face in the face lattice of a polyhedron

[gh-kliem]

We implement two methods that look up a face in the face lattice of a polyhedron:

• meet_of_Vrep -- the smallest face containing specified Vrepresentatives
• join_of_Hrep -- the largest face contained specified facets

This allows an easy answer for ​​https://ask.sagemath.org/question/34485/what-is-the-most-efficient-way-to-look-up-a-face-in-the-face-lattice-of-a-polyhedron/#50965

[gh-kliem]

Last 10 new commits:

​295039a	documentation
​d36da4a	coverage and small improvement
​2d0f0d9	method `reset` for the face iterator
​6d99a4c	typo
​5711f6b	method `ignore_subsets`
​f6633bd	method current and fix for reset
​e8c17c3	join_of_Vrep and meet_of_facets for face iterator
​b371a73	expose in combinatorial_polyhedron
​954b5c8	raise index error for index error
​521f9e0	expose in Polyhedron_base

[gh-kliem]

There are failing tests. I would wait until the dependices are taken care of to make things work again.

[gh-kliem]

New commits:

​19d2742	join_of_Vrep and meet_of_facets for face iterator
​c5c17ff	account for changes in the newest develop
​7a4b950	expose in combinatorial_polyhedron
​f108b06	raise index error for index error; fix doctests
​a69ba9f	expose in Polyhedron_base
​a6f53af	fix doctests

[gh-kliem]

I want to redesign some of the setup before making everything more complicated.

More precisely, creating strucutures face_struct and faces_list_struct that take care of the details. Then this overly long argument list of get_next_level will just reduce to three arguments or so. This would also make future changes easier including the transition to bitsets.pxi.

[gh-kliem]

New commits:

​2a14433	join_of_Vrep and meet_of_facets for face iterator
​530d85c	expose in combinatorial_polyhedron
​73a8be1	raise index error for index error; fix doctests
​c7ce411	expose in Polyhedron_base
​7e0a25e	fix doctests

[mkoeppe]

Some quick comments:

• "The offset is taken correctly" - what's that?
• Change "reseted" to "reset"
• "In case of an unbounded polyhedron" -> "In the case ..."

[gh-kliem]

I discovered a bug by accident:

sage: P = polytopes.permutahedron(3, backend='cdd')                                                                                                                                 
sage: [x.ambient_Hrepresentation() for x in P.facets()]                                                                                                                             
[(An inequality (1, 1, 0) x - 3 >= 0, An inequality (1, 0, 1) x - 3 >= 0),
 (An inequality (1, 1, 0) x - 3 >= 0, An inequality (1, 0, 0) x - 1 >= 0),
 (An inequality (1, 1, 0) x - 3 >= 0, An inequality (0, 1, 1) x - 3 >= 0),
 (An inequality (1, 1, 0) x - 3 >= 0, An inequality (0, 1, 0) x - 1 >= 0),
 (An inequality (1, 1, 0) x - 3 >= 0, An inequality (0, 0, 1) x - 1 >= 0),
 (An inequality (1, 1, 0) x - 3 >= 0, An equation (1, 1, 1) x - 6 == 0)]

The problem is that the backend cdd doesn't put equations in a stable place, which my code assumed. It depends whether you initialize from Vrep or from Hrep I guess.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​69740e7

improved documentation

[mkoeppe]

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​84d05ce

fix a variable name and add a doctest

[gh-kliem]

rebased

---

New commits:

​901b789	join_of_Vrep and meet_of_facets for face iterator
​2a5a782	expose in combinatorial_polyhedron
​6cecfd7	raise index error for index error; fix doctests
​9a66e32	expose in Polyhedron_base
​977c088	fix doctests
​15e2c58	improved documentation
​deca8e0	fix a variable name and add a doctest

[gh-kliem]

From ​https://trac.sagemath.org/ticket/30469#comment:25

• :meth:~sage.geometry.polyhedron.base.Polyhedron_base.meet_of_facets | largest face contained specified facets` --> contained in the specified facets
• Polyhedron_base.join_of_Vrep(self, *Vrepresentatives): ... in case of... --> "in the case of" * 2
• FaceIterator_base(SageObject).join_of_Vrep(self, *indices): .. NOTE:: ... may not te well defined. --> "be"
• Why is it meet_of_facets, but not meet_of_Hrep? If rays/lines are allowed on the primal side, why doesn't it make sense to consider equations on the dual side?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​ac9ff1e

typos

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​e9769ab

another typo

[gh-kliem]

Replying to gh-kliem:

From ​https://trac.sagemath.org/ticket/30469#comment:25

• :meth:~sage.geometry.polyhedron.base.Polyhedron_base.meet_of_facets | largest face contained specified facets` --> contained in the specified facets
• Polyhedron_base.join_of_Vrep(self, *Vrepresentatives): ... in case of... --> "in the case of" * 2
• FaceIterator_base(SageObject).join_of_Vrep(self, *indices): .. NOTE:: ... may not te well defined. --> "be"

All fixed.

• Why is it meet_of_facets, but not meet_of_Hrep? If rays/lines are allowed on the primal side, why doesn't it make sense to consider equations on the dual side?

Lines are basically ignored as are equations.

However, rays do make a difference. This is how I came up with it 12 months ago. join_of_Vrep is just better than join_of_vertices_and_rays, isn't it?

In the end, I don't know what is best.

[yzh]

I'm very new to the topic. I just started learning about combinatorial polyhedra today, and haven't read all the code changes on this ticket yet. But I'm wondering if it makes more sense to start by faces of the (possibly unbounded) polyhedron as input. Let

P = Polyhedron(vertices=[[1,0], [0,1]], rays=[[1,1]])
sage: V = P.Vrepresentation()
sage: V                                                                        
(A vertex at (0, 1), A vertex at (1, 0), A ray in the direction (1, 1))

Since "A ray in the direction (1, 1)" is not a face of P, V[2] shouldn't be taken as input for the look-up. Instead, the input could be something like [(1, 2), 0], where each element (= a single index or a tuple of indices) in the list must describe an actual face of the polyhedron.

For example, we can define the functions join_of_faces and meet_of_faces with input

• faces: list (iterator) of faces (not necessarily Vrep or facets), each face is described by a tuple of indices. Here index can be more flexible, such as P.Vrepresentation(4), etc.
• an argument telling whether the indices are primal or dual ones.

Just my two cents.

[gh-kliem]

As for the combinatorial polytope: It's just a coatom(ist)ic and atom(ist)ic lattice with some extra properties. Hence the name join_of and meet_of. As the lattice is coatom(ist)ic and atom(ist)ic the concept of meet_of_faces and join_of_faces is redundant. To compute the meet of faces it suffices to compute the meet of the coatoms/facets. As it is only one depth-first lookup the computation time of doing this is usually neglectable.

The unbounded polyhedron (no lines/linear subspaces assumed) is a bit different. It is basically a polytope with marked far face. The face lattice is given by [0, 1] \ [0, F]. Where [0,1] denotes the lattice of the bounded polytope and F is the far face.

The rays are vertices of that bounded polytope and thus it makes sense to include them in the notion of a join. Computing the join is always well defined with the exception for the case where you only include elements of the far face, e.g. rays.

However, I don't have a strong opinion on a name scheme. In the combinatorial polyhedron module it should be join and meet I think, but in polyhedron/base.py it might as well be look_up_face or find_face.

[yzh]

Given two faces (as polyhedra), if I'm interested in finding the largest face (as polyhedron) that is contained in both faces, is it better to call intersection or to apply join_of_facets?

[gh-kliem]

If you want the largest face contained in both of them, you need to do the intersection, which translates to the meet operation in the face lattice.

I would propose not defining the meet of the faces but instead the meet of the coatoms/facets. If you give me the faces as polyhedron, I need to first determine how your facet-description of the facets corresponds to the facets of the original polyhedron.

Meet should be much faster, because all the information is already encoded in the combinatorial polyhedron. It is just a simple depth first search in a lattice, doing a couple of intersections and subset checks of bitsets.

But it also depends on the kind of output you want.

[gh-kliem]

Actually what you describe is really easy to do, but neither implemented nor intended for in this ticket:

In your case you know (at least in the bounded case) that the intersection of both polytopes is exactly the convex hull of the intersection of the vertices. That is of course by far not true in the general case.

To determine the equations is really easy. It's just computing the kernel of a matrix. The possibly redundant inequalities are given by the union of the inequalities. So the only thing left to do (at least with backend field) is to figure out, which of the inequalities is redundant.

But this is exactly the functionality that face_as_combinatorial_polyhedron tries to implement as well.

[yzh]

Thanks!

Another possibly irrelevant question: for an empty face, does its ambient_H_indices or ambient_Hrepresentation mean anything?

P =  Polyhedron(rays=[(1,0),(0,1)])
F = P.faces(-1)[0]
F.ambient_H_indices()  #  returns  (0, 1)
F.ambient_Hrepresentation() #  (An inequality (1, 0) x + 0 >= 0, An inequality (0, 1) x + 0 >= 0)
P.faces(0)[0].ambient_Hrepresentation() # same output as above

[gh-kliem]

Stupid corner cases.

Yes, usually it ambient_H_indices or ambient_Hrepresentation does mean something for the empty face as well. But not in this case. There is at least two different definitions for faces:

• An intersection of hyperplanes corresponding to non-redundant inequalities.
• An intersection with a single hyperplane which corresponds to some inequality that the polyhedron satisfies (and maybe add the polyhedron as a face itself).

Usually they agree. But not in some case. In your example, there are two inequalities that are not redundant. They define those two unbounded 1-dimension faces, the 2-dimensional face is the empty intersection and the vertex. The empty face cannot be represented by those inequalities.

Different story for the second definition.

If you have at least two vertices, then things are fine. Sometimes cones are so much nicer...

That is why the CombinatorialPolyhedron itself just ignores the empty face and the full-dimensional face for the face iterator.

[yzh]

I'm reading base.py and find the following corner case:

sage: P = Polyhedron(vertices=[[1,0,0], [0,1,0]], rays=[[1,1,0]])               
sage: Q = Polyhedron(vertices=[[0,0],[0,1],[1,1],[1,0]])
sage: a, b, c = Q.Hrepresentation()[:3]
sage: P.meet_of_facets(b, c)
sage: P.meet_of_facets(a, b)

That is, join_of_Vrep and meet_of_facets probably need to check if v and facet given as PolyhedronFace or V/H-representation actually correspond to self.

Some typos:

• Line 6960: empty line but ::
• Line 7025: maybe you meant that the index cannot correspond to an equation in the Hrepresentation.
• Line 7030: ValueError: 0 is not a facet --> An equation (1, 1, 1, 1, 1) x - 15 == 0 is not a facet?

[yzh]

In /src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx, do the methods join_of_Vrep and meet_of_facets of CombinatorialPolyhedron have input indices corresponding to Vrepresentation/Hrepresentation or .vertices()/.facets()?

I guess it is the latter, but the documentation doesn't say. (And the name .._Vrepr somehow suggests the former?)

I also got a SIGSEGV when trying the following

sage: Q = Polyhedron(lines=[[1,1]])                                             
sage: CQ = CombinatorialPolyhedron(Q)
sage: CQ.join_of_Vrep(0)

[gh-kliem]

Replying to yzh:

In /src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx, do the methods join_of_Vrep and meet_of_facets of CombinatorialPolyhedron have input indices corresponding to Vrepresentation/Hrepresentation or .vertices()/.facets()?

Same thing.

I guess it is the latter, but the documentation doesn't say. (And the name .._Vrepr somehow suggests the former?)

I also got a SIGSEGV when trying the following

sage: Q = Polyhedron(lines=[[1,1]])                                             
sage: CQ = CombinatorialPolyhedron(Q)
sage: CQ.join_of_Vrep(0)

Bad.

[gh-kliem]

Replying to yzh:

I'm reading base.py and find the following corner case:

sage: P = Polyhedron(vertices=[[1,0,0], [0,1,0]], rays=[[1,1,0]])               
sage: Q = Polyhedron(vertices=[[0,0],[0,1],[1,1],[1,0]])
sage: a, b, c = Q.Hrepresentation()[:3]
sage: P.meet_of_facets(b, c)
sage: P.meet_of_facets(a, b)

That is, join_of_Vrep and meet_of_facets probably need to check if v and facet given as PolyhedronFace or V/H-representation actually correspond to self.

Some typos:

• Line 6960: empty line but ::

This is common practice in sage. I don't know if it is good practice. It just means that this is another example explained by the previous comment.

• Line 7025: maybe you meant that the index cannot correspond to an equation in the Hrepresentation.
• Line 7030: ValueError: 0 is not a facet --> An equation (1, 1, 1, 1, 1) x - 15 == 0 is not a facet?

Everything else will be fixed with the next commit.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​bec98df

improved documentation and check for elements to be of the correct polyhedron

[yzh]

Replying to gh-kliem:

Replying to yzh:

In /src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx, do the methods join_of_Vrep and meet_of_facets of CombinatorialPolyhedron have input indices corresponding to Vrepresentation/Hrepresentation or .vertices()/.facets()?

Same thing.

They are different in some cases, such as

sage: CQ.Hrepresentation()                                                      
(An equation (1, -1) x + 0 == 0,)
sage: CQ.facets()                                                               
()

[gh-kliem]

Ok. Forgot that.

[yzh]

This is certainly detail, but it bothers me that Polyhedron.join_of_Vrep can take a line as Vrep (it simply ignores the line) whereas Polyhedron.meet_of_facets would raise an error for equation or index corresponding to an equation. It is not symmetric. Also, I would give facet.dim() + 1 == self.dim() a separate error message.

[gh-kliem]

Replying to yzh:

Replying to gh-kliem:

Replying to yzh:

In /src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx, do the methods join_of_Vrep and meet_of_facets of CombinatorialPolyhedron have input indices corresponding to Vrepresentation/Hrepresentation or .vertices()/.facets()?

Same thing.

They are different in some cases, such as

sage: CQ.Hrepresentation()                                                      
(An equation (1, -1) x + 0 == 0,)
sage: CQ.facets()                                                               
()

The whole equation thing is so messed up, I don't even know where to start.

It was never a good decision to permit those in the first place.

Sometimes it's first the equations and then the inequalities, sometimes it's the other way around.

Maybe have same consistently last throughout the module? What do you think.

I have the feeling that there are some really bad design decisions that are haunting me.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​bfb08f5

typo

[gh-kliem]

I mean the ambient_H_indices of faces are not matching up the Hrepresentation. That is one thing that is annoying.

[gh-kliem]

Replying to yzh:

This is certainly detail, but it bothers me that Polyhedron.join_of_Vrep can take a line as Vrep (it simply ignores the line) whereas Polyhedron.meet_of_facets would raise an error for equation or index corresponding to an equation. It is not symmetric. Also, I would give facet.dim() + 1 == self.dim() a separate error message.

As stated above (at least I think I said this). It don't have a strong opinion on this. So you are proposing to just ignore equations? Or raise an error in both cases?

[yzh]

Indeed, the indexing is very confusing, even when there is no equations.

sage: P = Polyhedron(vertices=[[1,0], [0,1]], rays=[[1,1]], backend='field')    
sage: P.Hrepresentation()                                                       
(An inequality (1/2, -1/2) x + 1/2 >= 0,
 An inequality (-1/2, 1/2) x + 1/2 >= 0,
 An inequality (1/2, 1/2) x - 1/2 >= 0)
sage: P.facets()                                                                
(A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 2 vertices,
 A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray,
 A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray)
sage: CP = CombinatorialPolyhedron(P)                                           
sage: CP.Hrepresentation()                                                      
(An inequality (1/2, -1/2) x + 1/2 >= 0,
 An inequality (-1/2, 1/2) x + 1/2 >= 0,
 An inequality (1/2, 1/2) x - 1/2 >= 0)
sage: CP.facets()                                                               
((A vertex at (0, 1), A ray in the direction (1, 1)),
 (A vertex at (1, 0), A ray in the direction (1, 1)),
 (A vertex at (1, 0), A vertex at (0, 1)))

P.Hrepresentation() and P.facets() are not ordered in the same way. Neither do P.facets() and CP.facets(). (Please ignore the above code block. I think I understand it now.)

Let's focus on the CombinatorialPolyhedron first.

As you mentioned, it is weird that a face of a combinatorial polyhedron has ambient_H_indicies() not matching ambient_Hrepresentation() in length. Example with equation:

sage: P = Polyhedron(vertices=[[1,0,0], [0,1,0]], rays=[[1,1,0]], backend='field')                                                                ​
sage: P.Hrepresentation()                                                       
(An inequality (1/2, -1/2, 0) x + 1/2 >= 0,
​An inequality (-1/2, 1/2, 0) x + 1/2 >= 0,
​An inequality (1/2, 1/2, 0) x - 1/2 >= 0,
​An equation (0, 0, 1) x + 0 == 0)
sage: P.facets()                                                                
(A 1-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 2 vertices,
​A 1-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 1 vertex and 1 ray,
​A 1-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 1 vertex and 1 ray)
sage: CP = CombinatorialPolyhedron(P)                                           
sage: CP.Hrepresentation()                                                      
(An equation (0, 0, 1) x + 0 == 0,
​An inequality (1/2, -1/2, 0) x + 1/2 >= 0,
​An inequality (-1/2, 1/2, 0) x + 1/2 >= 0,
​An inequality (1/2, 1/2, 0) x - 1/2 >= 0)
sage: CP.facets()                                                               
((A vertex at (0, 1, 0), A ray in the direction (1, 1, 0)),
​(A vertex at (1, 0, 0), A ray in the direction (1, 1, 0)),
​(A vertex at (1, 0, 0), A vertex at (0, 1, 0)))
sage: list(P.face_generator())                                                  
[A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 2 vertices and 1 ray,
​A -1-dimensional face of a Polyhedron in QQ^3,
​A 1-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 2 vertices,
​A 1-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 1 vertex and 1 ray,
​A 1-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 1 vertex and 1 ray,
​A 0-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 1 vertex,
​A 0-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 1 vertex]
sage: list(CP.face_iter())                                                      
[A 1-dimensional face of a 2-dimensional combinatorial polyhedron,
​A 1-dimensional face of a 2-dimensional combinatorial polyhedron,
​A 1-dimensional face of a 2-dimensional combinatorial polyhedron,
​A 0-dimensional face of a 2-dimensional combinatorial polyhedron,
​A 0-dimensional face of a 2-dimensional combinatorial polyhedron]
sage: f = list(P.face_generator())[2]
sage: f.ambient_Hrepresentation()                                       
(An equation (0, 0, 1) x + 0 == 0, An inequality (1/2, 1/2, 0) x - 1/2 >= 0)
sage: f.ambient_H_indices()  # not sorted                                         
(3, 2)
sage: cf = next(CP.face_iter())
sage: cf.ambient_Hrepresentation()                                      
(An inequality (1/2, 1/2, 0) x - 1/2 >= 0, An equation (0, 0, 1) x + 0 == 0)
sage: cf.n_ambient_Hrepresentation()     # discard equation                          
1
sage: cf.ambient_H_indices()             # discard equation
(2,)

In my opinion, even through the lines and the equations are not used in the computation of a Combinatorial Polyhedron, they should be included in the ambient_V/Hrepresentaion and ambient_V/H_indices. It would be hard to choose between equations+inequalities or inequalities+equations for Polyhedron_base since the order is not consistent throughout different backends. However, for CombinatorialPolyhedron, I think it would make sense to fix the order, say, to inequalities+equations. (Note that currently it is equations+inequalities, which seems bad to me.)

Using this order with inequalities first, there will be no strange offset treatment needed for meet_of_facets. We just require that the indices of CP.meet_of_facets(*indices) cannot exceed CP.n_facets(). Raise error otherwise.

Similarly, in CP.Vrepresentation(), the order is vertices+rays+lines (as in the code right now). Thus, for CP.join_of_Vrep(*indices), we just require that indices cannot exceed CP.n_vertices() + CP.n_rays(). Raise error otherwise.

Perhaps, it would be useful to have CP.rays/lines/inequalities/equations and cf.vertices/rays/lines/inequalities/equations just for output? Not sure about this, but at least CP.n_rays() seems needed.

Expected returns

sage: f.ambient_Hrepresentation()                                       
(An inequality (1/2, 1/2, 0) x - 1/2 >= 0, An equation (0, 0, 1) x + 0 == 0)
sage: f.ambient_H_indices()                                    
(2, 3)
sage: cf.ambient_Hrepresentation()                                      
(An inequality (1/2, 1/2, 0) x - 1/2 >= 0, An equation (0, 0, 1) x + 0 == 0)
sage: cf.n_ambient_Hrepresentation()
sage: cf.ambient_H_indices()
(2, 3)

[yzh]

Now consider Polyhedron_base.

To be honest, I often forget what join_of_.. and meet_of_... mean in this class. I propose to rename them. How do def least_common_superface_of_Vrep(self, *Vrepresentatives) and def greatest_common_subface_of_Hrep(self, *Hrepresentatives) sound? I would prefer that the two functions are sort of symmetric. Lines and equations are allowed, but are simply ignored when passing to self.face_generator().

In the Hrepresentation, ppl backend puts equations first, and field, cdd, normaliz put inequalities first. (I was not able to install polymake backend.) Thus, the offset treatment seems needed.

[gh-kliem]

Note that equations would be ignored explicitly and lines implicitly. But for the user it's the same, I guess.

As mentioned before, CombinatorialPolyhedron has some difficult design decisions. Equations are almost everywhere ignored. They are just added in the end to every face etc. Lines are not everywhere ignored. They really belong to the Vrepresentation of the CombinatorialPolyhedron and are part of the indices.

A reason for this is that I wanted to have the dimensions of the faces match up with the original polyhedron. So lines are needed to push the dimension up. Is this a bad choice? I don't know exactly, however I would like to wait with this until someone really wants to work with polyhedra with lines.

One problem is the following:

Maybe we could treat lines and equations more symmetric at some point. I don't know what would be best.

I would prefer to keep join_of_ etc. However, it does not hurt to add the aliases you suggested. They do make sense. Maybe meet_of_Hrep is really better for Polyhedron_base.

CombinatorialPolyhedron.facets and Polyhedron_base.facets are really different things unfortunatly. Polyhedron_base.facets just gives you the facets as faces in the order given by the face iterator (reverse). CombinatorialPolyhedron.facets gives you the facets as Vrepresentation indices in the order of Hrepresentation. Maybe I should just document this.

The order of the Vrepresentation of CombinatorialPolyhedron is not fixed. This just happened because you use Polyhedron_base to initialize it. You can give it in any order you like and shuffle. Hence the strange method vertices that gives you all the elements of Vrepresentation that are not contained in the far face. (This is also historic. At some point I wanted to initialize CombinatorialPolyhedron just from the incidence matrix, which is not possible, because some unbounded polyhedra have the same incidence matrix, if you do not consider the far face).

CombinatorialPolyhedron does not distinguish rays and lines currently. Those are just things that happen to be in the Vrepresentation. They are kept as indices and both contained in the far face. They also cannot really be distinguished in sage, because Cone does not know lines at all. It just keeps to opposite rays:

sage: C = Cone([[1,0], [0,1], [-1,0]])   
sage: C.incidence_matrix()                                                                                                                                                          
[0]
[1]
[1]
sage: CP = CombinatorialPolyhedron(C)                                                                                                                                               
sage: CP.Vrepresentation()                                                                                                                                                          
(N(0, 1), N(1, 0), N(-1, 0), N(0, 0))

What's the number of rays in this case. Is it one or three?

[gh-kliem]

I created #31834.

Please have a look after I pushed the changes, whether this looks reasonable.

[yzh]

Replying to gh-kliem:

Equations are almost everywhere ignored. They are just added in the end to every face etc. Lines are not everywhere ignored. They really belong to the Vrepresentation of the CombinatorialPolyhedron and are part of the indices.

That makes sense to me. I often need to deal with polyhedra with equations, so #31834 would be very helpful to me.

I would prefer to keep join_of_ etc. However, it does not hurt to add the aliases you suggested. They do make sense. Maybe meet_of_Hrep is really better for Polyhedron_base.

Sounds good! Keep join_of_Vrep and meet_of_Hrep for the base class and make aliases.

CombinatorialPolyhedron.facets and Polyhedron_base.facets are really different things unfortunatly.

Yes, I got it after looking into the code closely. At first, I imagined that CombinatorialPolyhedron.facets returns CombinatorialFace objects. (The output format was described in the docstring but I was not careful in reading it.) What's the reason to use an unusual output format? Can the conversion to Vrep indices be deferred to where that's needed (such as in __reduce__)?

The order of the Vrepresentation of CombinatorialPolyhedron is not fixed.

I assumed that CombinatorialPolyhedron has Vrepresentation being "[vertices, rays, lines]" because this expression appears many times in combinatorial_polyhedron/base.pyx Does it make sense to fix the order, or is it a waste of time? In fact, .../polyhedron/representation.py defines INEQUALITY = 0, EQUATION = 1, VERTEX = 2, RAY = 3, LINE = 4. It's also used in #31702. Currently, ppl and cdd backends do not respect this order (creating a lot of trouble for my doctests), and the other backends field and normaliz respect this order.

CombinatorialPolyhedron does not distinguish rays and lines currently. Those are just things that happen to be in the Vrepresentation. They are kept as indices and both contained in the far face.

Ok. But don't we want a "canonical" representation (in terms of numbers at least) for CombinatorialPolyhedron? Note that the polyhedron base class has that. For instant, are the combinatorial polyhedron generated by two opposite rays and the other polyhedron generated by one line considered the same?

They also cannot really be distinguished in sage, because Cone does not know lines at all. It just keeps to opposite rays: ... What's the number of rays in this case. Is it one or three?

Although Cone allows for non-pointed cones, I believe that Cone and in particular the Fan class mainly deal with pointed (i.e. strictly convex) cones.

    When the cone in question is not strictly convex, the standard form for
    the "generating rays" of the linear subspace is "basis vectors and their
    negatives", as in the following example::

        sage: plane = Cone([(1,0), (0,1), (-1,-1)])
        sage: plane.rays()
        N( 0,  1),
        N( 0, -1),
        N( 1,  0),
        N(-1,  0)
        in 2-d lattice N

Because Cone has its rays in a "canonical" form, it would be possible for a CombinatorialPolyhedron constructed from a Cone to combine opposite rays to lines.

In my opinion, CP in your example should have one vertex one ray and one line. I don't know how to express the line in terms of N(?, ?), but writing the vertex as N(0, 0) looks weird to me too.

[gh-kliem]

Just as with cones, Polyhedron_base mainly deals with polytopes and strictly convex polyhedra.

Things with lines should be correct, but I don't think it is a very important use case. So if things are a bit ugly in the backend, it is fine unless someone really works with that stuff and wants high speed performance (or unless working around the bad design is more work then creating a better design and working with it).

[gh-kliem]

While I agree that #31834 should be a base for this ticket, I don't think we should make this ticket here dependent on the other fixes.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

​ddce44a	raise index error for index error; fix doctests
​1305ea2	expose in Polyhedron_base
​9322749	fix doctests
​11e4a50	improved documentation
​fa75eb0	fix a variable name and add a doctest
​2f66609	typos
​44fbabf	another typo
​f126c36	improved documentation and check for elements to be of the correct polyhedron
​87b2f4d	typo
​c1c14b5	ignore equations; make aliases

[gh-kliem]

Ok, I think I have your requests taken care of for now.

[yzh]

src/sage/geometry/polyhedron/base.py

• facets -> Hrepresentatives

  def meet_of_Hrep(self, *Hrepresentatives):
      r"""
      Return the largest face that is contained in ``facets``.
  

• Equations at the end. Need to fix #31834. Without the fix, there could be a bug (offset=0 but ambient_H_indices returns equations first) at facet = H_indices[0] if H_indices[0] >= offset else H_indices[-1]. After the fix, the expected return of the doctest is (An inequality (0, 0, 1) x - 1 >= 0, An equation (1, 1, 1) x - 6 == 0).

  def meet_of_Hrep(self, *Hrepresentatives):
      TESTS:
        sage: P = polytopes.permutahedron(3, backend='field'
        sage: P.meet_of_Hrep(0).ambient_Hrepresentation()
        (An equation (1, 1, 1) x - 6 == 0, An inequality (0, 0, 1) x - 1 >= 0)
  

• Error message raise ValueError("{} is not a Hrepresentative".format(facet))

  else:
      raise ValueError("{} is the index of an equation and not of an inequality".format(facet))
  

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

​909c1fc	expose in Polyhedron_base
​d6fefe4	fix doctests
​50d7252	improved documentation
​cb6256a	fix a variable name and add a doctest
​c838617	typos
​fd044af	another typo
​20e282e	improved documentation and check for elements to be of the correct polyhedron
​f3187b5	typo
​e0841a9	ignore equations; make aliases
​be6236f	small fixes

[gh-kliem]

Replying to yzh:

src/sage/geometry/polyhedron/base.py

• facets -> Hrepresentatives

  def meet_of_Hrep(self, *Hrepresentatives):
      r"""
      Return the largest face that is contained in ``facets``.
  

• Equations at the end. Need to fix #31834. Without the fix, there could be a bug (offset=0 but ambient_H_indices returns equations first) at facet = H_indices[0] if H_indices[0] >= offset else H_indices[-1]. After the fix, the expected return of the doctest is (An inequality (0, 0, 1) x - 1 >= 0, An equation (1, 1, 1) x - 6 == 0).

  def meet_of_Hrep(self, *Hrepresentatives):
      TESTS:
        sage: P = polytopes.permutahedron(3, backend='field'
        sage: P.meet_of_Hrep(0).ambient_Hrepresentation()
        (An equation (1, 1, 1) x - 6 == 0, An inequality (0, 0, 1) x - 1 >= 0)
  

• Error message raise ValueError("{} is not a Hrepresentative".format(facet))

  else:
      raise ValueError("{} is the index of an equation and not of an inequality".format(facet))
  

Thanks. Should be taken care of.

[yzh]

The segmentation fault persists

sage: P = Polyhedron(lines=[[1,1]])                                                                        
sage: C = CombinatorialPolyhedron(P)                                                                       
sage: C.Vrepresentation()                                                                                  
(A line in the direction (1, 1), A vertex at (0, 0))
sage: C.join_of_Vrep(0)   

[yzh]

sage: P.join_of_Vrep(0) 

also causes a segmentation fault.

[gh-kliem]

And meet_of_Hrep as well.

I'm currently trying to figure out how to resolve this. The FaceIterator was never initialized in that case, which makes life a bit annoying.

[yzh]

Replying to gh-kliem:

And meet_of_Hrep as well.

I'm currently trying to figure out how to resolve this. The FaceIterator was never initialized in that case, which makes life a bit annoying.

I guess it's only on the extreme corner case where both C.vertices() and C.facets() are empty.

[gh-kliem]

I'm actually surprised that Polyhedron([[0,0]]) works well. I did not expect this.

[yzh]

Replying to gh-kliem:

I'm actually surprised that Polyhedron([[0,0]]) works well. I did not expect this.

Polyhedron([[0,0]]) is not as special as Polyhedron(lines=[[1,1]]) because the former has a vertex in its combinatorial polyhedron. To fix the segmentation fault, can we just raise error in FaceIterator_base.meet_of_Hrep etc when i >= self._n_facets + self._n_equations?

[gh-kliem]

The special thing about it is that it has a far face. It tries to figure out, whether face is contained in the far face, before finding the face via the iterator.

It assumes that the first face of visited_all is the far face, which is not initalized, hence the segmentation fault.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3d74488

fix segementation fault

[yzh]

A weird question: why are the error messages different? There is nothing wrong, though

+            sage: it._meet_of_coatoms(-1)
+            Traceback (most recent call last):
+            ...
+            OverflowError: can't convert negative value to size_t
...
+            sage: it._join_of_atoms(-1)
+            Traceback (most recent call last):
+            ...
+            IndexError: atoms out of range

[yzh]

Another question that might sound silly: Why does if n_atoms == self.coatoms.n_atoms(): mean that the face is the universe?

[yzh]

I have trouble understanding the following code.

+        while self.structure.current_dimension != self.structure.dimension:
+            if face_issubset(face, self.structure.face):
+                if face_issubset(self.structure.face, face):
+                    # Found our face.
+                    return 0
+            else:
+                # The face is not a subface/supface of the current face.
+                self.ignore_subsets()
+
+            d = self.next_dimension()

When face_issubset(face, self.structure.face) is False, after self.ignore_subsets(), why does it move to the next dimension directly? Does it skip all other faces of the current dimension?

Sorry for asking so many questions.

[gh-kliem]

Replying to yzh:

Another question that might sound silly: Why does if n_atoms == self.coatoms.n_atoms(): mean that the face is the universe?

If you contain all the atoms there are, then you must be the universe. A face containing all of the Vrepresentation is the universe.

[gh-kliem]

Replying to yzh:

I have trouble understanding the following code.

+        while self.structure.current_dimension != self.structure.dimension:
+            if face_issubset(face, self.structure.face):
+                if face_issubset(self.structure.face, face):
+                    # Found our face.
+                    return 0
+            else:
+                # The face is not a subface/supface of the current face.
+                self.ignore_subsets()
+
+            d = self.next_dimension()

When face_issubset(face, self.structure.face) is False, after self.ignore_subsets(), why does it move to the next dimension directly? Does it skip all other faces of the current dimension?

Sorry for asking so many questions.

Maybe self.next_dimension() is not the best name. It goes to the next face and returns it's dimension. Might not change.

[gh-kliem]

Replying to yzh:

A weird question: why are the error messages different? There is nothing wrong, though

+            sage: it._meet_of_coatoms(-1)
+            Traceback (most recent call last):
+            ...
+            OverflowError: can't convert negative value to size_t
...
+            sage: it._join_of_atoms(-1)
+            Traceback (most recent call last):
+            ...
+            IndexError: atoms out of range

Because it tries to convert -1 to a non-negative integer (I think when assigning it to i, which is defined to be of size_t). It fails then already and not in the next step.

[yzh]

Can def _meet_of_coatoms use

if not all(i in range(n_coatoms) for i in indices):
           raise IndexError("coatoms out of range")

at the beginning to make the error message more consistent?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​c3f028b

check with a nice error for negative indices

[gh-kliem]

ok, changed

[yzh]

It's true that face is contained in the far face, but why does this cause an error?

+        elif not self._bounded and not self.coatoms.n_faces():
+            # Note: It is important to catch this and not to run ``face_issubset`` to prevent a segmentation fault.
+            if face_len_atoms(face):
+                # Contained in the far face, if it contains anything.
+                raise ValueError("the join is not well-defined")

Let P = Polyhedron(lines=[[1, 1]]), I expected to get the line itself for P.join_of_Vrep(0), P.join_of_Vrep(1) and P.join_of_Vrep(0, 1). I thought that the smallest face containing (A line in the direction (1, 1) and/or A vertex at (0, 0)) is well defined in this example.

One can compare the above example to the following, which runs as expected.

sage: P = Polyhedron(lines=[[1, 1]], rays=[(-1,1)])                                                        
sage: P.Vrepresentation()                                                                                  
(A line in the direction (1, 1),
 A ray in the direction (-1, 0),
 A vertex at (0, 0))
sage: P.join_of_Vrep(2)                                                                                    
A 1-dimensional face of a Polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 1 line
sage: P.join_of_Vrep(0)                                                                                    
A 1-dimensional face of a Polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 1 line
sage: P.join_of_Vrep(0, 2)                                                                                 
A 1-dimensional face of a Polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 1 line
sage: P.join_of_Vrep(1, 2)                                                                                 
A 2-dimensional face of a Polyhedron in ZZ^2 defined as the convex hull of 1 vertex, 1 ray, 1 line

[gh-kliem]

Replying to yzh:

It's true that face is contained in the far face, but why does this cause an error?

+        elif not self._bounded and not self.coatoms.n_faces():
+            # Note: It is important to catch this and not to run ``face_issubset`` to prevent a segmentation fault.
+            if face_len_atoms(face):
+                # Contained in the far face, if it contains anything.
+                raise ValueError("the join is not well-defined")

Let P = Polyhedron(lines=[[1, 1]]), I expected to get the line itself for P.join_of_Vrep(0), P.join_of_Vrep(1) and P.join_of_Vrep(0, 1). I thought that the smallest face containing (A line in the direction (1, 1) and/or A vertex at (0, 0)) is well defined in this example.

One can compare the above example to the following, which runs as expected.

sage: P = Polyhedron(lines=[[1, 1]], rays=[(-1,1)])                                                        
sage: P.Vrepresentation()                                                                                  
(A line in the direction (1, 1),
 A ray in the direction (-1, 0),
 A vertex at (0, 0))
sage: P.join_of_Vrep(2)                                                                                    
A 1-dimensional face of a Polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 1 line
sage: P.join_of_Vrep(0)                                                                                    
A 1-dimensional face of a Polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 1 line
sage: P.join_of_Vrep(0, 2)                                                                                 
A 1-dimensional face of a Polyhedron in ZZ^2 defined as the convex hull of 1 vertex and 1 line
sage: P.join_of_Vrep(1, 2)                                                                                 
A 2-dimensional face of a Polyhedron in ZZ^2 defined as the convex hull of 1 vertex, 1 ray, 1 line

Just copying self._far_face properly takes care of this and the segmentation fault as well.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​110228e

better check for the far face containment

[yzh]

It's still not all correct in the corner cases.

sage: P = Polyhedron(vertices=[[1,0]], lines=[[1,1]])                                                      
sage: P.Vrepresentation()                                                                                  
(A line in the direction (1, 1), A vertex at (1, 0))
sage: P.join_of_Vrep(0)                                                                                    
A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 line

and

sage: R = Polyhedron(vertices=[[1,0]], rays=[[1,1]])                                                       
sage: R.Vrepresentation()                                                                                  
(A vertex at (1, 0), A ray in the direction (1, 1))
sage: R.join_of_Vrep(1)                                                                                    
A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray

Expected to have ValueError: the join is not well-defined for both.

With one more vertex, it's correct though:

sage: Q = Polyhedron(vertices=[[1,0], [0,1]], lines=[[1,1]])                                               
sage: Q.Vrepresentation()                                                                                  
(A line in the direction (1, 1), A vertex at (-1, 0), A vertex at (1, 0))
sage: Q.join_of_Vrep(0)
ValueError: the join is not well-defined

[gh-kliem]

Replying to yzh:

It's still not all correct in the corner cases.

sage: P = Polyhedron(vertices=[[1,0]], lines=[[1,1]])                                                      
sage: P.Vrepresentation()                                                                                  
(A line in the direction (1, 1), A vertex at (1, 0))
sage: P.join_of_Vrep(0)                                                                                    
A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 line

and

sage: R = Polyhedron(vertices=[[1,0]], rays=[[1,1]])                                                       
sage: R.Vrepresentation()                                                                                  
(A vertex at (1, 0), A ray in the direction (1, 1))
sage: R.join_of_Vrep(1)                                                                                    
A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray

Expected to have ValueError: the join is not well-defined for both.

With one more vertex, it's correct though:

sage: Q = Polyhedron(vertices=[[1,0], [0,1]], lines=[[1,1]])                                               
sage: Q.Vrepresentation()                                                                                  
(A line in the direction (1, 1), A vertex at (-1, 0), A vertex at (1, 0))
sage: Q.join_of_Vrep(0)
ValueError: the join is not well-defined

It's consistent. Still one can argue about whether it should be changed.

What you are proposing is the following: Given a bunch of indices. If there are in the far face, throw an error.

Current behavior is the following (and I also wasn't aware of this, after all the ticket has been sitting there for a year): Given no indices, return the empty face. Otherwise, intersect all facets that contain all given indices. If that is contained in the far face, throw an error.

In the first of your cases, there is no facet. In the second case, the only facet does not contain the index. In the third case, both facets contain the line and thus the join is really not defined.

A different example where currently everything works out is the following:

sage: P = Polyhedron(vertices=[[0,1],[1,0]], rays=[[0,1],[1,0]])                                                                                                                    
sage: P.Vrepresentation()                                                                                                                                                           
(A ray in the direction (0, 1),
 A vertex at (0, 1),
 A ray in the direction (1, 0),
 A vertex at (1, 0))
sage: for i in range(4): 
....:     P.join_of_Vrep(i) 
....:                                                                                                                                                                               
A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray
A 0-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex
A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray
A 0-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex

There is always exactly one face containing each ray. It even makes sense when you consider the lattice of the cone (or add the far face). The ray is of course a face of it's own there. However, there is a unique smallest face in the original polyhedron that contains it.

I think the current behavior makes sense. You could give me a bunch of rays of a large polyhedron and ask me, what's the smallest face containing them. As long as this is well-defined, you can except to get an answer and not an error. It's like computing gcd and lcm. It tries to find an answer for you and only fails if those two numbers do not have a unique meet or join.

[yzh]

Actually, I'm not considering combinatorial polyhedra, far face, etc. at all. Instead, I'm looking from a geometric point of view.

In the examples above, I see "A line in the direction (1, 1)" as l = Polyhedron(lines=[(1,1)]), and "A ray in the direction (1, 1)" as r = Polyhedron(rays=[(1,1)]). Recall that for a geometric polyhedron (of Polyhedron_base class), it is a face of itself, and the empty polyhedron is also a face.

None of the faces of P (including itself) contains l. Similarly, none of the faces of Q (including itself) contains r. Thus, I expect an error in the first two examples.

None of the faces of Q contains l. Therefore, the third example returns an error.

[gh-kliem]

I don't think one should see "A line in the direction (1, 1)" as l = Polyhedron(lines=[(1,1)]), because whenever you don't specify vertices, but rays or lines, then the vertex (0,0) is added.

P contains the line and not l, because l is the line and the vertex (0,0).

No, the third example returns an error, because both edges of Q contain the line, but not a vertex. So it is not well-defined.

Edit: Maybe rather say 1-faces instead of edges. There are two 1-faces of Q that contain the lines, but no 0-face of Q contains the line.

[yzh]

Replying to gh-kliem:

I don't think one should see "A line in the direction (1, 1)" as l = Polyhedron(lines=[(1,1)])

Aha! That explains everything.

[gh-kliem]

Thank you.

That part with the implicit vertex is super confusing...

[gh-kliem]

Last 10 new commits:

​6e4fb22	fix a variable name and add a doctest
​45cffb0	typos
​117a15e	another typo
​5d2d329	improved documentation and check for elements to be of the correct polyhedron
​6ee3093	typo
​4cf2b76	ignore equations; make aliases
​40cb506	small fixes
​dae644e	fix segementation fault
​fd50a98	check with a nice error for negative indices
​644a803	better check for the far face containment

[gh-kliem]

Rebased on new version of #31834, which fixes trivial merge conflict because of whitespaces.

[gh-kliem]

Only merged and rebased on #31834.

---

Last 10 new commits:

​9cbd015	fix a variable name and add a doctest
​4612fe5	typos
​3b09b68	another typo
​fe428aa	improved documentation and check for elements to be of the correct polyhedron
​dd79e1c	typo
​56c66a4	ignore equations; make aliases
​80c80e5	small fixes
​9afb276	fix segementation fault
​ef5c95e	check with a nice error for negative indices
​99ddc2f	better check for the far face containment