Merge pull request #96 from JuliaPolyhedra/bl/liftedconcrete

Concrete type for liftedrep

Author: blegat

src/Polyhedra.jl

@@ -41,6 +41,8 @@ end
 similar_type(::Type{<:Point}, ::FullDim{N}, ::Type{T}) where {N, T} = Point{N,T}
 similar_type(::Type{<:Vec}, ::FullDim{N}, ::Type{T}) where {N, T} = Vec{N,T}
 
+emptymatrix(::Type{MT}, m, n) where {MT<:AbstractMatrix} = MT(m, n)
+emptymatrix(::Type{SparseMatrixCSC{T, Int}}, m, n) where T = spzeros(T, m, n)
 similar_type(::Type{<:Matrix}, ::Type{T}) where T = Matrix{T}
 similar_type(::Type{SparseMatrixCSC{S, I}}, ::Type{T}) where {S, I, T} = SparseMatrixCSC{T, I}
 arraytype(::Type{<:AbstractSparseArray{T}}) where T = SparseVector{T, Int}

src/doubledescription.jl

@@ -68,7 +68,7 @@ end
 
 function doubledescription(v::VRepresentation{N, T}) where {N, T}
     checkvconsistency(v)
-    lv = LiftedVRepresentation{N, T}(v)
+    lv = LiftedVRepresentation{N, T, Matrix{T}}(v)
     R = -lv.R
     vl = doubledescription(MixedMatHRep{N+1, T}(R, zeros(T, size(R, 1)), lv.linset))
     LiftedHRepresentation{N, T}(vl.R, vl.Rlinset)

src/liftedrep.jl

@@ -3,32 +3,36 @@ export LiftedHRepresentation, LiftedVRepresentation
 # H-Represenation
 
 # No copy since I do not modify anything and a copy is done when building a polyhedron
-mutable struct LiftedHRepresentation{N, T} <: MixedHRep{N, T}
+mutable struct LiftedHRepresentation{N, T, MT<:AbstractMatrix{T}} <: MixedHRep{N, T}
     # Ax >= 0, it is [b -A] * [z; x] where z = 1
-    A::AbstractMatrix{T}
+    A::MT
     linset::IntSet
 
-    function LiftedHRepresentation{N, T}(A::AbstractMatrix, linset::IntSet=IntSet()) where {N, T}
+    function LiftedHRepresentation{N, T, MT}(A::MT, linset::IntSet=IntSet()) where {N, T, MT}
         if !isempty(linset) && last(linset) > size(A, 1)
             error("The elements of linset should be between 1 and the number of rows of A")
         end
         if size(A, 2) != N+1
             error("dimension does not match")
         end
-        new{N, T}(A, linset)
+        new{N, T, MT}(A, linset)
     end
 end
 
-similar_type(::Type{<:LiftedHRepresentation}, ::FullDim{N}, ::Type{T}) where {N,T} = LiftedHRepresentation{N,T}
-arraytype(p::Union{LiftedHRepresentation{N, T}, Type{LiftedHRepresentation{N, T}}}) where {N, T} = Vector{T}
+similar_type(::Type{LiftedHRepresentation{M, S, MT}}, ::FullDim{N}, ::Type{T}) where {M, S, N, T, MT} = LiftedHRepresentation{N, T, similar_type(MT, T)}
+arraytype(p::Union{LiftedHRepresentation{N, T, MT}, Type{LiftedHRepresentation{N, T, MT}}}) where {N, T, MT} = arraytype(MT)
 
