Merge pull request #274 from JuliaReach/schillic/272

#272 - Lazy discrete post operator

src/ReachSets/DiscretePost/DiscretePost.jl

@@ -5,12 +5,34 @@ Abstract supertype of all discrete post operators.
 
 ### Notes
 
-All discrete post operators should provide the following methods, in addition
+All discrete post operators should provide the following method, in addition
 to those provided for general post operators:
 ```julia
 tube⋂inv(op::DiscretePost, reach_tube, invariant, Rsets, start_interval)
-cluster(op::DiscretePost, reach_sets)
-isfixpoint(op::DiscretePost, reach_set, passed_list, loc_id)
 ```
 """
 abstract type DiscretePost <: PostOperator end
+
+function cluster(op::DiscretePost, reach_sets)
+    # TODO apply some clustering
+    return reach_sets
+end
+
+function isfixpoint(op::DiscretePost,
+                    reach_set::ReachSet{LazySet{N}, N},
+                    passed_list,
+                    loc_id
+                   ) where N
+    if isassigned(passed_list, loc_id)
+        for other_reach_set in passed_list[loc_id]
+            if reach_set.X ⊆ other_reach_set.X
+                info("found a fixpoint in some reach tube")
+                return true
+            end
+        end
+        return false
+    else
+        passed_list[loc_id] = Vector{ReachSet{LazySet{N}, N}}()
+        return false
+    end
+end

src/ReachSets/DiscretePost/LazyTextbookDiscretePost.jl

@@ -0,0 +1,112 @@
+# ==============================================================================
+# Textbook implementation of a discrete post operator, but with lazy operations.
+# ==============================================================================
+
+struct LazyTextbookDiscretePost <: DiscretePost
+    options::Options
+end
+
+LazyTextbookDiscretePost() = LazyTextbookDiscretePost(Options())
+
+function init(op::LazyTextbookDiscretePost, system, options_input)
+    options_input.dict[:n] = statedim(system, 1)
+
+    # solver-specific options (adds default values for unspecified options)
+    options = validate_solver_options_and_add_default_values!(options_input)
+
+    # Input -> Output variable mapping
+    options.dict[:inout_map] =
+        inout_map_reach(options[:partition], options[:blocks], options[:n])
+
+    # set up operator-specific options
+    # TODO temporarily use box approximation
+    @assert !haskey(op.options.dict, :overapproximation)
+    op.options.dict[:overapproximation] = Hyperrectangle
+
+    return options
+end
+
+function tube⋂inv!(op::LazyTextbookDiscretePost,
+                   reach_tube::Vector{<:ReachSet{<:LazySet{N}}},
+                   invariant,
+                   Rsets,
+                   start_interval
+                  ) where {N}
+    intersections = Vector{ReachSet{LazySet{N}, N}}()
+    for reach_set in reach_tube
+        R⋂I = Intersection(invariant, reach_set.X)
+        if isempty(R⋂I)
+            break
+        end
+        push!(intersections, ReachSet{LazySet{N}, N}(R⋂I,
+            reach_set.t_start + start_interval[1],
+            reach_set.t_end + start_interval[2]))
+    end
+
+    append!(Rsets, intersections)
+end
+
+function post(op::LazyTextbookDiscretePost,
+              system::InitialValueProblem{<:HybridSystem, <:LazySet{N}},
+              waiting_list,
+              passed_list,
+              source_loc_id,
+              tube⋂inv,
+              jumps
+             ) where {N}
+    HS = system.s
+    jumps += 1
+    for trans in out_transitions(HS, source_loc_id)
+        info("Considering transition: $trans")
+        target_loc_id = target(HS, trans)
+        target_loc = HS.modes[target(HS, trans)]
+        target_invariant = target_loc.X
+        trans_annot = HS.resetmaps[symbol(HS, trans)]
+        guard = trans_annot.X
+        assignment = trans_annot.A
+
+        # perform jumps
+        post_jump = Vector{ReachSet{LazySet{N}, N}}()
+        sizehint!(post_jump, length(tube⋂inv))
+        for reach_set in tube⋂inv
+            # check intersection with guard
+            R⋂G = Intersection(reach_set.X, guard)
+            if isempty(R⋂G)
+                continue
+            end
+            # apply assignment
+            A⌜R⋂G⌟ = LinearMap(assignment, R⋂G)
+            # intersect with target invariant
+            A⌜R⋂G⌟⋂I = Intersection(target_invariant, A⌜R⋂G⌟)
+            # overapproximate final set
+            res = overapproximate(A⌜R⋂G⌟⋂I, op.options[:overapproximation];
+                                  upper_bound=true)
+
+            # check if the final set is empty
+            if isempty(res)
+                continue
+            end
+
+            # store result
+            push!(post_jump, ReachSet{LazySet{N}, N}(res,
+                                                     reach_set.t_start,
+                                                     reach_set.t_end))
+        end
+
+        if (isempty(post_jump))
+            # transition can never be taken
+            continue
+        end
+
+        # apply clustering
+        clustered = cluster(op, post_jump)
+
+        # push new sets after jump (unless a fixpoint is detected)
+        for reach_set in clustered
+            if !isfixpoint(op, reach_set, passed_list, target_loc_id)
+                push!(passed_list[target_loc_id], reach_set)
+                push!(waiting_list, (target(HS, trans), reach_set, jumps))
+            end
+        end
+    end
+end

src/ReachSets/DiscretePost/TextbookDiscretePost.jl

@@ -60,20 +60,19 @@ function tube⋂inv!(op::TextbookDiscretePost,
         rs = reach_set.X
         # TODO temporary workaround for 1D sets
         if dim(rs) == 1
-            reach_tube = VPolytope(vertices_list(
+            rs_converted = VPolytope(vertices_list(
                 Approximations.overapproximate(rs, LazySets.Interval)))
-        # TODO offer a lazy intersection here
         # TODO offer more options instead of taking the VPolytope intersection
         elseif rs isa CartesianProductArray
-            reach_tube = VPolytope(vertices_list(rs))
+            rs_converted = VPolytope(vertices_list(rs))
         else
             error("unsupported set type for reach tube: $(typeof(rs))")
         end
-        intersect = intersection(invariant, reach_tube)
-        if isempty(intersect)
+        R⋂I = intersection(invariant, rs_converted)
+        if isempty(R⋂I)
             break
         end
-        push!(intersections, ReachSet{LazySet{N}, N}(intersect,
+        push!(intersections, ReachSet{LazySet{N}, N}(R⋂I,
             reach_set.t_start + start_interval[1],
             reach_set.t_end + start_interval[2]))
     end
@@ -113,21 +112,22 @@ function post(op::TextbookDiscretePost,
         sizehint!(post_jump, length(tube⋂inv))
         for reach_set in tube⋂inv
             # check intersection with guard
-            cap = intersection(guard, VPolytope(vertices_list(reach_set.X)))
-            if isempty(cap)
+            R⋂G = intersection(guard, VPolytope(vertices_list(reach_set.X)))
+            if isempty(R⋂G)
                 continue
             end
             # apply assignment
-            asgn = linear_map(assignment, cap)
+            A⌜R⋂G⌟ = linear_map(assignment, R⋂G)
             # intersect with target invariant
-            res = intersection(target_invariant, asgn)
-            if isempty(res)
+            A⌜R⋂G⌟⋂I = intersection(target_invariant, A⌜R⋂G⌟)
+            if isempty(A⌜R⋂G⌟⋂I)
                 continue
             end
 
             # store result
-            push!(post_jump, ReachSet{LazySet{N}, N}(res, reach_set.t_start,
-                                                      reach_set.t_end))
+            push!(post_jump, ReachSet{LazySet{N}, N}(A⌜R⋂G⌟⋂I,
+                                                     reach_set.t_start,
+                                                     reach_set.t_end))
         end
 
         if (isempty(post_jump))
@@ -150,27 +150,3 @@ function post(op::TextbookDiscretePost,
         end
     end
 end
-
-function cluster(op::TextbookDiscretePost, reach_sets)
-    # TODO apply some clustering
-    return reach_sets
-end
-
-function isfixpoint(op::TextbookDiscretePost,
-                    reach_set::ReachSet{LazySet{N}, N},
-                    passed_list,
-                    loc_id
-                   ) where N
-    if isassigned(passed_list, loc_id)
-        for other_reach_set in passed_list[loc_id]
-            if reach_set.X ⊆ other_reach_set.X
-                info("found a fixpoint in some reach tube")
-                return true
-            end
-        end
-        return false
-    else
-        passed_list[loc_id] = Vector{ReachSet{LazySet{N}, N}}()
-        return false
-    end
-end

src/ReachSets/ReachSets.jl

@@ -202,8 +202,10 @@ include("ContinuousPost/BFFPSV18.jl")
 export BFFPSV18
 
 include("DiscretePost/TextbookDiscretePost.jl")
+include("DiscretePost/LazyTextbookDiscretePost.jl")
 
-export TextbookDiscretePost
+export TextbookDiscretePost,
+       LazyTextbookDiscretePost
 
 # =========================
 # External reach interface

test/REQUIRE

@@ -1,3 +1,4 @@
 CDDLib
 Compat
+Optim
 Polyhedra

test/Reachability/unit_solve_hybrid.jl

@@ -1,2 +1,4 @@
+using Optim
+
 include("unit_solve_hybrid_bouncing_ball.jl")
 include("unit_solve_hybrid_thermostat.jl")

test/Reachability/unit_solve_hybrid_bouncing_ball.jl

@@ -37,4 +37,8 @@ input_options = Options(:mode=>"reach");
 problem_options = Options(:vars=>[1,2], :T=>5.0, :δ=>0.1, :plot_vars=>[1, 2],
                           :max_jumps=>1, :verbosity=>1);
 options_input = merge(problem_options, input_options);
+options_copy = Options(copy(options_input.dict))
 sol = solve(system, options_input);
+
+sol = solve(system, options_copy, Reachability.BFFPSV18(),
+            Reachability.ReachSets.LazyTextbookDiscretePost());