#668 - Replace inout_map by low-dimensional flowpipes

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

src/ReachSets/ContinuousPost/BFFPS19/reach.jl

@@ -42,7 +42,7 @@ Interface to reachability algorithms for an LTI system.
 
 ### Output
 
-A sequence of [`ReachSet`](@ref)s.
+A sequence of [`SparseReachSet`](@ref)s.
 
 ### Notes
 
@@ -53,7 +53,7 @@ is inferred from the first parameter of the system (resp. lazy set), ie.
 function reach_mixed(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}},
                               IVP{<:CLCDS{NUM}, <:LazySet{NUM}}},
                options::TwoLayerOptions
-              )::Vector{<:ReachSet} where {NUM <: Real}
+              ) where {NUM <: Real}
     # optional matrix conversion
     problem = matrix_conversion(problem, options)
 
@@ -87,16 +87,19 @@ function reach_mixed(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}},
         end
     end
 
+    # compute dimensions
+    dimensions = compute_dimensions(partition, blocks)
+
     # determine output function: linear map with the given matrix
     output_function = options[:output_function] != nothing ?
         (x -> options[:output_function] * x) :
         nothing
 
     # preallocate output vector
     if output_function == nothing
-        res_type = ReachSet{CartesianProductArray{NUM, LazySet{NUM}}}
+        res_type = SparseReachSet{CartesianProductArray{NUM, LazySet{NUM}}}
     else
-        res_type = ReachSet{Hyperrectangle{NUM}}
+        res_type = SparseReachSet{Hyperrectangle{NUM}}
     end
     res = (N == nothing) ? Vector{res_type}() : Vector{res_type}(undef, N)
 
@@ -177,6 +180,7 @@ function reach_mixed(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}},
     # add mode-specific block(s) argument
     push!(args, blocks)
     push!(args, partition)
+    push!(args, dimensions)
 
     # time step
     push!(args, options[:δ])
@@ -230,7 +234,7 @@ end
 function reach_mixed(problem::Union{IVP{<:CLCS{NUM}, <:LazySet{NUM}},
                               IVP{<:CLCCS{NUM}, <:LazySet{NUM}}},
                options::TwoLayerOptions
-              )::Vector{<:ReachSet} where {NUM <: Real}
+              ) where {NUM <: Real}
     info("Time discretization...")
     Δ = @timing discretize(problem, options[:δ],
         algorithm=options[:discretization], exp_method=options[:exp_method],

src/ReachSets/ContinuousPost/BFFPS19/reach_blocks.jl

@@ -72,13 +72,14 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
                        output_function::Union{Function, Nothing},
                        blocks::AbstractVector{Int},
                        partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
+                       dimensions::Vector{Int},
                        δ::NUM,
                        guards_proj::Vector{<:LazySet{NUM}},
                        block_options,
                        vars::Vector{Int},
                        termination::Function,
                        progress_meter::Union{Progress, Nothing},
-                       res::Vector{<:ReachSet}
+                       res::Vector{<:SparseReachSet}
                        )::Tuple{Int, Bool} where {NUM}
     ProgressMeter.update!(progress_meter, 1)
     array = CartesianProductArray(Xhat0[blocks])
@@ -88,7 +89,8 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
     t0 = zero(δ)
     t1 = δ
 
-    diff_blocks = setdiff(1:length(partition),blocks)
+    diff_blocks = setdiff(1:length(partition), blocks)
+    all_dimensions = 1:n
 
     b = length(blocks)
     bd = length(diff_blocks)
@@ -99,8 +101,8 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
             box_approximation(output_function(array_d))
 
     terminate, skip, X_store = termination_X(1, X_store, t0, X_store_d, blocks,
-                                diff_blocks, block_options, termination)
-    store!(res, 1, X_store, t0, t1, N)
+                                             diff_blocks, block_options, termination)
+    store!(res, 1, X_store, t0, t1, dimensions, N)
     if terminate
         return 1, skip
     end
@@ -127,25 +129,28 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int},
     k = 2
     @inbounds while true
         ProgressMeter.update!(progress_meter, k)
-        X_store = deco_post_sparse(b, blocks, Whatk_blocks, partition,
-                                          ϕpowerk, Xhatk, Xhat0, output_function, overapproximate)
+        X_store = deco_post_sparse(b, blocks, Whatk_blocks, partition, ϕpowerk,
+                                   Xhatk, Xhat0, output_function, overapproximate)
 
         t0 = t1
         t1 += δ
 
         terminate, skip, reach_set_intersected = termination(k, X_store, t0)
         if reach_set_intersected isa EmptySet
-            store!(res, k, X_store, t0, t1, N)
+            store!(res, k, X_store, t0, t1, dimensions, N)
             break
         end
-        if !(isdisjoint(reach_set_intersected, UnionSetArray(guards_proj)))
+        if isdisjoint(reach_set_intersected, UnionSetArray(guards_proj))
+            # keep sparse reach set
+            store!(res, k, X_store, t0, t1, all_dimensions, N)
+        else
+            # compute full-dimensional reach set
             X_store_d = deco_post_sparse(bd, diff_blocks, Whatk_diff_blocks, partition,
-                                              ϕpowerk, Xhatk_d, Xhat0, output_function, overapproximate)
+                                         ϕpowerk, Xhatk_d, Xhat0, output_function, overapproximate)
             X_store = getX_store(reach_set_intersected, X_store_d, block_options, blocks, diff_blocks)
+            store!(res, k, X_store, t0, t1, dimensions, N)
         end
 
-        store!(res, k, X_store, t0, t1, N)
-
         if terminate
             break
         end
@@ -181,13 +186,14 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM},
                        output_function::Union{Function, Nothing},
                        blocks::AbstractVector{Int},
                        partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
