Polyhedron.integral_points_count() gives inconsistent answers on slices

[Mark Bell]

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

[Simon King]

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]

I opened trac ticket #24152 for the bug in converting the polyhedron
into polymake.

[Mark Bell]

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]

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.

[Matthias Koeppe]

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]

Looks like to_linear_program() handles equations wrong!

 

[Matthias Koeppe]

Ticket at https://trac.sagemath.org/ticket/24154 

[Matthias Koeppe]

Ticket needs review. 

[Mark Bell]

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.

[Matthias Koeppe]

Can you show the output of the crash on CoCalc?

[William Stein]

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]

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

Available from http://www.math.ucdavis.edu/~latte/

 

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]

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]

I can't reproduce this locally, but you could see if passing 

--triangulation=cddlib

to LattE helps. 

[Mark Bell]

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]

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]

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.