Polytope volume function engines produce different results

[sagetrac-wishcow]

The volume currently does not distinguish between ambient and induced volume. This should be changed.

Different engine yield very different results, and example:
The lrs engine for polytope volume calculates volume respective to the dimension of the polytope,
while the auto engine calculates volume respective to the dimension of the ambient space.

Example:

m3 = matrix(ZZ, [[0,0,0],[0,0,1]])
p = Polyhedron(m3)
p.volume(engine="lrs")
p.volume()

1.0
0

(Note: The lrs allows calculation of volumes of facets without reducing the dimension of the ambient space, but it uses non-isometrical projections!)

The suggested resolution adds a parameter measure that essentially behaves as follows:

sage: P = Polyhedron([[0, 0], [1, 1]])
sage: P.volume()
0
sage: P.volume(measure='induced')                                                                                                                                   
sqrt(2)
sage: P.volume(measure='induced_rational') # optional -- latte_int
1

Depends on #20887
Depends on #22804

CC: @sagetrac-jakobkroeker @mkoeppe @videlec @jplab @mforets @mo271 @seblabbe

Component: geometry

Keywords: polytope, volume, lrs_volume, days88

Stopgaps: wrongAnswerMarker

Author: Moritz Firsching

Branch/Commit: f5d47cf

Reviewer: Jean-Philippe Labbé

Issue created by migration from https://trac.sagemath.org/ticket/16045

[sagetrac-jakobkroeker]

Stopgaps: wrongAnswerMarker

[sagetrac-jakobkroeker]

comment:1

what is a rational polytope? defined over rationals?

At least I could not get from documentation that using engine="lrs"
implies calculates volume respective to the dimension of the polytope
and using 'auto' implies calculating volume respective to the dimension of the ambient space,
so for me it is a wrong answer issue.

[jplab]

comment:2

This is actually very useful to have this. (I was actually looking for such a function just now).

It should perhaps be better documented with an optional parameter to switch this on and off and let lrs compute the volume of faces.

[mkoeppe]

comment:3

I agree, a keyword parameter that controls the measure used for lower-dim faces would be useful.
There are currently 3 options:

• 0
• answer according to standard lower-dim measure
• LattE always returns a rational answer, by normalizing with respect to the intersection lattice.

[mo271]

comment:4

A few questions:

• What should be the name for the parameter (which is essentially boolean) that controls what measure to be used?
  I can think of inherent, induced, relative, ambient... (I think instrinsic would not be a good idea.)

• Should the default for the parameter (if none is given) depend on the engine chosen?

• What should happen if the engine is not able to compute the inherent volume?

[videlec]

comment:6

Replying to @mo271:

A few questions:

• What should be the name for the parameter (which is essentially boolean) that controls what measure to be used? I can think of inherent, induced, relative, ambient... (I think instrinsic would not be a good idea.)

It is a troolean (at least for rational polytopes) as mentioned in comment:3 of Marcelo.

• Should the default for the parameter (if none is given) depend on the engine chosen?

No.

• What should happen if the engine is not able to compute the inherent volume?

A ValueError or a RuntimeError.

[mo271]

comment:8

So how about the following solve the problem:

We introduce a parameter "measure" with the following options:

• "ambient" (the default?!)
• "induced"
• "latte_renormalized" (or something that Matthias likes)

I suppose it would be wisest to catch measure="ambient" if the polytope is not full-dimensional and return "0" for all the engines (Even if the engines only can compute the induced volume of the polytope)

Should it could work like:

sage: m3 = matrix(ZZ, [[0,0,0],[0,0,1]])
sage: p = Polyhedron(m3)
sage: p.volume()
0
sage: p.volume(engine="lrs", measure="induced")
1.0
sage: p.volume(engine="lrs", measure="ambient")
0
sage: p.volume(measure="induced")
ValueError or a RuntimeError: "Please choose a different engine."

