Opened 7 years ago

Closed 7 years ago

# A lot of polytope constructors are broken

Reported by: Owned by: vdelecroix major sage-6.7 geometry ncohen Vincent Delecroix Nathann Cohen N/A a58da00 a58da00e4907251c2ab47d7af1109da31901ee36 #18211

A lot of polytopes constructors in sage.geometry.polyhedron.library. For example

1. The Rhombicuboctahedron does not return what it should (as Python division was thought as Sage integer divisions). There is in the code
verts = [ [-3/2, -1/2, -1/2], [-3/2, -1/2, 1/2], [-3/2, 1/2, -1/2],
...

1. The great_rhombicuboctahedron is defined over QQ but it should be defined over QQ[sqrt(2)]! There are in two places rough approximation of sqrt(2) in the code
v1 = QQ(131739771357/54568400000)   # ~ sqrt(2) + 1 but actually 2 due to Python division
v2 = QQ(104455571357/27284200000)   # ~ 2 sqrt(2) but actually 3 due to Python division


Instead, we should use the base_ring argument (with appropriate defaults) and use base_ring(2).sqrt() instead.

1. The functions unrelated to construction of Polytopes are moved out of the class Polytopes:
• Polytopes.orthonormal_1, Polytopes.project_1 will be renamed respectively zero_sum_projection and project_points
• the one line Polytopes._pfunc is just removed

While we're at it, remove the deprecations from #11634.

During the process I discovered two annoying bugs:

• #18214 Bug in volume computation of polyhedron
• #18220 Bug when creating a polyhedron with coefficients in RR

Moreover, we can get a great speed up with the following because many polytopes have now coordinates in a quadratic number fields:

• #18215: Huge speed up for hash of quadratic number field elements
• #18226: Native _mpfr_ method on quadratic number field elements

### comment:1 in reply to: ↑ description Changed 7 years ago by jdemeyer

• Description modified (diff)

### comment:3 Changed 7 years ago by vdelecroix

• Authors set to Vincent Delecroix

I have a big commit to push in few minutes...

### comment:4 Changed 7 years ago by vdelecroix

• Branch set to u/vdelecroix/18123
• Description modified (diff)
• Status changed from new to needs_review

New commits:

 ​5110cc9 Trac 18211: use LattE to compute the Ehrhart polynomial ​cf2a10f Trac 18211: another example using Birkhoff_polytope ​ca8da5a Trac 18211: fixed spelling LattE integral*e* ​d52053c Trac 18211: check error + options ​0f5122b Trac 18211: fix command line arguments ​b877932 Trac 18211: fix arguments + doctests + cwd=SAGE_TMP ​3b30bc8 Trac 18213: rewrite the polyhedron library ​e598e28 Trac 18123: fix doctests

### comment:5 follow-up: ↓ 6 Changed 7 years ago by ncohen

Depends on #18211?

### comment:6 in reply to: ↑ 5 Changed 7 years ago by vdelecroix

• Branch changed from u/vdelecroix/18123 to u/vdelecroix/18213
• Commit changed from e598e289689ad8c294500303c89d5c1a0de2a060 to c44b53efbe0fa899951481b96c4ab0a15ebb88e3
• Dependencies set to #18211

Depends on #18211?

Yes

New commits:

 ​c44b53e Trac 18213: fix doctests

### comment:7 Changed 7 years ago by vdelecroix

• Description modified (diff)

### comment:8 Changed 7 years ago by git

• Commit changed from c44b53efbe0fa899951481b96c4ab0a15ebb88e3 to 6e7cf214ea1ee23be11060e0327d20e2ff1a35c8

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

 ​6e7cf21 trac #18211: Seealsos

### comment:9 Changed 7 years ago by git

• Commit changed from 6e7cf214ea1ee23be11060e0327d20e2ff1a35c8 to eef42afa3e3229534f2ae2d3303fbb176eaecdd8

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

 ​227e14a Trac 18213: Seealsos ​eef42af Trac 18213: documentation

### comment:10 Changed 7 years ago by vdelecroix

• Description modified (diff)

### comment:11 follow-up: ↓ 12 Changed 7 years ago by jdemeyer

The golden_ratio() can obviously be defined in terms of sqrt(5), so I don't see the need for that. I'm also really wondering whether you need the custom sqrt() function, doesn't QuadraticField(n) work?

### comment:12 in reply to: ↑ 11 ; follow-up: ↓ 13 Changed 7 years ago by vdelecroix

• Status changed from needs_review to needs_work

The golden_ratio() can obviously be defined in terms of sqrt(5), so I don't see the need for that. I'm also really wondering whether you need the custom sqrt() function, doesn't QuadraticField(n) work?

Right. The problem of embedding should be done elsewhere (my version embedds in QQbar whereas the generic QuadraticField(n) embedds into RLF).

### comment:13 in reply to: ↑ 12 Changed 7 years ago by vdelecroix

The golden_ratio() can obviously be defined in terms of sqrt(5), so I don't see the need for that. I'm also really wondering whether you need the custom sqrt() function, doesn't QuadraticField(n) work?

Right. The problem of embedding should be done elsewhere (my version embedds in QQbar whereas the generic QuadraticField(n) embedds into RLF).

Here is one reason. The generator name of QuadraticField(2) is a and it better be sqrt2 in this case. For the golden ratio, the advantage is that its get printed as phi. In particular you have

sage: polytopes.icosahedron().volume()
5/6*phi + 5/6


The answer 5/6*a + 5/6 would have been terrible here. I have nothing against using phi = (sqrt5+1)/2. So I will at least remove this one. For the quadratic fields, one option is to globally make this change in another ticket. What do you think?

### comment:14 Changed 7 years ago by git

• Commit changed from eef42afa3e3229534f2ae2d3303fbb176eaecdd8 to c1804d5879eedab9140d07096bfd1ad29ae75fa0

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

 ​c1804d5 Trac 18213: (review) remove sqrt/golden_ratio + doc

### comment:15 Changed 7 years ago by vdelecroix

• Status changed from needs_work to needs_review

All right. I removed the sqrt and golden_ratio function. Needs review!

Vincent

### comment:16 follow-up: ↓ 17 Changed 7 years ago by ncohen

Vincent: the diff of your main commit is very, very very hard to read. Could you split in order to make it easier to understand?

### comment:17 in reply to: ↑ 16 Changed 7 years ago by vdelecroix

Vincent: the diff of your main commit is very, very very hard to read. Could you split in order to make it easier to understand?

Not sure. It is basically a complete rewrite. I added a ton of documentation, and modify many functions to fit with the actual definition!

### comment:18 follow-up: ↓ 19 Changed 7 years ago by ncohen

Okay. Then could you at least change the description of this ticket to explain all changes that you made? (you apparently reated new functions, for example.. Why?)

### comment:19 in reply to: ↑ 18 Changed 7 years ago by vdelecroix

• Description modified (diff)

Okay. Then could you at least change the description of this ticket to explain all changes that you made? (you apparently reated new functions, for example.. Why?)

done in the ticket description.

### comment:20 follow-up: ↓ 21 Changed 7 years ago by ncohen

About golden_ratio and sqrt: shouldn't they be in number_field_element_quadratic instead ?

About zero_sum_projection what about having it in matrix.<tab> ?

Nathann

### comment:21 in reply to: ↑ 20 Changed 7 years ago by vdelecroix

About golden_ratio and sqrt: shouldn't they be in number_field_element_quadratic instead ?

Removed in ​c1804d5.

About zero_sum_projection what about having it in matrix.<tab> ?

It is very specific. I am not sure how I would call it in matrix.<tab>.

### comment:22 Changed 7 years ago by git

• Commit changed from c1804d5879eedab9140d07096bfd1ad29ae75fa0 to 524f00ea3d1d3734a9915cd6052d0bf677fdb6e3

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

 ​a9ce461 trac #18211: print stderr when verbose=False ​9307916 Trac 18211: "in in" -> "in" ​a6d03cd Trac 18211: any option to the command (using **kwds) ​c935177 Trac 18211: handle bad output from count ​7215956 Trac 18211: write explicitely the options in function declaration ​74ef0f1 trac #18211: Nothing important ​ec0e946 Trac 18213: merge sage-6.7.beta1 ​0b87ff0 Trac 18213: rewrite the polyhedron library ​2ddcdef Trac 18213: fix doctests ​524f00e Trac 18213: Seealsos

### comment:23 Changed 7 years ago by ncohen

O_o

### comment:24 Changed 7 years ago by vdelecroix

Cleaner version based on sage-6.7.beta1 and the complete #18211.

### comment:25 follow-up: ↓ 26 Changed 7 years ago by ncohen

Here are some comments. Please understand that doing everything at once makes things *MUCH* harder to understand for the reviewer.

Improving the implementation and improving the doc could have been two separate commits. The work that you do on each function could have been a separate commit. Removing deprecations could have been a separate commit. There are also things that anybody can review (only code) and things for which you need to be used to the tools you deal with.

• The projection is not the orthogona projection I expected:
sage: sage: zero_sum_projection(3)*vector([1,0,0])
(0.7071067811865475, 0.4082482904638631)

You can either change the matrix or emphasize in the docstring that the projection is not unique. Something like "returns a (d-1)xd matrix of rank d whose kernel contain (1,...,1). At first sight, I expected the function to return a matrix of dimensions d\times d. I updated the docstring on this matter.
• I do not understand the paragraph in the doc of project_points but that's probably because I don't know the maths behind.
• Compatible?
This projection is the only one compatible with the restriction to the
first coordinates

• You should probably add a # tol after this
Its volume is \sqrt{d+1} / d!

Or even better, check that the difference between the two is 0 (still, with some tolerance)
• At this point I still do not know if it is very wise to use this definition of the projection, but it would be cool if every function which uses this projection (like !Simplex) could redirect toward the matrix function. This way people would have a chance to learn what exactly they get as a result.
sage: polytopes.icosahedron(base_ring=QQ)
...
TypeError: unable to convert 1/4*sqrt(5) + 1/4 to a rational


I added a small commit at public/18213

Nathann

Last edited 7 years ago by ncohen (previous) (diff)

### comment:26 in reply to: ↑ 25 ; follow-up: ↓ 28 Changed 7 years ago by vdelecroix

Here are some comments. Please understand that doing everything at once makes things *MUCH* harder to understand for the reviewer.

I am very sorry.

• The projection is not the orthogona projection I expected:

Which one did you expected? I improved the documentation.

• At this point I still do not know if it is very wise to use this definition of the projection, but it would be cool if every function which uses this projection (like !Simplex) could redirect toward the matrix function. This way people would have a chance to learn what exactly they get as a result.

I see. But I am really not confortable in moving this to matrix.<tab> as well. One possibility is to move it back to polytopes.<tab>. What do you think?

sage: polytopes.icosahedron(base_ring=QQ)
...
TypeError: unable to convert 1/4*sqrt(5) + 1/4 to a rational


Is it not clear enough? The coordinates of the icosahedron need to be defined in QQ[sqrt(5)], hence the error. The following is ok (but very slow, ~3secs on my computer)

sage: I = polytopes.icosahedron(base_ring=AA)


Apart adding your example to the documentation, I do not see much what I can do (see the last commit). A deprecation?

Vincent

### comment:27 Changed 7 years ago by git

• Commit changed from 524f00ea3d1d3734a9915cd6052d0bf677fdb6e3 to fac3cffdb326dad1d3c61045d1c33044c99b8b32

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

 ​7eb2326 trac #18213: Review ​fac3cff Trac 18213: (review) more documentation

### comment:28 in reply to: ↑ 26 ; follow-up: ↓ 31 Changed 7 years ago by ncohen

Hello,

Which one did you expected? I improved the documentation.

Well, I expected an orthogonal projection on this hyperplane. And I probably expected the base of the hyperplane to be the projection of d-1 vectors of the base from the original space I guess.

• At this point I still do not know if it is very wise to use this definition of the projection, but it would be cool if every function which uses this projection (like !Simplex) could redirect toward the matrix function. This way people would have a chance to learn what exactly they get as a result.

I see. But I am really not confortable in moving this to matrix.<tab> as well. One possibility is to move it back to polytopes.<tab>. What do you think?

All I was suggesting in the message above was to add a link in the doc. The message you added in project_points is perfect

The projection used is the matrix given by :func:zero_sum_projection.


To me it should appear whenever there is a 'project' argument in the function.

Is it not clear enough? The coordinates of the icosahedron need to be defined in QQ[sqrt(5)], hence the error.

Precisely: Could you add somewhere in the doc that the field must contain QQ(sqrt(5))?

Nathann

### comment:29 follow-up: ↓ 30 Changed 7 years ago by ncohen

the function dodecahedron does nothing with its base_ring argument. Also, why should the projection function return 'None' instead of '[]' when it is called with no argument?

Nathann

### comment:30 in reply to: ↑ 29 ; follow-up: ↓ 32 Changed 7 years ago by vdelecroix

the function dodecahedron does nothing with its base_ring argument. Also, why should the projection function return 'None' instead of '[]' when it is called with no argument?

You corrected it in 7eb2326

### comment:31 in reply to: ↑ 28 Changed 7 years ago by vdelecroix

Hello,

Which one did you expected? I improved the documentation.

Well, I expected an orthogonal projection on this hyperplane. And I probably expected the base of the hyperplane to be the projection of d-1 vectors of the base from the original space I guess.

The description is explicit in the doc (commit ​fac3cff)

• At this point I still do not know if it is very wise to use this definition of the projection, but it would be cool if every function which uses this projection (like !Simplex) could redirect toward the matrix function. This way people would have a chance to learn what exactly they get as a result.

I see. But I am really not confortable in moving this to matrix.<tab> as well. One possibility is to move it back to polytopes.<tab>. What do you think?

All I was suggesting in the message above was to add a link in the doc. The message you added in project_points is perfect

The projection used is the matrix given by :func:zero_sum_projection.


To me it should appear whenever there is a 'project' argument in the function.

Will do!

Is it not clear enough? The coordinates of the icosahedron need to be defined in QQ[sqrt(5)], hence the error.

Precisely: Could you add somewhere in the doc that the field must contain QQ(sqrt(5))?

There is one example at the end of the doc (from fac3cff)

### comment:32 in reply to: ↑ 30 Changed 7 years ago by ncohen

You corrected it in 7eb2326

Arggggg... I was looking at the original commit again. I was wondering why it has been removed >_<

The description is explicit in the doc (commit ​fac3cff)

Yes yes now it's good!

There is one example at the end of the doc (from fac3cff)

It happens several times though. Could you add it in the doc of base_ring? You already talk about this field there. It could be something like "If set to None then it will be QQ[\phi], otherwise it must be a field that contains it".

Nathann

### comment:33 Changed 7 years ago by ncohen

The 600-cell code would look a bit better if you were building all points at once instead of interrupting the flow to decide which ring you should work with.

### comment:34 Changed 7 years ago by ncohen

Funny. There is no automorphism_group function for polytopes O_o

### comment:35 Changed 7 years ago by ncohen

        OUTPUT:

EXAMPLES::


### comment:36 Changed 7 years ago by vdelecroix

Funny. There is no automorphism_group function for polytopes O_o

Yes! This is very sad... But you have to be careful about what it means. There is the combinatorial automorphism group and the isometry group. I have no idea how to implement the latter efficiently.

### comment:37 Changed 7 years ago by ncohen

Most probably the ugliest thing I ever saw

sage: Sequence([Graph()]).universe()
<class 'sage.graphs.graph.Graph'>
sage: Sequence([Graph(),1]).universe()
Category of objects
sage: Sequence([1]).universe()
Integer Ring


### comment:38 Changed 7 years ago by ncohen

HMmmmmmmm... I don't exactly know what the isometry group is, but let's try anyway: what about defining a complete graph on your points, in which each edge has a color associated to its length. Wouldn't the automorphism group of that be what you want?

### comment:39 Changed 7 years ago by ncohen

Okayyyyyyyyyyyyyy. End of the review. Nothing else to add ;-)

