Polyhedron.volume() ZeroDivisionError

edit

The following answer is only a workaround in a special case, not a general solution. It seems that the implemented method volume is buggy with a good rate, the bug is related to the method triangulate of a polyhedron!

The SHORT answer:

Use the following function, which splits the polyhedron into parts, then computes their volume, the bug shows when trying to triangulate! This is only important message of this answer!

def triangulated_volume( P ):
    """compute the triangulation of P, compute the volume as a sum of the pieces.
    """
    d = P.dimension()
    vertices = [ list(v) for v in P.vertices() ]
    return sum( [ abs( matrix( d+1, d+1, [ [ 1, ] + vertices[ t ] for t in T ] ).det() )
                  for T in P.triangulate() ] ) / factorial(d)

Alternatively, one can try to split into pyramids based in origin, which is an interior point in our case. But this fails relatively often!

def my_volume( P, interior_point = vector( [0,0,0] ) ):
    """compute the faces of P, use the interior point to split it
    into faces pyramids, compute the volume of each pyramid, sum them up.
    """
    return sum( [ Polyhedron( [ interior_point, ]
                              + [ vector( v )
                                  for face in P.faces(2)
                                  for v in face.vertices() ] ).volume() ] )

Finally, starting with a vector c with components $(x,y,z)$, $0 \leq x \leq y\leq z$, for example ( 0.1, 1.0, 1.1 ), and associating the polyhedron for it as in the posted question, then the volume is given by the formula:

4*x^2*y - 8/3*y^3 - 8*x^2*z + 4*y^2*z + 4*y*z^2 + 4/3*z^3

Of course, this is always working! That is all.

LOOONG answer:

This goes deeper (and also far away), but it may be relevant.

(And there was a lot of geometric fun while studying the object. Thanks for the question!)

(Moreover, in m.h.o. it is important to support computer algebra system with illustrative examples. We are in a century of all possibilities for the software, by the "good", inspiring examples are missing.)

The plan of the longer answer is as follows.

1. Confirm the bug and make a statistic of errors for some randomly generated data.
2. Compute the volume starting from the vector vector( QQ, [ 13/17, 14/17, 19/17 ] ) as in the posted question. No error.
3. Working over RR, or with a random example runs into a ZeroDivisionError.
4. Define an ad-hoc version of the volume, that works in our case. The one in the short answer. This ad-hoc function my_volume has an implementation based on the method volume() of a special class of polyhedra, pyramids, where it seems to work well, and the fact we already know an interior point for our polyhedra.
5. Check that the same result is obtained for volume and my_volume and triangulated_volume over QQ.
6. Find the formula for the volume of the specific polyhedra associated to the group $G$ and a starting vector. This is a mathematical issue, but visualization and computations are supported by sage, and sage is a pretty good mean to put the hands on the shape of the polyhedron. (Without it, the trial and error cycles would have been too many in my case.)
7. Check that my_volume (over QQ) and triangulated_volume always gets this formula.

ONE: We start by confirming the error:

import random
import traceback

A = diagonal_matrix( 3, [-1,-1, 1] )
B = diagonal_matrix( 3, [ 1,-1,-1] )
s, t = SymmetricGroup(3).gens()

GROUP = MatrixGroup( [ A, B, s.matrix(), t.matrix() ] )
print "The considered group has order %s" % GROUP.order()

G = GROUP.list()
G = list( G )    # since the list method above gives me a tuple...
G.sort()

def vol( c ):
    return Polyhedron( [ g * c for g in G ] ).volume()

vectors  = [ vector( QQ, [ randint( 5, 15 ) for _ in range(3) ] )
             for trial in range( 10 ) ]
vectors += [ vector( RR, [ random.uniform( 5, 15 ) for _ in range(3) ] )
             for trial in range( 10 ) ]

for v in vectors:
    try:
        volv = vol(v)
        print "%s has associated volume %s" % ( v, vol( v ) )
    except ZeroDivisionError:
        # traceback.print_exc()
        print "%s [got ZeroDivisionError]" % v

We let this run and obtain this time:

