Opened 5 years ago

Polyhedron.get_integral_point — at Version 8

Reported by: Owned by: Mark Bell major sage-8.2 geometry Polyhedron, integral_points Vincent Delecroix, Thierry Monteil, Matthias Köppe Mark Bell N/A u/mcbell/polyhedron_get_integral_point 82a64d2f6eb533c16413c681a8d89c130283fa1c

This patch adds a method Polyhedron.get_integral_point(index) which returns the nth integral point in the polyhedron. It is equivalent to sorted(Polyhedron.integral_points())[index]. However when Polyhedron.integral_points_count() does not need to enumerate all of the integral points, for example in rational Polyhedra, neither does this method. Hence it can be significantly faster.

This method is useful for performing random sampling of integral points since it allows points to be chosen uniformly at random via:

index = randint(0, P.integral_points_count())
point = P.get_integral_point(index)

comment:1 Changed 5 years ago by Mark Bell

Branch: → u/mcbell/polyhedron_get_integral_point

comment:2 Changed 5 years ago by Mark Bell

Authors: → Mark Bell Vincent Delecroix Thierry Monteil Matthias Köppe added → d3b9b396c1dc9d6ca8c9478b01035454ebbe1f74 PLEASE CHANGE → geometry modified (diff) Polyhedron integral_points added new → needs_review PLEASE CHANGE → enhancement

comment:3 Changed 5 years ago by Vincent Delecroix

Mark,

As we discussed in Warwick, this method is rather unnatural (integer points are not naturally indexed) and is not needed for random sampling. To my mind, it would make more sense to have a new method .random_integral_point implemented along the same lines.

comment:4 follow-up:  5 Changed 5 years ago by Mark Bell

Sure, that would be a simple enough thing to change - should I open a new ticket?

comment:5 in reply to:  4 ; follow-up:  7 Changed 5 years ago by Vincent Delecroix

Sure, that would be a simple enough thing to change - should I open a new ticket?

Though, I do not think I was right: e.g. integral points are naturally indexed for the lexicographical ordering. In the description, you claim to follow the ordering of integral_points. Is it always lexicographic? I could not find any mention of ordering in the documentation. Note in particular that integral_points claims to have two implementations

Uses either the naive algorithm (iterate over a rectangular
bounding box) or triangulation + Smith form.

The ordering is independent of the method?

BTW, your branch u/mcbell/polyhedron_get_integral_point does not exist. Is that normal?

About get_integral_point vs random_integral_point you are free to do whatever you want with your ticket :-)

comment:6 Changed 5 years ago by git

Commit: d3b9b396c1dc9d6ca8c9478b01035454ebbe1f74 → 82a64d2f6eb533c16413c681a8d89c130283fa1c

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

 ​e020df5 Initial version of get_integral_point method. ​ca57e15 Added docstring and more comments. ​f771506 Fixed docstring whitespace. ​626e399 Fixed line widths and fencepost error in upper bound. ​c46460d Added sorting comments to docstring and random_integral_point method. ​82a64d2 Added missing import.

comment:7 in reply to:  5 Changed 5 years ago by Mark Bell

Thanks for the great comments. I think I've now managed to push my branch onto trac so you should be able to see the code.

Though, I do not think I was right: e.g. integral points are naturally indexed for the lexicographical ordering. In the description, you claim to follow the ordering of integral_points. Is it always lexicographic? I could not find any mention of ordering in the documentation. Note in particular that integral_points claims to have two implementations

Uses either the naive algorithm (iterate over a rectangular
bounding box) or triangulation + Smith form.

The ordering is independent of the method?

Ah, you are correct, you are not guaranteed to get lexicographical ordering when the triangulation method is used. In part this is caused by the fact that the points that are found in the different triangles are just merged together using a set. Perhaps L5288 of src/sage/geometrypolyhedron/base.py should become "return tuple(sorted(points))" to make the method deterministic. So the correct thing to say is that this method returns

sorted(self.integral_points())[index]

I've updated the docstring and above explanation to reflect this.

My new suggestion, which is in the latest commits, is to also add a Polyhedron.random_integral_point() method that just uses the snippet I posted originally.

comment:8 Changed 5 years ago by Mark Bell

Description: modified (diff)
Note: See TracTickets for help on using tickets.