Speed up with MutableArithmetics

perf/runbenchmarks.jl

@@ -39,24 +39,33 @@ using StaticArrays
 #   samples:          9809
 #   evals/sample:     1
 
-let
-    h = hrep([HalfSpace([ 1.,  1], 1),
-              HalfSpace([ 1., -1], 0),
-              HalfSpace([-1.,  0], 0)])
+using BenchmarkTools
+using Polyhedra
+
+function vector(T)
+    h = hrep([HalfSpace(T[ 1,  1], T(1)),
+              HalfSpace(T[ 1, -1], T(0)),
+              HalfSpace(T[-1,  0], T(0))])
     display(@benchmark(doubledescription($h)))
 end
 
-let
-    h = hrep([HalfSpace((@SVector [ 1,  1]), 1),
-              HalfSpace((@SVector [ 1, -1]), 0),
-              HalfSpace((@SVector [-1,  0]), 0)])
+using StaticArrays
+function static(T)
+    h = hrep([HalfSpace((@SVector T[ 1,  1]), T(1)),
+              HalfSpace((@SVector T[ 1, -1]), T(0)),
+              HalfSpace((@SVector T[-1,  0]), T(0))])
     display(@benchmark(doubledescription($h)))
+    @profview_allocs for _ in 1:10000
+        doubledescription(h)
+    end
 end
 
-let
-    h = hrep([ 1  1
-               1 -1
-              -1  0],
-             [1., 0, 0])
+function matrix(T)
+    h = hrep(T[
+        1  1
+        1 -1
+       -1  0],
+        T[1, 0, 0],
+    )
     display(@benchmark(doubledescription($h)))
 end

src/Polyhedra.jl

@@ -9,8 +9,7 @@ using LinearAlgebra
 #      like RowEchelon.jl (slow) or LU (faster).
 import GenericLinearAlgebra
 
-import MutableArithmetics
-const MA = MutableArithmetics
+import MutableArithmetics as MA
 
 import MathOptInterface as MOI
 import MathOptInterface.Utilities as MOIU

src/center.jl

@@ -19,7 +19,7 @@ function maximum_radius_with_center(h::HRep{T}, center) where T
                     return zero(T)
                 end
             else
-                new_radius = (hs.β - center ⋅ hs.a) / n
+                new_radius = (hs.β - _dot(center, hs.a)) / n
                 if radius === nothing
                     radius = new_radius
                 else

src/comp.jl

@@ -1,27 +1,34 @@
+# We use `Zero` so that it is allocation-free, while `zero(Rational{BigInt})` would be costly
+_default_tol(::Type{<:Union{Integer, Rational}}) = MA.Zero()
+# `rtoldefault(BigFloat)` is quite slow and allocates a lot so we hope that
+# `tol` will be called at some upper level function and passed in here instead
+# of created everytime
+_default_tol(::Type{T}) where {T<:AbstractFloat} = Base.rtoldefault(T)
+
 isapproxzero(x::T; kws...) where {T<:Real} = x == zero(T)
-isapproxzero(x::T; ztol=Base.rtoldefault(T)) where {T<:AbstractFloat} = abs(x) < ztol
-isapproxzero(x::AbstractVector{T}; kws...) where {T<:Real} = isapproxzero(maximum(abs.(x)); kws...)
+isapproxzero(x::T; tol=Base.rtoldefault(T)) where {T<:AbstractFloat} = abs(x) < tol
+isapproxzero(x::AbstractVector{T}; kws...) where {T<:Real} = isapproxzero(maximum(abs, x); kws...)
 #isapproxzero(x::Union{SymPoint, Ray, Line}; kws...) = isapproxzero(coord(x); kws...)
 isapproxzero(x::Union{Ray, Line}; kws...) = isapproxzero(coord(x); kws...)
 isapproxzero(h::HRepElement; kws...) = isapproxzero(h.a; kws...) && isapproxzero(h.β; kws...)
 
-_isapprox(x::Union{T, AbstractArray{T}}, y::Union{T, AbstractArray{T}}) where {T<:Union{Integer, Rational}} = x == y
+_isapprox(x::Union{T, AbstractArray{T}}, y::Union{T, AbstractArray{T}}; tol) where {T<:Union{Integer, Rational}} = x == y
 # I check with zero because isapprox(0, 1e-100) is false...
 # but isapprox(1e-100, 2e-100) should be false
