diff --git a/src/Properties/check_blocks.jl b/src/Properties/check_blocks.jl index fe71ba33..e30d115e 100644 --- a/src/Properties/check_blocks.jl +++ b/src/Properties/check_blocks.jl @@ -33,10 +33,10 @@ The first time index where the property is violated, and 0 if the property is sa ϕpowerk_πbi[:, bj] @inline block(ϕpowerk_πbi::SparseMatrixCSC, bj::Int) = ϕpowerk_πbi[:, [bj]] -# sparse, with input +# sparse function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, Xhat0::Vector{<:LazySet{NUM}}, - U::ConstantNonDeterministicInput, + U::Union{ConstantNonDeterministicInput, Void}, overapproximate::Function, n::Int, N::Int, @@ -52,71 +52,19 @@ function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, b = length(blocks) Xhatk = Vector{LazySet{NUM}}(b) - Whatk = Vector{LazySet{NUM}}(b) - - inputs = next_set(U) - @inbounds for i in 1:b - bi = partition[blocks[i]] - Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) - end ϕpowerk = copy(ϕ) - k = 2 - p = Progress(N-1, 1, "Computing successors ") - @inbounds while true - update!(p, k) - for i in 1:b - bi = partition[blocks[i]] - Xhatk_bi = ZeroSet(length(bi)) - for (j, bj) in enumerate(partition) - block = ϕpowerk[bi, bj] - if findfirst(block) != 0 - Xhatk_bi = Xhatk_bi + block * Xhat0[j] - end - end - Xhatk[i] = Xhatk_bi + Whatk[i] - end - if !check_property(CartesianProductArray(Xhatk), prop) - return k - elseif k == N - break - end - - for i in 1:b + if U != nothing + Whatk = Vector{LazySet{NUM}}(b) + inputs = next_set(U) + @inbounds for i in 1:b bi = partition[blocks[i]] - Whatk[i] = - overapproximate(blocks[i], Whatk[i] + row(ϕpowerk, bi) * inputs) + Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) end - ϕpowerk = ϕpowerk * ϕ - k += 1 end - return 0 -end - - -# sparse, no input -function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, - Xhat0::Vector{<:LazySet{NUM}}, - n::Int, - N::Int, - blocks::AbstractVector{Int}, - partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, - prop::Property - )::Int where {NUM} - if !check_property(CartesianProductArray(Xhat0[blocks]), prop) - return 1 - elseif N == 1 - return 0 - end - - b = length(blocks) - Xhatk = Vector{LazySet{NUM}}(b) - - ϕpowerk = copy(ϕ) - k = 2 - p = Progress(N-1, 1, "Computing successors ") + p = Progress(N, 1, "Computing successors ") @inbounds while true update!(p, k) for i in 1:b @@ -128,14 +76,23 @@ function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, Xhatk_bi = Xhatk_bi + block * Xhat0[j] end end - Xhatk[i] = Xhatk_bi + Xhatk[i] = (U == nothing ? Xhatk_bi : Xhatk_bi + Whatk[i]) end + if !check_property(CartesianProductArray(Xhatk), prop) return k elseif k == N break end + if U != nothing + for i in 1:b + bi = partition[blocks[i]] + Whatk[i] = overapproximate(blocks[i], + Whatk[i] + row(ϕpowerk, bi) * inputs) + end + end + ϕpowerk = ϕpowerk * ϕ k += 1 end @@ -143,11 +100,10 @@ function check_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, return 0 end - -# dense, with input +# dense function check_blocks!(ϕ::AbstractMatrix{NUM}, Xhat0::Vector{<:LazySet{NUM}}, - U::ConstantNonDeterministicInput, + U::Union{ConstantNonDeterministicInput, Void}, overapproximate::Function, n::Int, N::Int, @@ -163,74 +119,21 @@ function check_blocks!(ϕ::AbstractMatrix{NUM}, b = length(blocks) Xhatk = Vector{LazySet{NUM}}(b) - Whatk = Vector{LazySet{NUM}}(b) - - inputs = next_set(U) - @inbounds for i in 1:b - bi = partition[blocks[i]] - Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) - end ϕpowerk = copy(ϕ) ϕpowerk_cache = similar(ϕ) - arr_length = length(partition) + 1 - k = 2 - p = Progress(N-1, 1, "Computing successors ") - @inbounds while true - update!(p, k) - for i in 1:b + if U != nothing + Whatk = Vector{LazySet{NUM}}(b) + inputs = next_set(U) + @inbounds for i in 1:b bi = partition[blocks[i]] - arr = Vector{LazySet{NUM}}(arr_length) - for (j, bj) in enumerate(partition) - arr[j] = ϕpowerk[bi, bj] * Xhat0[j] - end - arr[arr_length] = Whatk[i] - Xhatk[i] = MinkowskiSumArray(arr) - end - if !check_property(CartesianProductArray(Xhatk), prop) - return k - elseif k == N - break + Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) end - - for i in 1:b - bi = partition[blocks[i]] - Whatk[i] = - overapproximate(blocks[i], Whatk[i] + row(ϕpowerk, bi) * inputs) - end - A_mul_B!(ϕpowerk_cache, ϕpowerk, ϕ) - copy!(ϕpowerk, ϕpowerk_cache) - k += 1 end - return 0 -end - - -# dense, no input -function check_blocks!(ϕ::AbstractMatrix{NUM}, - Xhat0::Vector{<:LazySet{NUM}}, - n::Int, - N::Int, - blocks::AbstractVector{Int}, - partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, - prop::Property - )::Int where {NUM} - if !check_property(CartesianProductArray(Xhat0[blocks]), prop) - return 1 - elseif N == 1 - return 0 - end - - b = length(blocks) - Xhatk = Vector{LazySet{NUM}}(b) - - ϕpowerk = copy(ϕ) - ϕpowerk_cache = similar(ϕ) - - arr_length = length(partition) + arr_length = (U == nothing) ? length(partition) : length(partition) + 1 k = 2 - p = Progress(N-1, 1, "Computing successors ") + p = Progress(N, 1, "Computing successors ") @inbounds while true update!(p, k) for i in 1:b @@ -239,14 +142,26 @@ function check_blocks!(ϕ::AbstractMatrix{NUM}, for (j, bj) in enumerate(partition) arr[j] = ϕpowerk[bi, bj] * Xhat0[j] end + if U != nothing + arr[arr_length] = Whatk[i] + end Xhatk[i] = MinkowskiSumArray(arr) end + if !check_property(CartesianProductArray(Xhatk), prop) return k elseif k == N break end + if U != nothing + for i in 1:b + bi = partition[blocks[i]] + Whatk[i] = overapproximate(blocks[i], + Whatk[i] + row(ϕpowerk, bi) * inputs) + end + end + A_mul_B!(ϕpowerk_cache, ϕpowerk, ϕ) copy!(ϕpowerk, ϕpowerk_cache) k += 1 @@ -255,10 +170,11 @@ function check_blocks!(ϕ::AbstractMatrix{NUM}, return 0 end - -# lazymexp, no input +# lazy_expm function check_blocks!(ϕ::SparseMatrixExp{NUM}, Xhat0::Vector{<:LazySet{NUM}}, + U::Union{ConstantNonDeterministicInput, Void}, + overapproximate::Function, n::Int, N::Int, blocks::AbstractVector{Int}, @@ -273,68 +189,20 @@ function check_blocks!(ϕ::SparseMatrixExp{NUM}, b = length(blocks) Xhatk = Vector{LazySet{NUM}}(b) - ϕpowerk = SparseMatrixExp(copy(ϕ.M)) - arr_length = length(partition) - k = 2 - p = Progress(N-1, 1, "Computing successors ") - @inbounds while true - update!(p, k) - for i in 1:b + if U != nothing + Whatk = Vector{LazySet{NUM}}(b) + inputs = next_set(U) + @inbounds for i in 1:b bi = partition[blocks[i]] - arr = Vector{LazySet{NUM}}(arr_length) - ϕpowerk_πbi = get_rows(ϕpowerk, bi) - for (j, bj) in enumerate(partition) - arr[j] = block(ϕpowerk_πbi, bj) * Xhat0[j] - end - Xhatk[i] = MinkowskiSumArray(arr) - end - if !check_property(CartesianProductArray(Xhatk), prop) - return k - elseif k == N - break + Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) end - - ϕpowerk.M .= ϕpowerk.M + ϕ.M - k += 1 end - return 0 -end - - -# lazymexp, with input -function check_blocks!(ϕ::SparseMatrixExp{NUM}, - Xhat0::Vector{<:LazySet{NUM}}, - U::ConstantNonDeterministicInput, - overapproximate::Function, - n::Int, - N::Int, - blocks::AbstractVector{Int}, - partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, - prop::Property - )::Int where {NUM} - if !check_property(CartesianProductArray(Xhat0[blocks]), prop) - return 1 - elseif N == 1 - return 0 - end - - b = length(blocks) - Xhatk = Vector{LazySet{NUM}}(b) - Whatk = Vector{LazySet{NUM}}(b) - - inputs = next_set(U) - @inbounds for i in 1:b - bi = partition[blocks[i]] - Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) - end - ϕpowerk = SparseMatrixExp(copy(ϕ.M)) - - arr_length = length(partition) + 1 + arr_length = (U == nothing) ? length(partition) : length(partition) + 1 k = 2 - p = Progress(N-1, 1, "Computing successors ") + p = Progress(N, 1, "Computing successors ") @inbounds while true update!(p, k) for i in 1:b @@ -344,21 +212,27 @@ function check_blocks!(ϕ::SparseMatrixExp{NUM}, for (j, bj) in enumerate(partition) arr[j] = block(ϕpowerk_πbi, bj) * Xhat0[j] end - arr[arr_length] = Whatk[i] + if U != nothing + arr[arr_length] = Whatk[i] + end Xhatk[i] = MinkowskiSumArray(arr) end + if !check_property(CartesianProductArray(Xhatk), prop) return k elseif k == N break end - for i in 1:b - bi = partition[blocks[i]] - ϕpowerk_πbi = get_rows(ϕpowerk, bi) - Whatk[i] = - overapproximate(blocks[i], Whatk[i] + ϕpowerk_πbi * inputs) + if U != nothing + for i in 1:b + bi = partition[blocks[i]] + ϕpowerk_πbi = get_rows(ϕpowerk, bi) + Whatk[i] = + overapproximate(blocks[i], Whatk[i] + ϕpowerk_πbi * inputs) + end end + ϕpowerk.M .= ϕpowerk.M + ϕ.M k += 1 end diff --git a/src/Properties/check_property.jl b/src/Properties/check_property.jl index c460e7f8..ecbd5b2a 100644 --- a/src/Properties/check_property.jl +++ b/src/Properties/check_property.jl @@ -92,20 +92,19 @@ function check_property(S::AbstractSystem, end push!(args, Xhat0) - if !assume_homogeneous - push!(args, S.U) - - # overapproximation function (with or without iterative refinement) - if haskey(kwargs_dict, :block_types_iter) - block_types_iter = block_to_set_map(kwargs_dict[:block_types_iter]) - push!(args, (i, x) -> block_types_iter[i] == HPolygon ? - overapproximate(x, HPolygon, ε_iter) : - overapproximate(x, block_types_iter[i])) - elseif ε_iter < Inf - push!(args, (i, x) -> overapproximate(x, set_type_iter, ε_iter)) - else - push!(args, (i, x) -> overapproximate(x, set_type_iter)) - end + # inputs + push!(args, assume_homogeneous ? nothing : S.U) + + # overapproximation function (with or without iterative refinement) + if haskey(kwargs_dict, :block_types_iter) + block_types_iter = block_to_set_map(kwargs_dict[:block_types_iter]) + push!(args, (i, x) -> block_types_iter[i] == HPolygon ? + overapproximate(x, HPolygon, ε_iter) : + overapproximate(x, block_types_iter[i])) + elseif ε_iter < Inf + push!(args, (i, x) -> overapproximate(x, set_type_iter, ε_iter)) + else + push!(args, (i, x) -> overapproximate(x, set_type_iter)) end # ambient dimension diff --git a/src/ReachSets/reach.jl b/src/ReachSets/reach.jl index b0a8ccbd..f708c2cb 100644 --- a/src/ReachSets/reach.jl +++ b/src/ReachSets/reach.jl @@ -95,9 +95,8 @@ function reach(S::AbstractSystem, end push!(args, Xhat0) - if !assume_homogeneous - push!(args, S.U) - end + # inputs + push!(args, assume_homogeneous ? nothing : S.U) # overapproximation function (with or without iterative refinement) if haskey(kwargs_dict, :block_types_iter) diff --git a/src/ReachSets/reach_blocks.jl b/src/ReachSets/reach_blocks.jl index 4f740b45..f9ca242b 100644 --- a/src/ReachSets/reach_blocks.jl +++ b/src/ReachSets/reach_blocks.jl @@ -36,10 +36,10 @@ nondeterministic inputs. ϕpowerk_πbi[:, bj] @inline block(ϕpowerk_πbi::SparseMatrixCSC, bj::Int) = ϕpowerk_πbi[:, [bj]] -# sparse, with input +# sparse function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, Xhat0::Vector{<:LazySet{NUM}}, - U::ConstantNonDeterministicInput, + U::Union{ConstantNonDeterministicInput, Void}, overapproximate::Function, n::Int, N::Int, @@ -54,69 +54,17 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, b = length(blocks) Xhatk = Vector{LazySet{NUM}}(b) - Whatk = Vector{LazySet{NUM}}(b) - - inputs = next_set(U) - @inbounds for i in 1:b - bi = partition[blocks[i]] - Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) - end ϕpowerk = copy(ϕ) - k = 2 - p = Progress(N, 1, "Computing successors ") - @inbounds while true - update!(p, k) - for i in 1:b + if U != nothing + Whatk = Vector{LazySet{NUM}}(b) + inputs = next_set(U) + @inbounds for i in 1:b bi = partition[blocks[i]] - Xhatk_bi = ZeroSet(length(bi)) - for (j, bj) in enumerate(partition) - 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]) - end - res[k] = CartesianProductArray(copy(Xhatk)) - - if k == N - break + Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) end - - for i in 1:b - bi = partition[blocks[i]] - Whatk[i] = - overapproximate(blocks[i], Whatk[i] + row(ϕpowerk, bi) * inputs) - end - ϕpowerk = ϕpowerk * ϕ - k += 1 end - return nothing -end - - -# sparse, no input -function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, - Xhat0::Vector{<:LazySet{NUM}}, - overapproximate::Function, - n::Int, - N::Int, - blocks::AbstractVector{Int}, - partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, - res::Vector{CartesianProductArray{NUM}} - )::Void where {NUM} - res[1] = CartesianProductArray(Xhat0[blocks]) - if N == 1 - return nothing - end - - b = length(blocks) - Xhatk = Vector{LazySet{NUM}}(b) - - ϕpowerk = copy(ϕ) - k = 2 p = Progress(N, 1, "Computing successors ") @inbounds while true @@ -130,7 +78,8 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, Xhatk_bi = Xhatk_bi + block * Xhat0[j] end end - Xhatk[i] = overapproximate(blocks[i], Xhatk_bi) + Xhatk[i] = overapproximate(blocks[i], + U == nothing ? Xhatk_bi : Xhatk_bi + Whatk[i]) end res[k] = CartesianProductArray(copy(Xhatk)) @@ -138,6 +87,14 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, break end + if U != nothing + for i in 1:b + bi = partition[blocks[i]] + Whatk[i] = overapproximate(blocks[i], + Whatk[i] + row(ϕpowerk, bi) * inputs) + end + end + ϕpowerk = ϕpowerk * ϕ k += 1 end @@ -145,11 +102,10 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, return nothing end - -# dense, with input +# dense function reach_blocks!(ϕ::AbstractMatrix{NUM}, Xhat0::Vector{<:LazySet{NUM}}, - U::ConstantNonDeterministicInput, + U::Union{ConstantNonDeterministicInput, Void}, overapproximate::Function, n::Int, N::Int, @@ -164,72 +120,20 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM}, b = length(blocks) Xhatk = Vector{LazySet{NUM}}(b) - Whatk = Vector{LazySet{NUM}}(b) - - inputs = next_set(U) - @inbounds for i in 1:b - bi = partition[blocks[i]] - Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) - end ϕpowerk = copy(ϕ) ϕpowerk_cache = similar(ϕ) - arr_length = length(partition) + 1 - arr = Vector{LazySet{NUM}}(arr_length) - k = 2 - p = Progress(N, 1, "Computing successors ") - @inbounds while true - update!(p, k) - for i in 1:b - bi = partition[blocks[i]] - for (j, bj) in enumerate(partition) - arr[j] = ϕpowerk[bi, bj] * Xhat0[j] - end - arr[arr_length] = Whatk[i] - Xhatk[i] = overapproximate(blocks[i], MinkowskiSumArray(arr)) - end - res[k] = CartesianProductArray(copy(Xhatk)) - - if k == N - break - end - - for i in 1:b + if U != nothing + Whatk = Vector{LazySet{NUM}}(b) + inputs = next_set(U) + @inbounds for i in 1:b bi = partition[blocks[i]] - Whatk[i] = - overapproximate(blocks[i], Whatk[i] + row(ϕpowerk, bi) * inputs) + Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) end - A_mul_B!(ϕpowerk_cache, ϕpowerk, ϕ) - copy!(ϕpowerk, ϕpowerk_cache) - k += 1 - end - - return nothing -end - - -# dense, no input -function reach_blocks!(ϕ::AbstractMatrix{NUM}, - Xhat0::Vector{<:LazySet{NUM}}, - overapproximate::Function, - n::Int, - N::Int, - blocks::AbstractVector{Int}, - partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, - res::Vector{CartesianProductArray{NUM}} - )::Void where {NUM} - res[1] = CartesianProductArray(Xhat0[blocks]) - if N == 1 - return nothing end - b = length(blocks) - Xhatk = Vector{LazySet{NUM}}(b) - - ϕpowerk = copy(ϕ) - ϕpowerk_cache = similar(ϕ) - - arr = Vector{LazySet{NUM}}(length(partition)) + arr_length = (U == nothing) ? length(partition) : length(partition) + 1 + arr = Vector{LazySet{NUM}}(arr_length) k = 2 p = Progress(N, 1, "Computing successors ") @inbounds while true @@ -239,6 +143,9 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM}, for (j, bj) in enumerate(partition) arr[j] = ϕpowerk[bi, bj] * Xhat0[j] end + if U != nothing + arr[arr_length] = Whatk[i] + end Xhatk[i] = overapproximate(blocks[i], MinkowskiSumArray(arr)) end res[k] = CartesianProductArray(copy(Xhatk)) @@ -247,18 +154,27 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM}, break end + if U != nothing + for i in 1:b + bi = partition[blocks[i]] + Whatk[i] = overapproximate(blocks[i], + Whatk[i] + row(ϕpowerk, bi) * inputs) + end + end + A_mul_B!(ϕpowerk_cache, ϕpowerk, ϕ) copy!(ϕpowerk, ϕpowerk_cache) + k += 1 end return nothing end - -# lazymexp, no input +# lazy_expm function reach_blocks!(ϕ::SparseMatrixExp{NUM}, Xhat0::Vector{<:LazySet{NUM}}, + U::Union{ConstantNonDeterministicInput, Void}, overapproximate::Function, n::Int, N::Int, @@ -273,66 +189,17 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM}, b = length(blocks) Xhatk = Vector{LazySet{NUM}}(b) - ϕpowerk = SparseMatrixExp(copy(ϕ.M)) - k = 2 - p = Progress(N, 1, "Computing successors ") - @inbounds while true - update!(p, k) - for i in 1:b + if U != nothing + Whatk = Vector{LazySet{NUM}}(b) + inputs = next_set(U) + @inbounds for i in 1:b bi = partition[blocks[i]] - ϕpowerk_πbi = get_rows(ϕpowerk, bi) - Xhatk_bi = ZeroSet(length(bi)) - for (j, bj) in enumerate(partition) - πbi = block(ϕpowerk_πbi, bj) - if findfirst(πbi) != 0 - Xhatk_bi = Xhatk_bi + πbi * Xhat0[j] - end - end - Xhatk[i] = overapproximate(blocks[i], Xhatk_bi) - end - res[k] = CartesianProductArray(copy(Xhatk)) - - if k == N - break + Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) end - - ϕpowerk.M .= ϕpowerk.M + ϕ.M - k += 1 - end - - return nothing -end - - -# lazymexp, with input -function reach_blocks!(ϕ::SparseMatrixExp{NUM}, - Xhat0::Vector{<:LazySet{NUM}}, - U::ConstantNonDeterministicInput, - overapproximate::Function, - n::Int, - N::Int, - blocks::AbstractVector{Int}, - partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, - res::Vector{CartesianProductArray{NUM}} - )::Void where {NUM} - res[1] = CartesianProductArray(Xhat0[blocks]) - if N == 1 - return nothing end - b = length(blocks) - Xhatk = Vector{LazySet{NUM}}(b) - Whatk = Vector{LazySet{NUM}}(b) - - inputs = next_set(U) - @inbounds for i in 1:b - bi = partition[blocks[i]] - Whatk[i] = overapproximate(blocks[i], proj(bi, n) * inputs) - end - ϕpowerk = SparseMatrixExp(copy(ϕ.M)) - k = 2 p = Progress(N, 1, "Computing successors ") @inbounds while true @@ -347,9 +214,12 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM}, Xhatk_bi = Xhatk_bi + πbi * Xhat0[j] end end - Xhatk[i] = overapproximate(blocks[i], Xhatk_bi + Whatk[i]) - Whatk[i] = - overapproximate(blocks[i], Whatk[i] + ϕpowerk_πbi * inputs) + Xhatk[i] = overapproximate(blocks[i], + U == nothing ? Xhatk_bi : Xhatk_bi + Whatk[i]) + if U != nothing + Whatk[i] = + overapproximate(blocks[i], Whatk[i] + ϕpowerk_πbi * inputs) + end end res[k] = CartesianProductArray(copy(Xhatk))