The considered group has order 24
(5, 6, 14) has associated volume 22808/3
(12, 15, 15) has associated volume 13860
(5, 6, 9) has associated volume 2436
(8, 8, 6) has associated volume 6784/3
(13, 9, 12) has associated volume 28156/3
(7, 9, 7) has associated volume 5800/3
(7, 9, 9) has associated volume 3096
(12, 9, 10) has associated volume 16984/3
(7, 7, 7) has associated volume 2744/3
(14, 8, 5) has associated volume 30448/3
(12.9589728317627, 5.64395666109784, 9.16678680046540) [got ZeroDivisionError]
(7.01088675869899, 9.62861840238589, 7.41339518763128) [got ZeroDivisionError]
(7.71732224544219, 9.83842107910933, 13.1338440570290) [got ZeroDivisionError]
(13.1463788341743, 13.4704614482142, 11.4765415343236) [got ZeroDivisionError]
(5.68466999612573, 14.0075400171660, 14.7885047686280) [got ZeroDivisionError]
(9.13414931152926, 11.6321442939304, 10.2356341843752) [got ZeroDivisionError]
(8.20989707720079, 8.72810914926507, 9.28924497131092) has associated volume 6534.91775036
(10.3005535665733, 11.1392338434077, 8.80644198437511) [got ZeroDivisionError]
(9.92815378030292, 7.91391408097008, 12.3418545472610) [got ZeroDivisionError]
(5.66242280708254, 13.8289596505955, 12.1218231049418) has associated volume 44714.4988912

So we recover the error in most randomly generated cases.

The "group" G (a list) will be used in the sequel.

TWO: Recall that for the vector( QQ, [ 13/17, 14/17, 19/17 ] ) we can compute the volume without error:

c = vector( QQ, [ 13/17, 14/17, 19/17 ] )
print Polyhedron( [ g * c for g in G ] ).volume()

This gives:

20716/4913

THREE: The ZeroDivisionError occurs over RR...

c = vector( RR, [ 13/17., 14/17., 19/17. ] )
print "c =", c
try:
    P = Polyhedron( [ g * c for g in G ] )
    vol = P.volume()
    print "The volume is ", vol
except ZeroDivisionError:
    print "... sorry, ZeroDivisionError"

This gives:

c = (0.764705882352941, 0.823529411764706, 1.11764705882353)
... sorry, ZeroDivisionError

And further:

try:
    P.triangulate()
except ZeroDivisionError:
    print "... sorry, ZeroDivisionError during P.triangulate()"

and the expected echo is comming:

... sorry, ZeroDivisionError during P.triangulate()

FOUR & FIVE: Let us define our own functions my_volume and triangulated_volume as in the short answer, based on the fact that the volume method seems to work better for pyramids, and always -- avoiding it via det for the computed triangulation. Then we can compare in a random case:

c = vector( [ 12034702141, 13274012341, 230946121123212341 ] )
P = Polyhedron( [ g * c for g in G ] )
print "The volume of P = %s" % P.volume() 
print "My  volume of P = %s" % my_volume(P)
print "Triangulated... = %s" % triangulated_volume(P)

This gives:

The volume of P = 16423693391061491154712244902262907888642459299750456
My  volume of P = 16423693391061491154712244902262907888642459299750456
Triangulated... = 16423693391061491154712244902262907888642459299750456

SIX: Let us start with a general vector c of the form $(x,y,z)$. Without restricting the generality, we may and do assume $0 < x\leq y \leq z$ . Then the faces of the associated polyhedron are...

• eight hexagons, four plus four, and six rectangles.

