Schlegel projection breaks convexity

[jipilab]

The documentation string of .schlegel_projection reads:

   Return the Schlegel projection.

   * The polyhedron is translated such that its "center()" is at the
     origin.

   * The vertices are then normalized to the unit sphere

   * The normalized points are stereographically projected from a
     point slightly outside of the sphere.

When normalizing to the unit sphere this (potentially) completely breaks the convexity of the object.

Minimal example:

sage: fcube = polytopes.hypercube(4)
sage: tfcube = fcube.face_truncation(fcube.faces(0)[0])
sage: sp = tfcube.schlegel_projection()
sage: sp.plot()
Launched html viewer for Graphics3d Object

The pentagons are not planar although they should in a schlegel diagram.

The scaling to the unit sphere should be removed to preserve convexity.

This ticket fixes the projection while de-duplicating some code.

FOLLOW-UP: Fix .plot() and .show() to have a better behaviour of smaller dimensional objects.

[gh-kliem]

Is this the bug you showed during your presentation at sd109?

[jipilab]

New commits:

​26c03d1	Fixed chirality in stereographic projection
​bed13a6	examples and pep8

[gh-kliem]

Replying to gh-kliem:

Is this the bug you showed during your presentation at sd109?

?

[jipilab]

Replying to gh-kliem:

Replying to gh-kliem:

Is this the bug you showed during your presentation at sd109?

?

No, that one had to do with plotting in 3d and having missing faces or so. This one is really different: it was providing a wrong output.

Further, the method is improved and made more user-friendly, one just needs to be fed a facet and a positive number to get a schlegel diagram. Close to 0, the projection point is very close to the barycenter of the facet, otherwise, getting away in the direction of the representative point of the locus of points that see that facet.

Given these 2 things, this provides a projection of the polyhedron (think typically of something in 4d) into the affine hull of the facet, and that image is then transformed orthonormally into 3d and then the picture is drawn.

[jipilab]

Ach...

sage: cp = polytopes.cyclic_polytope(4,8).polar()                                                                                                                                                                  
sage: cp.schlegel_projection()                                                                                                                                                                                     
Traceback (most recent call last)
[snip]

~/sage/local/lib/python3.8/site-packages/sage/geometry/polyhedron/plot.py in __call__(self, x)
    302         segment = Polyhedron(vertices=[vector(RDF, x),self.projection_point])
    303         # The intersection of the segment with the facet
--> 304         preimage = (segment & self.facet).vertices()[0].vector()
    305         # The transformation matrix acts on the right:
    306         return preimage*self.A + self.b

IndexError: tuple index out of range

That line is an overkill... I should probably change it to something a bit simpler than taking the intersection.

The problem occurs prior to this, this is an artefact from the older version, somehow it is already converted to RDF before and then no intersection is found.

This example is likely to be one the extreme side, but it is good to fix this now.

[gh-LaisRast]

Replying to jipilab:

Ach...

sage: cp = polytopes.cyclic_polytope(4,8).polar()                                                                                                                                                                  
sage: cp.schlegel_projection()                                                                                                                                                                                     
Traceback (most recent call last)
[snip]

~/sage/local/lib/python3.8/site-packages/sage/geometry/polyhedron/plot.py in __call__(self, x)
    302         segment = Polyhedron(vertices=[vector(RDF, x),self.projection_point])
    303         # The intersection of the segment with the facet
--> 304         preimage = (segment & self.facet).vertices()[0].vector()
    305         # The transformation matrix acts on the right:
    306         return preimage*self.A + self.b

IndexError: tuple index out of range

That line is an overkill... I should probably change it to something a bit simpler than taking the intersection.

In Ziegler's book (Lectures on Polytopes) page 133, there is a nice formula for computing the intersection point.

[jipilab]

Replying to gh-LaisRast:

In Ziegler's book (Lectures on Polytopes) page 133, there is a nice formula for computing the intersection point.

Thanks for the pointer, that will do. There is still the issue that at that moment, one deals already with RDF, but having a direct formula should improve it. (There is still the danger of dividing by 0...)

[git]

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

​2b6e5c5

Changed method to get projection

[jipilab]

I still get an error at line 917 in base.py, there is a NaN... which I believe is coming from this division. Hmm. I'll see what can be done.

[git]

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

​4342b52

Fix doctests

[jipilab]

That should do it.

[gh-kliem]

• # There is no 4-d screen, we must project down to 3d

Either 3-d or 4d I would say.

• +        if not locus_polyhedron.relative_interior_contains(projection_point):
  +            raise ValueError("the chosen position is too large")
  

  If you enter position=-1, this is the error you see, which is somewhat strange. (Of course bullshit gives bullshit, but maybe we can be nice here.)

[gh-kliem]

There is a trailing space in the preimage line.

[gh-kliem]

• The biggest trouble with schlegel is that the default position now seems to be much too close !(?)

Take a look at at

polytopes.six_hundred_cell().plot() etc.

The default position that works good for stacking, seems to be too close for schlegel projection.

• Even worse: polytopes.permutahedron(4).show() isn't orthonormally projected anymore.

[git]

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

​c3df1d5

Add some Errors and easy fixes

[jipilab]

Replying to gh-kliem:

• The biggest trouble with schlegel is that the default position now seems to be much too close !(?)

Take a look at at

polytopes.six_hundred_cell().plot() etc.

The default position that works good for stacking, seems to be too close for schlegel projection.

Defaulting behavior will never be perfect.

• Even worse: polytopes.permutahedron(4).show() isn't orthonormally projected anymore.