Very good job. This code needed to be cleaned, and it is much better now.

But please, next time: smaller patches. You waste your reviewers' health :-P

Nathann

### comment:40 Changed 7 years ago by vdelecroix

HMmmmmmmm... I don't exactly know what the isometry group is, but let's try anyway: what about defining a complete graph on your points, in which each edge has a color associated to its length. Wouldn't the automorphism group of that be what you want?

That is a valid definition. But the standard one is the set of orthogonal matrices that preserve the polyhedron. And it would be nice for isometry_group to be a finite matrix group (and not permutations)!

### comment:41 Changed 7 years ago by git

• Commit changed from fac3cffdb326dad1d3c61045d1c33044c99b8b32 to 1a82754ea33046384dbefc6045de5bef0e0a2261

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

 ​21e70a7 Trac 18213: review commit ​1a82754 Trac 18213: doctest fix

### comment:42 Changed 7 years ago by vdelecroix

Patchbot should be happy after that.

### comment:43 follow-up: ↓ 48 Changed 7 years ago by ncohen

Some broken doctests in the patchbot's report. Also, can you be sure that the sqrt(6.) / factorial(5) are not subject to numerical noise? There is no #tol on them right now.

Nathann

### comment:44 Changed 7 years ago by ncohen

Hmmmmmm.. I can't seem to find on google an implementation of an isometry group function anywhere O_o