-_isapprox(x, y) = (isapproxzero(x) ? isapproxzero(y) : (isapproxzero(y) ? isapproxzero(x) : isapprox(x, y)))
+_isapprox(x, y; tol) = (isapproxzero(x; tol) ? isapproxzero(y; tol) : (isapproxzero(y; tol) ? isapproxzero(x; tol) : isapprox(x, y; rtol = tol)))
 
 #Base.isapprox(r::SymPoint, s::SymPoint) = _isapprox(coord(r), coord(s)) || _isapprox(coord(r), -coord(s))
 function _scaleray(r::Union{Line, Ray}, s::Union{Line, Ray})
     cr = coord(r)
     cs = coord(s)
     cr * sum(abs.(cs)), cs * sum(abs.(cr))
 end
-function Base.isapprox(r::Line, s::Line)
+function _isapprox(r::Line, s::Line; tol)
     rs, ss = _scaleray(r, s)
-    _isapprox(rs, ss) || _isapprox(rs, -ss)
+    _isapprox(rs, ss; tol) || _isapprox(rs, -ss; tol)
 end
-function Base.isapprox(r::Ray, s::Ray)
-    _isapprox(_scaleray(r, s)...)
+function _isapprox(r::Ray, s::Ray; tol)
+    _isapprox(_scaleray(r, s)...; tol)
 end
 # TODO check that Vec in GeometryTypes also does that
 function _scalehp(h1, h2)
@@ -33,28 +40,28 @@ function Base.:(==)(h1::HyperPlane, h2::HyperPlane)
     (a1, β1), (a2, β2) = _scalehp(h1, h2)
     (a1 == a2 && β1 == β2) || (a1 == -a2 && β1 == -β2)
 end
-function Base.isapprox(h1::HyperPlane, h2::HyperPlane)
+function _isapprox(h1::HyperPlane, h2::HyperPlane; tol)
     (a1, β1), (a2, β2) = _scalehp(h1, h2)
-    (_isapprox(a1, a2) && _isapprox(β1, β2)) || (_isapprox(a1, -a2) && _isapprox(β1, -β2))
+    (_isapprox(a1, a2; tol) && _isapprox(β1, β2; tol)) || (_isapprox(a1, -a2; tol) && _isapprox(β1, -β2; tol))
 end
 function Base.:(==)(h1::HalfSpace, h2::HalfSpace)
     (a1, β1), (a2, β2) = _scalehp(h1, h2)
     a1 == a2 && β1 == β2
 end
-function Base.isapprox(h1::HalfSpace, h2::HalfSpace)
+function _isapprox(h1::HalfSpace, h2::HalfSpace; tol)
     (a1, β1), (a2, β2) = _scalehp(h1, h2)
-    _isapprox(a1, a2) && _isapprox(β1, β2)
+    _isapprox(a1, a2; tol) && _isapprox(β1, β2; tol)
 end
 
-_lt(x::T, y::T) where {T<:Real} = x < y
-_lt(x::T, y::T) where {T<:AbstractFloat} = x < y && !_isapprox(x, y)
-_lt(x::S, y::T) where {S<:Real,T<:Real} = _lt(promote(x, y)...)
-_gt(x::S, y::T) where {S<:Real, T<:Real} = _lt(y, x)
-_leq(x::T, y::T) where {T<:Real} = x <= y
-_leq(x::T, y::T) where {T<:AbstractFloat} = x <= y || _isapprox(x, y)
-_leq(x::S, y::T) where {S<:Real,T<:Real} = _leq(promote(x, y)...)
-_geq(x::T, y::T) where {T<:Real} = _leq(y, x)
-_pos(x::T) where {T<:Real} = _gt(x, zero(T))
-_neg(x::T) where {T<:Real} = _lt(x, zero(T))
-_nonneg(x::T) where {T<:Real} = _geq(x, zero(T))
-_nonpos(x::T) where {T<:Real} = _leq(x, zero(T))
+_lt(x::T, y::T; tol) where {T<:Real} = x < y
+_lt(x::T, y::T; tol) where {T<:AbstractFloat} = x < y && !_isapprox(x, y; tol)
+_lt(x::S, y::T; tol) where {S<:Real,T<:Real} = _lt(promote(x, y)...; tol)
+_gt(x::S, y::T; tol) where {S<:Real, T<:Real} = _lt(y, x; tol)
+_leq(x::T, y::T; tol) where {T<:Real} = x <= y
+_leq(x::T, y::T; tol) where {T<:AbstractFloat} = x <= y || _isapprox(x, y; tol)
+_leq(x::S, y::T; tol) where {S<:Real,T<:Real} = _leq(promote(x, y)...; tol)
+_geq(x::T, y::T; tol) where {T<:Real} = _leq(y, x; tol)
+_pos(x::T; tol) where {T<:Real} = _gt(x, zero(T); tol)
+_neg(x::T; tol) where {T<:Real} = _lt(x, zero(T); tol)
+_nonneg(x::T; tol) where {T<:Real} = _geq(x, zero(T); tol)
+_nonpos(x::T; tol) where {T<:Real} = _leq(x, zero(T); tol)

