Opened 6 years ago

Closed 6 years ago

#18213 closed defect (fixed)

A lot of polytope constructors are broken

Reported by: vdelecroix Owned by:
Priority: major Milestone: sage-6.7
Component: geometry Keywords:
Cc: ncohen Merged in:
Authors: Vincent Delecroix Reviewers: Nathann Cohen
Report Upstream: N/A Work issues:
Branch: a58da00 (Commits) Commit: a58da00e4907251c2ab47d7af1109da31901ee36
Dependencies: #18211 Stopgaps:

Description (last modified by 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],
    ...
    
  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

Change History (52)

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

  • Description modified (diff)

comment:2 Changed 6 years ago by ncohen

  • Cc ncohen added

comment:3 Changed 6 years ago by vdelecroix

  • Authors set to Vincent Delecroix

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

comment:4 Changed 6 years ago by vdelecroix

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

New commits:

5110cc9Trac 18211: use LattE to compute the Ehrhart polynomial
cf2a10fTrac 18211: another example using Birkhoff_polytope
ca8da5aTrac 18211: fixed spelling LattE integral*e*
d52053cTrac 18211: check error + options
0f5122bTrac 18211: fix command line arguments
b877932Trac 18211: fix arguments + doctests + cwd=SAGE_TMP
3b30bc8Trac 18213: rewrite the polyhedron library
e598e28Trac 18123: fix doctests

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

Depends on #18211?

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

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

Replying to ncohen:

Depends on #18211?

Yes


New commits:

c44b53eTrac 18213: fix doctests

comment:7 Changed 6 years ago by vdelecroix

  • Description modified (diff)

comment:8 Changed 6 years ago by git

  • Commit changed from c44b53efbe0fa899951481b96c4ab0a15ebb88e3 to 6e7cf214ea1ee23be11060e0327d20e2ff1a35c8

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

6e7cf21trac #18211: Seealsos

comment:9 Changed 6 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:

227e14aTrac 18213: Seealsos
eef42afTrac 18213: documentation

comment:10 Changed 6 years ago by vdelecroix

  • Description modified (diff)

comment:11 follow-up: Changed 6 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: Changed 6 years ago by vdelecroix

  • Status changed from needs_review to needs_work

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

comment:13 in reply to: ↑ 12 Changed 6 years ago by 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?

comment:14 Changed 6 years ago by git

  • Commit changed from eef42afa3e3229534f2ae2d3303fbb176eaecdd8 to c1804d5879eedab9140d07096bfd1ad29ae75fa0

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

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

comment:15 Changed 6 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: Changed 6 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 6 years ago by 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!

comment:18 follow-up: Changed 6 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 6 years ago by vdelecroix

  • Description modified (diff)

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.

comment:20 follow-up: Changed 6 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 6 years ago by 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>.

comment:22 Changed 6 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:

a9ce461trac #18211: print stderr when verbose=False
9307916Trac 18211: "in in" -> "in"
a6d03cdTrac 18211: any option to the command (using **kwds)
c935177Trac 18211: handle bad output from count
7215956Trac 18211: write explicitely the options in function declaration
74ef0f1trac #18211: Nothing important
ec0e946Trac 18213: merge sage-6.7.beta1
0b87ff0Trac 18213: rewrite the polyhedron library
2ddcdefTrac 18213: fix doctests
524f00eTrac 18213: Seealsos

comment:23 Changed 6 years ago by ncohen

O_o

comment:24 Changed 6 years ago by vdelecroix

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

comment:25 follow-up: Changed 6 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.
  • 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

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

comment:26 in reply to: ↑ 25 ; follow-up: Changed 6 years ago by 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

comment:27 Changed 6 years ago by git

  • Commit changed from 524f00ea3d1d3734a9915cd6052d0bf677fdb6e3 to fac3cffdb326dad1d3c61045d1c33044c99b8b32

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

7eb2326trac #18213: Review
fac3cffTrac 18213: (review) more documentation

comment:28 in reply to: ↑ 26 ; follow-up: Changed 6 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: Changed 6 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: Changed 6 years ago by 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

comment:31 in reply to: ↑ 28 Changed 6 years ago by 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)

comment:32 in reply to: ↑ 30 Changed 6 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 6 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 6 years ago by ncohen

Funny. There is no automorphism_group function for polytopes O_o

comment:35 Changed 6 years ago by ncohen

        OUTPUT: 

        EXAMPLES::

comment:36 Changed 6 years ago by 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.

comment:37 Changed 6 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 6 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 6 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 6 years ago by 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)!

comment:41 Changed 6 years ago by git

  • Commit changed from fac3cffdb326dad1d3c61045d1c33044c99b8b32 to 1a82754ea33046384dbefc6045de5bef0e0a2261

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

21e70a7Trac 18213: review commit
1a82754Trac 18213: doctest fix

comment:42 Changed 6 years ago by vdelecroix

Patchbot should be happy after that.

comment:43 follow-up: Changed 6 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 6 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 6 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 6 years ago by ncohen

  • Status changed from needs_review to needs_work

comment:47 Changed 6 years ago by git

  • Commit changed from 1a82754ea33046384dbefc6045de5bef0e0a2261 to a58da00e4907251c2ab47d7af1109da31901ee36

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

120c1e9Trac 18213: more abs tol in doctests
d09d8a5Trac 18211: forgot a # optional in the doctest
a58da00Trac 18213: merge failed doctest from #18211

comment:48 in reply to: ↑ 43 Changed 6 years ago by 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.

comment:49 Changed 6 years ago by vdelecroix

  • Status changed from needs_work to needs_review

comment:50 Changed 6 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 6 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 6 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.