-LiftedHRepresentation(A::AbstractMatrix{T}, linset::IntSet=IntSet()) where {T <: Real} = LiftedHRepresentation{size(A,2)-1,T}(A, linset)
-LiftedHRepresentation(h::HRepresentation{N,T}) where {N,T} = LiftedHRepresentation{N,T}(h)
+LiftedHRepresentation{N, T}(A::AbstractMatrix{T}, linset::IntSet=IntSet()) where {N, T} = LiftedHRepresentation{N, T, typeof(A)}(A, linset)
+LiftedHRepresentation{N, T}(A::AbstractMatrix, linset::IntSet=IntSet()) where {N, T} = LiftedHRepresentation{N, T}(AbstractMatrix{T}(A), linset)
+LiftedHRepresentation(A::AbstractMatrix{T}, linset::IntSet=IntSet()) where T = LiftedHRepresentation{size(A, 2) - 1, T}(A, linset)
+function LiftedHRepresentation(h::HRepresentation{N, T}) where {N, T}
+    LiftedHRepresentation{N, T, arraytype(h) <: AbstractSparseVector ? SparseMatrixCSC{T, Int} : Matrix{T}}(h)
+end
 
-function LiftedHRepresentation{N, T}(hyperplanes::ElemIt{<:HyperPlane{N, T}}, halfspaces::ElemIt{<:HalfSpace{N, T}}) where {N, T}
+function LiftedHRepresentation{N, T, MT}(hyperplanes::ElemIt{<:HyperPlane{N, T}}, halfspaces::ElemIt{<:HalfSpace{N, T}}) where {N, T, MT}
     nhyperplane = length(hyperplanes)
     nhrep = nhyperplane + length(halfspaces)
-    A = Matrix{T}(nhrep, N+1)
+    A = emptymatrix(MT, nhrep, N+1)
     linset = IntSet(1:nhyperplane)
     for (i, h) in enumerate(hyperplanes)
         A[i,2:end] = -h.a
@@ -49,34 +53,36 @@ Base.get(hrep::LiftedHRepresentation{N, T}, idx::HIndex{N, T}) where {N, T} = va
 
 # V-Representation
 
-mutable struct LiftedVRepresentation{N,T} <: MixedVRep{N,T}
-    R::AbstractMatrix{T} # each row is a vertex if the first element is 1 and a ray otherwise
+mutable struct LiftedVRepresentation{N, T, MT<:AbstractMatrix{T}} <: MixedVRep{N, T}
+    R::MT # each row is a vertex if the first element is 1 and a ray otherwise
     linset::IntSet
 
-    function LiftedVRepresentation{N, T}(R::AbstractMatrix, linset::IntSet=IntSet([])) where {N, T}
+    function LiftedVRepresentation{N, T, MT}(R::MT, linset::IntSet=IntSet([])) where {N, T, MT}
         if length(R) > 0 && size(R, 2) != N+1
             error("dimension does not match")
         end
         if !isempty(linset) && last(linset) > size(R, 1)
             error("The elements of linset should be between 1 and the number of rows of R")
         end
-        new{N, T}(R, linset)
+        new{N, T, MT}(R, linset)
     end
 end
 
-similar_type(::Type{<:LiftedVRepresentation}, ::FullDim{N}, ::Type{T}) where {N,T} = LiftedVRepresentation{N,T}
-arraytype(p::Union{LiftedVRepresentation{N, T}, Type{LiftedVRepresentation{N, T}}}) where {N, T} = Vector{T}
+similar_type(::Type{LiftedVRepresentation{M, S, MT}}, ::FullDim{N}, ::Type{T}) where {M, S, N, T, MT} = LiftedVRepresentation{N, T, similar_type(MT, T)}
+arraytype(p::Union{LiftedVRepresentation{N, T, MT}, Type{LiftedVRepresentation{N, T, MT}}}) where {N, T, MT} = arraytype(MT)
 
