Opened 7 years ago
Closed 4 years ago
#16045 closed defect (fixed)
Polytope volume function engines produce different results
Reported by:  wishcow  Owned by:  

Priority:  minor  Milestone:  sage8.1 
Component:  geometry  Keywords:  polytope, volume, lrs_volume, days88 
Cc:  jakobkroeker, mkoeppe, vdelecroix, jipilab, mforets, moritz, slabbe  Merged in:  
Authors:  Moritz Firsching  Reviewers:  JeanPhilippe Labbé 
Report Upstream:  N/A  Work issues:  
Branch:  f5d47cf (Commits, GitHub, GitLab)  Commit:  f5d47cf5576ca2f752aaac47f859b2c997535af3 
Dependencies:  #20887 #22804  Stopgaps:  wrongAnswerMarker 
Description (last modified by )
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 nonisometrical 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
Change History (47)
comment:1 Changed 5 years ago by
 Cc jakobkroeker added
 Stopgaps set to wrongAnswerMarker
comment:2 Changed 5 years ago by
 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
 Cc mforets added
I agree, a keyword parameter that controls the measure used for lowerdim faces would be useful. There are currently 3 options:
 0
 answer according to standard lowerdim measure
 LattE always returns a rational answer, by normalizing with respect to the intersection lattice.
comment:4 followup: ↓ 6 Changed 5 years ago by
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
 Cc moritz added
comment:6 in reply to: ↑ 4 Changed 5 years ago by
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 thinkinstrinsic
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
 Cc slabbe added
comment:8 followup: ↓ 9 Changed 5 years ago by
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 fulldimensional 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
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 fulldimensional 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
minor remark: ticket #20887 also touches the polyhedron volume function
comment:11 Changed 5 years ago by
+1 on "induced_rational"
comment:12 Changed 4 years ago by
 Branch set to u/moritz/16045
comment:13 Changed 4 years ago by
 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 3dimensional 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 'notfulldimensional example: ', pvol ....: print 'fulldimensional example: ', cvol ....: print '' ....: measure = ambient engine = auto notfulldimensional example: 0 fulldimensional example: 8  measure = ambient engine = internal notfulldimensional example: 0 fulldimensional example: 8  measure = ambient engine = TOPCOM notfulldimensional example: 0 fulldimensional example: 8  measure = ambient engine = lrs notfulldimensional example: 0 fulldimensional example: 8.0  measure = ambient engine = latte notfulldimensional example: 0 fulldimensional example: 8  measure = induced engine = auto notfulldimensional example: 16.0 fulldimensional example: 8.0  measure = induced engine = internal notfulldimensional example: The induced measure can only be computed with the engine set to `auto` or `lrs` fulldimensional example: 8.0  measure = induced engine = TOPCOM notfulldimensional example: The induced measure can only be computed with the engine set to `auto` or `lrs` fulldimensional example: 8.0  measure = induced engine = lrs notfulldimensional example: 16.0 fulldimensional example: 8.0  measure = induced engine = latte notfulldimensional example: The induced measure can only be computed with the engine set to `auto` or `lrs` fulldimensional example: 8.0  measure = induced_rational engine = auto notfulldimensional example: 16 fulldimensional example: 8  measure = induced_rational engine = internal notfulldimensional example: The induced rational measure can only be computed with the engine set to `auto` or `latte` fulldimensional example: 8  measure = induced_rational engine = TOPCOM notfulldimensional example: The induced rational measure can only be computed with the engine set to `auto` or `latte` fulldimensional example: 8  measure = induced_rational engine = lrs notfulldimensional example: The induced rational measure can only be computed with the engine set to `auto` or `latte` fulldimensional example: 8  measure = induced_rational engine = latte notfulldimensional example: 16 fulldimensional 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 lowdim 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

comment:14 Changed 4 years ago by
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.
comment:15 Changed 4 years ago by
If you merge commits from #20887 here you should set it as a dependency.
comment:16 Changed 4 years ago by
Note that it's also defined for rational polytopes; and more generally for any (possibly irrational) translation of those.
comment:17 Changed 4 years ago by
 Dependencies set to #20887
comment:18 Changed 4 years ago by
 Status changed from needs_info to needs_work
comment:19 Changed 4 years ago by
I noticed that "lrs" is in fact not computing the Lebesgue measure of the affine hull. Form the lrs webpage:
If the option is applied to a Vrepresentation 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 4 years ago by
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 followup: ↓ 22 Changed 4 years ago by
One thing that would be useful to have, would be a function, that takes a non fulldimensional poytope and returns a fulldimensional 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 1dimensional 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 3dimensional polyhedron in ZZ^4 defined as the convex hull of 4 vertices sage: s.affine_hull() A 3dimensional 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 1dimensional 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 3dimensional polyhedron in ZZ^4 defined as the convex hull of 24 vertices sage: affine_hull(P) A 3dimensional 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 fulldimensional polytope in its affine hull might lead to nonrational 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 fulldimenional polytopes here.
What do you think?
comment:22 in reply to: ↑ 21 ; followup: ↓ 23 Changed 4 years ago by
Replying to moritz:
One thing that would be useful to have, would be a function, that takes a non fulldimensional poytope and returns a fulldimensional 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 4 years ago by
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
 Dependencies changed from #20887 to #20887 #22804
comment:25 Changed 4 years ago by
 Commit changed from 2e2764eed2b78dde93138f613ca98712ecc9d3e4 to 7a728f10bf75e2e1233a2901e6b831336b5f38f2
