From 98c42a727988fc9639891b3295db44b81c6b45df Mon Sep 17 00:00:00 2001 From: lxvm Date: Fri, 29 Dec 2023 17:13:12 -0800 Subject: [PATCH 1/3] initial commit: new tests and temporary bugfix --- src/triangulation.jl | 2 +- test/runtests.jl | 1 + test/volume.jl | 147 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 149 insertions(+), 1 deletion(-) create mode 100644 test/volume.jl diff --git a/src/triangulation.jl b/src/triangulation.jl index a7335c2d..1dddc5c5 100644 --- a/src/triangulation.jl +++ b/src/triangulation.jl @@ -44,7 +44,7 @@ function triangulation_indices(p::Polyhedron) incident_idx = Dict(h => Set(incidentpointindices(p, h)) for h in h_idx) is_weak_adjacent = Dict((hi, hj) => !isempty(incident_idx[hi] ∩ incident_idx[hj]) for hi in h_idx for hj in h_idx) _triangulation(Δs, Δ, v_idx, h_idx, incident_idx, is_weak_adjacent, fulldim(p)) - return Δs + return unique(Δs) end function triangulation(p::Polyhedron) return map(Δ -> vrep(get.(p, Δ)), triangulation_indices(p)) diff --git a/test/runtests.jl b/test/runtests.jl index e164dcc7..0af2b70d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -38,3 +38,4 @@ for (arith, T) in (("floating point", Float64), ("exact", Rational{BigInt})) end include("../examples/test_examples.jl") +include("volume.jl") diff --git a/test/volume.jl b/test/volume.jl new file mode 100644 index 00000000..5dad3b50 --- /dev/null +++ b/test/volume.jl @@ -0,0 +1,147 @@ +using Test +using Polyhedra +using CDDLib + +# simple algorithm for calculating area of polygons +# requires vertices to be sorted (counter)clockwise +function shoelace(verts::AbstractMatrix{<:Real}) + @assert size(verts, 1) == 2 "shoelace only works for polygons" + xs = verts[begin,:] + ys = verts[end,:] + A = (ys[end]+ys[begin])*(xs[end]-xs[begin]) + for i in axes(verts,2)[begin:end-1] + A += (ys[i]+ys[i+1])*(xs[i]-xs[i+1]) + end + 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) + +# 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) + +# issue 249 + +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) + +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) + +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 + ]),lib) +# reportedly solution taken from https://doi.org/10.1016/j.dam.2018.10.038 +@test volume(L3) == 1//180 + +# 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 + +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 + +ϕ² = ϕ^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 From 76582d1ee26153c77035785491ed4add388878dc Mon Sep 17 00:00:00 2001 From: lxvm Date: Tue, 30 Jan 2024 23:21:45 -0500 Subject: [PATCH 2/3] add CDDLib as test dependency --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index 3b9d50da..475ef2ed 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [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" From a8671132afe8af7ebd517dfb31930991b11a6d02 Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 31 Jan 2024 09:21:51 -0500 Subject: [PATCH 3/3] apply suggestions --- test/Project.toml | 1 - test/volume.jl | 244 ++++++++++++++++++++++++---------------------- 2 files changed, 130 insertions(+), 115 deletions(-) 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