A lot of polytope constructors are broken

[videlec]

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

Depends on #18211

CC: @nathanncohen

Component: geometry

Author: Vincent Delecroix

Branch/Commit: a58da00

Reviewer: Nathann Cohen

Issue created by migration from https://trac.sagemath.org/ticket/18213

[jdemeyer]

Description changed:

--- 
+++ 
@@ -7,10 +7,12 @@
 ...

-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 +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
-v2 = QQ(104455571357/27284200000)   # ~ 2 sqrt(2)
+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

• Now that Sage supports Polyhedron over number fields, we should use them! +Instead, we should use the base_ring argument (with appropriate defaults) and use base_ring(2).sqrt() instead.

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

[videlec]

Author: Vincent Delecroix

[videlec]

comment:3

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

[videlec]

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

[videlec]

Branch: u/vdelecroix/18123

[videlec]

Commit: e598e28

[videlec]

Description changed:

--- 
+++ 
@@ -1,4 +1,4 @@
-A lot of polytopes constructors in `sage.geometry.polyhedron.library`.
+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

@@ -15,4 +15,4 @@

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

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

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:5

Depends on #18211?

[videlec]

Dependencies: #18211

[videlec]

comment:6

Replying to @nathanncohen:

Depends on #18211?

Yes

---

New commits:

c44b53e

Trac 18213: fix doctests

[videlec]

Changed commit from e598e28 to c44b53e

[videlec]

Changed branch from u/vdelecroix/18123 to u/vdelecroix/18213

[videlec]

Description changed:

--- 
+++ 
@@ -16,3 +16,9 @@
 Instead, we should use the `base_ring` argument (with appropriate defaults) and use `base_ring(2).sqrt()` instead.

 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

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

6e7cf21

trac #18211: Seealsos

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from c44b53e to 6e7cf21

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

227e14a	Trac 18213: Seealsos
eef42af	Trac 18213: documentation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 6e7cf21 to eef42af

[videlec]

Description changed:

--- 
+++ 
@@ -22,3 +22,4 @@
 - #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

[jdemeyer]

comment:11

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?

[videlec]

comment:12

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

[videlec]

comment:13

Replying to @videlec:

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?

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from eef42af to c1804d5

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

c1804d5

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

[videlec]

comment:15

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

Vincent

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:16

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?

[videlec]

comment:17

Replying to @nathanncohen:

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!

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:18

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

[videlec]

Description changed:

--- 
+++ 
@@ -15,6 +15,10 @@

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:

[videlec]

comment:19

Replying to @nathanncohen:

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.

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:20

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

[videlec]

comment:21

Replying to @nathanncohen:

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

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from c1804d5 to 524f00e

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:23

O_o

[videlec]

comment:24

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

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:25

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

[videlec]

comment:26

Replying to @nathanncohen:

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

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 524f00e to fac3cff

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

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

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:28

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

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:29

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

[videlec]

comment:30

Replying to @nathanncohen:

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

[videlec]

comment:31

Replying to @nathanncohen:

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)

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:32

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

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:33

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.

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:34

Funny. There is no automorphism_group function for polytopes O_o

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:35

        OUTPUT: 

        EXAMPLES::

[videlec]

comment:36

Replying to @nathanncohen:

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.

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:37

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

[6bdad4c1-1e26-4f2f-a442-a01a2292c181]

comment:38

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?