comment:26 Changed 4 years ago by
 Commit changed from 7a728f10bf75e2e1233a2901e6b831336b5f38f2 to 1f64267248f427828b1d0225b947086b8540e619
Branch pushed to git repo; I updated commit sha1. New commits:
1f64267  added one forgotten line in a doctest

comment:27 followup: ↓ 28 Changed 4 years ago by
 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"?)
comment:28 in reply to: ↑ 27 Changed 4 years ago by
Replying to moritz:
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"?)
Volume of the line segment between (a,0)
and (0,b)
are already good examples. Then you can do that for the triangle (a,0,0)
, (0,b,0)
and (0,0,c)
. And you can consider polytopes obtained by permuting the entries of a vector like
sage: Polyhedron(vertices=[(0,0,1,1),(0,1,1,0),(1,1,0,0)]) A 2dimensional polyhedron in ZZ^4 defined as the convex hull of 3 vertices
comment:29 Changed 4 years ago by
 Commit changed from 1f64267248f427828b1d0225b947086b8540e619 to 545507350c90451a2532b5d082fb09ecf2f031fb
Branch pushed to git repo; I updated commit sha1. New commits:
5455073  doctest for 'induced_rational'

comment:30 Changed 4 years ago by
 Status changed from needs_info to needs_review
comment:31 Changed 4 years ago by
 Reviewers set to JeanPhilippe Labbé
 Status changed from needs_review to needs_work
Hi Moritz,
Small comments:
 In the doc:
measure.Different engines
(add space)  Further, right after this sentence, could there be a short sentence, maybe as a note, to say what "for lattice polytope" means exactly?.
 For the example,
P = Polyhedron([[0, 0], [1, 1]])
I would also add the output with the parameterinduced_rational
. This will make it contrast right away, rather than waiting for polyhedronI
. (It is minor but very instructive).  Perhaps name the second polyhedron
P
differently. if polyhedron is actually fulldimensional, return volume with ambient measuren
(remove then
at the end).
You have some failed tests:
********************************************************************** File "src/sage/geometry/polyhedron/base.py", line 4373, in sage.geometry.polyhedron.base.Polyhedron_base.volume Failed example: w = Dinexact.faces(2)[0].as_polyhedron().volume(measure='induced', engine='internal'); w Expected: 1.534062710738235 Got: 1.5340627107382352
Failed example: I.volume(measure='induced_rational') Exception raised: NotImplementedError: You must install the optional latte_int package for this function to work.
(add # optional  latte_int
)
comment:32 Changed 4 years ago by
 Commit changed from 545507350c90451a2532b5d082fb09ecf2f031fb to 1a8cefc5171c26c4fa27a198a48dd6a661fd59b0
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
e5a01d0  add measure option

9350ef2  using new affine_hull function to provide induced measure

bddd3ed  added one forgotten line in a doctest

a4e62a9  doctest for 'induced_rational'

1a8cefc  small improvements suggested by JP

comment:33 Changed 4 years ago by
 Status changed from needs_work to needs_info
Thanks for the comments, JP!
I fixed everything, except perhaps the following; where I was not quite sure what is going on. (Perhaps somehow different machines have different ideas what a real double is..)
w = Dinexact.faces(2)[0].as_polyhedron().volume(measure='induced', engine='internal'); w Expected: 1.534062710738235 Got: 1.5340627107382352
I hope I managed to add # optional  latte_int
in all relevant places: I didn't test it without latte_int
installed.
comment:34 Changed 4 years ago by
Use # abs tol marker in the comment. See Special Markup to Influence Doctests in the documentation.
comment:35 Changed 4 years ago by
 Commit changed from 1a8cefc5171c26c4fa27a198a48dd6a661fd59b0 to fa46ceef2634458cd6dfc4146feb65f1945ea5a2
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

comment:36 Changed 4 years ago by
 Status changed from needs_info to needs_review
Thanks, Sébastian! I added the marker..
comment:37 Changed 4 years ago by
 Keywords days88 added
comment:38 followups: ↓ 40 ↓ 41 Changed 4 years ago by
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.
comment:39 Changed 4 years ago by
 Status changed from needs_review to needs_work
comment:40 in reply to: ↑ 38 ; followup: ↓ 42 Changed 4 years ago by
Replying to jipilab:
set(PI.Vrepresentation()) # optional  pynormalizQuestion: 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.
comment:41 in reply to: ↑ 38 ; followup: ↓ 44 Changed 4 years ago by
 Status changed from needs_work to needs_review
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:
https://trac.sagemath.org/ticket/23586
see also:
https://trac.sagemath.org/ticket/23637
and
comment:42 in reply to: ↑ 40 Changed 4 years ago by
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...
comment:43 Changed 4 years ago by
 Branch changed from u/moritz/16045 to u/jipilab/16045
comment:44 in reply to: ↑ 41 Changed 4 years ago by
 Commit changed from fa46ceef2634458cd6dfc4146feb65f1945ea5a2 to f5d47cf5576ca2f752aaac47f859b2c997535af3
 Description modified (diff)
 Status changed from needs_review to positive_review
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:
https://trac.sagemath.org/ticket/23586
see also:
https://trac.sagemath.org/ticket/23637
and
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

comment:45 Changed 4 years ago by
Thanks for the review, JP! I don't think we should put a dependency for the pynormalizsetthing.
comment:46 Changed 4 years ago by
 Milestone set to sage8.1
comment:47 Changed 4 years ago by
 Branch changed from u/jipilab/16045 to f5d47cf5576ca2f752aaac47f859b2c997535af3
 Resolution set to fixed
 Status changed from positive_review to closed
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.