From ef60e2e17b6134e417651f1382319c4b40665376 Mon Sep 17 00:00:00 2001 From: Unknown Date: Sun, 2 Jun 2019 01:46:13 +1000 Subject: [PATCH 01/15] Added overapproximation of Lazy linear map of cartesian product array --- src/Approximations/overapproximate.jl | 55 +++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 9f54aa77f3..9466c5b39a 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -516,3 +516,58 @@ function overapproximate(cap::Intersection{N, return overapproximate(swap(cap), dir; kwargs...) end +""" + overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, <:LazySet{N}}}, + dir::Type{<:AbstractDirections})::CartesianProductArray{N, HPolytope{N}} where {N} + +Overapproximating a lazy linear map of cartesian product array with template directions for each block. + +### Input + +- `lm` -- lazy linear map of cartesian product array +- `dir` -- (concrete) direction representation + +### Output + +An `CartesianProductArray` with overapproximation of the each block with the directions from `dir`. +""" +function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, <:LazySet{N}}}, + dir::Type{<:AbstractDirections})::CartesianProductArray{N, HPolytope{N}} where {N} + + if !(size(lm.M, 2) == dim(lm.X)) + error("Matrix needs to be commensurate with the cartesian product") + end + + cp = lm.X + M = lm.M + + array = Vector{HPolytope{N}} + sizehint!(array,length(cp.array)) + col_st_index, col_end_ind = 0, 0 + + for bi in 1:length(cp.array) + col_st_index = col_end_ind + 1 + n = dim(cp.array[bi]) + col_end_ind += n + + start_ind, end_ind = 1, n + matrices = Vector{Matrix{N}}() + + while end_ind <= size(M, 2) + push!(matrices, M[col_st_index : col_end_ind, start_ind : end_ind]) + start_ind = end_ind + 1 + end_ind += n + end + + v_block = MinkowskiSumArray() + + for m in matrices + push!(v_block.array,LinearMap(m, cp.array[bi])) + end + + push!(array, Approximations.overapproximate(v_block, dir{N}(dim(v_block)))) + end + + result = CartesianProductArray{N, HPolytope{N}}(array) + return result +end From ecabee56473f9de522434a8d71c915fda6a50a96 Mon Sep 17 00:00:00 2001 From: Unknown Date: Tue, 4 Jun 2019 14:45:14 +1000 Subject: [PATCH 02/15] updated overaproximation of decomposed linear map according to Christian comments --- src/Approximations/overapproximate.jl | 77 +++++++++++++++++++-------- 1 file changed, 56 insertions(+), 21 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 9466c5b39a..f9384630f0 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -525,14 +525,15 @@ Overapproximating a lazy linear map of cartesian product array with template dir ### Input - `lm` -- lazy linear map of cartesian product array -- `dir` -- (concrete) direction representation +- Type{<:CartesianProductArray} -- type to specify decomposed overapproximation +- `oa` -- (concrete) direction representation ### Output An `CartesianProductArray` with overapproximation of the each block with the directions from `dir`. """ -function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, <:LazySet{N}}}, - dir::Type{<:AbstractDirections})::CartesianProductArray{N, HPolytope{N}} where {N} +function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, + <:LazySet{N}}}, ::Type{<:CartesianProductArray}, oa::Union{Type{<:AbstractDirections},Type{<:Hyperrectangle}, Type{<:Interval}}=Hyperrectangle) where {N} if !(size(lm.M, 2) == dim(lm.X)) error("Matrix needs to be commensurate with the cartesian product") @@ -541,33 +542,67 @@ function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, <:LazySet{N cp = lm.X M = lm.M - array = Vector{HPolytope{N}} - sizehint!(array,length(cp.array)) col_st_index, col_end_ind = 0, 0 - for bi in 1:length(cp.array) - col_st_index = col_end_ind + 1 - n = dim(cp.array[bi]) - col_end_ind += n + return decomposed_overapproximation(cp, M, col_st_index, col_end_ind, oa) +end - start_ind, end_ind = 1, n - matrices = Vector{Matrix{N}}() +#template directions +function decomposed_overapproximation(cp::CartesianProductArray{N,<:LazySet{N}}, M::AbstractArray{N}, + col_st_index::Int, col_end_ind::Int, oa::Type{<:AbstractDirections}) where {N} - while end_ind <= size(M, 2) - push!(matrices, M[col_st_index : col_end_ind, start_ind : end_ind]) - start_ind = end_ind + 1 - end_ind += n - end + array = Vector{HPolytope{N}}() - v_block = MinkowskiSumArray() + sizehint!(array,length(cp.array)) - for m in matrices - push!(v_block.array,LinearMap(m, cp.array[bi])) - end + for bi in cp.array + v_block, col_st_index, col_end_ind = overapproximate_block(bi, col_st_index, col_end_ind, M) - push!(array, Approximations.overapproximate(v_block, dir{N}(dim(v_block)))) + push!(array, Approximations.overapproximate(v_block, oa{N}(dim(v_block)))) end result = CartesianProductArray{N, HPolytope{N}}(array) return result end + +#Interval Hull +function decomposed_overapproximation(cp::CartesianProductArray{N,<:LazySet{N}}, M::AbstractArray{N}, + col_st_index::Int, col_end_ind::Int, oa::Union{Type{<:Hyperrectangle}, Type{<:Interval}}) where {N} + + array = Vector{oa{N}}() + + sizehint!(array,length(cp.array)) + + for bi in cp.array + v_block, col_st_index, col_end_ind = overapproximate_block(bi, col_st_index, col_end_ind, M) + + push!(array, Approximations.overapproximate(v_block, oa)) + end + + result = CartesianProductArray{N, oa{N}}(array) + return result +end + +#overapproximation of linear map to each block +function overapproximate_lm_block(bi:: LazySet{N}, col_st_index::Int, col_end_ind::Int, M::AbstractArray{N}) where {N} + col_st_index = col_end_ind + 1 + n = dim(bi) + col_end_ind += n + + start_ind, end_ind = 1, n + matrices = Vector{Matrix{N}}() + + while end_ind <= size(M, 2) + push!(matrices, M[col_st_index : col_end_ind, start_ind : end_ind]) + start_ind = end_ind + 1 + end_ind += n + end + + v_block = MinkowskiSumArray() + + for m in matrices + push!(v_block.array,LinearMap(m, bi)) + end + + return v_block, col_st_index, col_end_ind +end From b98b789116203d1b17f12e039b4f74333ca97aa0 Mon Sep 17 00:00:00 2001 From: Unknown Date: Tue, 4 Jun 2019 16:31:30 +1000 Subject: [PATCH 03/15] Added using @view --- src/Approximations/overapproximate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index f9384630f0..e945117c1f 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -593,7 +593,7 @@ function overapproximate_lm_block(bi:: LazySet{N}, col_st_index::Int, col_end_in matrices = Vector{Matrix{N}}() while end_ind <= size(M, 2) - push!(matrices, M[col_st_index : col_end_ind, start_ind : end_ind]) + push!(matrices, @view M[col_st_index : col_end_ind, start_ind : end_ind]) start_ind = end_ind + 1 end_ind += n end From cf01ae601496db2133f9b7b502bac1e4e533dbc2 Mon Sep 17 00:00:00 2001 From: Unknown Date: Tue, 4 Jun 2019 18:50:19 +1000 Subject: [PATCH 04/15] updated according to comments --- src/Approximations/overapproximate.jl | 67 +++++++++++---------------- 1 file changed, 26 insertions(+), 41 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index e945117c1f..5fd92e1345 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -517,8 +517,8 @@ function overapproximate(cap::Intersection{N, end """ - overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, <:LazySet{N}}}, - dir::Type{<:AbstractDirections})::CartesianProductArray{N, HPolytope{N}} where {N} + overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, + <:LazySet{N}}}, ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} Overapproximating a lazy linear map of cartesian product array with template directions for each block. @@ -533,76 +533,61 @@ Overapproximating a lazy linear map of cartesian product array with template dir An `CartesianProductArray` with overapproximation of the each block with the directions from `dir`. """ function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, - <:LazySet{N}}}, ::Type{<:CartesianProductArray}, oa::Union{Type{<:AbstractDirections},Type{<:Hyperrectangle}, Type{<:Interval}}=Hyperrectangle) where {N} + <:LazySet{N}}}, ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} - if !(size(lm.M, 2) == dim(lm.X)) - error("Matrix needs to be commensurate with the cartesian product") - end + @assert !(size(lm.M, 2) == dim(lm.X)) "matrix needs to be commensurate with the cartesian product" cp = lm.X M = lm.M - col_st_index, col_end_ind = 0, 0 + col_end_ind = 0 return decomposed_overapproximation(cp, M, col_st_index, col_end_ind, oa) end #template directions -function decomposed_overapproximation(cp::CartesianProductArray{N,<:LazySet{N}}, M::AbstractArray{N}, - col_st_index::Int, col_end_ind::Int, oa::Type{<:AbstractDirections}) where {N} - - array = Vector{HPolytope{N}}() +function decomposed_overapproximation(cp::CartesianProductArray{N,<:LazySet{N}}, M::AbstractMatrix{N}, + col_end_ind::Int, oa) where {N} + array = allocate_result(oa, N) sizehint!(array,length(cp.array)) for bi in cp.array - v_block, col_st_index, col_end_ind = overapproximate_block(bi, col_st_index, col_end_ind, M) - - push!(array, Approximations.overapproximate(v_block, oa{N}(dim(v_block)))) + v_block, col_end_ind = overapproximate_block(bi, col_end_ind, M) + push!(array, Approximations.overapproximate(v_block, oa)) end - result = CartesianProductArray{N, HPolytope{N}}(array) + result = CartesianProductArray{N}(array) return result end -#Interval Hull -function decomposed_overapproximation(cp::CartesianProductArray{N,<:LazySet{N}}, M::AbstractArray{N}, - col_st_index::Int, col_end_ind::Int, oa::Union{Type{<:Hyperrectangle}, Type{<:Interval}}) where {N} - - array = Vector{oa{N}}() - - sizehint!(array,length(cp.array)) - - for bi in cp.array - v_block, col_st_index, col_end_ind = overapproximate_block(bi, col_st_index, col_end_ind, M) +function allocate_result(oa, N) + return Vector{LazySet{N}}() +end - push!(array, Approximations.overapproximate(v_block, oa)) - end +function allocate_result(oa::Type{<:LazySet}, N) + return Vector{oa{N}}() +end - result = CartesianProductArray{N, oa{N}}(array) - return result +function allocate_result(oa::Type{<:AbstractDirections}, N) + return Vector{HPolytope{N}}() end #overapproximation of linear map to each block -function overapproximate_lm_block(bi:: LazySet{N}, col_st_index::Int, col_end_ind::Int, M::AbstractArray{N}) where {N} +function overapproximate_lm_block(bi::LazySet{N}, col_end_ind::Int, M::AbstractMatrix{N}) where {N} col_st_index = col_end_ind + 1 n = dim(bi) col_end_ind += n - start_ind, end_ind = 1, n + row_start_ind, row_end_ind = 1, n matrices = Vector{Matrix{N}}() - while end_ind <= size(M, 2) - push!(matrices, @view M[col_st_index : col_end_ind, start_ind : end_ind]) - start_ind = end_ind + 1 - end_ind += n - end - v_block = MinkowskiSumArray() - - for m in matrices - push!(v_block.array,LinearMap(m, bi)) + while row_end_ind <= size(M, 2) + push!(v_block.array, LinearMap(M[col_st_index : col_end_ind, row_start_ind : row_end_ind], bi)) + row_start_ind = row_end_ind + 1 + row_end_ind += n end - return v_block, col_st_index, col_end_ind + return v_block, col_end_ind end From 5b97aec6e0a09f737cc8aa0010ab7604c8556d35 Mon Sep 17 00:00:00 2001 From: Unknown Date: Tue, 4 Jun 2019 18:57:12 +1000 Subject: [PATCH 05/15] small fixes --- src/Approximations/overapproximate.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 5fd92e1345..8101813f2e 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -535,14 +535,14 @@ An `CartesianProductArray` with overapproximation of the each block with the dir function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, <:LazySet{N}}}, ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} - @assert !(size(lm.M, 2) == dim(lm.X)) "matrix needs to be commensurate with the cartesian product" + @assert size(lm.M, 2) == dim(lm.X) "matrix needs to be commensurate with the cartesian product" cp = lm.X M = lm.M col_end_ind = 0 - return decomposed_overapproximation(cp, M, col_st_index, col_end_ind, oa) + return decomposed_overapproximation(cp, M, col_end_ind, oa) end #template directions @@ -553,11 +553,11 @@ function decomposed_overapproximation(cp::CartesianProductArray{N,<:LazySet{N}}, sizehint!(array,length(cp.array)) for bi in cp.array - v_block, col_end_ind = overapproximate_block(bi, col_end_ind, M) + v_block, col_end_ind = overapproximate_lm_block(bi, col_end_ind, M) push!(array, Approximations.overapproximate(v_block, oa)) end - result = CartesianProductArray{N}(array) + result = CartesianProductArray(array) return result end From 80ab0cb705d44df1d63b4605188cd0d5488b9bde Mon Sep 17 00:00:00 2001 From: Unknown Date: Tue, 4 Jun 2019 23:43:19 +1000 Subject: [PATCH 06/15] Correct minkowski sum addition --- src/Approximations/overapproximate.jl | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 8101813f2e..e93b9bfa60 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -549,15 +549,22 @@ end function decomposed_overapproximation(cp::CartesianProductArray{N,<:LazySet{N}}, M::AbstractMatrix{N}, col_end_ind::Int, oa) where {N} array = allocate_result(oa, N) - sizehint!(array,length(cp.array)) + min_array = Dict{Int,MinkowskiSumArray{N}}() + for i in 1:size(M, 2) + min_array[i] = MinkowskiSumArray() + end + for bi in cp.array - v_block, col_end_ind = overapproximate_lm_block(bi, col_end_ind, M) - push!(array, Approximations.overapproximate(v_block, oa)) + lms, col_end_ind = overapproximate_lm_block(bi, col_end_ind, M) + for i in 1:length(lms) + push!(min_array[i].array, lms[i]) + end end - result = CartesianProductArray(array) + result = CartesianProductArray([Approximations.overapproximate(min_array[key], oa) + for key in sort(collect(keys(min_array)))]) return result end @@ -578,16 +585,13 @@ function overapproximate_lm_block(bi::LazySet{N}, col_end_ind::Int, M::AbstractM col_st_index = col_end_ind + 1 n = dim(bi) col_end_ind += n - row_start_ind, row_end_ind = 1, n - matrices = Vector{Matrix{N}}() - - v_block = MinkowskiSumArray() + lms = Vector{LazySet}() while row_end_ind <= size(M, 2) - push!(v_block.array, LinearMap(M[col_st_index : col_end_ind, row_start_ind : row_end_ind], bi)) + push!(lms, LinearMap(M[row_start_ind : row_end_ind, col_st_index : col_end_ind], bi)) row_start_ind = row_end_ind + 1 row_end_ind += n end - return v_block, col_end_ind + return lms, col_end_ind end From c7fc98d5bf66c0d8f60737437e0f64d296fee560 Mon Sep 17 00:00:00 2001 From: Unknown Date: Wed, 5 Jun 2019 01:10:10 +1000 Subject: [PATCH 07/15] cleaned up the code --- src/Approximations/overapproximate.jl | 49 +++++++++++---------------- 1 file changed, 20 insertions(+), 29 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index e93b9bfa60..70d599f2d5 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -546,26 +546,33 @@ function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, end #template directions -function decomposed_overapproximation(cp::CartesianProductArray{N,<:LazySet{N}}, M::AbstractMatrix{N}, - col_end_ind::Int, oa) where {N} +function decomposed_overapproximation(cpa::CartesianProductArray{N,<:LazySet{N}}, + M::AbstractMatrix{N}, + col_end_ind::Int, oa) where {N} array = allocate_result(oa, N) - sizehint!(array,length(cp.array)) + sizehint!(array,length(cpa.array)) - min_array = Dict{Int,MinkowskiSumArray{N}}() for i in 1:size(M, 2) - min_array[i] = MinkowskiSumArray() + ms = overapproximate_row_blocks(cpa, M, i) + push!(array, overapproximate(ms, oa)) end - for bi in cp.array - lms, col_end_ind = overapproximate_lm_block(bi, col_end_ind, M) - for i in 1:length(lms) - push!(min_array[i].array, lms[i]) - end + result = CartesianProductArray(array) + return result +end + +#overapproximation of linear map to each block +function overapproximate_row_blocks(cpa::CartesianProductArray{N,<:LazySet{N}}, + M::AbstractMatrix{N}, i::Int) where {N} + col_start_ind, col_end_ind = 1, 0 + h_min_sum = MinkowskiSumArray() + for bi in cpa.array + col_end_ind += dim(bi) + push!(h_min_sum.array, LinearMap(M[i : i, col_start_ind : col_end_ind], bi)) + col_start_ind = col_end_ind + 1 end - result = CartesianProductArray([Approximations.overapproximate(min_array[key], oa) - for key in sort(collect(keys(min_array)))]) - return result + return h_min_sum end function allocate_result(oa, N) @@ -579,19 +586,3 @@ end function allocate_result(oa::Type{<:AbstractDirections}, N) return Vector{HPolytope{N}}() end - -#overapproximation of linear map to each block -function overapproximate_lm_block(bi::LazySet{N}, col_end_ind::Int, M::AbstractMatrix{N}) where {N} - col_st_index = col_end_ind + 1 - n = dim(bi) - col_end_ind += n - row_start_ind, row_end_ind = 1, n - lms = Vector{LazySet}() - while row_end_ind <= size(M, 2) - push!(lms, LinearMap(M[row_start_ind : row_end_ind, col_st_index : col_end_ind], bi)) - row_start_ind = row_end_ind + 1 - row_end_ind += n - end - - return lms, col_end_ind -end From 8742462bcce317bed81cacdceab16658c3b03a0e Mon Sep 17 00:00:00 2001 From: Unknown Date: Wed, 5 Jun 2019 01:23:36 +1000 Subject: [PATCH 08/15] added tests --- test/unit_overapproximate.jl | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/test/unit_overapproximate.jl b/test/unit_overapproximate.jl index e318ea8273..ac1a335822 100644 --- a/test/unit_overapproximate.jl +++ b/test/unit_overapproximate.jl @@ -47,7 +47,7 @@ for N in [Float64, Rational{Int}, Float32] h2 = Hyperrectangle(N[2.5, 4.5], N[1/2, 1/2]) H = overapproximate(h1 × h2, Hyperrectangle) # defaults to convert method @test low(H) == N[0, 2, 4] && high(H) == N[1, 3, 5] - + # overapproximation of the lazy linear map of a hyperrectangular set H = Hyperrectangle(N[0, 0], N[1/2, 1]) M = Diagonal(N[2, 2]) @@ -129,4 +129,22 @@ for N in [Float64, Float32] Y_zonotope = overapproximate(Y, Zonotope) # overapproximate with a zonotope @test Y_polygon ⊆ Y_zonotope @test !(Y_zonotope ⊆ Y_polygon) + + #Decomposed approximation of lazy linear map of CartesianProductArray + i1 = Interval([0, 1]) + h = Hyperrectangle(low=[3., 3.], high=[5., 5.]) + M = [1. 2. 3.; 4. 5. 6.; 7. 8. 9.] + cpa = CartesianProductArray([i1, h]) + lm = M * cpa + + oa = overapproximate(lm, Hyperrectangle) + oa_box = overapproximate(lm, Approximations.BoxDirections) + d_oa_d_hp = overapproximate(lm, CartesianProductArray) + d_oa_d_box = overapproximate(lm, CartesianProductArray, Approximations.BoxDirections) + oa_d_hp = overapproximate(d_oa_d_hp) + oa_d_box = overapproximate(d_oa_d_box, Approximations.BoxDirections) + + @test oa == oa_d_hp + @test oa_box == oa_d_box + end From 54d73fdd2e440f41a81df078c310b1717874e3df Mon Sep 17 00:00:00 2001 From: Unknown Date: Wed, 5 Jun 2019 04:11:41 +1000 Subject: [PATCH 09/15] updated linear map to use correct rows and updated tests --- src/Approximations/overapproximate.jl | 32 +++++++++------------ test/unit_overapproximate.jl | 41 ++++++++++++++++----------- 2 files changed, 38 insertions(+), 35 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 70d599f2d5..25a1012d89 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -520,55 +520,51 @@ end overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, <:LazySet{N}}}, ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} -Overapproximating a lazy linear map of cartesian product array with template directions for each block. +Decompose a lazy linear map of a cartesian product array. ### Input - `lm` -- lazy linear map of cartesian product array -- Type{<:CartesianProductArray} -- type to specify decomposed overapproximation -- `oa` -- (concrete) direction representation +- `CartesianProductArray` -- type for dispatch +- `oa` -- approximation option for decomposition ### Output -An `CartesianProductArray` with overapproximation of the each block with the directions from `dir`. +A `CartesianProductArray` representing the decomposed linear map. """ function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, <:LazySet{N}}}, ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} @assert size(lm.M, 2) == dim(lm.X) "matrix needs to be commensurate with the cartesian product" - cp = lm.X + cpa = lm.X M = lm.M col_end_ind = 0 - return decomposed_overapproximation(cp, M, col_end_ind, oa) -end - -#template directions -function decomposed_overapproximation(cpa::CartesianProductArray{N,<:LazySet{N}}, - M::AbstractMatrix{N}, - col_end_ind::Int, oa) where {N} array = allocate_result(oa, N) sizehint!(array,length(cpa.array)) - for i in 1:size(M, 2) - ms = overapproximate_row_blocks(cpa, M, i) + i, n = 1, 0 + for bi in cpa.array + n += dim(bi) + ms = blocks_linear_map(cpa, M, i, n) push!(array, overapproximate(ms, oa)) + i += n end result = CartesianProductArray(array) return result end -#overapproximation of linear map to each block -function overapproximate_row_blocks(cpa::CartesianProductArray{N,<:LazySet{N}}, - M::AbstractMatrix{N}, i::Int) where {N} + +function blocks_linear_map(cpa::CartesianProductArray{N,<:LazySet{N}}, + M::AbstractMatrix{N}, row_start_ind::Int, row_end_ind::Int) where {N} col_start_ind, col_end_ind = 1, 0 h_min_sum = MinkowskiSumArray() for bi in cpa.array col_end_ind += dim(bi) - push!(h_min_sum.array, LinearMap(M[i : i, col_start_ind : col_end_ind], bi)) + push!(h_min_sum.array, LinearMap(M[row_start_ind : row_end_ind, col_start_ind : col_end_ind], bi)) col_start_ind = col_end_ind + 1 end diff --git a/test/unit_overapproximate.jl b/test/unit_overapproximate.jl index ac1a335822..4e335bbe46 100644 --- a/test/unit_overapproximate.jl +++ b/test/unit_overapproximate.jl @@ -53,6 +53,30 @@ for N in [Float64, Rational{Int}, Float32] M = Diagonal(N[2, 2]) OA = overapproximate(M*H, Hyperrectangle) @test OA isa Hyperrectangle && OA.center == N[0, 0] && OA.radius == N[1, 2] + + #overapproximation of Minkowski sum of linear maps for each block in the row block + i1 = Interval(N[0, 1]) + h = Hyperrectangle(low=N[3, 4], high=N[5, 7]) + M = N[1 2 3; 4 5 6; 7 8 9] + cpa = CartesianProductArray([i1, h]) + lm = M * cpa + + oa = Approximations.overapproximate(lm, Hyperrectangle) + oa_box = Approximations.overapproximate(lm, Approximations.BoxDirections) + d_oa_d_hp = Approximations.overapproximate(lm, CartesianProductArray) + d_oa_d_box = Approximations.overapproximate(lm, CartesianProductArray, Approximations.BoxDirections) + oa_d_hp = Approximations.overapproximate(d_oa_d_hp) + oa_d_box = Approximations.overapproximate(d_oa_d_box, Approximations.BoxDirections) + + @test oa == oa_d_hp + @test oa_box == oa_d_box + + for (oax, set_type) in [(d_oa_d_hp, Hyperrectangle), (d_oa_d_box, HPolytope)] + @test oax isa CartesianProductArray + arr = oax.array + @test length(arr) == 2 && dim(arr[1]) == 1 && dim(arr[2]) == 2 + @test all(X -> X isa set_type, arr) + end end # tests that do not work with Rational{Int} @@ -130,21 +154,4 @@ for N in [Float64, Float32] @test Y_polygon ⊆ Y_zonotope @test !(Y_zonotope ⊆ Y_polygon) - #Decomposed approximation of lazy linear map of CartesianProductArray - i1 = Interval([0, 1]) - h = Hyperrectangle(low=[3., 3.], high=[5., 5.]) - M = [1. 2. 3.; 4. 5. 6.; 7. 8. 9.] - cpa = CartesianProductArray([i1, h]) - lm = M * cpa - - oa = overapproximate(lm, Hyperrectangle) - oa_box = overapproximate(lm, Approximations.BoxDirections) - d_oa_d_hp = overapproximate(lm, CartesianProductArray) - d_oa_d_box = overapproximate(lm, CartesianProductArray, Approximations.BoxDirections) - oa_d_hp = overapproximate(d_oa_d_hp) - oa_d_box = overapproximate(d_oa_d_box, Approximations.BoxDirections) - - @test oa == oa_d_hp - @test oa_box == oa_d_box - end From ece7d747f876ca55647289eaebd84215ff3f43fe Mon Sep 17 00:00:00 2001 From: Unknown Date: Wed, 5 Jun 2019 05:23:36 +1000 Subject: [PATCH 10/15] fixed tests and type instability --- src/Approximations/overapproximate.jl | 2 +- test/unit_overapproximate.jl | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 25a1012d89..a92bac1140 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -561,7 +561,7 @@ end function blocks_linear_map(cpa::CartesianProductArray{N,<:LazySet{N}}, M::AbstractMatrix{N}, row_start_ind::Int, row_end_ind::Int) where {N} col_start_ind, col_end_ind = 1, 0 - h_min_sum = MinkowskiSumArray() + h_min_sum = MinkowskiSumArray(length(cpa.array), N) for bi in cpa.array col_end_ind += dim(bi) push!(h_min_sum.array, LinearMap(M[row_start_ind : row_end_ind, col_start_ind : col_end_ind], bi)) diff --git a/test/unit_overapproximate.jl b/test/unit_overapproximate.jl index 4e335bbe46..2362250d31 100644 --- a/test/unit_overapproximate.jl +++ b/test/unit_overapproximate.jl @@ -61,12 +61,12 @@ for N in [Float64, Rational{Int}, Float32] cpa = CartesianProductArray([i1, h]) lm = M * cpa - oa = Approximations.overapproximate(lm, Hyperrectangle) - oa_box = Approximations.overapproximate(lm, Approximations.BoxDirections) - d_oa_d_hp = Approximations.overapproximate(lm, CartesianProductArray) - d_oa_d_box = Approximations.overapproximate(lm, CartesianProductArray, Approximations.BoxDirections) - oa_d_hp = Approximations.overapproximate(d_oa_d_hp) - oa_d_box = Approximations.overapproximate(d_oa_d_box, Approximations.BoxDirections) + oa = overapproximate(lm, Hyperrectangle) + oa_box = overapproximate(lm, Approximations.BoxDirections) + d_oa_d_hp = overapproximate(lm, CartesianProductArray) + d_oa_d_box = overapproximate(lm, CartesianProductArray, Approximations.BoxDirections) + oa_d_hp = overapproximate(d_oa_d_hp) + oa_d_box = overapproximate(d_oa_d_box, Approximations.BoxDirections) @test oa == oa_d_hp @test oa_box == oa_d_box From 80eed63c8094a9889e1cfc648c05c0ea9ab4b2ad Mon Sep 17 00:00:00 2001 From: Unknown Date: Wed, 5 Jun 2019 17:55:26 +1000 Subject: [PATCH 11/15] small refactoring according to Christian's recommendations --- src/Approximations/overapproximate.jl | 36 +++++++++++++-------------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index a92bac1140..e12f8c7376 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -517,8 +517,8 @@ function overapproximate(cap::Intersection{N, end """ - overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, - <:LazySet{N}}}, ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} + overapproximate(lm::LinearMap{N, <:CartesianProductArray}, + ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} Decompose a lazy linear map of a cartesian product array. @@ -532,25 +532,21 @@ Decompose a lazy linear map of a cartesian product array. A `CartesianProductArray` representing the decomposed linear map. """ -function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, - <:LazySet{N}}}, ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} - - @assert size(lm.M, 2) == dim(lm.X) "matrix needs to be commensurate with the cartesian product" +function overapproximate(lm::LinearMap{N, <:CartesianProductArray}, + ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} cpa = lm.X M = lm.M - col_end_ind = 0 - array = allocate_result(oa, N) - sizehint!(array,length(cpa.array)) + sizehint!(array, length(cpa.array)) - i, n = 1, 0 + row_start_ind, row_end_ind = 1, 0 for bi in cpa.array - n += dim(bi) - ms = blocks_linear_map(cpa, M, i, n) + row_end_ind += dim(bi) + ms = blocks_linear_map(cpa, M, row_start_ind : row_end_ind) push!(array, overapproximate(ms, oa)) - i += n + row_start_ind += row_end_ind end result = CartesianProductArray(array) @@ -558,17 +554,19 @@ function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N, end -function blocks_linear_map(cpa::CartesianProductArray{N,<:LazySet{N}}, - M::AbstractMatrix{N}, row_start_ind::Int, row_end_ind::Int) where {N} +function blocks_linear_map(cpa::CartesianProductArray, + M::AbstractMatrix{N}, + row_range::UnitRange{Int64}) where {N} col_start_ind, col_end_ind = 1, 0 - h_min_sum = MinkowskiSumArray(length(cpa.array), N) + array = Vector{LazySet{N}}() + sizehint!(array, length(cpa.array)) for bi in cpa.array col_end_ind += dim(bi) - push!(h_min_sum.array, LinearMap(M[row_start_ind : row_end_ind, col_start_ind : col_end_ind], bi)) - col_start_ind = col_end_ind + 1 + push!(array, LinearMap(M[row_range, col_start_ind : col_end_ind], bi)) + col_start_ind += col_end_ind end - return h_min_sum + return MinkowskiSumArray(array) end function allocate_result(oa, N) From 13e7b9c860cce06a9a4921f71298e99d0e2b7d68 Mon Sep 17 00:00:00 2001 From: kostakoida Date: Wed, 5 Jun 2019 18:00:12 +1000 Subject: [PATCH 12/15] Update src/Approximations/overapproximate.jl Co-Authored-By: Christian Schilling --- src/Approximations/overapproximate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index e12f8c7376..5e8b1e3d87 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -556,7 +556,7 @@ end function blocks_linear_map(cpa::CartesianProductArray, M::AbstractMatrix{N}, - row_range::UnitRange{Int64}) where {N} + row_range::UnitRange{Int}) where {N} col_start_ind, col_end_ind = 1, 0 array = Vector{LazySet{N}}() sizehint!(array, length(cpa.array)) From e20224a083c2e5e67c4b8c61edb44ffdf3eeff54 Mon Sep 17 00:00:00 2001 From: Unknown Date: Thu, 6 Jun 2019 15:50:32 +1000 Subject: [PATCH 13/15] added new tests and refactoring added functionality to keep the original structure and preallocate corresponding set --- src/Approximations/overapproximate.jl | 95 ++++++++++++++++----------- test/unit_overapproximate.jl | 22 ++++++- 2 files changed, 77 insertions(+), 40 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 5e8b1e3d87..fe0052f13d 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -334,7 +334,7 @@ function overapproximate(X::LazySet{N}, end """ - overapproximate(S::LazySet{N}, ::Type{Interval}) where {N<:Real} + overapproximate(S::LazySet{N}, ::Union{Type{Interval}, Type{Interval{N}}}) where {N<:Real} Return the overapproximation of a real unidimensional set with an interval. @@ -352,11 +352,10 @@ An interval. The method relies on the exact conversion to `Interval`. Two support function evaluations are needed in general. """ -function overapproximate(S::LazySet{N}, ::Type{Interval}) where {N<:Real} +function overapproximate(S::LazySet{N}, ::Union{Type{Interval}, Type{Interval{N}}}) where {N<:Real} @assert dim(S) == 1 "cannot overapproximate a $(dim(S))-dimensional set with an `Interval`" return convert(Interval, S) end - function overapproximate_cap_helper(X::LazySet{N}, # compact set P::AbstractPolyhedron{N}, # polyhedron dir::AbstractDirections{N}; @@ -517,66 +516,84 @@ function overapproximate(cap::Intersection{N, end """ - overapproximate(lm::LinearMap{N, <:CartesianProductArray}, - ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} + overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}}, + ::Type{CartesianProductArray{N, O}}) where {N, T<:LazySet{N}, O<:LazySet{N}} -Decompose a lazy linear map of a cartesian product array. +Decompose a lazy linear map of a cartesian product array keeping original block structure. ### Input - `lm` -- lazy linear map of cartesian product array - `CartesianProductArray` -- type for dispatch -- `oa` -- approximation option for decomposition ### Output A `CartesianProductArray` representing the decomposed linear map. """ -function overapproximate(lm::LinearMap{N, <:CartesianProductArray}, - ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N} +function overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}}, + ::Type{CartesianProductArray{N, O}}) where {N, T<:LazySet{N}, O<:LazySet{N}} + + M, cpa = lm.M, lm.X + array = Vector{O}(undef, length(cpa.array)) + return _overapproximate_lm_cpa!(M, cpa, array, O) +end + +""" + overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}}, + ::Type{CartesianProductArray{N, HPolytope{N}}}, + dir::AbstractDirections{N}) where {N, T<:LazySet{N}} + +Decompose a lazy linear map of a cartesian product array. + +### Input + +- `lm` -- lazy linear map of cartesian product array +- `CartesianProductArray` -- type for dispatch +- `dir` -- approximation option for decomposition + +### Output - cpa = lm.X - M = lm.M +A `CartesianProductArray` representing the decomposed linear map. +""" +function overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}}, + ::Type{<:CartesianProductArray}, + dir::Type{<:AbstractDirections}) where {N, T<:LazySet{N}} + M, cpa = lm.M, lm.X + array = Vector{HPolytope{N}}(undef, length(cpa.array)) + return _overapproximate_lm_cpa!(M, cpa, array, dir) +end - array = allocate_result(oa, N) - sizehint!(array, length(cpa.array)) +#same for <:AbstractHyperrectangle +function overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}}, + ::Type{<:CartesianProductArray}, + oa::Type{<:AbstractHyperrectangle}) where {N, T<:LazySet{N}} + M, cpa = lm.M, lm.X + array = Vector{oa{N}}(undef, length(cpa.array)) + return _overapproximate_lm_cpa!(M, cpa, array, oa) +end +function _overapproximate_lm_cpa!(M, cpa, array, O) row_start_ind, row_end_ind = 1, 0 - for bi in cpa.array + @inbounds for (i, bi) in enumerate(cpa.array) row_end_ind += dim(bi) ms = blocks_linear_map(cpa, M, row_start_ind : row_end_ind) - push!(array, overapproximate(ms, oa)) - row_start_ind += row_end_ind + array[i] = overapproximate(ms, O) + row_start_ind = row_end_ind + 1 end - - result = CartesianProductArray(array) - return result + return CartesianProductArray(array) end - -function blocks_linear_map(cpa::CartesianProductArray, +function blocks_linear_map(cpa::CartesianProductArray{N, T}, M::AbstractMatrix{N}, - row_range::UnitRange{Int}) where {N} + row_range::UnitRange{Int}) where {N, T<:LazySet{N}} + array = Vector{LinearMap{N, <:T}}(undef, length(cpa.array)) + col_start_ind, col_end_ind = 1, 0 - array = Vector{LazySet{N}}() - sizehint!(array, length(cpa.array)) - for bi in cpa.array + @inbounds for (i, bi) in enumerate(cpa.array) col_end_ind += dim(bi) - push!(array, LinearMap(M[row_range, col_start_ind : col_end_ind], bi)) - col_start_ind += col_end_ind + array[i] = LinearMap(M[row_range, col_start_ind : col_end_ind], bi) + col_start_ind = col_end_ind + 1 end return MinkowskiSumArray(array) end - -function allocate_result(oa, N) - return Vector{LazySet{N}}() -end - -function allocate_result(oa::Type{<:LazySet}, N) - return Vector{oa{N}}() -end - -function allocate_result(oa::Type{<:AbstractDirections}, N) - return Vector{HPolytope{N}}() -end diff --git a/test/unit_overapproximate.jl b/test/unit_overapproximate.jl index 2362250d31..d524332fdb 100644 --- a/test/unit_overapproximate.jl +++ b/test/unit_overapproximate.jl @@ -63,7 +63,7 @@ for N in [Float64, Rational{Int}, Float32] oa = overapproximate(lm, Hyperrectangle) oa_box = overapproximate(lm, Approximations.BoxDirections) - d_oa_d_hp = overapproximate(lm, CartesianProductArray) + d_oa_d_hp = overapproximate(lm, CartesianProductArray{N, Hyperrectangle{N}}) d_oa_d_box = overapproximate(lm, CartesianProductArray, Approximations.BoxDirections) oa_d_hp = overapproximate(d_oa_d_hp) oa_d_box = overapproximate(d_oa_d_box, Approximations.BoxDirections) @@ -77,6 +77,26 @@ for N in [Float64, Rational{Int}, Float32] @test length(arr) == 2 && dim(arr[1]) == 1 && dim(arr[2]) == 2 @test all(X -> X isa set_type, arr) end + + i1 = Interval(N[0, 1]) + i2 = Interval(N[2, 3]) + i3 = Interval(N[1, 4]) + M = N[1 2; 0 1] + cpa = CartesianProductArray([i1, i2]) + lm = M * cpa + d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}}) + oa = overapproximate(lm, 1e-8) + @test oa ⊆ d_oa + + cpa = CartesianProductArray([i1, i2, i3]) + M = N[1 2 0; 0 1 0; 0 1 1] + lm = M * cpa + d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}}) + oa = overapproximate(lm) + @test overapproximate(d_oa) == oa + @test typeof(d_oa) == CartesianProductArray{N, Interval{N}} + + end # tests that do not work with Rational{Int} From 16f4a714ef204cd2a4690ea94a0a209f30bb5748 Mon Sep 17 00:00:00 2001 From: Unknown Date: Thu, 6 Jun 2019 16:34:15 +1000 Subject: [PATCH 14/15] fixed tests and function signature --- src/Approximations/overapproximate.jl | 4 ++-- test/unit_overapproximate.jl | 20 +++++++++++--------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index fe0052f13d..0861e1715a 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -334,7 +334,7 @@ function overapproximate(X::LazySet{N}, end """ - overapproximate(S::LazySet{N}, ::Union{Type{Interval}, Type{Interval{N}}}) where {N<:Real} + overapproximate(S::LazySet{N}, ::Type{<:Interval}) where {N<:Real} Return the overapproximation of a real unidimensional set with an interval. @@ -352,7 +352,7 @@ An interval. The method relies on the exact conversion to `Interval`. Two support function evaluations are needed in general. """ -function overapproximate(S::LazySet{N}, ::Union{Type{Interval}, Type{Interval{N}}}) where {N<:Real} +function overapproximate(S::LazySet{N}, ::Type{<:Interval}) where {N<:Real} @assert dim(S) == 1 "cannot overapproximate a $(dim(S))-dimensional set with an `Interval`" return convert(Interval, S) end diff --git a/test/unit_overapproximate.jl b/test/unit_overapproximate.jl index d524332fdb..66751742b4 100644 --- a/test/unit_overapproximate.jl +++ b/test/unit_overapproximate.jl @@ -81,13 +81,6 @@ for N in [Float64, Rational{Int}, Float32] i1 = Interval(N[0, 1]) i2 = Interval(N[2, 3]) i3 = Interval(N[1, 4]) - M = N[1 2; 0 1] - cpa = CartesianProductArray([i1, i2]) - lm = M * cpa - d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}}) - oa = overapproximate(lm, 1e-8) - @test oa ⊆ d_oa - cpa = CartesianProductArray([i1, i2, i3]) M = N[1 2 0; 0 1 0; 0 1 1] lm = M * cpa @@ -95,8 +88,7 @@ for N in [Float64, Rational{Int}, Float32] oa = overapproximate(lm) @test overapproximate(d_oa) == oa @test typeof(d_oa) == CartesianProductArray{N, Interval{N}} - - + end # tests that do not work with Rational{Int} @@ -174,4 +166,14 @@ for N in [Float64, Float32] @test Y_polygon ⊆ Y_zonotope @test !(Y_zonotope ⊆ Y_polygon) + #decomposed linear map approximation + i1 = Interval(N[0, 1]) + i2 = Interval(N[2, 3]) + M = N[1 2; 0 1] + cpa = CartesianProductArray([i1, i2]) + lm = M * cpa + d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}}) + oa = overapproximate(lm, OctDirections) + @test oa ⊆ d_oa + end From 03823d1511117d81b1b2763c37e8d969dc4134c2 Mon Sep 17 00:00:00 2001 From: Unknown Date: Thu, 6 Jun 2019 18:15:10 +1000 Subject: [PATCH 15/15] updated tests due to sparse --- test/unit_overapproximate.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/unit_overapproximate.jl b/test/unit_overapproximate.jl index 66751742b4..c988e90de1 100644 --- a/test/unit_overapproximate.jl +++ b/test/unit_overapproximate.jl @@ -88,7 +88,7 @@ for N in [Float64, Rational{Int}, Float32] oa = overapproximate(lm) @test overapproximate(d_oa) == oa @test typeof(d_oa) == CartesianProductArray{N, Interval{N}} - + end # tests that do not work with Rational{Int} @@ -165,7 +165,9 @@ for N in [Float64, Float32] Y_zonotope = overapproximate(Y, Zonotope) # overapproximate with a zonotope @test Y_polygon ⊆ Y_zonotope @test !(Y_zonotope ⊆ Y_polygon) +end +for N in [Float64] # due to sparse vectors: a = sparse(Float32[1 -1; 1 1];); a \ Float32[4, 10] #decomposed linear map approximation i1 = Interval(N[0, 1]) i2 = Interval(N[2, 3]) @@ -175,5 +177,4 @@ for N in [Float64, Float32] d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}}) oa = overapproximate(lm, OctDirections) @test oa ⊆ d_oa - end