Skip to content

Commit

Permalink
Merge pull request #119 from JuliaReach/schillic/111
Browse files Browse the repository at this point in the history
#111 - Fix slowdowns with sparse models and 1D blocks
  • Loading branch information
schillic authored Mar 26, 2018
2 parents ee68c59 + a9c2510 commit 49d4d4d
Show file tree
Hide file tree
Showing 5 changed files with 116 additions and 42 deletions.
48 changes: 28 additions & 20 deletions src/Properties/check_blocks.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,15 @@ OUTPUT:
The first time index where the property is violated, and 0 if the property is satisfied.
=#

# helper function
@inline G0(bi::AbstractVector{Int}, n::Int) =
sparse(1:length(bi), bi, ones(length(bi)), length(bi), n)
# helper functions
@inline proj(bi::UnitRange{Int}, n::Int) =
sparse(1:length(bi), bi, ones(length(bi)), length(bi), n)
@inline proj(bi::Int, n::Int) = sparse([1], [bi], ones(1), 1, n)
@inline row(ϕpowerk::AbstractMatrix, bi::UnitRange{Int}) = ϕpowerk[bi, :]
@inline row(ϕpowerk::AbstractMatrix, bi::Int) = ϕpowerk[[bi], :]
@inline block(ϕpowerk_πbi::SparseMatrixCSC, bi::UnitRange{Int}) =
ϕpowerk_πbi[:, bj]
@inline block(ϕpowerk_πbi::SparseMatrixCSC, bi::Int) = ϕpowerk_πbi[:, [bj]]

# sparse, with input
function check_blocks!::SparseMatrixCSC{NUM, Int},
Expand All @@ -35,7 +41,7 @@ function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
prop::Property
)::Int where {NUM}
if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
Expand All @@ -51,7 +57,7 @@ function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
inputs = next_set(U)
@inbounds for i in 1:b
bi = partition[blocks[i]]
Whatk[i] = overapproximate(blocks[i], G0(bi, n) * inputs)
Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs)
end
ϕpowerk = copy(ϕ)

