Opened 2 years ago

Closed 23 months ago

#24451 closed enhancement (fixed)

Polyhedron.get_integral_point

Reported by: mcbell Owned by:
Priority: major Milestone: sage-8.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 mcbell)

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 mcbell

  • Branch set to u/mcbell/polyhedron_get_integral_point

comment:2 Changed 2 years ago by mcbell

  • Authors set to Mark Bell
  • 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 vdelecroix

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: Changed 2 years ago by mcbell

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

comment:5 in reply to: ↑ 4 ; follow-up: Changed 2 years ago by vdelecroix

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 git

  • Commit changed from d3b9b396c1dc9d6ca8c9478b01035454ebbe1f74 to 82a64d2f6eb533c16413c681a8d89c130283fa1c

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

e020df5Initial version of get_integral_point method.
ca57e15Added docstring and more comments.
f771506Fixed docstring whitespace.
626e399Fixed line widths and fencepost error in upper bound.
c46460dAdded sorting comments to docstring and random_integral_point method.
82a64d2Added missing import.

comment:7 in reply to: ↑ 5 Changed 2 years ago by mcbell

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 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 2 years ago by mcbell

  • Description modified (diff)

comment:9 Changed 2 years ago by vdelecroix

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

  • Commit changed from 82a64d2f6eb533c16413c681a8d89c130283fa1c to 11aabdd974fba81accf9966584a159fba8621c2e

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

11aabddAdded check to make sure there are integral points before sampling and doctests for the different errors.

comment:11 Changed 2 years ago by vdelecroix

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 git

  • Commit changed from 11aabdd974fba81accf9966584a159fba8621c2e to 4e76f7d28eb839bac153ea682e50707c5bdb4754

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

4e76f7dUsing current_randstate() instead of random.

comment:13 Changed 2 years ago by git

  • Commit changed from 4e76f7d28eb839bac153ea682e50707c5bdb4754 to 99683415186c9890bc065fe0883f383fd2bd0019

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

9968341Fixed upper and lower bounds to be integers.

comment:14 Changed 2 years ago by vdelecroix

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: Changed 2 years ago by vdelecroix

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 2 years ago by vdelecroix

  • Reviewers set to Vincent Delecroix

comment:17 in reply to: ↑ 15 Changed 2 years ago by mcbell

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

  • Commit changed from 99683415186c9890bc065fe0883f383fd2bd0019 to 5aea167963a784f90cdd6a576f7f0e446b288666

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

911bfdcReturned point is now immutable.
73af991More explict error cases.
5aea167Rewitten to repeatedly intersect the polyhedron with halfspaces.

comment:19 Changed 2 years ago by vdelecroix

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

  • Commit changed from 5aea167963a784f90cdd6a576f7f0e446b288666 to ff28e93eda39a2ff8a83dc8cf68af3d36f5d7976

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

ff28e93Construct halfspaces using self.parent() to stay in the same ring.

comment:21 Changed 2 years ago by vdelecroix

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 2 years ago by vdelecroix (previous) (diff)

comment:22 Changed 2 years ago by git

  • Commit changed from ff28e93eda39a2ff8a83dc8cf68af3d36f5d7976 to 8d3f64162697d0580f8b7657ace4ee703a27bd03

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

7a0b647Merge branch 'develop' into t/24451/polyhedron_get_integral_point
8d3f641Using correct syntax for creating halfspaces via parent.

comment:23 Changed 2 years ago by vdelecroix

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 git

  • Commit changed from 8d3f64162697d0580f8b7657ace4ee703a27bd03 to ed313737d7579954968b5feee89ada4901cdd5e3

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

ed31373Reworked docstring.

comment:25 Changed 2 years ago by git

  • Commit changed from ed313737d7579954968b5feee89ada4901cdd5e3 to f8350785ffca7f2078cc1fd6537c77288eb05eb0

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

f835078Added code fragment markers to docstring.

comment:26 Changed 2 years ago by git

  • Commit changed from f8350785ffca7f2078cc1fd6537c77288eb05eb0 to f7e095abe7c5b5e56a5e4b18df991b37a841e68a

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

f7e095aActually passing on kwds.

comment:27 Changed 2 years ago by vdelecroix

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

comment:28 Changed 2 years ago by git

  • Commit changed from f7e095abe7c5b5e56a5e4b18df991b37a841e68a to f585b5ec67fd3048ab35c7b5cd800306563ebeed

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

f585b5eUsing **kwds in get_integral_point and doctesting these.

comment:29 Changed 2 years ago by vdelecroix

typo in the doc: it bisects the the upper

comment:30 Changed 2 years ago by git

  • Commit changed from f585b5ec67fd3048ab35c7b5cd800306563ebeed to ed9bdefe2b193843cba5fe794b0dfd1e3cdcc2b9

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

ed9bdefDocstring typo.

comment:31 Changed 2 years ago by vdelecroix

  • Status changed from needs_review to positive_review

let it go

comment:32 Changed 2 years ago by vbraun

  • Status changed from positive_review to 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 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 git

  • Commit changed from ed9bdefe2b193843cba5fe794b0dfd1e3cdcc2b9 to 73fbd0571a629071bceb1300d699f0b826e983c1

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

73fbd05Added # optional - latte_int markers to doctests that require latte-int.

comment:34 Changed 23 months ago by mcbell

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 tscrim

  • 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 sage-devel within the past month or so about the progress of this.

comment:36 Changed 23 months ago by vbraun

  • Branch changed from u/mcbell/polyhedron_get_integral_point to 73fbd0571a629071bceb1300d699f0b826e983c1
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.