#30015 closed defect (fixed)

Schlegel projection breaks convexity

Reported by: jipilab Owned by:
Priority: major Milestone: sage-9.3
Component: geometry Keywords: polytope, schlegel
Cc: gh-kliem, gh-LaisRast Merged in:
Authors: Jean-Philippe Labbé Reviewers: Jonathan Kliem
Report Upstream: N/A Work issues:
Branch: f10d571 (Commits, GitHub, GitLab) Commit: f10d5714bc757db8ec8910359931fb1bcb0dbb29
Dependencies: Stopgaps:

Status badges

Description (last modified by 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.

Change History (36)

comment:1 follow-up: Changed 17 months ago by gh-kliem

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

comment:2 Changed 14 months ago by jipilab

  • Authors set to Jean-Philippe Labbé
  • Branch set to u/jipilab/schlegel
  • Commit set to bed13a6539a1750c63250faee16342bf33bed8a1
  • Description modified (diff)
  • Status changed from new to needs_review

New commits:

26c03d1Fixed chirality in stereographic projection
bed13a6examples and pep8

comment:3 in reply to: ↑ 1 ; follow-up: Changed 14 months ago by gh-kliem

Replying to gh-kliem:

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

?

comment:4 in reply to: ↑ 3 Changed 14 months ago by 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.

comment:5 follow-up: Changed 14 months ago by jipilab

  • Status changed from needs_review to needs_work

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.

Last edited 14 months ago by jipilab (previous) (diff)

comment:6 in reply to: ↑ 5 ; follow-up: Changed 14 months ago by 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.

comment:7 in reply to: ↑ 6 Changed 14 months ago by 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...)

Last edited 14 months ago by jipilab (previous) (diff)

comment:8 Changed 14 months ago by git

  • Commit changed from bed13a6539a1750c63250faee16342bf33bed8a1 to 2b6e5c5ebe3cbda4a6c33d44342eeb55a351ccb3

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

2b6e5c5Changed method to get projection

comment:9 Changed 14 months ago by 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.

comment:10 Changed 14 months ago by git

  • Commit changed from 2b6e5c5ebe3cbda4a6c33d44342eeb55a351ccb3 to 4342b529d11d548b4fca2ea229c927d2b364ab3e

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

4342b52Fix doctests

comment:11 Changed 14 months ago by jipilab

  • Status changed from needs_work to needs_review

That should do it.

comment:12 Changed 14 months ago by jipilab

  • Status changed from needs_review to needs_work

comment:13 Changed 14 months ago by jipilab

  • Status changed from needs_work to needs_review

comment:14 Changed 13 months ago by mkoeppe

  • Milestone changed from sage-9.2 to sage-9.3

comment:15 Changed 13 months ago by 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.)

comment:16 Changed 13 months ago by gh-kliem

There is a trailing space in the preimage line.

comment:17 follow-up: Changed 13 months ago by 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.

comment:18 Changed 13 months ago by gh-kliem

  • Status changed from needs_review to needs_work

comment:19 Changed 13 months ago by git

  • Commit changed from 4342b529d11d548b4fca2ea229c927d2b364ab3e to c3df1d5855585ee893a86fbcd85d341fe4d2d58d

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

c3df1d5Add some Errors and easy fixes

comment:20 in reply to: ↑ 17 ; follow-up: Changed 13 months ago by jipilab

  • Description modified (diff)
  • Status changed from needs_work to needs_review

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...

Last edited 13 months ago by jipilab (previous) (diff)

comment:21 Changed 13 months ago by 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

    diff --git a/src/sage/geometry/polyhedron/base.py b/src/sage/geometry/polyhedron/base.py
    index 59eeb0aa41..f3551f6045 100644
    a b class Polyhedron_base(Element): 
    10401040            elif polyhedron.dimension() == 4:
    10411041                # There is no 4-d screen, we must project down to 3d
    10421042                return polyhedron.schlegel_projection(position=position)
     1043            elif polyhedron.dim() <= 3:
     1044                return polyhedron.affine_hull_projection(orthonormal=True, extend=True).projection()
    10431045            else:
    10441046                return polyhedron.projection()

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

comment:22 in reply to: ↑ 20 Changed 13 months ago by 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.

comment:23 Changed 13 months ago by 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.

comment:24 Changed 13 months ago by git

  • Commit changed from c3df1d5855585ee893a86fbcd85d341fe4d2d58d to dae3bbfb2f6d65eeb5f012b064d05e130a5d4f82

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

dae3bbfFixed default behavior and tweaks

comment:25 Changed 13 months ago by 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)

comment:26 Changed 13 months ago by git

  • Commit changed from dae3bbfb2f6d65eeb5f012b064d05e130a5d4f82 to ec0edd73b6bf6e9d5ca52fb50a4cf04fe02c9983

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

ec0edd7Fixed higher dim/dirt

comment:27 Changed 13 months ago by 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...)

comment:28 Changed 13 months ago by 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)
Last edited 13 months ago by gh-kliem (previous) (diff)

comment:29 Changed 13 months ago by gh-kliem

  • Reviewers set to Jonathan Kliem

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

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

comment:30 follow-up: Changed 13 months ago by 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.

comment:31 Changed 13 months ago by 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.

comment:32 Changed 12 months ago by git

  • Commit changed from ec0edd73b6bf6e9d5ca52fb50a4cf04fe02c9983 to f10d5714bc757db8ec8910359931fb1bcb0dbb29

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

f10d571Schoenheit fixes + test

comment:33 Changed 12 months ago by jipilab

  • Status changed from needs_review to positive_review

comment:34 in reply to: ↑ 30 Changed 12 months ago by 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.

comment:35 Changed 12 months ago by 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::

comment:36 Changed 12 months ago by vbraun

  • Branch changed from u/jipilab/schlegel to f10d5714bc757db8ec8910359931fb1bcb0dbb29
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.