src/doubledescription.jl

@@ -19,7 +19,7 @@ function dualfullspace(rep::Representation{T}) where T
 end
 
 """
-    doubledescription(h::HRepresentation)
+    doubledescription(h::HRepresentation; tol)
 
 Computes the V-representation of the polyhedron represented by `h` using the Double-Description algorithm [MRTT53, FP96].
 It maintains a list of the points, rays and lines that are in the resulting representation
@@ -305,44 +305,44 @@ function set_in!(data, I, index)
         end
     end
 end
-function add_element!(data, k, el, tight)
+function add_element!(data, k, el, tight; tol)
     cutoff = nothing
     for i in reverse(1:k)
         if data.cutline[i] !== nothing
             el = line_project(el, data.cutline[i], data.halfspaces[i])
             # The line is in all halfspaces from `i+1` up so projecting with it does not change it.
             push!(tight, i)
-            return add_adjacent_element!(data, i, k, el, data.lineray[i], tight)
+            return add_adjacent_element!(data, i, k, el, data.lineray[i], tight; tol)
         end
         # Could avoid computing `dot` twice between `el` and the halfspace here.
-        if !(el in data.halfspaces[i])
+        if !in(el, data.halfspaces[i]; tol)
             cutoff = i
             break
         end
-        if el in hyperplane(data.halfspaces[i])
+        if in(el, hyperplane(data.halfspaces[i]); tol)
             push!(tight, i)
         end
     end
     index = add_index!(data, cutoff, el, tight)
     set_in!(data, (index.cutoff + 1):k, index)
     return index
 end
-function project_onto_affspace(data, offset, el, hyperplanes)
+function project_onto_affspace(data, offset, el, hyperplanes; tol)
     for i in reverse(eachindex(hyperplanes))
         line = data.cutline[offset + i]
         h = hyperplanes[i]
         if line !== nothing
             el = line_project(el, line, h)
-        elseif !(el in h)
+        elseif !in(el, h; tol)
             # e.g. 0x1 + 0x2 = 1 or -1 = x1 = 1, ...
             return nothing
         end
     end
     return el
 end
-function add_adjacent_element!(data, i, k, el, parent, tight)
-    index = add_element!(data, i - 1, el, tight)
-    add_intersection!(data, index, parent, nothing, i - 1)
+function add_adjacent_element!(data, i, k, el, parent, tight; tol)
+    index = add_element!(data, i - 1, el, tight; tol)
+    add_intersection!(data, index, parent, nothing, i - 1; tol)
     set_in!(data, i:k, index)
     return index
 end
@@ -377,11 +377,11 @@ Combine `el1` and `el2` which are on two different sides of the halfspace `h` so
 that the combination is in the corresponding hyperplane.
 """
 function combine(h::HalfSpace, el1::VRepElement, el2::VRepElement)
-    combine(h.β, el1, h.a ⋅ el1, el2, h.a ⋅ el2)
+    combine(h.β, el1, _dot(h.a, el1), el2, _dot(h.a, el2))
 end
 
 """
-    add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1)
+    add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1; tol)
 
 Add the intersection of the elements of indices `idx1` and `idx2` that belong to
 the hyperplane corresponding the halfspace of index `hp_idx` (see 3.2 (ii) of
@@ -396,9 +396,9 @@ to avoid adding redundant elements.
 Double description method revisited
 *Combinatorics and computer science*, *Springer*, **1996**, 91-111
 """
-function add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1)
+function add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1; tol)
     if idx1.cutoff > idx2.cutoff
-        return add_intersection!(data, idx2, idx1, hp_idx, hs_idx)
+        return add_intersection!(data, idx2, idx1, hp_idx, hs_idx; tol)
     end
     i = idx2.cutoff
     # Condition (c_k) in [FP96]
@@ -418,31 +418,35 @@ function add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1)
     #      in case `rj, rj'` have inherited adjacency
     tight = tight_halfspace_indices(data, idx1) ∩ tight_halfspace_indices(data, idx2)
     push!(tight, i)
-    return add_adjacent_element!(data, i, i, newel, idx1, tight)
+    return add_adjacent_element!(data, i, i, newel, idx1, tight; tol)
 end
 
+# `MA.operate` is going to redirect to an implementation that
+# exploits the mutable API of `BigFloat` of `BigInt`
+_dot(a, b) = MA.operate(LinearAlgebra.dot, a, b)
+
 _shift(el::AbstractVector, line::Line) = el + Polyhedra.coord(line)
 _shift(el::Line, line::Line) = el + line
 _shift(el::Ray, line::Line) = el + Polyhedra.Ray(Polyhedra.coord(line))
 function _λ_proj(r::VStruct, line::Line, h::HRepElement)
     # Line or ray `r`:
     # (r + λ * line) ⋅ h.a == 0
     # λ = -r ⋅ h.a / (line ⋅ h.a)