• One hexagon of the first four group has the vertices $\sigma(x,y,z)$, $\sigma$ in S(3) being a permutation that acts on the components. The six vertices live in the $(X,Y,Z)$ plane $X+Y+Z=x+y+z$, the hexagon $$(x,y,z)\to(x,z,y)\to(y,z,x)\to(z,y,x)\to(z,x,y)\to(y,x,z)\to(x,y,z)$$ (only $(yz)$ or $(xy)$ appear as permutation for each side)
  is centered in $$G=(g,g,g)\ , \ g=\frac 13(x+y+z)\ ,$$ $OG$ is perpendicular to the plane and of length $ g\sqrt 3 $, the sides have length either $a=(z-y)\sqrt2$ or $b=(y-x)\sqrt2$, they alternate $a,b,a,b,a,b$, and we can draw the full lines supporting the segments $a,a,a $ inside to get an equilateral triangle with the side $a+2b$. The area of the hexagon is thus $$\frac{\sqrt3}4(\ (a+2b)^2-3b^2\ ) = \frac{\sqrt3}4(\ a^2+4ab+b^2\ ) = \frac{\sqrt3}2(\ (z-y)^2+4(z-y)(y-x)+(y-x)^2\ ) \ .$$ We know $OG$, and thus the volume: $$ \frac 16(x+y+z)\left[\ (z-y)^2+4(z-y)(y-x)+(y-x)^2\ \right] \ .$$

  The eight planes have equations $X+Y+Z=\pm x\pm y\pm z$, and four out of the eight hexagons are of the same shape as computed above.

• One of the other four hexagons has the vertices $$(-x,y,-z)\to(-x,z,-y)\to(y,z,x)\to(z,y,x)\to(z,-x,-y)\to(y,-x,-z)\dots $$ (twisted permutations), and is centered in $$H=(h,h,-h)\ , \ h=\frac 13(y+z-x)\ ,$$ $OH$ is perpendicular to the plane and of length $ h\sqrt 3 $, the sides have length either $a=(z-y)\sqrt2$ or $b=(y+x)\sqrt2$, they alternate $a,b,a,b,a,b$, and we can draw the full lines supporting the segments $a,a,a $ inside to get an equilateral triangle with the side $a+2b$. The area of the hexagon is thus $$\frac{\sqrt3}4(\ (a+2b)^2-3b^2\ ) = \frac{\sqrt3}4(\ a^2+4ab+b^2\ ) = \frac{\sqrt3}2(\ (z-y)^2+4(z-y)(y+x)+(y+x)^2\ ) \ .$$ We know $OH$, and thus the volume: $$ \frac 16(-x+y+z)\left[\ (z-y)^2+4(z-y)(y+x)+(y+x)^2\ \right] \ .$$

• The six rectangles are of the same shape now, one of them has vertices ${z}\times{\ (x,y),(y,x),(-x,-y),(-y,-x)\ }$, the other being obtained by (the action of $G$, so by) changing the component where $z$ is placed in, and some signs. The corresponding six pyramids of the same shape have each height $z$, and base area $2(y^2-x^2)$. Each volume is thus $$\frac 23z(y+x)(y-x)\ .$$

• This gives a formula for the total volume.

For this we define:

def V1(x,y,z):    return 2/3*z*(y-x)*(y+x)
def V2(x,y,z):    return 1/6*( x+y+z)*((z-y)^2+4*(z-y)*(y-x)+(y-x)^2)
def V3(x,y,z):    return 1/6*(-x+y+z)*((z-y)^2+4*(z-y)*(y+x)+(y+x)^2)
def V (x,y,z):    return 6*V1(x,y,z) + 4*V2(x,y,z) + 4*V3(x,y,z)

def Vinfo(x,y,z): print V1(x,y,z), V2(x,y,z), V3(x,y,z), V(x,y,z)

And the formula is:

var( 'x,y,z' );
V(x,y,z).factor()

This gives:

4*x^2*y - 8/3*y^3 - 8*x^2*z + 4*y^2*z + 4*y*z^2 + 4/3*z^3

Let us verify this formula. It may seem useful to get an explicit example an require from sage all information on all computed volumes. It is simple to rewrite my_volume to get also the information during the computations, for instance:

def info_my_volume( P, interior_point = vector( [0,0,0] ) ):
    """compute the faces of P, use the interior point to split it
    into faces pyramids, compute the volume of each pyramid, sum them up.
    """
    vol = 0
    for face in P.faces(2):
        Q = Polyhedron( [ interior_point, ]
                        + [ vector( v )for v in face.vertices() ] )
        print "Face: %s %s vol = %s" % ( face
                                         , [ list(v) for v in face.vertices() ]
                                         , Q.volume() )
        vol += Q.volume()
    return vol

