From da0316b0151d3d004119ff719e9214b9910916c1 Mon Sep 17 00:00:00 2001 From: schillic Date: Mon, 4 Mar 2019 22:29:26 +0100 Subject: [PATCH] unify block_options --- .../ContinuousPost/BFFPSV18/BFFPSV18.jl | 190 ++++++------------ .../ContinuousPost/BFFPSV18/check_property.jl | 48 +---- .../ContinuousPost/BFFPSV18/reach.jl | 78 +++---- src/ReachSets/ReachSets.jl | 3 +- src/Utils/Utils.jl | 76 +------ test/Reachability/solve_continuous.jl | 10 +- 6 files changed, 118 insertions(+), 287 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl index ac7a9560..4db550f0 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl @@ -55,9 +55,6 @@ function options_BFFPSV18() info="use an analysis for sparse discretized matrices?"), # reachability options - OptionSpec(:lazy_X0, false, domain=Bool, - info="keep the discretized and decomposed initial states a lazy " * - "set?"), OptionSpec(:lazy_inputs_interval, lazy_inputs_interval_always, domain=Union{Int, Function}, domain_check=(v -> !(v isa Int) || v >= -1), @@ -66,49 +63,14 @@ function options_BFFPSV18() "predicate over indices; the default corresponds to ``-1``"), # approximation options - OptionSpec(:block_types, nothing, domain=Union{Nothing, - Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}}, - info="short hand to set ':block_types_init' and " * - "':block_types_iter'"), - OptionSpec(:block_types_init, nothing, domain=Union{Nothing, - Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}}, - info="set type for the approximation of the initial states for " * - "each block"), - OptionSpec(:block_types_iter, nothing, domain=Union{Nothing, - Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}}, - info="set type for the approximation of the states ``X_k``, " * - "``k>0``, for each block"), - OptionSpec(:ε, Inf, domain=Float64, domain_check=(v -> v > 0.), - info="short hand to set `:ε_init` and `:ε_iter`"), - OptionSpec(:ε_init, Inf, domain=Float64, domain_check=(v -> v > 0.), - info="error bound for the approximation of the initial states" * + OptionSpec(:block_options, nothing, domain=Any, + info="short hand to set ':block_options_init' and " * + "':block_options_iter'"), + OptionSpec(:block_options_init, nothing, domain=Any, + info="option for the approximation of the initial states " * "(during decomposition)"), - OptionSpec(:ε_iter, Inf, domain=Float64, domain_check=(v -> v > 0.), - info="error bound for the approximation of the states ``X_k``, " * - "``k>0``"), - OptionSpec(:set_type, Hyperrectangle, domain=Union{Type{HPolygon}, - Type{Hyperrectangle}, Type{LazySets.Interval}}, - info="short hand to set `:set_type_init` and `:set_type_iter`"), - OptionSpec(:set_type_init, Hyperrectangle, domain=Union{Type{HPolygon}, - Type{Hyperrectangle}, Type{LazySets.Interval}}, - info="set type for the approximation of the initial states" * - "(during decomposition)"), - OptionSpec(:set_type_iter, Hyperrectangle, domain=Union{Type{HPolygon}, - Type{Hyperrectangle}, Type{LazySets.Interval}}, - info="set type for the approximation of the states ``X_k``, " * - "``k>0``"), - OptionSpec(:template_directions, :nothing, domain=Symbol, - domain_check=(v::Symbol -> v in [:box, :oct, :boxdiag, :nothing]), - info="short hand to set `template_directions_init` and " * - "`template_directions_iter`"), - OptionSpec(:template_directions_init, :nothing, domain=Symbol, - domain_check=(v::Symbol -> v in [:box, :oct, :boxdiag, :nothing]), - info="directions to use for the approximation of the initial " * - "states (during decomposition)"), - OptionSpec(:template_directions_iter, :nothing, domain=Symbol, - domain_check=(v::Symbol -> v in [:box, :oct, :boxdiag, :nothing]), - info="directions to use for the approximation of the states " * - "``X_k``, ``k>0``, for each block"), + OptionSpec(:block_options_iter, nothing, domain=Any, + info="option for the approximation of the states ``X_k``, ``k>0``"), # convenience options OptionSpec(:assume_homogeneous, false, domain=Bool, @@ -133,76 +95,15 @@ function normalization_BFFPSV18!(𝑂::TwoLayerOptions) end end - # :block_types options - block_types = nothing - dict_type = Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}} - if !haskey_specified(𝑂, :block_types) && haskey(𝑂, :set_type) && - haskey_specified(𝑂, :partition) - 𝑂.specified[:block_types] = dict_type(𝑂[:set_type] => copy(𝑂[:partition])) - end - if !haskey_specified(𝑂, :block_types_init) && block_types != nothing - 𝑂.specified[:block_types_init] = block_types - end - if !haskey_specified(𝑂, :block_types_iter) && block_types != nothing - 𝑂.specified[:block_types_iter] = block_types - end - - # :ε, :set_type, and :template_directions options - ε = 𝑂[:ε] - if haskey_specified(𝑂, :set_type) - # use the provided set type - set_type = 𝑂[:set_type] - elseif ε < Inf - # use polygons - set_type = HPolygon - 𝑂[:set_type] = HPolygon - else - # use hyperrectangles - set_type = 𝑂[:set_type] - end - # - if !haskey_specified(𝑂, :ε_init) - 𝑂.specified[:ε_init] = - (haskey_specified(𝑂, :set_type_init) && 𝑂[:set_type_init] == HPolygon) || - (!haskey_specified(𝑂, :set_type_init) && set_type == HPolygon) ? - ε : - Inf - end - # - if !haskey_specified(𝑂, :set_type_init) - 𝑂.specified[:set_type_init] = 𝑂[:ε_init] < Inf ? HPolygon : set_type - end - # - if !haskey_specified(𝑂, :template_directions_init) - 𝑂.specified[:template_directions_init] = - haskey_specified(𝑂, :template_directions_init) ? - 𝑂[:template_directions_init] : - haskey_specified(𝑂, :template_directions) ? - 𝑂[:template_directions] : - :nothing - end - # - if !haskey_specified(𝑂, :ε_iter) - 𝑂.specified[:ε_iter] = - (haskey_specified(𝑂, :set_type_iter) && 𝑂[:set_type_iter] == HPolygon) || - (!haskey_specified(𝑂, :set_type_iter) && set_type == HPolygon) ? - ε : - Inf - end - # - if !haskey_specified(𝑂, :set_type_iter) - 𝑂.specified[:set_type_iter] = 𝑂[:ε_iter] < Inf ? HPolygon : set_type + # :block_options options: use convenience option for '_init' and '_iter' + if !haskey_specified(𝑂, :block_options_init) && + haskey_specified(𝑂, :block_options) + 𝑂.specified[:block_options_init] = 𝑂[:block_options] end - # - if !haskey_specified(𝑂, :template_directions_iter) - 𝑂.specified[:template_directions_iter] = - haskey_specified(𝑂, :template_directions_iter) ? - 𝑂[:template_directions_iter] : - haskey_specified(𝑂, :template_directions) ? - 𝑂[:template_directions] : - :nothing + if !haskey_specified(𝑂, :block_options_iter) && + haskey_specified(𝑂, :block_options) + 𝑂.specified[:block_options_iter] = 𝑂[:block_options] end - # nothing end @@ -214,27 +115,36 @@ function validation_BFFPSV18(𝑂) "':lazy_expm' with deactivated option ':lazy_expm_discretize'")) end - # block_types - if haskey_specified(𝑂, :block_types) - for (key, value) in 𝑂[:block_types] - if !(key <: LazySet) - throw(DomainError(key, "the keys of the `:block_types` " * - "dictionary should be lazy sets")) - elseif !(typeof(value) <: AbstractVector{<:AbstractVector{Int}}) - throw(DomainError(value, "the values of the `:block_types` " * - "dictionary should be vectors of " * - "vectors")) + # block_options + for b_options in [:block_options, :block_options_init, :block_options_iter] + if haskey_specified(𝑂, b_options) + bo = 𝑂[b_options] + if bo isa Type{<:LazySet} || bo isa Type{<:AbstractDirections} + # uniform options + elseif bo isa Symbol + # template directions + option = get(Utils.template_direction_symbols, bo, nothing) + if option == nothing + throw(DomainError(key, "if the `$b_options` option is a " * + "Symbol, it must be one of " * + "$(keys(Utils.template_direction_symbols))")) + end + 𝑂.specified[b_options] = option + elseif bo isa AbstractVector || bo isa Dict{Int, Any} + # mapping + elseif bo isa Real || bo isa Pair{<:UnionAll, <:Real} + ε = bo isa Real ? bo : bo[2] + if ε <= 0 + throw(DomainError(key, "the `$b_options` option must be " * + "positive")) + end + else + throw(DomainError(key, "the `$b_options` option does not accept " * + "$bo")) end end end - # ε-close approximation - if (𝑂[:ε_init] < Inf && 𝑂[:set_type_init] != HPolygon) || - (𝑂[:ε_iter] < Inf && 𝑂[:set_type_iter] != HPolygon) - throw(DomainError("ε-close approximation is only supported with the " * - "set type 'HPolygon'")) - end - nothing end @@ -312,6 +222,17 @@ function init!(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂::Options) # list of all interesting block indices in the partition 𝑂validated[:blocks] = compute_blocks(𝑂validated[:vars], 𝑂validated[:partition]) + # :block_options_init & :block_options_iter options: + # set default according to :partition + if !haskey_specified(𝒫.options, :block_options_init) + 𝑂validated[:block_options_init] = + compute_default_block_options(𝑂validated[:partition]) + end + if !haskey_specified(𝒫.options, :block_options_iter) + 𝑂validated[:block_options_iter] = + compute_default_block_options(𝑂validated[:partition]) + end + # Input -> Output variable mapping 𝑂validated[:inout_map] = inout_map_reach(𝑂validated[:partition], 𝑂validated[:blocks], 𝑂validated[:n]) @@ -410,3 +331,12 @@ function compute_blocks(vars, partition) sizehint!(blocks, length(blocks)) return blocks end + +function compute_default_block_options(partition) + # use Interval for 1D blocks and Hyperrectangle otherwise + block_options = Vector{Type{<:LazySet}}(undef, length(partition)) + for (i, block) in enumerate(partition) + block_options[i] = length(block) == 1 ? Interval : Hyperrectangle + end + return block_options +end diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl index 4046b27d..d5ca9e33 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl @@ -35,13 +35,7 @@ function check_property(S::IVP{<:AbstractDiscreteSystem}, n = statedim(S) blocks = options[:blocks] partition = convert_partition(options[:partition]) - dir = template_direction_symbols[options[:template_directions_init]] - block_sizes = compute_block_sizes(partition) N = ceil(Int, options[:T] / options[:δ]) - ε_init = options[:ε_init] - set_type_init = options[:set_type_init] - ε_iter = options[:ε_iter] - set_type_iter = options[:set_type_iter] # Cartesian decomposition of the initial set if length(partition) == 1 && length(partition[1]) == n @@ -50,22 +44,8 @@ function check_property(S::IVP{<:AbstractDiscreteSystem}, else info("- Decomposing X0") @timing begin - if options[:lazy_X0] - Xhat0 = array(decompose_helper(S.x0, block_sizes, n)) - elseif dir != nothing - Xhat0 = array(decompose(S.x0, directions=dir, - blocks=block_sizes)) - elseif options[:block_types_init] != nothing && - !isempty(options[:block_types_init]) - Xhat0 = array(decompose(S.x0, ε=ε_init, - block_types=options[:block_types_init])) - elseif set_type_init == LazySets.Interval - Xhat0 = array(decompose(S.x0, set_type=set_type_init, ε=ε_init, - blocks=ones(Int, n))) - else - Xhat0 = array(decompose(S.x0, set_type=set_type_init, ε=ε_init, - blocks=block_sizes)) - end + Xhat0 = array(decompose(S.x0, options[:partition], + options[:block_options_init])) end end @@ -88,21 +68,15 @@ function check_property(S::IVP{<:AbstractDiscreteSystem}, end push!(args, U) - # raw overapproximation function - dir = template_direction_symbols[options[:template_directions_iter]] - if dir != nothing - overapproximate_fun = - (i, x) -> overapproximate(x, dir(length(partition[i]))) - elseif options[:block_types_iter] != nothing - block_types_iter = block_to_set_map(options[:block_types_iter]) - overapproximate_fun = (i, x) -> block_types_iter[i] == HPolygon ? - overapproximate(x, HPolygon, ε_iter) : - overapproximate(x, block_types_iter[i]) - elseif ε_iter < Inf - overapproximate_fun = - (i, x) -> overapproximate(x, set_type_iter, ε_iter) + # overapproximation function for states + block_options_iter = options[:block_options_iter] + if block_options_iter isa AbstractVector || + block_options_iter isa Dict{Int, Any} + # individual overapproximation options per block + overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter[i]) else - overapproximate_fun = (i, x) -> overapproximate(x, set_type_iter) + # uniform overapproximation options for each block + overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter) end # overapproximate function for inputs @@ -117,7 +91,7 @@ function check_property(S::IVP{<:AbstractDiscreteSystem}, end # further sets of the series function _f(k, i, x::MinkowskiSum{NUM, <:CacheMinkowskiSum}) where NUM - if ε_iter == Inf + if has_constant_directions(block_options_iter, i) # forget sets if we do not use epsilon-close approximation forget_sets!(x.X) end diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl index 974bbf2c..ec8abe47 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl @@ -1,7 +1,9 @@ -import LazySets: CacheMinkowskiSum, +using ..Utils: LDS, CLCDS + +using LazySets: CacheMinkowskiSum, isdisjoint -import ..Utils: LDS, CLCDS +import LazySets.Approximations: overapproximate """ reach(S, invariant, options) @@ -48,14 +50,7 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}}, n = statedim(S) blocks = options[:blocks] partition = convert_partition(options[:partition]) - dir = template_direction_symbols[options[:template_directions_init]] - block_sizes = compute_block_sizes(partition) N = ceil(Int, options[:T] / options[:δ]) - ε_init = options[:ε_init] - set_type_init = options[:set_type_init] - ε_iter = options[:ε_iter] - set_type_iter = options[:set_type_iter] - # Cartesian decomposition of the initial set if length(partition) == 1 && length(partition[1]) == n @@ -64,22 +59,8 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}}, else info("- Decomposing X0") @timing begin - if options[:lazy_X0] - Xhat0 = array(decompose_helper(S.x0, block_sizes, n)) - elseif dir != nothing - Xhat0 = array(decompose(S.x0, directions=dir, - blocks=block_sizes)) - elseif options[:block_types_init] != nothing && - !isempty(options[:block_types_init]) - Xhat0 = array(decompose(S.x0, ε=ε_init, - block_types=options[:block_types_init])) - elseif set_type_init == LazySets.Interval - Xhat0 = array(decompose(S.x0, set_type=set_type_init, ε=ε_init, - blocks=ones(Int, n))) - else - Xhat0 = array(decompose(S.x0, set_type=set_type_init, ε=ε_init, - blocks=block_sizes)) - end + Xhat0 = array(decompose(S.x0, options[:partition], + options[:block_options_init])) end end @@ -114,19 +95,14 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}}, push!(args, U) # overapproximation function for states - dir = template_direction_symbols[options[:template_directions_iter]] - if dir != nothing - overapproximate_fun = (i, x) -> overapproximate(x, dir(length(partition[i]))) - elseif options[:block_types_iter] != nothing - block_types_iter = block_to_set_map(options[:block_types_iter]) - overapproximate_fun = (i, x) -> (block_types_iter[i] == HPolygon) ? - overapproximate(x, HPolygon, ε_iter) : - overapproximate(x, block_types_iter[i]) - elseif ε_iter < Inf - overapproximate_fun = - (i, x) -> overapproximate(x, set_type_iter, ε_iter) + block_options_iter = options[:block_options_iter] + if block_options_iter isa AbstractVector || + block_options_iter isa Dict{Int, Any} + # individual overapproximation options per block + overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter[i]) else - overapproximate_fun = (i, x) -> overapproximate(x, set_type_iter) + # uniform overapproximation options for each block + overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter) end push!(args, overapproximate_fun) @@ -142,7 +118,7 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}}, end # further sets of the series function _f(k, i, x::MinkowskiSum{NUM, <:CacheMinkowskiSum}) where NUM - if ε_iter == Inf + if has_constant_directions(block_options_iter, i) # forget sets if we do not use epsilon-close approximation forget_sets!(x.X) end @@ -252,3 +228,29 @@ function termination_inv_N(N, inv, k, set, t0) return (false, false) end end + +function overapproximate(X::LazySet, pair::Pair) + return overapproximate(X, pair[1], pair[2]) +end + +function has_constant_directions(block_options::AbstractVector, i::Int) + return has_constant_directions(block_options[i]) +end + +function has_constant_directions(block_options::Dict{<:UnionAll, <:Real}, + i::Int) + return has_constant_directions(block_options[i]) +end + +function has_constant_directions(block_options::Pair{<:UnionAll, <:Real}, + i::Int) + return has_constant_directions(block_options[2]) +end + +function has_constant_directions(block_options::Real, i::Int) + return ε == Inf +end + +function has_constant_directions(block_options, i::Int) + return true +end diff --git a/src/ReachSets/ReachSets.jl b/src/ReachSets/ReachSets.jl index e1f4e404..3c49d29d 100644 --- a/src/ReachSets/ReachSets.jl +++ b/src/ReachSets/ReachSets.jl @@ -17,7 +17,8 @@ include("../compat.jl") using LazySets.Approximations: symmetric_interval_hull, decompose, overapproximate, - box_approximation + box_approximation, + AbstractDirections using Reachability: @timing, Options, OptionSpec, TwoLayerOptions, AbstractOptions, validate_and_wrap_options, print_option_spec, diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index 42ca507d..ae869896 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -20,9 +20,7 @@ export print_sparsity, # Block Structure export @block_id, add_dimension, - block_to_set_map, - convert_partition, - compute_block_sizes + convert_partition # Usability export @filename_to_png, @@ -33,9 +31,6 @@ export template_direction_symbols, matrix_conversion, matrix_conversion_lazy_explicit -# temporary helper function -export decompose_helper - # Extension of MathematicalSystems for use inside Reachability.jl include("systems.jl") @@ -327,35 +322,6 @@ macro relpath(name::String) end end -""" - block_to_set_map(dict::Dict{Type{<:LazySet}, - AbstractVector{<:AbstractVector{Int}}}) - -Invert a map (set type -> block structure) to a map (block index -> set type). - -### Input - -- `dict` -- map (set type -> block structure) - -### Ouput - -Vector mapping a block index to the respective set type. -""" -function block_to_set_map(dict::Dict{Type{<:LazySet}, - AbstractVector{<:AbstractVector{Int}}}) - # we assume that blocks do not overlap - set_type = Vector{Type{<:LazySet}}() - initial_block_indices = Vector{Int}() - @inbounds for (key, val) in dict - for bi in val - push!(set_type, key) - push!(initial_block_indices, bi[1]) - end - end - s = sortperm(initial_block_indices) - return set_type[s] -end - """ convert_partition(partition::AbstractVector{<:AbstractVector{Int}})::Union{ Vector{Int}, Vector{UnitRange{Int}}, Vector{Union{UnitRange{Int}, Int}}} @@ -413,46 +379,6 @@ template_direction_symbols = Dict( :boxdiag => Approximations.BoxDiagDirections ) -""" - compute_block_sizes(partition::Union{Vector{Int}, Vector{UnitRange{Int}}} - )::Vector{Int} - -Conversion of the blocks representation from partition to block sizes. - -### Input - -- `partition` -- partition, a vector of block index ranges - -### Output - -A vector where the ``i``th entry contains the size of the ``i``th block. -""" -function compute_block_sizes(partition::Union{Vector{Int}, Vector{UnitRange{Int}}} - )::Vector{Int} - res = Vector{Int}(undef, length(partition)) - for (i, block) in enumerate(partition) - res[i] = length(block) - end - return res -end - -""" -This is a temporary decomposition function that only applies a projection and -keeps the sets lazy. -Eventually this should be handled in LazySets. -""" -function decompose_helper(S::LazySet{N}, blocks::AbstractVector{Int}, - n::Int=dim(S)) where {N} - result = Vector{LazySet{N}}(undef, length(blocks)) - block_start = 1 - @inbounds for (i, bi) in enumerate(blocks) - M = sparse(1:bi, block_start:(block_start + bi - 1), ones(N, bi), bi, n) - result[i] = M * S - block_start += bi - end - return CartesianProductArray(result) -end - # sparse/dense matrix conversion function matrix_conversion(Δ, options; A_passed=nothing) if A_passed == nothing diff --git a/test/Reachability/solve_continuous.jl b/test/Reachability/solve_continuous.jl index a429ccae..0354274c 100644 --- a/test/Reachability/solve_continuous.jl +++ b/test/Reachability/solve_continuous.jl @@ -65,8 +65,7 @@ s = solve(ContinuousSystem(A, X0), # template directions (eg. :box, :oct, :octbox) s = solve(ContinuousSystem(A, X0), Options(:T=>0.1, :ε_proj=>1e-5, :set_type_proj=>HPolygon), - op=BFFPSV18(:vars=>[1,3], :partition=>[1:4], - :template_directions => :oct)) + op=BFFPSV18(:vars=>[1,3], :partition=>[1:4], :block_options => :oct)) s = solve(ContinuousSystem(A, X0), Options(:T=>0.1), @@ -92,16 +91,15 @@ s = solve(ContinuousSystem(sparse(A), X0), Options(:T=>0.1), op=BFFPSV18(:vars=>[1,3], :partition=>[[i] for i in 1:4], :lazy_expm=>true, :lazy_expm_discretize=>true, - :set_type=>Interval)) + :block_options=>Interval)) # =============================== # System with an odd dimension # =============================== A = randn(5, 5); X0 = BallInf(ones(5), 0.1) s = solve(ContinuousSystem(A, X0), Options(:T=>0.1), - op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4, [5]], :block_types=> - Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}( - HPolygon=>[1:2, 3:4], Interval=>[[5]]))) + op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4, [5]], :block_options=> + [HPolygon, HPolygon, Interval])) # =============================== # System with an odd dimension