JuliaReach / LazySets.jl

Scalable symbolic-numeric set computations in Julia
https://juliareach.github.io/LazySets.jl/
Other
226 stars 32 forks source link

Elimination with mixed types #1952

Closed mforets closed 3 years ago

mforets commented 4 years ago

If M is a dense matrix and if P is a polyhedron whose constraints are defined by sparse arrays, elimination fails because Phrep = Polyhedra.hrep(y_eq_Mx, Ax_leq_b) gets mixed types.

For example,

using LazySets, Polyhedra, CDDLib, SparseArrays

p = HPolyhedron([LazySets.HalfSpace(sprand(5, 1.0), 1.0) for _ in 1:10]);

m = rand(5, 5)

linear_map(m, p, algorithm="elimination")

gives

julia> linear_map(m, p, algorithm="elimination")
 = = enter = =
typeof(M) = Array{Float64,2}
typeof(P) = HPolyhedron{Float64}
typeof(P.constraints[1]) = LazySets.HalfSpace{Float64,SparseVector{Float64,Int64}}
typeof(Ax_leq_b) = Array{Polyhedra.HalfSpace{Float64,SparseVector{Float64,Int64}},1}
typeof(y_eq_Mx) = Array{HyperPlane{Float64,Array{Float64,1}},1}

 = = end = =
ERROR: MethodError: no method matching Polyhedra.Intersection(::Int64, ::Array{HyperPlane{Float64,Array{Float64,1}},1}, ::Array{Polyhedra.HalfSpace{Float64,Spars
eVector{Float64,Int64}},1})
Closest candidates are:
  Polyhedra.Intersection(::Union{Int64, StaticArrays.Size}, ::Union{AbstractArray{HyperPlane{T,AT},1}, Polyhedra.AbstractRepIterator{#s17,HyperPlane{T,AT}} where
 #s17, Polyhedra.AllRepIterator{#s18,HyperPlane{T,AT},LinElemT,LRT,RT} where RT<:Polyhedra.AbstractRepIterator{#s18,HyperPlane{T,AT}} where LRT<:Polyhedra.Abstra
ctRepIterator{#s18,LinElemT} where LinElemT where #s18}, ::Union{AbstractArray{Polyhedra.HalfSpace{T,AT},1}, Polyhedra.AbstractRepIterator{#s17,Polyhedra.HalfSpa
ce{T,AT}} where #s17, Polyhedra.AllRepIterator{#s18,Polyhedra.HalfSpace{T,AT},LinElemT,LRT,RT} where RT<:Polyhedra.AbstractRepIterator{#s18,Polyhedra.HalfSpace{T
,AT}} where LRT<:Polyhedra.AbstractRepIterator{#s18,LinElemT} where LinElemT where #s18}) where {T, AT} at /home/mforets/.julia/dev/Polyhedra/src/vecrep.jl:68
Stacktrace:
 [1] #hrep#70(::Int64, ::typeof(hrep), ::Array{HyperPlane{Float64,Array{Float64,1}},1}, ::Array{Polyhedra.HalfSpace{Float64,SparseVector{Float64,Int64}},1}) at /
home/mforets/.julia/dev/Polyhedra/src/vecrep.jl:30
 [2] hrep(::Array{HyperPlane{Float64,Array{Float64,1}},1}, ::Array{Polyhedra.HalfSpace{Float64,SparseVector{Float64,Int64}},1}) at /home/mforets/.julia/dev/Polyh
edra/src/vecrep.jl:30
 [3] _linear_map_hrep(::Array{Float64,2}, ::HPolyhedron{Float64}, ::LazySets.LinearMapElimination{CDDLib.Library,BlockElimination}) at /home/mforets/.julia/dev/L
azySets/src/Interfaces/AbstractPolyhedron_functions.jl:736
 [4] _linear_map_hrep_helper(::Array{Float64,2}, ::HPolyhedron{Float64}, ::LazySets.LinearMapElimination{CDDLib.Library,BlockElimination}) at /home/mforets/.juli
a/dev/LazySets/src/Interfaces/AbstractPolyhedron_functions.jl:661
 [5] #linear_map#33(::String, ::Bool, ::Float64, ::Nothing, ::Nothing, ::Nothing, ::typeof(linear_map), ::Array{Float64,2}, ::HPolyhedron{Float64}) at /home/mfor
ets/.julia/dev/LazySets/src/Interfaces/AbstractPolyhedron_functions.jl:618
 [6] (::getfield(LazySets, Symbol("#kw##linear_map")))(::NamedTuple{(:algorithm,),Tuple{String}}, ::typeof(linear_map), ::Array{Float64,2}, ::HPolyhedron{Float64
}) at ./none:0
 [7] top-level scope at REPL[19]:1
mforets commented 4 years ago

The fix would be to force the same vector type in

    Ax_leq_b = [Polyhedra.HalfSpace(vcat(zeros(N, m), c.a), c.b) for c in constraints_list(P)]
    y_eq_Mx = [Polyhedra.HyperPlane(vcat(-id_m[i, :], M[i, :]), zero(N)) for i in 1:m]

This issue is basicallty the same as in "lift", see https://github.com/JuliaReach/LazySets.jl/issues/1942.

schillic commented 3 years ago

This works in master now.