Maybe it would be ok to choose a differnt engine for the last command automatically.

What do you think?

[videlec]

comment:9

Replying to @mo271:

So how about the following solve the problem:

We introduce a parameter "measure" with the following options:

• "ambient" (the default?!)
• "induced"
• "latte_renormalized" (or something that Matthias likes)

Note that "latte_renormalized" only make sense for polyhedron defined over Q. I guess it would be better "induced_rational" rather than any reference to latte. This renormalization is very natural.

I suppose it would be wisest to catch measure="ambient" if the polytope is not full-dimensional and return "0" for all the engines (Even if the engines only can compute the induced volume of the polytope)

Should it could work like:

sage: m3 = matrix(ZZ, [[0,0,0],[0,0,1]])
sage: p = Polyhedron(m3)
sage: p.volume()
0
sage: p.volume(engine="lrs", measure="induced")
1.0
sage: p.volume(engine="lrs", measure="ambient")
0
sage: p.volume(measure="induced")
ValueError or a RuntimeError: "Please choose a different engine."

Maybe it would be ok to choose a differnt engine for the last command automatically.

+1. I also do not like the behavior of the last command. The default choice of the engine should take into account the measure argument (which should be firt BTW).

def volume(measure="ambient", engine=None):
    r"""
    Return the volume (or relative volume) of the polytope
    """

[mforets]

comment:10

minor remark: ticket #20887 also touches the polyhedron volume function

[mkoeppe]

comment:11

+1 on "induced_rational"

[mo271]

Branch: u/moritz/16045

[mo271]

Author: Moritz Firsching

[mo271]

Commit: 2e2764e

[mo271]

comment:13

I am building on #20887 and introduced a measure option.

Here is a test (make sure ou have latte_int, topcom and lrs installed)

sage: C = polytopes.cube()
sage: P = polytopes.permutahedron(4)
sage: P
A 3-dimensional polyhedron in ZZ^4 defined as the convex hull of 24 vertices
sage: measures = ["ambient", "induced", "induced_rational"]
sage: engines = ["auto", "internal", "TOPCOM", "lrs", "latte"]
sage: for m in measures:
....:     for e in engines:
....:         try:
....:             pvol = P.volume(measure=m, engine=e)
....:         except TypeError as err:
....:             pvol = err
....:         try:
....:             cvol = C.volume(measure=m, engine=e)
....:         except Exception as err:
....:             cvol = err
....:         print 'measure =',m, 'engine =',e 
....:         print 'not-fulldimensional example: ', pvol
....:         print 'full-dimensional example: ', cvol
....:         print '--------------'
....:         
measure = ambient engine = auto
not-fulldimensional example:  0
full-dimensional example:  8
--------------
measure = ambient engine = internal
not-fulldimensional example:  0
full-dimensional example:  8
--------------
measure = ambient engine = TOPCOM
not-fulldimensional example:  0
full-dimensional example:  8
--------------
measure = ambient engine = lrs
not-fulldimensional example:  0
full-dimensional example:  8.0
--------------
measure = ambient engine = latte
not-fulldimensional example:  0
full-dimensional example:  8
--------------
measure = induced engine = auto
not-fulldimensional example:  16.0
full-dimensional example:  8.0
--------------
measure = induced engine = internal
not-fulldimensional example:  The induced measure can only be computed with the engine set to `auto` or `lrs`
full-dimensional example:  8.0
--------------
measure = induced engine = TOPCOM
not-fulldimensional example:  The induced measure can only be computed with the engine set to `auto` or `lrs`
full-dimensional example:  8.0
--------------
measure = induced engine = lrs
not-fulldimensional example:  16.0
full-dimensional example:  8.0
--------------
measure = induced engine = latte
not-fulldimensional example:  The induced measure can only be computed with the engine set to `auto` or `lrs`
full-dimensional example:  8.0
--------------
measure = induced_rational engine = auto
not-fulldimensional example:  16
full-dimensional example:  8
--------------
measure = induced_rational engine = internal
not-fulldimensional example:  The induced rational measure can only be computed with the engine set to `auto` or `latte`
full-dimensional example:  8
--------------
measure = induced_rational engine = TOPCOM
not-fulldimensional example:  The induced rational measure can only be computed with the engine set to `auto` or `latte`
full-dimensional example:  8
--------------
measure = induced_rational engine = lrs
not-fulldimensional example:  The induced rational measure can only be computed with the engine set to `auto` or `latte`
full-dimensional example:  8
--------------
measure = induced_rational engine = latte
not-fulldimensional example:  16
full-dimensional example:  8
--------------