Expand All @@ -63,8 +69,9 @@ function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
bi = partition[blocks[i]]
Xhatk_bi = ZeroSet(length(bi))
for (j, bj) in enumerate(partition)
if findfirst(ϕpowerk[bi, bj]) != 0
Xhatk_bi = Xhatk_bi + ϕpowerk[bi, bj] * Xhat0[j]
block = ϕpowerk[bi, bj]
if findfirst(block) != 0
Xhatk_bi = Xhatk_bi + block * Xhat0[j]
end
end
Xhatk[i] = Xhatk_bi + Whatk[i]
Expand All @@ -78,7 +85,7 @@ function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
for i in 1:b
bi = partition[blocks[i]]
Whatk[i] =
overapproximate(blocks[i], Whatk[i] + ϕpowerk[bi, :] * inputs)
overapproximate(blocks[i], Whatk[i] + row(ϕpowerk, bi) * inputs)
end
ϕpowerk = ϕpowerk * ϕ
k += 1
Expand All @@ -94,7 +101,7 @@ function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
prop::Property
)::Int where {NUM}
if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
Expand All @@ -116,8 +123,9 @@ function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
bi = partition[blocks[i]]
Xhatk_bi = ZeroSet(length(bi))
for (j, bj) in enumerate(partition)
if findfirst(ϕpowerk[bi, bj]) != 0
Xhatk_bi = Xhatk_bi + ϕpowerk[bi, bj] * Xhat0[j]
block = ϕpowerk[bi, bj]
if findfirst(block) != 0
Xhatk_bi = Xhatk_bi + block * Xhat0[j]
end
end
Xhatk[i] = Xhatk_bi
Expand All @@ -144,7 +152,7 @@ function check_blocks!(ϕ::AbstractMatrix{NUM},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
prop::Property
)::Int where {NUM}
if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
Expand All @@ -160,7 +168,7 @@ function check_blocks!(ϕ::AbstractMatrix{NUM},
inputs = next_set(U)
@inbounds for i in 1:b
bi = partition[blocks[i]]
Whatk[i] = overapproximate(blocks[i], G0(bi, n) * inputs)
Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs)
end
ϕpowerk = copy(ϕ)
ϕpowerk_cache = similar(ϕ)
Expand Down Expand Up @@ -188,7 +196,7 @@ function check_blocks!(ϕ::AbstractMatrix{NUM},
for i in 1:b
bi = partition[blocks[i]]
Whatk[i] =
overapproximate(blocks[i], Whatk[i] + ϕpowerk[bi, :] * inputs)
overapproximate(blocks[i], Whatk[i] + row(ϕpowerk, bi) * inputs)
end
A_mul_B!(ϕpowerk_cache, ϕpowerk, ϕ)
copy!(ϕpowerk, ϕpowerk_cache)
Expand All @@ -205,7 +213,7 @@ function check_blocks!(ϕ::AbstractMatrix{NUM},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
prop::Property
)::Int where {NUM}
if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
Expand Down Expand Up @@ -254,7 +262,7 @@ function check_blocks!(ϕ::SparseMatrixExp{NUM},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
prop::Property
)::Int where {NUM}
if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
Expand All @@ -278,7 +286,7 @@ function check_blocks!(ϕ::SparseMatrixExp{NUM},
arr = Vector{LazySet{NUM}}(arr_length)
ϕpowerk_πbi = get_rows(ϕpowerk, bi)
for (j, bj) in enumerate(partition)
arr[j] = ϕpowerk_πbi[:, bj] * Xhat0[j]
arr[j] = block(ϕpowerk_πbi, bj) * Xhat0[j]
end
Xhatk[i] = MinkowskiSumArray(arr)
end
Expand All @@ -304,7 +312,7 @@ function check_blocks!(ϕ::SparseMatrixExp{NUM},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
prop::Property
)::Int where {NUM}
if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
Expand All @@ -320,7 +328,7 @@ function check_blocks!(ϕ::SparseMatrixExp{NUM},
inputs = next_set(U)
@inbounds for i in 1:b
bi = partition[blocks[i]]
Whatk[i] = overapproximate(blocks[i], G0(bi, n) * inputs)
Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs)
end
ϕpowerk = SparseMatrixExp(copy.M))

Expand All @@ -334,7 +342,7 @@ function check_blocks!(ϕ::SparseMatrixExp{NUM},
arr = Vector{LazySet{NUM}}(arr_length)
ϕpowerk_πbi = get_rows(ϕpowerk, bi)
for (j, bj) in enumerate(partition)
arr[j] = ϕpowerk_πbi[:, bj] * Xhat0[j]
arr[j] = block(ϕpowerk_πbi, bj) * Xhat0[j]
end
arr[arr_length] = Whatk[i]
Xhatk[i] = MinkowskiSumArray(arr)
Expand Down
2 changes: 1 addition & 1 deletion src/Properties/check_property.jl
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ function check_property(S::AbstractSystem,
# add mode-specific block(s) argument
if algorithm == "explicit"
push!(args, kwargs_dict[:blocks])
push!(args, kwargs_dict[:partition])
push!(args, convert_partition(kwargs_dict[:partition]))
algorithm_backend = "explicit_blocks"
else
error("Unsupported algorithm: ", algorithm)
Expand Down
2 changes: 1 addition & 1 deletion src/ReachSets/reach.jl
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ function reach(S::AbstractSystem,
# preallocate output vector and add mode-specific block(s) argument
if algorithm == "explicit"
push!(args, kwargs_dict[:blocks])
push!(args, kwargs_dict[:partition])
push!(args, convert_partition(kwargs_dict[:partition]))
res = Vector{CartesianProductArray{numeric_type}}(N)
algorithm_backend = "explicit_blocks"
else
Expand Down
46 changes: 27 additions & 19 deletions src/ReachSets/reach_blocks.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,15 @@ It is obtained by reachability computation of a discrete affine system with
nondeterministic inputs.
=#

# helper function
@inline G0(bi::AbstractVector{Int}, n::Int) =
# helper functions
@inline proj(bi::UnitRange{Int}, n::Int) =
sparse(1:length(bi), bi, ones(length(bi)), length(bi), n)
@inline proj(bi::Int, n::Int) = sparse([1], [bi], ones(1), 1, n)
@inline row(ϕpowerk::AbstractMatrix, bi::UnitRange{Int}) = ϕpowerk[bi, :]
@inline row(ϕpowerk::AbstractMatrix, bi::Int) = ϕpowerk[[bi], :]
@inline block(ϕpowerk_πbi::SparseMatrixCSC, bi::UnitRange{Int}) =
ϕpowerk_πbi[:, bj]
@inline block(ϕpowerk_πbi::SparseMatrixCSC, bi::Int) = ϕpowerk_πbi[:, [bj]]

# sparse, with input
function reach_blocks!::SparseMatrixCSC{NUM, Int},
Expand All @@ -38,7 +44,7 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
res::Vector{CartesianProductArray{NUM}}
)::Void where {NUM}
res[1] = CartesianProductArray(Xhat0[blocks])
Expand All @@ -53,7 +59,7 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
inputs = next_set(U)
@inbounds for i in 1:b
bi = partition[blocks[i]]
Whatk[i] = overapproximate(blocks[i], G0(bi, n) * inputs)
Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs)
end
ϕpowerk = copy(ϕ)

Expand All @@ -65,8 +71,9 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
bi = partition[blocks[i]]
Xhatk_bi = ZeroSet(length(bi))
for (j, bj) in enumerate(partition)
if findfirst(ϕpowerk[bi, bj]) != 0
Xhatk_bi = Xhatk_bi + ϕpowerk[bi, bj] * Xhat0[j]
block = ϕpowerk[bi, bj]
if findfirst(block) != 0
Xhatk_bi = Xhatk_bi + block * Xhat0[j]
end
end
Xhatk[i] = overapproximate(blocks[i], Xhatk_bi + Whatk[i])
Expand All @@ -80,7 +87,7 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
for i in 1:b
bi = partition[blocks[i]]
Whatk[i] =
overapproximate(blocks[i], Whatk[i] + ϕpowerk[bi, :] * inputs)
overapproximate(blocks[i], Whatk[i] + row(ϕpowerk, bi) * inputs)
end
ϕpowerk = ϕpowerk * ϕ
k += 1
Expand All @@ -97,7 +104,7 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
res::Vector{CartesianProductArray{NUM}}
)::Void where {NUM}
res[1] = CartesianProductArray(Xhat0[blocks])
Expand All @@ -118,8 +125,9 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
bi = partition[blocks[i]]
Xhatk_bi = ZeroSet(length(bi))
for (j, bj) in enumerate(partition)
if findfirst(ϕpowerk[bi, bj]) != 0
Xhatk_bi = Xhatk_bi + ϕpowerk[bi, bj] * Xhat0[j]
block = ϕpowerk[bi, bj]
if findfirst(block) != 0
Xhatk_bi = Xhatk_bi + block * Xhat0[j]
end
end
Xhatk[i] = overapproximate(blocks[i], Xhatk_bi)
Expand All @@ -146,7 +154,7 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
res::Vector{CartesianProductArray{NUM}}
)::Void where {NUM}
res[1] = CartesianProductArray(Xhat0[blocks])
Expand All @@ -161,7 +169,7 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM},
inputs = next_set(U)
@inbounds for i in 1:b
bi = partition[blocks[i]]
Whatk[i] = overapproximate(blocks[i], G0(bi, n) * inputs)
Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs)
end
ϕpowerk = copy(ϕ)
ϕpowerk_cache = similar(ϕ)
Expand Down Expand Up @@ -189,7 +197,7 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM},
for i in 1:b
bi = partition[blocks[i]]
Whatk[i] =
overapproximate(blocks[i], Whatk[i] + ϕpowerk[bi, :] * inputs)
overapproximate(blocks[i], Whatk[i] + row(ϕpowerk, bi) * inputs)
end
A_mul_B!(ϕpowerk_cache, ϕpowerk, ϕ)
copy!(ϕpowerk, ϕpowerk_cache)
Expand All @@ -207,7 +215,7 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
res::Vector{CartesianProductArray{NUM}}
)::Void where {NUM}
res[1] = CartesianProductArray(Xhat0[blocks])
Expand Down Expand Up @@ -255,7 +263,7 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
res::Vector{CartesianProductArray{NUM}}
)::Void where {NUM}
res[1] = CartesianProductArray(Xhat0[blocks])
Expand All @@ -277,7 +285,7 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM},
ϕpowerk_πbi = get_rows(ϕpowerk, bi)
Xhatk_bi = ZeroSet(length(bi))
for (j, bj) in enumerate(partition)
πbi = ϕpowerk_πbi[:, bj]
πbi = block(ϕpowerk_πbi, bj)
if findfirst(πbi) != 0
Xhatk_bi = Xhatk_bi + πbi * Xhat0[j]
end
Expand Down Expand Up @@ -306,7 +314,7 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM},
n::Int,
N::Int,
blocks::AbstractVector{Int},
partition::AbstractVector{<:AbstractVector{Int}},
partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
res::Vector{CartesianProductArray{NUM}}
)::Void where {NUM}
res[1] = CartesianProductArray(Xhat0[blocks])
Expand All @@ -321,7 +329,7 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM},
inputs = next_set(U)
@inbounds for i in 1:b
bi = partition[blocks[i]]
Whatk[i] = overapproximate(blocks[i], G0(bi, n) * inputs)
Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs)
end
ϕpowerk = SparseMatrixExp(copy.M))

Expand All @@ -334,7 +342,7 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM},
ϕpowerk_πbi = get_rows(ϕpowerk, bi)
Xhatk_bi = ZeroSet(length(bi))
for (j, bj) in enumerate(partition)
πbi = ϕpowerk_πbi[:, bj]
πbi = block(ϕpowerk_πbi, bj)
if findfirst(πbi) != 0
Xhatk_bi = Xhatk_bi + πbi * Xhat0[j]
end
Expand Down
Loading

0 comments on commit 49d4d4d

Please sign in to comment.