Polyhedron.integral_points_count() gives inconsistent answers on slices

122 views
Skip to first unread message

Mark Bell

unread,
Nov 3, 2017, 9:52:52 AM11/3/17
to sage-devel
I want to count the number of integral points inside of a bounded polyhedron. The polyhedron is defined by the following system of equations:

> eqns = [
[-100, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, -1, 0, 0, 0, 0, -1, 1, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, -1, 0, 0, 1, 1, -1, 0, 0],
[0, 0, 0, 0, -1, 0, -1, 1, 0, 1, 0, 0, 0],
[0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, -1],
[0, 0, -1, -1, 0, 1, 0, 0, 0, 0, 0, 1, 0],
[0, -1, 1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0]
]

and the following system of inequalities:

> ieqs = [
[-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[-1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
[-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
[-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
]

And so I did:
> Polyhedron(eqns=eqns, ieqs=ieqs).integral_points_count()
1260

However, if I look at all of the slices obtained by fixing the first variable x_0 (which must be 1 <= x_0 <= 100 by the first equation and inequality) then I get the following numbers of integral points in each slice:

> [Polyhedron(eqns=eqns+[[-k, 1, 0,0,0,0,0,0,0,0,0,0,0]], ieqs=ieqs).integral_points_count() for k in range(1, 101)]
[84, 84, 84, ...., 0, 0, 0]

Which is fewer points in total:

> sum(_)
1197

Can anyone reproduce this? Is this a bug in Sage or LattE?

I'm using Sage 8.0 with LattE installed via "sage -i latte_int" and have attached the full example code.

Thanks,
Mark
scratch.py

Simon King

unread,
Nov 3, 2017, 10:37:13 AM11/3/17
to sage-...@googlegroups.com
Hi!

On 2017-11-03, 'Mark Bell' via sage-devel <sage-...@googlegroups.com> wrote:
> And so I did:
>> Polyhedron(eqns=eqns, ieqs=ieqs).integral_points_count()
> 1260
>
> However, if I look at all of the slices obtained by fixing the first
> variable x_0 (which must be 1 <= x_0 <= 100 by the first equation and
> inequality) then I get the following numbers of integral points in each
> slice:
>
>> [Polyhedron(eqns=eqns+[[-k, 1, 0,0,0,0,0,0,0,0,0,0,0]],
> ieqs=ieqs).integral_points_count() for k in range(1, 101)]
> [84, 84, 84, ...., 0, 0, 0]
>
> Which is fewer points in total:
>
>> sum(_)
> 1197

Sorry that I am not directly answering your question. What does polymake
have to say about that polyhedron?

First of all, there seems to be some bug in the polymake interface
(that I authored, so: Sorry...):
sage: P = Polyhedron(eqns = eqns, ieqs=ieqs)
sage: PP = polymake(P)

The conversion of P into a polymake polytope fails first (that's the bug),
but when one tries again, there is no complaint.

sage: PP.VERTICES
1 1 1 46/3 1 49/3 46/3 0 49/3 49/3 49/3 0 1
1 45/2 45/2 1 1 47/2 1 0 47/2 2 2 0 1
1 45/2 1 1 1 2 1 0 47/2 2 47/2 0 45/2
1 1 1 1 45/2 2 1 0 2 47/2 47/2 0 45/2
1 1 45/2 1 45/2 47/2 1 0 2 47/2 2 0 1
sage: P.vertices_list() # Sage agrees with Polymake
[[1, 1, 46/3, 1, 49/3, 46/3, 0, 49/3, 49/3, 49/3, 0, 1],
[45/2, 45/2, 1, 1, 47/2, 1, 0, 47/2, 2, 2, 0, 1],
[45/2, 1, 1, 1, 2, 1, 0, 47/2, 2, 47/2, 0, 45/2],
[1, 1, 1, 45/2, 2, 1, 0, 2, 47/2, 47/2, 0, 45/2],
[1, 45/2, 1, 45/2, 47/2, 1, 0, 2, 47/2, 2, 0, 1]]
sage: PP.N_LATTICE_POINTS # Sage agrees with Polymake again.
1260

Best regards,
Simon

Simon King

unread,
Nov 3, 2017, 10:45:21 AM11/3/17
to sage-...@googlegroups.com
I opened trac ticket #24152 for the bug in converting the polyhedron
into polymake.

Mark Bell

unread,
Nov 3, 2017, 11:19:35 AM11/3/17
to sage-devel
Great suggestion. Could you ask polymake to compute what the number of lattice points in each x_0 = k (for k = 1, ..., 100) slice? Does it match:
    [84, 84, 84, 81, 81, 81, 75, 75, 75, 66, 66, 66, 54, 54, 54, 39, 39, 39, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Mark Bell

unread,
Nov 3, 2017, 4:41:33 PM11/3/17
to sage-devel
So when k = 19, sage.interfaces.latte.count and P.integral_points_count() give different answers. This seems to be the cause of the difference. Here is some more code showing the difference.
scratch.py

Matthias Koeppe

unread,
Nov 3, 2017, 8:06:13 PM11/3/17
to sage-devel
On Friday, November 3, 2017 at 1:41:33 PM UTC-7, Mark Bell wrote:
So when k = 19, sage.interfaces.latte.count and P.integral_points_count() give different answers. This seems to be the cause of the difference. Here is some more code showing the difference.

The bug appears to be in the preprocessing code of integral_points_count().
If you use preprocess=False, the number of lattice points is correct (19).

 

Matthias Koeppe

unread,
Nov 3, 2017, 8:22:16 PM11/3/17
to sage-devel
Looks like to_linear_program() handles equations wrong!
 

Matthias Koeppe

unread,
Nov 3, 2017, 8:44:43 PM11/3/17
to sage-devel

Matthias Koeppe

unread,
Nov 3, 2017, 11:07:48 PM11/3/17
to sage-devel
Ticket needs review. 

Mark Bell

unread,
Nov 4, 2017, 6:46:01 PM11/4/17
to sage-devel
Thanks! 

I'm not sure if this is the right place to add this but the attached example also causes LattE to crash, but I only see this behaviour when I run this code on CoCalc. When I run the script on my local machine with latte_int installed it completes and prints out 19958. Surely this isn't being caused by the same bug.
error.py

Matthias Koeppe

unread,
Nov 4, 2017, 9:27:48 PM11/4/17
to sage-devel
Can you show the output of the crash on CoCalc?

William Stein

unread,
Nov 4, 2017, 10:05:31 PM11/4/17
to sage-devel
On Sat, Nov 4, 2017 at 6:27 PM, Matthias Koeppe
<matthia...@gmail.com> wrote:
> Can you show the output of the crash on CoCalc?

Also, the crash could just be running out of memory...

>
>
> On Saturday, November 4, 2017 at 3:46:01 PM UTC-7, Mark Bell wrote:
>>
>> Thanks!
>>
>> I'm not sure if this is the right place to add this but the attached
>> example also causes LattE to crash, but I only see this behaviour when I run
>> this code on CoCalc. When I run the script on my local machine with
>> latte_int installed it completes and prints out 19958. Surely this isn't
>> being caused by the same bug.
>
> --
> You received this message because you are subscribed to the Google Groups
> "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sage-devel+...@googlegroups.com.
> To post to this group, send email to sage-...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sage-devel.
> For more options, visit https://groups.google.com/d/optout.



--
William (http://wstein.org)

Mark Bell

unread,
Nov 5, 2017, 4:46:30 AM11/5/17
to sage-devel
Here is the full CoCalc output when I run "sage error.py" in a terminal, either via the web interface or via sshing into the machine:

Traceback (most recent call last):
  File "error.py", line 42, in <module>
    print(P.integral_points_count(preprocess=False))
  File "/ext/sage/sage-8.0/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base_QQ.py", line 218, in integral_points_count
    **kwds)
  File "/ext/sage/sage-8.0/local/lib/python2.7/site-packages/sage/interfaces/latte.py", line 159, in count
    raise RuntimeError("LattE integrale program failed (exit code {})".format(ret_code) + err.strip())
RuntimeError: LattE integrale program failed (exit code -6):
This is LattE integrale 1.7.3
 
Invocation: count '--redundancy-check=none' --cdd /dev/stdin
Warning: Not performing check for empty polytope, because it is unimplemented for the CDD-style input format.
size = 22 x 19
Number Type = rational
Ax <= b, given as (b|-A):
=========================
[-1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[-1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[-1 1 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[-1 0 1 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[-1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[-1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0]
[97 -2 -8 -1 -7 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[-1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[97 -2 -5 -1 -4 0 -3 0 0 0 0 0 0 0 0 0 0 0 0]

Ax = b, given as (b|-A):
========================
[0 1 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1]
[0 0 1 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 -1 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0]
[0 0 1 -1 2 0 0 0 0 0 0 0 0 0 0 -1 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 1 0 0 0 0 0 0 -1 0 0 0 0 0]
[-100 2 8 1 7 0 0 0 0 0 0 0 3 0 0 0 0 0 0]
[0 0 1 0 1 0 0 0 0 0 0 -1 0 0 0 0 0 0 0]
[0 0 1 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0 0]
[0 0 1 1 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0]
[-100 2 5 1 4 0 3 0 3 0 0 0 0 0 0 0 0 0 0]
[-100 2 5 1 4 0 3 3 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0]

Computing hermitean normal form.
sol:
[0 0 0 0 0 0 -100 0 0 -100 0 0 0 0 0 0 0 0]
Particular solution:
Basis:
New inequalities:
[99 -3 -3 -5 -5 -1 0]
[-101 3 3 5 6 2 0]
[-1 0 0 0 0 1 0]
[-1 0 0 0 1 0 0]
[-1 0 0 1 0 0 0]
[-1 1 0 0 0 0 0]
[-3 3 3 -3 -3 0 0]
[99 -3 -3 -4 -5 -2 0]
[-3 0 3 0 0 0 0]
Time for reading and preprocessing: 0 sec
The polytope has 12 vertices.
Time for computing vertices and supporting cones: 0 sec
Time for dualizing general cones: 0 sec
Dualizing all cones...All cones are now dualized.
Time for dualizing general cones: 0 sec
decomposeCones_Single: Decomposing all cones. (Memory Save on)
12 cones total to be done!decomposeCones_Single: degree = 1
Number of cones: 12
Triangulating cone... done.
count: barvinok/barvinok.cpp:572: int barvinokDecomposition_Single(listCone*, Single_Cone_Parameters*): Assertion `num_rays == Parameters->Number_of_Variables'
failed.

Mark Bell

unread,
Nov 5, 2017, 5:33:01 AM11/5/17
to sage-devel
Thanks for the suggestion. That might be the case but if I run "/usr/bin/time -v sage error.py" then I get 

Command exited with non-zero status 1
Command being timed: "sage error.py"
User time (seconds): 2.22
System time (seconds): 0.67
Percent of CPU this job got: 86%
Elapsed (wall clock) time (h:mm:ss or m:ss): 0:03.35
Average shared text size (kbytes): 0
Average unshared data size (kbytes): 0
Average stack size (kbytes): 0
Average total size (kbytes): 0
Maximum resident set size (kbytes): 200116
Average resident set size (kbytes): 0
Major (requiring I/O) page faults: 0
Minor (reclaiming a frame) page faults: 81745
Voluntary context switches: 7875
Involuntary context switches: 191
Swaps: 0
File system inputs: 22
File system outputs: 4
Socket messages sent: 0
Socket messages received: 0
Signals delivered: 0
Page size (bytes): 4096
Exit status: 1

So it appears to be using ~200 MB of the 1000 MB shared RAM memory limit. Is there a better way to determine the resources being used / if I am hitting this limit?

Matthias Koeppe

unread,
Nov 5, 2017, 5:01:22 PM11/5/17
to sage-devel
I can't reproduce this locally, but you could see if passing 

--triangulation=cddlib

to LattE helps. 

Mark Bell

unread,
Nov 5, 2017, 6:02:01 PM11/5/17
to sage-devel
Excellent! Adding the option prevents the crash. So on CoCalc (and locally) I can do:
   print(P.integral_points_count(preprocess=False, triangulation='cddlib'))
And it prints out 19958.

Without this extra option this example crashes on CoCalc but runs fine locally. What is this this option doing?

Dima Pasechnik

unread,
Nov 6, 2017, 2:25:58 AM11/6/17
to sage-devel
lattE needs to compute a triangulation of the polytope; by default it uses 4ti2, but can instead use cddlib. Hence this option. See
https://www.math.ucdavis.edu/~latte/software/packages/latte_current/manual_v1.7.2.pdf

Mark Bell

unread,
Nov 6, 2017, 7:12:50 AM11/6/17
to sage-devel
Ah, I see. Is it possible that by doing "sage -i latte_int" I have installed a different version of 4ti2 than the one that CoCalc uses? I'm still not sure why (without this option) I see different behaviours on my machine and on CoCalc. If it matters I have Sage 8.0 running on Mac OS Sierra 10.12.6.
Reply all
Reply to author
Forward
0 new messages