Merge pull request #3308 from JuliaReach/schillic/993

#993 - More efficient support vector for UnionSet and UnionSetArray

src/LazyOperations/ConvexHullArray.jl

@@ -110,18 +110,7 @@ A support vector in the given direction.
 """
 function σ(d::AbstractVector, cha::ConvexHullArray)
     @assert !isempty(cha.array) "an empty convex hull is not allowed"
-    svec = d
-    N = eltype(d)
-    rmax = N(-Inf)
-    for chi in cha.array
-        si = σ(d, chi)
-        ri = dot(d, si)
-        if ri > rmax
-            rmax = ri
-            svec = si
-        end
-    end
-    return svec
+    return _σ_union(d, array(cha))
 end
 
 """

src/LazyOperations/UnionSet.jl

@@ -124,8 +124,8 @@ function σ(d::AbstractVector, cup::UnionSet; algorithm="support_vector")
         return dot(d, σX) > dot(d, σY) ? σX : σY
 
     elseif algorithm == "support_function"
-        m = argmax([ρ(d, X), ρ(d, Y)])
-        return m == 1 ? σ(d, X) : σ(d, Y)
+        ρX, ρY = ρ(d, X), ρ(d, Y)
+        return ρX > ρY ? σ(d, X) : σ(d, Y)
 
     else
         error("algorithm $algorithm for the support vector of a `UnionSet` is" *

src/LazyOperations/UnionSetArray.jl

@@ -109,25 +109,36 @@ support function of each set directly and then calls only the support vector of
 one of the ``Xᵢ``.
 """
 function σ(d::AbstractVector, cup::UnionSetArray; algorithm="support_vector")
-   A = array(cup)
+   arr = array(cup)
 
    if algorithm == "support_vector"
-       σarray = map(Xi -> σ(d, Xi), A)
-       ρarray = map(vi -> dot(d, vi), σarray)
-       m = argmax(ρarray)
-       return σarray[m]
+       return _σ_union(d, arr)
 
    elseif algorithm == "support_function"
-       ρarray = map(Xi -> ρ(d, Xi), A)
-       m = argmax(ρarray)
-       return σ(d, A[m])
+       m = argmax(i -> ρ(d, @inbounds arr[i]), eachindex(arr))
+       return σ(d, arr[m])
 
    else
        error("algorithm $algorithm for the support vector of a " *
              "`UnionSetArray` is unknown")
    end
 end
 
+function _σ_union(d::AbstractVector, sets)
+    σmax = d
+    N = eltype(d)
+    ρmax = N(-Inf)
+    for Xi in sets
+        σX = σ(d, Xi)
+        ρX = dot(d, σX)
+        if ρX > ρmax
+            ρmax = ρX
+            σmax = σX
+        end
+    end
+    return σmax
+end
+
 """
    ρ(d::AbstractVector, cup::UnionSetArray)
 
@@ -149,8 +160,7 @@ The support function of the union of a finite number of sets ``X₁, X₂, ...``
 can be obtained as the maximum of ``ρ(d, X₂), ρ(d, X₂), ...``.
 """
 function ρ(d::AbstractVector, cup::UnionSetArray)
-   A = array(cup)
-   return maximum(Xi -> ρ(d, Xi), A)
+   return maximum(Xi -> ρ(d, Xi), array(cup))
 end
 
 """