Speed up with MutableArithmetics (#343)

* Speed up with MutableArithmetics

* Fix tol

* Unifty tol

* Fixes

* Add missing tol

* Fixes

* Bump Julia to v1.10

* Fix doctest output

.github/workflows/ci.yml

@@ -24,7 +24,7 @@ jobs:
           - version: '1'
             os: ubuntu-latest
             arch: x86
-          - version: '1.6'
+          - version: '1.10'
             os: ubuntu-latest
             arch: x64
     steps:

Project.toml

@@ -27,11 +27,11 @@ PolyhedraRecipesBaseExt = "RecipesBase"
 [compat]
 GenericLinearAlgebra = "0.2, 0.3"
 GeometryBasics = "0.2, 0.3, 0.4, 0.5"
-JuMP = "0.23, 1"
+JuMP = "1"
 LinearAlgebra = "1.6"
 MathOptInterface = "1"
 MutableArithmetics = "1"
-RecipesBase = "0.7, 0.8, 1.0"
+RecipesBase = "1.0"
 SparseArrays = "1.6"
-StaticArrays = "0.12, 1.0"
-julia = "1.6"
+StaticArrays = "1.0"
+julia = "1.10"

docs/src/optimization.md

@@ -40,11 +40,11 @@ julia> using JuMP
 
 julia> model = Model()
 A JuMP Model
-Feasibility problem with:
-Variables: 0
-Model mode: AUTOMATIC
-CachingOptimizer state: NO_OPTIMIZER
-Solver name: No optimizer attached.
+├ solver: none
+├ objective_sense: FEASIBILITY_SENSE
+├ num_variables: 0
+├ num_constraints: 0
+└ Names registered in the model: none
 
 julia> @variable(model, λ[1:2])
 2-element Vector{VariableRef}:

ext/PolyhedraGeometryBasicsExt.jl

@@ -90,15 +90,15 @@ function _isdup(zray, triangles)
 end
 _isdup(poly, hidx, triangles) = _isdup(get(poly, hidx).a, triangles)
 
-function decompose_plane!(triangles::Vector, d::FullDim, zray, incident_points, incident_lines, incident_rays, exit_point::Function, counterclockwise::Function, rotate::Function)
+function decompose_plane!(triangles::Vector, d::FullDim, zray, incident_points, incident_lines, incident_rays, exit_point::Function, counterclockwise::Function, rotate::Function; tol)
     # xray should be the rightmost ray
     xray = nothing
     # xray should be the leftmost ray
     yray = nothing
     isapproxzero(zray) && return
 
     # Checking rays
-    hull, lines, rays = _planar_hull(d, incident_points, incident_lines, incident_rays, counterclockwise, rotate)
+    hull, lines, rays = _planar_hull(d, incident_points, incident_lines, incident_rays, counterclockwise, rotate; tol)
     if isempty(lines)
         if length(hull) + length(rays) < 3
             return
@@ -190,24 +190,24 @@ function fulldecompose(triangles::Vector{_Tri{N,T}}) where {N,T}
     return pts, faces, ns
 end
 
-function fulldecompose(poly_geom::Mesh, ::Type{T}) where T
+function fulldecompose(poly_geom::Mesh, ::Type{T}; tol = Polyhedra._default_tol(T)) where T
     poly = poly_geom.polyhedron
     exit_point = scene(poly_geom, T)
     triangles = _Tri{2,T}[]
     z = StaticArrays.SVector(zero(T), zero(T), one(T))
-    decompose_plane!(triangles, FullDim(poly), z, collect(points(poly)), lines(poly), rays(poly), exit_point, counterclockwise, rotate)
+    decompose_plane!(triangles, FullDim(poly), z, collect(points(poly)), lines(poly), rays(poly), exit_point, counterclockwise, rotate; tol)
     return fulldecompose(triangles)
 end
 
-function fulldecompose(poly_geom::Mesh{3}, ::Type{T}) where T
+function fulldecompose(poly_geom::Mesh{3}, ::Type{T}; tol = Polyhedra._default_tol(T)) where T
     poly = poly_geom.polyhedron
     exit_point = scene(poly_geom, T)
     triangles = _Tri{3,T}[]
     function decompose_plane(hidx)
         zray = get(poly, hidx).a
         counterclockwise(a, b) = dot(cross(a, b), zray)
         rotate(r) = cross(zray, r)
-        decompose_plane!(triangles, FullDim(poly), zray, incidentpoints(poly, hidx), incidentlines(poly, hidx), incidentrays(poly, hidx), exit_point, counterclockwise, rotate)
+        decompose_plane!(triangles, FullDim(poly), zray, incidentpoints(poly, hidx), incidentlines(poly, hidx), incidentrays(poly, hidx), exit_point, counterclockwise, rotate; tol)
     end
     for hidx in eachindex(hyperplanes(poly))
         decompose_plane(hidx)

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,40 @@
+# 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)
+
+_default_tol(::Type{T}, ::Type{U}) where {T,U} = _default_tol(promote_type(T, U))
+
+# We should **not** define `isless(::Float64, ::MutableArithmetics.Zero)`.
+# When we get such `MethodError`, it a nice error that always indicate some underlying issue
+# so we should fix that underlying issue, not add this method that would just hide bugs.
+
 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 +46,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)