### comment:45 Changed 7 years ago by jkeitel

Hi Vincent and Nathann

Just a quick comment, as I'm not sure how useful this will be for general polyhedra: There is a method for lattice polytopes: LatticePolytope_PPL (in sage.geometry.polyhedron.ppl_lattice_polytope) has a method lattice_automorphism_group()

Cheers, Jan

### comment:46 Changed 7 years ago by ncohen

• Status changed from needs_review to needs_work

### comment:47 Changed 7 years ago by git

• Commit changed from 1a82754ea33046384dbefc6045de5bef0e0a2261 to a58da00e4907251c2ab47d7af1109da31901ee36

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

 ​120c1e9 Trac 18213: more abs tol in doctests ​d09d8a5 Trac 18211: forgot a # optional in the doctest ​a58da00 Trac 18213: merge failed doctest from #18211

### comment:48 in reply to: ↑ 43 Changed 7 years ago by vdelecroix

Some broken doctests in the patchbot's report. Also, can you be sure that the sqrt(6.) / factorial(5) are not subject to numerical noise? There is no #tol on them right now.

### comment:49 Changed 7 years ago by vdelecroix

• Status changed from needs_work to needs_review

### comment:50 Changed 7 years ago by ncohen

• Reviewers set to Nathann Cohen
• Status changed from needs_review to positive_review

Okayyyyyyy. Then positive_review!

Nathann

### comment:51 Changed 7 years ago by vdelecroix

• Summary changed from A lot of polytopes constructor are broken to A lot of polytope constructors are broken

### comment:52 Changed 7 years ago by vbraun

• Branch changed from u/vdelecroix/18213 to a58da00e4907251c2ab47d7af1109da31901ee36
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.