diff --git a/test/Project.toml b/test/Project.toml index 475ef2ed..3b9d50da 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,5 @@ [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" -CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60" GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" diff --git a/test/volume.jl b/test/volume.jl index 5dad3b50..4954b681 100644 --- a/test/volume.jl +++ b/test/volume.jl @@ -1,6 +1,5 @@ using Test using Polyhedra -using CDDLib # simple algorithm for calculating area of polygons # requires vertices to be sorted (counter)clockwise @@ -15,133 +14,150 @@ function shoelace(verts::AbstractMatrix{<:Real}) A = abs(A) A isa AbstractFloat ? A/2 : A//2 end -isvertsapprox = (verts, points) -> all(any(isapprox(p, v) for v in verts) for p in points) +isvertsapprox(verts, points) = all(any(isapprox(p, v) for v in verts) for p in points) # issue 285: area of square [-0.5, -0.5] x [0.5, 0.5] - -ineq = [ - HalfSpace([-2.0, -2.0], 4.0), - HalfSpace([-2.0, -1.0], 2.5), - HalfSpace([-2.0, 0.0], 2.0), - HalfSpace([-2.0, 1.0], 2.5), - HalfSpace([-2.0, 2.0], 4.0), - HalfSpace([-1.0, -2.0], 2.5), - HalfSpace([-1.0, -1.0], 1.0), - HalfSpace([-1.0, 0.0], 0.5), - HalfSpace([-1.0, 1.0], 1.0), - HalfSpace([-1.0, 2.0], 2.5), - HalfSpace([0.0, -2.0], 2.0), - HalfSpace([0.0, -1.0], 0.5), - HalfSpace([0.0, 0.0], 0.0), - HalfSpace([0.0, 1.0], 0.5), - HalfSpace([0.0, 2.0], 2.0), - HalfSpace([1.0, -2.0], 2.5), - HalfSpace([1.0, -1.0], 1.0), - HalfSpace([1.0, 0.0], 0.5), - HalfSpace([1.0, 1.0], 1.0), - HalfSpace([1.0, 2.0], 2.5), - HalfSpace([2.0, -2.0], 4.0), - HalfSpace([2.0, -1.0], 2.5), - HalfSpace([2.0, 0.0], 2.0), - HalfSpace([2.0, 1.0], 2.5), - HalfSpace([2.0, 2.0], 4.0), -] -square = polyhedron(reduce(intersect, ineq)) -sqverts = [-1 1 1 -1; -1 -1 1 1]/2 -@test isvertsapprox(eachcol(sqverts), points(square)) -@test volume(square) ≈ shoelace(sqverts) +function check_vol_issue285(lib) + ineq = [ + HalfSpace([-2.0, -2.0], 4.0), + HalfSpace([-2.0, -1.0], 2.5), + HalfSpace([-2.0, 0.0], 2.0), + HalfSpace([-2.0, 1.0], 2.5), + HalfSpace([-2.0, 2.0], 4.0), + HalfSpace([-1.0, -2.0], 2.5), + HalfSpace([-1.0, -1.0], 1.0), + HalfSpace([-1.0, 0.0], 0.5), + HalfSpace([-1.0, 1.0], 1.0), + HalfSpace([-1.0, 2.0], 2.5), + HalfSpace([0.0, -2.0], 2.0), + HalfSpace([0.0, -1.0], 0.5), + HalfSpace([0.0, 0.0], 0.0), + HalfSpace([0.0, 1.0], 0.5), + HalfSpace([0.0, 2.0], 2.0), + HalfSpace([1.0, -2.0], 2.5), + HalfSpace([1.0, -1.0], 1.0), + HalfSpace([1.0, 0.0], 0.5), + HalfSpace([1.0, 1.0], 1.0), + HalfSpace([1.0, 2.0], 2.5), + HalfSpace([2.0, -2.0], 4.0), + HalfSpace([2.0, -1.0], 2.5), + HalfSpace([2.0, 0.0], 2.0), + HalfSpace([2.0, 1.0], 2.5), + HalfSpace([2.0, 2.0], 4.0), + ] + square = polyhedron(reduce(intersect, ineq), lib) + sqverts = [-1 1 1 -1; -1 -1 1 1]/2 + @assert isvertsapprox(eachcol(sqverts), points(square)) + return volume(square) ≈ shoelace(sqverts) +end # issue 249 +function check_vol_issue249_1(lib) + poly = polyhedron(HalfSpace([-1.0, 0.0], 0.0) ∩ + HalfSpace([0.0, -1.0], 0.0) ∩ + HalfSpace([1.0, 0.0], 1.0) ∩ + HalfSpace([0.0, 1.0], 1.0) ∩ + HalfSpace([-0.2, -0.8], -0.0) ∩ + HalfSpace([0.6, 0.4], 0.6), lib) -poly = polyhedron(HalfSpace([-1.0, 0.0], 0.0) ∩ - HalfSpace([0.0, -1.0], 0.0) ∩ - HalfSpace([1.0, 0.0], 1.0) ∩ - HalfSpace([0.0, 1.0], 1.0) ∩ - HalfSpace([-0.2, -0.8], -0.0) ∩ - HalfSpace([0.6, 0.4], 0.6)) - -verts = [1/3 1 0 0; 1 0 0 1] -@test isvertsapprox(eachcol(verts), points(poly)) -@test volume(poly) ≈ shoelace(verts) + verts = [1/3 1 0 0; 1 0 0 1] + @assert isvertsapprox(eachcol(verts), points(poly)) + return volume(poly) ≈ shoelace(verts) +end -poly2 = polyhedron(# HalfSpace([-1.0, 0.0], 0.0) ∩ +function check_vol_issue249_2(lib) + poly2 = polyhedron(# HalfSpace([-1.0, 0.0], 0.0) ∩ HalfSpace([0.0, -1.0], 0.0) ∩ # Comment here to get correct area HalfSpace([1.0, 0.0], 1.0) ∩ HalfSpace([0.0, 1.0], 1.0) ∩ HalfSpace([-0.6, -0.4], -0.6) ∩ - HalfSpace([0.2, 0.8], 0.95)) -verts2 = [1/3 3/4 1 1; 1 1 15/16 0] -@test isvertsapprox(eachcol(verts2), points(poly2)) -@test volume(poly2) ≈ shoelace(verts2) + HalfSpace([0.2, 0.8], 0.95), lib) + verts2 = [1/3 3/4 1 1; 1 1 15/16 0] + @assert isvertsapprox(eachcol(verts2), points(poly2)) + return volume(poly2) ≈ shoelace(verts2) +end -lib=CDDLib.Library(:exact) -L3 = polyhedron(vrep([ - 0 0 0 0 0 0 - 0 0 1 0 0 0 - 0 1 0 0 0 0 - 0 1 1 0 0 1 - 1 0 0 0 0 0 - 1 0 1 0 1 0 - 1 1 0 1 0 0 - 1 1 1 1 1 1 +function check_vol_issue249_3(lib) + L3 = polyhedron(vrep([ + 0 0 0 0 0 0 + 0 0 1 0 0 0 + 0 1 0 0 0 0 + 0 1 1 0 0 1 + 1 0 0 0 0 0 + 1 0 1 0 1 0 + 1 1 0 1 0 0 + 1 1 1 1 1 1 ]),lib) -# reportedly solution taken from https://doi.org/10.1016/j.dam.2018.10.038 -@test volume(L3) == 1//180 - + # reportedly solution taken from https://doi.org/10.1016/j.dam.2018.10.038 + return volume(L3) == 1//180 +end # examples similar to Cohen-Hickey paper, table 1 # https://doi.org/10.1145/322139.322141 -s5 = polyhedron(vrep([ - 0 0 0 0 0 - 0 0 0 0 1 - 0 0 0 1 0 - 0 0 1 0 0 - 0 1 0 0 0 - 1 0 0 0 0 - ]), lib) -@test volume(s5) == 1//120 +function check_vol_simplex(n, lib) + s = polyhedron(vrep([i==j for i in 0:n, j in 1:n]), lib) + return volume(s) == 1//factorial(n) +end -flib = CDDLib.Library() -ϕ = (1 + √5)/2 -isoc = polyhedron(vrep([ - 0 1 ϕ - 0 1 -ϕ - 0 -1 ϕ - 0 -1 -ϕ - 1 ϕ 0 - 1 -ϕ 0 - -1 ϕ 0 - -1 -ϕ 0 - ϕ 0 1 - ϕ 0 -1 - -ϕ 0 1 - -ϕ 0 -1 - ]), flib) -# https://en.wikipedia.org/wiki/Regular_icosahedron -@test volume(isoc) ≈ (5/12)*(3+√5)*2^3 +function check_vol_isocahedron(lib) + ϕ = (1 + √5)/2 + isoc = polyhedron(vrep([ + 0 1 ϕ + 0 1 -ϕ + 0 -1 ϕ + 0 -1 -ϕ + 1 ϕ 0 + 1 -ϕ 0 + -1 ϕ 0 + -1 -ϕ 0 + ϕ 0 1 + ϕ 0 -1 + -ϕ 0 1 + -ϕ 0 -1 + ]), lib) + # https://en.wikipedia.org/wiki/Regular_icosahedron + return volume(isoc) ≈ (5/12)*(3+√5)*2^3 +end + +function check_vol_dodecahedron(lib) + ϕ = (1 + √5)/2 + ϕ² = ϕ^2 + dodec = polyhedron(vrep([ + ϕ ϕ ϕ + ϕ ϕ -ϕ + ϕ -ϕ ϕ + ϕ -ϕ -ϕ + -ϕ ϕ ϕ + -ϕ ϕ -ϕ + -ϕ -ϕ ϕ + -ϕ -ϕ -ϕ + 0 1 ϕ² + 0 1 -ϕ² + 0 -1 ϕ² + 0 -1 -ϕ² + 1 ϕ² 0 + 1 -ϕ² 0 + -1 ϕ² 0 + -1 -ϕ² 0 + ϕ² 0 1 + ϕ² 0 -1 + -ϕ² 0 1 + -ϕ² 0 -1 + ]), lib) + # https://en.wikipedia.org/wiki/Regular_dodecahedron + return volume(dodec) ≈ (1/4)*(15+7*√5)*2^3 +end -ϕ² = ϕ^2 -dodec = polyhedron(vrep([ - ϕ ϕ ϕ - ϕ ϕ -ϕ - ϕ -ϕ ϕ - ϕ -ϕ -ϕ - -ϕ ϕ ϕ - -ϕ ϕ -ϕ - -ϕ -ϕ ϕ - -ϕ -ϕ -ϕ - 0 1 ϕ² - 0 1 -ϕ² - 0 -1 ϕ² - 0 -1 -ϕ² - 1 ϕ² 0 - 1 -ϕ² 0 - -1 ϕ² 0 - -1 -ϕ² 0 - ϕ² 0 1 - ϕ² 0 -1 - -ϕ² 0 1 - -ϕ² 0 -1 - ]), flib) -# https://en.wikipedia.org/wiki/Regular_dodecahedron -@test volume(dodec) ≈ (1/4)*(15+7*√5)*2^3 +@testset "volumes" begin + lib_float = DefaultLibrary{Float64}() + lib_exact = DefaultLibrary{Int}() + @testset "simplex" for n in 1:5 + @test check_vol_simplex(n, lib_exact) + end + @test check_vol_issue285(lib_float) + @test check_vol_issue249_1(lib_float) + @test check_vol_issue249_2(lib_float) + @test check_vol_issue249_3(lib_exact) + @test check_vol_isocahedron(lib_float) + @test check_vol_dodecahedron(lib_float) +end