The options are explained as

        - ``measure`` -- string. The measure to use. Allowed values are:

          * ``ambient``: Lebesgue measure of ambient space (volume)
          * ``induced``: Lebesgue measure of affine hull (relative volume)
          * ``induced_rational``: TODO: explanation (only makes sense for integer polytopes?!)

        - ``engine`` -- string. The backend to use. Allowed values are:

          * ``'auto'`` (default): choose engine according to measure
          * ``'internal'``: see :meth:`triangulate`.
          * ``'TOPCOM'``: see :meth:`triangulate`.
          * ``'lrs'``: use David Avis's lrs program (optional).
          * ``'latte'``: use LattE integrale program (optional).

What could I write for induced_rational?

I have not yet added doctests.

Comments welcome!

---

New commits:

17911f7	Added integrate method to Polyhedron base.py
91d885e	Fix the docstrings, and enhance integrate method.
507aea5	Add the new volume engine Polyhedron.
a563192	fix a bug in _volume_latte
9ac8279	reject RDF, but add an example
ec15897	add test for low-dim polytope and fix docstring typo
1cf59f4	fix last doctest in integrate (exception mismatch)
dc544fa	fix optional tag in integrate test
31ea148	Merge branch 'u/mforets/20887' of git://trac.sagemath.org/sage into 16045
2e2764e	add measure option

[videlec]

comment:14

What could I write for induced_rational?

(I think this is describe in the latte manual) Given the polytope P in Z^d it spans an affine space E. The ambient measure is the only scaling of the Lebesgue measure so that the lattice E \cap Z^d has covolume 1 in E.

[videlec]

comment:15

If you merge commits from #20887 here you should set it as a dependency.

[mkoeppe]

comment:16

Note that it's also defined for rational polytopes; and more generally for any (possibly irrational) translation of those.

[mo271]

Dependencies: #20887

[mo271]

comment:19

I noticed that "lrs" is in fact not computing the Lebesgue measure of the affine hull. Form the lrs web-page:

If the option is applied to a V-representation of a polytope that is not full dimensional, the volume of a projected polytope is computed. The projection used is to the lexicographically smallest coordinate subspace, see Avis, Fukuda, Picozzi (2002).

In fact, this shows already in one of the old doctests:

sage: I = Polyhedron([(0,0), (1,1)])
sage: I.volume(engine='lrs') #optional - lrslib
1.0

The result with the induced measure should be sqrt(2) (or 1.414213562373095 I guess)

So this should not be used to compute the volume of facets, as mentioned in comment:2 !

In any case I still think it would be useful to have the option for this measure, except that at the moment, none of the engines is able to actually compute appropriate induced volumes.

Perhaps it would be useful to open a ticket to ask for a function, that does compute the induced measure?! A naive implementation could choose an affine basis of vertices and then some isometry.

[jplab]

comment:20

It would be nice to have an explanation of these subtleties in the doc of the function (as shown by Moritz). Consequently, the issue described in the ticket would be resolved so as to address the reason why it returns different results.

IMHO, I guess it is okay to have one method returning different results depending on the value of the parameters. So +1 to have an optional parameter which specifies the measure (and hence the engine in some cases) as Vincent mentionned.

