Opened 3 years ago
Closed 2 years ago
#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: |
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 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: 3 Changed 3 years ago by
comment:2 Changed 2 years ago by
Authors: | → Jean-Philippe Labbé |
---|---|
Branch: | → u/jipilab/schlegel |
Commit: | → bed13a6539a1750c63250faee16342bf33bed8a1 |
Description: | modified (diff) |
Status: | new → needs_review |
comment:3 follow-up: 4 Changed 2 years ago by
comment:4 Changed 2 years ago by
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: 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/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.
comment:6 follow-up: 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/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 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 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...)
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: | sage-9.2 → sage-9.3 |
---|
comment:15 Changed 2 years ago by
# There is no 4-d screen, we must project down to 3d
Either
3-d
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 follow-up: 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 follow-up: 22 Changed 2 years ago by
Description: | modified (diff) |
---|---|
Status: | needs_work → 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...
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 4-d 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 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 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 follow-up: 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 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 2 years ago by
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 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?