-LiftedVRepresentation(R::AbstractMatrix{T}, linset::IntSet=IntSet()) where {T <: Real} = LiftedVRepresentation{size(R,2)-1,T}(R, linset)
-LiftedVRepresentation(v::VRepresentation{N,T}) where {N,T} = LiftedVRepresentation{N,T}(v)
+LiftedVRepresentation{N, T}(R::AbstractMatrix{T}, linset::IntSet=IntSet()) where {N, T} = LiftedVRepresentation{N, T, typeof(R)}(R, linset)
+LiftedVRepresentation{N, T}(R::AbstractMatrix, linset::IntSet=IntSet()) where {N, T} = LiftedVRepresentation{N, T}(AbstractMatrix{T}(R), linset)
+LiftedVRepresentation(R::AbstractMatrix{T}, linset::IntSet=IntSet()) where T = LiftedVRepresentation{size(R, 2) - 1, T}(R, linset)
+LiftedVRepresentation(v::VRepresentation{N,T}) where {N, T} = LiftedVRepresentation{N, T, Matrix{T}}(v)
 
-function LiftedVRepresentation{N, T}(vits::VIt{N, T}...) where {N, T}
+function LiftedVRepresentation{N, T, MT}(vits::VIt{N, T}...) where {N, T, MT}
     points, lines, rays = fillvits(FullDim{N}(), vits...)
     npoint = length(points)
     nline = length(lines)
     nray = length(rays)
     nvrep = npoint + nline + nray
-    R = Matrix{T}(nvrep, N+1)
+    R = emptymatrix(MT, nvrep, N+1)
     linset = IntSet()
     function _fill(offset, z, ps)
         for (i, p) in enumerate(ps)

src/lphrep.jl