-    return -r ⋅ h.a / (line ⋅ h.a)
+    return -_dot(r, h.a) / _dot(line, h.a)
 end
 function _λ_proj(x::AbstractVector, line::Line, h::HRepElement)
     # Point `x`:
     # (x + λ * line) ⋅ h.a == h.β
     # λ = (h.β - x ⋅ h.a) / (line ⋅ h.a)
-    return (h.β - x ⋅ h.a) / (line ⋅ h.a)
+    return (h.β - _dot(x, h.a)) / _dot(line, h.a)
 end
 function line_project(el, line, h)
     λ = _λ_proj(el, line, h)
     return Polyhedra.simplify(_shift(el, λ * line))
 end
-function hline(data, line::Line, i, h)
-    value = h.a ⋅ line
-    if !Polyhedra.isapproxzero(value)
+function hline(data, line::Line, i, h; tol)
+    value = _dot(h.a, line)
+    if !Polyhedra.isapproxzero(value; tol)
         if data.cutline[i] === nothing
             # Make `lineray` point inward
             data.cutline[i] = value > 0 ? -line : line
@@ -456,20 +460,20 @@ function hline(data, line::Line, i, h)
 end
 
 # TODO remove solver arg `_`, it is kept to avoid breaking code
-function doubledescription(hr::HRepresentation, _ = nothing)
+function doubledescription(hr::HRepresentation{T}, solver = nothing; tol = _default_tol(T)) where {T}
     v = Polyhedra.dualfullspace(hr)
     hps = Polyhedra.lazy_collect(hyperplanes(hr))
     hss = Polyhedra.lazy_collect(halfspaces(hr))
     data = DoubleDescriptionData{pointtype(v), raytype(v), linetype(v)}(fulldim(hr), hps, hss)
     for line in lines(v)
         cut = false
         for i in reverse(eachindex(hps))
-            cut, line = hline(data, line, nhalfspaces(hr) + i, hps[i])
+            cut, line = hline(data, line, nhalfspaces(hr) + i, hps[i]; tol)
             cut && break
         end
         if !cut
             for i in reverse(eachindex(hss))
-                cut, line = hline(data, line, i, hss[i])
+                cut, line = hline(data, line, i, hss[i]; tol)
                 cut && break
             end
         end
@@ -484,14 +488,14 @@ function doubledescription(hr::HRepresentation, _ = nothing)
         line = data.cutline[i]
         if line !== nothing
             ray = Polyhedra.Ray(Polyhedra.coord(line))
-            data.lineray[i] = add_element!(data, i - 1, ray, BitSet((i + 1):length(hss)))
+            data.lineray[i] = add_element!(data, i - 1, ray, BitSet((i + 1):length(hss)); tol)
         end
     end
     @assert isone(npoints(v))
     # Add the origin
-    orig = project_onto_affspace(data, nhalfspaces(hr), first(points(v)), hps)
+    orig = project_onto_affspace(data, nhalfspaces(hr), first(points(v)), hps; tol)
     if orig !== nothing
-        add_element!(data, nhalfspaces(hr), orig, resized_bitset(data))
+        add_element!(data, nhalfspaces(hr), orig, resized_bitset(data); tol)
     end
     for i in reverse(eachindex(hss))
         if data.cutline[i] === nothing && isempty(data.cutpoints[i]) && isempty(data.cutrays[i])
@@ -507,11 +511,11 @@ function doubledescription(hr::HRepresentation, _ = nothing)
             # Catches new adjacent rays, see 3.2 (ii) of [FP96]
             for p1 in data.pin[i], p2 in data.pin[i]
                 if p1.cutoff < p2.cutoff
-                    add_intersection!(data, p1, p2, i)
+                    add_intersection!(data, p1, p2, i; tol)
                 end
             end
             for p in data.pin[i], r in data.rin[i]
-                add_intersection!(data, p, r, i)
+                add_intersection!(data, p, r, i; tol)
             end
         end
         deleteat!(data.cutpoints, i)
@@ -521,7 +525,7 @@ function doubledescription(hr::HRepresentation, _ = nothing)
                 # We encounter both `r1, r2` and `r2, r1`.
                 # Break this symmetry with:
                 if r1.cutoff < r2.cutoff
-                    add_intersection!(data, r1, r2, i)
+                    add_intersection!(data, r1, r2, i; tol)
                 end
             end
         end