Opened 2 years ago
Closed 23 months ago
#24451 closed enhancement (fixed)
Polyhedron.get_integral_point
Reported by:  mcbell  Owned by:  

Priority:  major  Milestone:  sage8.2 
Component:  geometry  Keywords:  Polyhedron, integral_points 
Cc:  vdelecroix, tmonteil, mkoeppe  Merged in:  
Authors:  Mark Bell  Reviewers:  Vincent Delecroix, Travis Scrimshaw 
Report Upstream:  N/A  Work issues:  
Branch:  73fbd05 (Commits)  Commit:  73fbd0571a629071bceb1300d699f0b826e983c1 
Dependencies:  Stopgaps: 
Description (last modified by )
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)
Change History (36)
comment:1 Changed 2 years ago by
 Branch set to u/mcbell/polyhedron_get_integral_point
comment:2 Changed 2 years ago by
 Cc vdelecroix tmonteil mkoeppe added
 Commit set to d3b9b396c1dc9d6ca8c9478b01035454ebbe1f74
 Component changed from PLEASE CHANGE to geometry
 Description modified (diff)
 Keywords Polyhedron integral_points added
 Status changed from new to needs_review
 Type changed from PLEASE CHANGE to enhancement
comment:3 Changed 2 years ago by
comment:4 followup: ↓ 5 Changed 2 years ago by
Sure, that would be a simple enough thing to change  should I open a new ticket?
comment:5 in reply to: ↑ 4 ; followup: ↓ 7 Changed 2 years ago by
Replying to mcbell:
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 2 years ago by
 Commit changed from d3b9b396c1dc9d6ca8c9478b01035454ebbe1f74 to 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 2 years ago by
Replying to vdelecroix:
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 thatintegral_points
claims to have two implementationsUses 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 2 years ago by
 Description modified (diff)
comment:9 Changed 2 years ago by
Concerning, P.integral_points()
, I think that it should even be an iterator (as most things in Python3). Hence I am against sorting. If the user needs a sorted list, she can just apply a sort on the result.
For polytopes "with small interior", I believe that it is much faster to compute the set of points first and then do random sampling (while keeping the list in some cache). But let us keep that for later on.
For more concrete comments: when the polytope has no integral points the method random_integral_point
should raise an EmptySetError
(from sage.categories.sets_cat
). And in both methods, the error have to be doctested. It should look like
sage: my_polytope_without_integral_points.random_integral_point() Traceback (most recent call last): ... EmptySetError: errormessage
(you can have a look around in the code)
comment:10 Changed 2 years ago by
 Commit changed from 82a64d2f6eb533c16413c681a8d89c130283fa1c to 11aabdd974fba81accf9966584a159fba8621c2e
Branch pushed to git repo; I updated commit sha1. New commits:
11aabdd  Added check to make sure there are integral points before sampling and doctests for the different errors.

comment:11 Changed 2 years ago by
For random sampling you should not use the directly the random
module but
sage: from sage.misc.randstate import current_randstate sage: random = current_randstate().python_random() sage: random.randint(0,5) 0 sage: random.randint(0,5) 5
The reason is that Sage has to take care of a lot of random generators coming from different libraries (Python, gmp, GAP, PARI/GP, etc). All of them are initialized through a unique function set_random_seed
. As a consequence, all the generators has to be taken from the current_randstate
.
You can compare
sage: set_random_seed(5); [current_randstate().python_random().randint(0,10) for _ in range(5)] [0, 5, 2, 3, 5] sage: set_random_seed(5); [current_randstate().python_random().randint(0,10) for _ in range(5)] [0, 5, 2, 3, 5]
versus
sage: import random sage: set_random_seed(5); [random.randint(0,10) for _ in range(5)] [9, 8, 3, 8, 2] sage: set_random_seed(5); [random.randint(0,10) for _ in range(5)] [9, 6, 1, 8, 5]
comment:12 Changed 2 years ago by
 Commit changed from 11aabdd974fba81accf9966584a159fba8621c2e to 4e76f7d28eb839bac153ea682e50707c5bdb4754
Branch pushed to git repo; I updated commit sha1. New commits:
4e76f7d  Using current_randstate() instead of random.

comment:13 Changed 2 years ago by
 Commit changed from 4e76f7d28eb839bac153ea682e50707c5bdb4754 to 99683415186c9890bc065fe0883f383fd2bd0019
