A lot of polytope constructors are broken

[vdelecroix]

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

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

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

[vdelecroix]

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

[vdelecroix]

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

[ncohen]

Depends on #18211?

[vdelecroix]

Replying to ncohen:

Depends on #18211?

Yes

---

New commits:

​c44b53e

Trac 18213: fix doctests

[git]

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

​6e7cf21

trac #18211: Seealsos

[git]

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

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

[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?

[vdelecroix]

Replying to 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?

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

[vdelecroix]

Replying to vdelecroix:

Replying to 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?

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?

[git]

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

​c1804d5

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

[vdelecroix]

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

Vincent

[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?

[vdelecroix]

Replying to 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?

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

[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?)

[vdelecroix]

Replying to 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?)

done in the ticket description.

[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

[vdelecroix]

Replying to ncohen:

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

[git]

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

[ncohen]

O_o

[vdelecroix]

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

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

• What about this error?

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

[vdelecroix]

Replying to ncohen:

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?

• What about this error?

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

[git]

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

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

[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

[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

[vdelecroix]

Replying to 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?

You corrected it in 7eb2326

[vdelecroix]

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

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)

[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

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

[ncohen]

Funny. There is no automorphism_group function for polytopes O_o

[ncohen]

        OUTPUT: 

        EXAMPLES::

[vdelecroix]

Replying to ncohen:

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.

[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

[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?

[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

[vdelecroix]

Replying to 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?

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)!

[git]

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

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

[vdelecroix]

Patchbot should be happy after that.

[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

[ncohen]

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

[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

[git]

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

[vdelecroix]

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

I don't know... Patchbot did not complain about this.

[ncohen]

Okayyyyyyy. Then positive_review!

Nathann