Opened 2 years ago
Closed 2 years ago
#30015 closed defect (fixed)
Schlegel projection breaks convexity
Reported by:  JeanPhilippe Labbé  Owned by:  

Priority:  major  Milestone:  sage9.3 
Component:  geometry  Keywords:  polytope, schlegel 
Cc:  ghkliem, Laith Rastanawi  Merged in:  
Authors:  JeanPhilippe Labbé  Reviewers:  Jonathan Kliem 
Report Upstream:  N/A  Work issues:  
Branch:  f10d571 (Commits, GitHub, GitLab)  Commit:  f10d5714bc757db8ec8910359931fb1bcb0dbb29 
Dependencies:  Stopgaps: 
Description (last modified by )
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 deduplicating some code.
FOLLOWUP: Fix .plot() and .show() to have a better behaviour of smaller dimensional objects.
Change History (36)
comment:1 followup: 3 Changed 2 years ago by
comment:2 Changed 2 years ago by
Authors:  → JeanPhilippe Labbé 

Branch:  → u/jipilab/schlegel 
Commit:  → bed13a6539a1750c63250faee16342bf33bed8a1 
Description:  modified (diff) 
Status:  new → needs_review 
comment:3 followup: 4 Changed 2 years ago by
comment:4 Changed 2 years ago by
Replying to ghkliem:
Replying to ghkliem:
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 userfriendly, 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 followup: 6 Changed 2 years ago by
Status:  needs_review → 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/sitepackages/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.
comment:6 followup: 7 Changed 2 years ago by
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/sitepackages/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 rangeThat 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 Changed 2 years ago by
Replying to ghLaisRast:
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...)
comment:8 Changed 2 years ago by
Commit:  bed13a6539a1750c63250faee16342bf33bed8a1 → 2b6e5c5ebe3cbda4a6c33d44342eeb55a351ccb3 

Branch pushed to git repo; I updated commit sha1. New commits:
2b6e5c5  Changed method to get projection

comment:9 Changed 2 years ago by
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 2 years ago by
Commit:  2b6e5c5ebe3cbda4a6c33d44342eeb55a351ccb3 → 4342b529d11d548b4fca2ea229c927d2b364ab3e 

Branch pushed to git repo; I updated commit sha1. New commits:
4342b52  Fix doctests

comment:12 Changed 2 years ago by
Status:  needs_review → needs_work 

comment:13 Changed 2 years ago by
Status:  needs_work → needs_review 

comment:14 Changed 2 years ago by
Milestone:  sage9.2 → sage9.3 

comment:15 Changed 2 years ago by
# There is no 4d screen, we must project down to 3d
Either
3d
or4d
I would say.

+ if not locus_polyhedron.relative_interior_contains(projection_point): + raise ValueError("the chosen position is too large")
If you enterposition=1
, this is the error you see, which is somewhat strange. (Of course bullshit gives bullshit, but maybe we can be nice here.)
comment:17 followup: 20 Changed 2 years ago by
 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 2 years ago by
Status:  needs_review → needs_work 

comment:19 Changed 2 years ago by
Commit:  4342b529d11d548b4fca2ea229c927d2b364ab3e → c3df1d5855585ee893a86fbcd85d341fe4d2d58d 

Branch pushed to git repo; I updated commit sha1. New commits:
c3df1d5  Add some Errors and easy fixes

comment:20 followup: 22 Changed 2 years ago by
Description:  modified (diff) 

Status:  needs_work → needs_review 
Replying to ghkliem:
 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...
comment:21 Changed 2 years ago by
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): 1040 1040 elif polyhedron.dimension() == 4: 1041 1041 # There is no 4d screen, we must project down to 3d 1042 1042 return polyhedron.schlegel_projection(position=position) 1043 elif polyhedron.dim() <= 3: 1044 return polyhedron.affine_hull_projection(orthonormal=True, extend=True).projection() 1043 1045 else: 1044 1046 return polyhedron.projection()
to this ticket? Then this won't be a regression in that way.
comment:22 Changed 2 years ago by
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 2 years ago by
Ok sure.
But again: it was not doing an orthonormal projection, otherwise the minimal nonworking 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 2 years ago by
Commit:  c3df1d5855585ee893a86fbcd85d341fe4d2d58d → dae3bbfb2f6d65eeb5f012b064d05e130a5d4f82 

Branch pushed to git repo; I updated commit sha1. New commits:
dae3bbf  Fixed default behavior and tweaks

comment:25 Changed 2 years ago by
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 2 years ago by
Commit:  dae3bbfb2f6d65eeb5f012b064d05e130a5d4f82 → ec0edd73b6bf6e9d5ca52fb50a4cf04fe02c9983 

Branch pushed to git repo; I updated commit sha1. New commits:
ec0edd7  Fixed higher dim/dirt

comment:27 Changed 2 years ago by
... 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 2 years ago by
 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)
comment:29 Changed 2 years ago by
Reviewers:  → 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 followup: 34 Changed 2 years ago by
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 2 years ago by
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 2 years ago by
Commit:  ec0edd73b6bf6e9d5ca52fb50a4cf04fe02c9983 → f10d5714bc757db8ec8910359931fb1bcb0dbb29 

Branch pushed to git repo; I updated commit sha1. New commits:
f10d571  Schoenheit fixes + test

comment:33 Changed 2 years ago by
Status:  needs_review → positive_review 

comment:34 Changed 2 years ago by
Replying to ghkliem:
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 2 years ago by
Typos, to be fixed here or in a followup 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 2 years ago by
Branch:  u/jipilab/schlegel → f10d5714bc757db8ec8910359931fb1bcb0dbb29 

Resolution:  → fixed 
Status:  positive_review → closed 
Is this the bug you showed during your presentation at sd109?