Then we have in a concrete example:

x,y,z = 1, 6, 9
c = vector( [ x,y,z ] )
P = Polyhedron( [ g*c for g in G ] )
vol = info_my_volume( P )
print "P.volume() =", P.volume()
print "Our volume =", vol

And this gives:

Face: <0,1,2,3> [[-9, -6, 1], [-9, -1, 6], [-9, 1, -6], [-9, 6, -1]] vol = 210
Face: <0,1,4,5,8,9> [[-9, -6, 1], [-9, -1, 6], [-6, -9, 1], [-6, -1, 9], [-1, -9, 6], [-1, -6, 9]] vol = 752/3
Face: <2,3,6,7,10,11> [[-9, 1, -6], [-9, 6, -1], [-6, 1, -9], [-6, 9, -1], [-1, 6, -9], [-1, 9, -6]] vol = 752/3
Face: <0,2,4,6,12,13> [[-9, -6, 1], [-9, 1, -6], [-6, -9, 1], [-6, 1, -9], [1, -9, -6], [1, -6, -9]] vol = 994/3
Face: <1,3,5,7,14,15> [[-9, -1, 6], [-9, 6, -1], [-6, -1, 9], [-6, 9, -1], [1, 6, 9], [1, 9, 6]] vol = 994/3
Face: <4,8,12,16> [[-6, -9, 1], [-1, -9, 6], [1, -9, -6], [6, -9, -1]] vol = 210
Face: <6,10,13,17> [[-6, 1, -9], [-1, 6, -9], [1, -6, -9], [6, -1, -9]] vol = 210
Face: <5,9,14,18> [[-6, -1, 9], [-1, -6, 9], [1, 6, 9], [6, 1, 9]] vol = 210
Face: <7,11,15,19> [[-6, 9, -1], [-1, 9, -6], [1, 9, 6], [6, 9, 1]] vol = 210
Face: <12,13,16,17,20,21> [[1, -9, -6], [1, -6, -9], [6, -9, -1], [6, -1, -9], [9, -6, -1], [9, -1, -6]] vol = 752/3
Face: <8,9,16,18,20,22> [[-1, -9, 6], [-1, -6, 9], [6, -9, -1], [6, 1, 9], [9, -6, -1], [9, 1, 6]] vol = 994/3
Face: <10,11,17,19,21,23> [[-1, 6, -9], [-1, 9, -6], [6, -1, -9], [6, 9, 1], [9, -1, -6], [9, 6, 1]] vol = 994/3
Face: <14,15,18,19,22,23> [[1, 6, 9], [1, 9, 6], [6, 1, 9], [6, 9, 1], [9, 1, 6], [9, 6, 1]] vol = 752/3
Face: <20,21,22,23> [[9, -6, -1], [9, -1, -6], [9, 1, 6], [9, 6, 1]] vol = 210
P.volume() = 3588
Our volume = 3588

The information agrees with:

sage: Vinfo(x,y,z)
210 752/3 994/3 3588

SEVEN: Even now after performing all computations and after some dozen P.show() calls for different shapes, i am not quite sure that the splitted computation always works, and the formula above is recovering the volume. Let us make some 100 tests over QQ.

okCounter = 0    # so far
def my_formula( x,y,z ):
    return 4*x^2*y - 8/3*y^3 - 8*x^2*z + 4*y^2*z + 4*y*z^2 + 4/3*z^3

NTRIALS = 100
for test in xrange(NTRIALS):
    c = [ randint( 1, 10000 ) for _ in range(3) ]
    c . sort()
    x, y, z = c
    c = vector( QQ, c )
    P = Polyhedron( [ g*c for g in G ] )

    vol1 = P.volume()
    vol2 = my_volume( P )
    vol3 = triangulated_volume( P )
    vol4 = triangulated_volume( P )
    # print "[%s] %s %s %s -> %s" % ( test, x, y, z, vol1 )
    if len( { vol1, vol2, vol3, vol4 } ) == 1:
        okCounter += 1    # ok, all answers coincide!

print "%s trials, okCounter = %s" % ( NTRIALS, okCounter )

And i've got this time:

100 trials, okCounter = 100