It was a mere coincidence that it was orthonormally projected by the method show (!). This indicates that if one has a lower dimensional object (<=3) in an ambient dimension which is strictly larger than 3, one would prefer to visualize (by default) the object in its affine hull, and that, orthonormally projected. This is what you mean?

In this case, this is a different ticket fixing the .show()/.plot() methods. This ticket fixes the schlegel projection which it does...

[gh-kliem]

Thank you for exposing the position argument.

Of course default behavior will never be perfect, but it is stupid to be stuck with an awful default and don't know how to change it.

If I'm now unhappy with the default, I can inspect the plot method to find the option I need.

polytopes.permutehdron(4).show() was orthonormally projected because the default for ambient dimension 4 is (previous to this ticket) schlegel projection. So it was orthonormal for the wrong reasons. Nevertheless, it worked.

How about adding

src/sage/geometry/polyhedron/base.py

@@ -1040 +1040 @@
             elif polyhedron.dimension() == 4:
                 # There is no 4-d screen, we must project down to 3d
                 return polyhedron.schlegel_projection(position=position)
+            elif polyhedron.dim() <= 3:
+                return polyhedron.affine_hull_projection(orthonormal=True, extend=True).projection()
             else:
                 return polyhedron.projection()

to this ticket? Then this won't be a regression in that way.

[gh-kliem]

Replying to jipilab:

[...]

In this case, this is a different ticket fixing the .show()/.plot() methods. This ticket fixes the schlegel projection which it does...

But you have introduced a change in this ticket, this is why I propose to "fix" it here. At least some temporary solution we can live with.

The default plotting method before this ticket for polytopes.permutahedron(4) was schlegel projection. You changed this to just projection.

[jipilab]

Ok sure.

But again: it was not doing an orthonormal projection, otherwise the minimal non-working example in the ticket description would not occur. It was a mere accident because the permutahedron realization given there is inscribed.

The above addition looks good to me. I would even prioritize it over polyhedron.dimension() == 4, because... I guess that whenever you have an object that is at most dimension 3, you would like to see "it" exactly how it is, and I actually would give an option to turn off the orthonormal if the user wants to. Then, if the object has dimension 4, it will do the schlegel projection (without respect at all to the ambient dimension, which make more sense to me).

I'll have another round of look at it.

[git]

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

​dae3bbf

Fixed default behavior and tweaks

[jipilab]

It should be all fixed now. I changed the default behavior for when the locus polyhedron is compact, to take the midpoint of the line segment. The rest is unchanged. Further, I exposed the option to remove orthonormality in the plotting...

Have a look at:

sage: p = polytopes.six_hundred_cell()                                                                                                                                                                             
sage: p.show(position=4,point={'size':0},frame=False)

[git]

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

​ec0edd7

Fixed higher dim/dirt

[jipilab]

... Ok. So ready for review again... It should be a small improvement on the plotting procedure. To me it makes more sense as of now. (see the internal def project function which decides what to do...)

[gh-kliem]

-        def project(polyhedron,ortho):
+        def project(polyhedron, ortho):

-        projection = project(self,orthonormal)
+        projection = project(self, orthonormal)

I propose the following doctest to show that this has been fixed::

sage: fcube = polytopes.hypercube(4)                                                                                                                                                
sage: tfcube = fcube.face_truncation(fcube.faces(0)[0])                                                                                                                             
sage: sp = tfcube.schlegel_projection()                                                                                                                                             
sage: for face in tfcube.faces(2): 
....:     vertices = face.ambient_Vrepresentation() 
....:     indices = [sp.coord_index_of(vector(x)) for x in vertices] 
....:     projected_vertices = [sp.transformed_coords[i] for i in indices] 
....:     assert Polyhedron(projected_vertices).dim() == 2 

This line here is hard to read due to the underscore. I would exchange it for something else:

ineq = [_ for _ in facet.ambient_Hrepresentation() if _.is_inequality()][0]

+        A, b = self.facet.as_polyhedron().affine_hull_projection(as_affine_map=True, orthonormal=True,extend=True)
-        A,b = self.facet.as_polyhedron().affine_hull_projection(as_affine_map=True, orthonormal=True, extend=True)

-        self.projection_point = vector(RDF,projection_point)
+        self.projection_point = vector(RDF, projection_point)

[gh-kliem]

Once done, you can put it on positive review on my behalf.

Thank you for fixing this. It's great to have this.

[gh-kliem]

With your ticket, how do I obtain a decent picture of polytopes.cyclic_polytope(4,10)?

Or maybe one should really choose better points on the moment curve for pictures. The pairwise euclidean distance of the vertices behaves awful.

[gh-kliem]

I guess you can basically forget the last comment. It is just to note that the current realization doesn't work for pictures.

Anyway, if you agree to the changes I suggested, I can also do them.

[git]

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

​f10d571

Schoenheit fixes + test

[jipilab]

Replying to gh-kliem:

With your ticket, how do I obtain a decent picture of polytopes.cyclic_polytope(4,10)?

Or maybe one should really choose better points on the moment curve for pictures. The pairwise euclidean distance of the vertices behaves awful.

There isn't really any good choice of parameters to get "a nice picture" of the cyclic polytope on the moment curve.

A better choice would be to essentially add the cyclic polytope construction on the trigonometric curve so that the coordinates all have the same "scale". ... another nice thing that would be easy to implement and that's not too hard.

[slelievre]

Typos, to be fixed here or in a follow-up ticket.

-        A different values of ``position`` changes the projection::
+        A different value of ``position`` changes the projection::

-        sees more than one facet resulting in a error::
+        sees more than one facet resulting in an error::