Possible bug in volume

[itsdfish]

Hi,

I asked about this issue on discourse, but did not receive a reply. I am opening this issue because the result from volume seems to be a bug unless I am missing something about the calculation. Here is a MWE:

Results

2D
Polyhedra.jl: 0.982
QHull.jl: 0.982
3D
Polyhedra.jl: 4.987
QHull.jl: 0.932

Code

using QHull, Polyhedra

points1 = rand(1000, 2) 
ph1 = polyhedron(vrep(points1))

hull1 = chull(points1)

println("2D")
println("Polyhedra.jl: $(round(volume(ph1), digits=3))")
println("QHull.jl: $(round(hull1.volume, digits=3))")


points2 = rand(1000, 3) 
ph2 = polyhedron(vrep(points2))

hull2 = chull(points2)

println("3D")
println("Polyhedra.jl: $(round(volume(ph2), digits=3))")
println("QHull.jl: $(round(hull2.volume, digits=3))")

Version Info

(@v1.7) pkg> st QHull Polyhedra
      Status `~/.julia/environments/v1.7/Project.toml`
  [67491407] Polyhedra v0.6.17
  [a8468747] QHull v0.2.2

julia> versioninfo()
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, haswell)
Environment:
  JULIA_CMDSTAN_HOME = /home/dfish/cmdstan
  JULIA_NUM_THREADS = 4
  JULIA_EDITOR = code

[blegat]

Thanks for reporting this. Can you reproduce it for a simpler explicit example rather than a random polytope of 1000 points ?

[itsdfish]

No problem. I might need some suggestions, as this is outside of my area of expertise. In the code above, tried to make a 1 X 1 square and a 1 X 1 X 1 cube. My understanding is that with a sufficient number of samples, the approximate volume would be 1 in each case. Would it be better to generate points along the perimeter instead?

[blegat]

You should get the correct volume if you give the 8 extreme vertices, does Polyhedra.volume return 1 in this case ?

[itsdfish]

It appears to work correctly in this simple case:

Code

using Polyhedra

vertices = [[i, j, k] for i in 0:1 for j in 0:1 for k in 0:1]

verticies = permutedims(hcat(vertices...))

p = polyhedron(vrep(vertices))

vol = Polyhedra.volume(p)

Result

julia> vol = Polyhedra.volume(p)
1//1

[lxvm]

@itsdfish @blegat I can produce an example with an error in the volume, but it is also possible to fix it using the library.

using SymmetryReduceBZ, Polyhedra
real_latvecs = [1 0; 0 1]
convention="ordinary"
atom_types=[0]
atom_pos = Array([0 0]')
coordinates = "Cartesian"
bzformat = "half-space"
makeprim=false
bz = SymmetryReduceBZ.Symmetry.calc_bz(real_latvecs,atom_types,atom_pos,coordinates,
    bzformat,makeprim,convention) # the square [-0.5,-0.5]x[0.5,0.5]
p = polyhedron(bz, SymmetryReduceBZ.Symmetry.CDDLib.Library())
volume(p) # 2.0, incorrect
removehredundancy!(p)
volume(p) # 1.0, correct

[itsdfish]

@lxvm, thanks for the work around!

[My-Kingdom-for-a-Cat]

It seems the issue is the same as in [#249]: the function triangulation_indices() returns redundant simplices. Removing duplicates as explained here yields the correct volume 1.0.