Opened 5 years ago

Closed 5 years ago

# Polyhedron.get_integral_point

Reported by: Owned by: Mark Bell major sage-8.2 geometry Polyhedron, integral_points Vincent Delecroix, Thierry Monteil, Matthias Köppe Mark Bell Vincent Delecroix, Travis Scrimshaw N/A 73fbd05 73fbd0571a629071bceb1300d699f0b826e983c1

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)

### comment:9 Changed 5 years ago by Vincent Delecroix

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: error-message
```

(you can have a look around in the code)

### comment:10 Changed 5 years ago by git

Commit: 82a64d2f6eb533c16413c681a8d89c130283fa1c → 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 5 years ago by Vincent Delecroix

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 5 years ago by git

Commit: 11aabdd974fba81accf9966584a159fba8621c2e → 4e76f7d28eb839bac153ea682e50707c5bdb4754

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

 ​4e76f7d `Using current_randstate() instead of random.`

### comment:13 Changed 5 years ago by git

Commit: 4e76f7d28eb839bac153ea682e50707c5bdb4754 → 99683415186c9890bc065fe0883f383fd2bd0019

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

 ​9968341 `Fixed upper and lower bounds to be integers.`

### comment:14 Changed 5 years ago by Vincent Delecroix

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 follow-up:  17 Changed 5 years ago by Vincent Delecroix

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, ..., dimension-1}:
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 5 years ago by Vincent Delecroix

Reviewers: → Vincent Delecroix

### comment:17 in reply to:  15 Changed 5 years ago by Mark Bell

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 with `extra_eqns` along the loop. In pseudo code

```P = self
for i in {0, ..., dimension-1}:
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 pseudo-code 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 5 years ago by git

Commit: 99683415186c9890bc065fe0883f383fd2bd0019 → 5aea167963a784f90cdd6a576f7f0e446b288666

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

 ​911bfdc `Returned point is now immutable.` ​73af991 `More explict error cases.` ​5aea167 `Rewitten to repeatedly intersect the polyhedron with halfspaces.`

### comment:19 Changed 5 years ago by Vincent Delecroix

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=[[guess-1] + [0] * i + [-1] + [0] * (D - i - 1)])
```

### comment:20 Changed 5 years ago by git

Commit: 5aea167963a784f90cdd6a576f7f0e446b288666 → 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 5 years ago by Vincent Delecroix

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 3-dimensional 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 0-dimensional 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 3-dimensional polyhedron in QQ^3 defined as the convex hull of 1 vertex, 1 ray, 2 lines
```
Last edited 5 years ago by Vincent Delecroix (previous) (diff)

### comment:22 Changed 5 years ago by git

Commit: ff28e93eda39a2ff8a83dc8cf68af3d36f5d7976 → 8d3f64162697d0580f8b7657ace4ee703a27bd03

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

 ​7a0b647 `Merge branch 'develop' into t/24451/polyhedron_get_integral_point` ​8d3f641 `Using correct syntax for creating halfspaces via parent.`

### comment:23 Changed 5 years ago by Vincent Delecroix

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 5 years ago by git

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

 ​ed31373 `Reworked docstring.`

### comment:25 Changed 5 years ago by git

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

 ​f835078 `Added code fragment markers to docstring.`

### comment:26 Changed 5 years ago by git

Commit: f8350785ffca7f2078cc1fd6537c77288eb05eb0 → f7e095abe7c5b5e56a5e4b18df991b37a841e68a

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

 ​f7e095a `Actually passing on kwds.`

### comment:27 Changed 5 years ago by Vincent Delecroix

Can you doctest the fact that `**kwds` are passed (for example you can pass a wrong keyword)?

### comment:28 Changed 5 years ago by git

Commit: f7e095abe7c5b5e56a5e4b18df991b37a841e68a → 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 5 years ago by Vincent Delecroix

typo in the doc: `it bisects the the upper`

### comment:30 Changed 5 years ago by git

Commit: f585b5ec67fd3048ab35c7b5cd800306563ebeed → ed9bdefe2b193843cba5fe794b0dfd1e3cdcc2b9

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

 ​ed9bdef `Docstring typo.`

### comment:31 Changed 5 years ago by Vincent Delecroix

Status: needs_review → positive_review

let it go

### comment:32 Changed 5 years ago by Volker Braun

Status: positive_review → needs_work
```sage -t --long --warn-long 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/site-packages/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/site-packages/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/site-packages/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/site-packages/sage/geometry/polyhedron/base_QQ.py", line 218, in integral_points_count
**kwds)
File "/mnt/disk/home/release/Sage/local/lib/python2.7/site-packages/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/site-packages/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/site-packages/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/site-packages/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/site-packages/sage/geometry/polyhedron/base_QQ.py", line 218, in integral_points_count
**kwds)
File "/mnt/disk/home/release/Sage/local/lib/python2.7/site-packages/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/site-packages/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/site-packages/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/site-packages/sage/geometry/polyhedron/base.py", line 5470, in random_integral_point
return self.get_integral_point(current_randstate().python_random().randint(0, count-1), **kwds)
File "/mnt/disk/home/release/Sage/local/lib/python2.7/site-packages/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/site-packages/sage/geometry/polyhedron/base_QQ.py", line 218, in integral_points_count
**kwds)
File "/mnt/disk/home/release/Sage/local/lib/python2.7/site-packages/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/site-packages/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/site-packages/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/site-packages/sage/geometry/polyhedron/base.py", line 5470, in random_integral_point
return self.get_integral_point(current_randstate().python_random().randint(0, count-1), **kwds)
File "/mnt/disk/home/release/Sage/local/lib/python2.7/site-packages/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/site-packages/sage/geometry/polyhedron/base_QQ.py", line 218, in integral_points_count
**kwds)
File "/mnt/disk/home/release/Sage/local/lib/python2.7/site-packages/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 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 5 years ago by git

Commit: ed9bdefe2b193843cba5fe794b0dfd1e3cdcc2b9 → 73fbd0571a629071bceb1300d699f0b826e983c1

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

 ​73fbd05 `Added # optional - latte_int markers to doctests that require latte-int.`

### comment:34 Changed 5 years ago by Mark Bell

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 5 years ago by Travis Scrimshaw

Reviewers: Vincent Delecroix → Vincent Delecroix, Travis Scrimshaw needs_work → positive_review

Right now there is not, but Erik Bray (embray) is currently working on developing this. IIRC he has a post on sage-devel within the past month or so about the progress of this.

### comment:36 Changed 5 years ago by Volker Braun

Branch: u/mcbell/polyhedron_get_integral_point → 73fbd0571a629071bceb1300d699f0b826e983c1 → fixed positive_review → closed
Note: See TracTickets for help on using tickets.