@@ -100,7 +100,7 @@ LPHRepresentation(rep::HRep{N, T}) where {N, T} = LPHRepresentation{N, T, _matty
 function LPHRepresentation{N, T, MT}(hyperplanes::ElemIt{<:HyperPlane{N, T}}, halfspaces::ElemIt{<:HalfSpace{N, T}}) where {N, T, MT}
     nhyperplane = length(hyperplanes)
     nhrep = nhyperplane + length(halfspaces)
-    A = MT <: AbstractSparseArray ? spzeros(eltype(MT), nhrep, N) : MT(nhrep, N)
+    A = emptymatrix(MT, nhrep, N)
     lb = Vector{T}(nhrep)
     ub = Vector{T}(nhrep)
     MPB.HighLevelInterface.warn_no_inf(T)

src/matrep.jl

@@ -48,7 +48,7 @@ MixedMatHRep(h::HRep{N,T}) where {N,T} = MixedMatHRep{N,T}(h)
 function MixedMatHRep{N, T, MT}(hyperplanes::HyperPlaneIt{N, T}, halfspaces::HalfSpaceIt{N, T}) where {N, T, MT}
     nhyperplane = length(hyperplanes)
     nhrep = nhyperplane + length(halfspaces)
-    A = MT <: AbstractSparseArray ? spzeros(eltype(MT), nhrep, N) : MT(nhrep, N)
+    A = emptymatrix(MT, nhrep, N)
     b = Vector{T}(nhrep)
     linset = IntSet(1:nhyperplane)
     for (i, h) in enumerate(hyperplanes)
@@ -117,8 +117,8 @@ function MixedMatVRep{N, T, MT}(vits::VIt{N, T}...) where {N, T, MT}
     npoint = length(points)
     nline = length(lines)
     nray = length(rays)
-    V = Matrix{T}(npoint, N)
-    R = Matrix{T}(nline + nray, N)
+    V = emptymatrix(MT, npoint, N)
+    R = emptymatrix(MT, nline + nray, N)
     Rlinset = IntSet()
     function _fill!(M, linset, offset, ps)
         for (i, p) in enumerate(ps)

test/doubledescription.jl

@@ -124,14 +124,16 @@
                 @testset "Exact" begin
                     v = conichull([1, 0],
                                   [0, 1])
-                    h = @inferred doubledescription(v)
+                    h = doubledescription(v)
+                    #h = @inferred doubledescription(v)
                     @test !hashyperplanes(h)
                     @test collect(halfspaces(h)) == [HalfSpace([0, -1], 0), HalfSpace([-1, 0], 0), HalfSpace([0, 0], 1)] # FIXME get rid of (0, 0) 1
                 end
                 @testset "Numerical" begin
                     v = conichull([1., 0.],
                                   [0., 1.])
-                    h = @inferred doubledescription(v)
+                    h = doubledescription(v)
+                    #h = @inferred doubledescription(v)
                     @test !hashyperplanes(h)
                     @test collect(halfspaces(h)) == [HalfSpace([0, -1], 0), HalfSpace([-1, 0], 0), HalfSpace([0, 0], 1)] # FIXME get rid of (0, 0) 1
                 end

test/representation.jl

@@ -57,8 +57,8 @@ Polyhedra.@subrepelem InconsistentVRep Ray rays
         @test ine.A !== A
         #@test linset(ine) == ls
         @test ine.linset !== ls
-        @test Polyhedra.similar_type(LiftedHRepresentation{2, Int}, Float64) == LiftedHRepresentation{2, Float64}
-        @test Polyhedra.similar_type(LiftedHRepresentation{2, Int}, FullDim{3}(), Float64) == LiftedHRepresentation{3, Float64}
+        @test Polyhedra.similar_type(LiftedHRepresentation{2, Int, Matrix{Int}}, Float64) == LiftedHRepresentation{2, Float64, Matrix{Float64}}
+        @test Polyhedra.similar_type(LiftedHRepresentation{2, Int, SparseMatrixCSC{Int, Int}}, FullDim{3}(), Float64) == LiftedHRepresentation{3, Float64, SparseMatrixCSC{Float64, Int}}
 
         A2 = [1 1; -1 0; 0 -1]
         b2 = [1, 0, 0]
@@ -79,8 +79,8 @@ Polyhedra.@subrepelem InconsistentVRep Ray rays
         #@test linset(ext) == Vlinset
         @test ext.linset !== Vlinset
 
-        @test Polyhedra.similar_type(LiftedVRepresentation{2, Int}, Float64) == LiftedVRepresentation{2, Float64}
-        @test Polyhedra.similar_type(LiftedVRepresentation{2, Int}, FullDim{3}(), Float64) == LiftedVRepresentation{3, Float64}
+        @test Polyhedra.similar_type(LiftedVRepresentation{2, Int, SparseMatrixCSC{Int, Int}}, Float64) == LiftedVRepresentation{2, Float64, SparseMatrixCSC{Float64, Int}}
+        @test Polyhedra.similar_type(LiftedVRepresentation{2, Int, Matrix{Int}}, FullDim{3}(), Float64) == LiftedVRepresentation{3, Float64, Matrix{Float64}}
     end
 
     @testset "eltype for some iterators is incorrect #7" begin
@@ -201,7 +201,7 @@ Polyhedra.@subrepelem InconsistentVRep Ray rays
         N = 5
         M = 10
         T = Int64
-        reps = [MixedMatHRep{N, T, Matrix{T}}, MixedMatVRep{N, T, Matrix{T}}, LiftedHRepresentation{N, T}, LiftedVRepresentation{N, T}]
+        reps = [MixedMatHRep{N, T, Matrix{T}}, MixedMatVRep{N, T, Matrix{T}}, LiftedHRepresentation{N, T, Matrix{T}}, LiftedVRepresentation{N, T, Matrix{T}}]
         for rep in reps
             changedrep = Polyhedra.similar_type(rep, FullDim{M}())
             @test fulldim(changedrep) == M
@@ -299,7 +299,7 @@ Polyhedra.@subrepelem InconsistentVRep Ray rays
     @testset "V-consistency with iterator constructor" begin
         T = Int
         AT = Vector{Int}
-        for VRepType in (Polyhedra.LiftedVRepresentation{2, T},
+        for VRepType in (Polyhedra.LiftedVRepresentation{2, T, Matrix{T}},
                          Polyhedra.MixedMatVRep{2, T, Matrix{T}},
                          Polyhedra.Hull{2, T, AT})
             @test_throws ErrorException VRepType(AT[], [Line([1, 2])])