From 5097fca62dc50e38bf86b50f941a0035072ab7a3 Mon Sep 17 00:00:00 2001 From: schillic Date: Wed, 11 Sep 2019 10:56:17 +0200 Subject: [PATCH] remove inout_map option --- .../ContinuousPost/BFFPS19/BFFPS19.jl | 3 - .../ContinuousPost/BFFPSV18/BFFPSV18.jl | 3 - .../BFFPSV18/inout_map_reach.jl | 44 ---- src/ReachSets/ReachSets.jl | 1 - src/ReachSets/project_reach.jl | 230 +++++------------- src/solve.jl | 11 - 6 files changed, 59 insertions(+), 233 deletions(-) delete mode 100644 src/ReachSets/ContinuousPost/BFFPSV18/inout_map_reach.jl diff --git a/src/ReachSets/ContinuousPost/BFFPS19/BFFPS19.jl b/src/ReachSets/ContinuousPost/BFFPS19/BFFPS19.jl index a0b9afdc..a8226667 100644 --- a/src/ReachSets/ContinuousPost/BFFPS19/BFFPS19.jl +++ b/src/ReachSets/ContinuousPost/BFFPS19/BFFPS19.jl @@ -215,9 +215,6 @@ function init!(š¯’«::BFFPS19, š¯‘†::AbstractSystem, š¯‘‚::Options) compute_default_block_options(š¯‘‚validated[:partition]) end - # Input -> Output variable mapping - š¯‘‚validated[:inout_map] = inout_map_reach(š¯‘‚validated[:partition], š¯‘‚validated[:blocks], š¯‘‚validated[:n]) - if š¯‘‚validated[:project_reachset] š¯‘‚validated[:output_function] = nothing else diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl index 2fb3667a..3a30e05d 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl @@ -222,9 +222,6 @@ function init!(š¯’«::BFFPSV18, š¯‘†::AbstractSystem, š¯‘‚::Options) compute_default_block_options(š¯‘‚validated[:partition]) end - # Input -> Output variable mapping - š¯‘‚validated[:inout_map] = inout_map_reach(š¯‘‚validated[:partition], š¯‘‚validated[:blocks], š¯‘‚validated[:n]) - if š¯‘‚validated[:project_reachset] š¯‘‚validated[:output_function] = nothing else diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/inout_map_reach.jl b/src/ReachSets/ContinuousPost/BFFPSV18/inout_map_reach.jl deleted file mode 100644 index b5454ad6..00000000 --- a/src/ReachSets/ContinuousPost/BFFPSV18/inout_map_reach.jl +++ /dev/null @@ -1,44 +0,0 @@ -""" - inout_map_reach(partition::AbstractVector{<:AbstractVector{Int}}, - blocks::AbstractVector{Int}, - n::Int) - -Compute a mapping from input dimension to output dimension. -The map returns `0` for all input dimensions not present in the output. - -### Input - -- `partition` -- block partition; elements are start and end indices of a block -- `blocks` -- list of all output block indices in the partition -- `n` -- total number of input dimensions - -### Output - -A vector mapping from input dimension (index) to output dimension (entry). -""" -function inout_map_reach(partition::AbstractVector{<:AbstractVector{Int}}, - blocks::AbstractVector{Int}, - n::Int) - inout = Vector{Int}(undef, n) - blocks_idx = 1 - next_block_idx = blocks[blocks_idx] - out_idx = 1 - for (i, bi) in enumerate(partition) - if next_block_idx == i - # block/variables present in output - inout[bi] = out_idx:(out_idx + length(bi) - 1) - out_idx += length(bi) - if blocks_idx < length(blocks) - blocks_idx += 1 - next_block_idx = blocks[blocks_idx] - else - next_block_idx = 0 - end - else - @assert next_block_idx > i || next_block_idx == 0 - # block/variables not present in output - inout[bi] .= 0 - end - end - return inout -end diff --git a/src/ReachSets/ReachSets.jl b/src/ReachSets/ReachSets.jl index a44fd4d5..79849b27 100644 --- a/src/ReachSets/ReachSets.jl +++ b/src/ReachSets/ReachSets.jl @@ -53,7 +53,6 @@ include("ContinuousPost/BFFPSV18/check_blocks.jl") include("ContinuousPost/BFFPSV18/check_property.jl") include("ContinuousPost/BFFPSV18/partitions.jl") include("ContinuousPost/BFFPSV18/compute_dimensions.jl") -include("ContinuousPost/BFFPSV18/inout_map_reach.jl") # "explicit" backends push!(available_algorithms_check, "explicit_blocks"=>Dict("func"=>check_blocks, diff --git a/src/ReachSets/project_reach.jl b/src/ReachSets/project_reach.jl index 934426ee..10463677 100644 --- a/src/ReachSets/project_reach.jl +++ b/src/ReachSets/project_reach.jl @@ -1,5 +1,5 @@ """ - project_reach(Rsets, vars, n, options) + project_reach(Rsets, vars, options) Projection of a reachability analysis result in 2D. @@ -7,7 +7,6 @@ Projection of a reachability analysis result in 2D. - `Rsets` -- reachable states representation - `vars` -- variables to plot; two-dimensional index vector -- `n` -- system dimension - `options` -- options ### Notes @@ -17,186 +16,83 @@ The `vars` argument is required even if the optional argument a dimension from this variable. """ function project_reach( - Rsets::Vector{<:AbstractReachSet{<:LazySets.LazySet{numeric_type}}}, + Rsets::Vector{<:AbstractReachSet{<:LazySets.LazySet{N}}}, vars::Vector{Int64}, - n::Int64, - options::AbstractOptions)::Vector{<:AbstractReachSet} where {numeric_type<:Real} - + options::AbstractOptions)::Vector{<:AbstractReachSet} where {N<:Real} # parse input - @assert(length(vars) == 2) - if n == 2 - return project_2d_reach(Rsets, vars, n, options) - end - - # first projection dimension + @assert length(vars) == 2 "we only support projection to two dimensions" xaxis = vars[1] - if xaxis == 0 - got_time = true - xaxis = n+1 # we add a new dimension for time - else - got_time = false - if (xaxis <= 0 || xaxis > n) - throw(DomainError("value $xaxis for X variable not allowed")) - end - end - - # build projection matrix - projection_matrix = options[:projection_matrix] - output_function = !options[:project_reachset] - m = got_time ? n+1 : n - if projection_matrix == nothing - # projection to a state variable - yaxis = vars[2] - if (yaxis <= 0 || yaxis > n) - throw(DomainError("value $yaxis for Y variable not allowed")) - end - projection_matrix = sparse([1, 2], [xaxis, yaxis], [1.0, 1.0], 2, m) - else - @assert(size(projection_matrix) == (1,n)) - # make vector a 2-row matrix - projection_matrix = sparse(fill(2, n), 1:n, projection_matrix[1, :], 2, m) - projection_matrix[1, xaxis] = 1.0 + yaxis = vars[2] + got_time = (xaxis == 0) + if xaxis < 0 + throw(DomainError("value $xaxis for x variable not allowed")) + elseif yaxis <= 0 + throw(DomainError("value $yaxis for y variable not allowed")) end - - # apply optional transformation to projection matrix - if options[:transformation_matrix] != nothing - transformation_matrix = options[:transformation_matrix] - if got_time - # add another dimension for time: block matrix [S 0; 0 1] - transformation_matrix = sparse(cat([1, 2], transformation_matrix, [1])) - end - projection_matrix = projection_matrix * transformation_matrix - end - - N = length(Rsets) + n = options[:n] + projection_matrix_high_dimensional = options[:projection_matrix] + transformation_matrix = options[:transformation_matrix] # allocate output and define overapproximation function Īµ = options[:Īµ_proj] if Īµ < Inf oa = x -> overapproximate(x, HPolygon, Īµ) - RsetsProj = Vector{ReachSet{HPolygon{numeric_type}}}(undef, N) + RsetsProj = Vector{ReachSet{HPolygon{N}}}(undef, length(Rsets)) else set_type = options[:set_type_proj] oa = x -> overapproximate(x, set_type) - RsetsProj = Vector{ReachSet{set_type{numeric_type}}}(undef, N) + RsetsProj = Vector{ReachSet{set_type{N}}}(undef, length(Rsets)) end - if got_time - @inbounds for i in 1:N - t0 = time_start(Rsets[i]) - t1 = time_end(Rsets[i]) - radius = (t1 - t0)/2.0 - RsetsProj[i] = ReachSet( - oa(projection_matrix * - CartesianProduct(set(Rsets[i]), - BallInf([t0 + radius], radius))), t0, t1) - end - else - @inbounds for i in 1:N - RsetsProj[i] = ReachSet( - oa(projection_matrix * set(Rsets[i])), - time_start(Rsets[i]), time_end(Rsets[i])) - end - end - - return RsetsProj -end - -""" - project_reach(Rsets, vars, n, options) - -This function projects a sequence of sets into the time variable, or can be -used to take a linear combination of the given variables. - -### Input - -- `Rsets` -- reachable states representation -- `vars` -- variables to plot; two-dimensional index vector -- `n` -- system dimension -- `options` -- options - -### Notes - -The input `Rsets` is an array of sets (instead of a `CartesianProductArray`). -This array contains the collection of reach sets in 2D. - -It is assumed that the variable given in vars belongs to the block computed -in the sequence of 2D sets `Rsets`. -""" -function project_2d_reach( - Rsets::Vector{<:AbstractReachSet{<:LazySets.LazySet{numeric_type}}}, - vars::Vector{Int64}, - n::Int64, - options::AbstractOptions)::Vector{<:AbstractReachSet} where {numeric_type<:Real} - - # first projection dimension - xaxis = vars[1] - if xaxis == 0 - got_time = true - xaxis = 3 # time is associated to dimension 3 - else - got_time = false - if (xaxis <= 0 || xaxis > n) - throw(DomainError()) + @inbounds for (i, rs) in enumerate(Rsets) + X = set(rs) + if projection_matrix_high_dimensional == nothing + # use a simple projection to state variables + if got_time + # projection to a single state variable + projection_matrix = sparse([1], [yaxis], [1.0], 1, n) + else + # projection to two state variables + projection_matrix = + sparse([1, 2], [xaxis, yaxis], [1.0, 1.0], 2, n) + end + else + @assert size(projection_matrix, 1) == 1 "currently we only " * + "support one-dimensional projection matrices" + if got_time + # create a 1-row matrix + projection_matrix = + sparse(fill(1, n), 1:n, projection_matrix[1, :], 1, n) + else + # create a 2-row matrix + projection_matrix = + sparse(fill(2, n), 1:n, projection_matrix[1, :], 2, n) + projection_matrix[1, xaxis] = 1.0 + end end - # map back to 2d block - xaxis = iseven(xaxis) ? 2 : 1 - end - # build projection matrix - projection_matrix = options[:projection_matrix] - output_function = !options[:project_reachset] - if projection_matrix == nothing - # projection to a state variable - yaxis = vars[2] - if (yaxis <= 0 || yaxis > n) - throw(DomainError()) + # apply optional transformation to projection matrix + if transformation_matrix != nothing + if got_time + # add another dimension for time: block matrix [S 0; 0 1] + transformation_matrix = + sparse(cat([1, 2], transformation_matrix, [1])) + end + projection_matrix = projection_matrix * transformation_matrix end - # map back to 2d block - yaxis = iseven(yaxis) ? 2 : 1 - m = got_time ? 3 : 2 - projection_matrix = sparse([1, 2], [xaxis, yaxis], [1.0, 1.0], 2, m) - elseif !output_function - error("projection matrix not allowed for this algorithm") - end - - N = length(Rsets) - # allocate output and define overapproximation function - Īµ = options[:Īµ_proj] - if Īµ < Inf - oa = x -> overapproximate(x, HPolygon, Īµ) - RsetsProj = Vector{ReachSet{HPolygon{numeric_type}}}(undef, N) - else - set_type = options[:set_type_proj] - oa = x -> overapproximate(x, set_type) - RsetsProj = Vector{ReachSet{set_type{numeric_type}}}(undef, N) - end + # project set + rs = project(rs, projection_matrix) + projected = set(rs) - if output_function - @inbounds for i in 1:N - t0 = time_start(Rsets[i]) - t1 = time_end(Rsets[i]) - radius = (t1 - t0)/2.0 - RsetsProj[i] = ReachSet( - oa(CartesianProduct(BallInf([t0 + radius], radius), - set(Rsets[i]))), t0, t1) - end - elseif got_time # x variable is 'time' - @inbounds for i in 1:N - t0 = time_start(Rsets[i]) - t1 = time_end(Rsets[i]) - radius = (t1 - t0)/2.0 - RsetsProj[i] = ReachSet( - oa(projection_matrix * - CartesianProduct(set(Rsets[i]), - BallInf([t0 + radius], radius))), t0, t1) - end - else - @inbounds for i in 1:N - RsetsProj[i] = ReachSet(oa(projection_matrix * set(Rsets[i])), - time_start(Rsets[i]), time_end(Rsets[i])) + # add time dimension + t0 = time_start(rs) + t1 = time_end(rs) + if got_time + time_interval = Interval(t0, t1) + projected = CartesianProduct(time_interval, projected) end + RsetsProj[i] = ReachSet(oa(projected), t0, t1) end return RsetsProj @@ -218,15 +114,7 @@ A projection matrix can be given in the options structure, or passed as a dictionary entry. """ function project(Rsets::Vector{<:AbstractReachSet}, options::AbstractOptions) - plot_vars = copy(options[:plot_vars]) - for i in 1:length(plot_vars) - if plot_vars[i] != 0 - plot_vars[i] = options[:inout_map][plot_vars[i]] - end - end - reduced_n = sum(x -> x != 0, options[:inout_map]) - output_function = !options[:project_reachset] - RsetsProj = project_reach(Rsets, plot_vars, reduced_n, options) + return project_reach(Rsets, options[:plot_vars], options) end project(reach_sol::ReachSolution) = project(reach_sol.Xk, reach_sol.options) diff --git a/src/solve.jl b/src/solve.jl index 0239ea56..d4be5be3 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -115,9 +115,6 @@ function solve!(system::InitialValueProblem{<:HybridSystem, options = init!(opD, HS, options_input) time_horizon = options[:T] max_jumps = options[:max_jumps] - if opC isa BFFPSV18 - inout_map = nothing - end property = options[:mode] == "check" ? options[:property] : nothing # waiting_list entries: @@ -172,10 +169,6 @@ function solve!(system::InitialValueProblem{<:HybridSystem, reach_tube = solve!(IVP(loc, set(X0)), options_copy, op=opC) - if opC isa BFFPSV18 - inout_map = reach_tube.options[:inout_map] # TODO temporary hack - end - # get the property for the current location property_loc = property isa Dict ? get(property, loc_id, nothing) : @@ -236,10 +229,6 @@ function solve!(system::InitialValueProblem{<:HybridSystem, end # Projection - if opC isa BFFPSV18 - options[:inout_map] = inout_map - end - if options[:project_reachset] || options[:projection_matrix] != nothing info("Projection...") RsetsProj = @timing project(Rsets, options)