Skip to content

Commit

Permalink
Trac #22497: generic latte_int interface to count
Browse files Browse the repository at this point in the history
We provide a file `sage/interfaces/latte.py` that provides a uniform
function `count` to query the LattE program `count` (`integrate` will be
dealt with in another ticket). It is used in this ticket to refactor
`ehrhart_polynomial` and `integral_points_count`.

Prerequisite for #18232.

URL: https://trac.sagemath.org/22497
Reported by: vdelecroix
Ticket author(s): Vincent Delecroix
Reviewer(s): Matthias Koeppe
  • Loading branch information
Release Manager authored and vbraun committed Mar 4, 2017
2 parents 1d34baf + d5ff154 commit 84abd08
Show file tree
Hide file tree
Showing 3 changed files with 219 additions and 88 deletions.
58 changes: 20 additions & 38 deletions src/sage/geometry/polyhedron/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -4483,7 +4483,7 @@ def bounding_box(self, integral=False):
box_min.append(min_coord)
return (tuple(box_min), tuple(box_max))

def integral_points_count(self, verbose=False):
def integral_points_count(self, verbose=False, use_Hrepresentation=False, **kwds):
r"""
Return the number of integral points in the polyhedron.
Expand All @@ -4494,6 +4494,13 @@ def integral_points_count(self, verbose=False):
- ``verbose`` (boolean; ``False`` by default) -- whether to display
verbose output.
- ``use_Hrepresentation`` - (boolean; ``False`` by default) -- whether
to send the H or V representation to LattE
.. SEEALSO::
:mod:`~sage.interfaces.latte` the interface to LattE interfaces
EXAMPLES::
sage: P = polytopes.cube()
Expand All @@ -4517,12 +4524,15 @@ def integral_points_count(self, verbose=False):
sage: Q.integral_points_count() # optional - latte_int
Traceback (most recent call last):
...
RuntimeError: LattE integrale failed (exit code 1) to execute...
...Parse error in CDD-style input file /dev/stdin
RuntimeError: LattE integrale program failed (exit code 1):
...
Invocation: count '--redundancy-check=none' --cdd /dev/stdin
...
Parse error in CDD-style input file /dev/stdin
sage: Q.integral_points_count(verbose=True) # optional - latte_int
Traceback (most recent call last):
...
RuntimeError: LattE integrale failed (exit code 1) to execute count --cdd /dev/stdin, see error message above
RuntimeError: LattE integrale program failed (exit code 1), see error message above
TESTS:
Expand All @@ -4538,40 +4548,12 @@ def integral_points_count(self, verbose=False):
if self.is_empty():
return 0

from subprocess import Popen, PIPE
from sage.misc.misc import SAGE_TMP
from sage.rings.integer import Integer

ine = self.cdd_Hrepresentation()
args = ['count', '--cdd', '/dev/stdin']

try:
# The cwd argument is needed because latte
# always produces diagnostic output files.
latte_proc = Popen(args,
stdin=PIPE, stdout=PIPE,
stderr=(None if verbose else PIPE),
cwd=str(SAGE_TMP))
except OSError:
from sage.misc.package import PackageNotFoundError
raise PackageNotFoundError('latte_int')

ans, err = latte_proc.communicate(ine)
ret_code = latte_proc.poll()
if ret_code:
if err is None:
err = ", see error message above"
else:
err = ":\n" + err
raise RuntimeError("LattE integrale failed (exit code {}) to execute {}".format(ret_code, ' '.join(args)) + err.strip())

try:
return Integer(ans.splitlines()[-1])
except IndexError:
# opening a file is slow (30e-6s), so we read the file
# numOfLatticePoints only in case of a IndexError above
with open(SAGE_TMP+'/numOfLatticePoints', 'r') as f:
return Integer(f.read())
from sage.interfaces.latte import count
return count(
self.cdd_Hrepresentation() if use_Hrepresentation else self.cdd_Vrepresentation(),
cdd=True,
verbose=verbose,
**kwds)

def integral_points(self, threshold=100000):
r"""
Expand Down
69 changes: 19 additions & 50 deletions src/sage/geometry/polyhedron/base_ZZ.py
Original file line number Diff line number Diff line change
Expand Up @@ -178,6 +178,10 @@ def ehrhart_polynomial(self, verbose=False, dual=None,
for lattice point enumeration (see
https://www.math.ucdavis.edu/~latte/).
.. SEEALSO::
:mod:`~sage.interfaces.latte` the interface to LattE integrale
EXAMPLES::
sage: P = Polyhedron(vertices=[(0,0,0),(3,3,3),(-3,2,1),(1,-1,-2)])
Expand Down Expand Up @@ -224,31 +228,32 @@ def ehrhart_polynomial(self, verbose=False, dual=None,
sage: p = P.ehrhart_polynomial(maxdet=5, verbose=True) # optional - latte_int
This is LattE integrale ...
...
Invocation: count --ehrhart-polynomial '--redundancy-check=none' '--maxdet=5' --cdd ...
Invocation: count --ehrhart-polynomial '--redundancy-check=none' --cdd '--maxdet=5' /dev/stdin
...
sage: p # optional - latte_int
1/2*t^2 + 3/2*t + 1
sage: p = P.ehrhart_polynomial(dual=True, verbose=True) # optional - latte_int
This is LattE integrale ...
...
Invocation: count --ehrhart-polynomial '--redundancy-check=none' --dual --cdd ...
Invocation: count --ehrhart-polynomial '--redundancy-check=none' --cdd --dual /dev/stdin
...
sage: p # optional - latte_int
1/2*t^2 + 3/2*t + 1
sage: p = P.ehrhart_polynomial(irrational_primal=True, verbose=True) # optional - latte_int
This is LattE integrale ...
...
Invocation: count --ehrhart-polynomial '--redundancy-check=none' --irrational-primal --cdd ...
Invocation: count --ehrhart-polynomial '--redundancy-check=none' --cdd --irrational-primal /dev/stdin
...
sage: p # optional - latte_int
1/2*t^2 + 3/2*t + 1
sage: p = P.ehrhart_polynomial(irrational_all_primal=True, verbose=True) # optional - latte_int
This is LattE integrale ...
...
Invocation: count --ehrhart-polynomial '--redundancy-check=none' --irrational-all-primal --cdd ...
Invocation: count --ehrhart-polynomial '--redundancy-check=none' --cdd --irrational-all-primal /dev/stdin
...
sage: p # optional - latte_int
1/2*t^2 + 3/2*t + 1
Expand All @@ -257,22 +262,17 @@ def ehrhart_polynomial(self, verbose=False, dual=None,
sage: P.ehrhart_polynomial(bim_bam_boum=19) # optional - latte_int
Traceback (most recent call last):
...
RuntimeError: LattE integrale failed with exit code 1 to execute...
RuntimeError: LattE integrale program failed (exit code 1):
...
Invocation: count --ehrhart-polynomial '--redundancy-check=none' --cdd '--bim-bam-boum=19' /dev/stdin
Unknown command/option --bim-bam-boum=19
"""
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
R = PolynomialRing(QQ, 't')
if self.is_empty():
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.rational_field import QQ
R = PolynomialRing(QQ, 't')
return R.zero()

from sage.misc.misc import SAGE_TMP
from subprocess import Popen, PIPE

ine = self.cdd_Hrepresentation()

args = ['count', '--ehrhart-polynomial']
if 'redundancy_check' not in kwds:
args.append('--redundancy-check=none')

# note: the options below are explicitely written in the function
# declaration in order to keep tab completion (see #18211).
kwds.update({
Expand All @@ -287,40 +287,9 @@ def ehrhart_polynomial(self, verbose=False, dual=None,
'triangulation' : triangulation,
'triangulation_max_height': triangulation_max_height})

for key,value in kwds.items():
if value is None or value is False:
continue

key = key.replace('_','-')
if value is True:
args.append('--{}'.format(key))
else:
args.append('--{}={}'.format(key, value))
args += ['--cdd', '/dev/stdin']

try:
# The cwd argument is needed because latte
# always produces diagnostic output files.
latte_proc = Popen(args,
stdin=PIPE, stdout=PIPE,
stderr=(None if verbose else PIPE),
cwd=str(SAGE_TMP))
except OSError:
from sage.misc.package import PackageNotFoundError
raise PackageNotFoundError('latte_int')

ans, err = latte_proc.communicate(ine)
ret_code = latte_proc.poll()
if ret_code:
if err is None:
err = ", see error message above"
else:
err = ":\n" + err
raise RuntimeError("LattE integrale failed with exit code {} to execute {}".format(ret_code, ' '.join(args)) + err.strip())

p = ans.splitlines()[-2]

return R(p)
from sage.interfaces.latte import count
ine = self.cdd_Hrepresentation()
return count(ine, cdd=True, ehrhart_polynomial=True, verbose=verbose, **kwds)

@cached_method
def polar(self):
Expand Down
180 changes: 180 additions & 0 deletions src/sage/interfaces/latte.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,180 @@
r"""
Generic interface to LattE integrale programs
"""
#*****************************************************************************
# Copyright (C) 2017 Vincent Delecroix <[email protected]>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
# http://www.gnu.org/licenses/
#*****************************************************************************

def count(arg, ehrhart_polynomial=False, multivariate_generating_function=False, raw_output=False, verbose=False, **kwds):
r"""
Call to the program count from LattE integrale
INPUT:
- ``arg`` -- a cdd or LattE description string
- ``ehrhart_polynomial``, ``multivariate_generating_function`` -- to
compute Ehrhart polynomial or multivariate generating function instead of
just counting points
- ``raw_output`` -- if ``True`` then return directly the output string from LattE
- For all other options of the count program, consult the LattE manual
OUTPUT:
Either a string (if ``raw_output`` if set to ``True``) or an integer (when
counting points), or a polynomial (if ``ehrhart_polynomial`` is set to
``True``) or a multivariate THING (if ``multivariate_generating_function``
is set to ``True``)
EXAMPLES::
sage: from sage.interfaces.latte import count
sage: P = 2 * polytopes.cube()
Counting integer points from either the H or V representation::
sage: count(P.cdd_Hrepresentation(), cdd=True) # optional - latte_int
125
sage: count(P.cdd_Vrepresentation(), cdd=True) # optional - latte_int
125
Ehrhart polynomial::
sage: count(P.cdd_Hrepresentation(), cdd=True, ehrhart_polynomial=True) # optional - latte_int
64*t^3 + 48*t^2 + 12*t + 1
Multivariate generating function currently only work with ``raw_output=True``::
sage: opts = {'cdd': True,
....: 'multivariate_generating_function': True,
....: 'raw_output': True}
sage: cddin = P.cdd_Hrepresentation()
sage: print(count(cddin, **opts)) # optional - latte_int
x[0]^2*x[1]^(-2)*x[2]^(-2)/((1-x[1])*(1-x[2])*(1-x[0]^(-1)))
+ x[0]^(-2)*x[1]^(-2)*x[2]^(-2)/((1-x[1])*(1-x[2])*(1-x[0]))
+ x[0]^2*x[1]^(-2)*x[2]^2/((1-x[1])*(1-x[0]^(-1))*(1-x[2]^(-1)))
+ x[0]^(-2)*x[1]^(-2)*x[2]^2/((1-x[1])*(1-x[0])*(1-x[2]^(-1)))
+ x[0]^2*x[1]^2*x[2]^(-2)/((1-x[2])*(1-x[0]^(-1))*(1-x[1]^(-1)))
+ x[0]^(-2)*x[1]^2*x[2]^(-2)/((1-x[2])*(1-x[0])*(1-x[1]^(-1)))
+ x[0]^2*x[1]^2*x[2]^2/((1-x[0]^(-1))*(1-x[1]^(-1))*(1-x[2]^(-1)))
+ x[0]^(-2)*x[1]^2*x[2]^2/((1-x[0])*(1-x[1]^(-1))*(1-x[2]^(-1)))
TESTS:
Testing raw output::
sage: from sage.interfaces.latte import count
sage: P = polytopes.cuboctahedron()
sage: cddin = P.cdd_Vrepresentation()
sage: count(cddin, cdd=True, raw_output=True) # optional - latte_int
'19'
sage: count(cddin, cdd=True, raw_output=True, ehrhart_polynomial=True) # optional - latte_int
' + 1 * t^0 + 10/3 * t^1 + 8 * t^2 + 20/3 * t^3'
sage: count(cddin, cdd=True, raw_output=True, multivariate_generating_function=True) # optional - latte_int
'x[0]^(-1)*x[1]^(-1)/((1-x[0]*x[2])*(1-x[0]^(-1)*x[1])*...x[0]^(-1)*x[2]^(-1)))\n'
Testing the ``verbose`` option::
sage: n = count(cddin, cdd=True, verbose=True, raw_output=True) # optional - latte_int
This is LattE integrale ...
...
Invocation: count '--redundancy-check=none' --cdd /dev/stdin
...
Total Unimodular Cones: ...
Maximum number of simplicial cones in memory at once: ...
<BLANKLINE>
**** The number of lattice points is: ****
Total time: ... sec
"""
from subprocess import Popen, PIPE
from sage.misc.misc import SAGE_TMP
from sage.rings.integer import Integer

args = ['count']
if ehrhart_polynomial and multivariate_generating_function:
raise ValueError
if ehrhart_polynomial:
args.append('--ehrhart-polynomial')
elif multivariate_generating_function:
args.append('--multivariate-generating-function')

if 'redundancy_check' not in kwds:
args.append('--redundancy-check=none')

for key,value in kwds.items():
if value is None or value is False:
continue

key = key.replace('_','-')
if value is True:
args.append('--{}'.format(key))
else:
args.append('--{}={}'.format(key, value))

if multivariate_generating_function:
from sage.misc.temporary_file import tmp_filename
filename = tmp_filename()
with open(filename, 'w') as f:
f.write(arg)
args += [filename]
else:
args += ['/dev/stdin']

try:
# The cwd argument is needed because latte
# always produces diagnostic output files.
latte_proc = Popen(args,
stdin=PIPE, stdout=PIPE,
stderr=(None if verbose else PIPE),
cwd=str(SAGE_TMP))
except OSError:
from sage.misc.package import PackageNotFoundError
raise PackageNotFoundError('latte_int')

ans, err = latte_proc.communicate(arg)
ret_code = latte_proc.poll()
if ret_code:
if err is None:
err = ", see error message above"
else:
err = ":\n" + err
raise RuntimeError("LattE integrale program failed (exit code {})".format(ret_code) + err.strip())


if ehrhart_polynomial:
ans = ans.splitlines()[-2]
if raw_output:
return ans
else:
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.rational_field import QQ
R = PolynomialRing(QQ, 't')
return R(ans)
elif multivariate_generating_function:
with open(filename + '.rat') as f:
ans = f.read()
if raw_output:
return ans
else:
raise NotImplementedError("there is no Sage object to handle multivariate series from LattE, use raw_output=True")
else:
ans = ans.splitlines()[-1]
if not ans:
# opening a file is slow (30e-6s), so we read the file
# numOfLatticePoints only in case of a IndexError above
with open(SAGE_TMP+'/numOfLatticePoints', 'r') as f:
ans = f.read()

if raw_output:
return ans
else:
from sage.rings.integer import Integer
return Integer(ans)

0 comments on commit 84abd08

Please sign in to comment.