Finally, add a TODO that would ask for the lebesgue measure which would be addressed in a separate ticket and explain what lrs is doing with a possible link to the webpage.

What do you think?

[mo271]

comment:21

One thing that would be useful to have, would be a function, that takes a non full-dimensional poytope and returns a full-dimensional one that is affinely equivalent to the original one. Note that the usual .affine_hull method does not achive that goal:

sage: p = Polyhedron([[0,0],[1,1]])
sage: p.affine_hull()
A 1-dimensional polyhedron in ZZ^1 defined as the convex hull of 2 vertices
sage: p.affine_hull().vertices()
(A vertex at (0), A vertex at (1))

The length of the affine hull is not sqrt(2) as one might expect. Another example would be

sage: s = polytopes.simplex()
sage: s
A 3-dimensional polyhedron in ZZ^4 defined as the convex hull of 4 vertices
sage: s.affine_hull()
A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices
sage: _.show()
Launched jmol viewer for Graphics3d Object

Here the simplex is not regular after using .affine_hull.

Consider this new affine_hull function:

def affine_hull(P, preserve_volume=True):
    if P.ambient_dim() == P.dim():
        return P
    Q = P.translation(-vector(P.vertices()[0]))
    v = Q.vertices()[0]
    assert list(v) == P.ambient_dim()*[0]
    import itertools
    M = Matrix([list(w) for w in itertools.islice(v.neighbors(), P.dim())])
    if P.base_ring() in [ZZ,QQ]:
        M = Matrix(AA,M)
    A = M.gram_schmidt(orthonormal = True)[0]
    return Polyhedron([A*vector(v) for v in Q.vertices()])

Using this gets the desired results:

sage: p = Polyhedron([[0,0],[1,1]])
sage: affine_hull(p)
A 1-dimensional polyhedron in AA^1 defined as the convex hull of 2 vertices
sage: affine_hull(p).vertices()
(A vertex at (0), A vertex at (1.414213562373095?))

It is actually useful to plot the permutahedron nicely:

sage: P = polytopes.permutahedron(4)
sage: P
A 3-dimensional polyhedron in ZZ^4 defined as the convex hull of 24 vertices
sage: affine_hull(P)
A 3-dimensional polyhedron in AA^3 defined as the convex hull of 24 vertices
sage: _.plot()
Launched jmol viewer for Graphics3d Object

(There is an inherent problem, that considering a rational polytope as full-dimensional polytope in its affine hull might lead to non-rational polyopes: basically you want to work in a field with square roots.)

In any case I propose to introduce an improved version of the .affine_hull method, perhaps with an parameter preserve_volume=True (or isometry=True) an then use this method to handle the cases of the volume of non full-dimenional polytopes here.

What do you think?

[videlec]

comment:22

Replying to @mo271:

One thing that would be useful to have, would be a function, that takes a non full-dimensional poytope and returns a full-dimensional one that is affinely equivalent to the original one.

[...]

What do you think?

I agree that it would be useful and essentially straightforward (ignoring details about base ring). You can open a new ticket and move the discussion there.

[mo271]

comment:23

Replying to @videlec:

I agree that it would be useful and essentially straightforward (ignoring details about base ring). You can open a new ticket and move the discussion there.

I opened #22804 and already added some code. This could eventually become a depency here (#16045) and we can end up with a '.volume' method that is more useful...

[mo271]

Changed dependencies from #20887 to #20887 #22804

[seblabbe]

comment:34

Use # abs tol marker in the comment. See Special Markup to Influence Doctests in the documentation.

[sagetrac-git]

Changed commit from 1a8cefc to fa46cee

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

b602072	add measure option
09fb035	using new affine_hull function to provide induced measure
5b68b13	added one forgotten line in a doctest
5f452f5	doctest for 'induced_rational'
b7b10da	small improvements suggested by JP
fa46cee	abs tol marker added