Branch pushed to git repo; I updated commit sha1. New commits:
9968341  Fixed upper and lower bounds to be integers.

comment:14 Changed 2 years ago by
For consistency, could you make the output of get_integral_point
an immutable vector
sage: P = polytopes.icosahedron() sage: P.integral_points()[0].is_mutable() False sage: P.get_integral_point(0).is_mutable() True
(Just use the method .set_immutable()
of vectors.)
comment:15 followup: ↓ 17 Changed 2 years ago by
Why do you consider the inequality [lower] + [0] * i + [1] + [0] * (D  i  1)
. Shouldn't this be redundant?
I believe that it would be faster to just add one equality at a time (we can possibly speed that up later on). Namely, I would rather keep the polyhedron self
intersected with extra_eqns
along the loop. In pseudo code
P = self for i in {0, ..., dimension1}: P_left = P intersection left_guess count_left = P_left.count() if count_left < index: P = P_left else: P = P intersection right_guess index = count_left
comment:16 Changed 2 years ago by
 Reviewers set to Vincent Delecroix
comment:17 in reply to: ↑ 15 Changed 2 years ago by
Replying to vdelecroix:
Why do you consider the inequality
[lower] + [0] * i + [1] + [0] * (D  i  1)
. Shouldn't this be redundant?
So this was needed precisely because of the fact that previously found inequalities were not added to the polyhedron.
I believe that it would be faster to just add one equality at a time (we can possibly speed that up later on). Namely, I would rather keep the polyhedron
self
intersected withextra_eqns
along the loop. In pseudo codeP = self for i in {0, ..., dimension1}: P_left = P intersection left_guess count_left = P_left.count() if count_left < index: P = P_left else: P = P intersection right_guess index = count_left
Sure, that would probably be even faster since the polyhedron will become smaller as the process goes on. I've just finished a new version based on your pseudocode which repeatedly defines P as the intersection of P with a halfspace. Although of course it includes a while loop to keep doing the subdivisions on the x_i coordinate until it is uniquely determined.
comment:18 Changed 2 years ago by
 Commit changed from 99683415186c9890bc065fe0883f383fd2bd0019 to 5aea167963a784f90cdd6a576f7f0e446b288666
comment:19 Changed 2 years ago by
With the current code, you might construct and intersect polyhedra on different base ring. Instead of Polyhedron(XXX)
it would be better to use
S = self.parent() S(ieqs=[[guess1] + [0] * i + [1] + [0] * (D  i  1)])
comment:20 Changed 2 years ago by
 Commit changed from 5aea167963a784f90cdd6a576f7f0e446b288666 to ff28e93eda39a2ff8a83dc8cf68af3d36f5d7976
Branch pushed to git repo; I updated commit sha1. New commits:
ff28e93  Construct halfspaces using self.parent() to stay in the same ring.

comment:21 Changed 2 years ago by
Apparently, there is a problem with the last commit
sage t long src/sage/geometry/polyhedron/base.py ********************************************************************** File "src/sage/geometry/polyhedron/base.py", line 5332, in sage.geometry.polyhedron.base.Polyhedron_base.get_integral_point Failed example: P.get_integral_point(1) Expected: (0, 0) Got: (0, 1) ********************************************************************** File "src/sage/geometry/polyhedron/base.py", line 5340, in sage.geometry.polyhedron.base.Polyhedron_base.get_integral_point Failed example: [Q.get_integral_point(i) for i in range(Q.integral_points_count())] == sorted(Q.integral_points()) Expected: True Got: False ********************************************************************** File "src/sage/geometry/polyhedron/base.py", line 5408, in sage.geometry.polyhedron.base.Polyhedron_base.random_integral_point Failed example: P.random_integral_point() in P.integral_points() Expected: True Got: False **********************************************************************
The reason is that the syntax is different (I was wrong in comment:19)
sage: P = Polyhedron(ieqs=[[0,1,0,0]]) sage: P A 3dimensional polyhedron in QQ^3 defined as the convex hull of 1 vertex, 1 ray, 2 lines sage: PP Polyhedra in QQ^3 sage: PP(ieqs=[[0,1,0,0]]) # this is wrong A 0dimensional polyhedron in QQ^3 defined as the convex hull of 1 vertex sage: Vrep = None sage: Hrep = [[[0,1,0,0]], []] sage: PP(Vrep, Hrep) # this is right A 3dimensional polyhedron in QQ^3 defined as the convex hull of 1 vertex, 1 ray, 2 lines
comment:22 Changed 2 years ago by
 Commit changed from ff28e93eda39a2ff8a83dc8cf68af3d36f5d7976 to 8d3f64162697d0580f8b7657ace4ee703a27bd03