+                       dimensions::Vector{Int},
                        δ::NUM,
                        guards_proj::Vector{<:LazySet{NUM}},
                        block_options,
                        vars::Vector{Int},
                        termination::Function,
                        progress_meter::Union{Progress, Nothing},
-                       res::Vector{<:ReachSet}
+                       res::Vector{<:SparseReachSet}
                        )::Tuple{Int, Bool} where {NUM}
     ProgressMeter.update!(progress_meter, 1)
     array = CartesianProductArray(Xhat0[blocks])
@@ -197,7 +203,8 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM},
     t0 = zero(δ)
     t1 = δ
 
-    diff_blocks = setdiff(1:length(partition),blocks)
+    diff_blocks = setdiff(1:length(partition), blocks)
+    all_dimensions = 1:n
 
     b = length(blocks)
     bd = length(diff_blocks)
@@ -208,8 +215,8 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM},
         box_approximation(output_function(array_d))
 
     terminate, skip, X_store = termination_X(1, X_store, t0, X_store_d, blocks,
-                                diff_blocks, block_options, termination)
-    store!(res, 1, X_store, t0, t1, N)
+                                             diff_blocks, block_options, termination)
+    store!(res, 1, X_store, t0, t1, dimensions, N)
     if terminate
         return 1, skip
     end
@@ -239,23 +246,27 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM},
     @inbounds while true
         ProgressMeter.update!(progress_meter, k)
         X_store = deco_post_dense(b, blocks, Whatk_blocks, partition, ϕpowerk,
-                    arr, arr_length, U, Xhat0, Xhatk, output_function, overapproximate)
+                                  arr, arr_length, U, Xhat0, Xhatk, output_function, overapproximate)
         t0 = t1
         t1 += δ
         terminate, skip, reach_set_intersected = termination(k, X_store, t0)
         if reach_set_intersected isa EmptySet
-            store!(res, k, X_store, t0, t1, N)
+            store!(res, k, X_store, t0, t1, dimensions, N)
             break
         end
-        if !(isdisjoint(reach_set_intersected, UnionSetArray(guards_proj)))
+
+        if isdisjoint(reach_set_intersected, UnionSetArray(guards_proj))
+            # keep sparse reach set
+            store!(res, k, X_store, t0, t1, all_dimensions, N)
+        else
+            # compute full-dimensional reach set
             X_store_d = deco_post_dense(bd, diff_blocks, Whatk_diff_blocks, partition, ϕpowerk, arr,
-                                              arr_length, U, Xhat0, Xhatk_d, output_function, overapproximate)
+                                        arr_length, U, Xhat0, Xhatk_d, output_function, overapproximate)
 
             X_store = getX_store(reach_set_intersected, X_store_d, block_options, blocks, diff_blocks)
+            store!(res, k, X_store, t0, t1, dimensions, N)
         end
 
-        store!(res, k, X_store, t0, t1, N)
-
         if terminate
             break
         end

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

src/ReachSets/ContinuousPost/BFFPSV18/compute_dimensions.jl

@@ -0,0 +1,15 @@
+function compute_dimensions(partition, blocks)
+    dimensions = Vector{Int}()
+    dims = 0
+    next_idx = 1
+    for block_idx in blocks
+        while next_idx < block_idx
+            # sum up dimensions of skipped blocks
+            dims += length(partition[next_idx])
+            next_idx += 1
+        end
+        append!(dimensions, partition[next_idx])
+        next_idx += 1
+    end
+    return dimensions
+end

src/ReachSets/ContinuousPost/BFFPSV18/reach.jl

@@ -17,7 +17,7 @@ Interface to reachability algorithms for an LTI system.
 
 ### Output
 
-A sequence of [`ReachSet`](@ref)s.
+A sequence of [`SparseReachSet`](@ref)s.
 
 ### Notes
 
@@ -31,7 +31,7 @@ is inferred from the first parameter of the system (resp. lazy set), ie.
 function reach(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}},
                               IVP{<:CLCDS{NUM}, <:LazySet{NUM}}},
                options::TwoLayerOptions
-              )::Vector{<:ReachSet} where {NUM <: Real}
+              ) where {NUM <: Real}
     # optional matrix conversion
     problem = matrix_conversion(problem, options)
 
@@ -66,16 +66,19 @@ function reach(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}},
         end
     end
 
+    # compute dimensions
+    dimensions = compute_dimensions(partition, blocks)
+
     # determine output function: linear map with the given matrix
     output_function = options[:output_function] != nothing ?
         (x -> options[:output_function] * x) :
         nothing
 
     # preallocate output vector
     if output_function == nothing
-        res_type = ReachSet{CartesianProductArray{NUM, LazySet{NUM}}}
+        res_type = SparseReachSet{CartesianProductArray{NUM, LazySet{NUM}}}
     else
-        res_type = ReachSet{Hyperrectangle{NUM}}
+        res_type = SparseReachSet{Hyperrectangle{NUM}}
     end
     res = (N == nothing) ? Vector{res_type}() : Vector{res_type}(undef, N)
 
@@ -156,6 +159,7 @@ function reach(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}},
     # add mode-specific block(s) argument
     push!(args, blocks)
     push!(args, partition)
+    push!(args, dimensions)
 
     # time step
     push!(args, options[:δ])
@@ -204,7 +208,7 @@ end
 function reach(problem::Union{IVP{<:CLCS{NUM}, <:LazySet{NUM}},
                               IVP{<:CLCCS{NUM}, <:LazySet{NUM}}},
                options::TwoLayerOptions
-              )::Vector{<:ReachSet} where {NUM <: Real}
+              ) where {NUM <: Real}
     info("Time discretization...")
     Δ = @timing discretize(problem, options[:δ],
         algorithm=options[:discretization], exp_method=options[:exp_method],