[mo271]

comment:36

Thanks, Sébastian!
I added the marker..

[jplab]

Changed keywords from polytope, volume, lrs_volume to polytope, volume, lrs_volume, days88

[jplab]

comment:38

There is a hidden failed test (probably got caught by chance by the bot):

sage -t --long src/sage/geometry/polyhedron/backend_normaliz.py
**********************************************************************
File "src/sage/geometry/polyhedron/backend_normaliz.py", line 412, in sage.geometry.polyhedron.backend_normaliz.Polyhedron_normaliz.integral_hull
Failed example:
    set(PI.Vrepresentation())                                # optional - pynormaliz
Expected:
    {A vertex at (-1, 0), A vertex at (0, 1), A ray in the direction (1, 0)}
Got:
    {A ray in the direction (1, 0), A vertex at (-1, 0), A vertex at (0, 1)}
**********************************************************************
File "src/sage/geometry/polyhedron/backend_normaliz.py", line 420, in sage.geometry.polyhedron.backend_normaliz.Polyhedron_normaliz.integral_hull
Failed example:
    set(PI.Vrepresentation())                                # optional - pynormaliz
Expected:
    {A vertex at (1, 0),
     A ray in the direction (1, 0),
     A line in the direction (1, -1)}
Got:
    {A line in the direction (1, -1),
     A ray in the direction (1, 0),
     A vertex at (1, 0)}
**********************************************************************
1 item had failures:
   2 of  10 in sage.geometry.polyhedron.backend_normaliz.Polyhedron_normaliz.integral_hull
    [100 tests, 2 failures, 4.39 s]

Question: I have "normaliz" in my optional installed packages but not "pynormaliz", is the above naming correct then?

I am aware that this is not the topic of this ticket but since the tests failed here, it is probably best to treat this here.

[mkoeppe]

comment:40

Replying to @jplab:

    set(PI.Vrepresentation())                                # optional - pynormaliz

Question: I have "normaliz" in my optional installed packages but not "pynormaliz", is the above naming correct then?

The only way that sage accesses normaliz is via pynormaliz, so this is correct.

[mo271]

comment:41

I am aware that this is not the topic of this ticket but since the tests failed here, it is probably best to treat this here.

I disagree. This should not be treated here, but a new ticket should be opened. If we proceed as you propose (and everybody else does the same), this random pynormaliz bug will be treated in every ticket that is around.

In fact it already seems to be fixed here:

#23586

see also:

#23637

and

sagemath/sagetrac-mirror@effcedb...cb2fcfd

[jplab]

comment:42

The only way that sage accesses normaliz is via pynormaliz, so this is correct.

My bad, I confused both and only installed normaliz. Oh well...

[jplab]

Changed branch from u/moritz/16045 to u/jipilab/16045

[jplab]

Changed commit from fa46cee to f5d47cf

[jplab]

comment:44

I disagree. This should not be treated here, but a new ticket should be opened. If we proceed as you propose (and everybody else does the same), this random pynormaliz bug will be treated in every ticket that is around.

In fact it already seems to be fixed here:

#23586

see also:

#23637

and

sagemath/sagetrac-mirror@effcedb...cb2fcfd

All right, good to know, thanks for looking it up!

Should we add a dependancy in such cases, or perhaps just wait...?

I corrected two indentations. This looks good to me now, tests pass on my machine with pynormaliz (!).

I extended the description of the ticket a bit, if you believe it should contain some more information you can add it, otherwise, I'd it is ready to go.

---

New commits:

2b0fb6d	Merge branch 'u/moritz/16045' of git://trac.sagemath.org/sage into 16045
f5d47cf	Corrected some indentations

[mo271]

comment:45

Thanks for the review, JP! I don't think we should put a dependency for the pynormaliz-set-thing.

[vbraun]

Changed branch from u/jipilab/16045 to f5d47cf