Merge pull request #3221 from JuliaReach/schillic/3203

#3203 - Make rand of HParallelotope return a nonempty set

src/Sets/HParallelotope.jl

@@ -344,16 +344,42 @@ A random parallelotope.
 ### Notes
 
 All numbers are normally distributed with mean 0 and standard deviation 1.
+
+### Algorithm
+
+The directions matrix and offset vector are created randomly. On average there
+is a good chance that this resulting set is empty. We then modify the offset to
+ensure non-emptiness.
+
+There is a chance that the resulting set represents an unbounded set. This
+implementation checks for that case and then samples a new set.
 """
 function rand(::Type{HParallelotope};
               N::Type{<:Real}=Float64,
               dim::Int=2,
               rng::AbstractRNG=GLOBAL_RNG,
               seed::Union{Int, Nothing}=nothing)
     rng = reseed(rng, seed)
-    D = randn(rng, N, dim, dim)
-    offset = randn(rng, N, 2 * dim)
-    return HParallelotope(D, offset)
+
+    while true
+        D = randn(rng, N, dim, dim)
+        offset = randn(rng, N, 2 * dim)
+
+        # make sure that the set is not empty:
+        # `offset[i] >= -offset[i-dim]` for all `i ∈ dim+1:2*dim`
+        @inbounds for i in dim+1:2*dim
+            offset[i] = -offset[i-dim] + abs(offset[i])
+        end
+
+        P = HParallelotope(D, offset)
+
+        # convert to polyhedron to check boundedness
+        Q = HPolyhedron(vcat(P.directions, -P.directions), P.offset)
+        if isbounded(Q)
+            return P
+        end
+        # set is unbounded; sample a new set in the next iteration
+    end
 end
 
 """