Opened 8 years ago

Last modified 4 years ago

#16045 closed defect

Polytope volume function engines produce different results — at Version 27

Reported by: wishcow Owned by:
Priority: minor Milestone: sage-8.1
Component: geometry Keywords: polytope, volume, lrs_volume, days88
Cc: jakobkroeker, mkoeppe, vdelecroix, jipilab, mforets, moritz, slabbe Merged in:
Authors: Moritz Firsching Reviewers:
Report Upstream: N/A Work issues:
Branch: u/moritz/16045 (Commits, GitHub, GitLab) Commit: 1f64267248f427828b1d0225b947086b8540e619
Dependencies: #20887 #22804 Stopgaps: wrongAnswerMarker

Status badges

Description (last modified by moritz)

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!)

Change History (27)

comment:1 Changed 5 years ago by jakobkroeker

  • Cc jakobkroeker added
  • Stopgaps set to wrongAnswerMarker

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.

comment:2 Changed 5 years ago by jipilab

  • Cc mkoeppe vdelecroix jipilab added

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.

comment:3 Changed 5 years ago by mkoeppe

  • Cc mforets added

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.

comment:4 follow-up: Changed 5 years ago by moritz

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?

comment:5 Changed 5 years ago by moritz

  • Cc moritz added

comment:6 in reply to: ↑ 4 Changed 5 years ago by vdelecroix

Replying to moritz:

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.

comment:7 Changed 5 years ago by slabbe

  • Cc slabbe added

comment:8 follow-up: Changed 5 years ago by moritz

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?

comment:9 in reply to: ↑ 8 Changed 5 years ago by vdelecroix

Replying to moritz:

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
    """

comment:10 Changed 5 years ago by mforets

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

comment:11 Changed 5 years ago by mkoeppe

+1 on "induced_rational"

comment:12 Changed 5 years ago by moritz

  • Branch set to u/moritz/16045

comment:13 Changed 5 years ago by moritz

  • Authors set to Moritz Firsching
  • Commit set to 2e2764eed2b78dde93138f613ca98712ecc9d3e4
  • Status changed from new to needs_info

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:

17911f7Added integrate method to Polyhedron base.py
91d885eFix the docstrings, and enhance integrate method.
507aea5Add the new volume engine Polyhedron.
a563192fix a bug in _volume_latte
9ac8279reject RDF, but add an example
ec15897add test for low-dim polytope and fix docstring typo
1cf59f4fix last doctest in integrate (exception mismatch)
dc544fafix optional tag in integrate test
31ea148Merge branch 'u/mforets/20887' of git://trac.sagemath.org/sage into 16045
2e2764eadd measure option

comment:14 Changed 5 years ago by vdelecroix

What could I write for induced_rational?

(I think this is describe in the latte manual) Given the polytope P in Zd 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.

Last edited 5 years ago by vdelecroix (previous) (diff)

comment:15 Changed 5 years ago by vdelecroix

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

comment:16 Changed 5 years ago by mkoeppe

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

comment:17 Changed 5 years ago by moritz

  • Dependencies set to #20887

comment:18 Changed 5 years ago by mkoeppe

  • Status changed from needs_info to needs_work

comment:19 Changed 5 years ago by moritz

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.

comment:20 Changed 5 years ago by jipilab

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?

comment:21 follow-up: Changed 5 years ago by moritz

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?

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

Replying to moritz:

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.

comment:23 in reply to: ↑ 22 Changed 5 years ago by moritz

Replying to vdelecroix:

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...

comment:24 Changed 4 years ago by moritz

  • Dependencies changed from #20887 to #20887 #22804

comment:25 Changed 4 years ago by git

  • Commit changed from 2e2764eed2b78dde93138f613ca98712ecc9d3e4 to 7a728f10bf75e2e1233a2901e6b831336b5f38f2

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

f14a9e8add measure option
7a728f1using new affine_hull function to provide induced measure

comment:26 Changed 4 years ago by git

  • Commit changed from 7a728f10bf75e2e1233a2901e6b831336b5f38f2 to 1f64267248f427828b1d0225b947086b8540e619

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

1f64267added one forgotten line in a doctest

comment:27 Changed 4 years ago by moritz

  • Description modified (diff)
  • Status changed from needs_work to needs_info

I now used the new affine hull function to provide induced volumes; comments welcome!

I would like to include a good example / doctest for the induced rational volume, any hints? (Also to improve the description ind the section "Input"?)

Note: See TracTickets for help on using tickets.