comment:23 Changed 2 years ago by
In the description, I would rather use "Return the ``index``th integral point."
as sorted(self.integral_points())[index]
is not very readable. Though you can use the code in the paragraph description below, and more importantly (as you already did) in the examples.
On the other hand, the OUTPUT section is intended to describe the output type (and is not mandatory). As your function is simple enough, I would just remove it. You can merge the paragraph describing the method and the one from the OUTPUT.
After that, I think the branch is good enough for a first version!
comment:24 Changed 2 years ago by
 Commit changed from 8d3f64162697d0580f8b7657ace4ee703a27bd03 to ed313737d7579954968b5feee89ada4901cdd5e3
Branch pushed to git repo; I updated commit sha1. New commits:
ed31373  Reworked docstring.

comment:25 Changed 2 years ago by
 Commit changed from ed313737d7579954968b5feee89ada4901cdd5e3 to f8350785ffca7f2078cc1fd6537c77288eb05eb0
Branch pushed to git repo; I updated commit sha1. New commits:
f835078  Added code fragment markers to docstring.

comment:26 Changed 2 years ago by
 Commit changed from f8350785ffca7f2078cc1fd6537c77288eb05eb0 to f7e095abe7c5b5e56a5e4b18df991b37a841e68a
Branch pushed to git repo; I updated commit sha1. New commits:
f7e095a  Actually passing on kwds.

comment:27 Changed 2 years ago by
Can you doctest the fact that **kwds
are passed (for example you can pass a wrong keyword)?
comment:28 Changed 2 years ago by
 Commit changed from f7e095abe7c5b5e56a5e4b18df991b37a841e68a to f585b5ec67fd3048ab35c7b5cd800306563ebeed
Branch pushed to git repo; I updated commit sha1. New commits:
f585b5e  Using **kwds in get_integral_point and doctesting these.

comment:29 Changed 2 years ago by
typo in the doc: it bisects the the upper
comment:30 Changed 2 years ago by
 Commit changed from f585b5ec67fd3048ab35c7b5cd800306563ebeed to ed9bdefe2b193843cba5fe794b0dfd1e3cdcc2b9
Branch pushed to git repo; I updated commit sha1. New commits:
ed9bdef  Docstring typo.

comment:32 Changed 2 years ago by
 Status changed from positive_review to needs_work
sage t long warnlong 78.8 src/sage/geometry/polyhedron/base.py ********************************************************************** File "src/sage/geometry/polyhedron/base.py", line 5369, in sage.geometry.polyhedron.base.Polyhedron_base.get_integral_point Failed example: Q.get_integral_point(0, explicit_enumeration_threshold=0, triangulation='cddlib') Exception raised: Traceback (most recent call last): File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 517, in _run self.compile_and_execute(example, compiler, test.globs) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 920, in compile_and_execute exec(compiled, globs) File "<doctest sage.geometry.polyhedron.base.Polyhedron_base.get_integral_point[7]>", line 1, in <module> Q.get_integral_point(Integer(0), explicit_enumeration_threshold=Integer(0), triangulation='cddlib') File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base.py", line 5386, in get_integral_point if not 0 <= index < self.integral_points_count(**kwds): File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base_QQ.py", line 218, in integral_points_count **kwds) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/interfaces/latte.py", line 150, in count raise PackageNotFoundError('latte_int') PackageNotFoundError: the package 'latte_int' was not found. You can install it by running 'sage i latte_int' in a shell ********************************************************************** File "src/sage/geometry/polyhedron/base.py", line 5371, in sage.geometry.polyhedron.base.Polyhedron_base.get_integral_point Failed example: Q.get_integral_point(0, explicit_enumeration_threshold=0, triangulation='cddlib', foo=True) Expected: Traceback (most recent call last): ... RuntimeError: ... Got: <BLANKLINE> Traceback (most recent call last): File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 517, in _run self.compile_and_execute(example, compiler, test.globs) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 920, in compile_and_execute exec(compiled, globs) File "<doctest sage.geometry.polyhedron.base.Polyhedron_base.get_integral_point[8]>", line 1, in <module> Q.get_integral_point(Integer(0), explicit_enumeration_threshold=Integer(0), triangulation='cddlib', foo=True) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base.py", line 5386, in get_integral_point if not 0 <= index < self.integral_points_count(**kwds): File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base_QQ.py", line 218, in integral_points_count **kwds) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/interfaces/latte.py", line 150, in count raise PackageNotFoundError('latte_int') PackageNotFoundError: the package 'latte_int' was not found. You can install it by running 'sage i latte_int' in a shell ********************************************************************** File "src/sage/geometry/polyhedron/base.py", line 5443, in sage.geometry.polyhedron.base.Polyhedron_base.random_integral_point Failed example: P.random_integral_point(explicit_enumeration_threshold=0, triangulation='cddlib') # random Exception raised: Traceback (most recent call last): File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 517, in _run self.compile_and_execute(example, compiler, test.globs) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 920, in compile_and_execute exec(compiled, globs) File "<doctest sage.geometry.polyhedron.base.Polyhedron_base.random_integral_point[3]>", line 1, in <module> P.random_integral_point(explicit_enumeration_threshold=Integer(0), triangulation='cddlib') # random File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base.py", line 5470, in random_integral_point return self.get_integral_point(current_randstate().python_random().randint(0, count1), **kwds) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base.py", line 5386, in get_integral_point if not 0 <= index < self.integral_points_count(**kwds): File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base_QQ.py", line 218, in integral_points_count **kwds) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/interfaces/latte.py", line 150, in count raise PackageNotFoundError('latte_int') PackageNotFoundError: the package 'latte_int' was not found. You can install it by running 'sage i latte_int' in a shell ********************************************************************** File "src/sage/geometry/polyhedron/base.py", line 5445, in sage.geometry.polyhedron.base.Polyhedron_base.random_integral_point Failed example: P.random_integral_point(explicit_enumeration_threshold=0, triangulation='cddlib', foo=7) Expected: Traceback (most recent call last): ... RuntimeError: ... Got: <BLANKLINE> Traceback (most recent call last): File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 517, in _run self.compile_and_execute(example, compiler, test.globs) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/doctest/forker.py", line 920, in compile_and_execute exec(compiled, globs) File "<doctest sage.geometry.polyhedron.base.Polyhedron_base.random_integral_point[4]>", line 1, in <module> P.random_integral_point(explicit_enumeration_threshold=Integer(0), triangulation='cddlib', foo=Integer(7)) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base.py", line 5470, in random_integral_point return self.get_integral_point(current_randstate().python_random().randint(0, count1), **kwds) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base.py", line 5386, in get_integral_point if not 0 <= index < self.integral_points_count(**kwds): File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/geometry/polyhedron/base_QQ.py", line 218, in integral_points_count **kwds) File "/mnt/disk/home/release/Sage/local/lib/python2.7/sitepackages/sage/interfaces/latte.py", line 150, in count raise PackageNotFoundError('latte_int') PackageNotFoundError: the package 'latte_int' was not found. You can install it by running 'sage i latte_int' in a shell ********************************************************************** 2 items had failures: 2 of 12 in sage.geometry.polyhedron.base.Polyhedron_base.get_integral_point 2 of 10 in sage.geometry.polyhedron.base.Polyhedron_base.random_integral_point [998 tests, 4 failures, 55.73 s]
comment:33 Changed 2 years ago by
 Commit changed from ed9bdefe2b193843cba5fe794b0dfd1e3cdcc2b9 to 73fbd0571a629071bceb1300d699f0b826e983c1
Branch pushed to git repo; I updated commit sha1. New commits:
73fbd05  Added # optional  latte_int markers to doctests that require latteint.

comment:34 Changed 23 months ago by
I think this updated version should have the doctests that require latte_int flagged correctly. However to test this I had to go into my install of sage and manually delete /build/pkgs/latte_int and ./local/bin/count. Is there an easier way to (temporarily) uninstall an optional sage package?
comment:35 Changed 23 months ago by
 Reviewers changed from Vincent Delecroix to Vincent Delecroix, Travis Scrimshaw
 Status changed from needs_work to positive_review
Right now there is not, but Erik Bray (embray) is currently working on developing this. IIRC he has a post on sagedevel within the past month or so about the progress of this.
comment:36 Changed 23 months ago by
 Branch changed from u/mcbell/polyhedron_get_integral_point to 73fbd0571a629071bceb1300d699f0b826e983c1
 Resolution set to fixed
 Status changed from positive_review to closed
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.