Merge pull request #633 from JuliaHomotopyContinuation/automatic-triangle-inequality

Automatic triangle inequality

Author: PBrdng

src/monodromy.jl

@@ -37,7 +37,7 @@ Base.@kwdef struct MonodromyOptions{D,GA<:Union{Nothing,GroupActions}}
     permutations::Bool = false
     # unique points options
     distance::D = EuclideanNorm()
-    triangle_inequality::Bool = true
+    triangle_inequality::Union{Nothing,Bool} = nothing
     unique_points_atol::Union{Nothing,Float64} = nothing
     unique_points_rtol::Union{Nothing,Float64} = nothing
     #
@@ -458,7 +458,7 @@ function MonodromySolver(
     unique_points = UniquePoints(
         x₀,
         1;
-        metric = options.distance,
+        distance = options.distance,
         group_actions = group_actions,
         triangle_inequality = options.triangle_inequality,
     )

src/unique_points.jl

@@ -127,12 +127,14 @@ Base.iterate(p::SymmetricGroup, s) = iterate(p.permutations, s)
     UniquePoints{T, Id, M}
 
 A data structure for assessing quickly whether a point is close to an indexed point as
-determined by the given distances function `M`. The distance function has to be a *metric*.
-The indexed points are only stored by their identifiers `Id`. `triangle_inequality` should be set to `true`, if the metric satisfies the triangle inequality. Otherwise, it should be set to `false`.
+determined by the given distance function `M`. 
+The indexed points are only stored by their identifiers `Id`. 
+`triangle_inequality` should be set to `true`, if the distance function satisfies the triangle inequality. 
+Otherwise, it should be set to `false`. If `triangle_inequality` is nothing the algorithm will try to detect whether the triangle is satisfied.
 
     UniquePoints(v::AbstractVector{T}, id::Id;
-                    metric = EuclideanNorm(),
-                    triangle_inequality = true,
+                    distance = EuclideanNorm(),
+                    triangle_inequality = nothing,
                     group_actions = nothing)
 
     
@@ -167,16 +169,35 @@ end
 function UniquePoints(
     v::AbstractVector,
     id;
-    metric = EuclideanNorm(),
-    triangle_inequality = true,
+    distance = EuclideanNorm(),
+    triangle_inequality = nothing,
     group_action = nothing,
     group_actions = isnothing(group_action) ? nothing : GroupActions(group_action),
 )
     if (group_actions isa Tuple) || (group_actions isa AbstractVector)
         group_actions = GroupActions(group_actions)
     end
 
-    tree = VoronoiTree(v, id; metric = metric, triangle_inequality = triangle_inequality)
+    d = distance
+
+    if isnothing(triangle_inequality)
+        if typeof(d) <: AbstractNorm
+            triangle_inequality = true
+        else
+            n = length(v)
+            v₁ = randn(ComplexF64, n)
+            v₂ = randn(ComplexF64, n)
+            v₃ = randn(ComplexF64, n)
+            if d(v₁, v₂) ≤ d(v₁, v₃) + d(v₃, v₂) &&
+               d(v₁, 4 .* v₁) ≤ d(v₁, 2 .* v₁) + d(2 .* v₁, 4 .* v₁)
+                triangle_inequality = true
+            else
+                triangle_inequality = false
+            end
+        end
+    end
+
+    tree = VoronoiTree(v, id; distance = d, triangle_inequality = triangle_inequality)
     UniquePoints(tree, group_actions, zeros(eltype(v), length(v)))
 end
 
@@ -257,7 +278,7 @@ function add!(
     atol::Float64 = 1e-14,
     rtol::Float64 = sqrt(eps()),
 ) where {T,Id,M,GA}
-    n = UP.tree.metric(v, UP.zero_vec)
+    n = UP.tree.distance(v, UP.zero_vec)
     rad = max(atol, rtol * n)
     add!(UP, v, id, rad)
 end
@@ -266,10 +287,10 @@ end
 ## Multiplicities ##
 ####################
 """
-    multiplicities(vectors; metric = EuclideanNorm(), atol = 1e-14, rtol = 1e-8, kwargs...)
+    multiplicities(vectors; distance = EuclideanNorm(), atol = 1e-14, rtol = 1e-8, kwargs...)
 
 Returns a `Vector{Vector{Int}}` `v`. Each vector `w` in 'v' contains all indices `i`,`j`
-such that `w[i]` and `w[j]` have `distance` at most `max(atol, rtol * metric(0,w[i]))`.
+such that `w[i]` and `w[j]` have `distance` at most `max(atol, rtol * distance(0,w[i]))`.
 The remaining `kwargs` are things that can be passed to [`UniquePoints`](@ref).
 
 ```julia-repl
@@ -289,19 +310,19 @@ julia> m = multiplicities(X, group_action = permutation)
 ```
 """
 multiplicities(v; kwargs...) = multiplicities(identity, v; kwargs...)
-function multiplicities(f::F, v; metric = EuclideanNorm(), kwargs...) where {F<:Function}
+function multiplicities(f::F, v; distance = EuclideanNorm(), kwargs...) where {F<:Function}
     isempty(v) && return Vector{Vector{Int}}()
-    _multiplicities(f, v, metric; kwargs...)
+    _multiplicities(f, v, distance; kwargs...)
 end
 function _multiplicities(
     f::F,
     V,
-    metric;
+    distance;
     atol::Float64 = 1e-14,
     rtol::Float64 = 1e-8,
     kwargs...,
 ) where {F<:Function}
-    unique_points = UniquePoints(f(first(V)), 1; metric = metric, kwargs...)
+    unique_points = UniquePoints(f(first(V)), 1; distance = distance, kwargs...)
     mults = Dict{Int,Vector{Int}}()
     for (i, vᵢ) in enumerate(V)
         wᵢ = f(vᵢ)
@@ -317,14 +338,14 @@ function _multiplicities(
     collect(values(mults))
 end
 """
-    unique_points(vectors; metric = EuclideanNorm(), atol = 1e-14, rtol = 1e-8, kwargs...)
+    unique_points(vectors; distance = EuclideanNorm(), atol = 1e-14, rtol = 1e-8, kwargs...)
 
-Returns all elements in `vector` for which two elements have `distance` at most `max(atol, rtol * metric(0,w[i]))`.
+Returns all elements in `vector` for which two elements have `distance` at most `max(atol, rtol * distance(0,w[i]))`.
 Note that the output can depend on the order of elements in vectors.
 The remaining `kwargs` are things that can be passed to [`UniquePoints`](@ref).
 """
-function unique_points(V; metric = EuclideanNorm(), atol = 1e-14, rtol = 1e-8, kwargs...)
-    unique_points = UniquePoints(first(V), 1; metric = metric, kwargs...)
+function unique_points(V; distance = EuclideanNorm(), atol = 1e-14, rtol = 1e-8, kwargs...)
+    unique_points = UniquePoints(first(V), 1; distance = distance, kwargs...)
     out = Vector{eltype(V)}()
     for (i, vᵢ) in enumerate(V)
         _, new_point = add!(unique_points, vᵢ, i; atol = atol, rtol = rtol)

src/voronoi_tree.jl

@@ -29,9 +29,9 @@ Base.length(node::VTNode) = node.nentries
 Base.isempty(node::VTNode) = length(node) == 0
 capacity(node::VTNode) = length(node.children)
 
-function compute_distances!(node, x, metric::M) where {M}
+function compute_distances!(node, x, distance::M) where {M}
     for j = 1:length(node)
-        node.distances[j] = (metric(x, view(node.values, :, j)), j)
+        node.distances[j] = (distance(x, view(node.values, :, j)), j)
     end
     node.distances
 end
@@ -40,7 +40,7 @@ function search_in_radius(
     node::VTNode{T,Id},
     x,
     tol::Real,
-    metric::M,
+    distance::M,
     triangle_inequality::Bool,
 ) where {T,Id,M}
     !isempty(node) || return nothing
@@ -62,7 +62,7 @@ function search_in_radius(
     m₁ = m₂ = m₃ = (Inf, 1)
 
     # we have a distances cache per thread
-    distances = compute_distances!(node, x, metric)
+    distances = compute_distances!(node, x, distance)
     for i = 1:n
         dᵢ = first(distances[i])
         # early exit
@@ -95,7 +95,7 @@ function search_in_radius(
     # Check smallest element first
     if isassigned(node.children, last(m₁))
         retid =
-            search_in_radius(node.children[last(m₁)], x, tol, metric, triangle_inequality)
+            search_in_radius(node.children[last(m₁)], x, tol, distance, triangle_inequality)
         if !isnothing(retid)
             # we rely on the distances for insertion, so place the smallest element first
             distances[1] = m₁
@@ -112,7 +112,7 @@ function search_in_radius(
     # Case 2) If we have a triangle in equality, we know m₂[1] - m₁[1] ≤ 2tol 
     if isassigned(node.children, last(m₂))
         retid =
-            search_in_radius(node.children[last(m₂)], x, tol, metric, triangle_inequality)
+            search_in_radius(node.children[last(m₂)], x, tol, distance, triangle_inequality)
         if !isnothing(retid)
             # we rely on the distances for insertion, so place the smallest element first
             distances[1] = m₁
@@ -128,7 +128,7 @@ function search_in_radius(
     # Since we know also the third element, let's check it
     if isassigned(node.children, last(m₃))
         retid =
-            search_in_radius(node.children[last(m₃)], x, tol, metric, triangle_inequality)
+            search_in_radius(node.children[last(m₃)], x, tol, distance, triangle_inequality)
         if !isnothing(retid)
             # we rely on the distances for insertion, so place at the first place the smallest element
             distances[1] = m₁
@@ -146,8 +146,13 @@ function search_in_radius(
         dᵢ, i = distances[k]
         if dᵢ - m₁[1] < 2tol || !triangle_inequality
             if isassigned(node.children, i)
-                retid =
-                    search_in_radius(node.children[i], x, tol, metric, triangle_inequality)
+                retid = search_in_radius(
+                    node.children[i],
+                    x,
+                    tol,
+                    distance,
+                    triangle_inequality,
+                )
                 if !isnothing(retid)
                     return retid::Id
                 end
@@ -165,7 +170,7 @@ function _insert!(
     node::VTNode{T,Id},
     v,
     id::Id,
-    metric::M;
+    distance::M;
     use_distances::Bool = false,
 ) where {T,Id,M}
     # if not filled so far, just add it to the current node
@@ -179,14 +184,14 @@ function _insert!(
     if use_distances
         dᵢ, minᵢ = first(node.distances)
     else
-        compute_distances!(node, v, metric)
+        compute_distances!(node, v, distance)
         dᵢ, minᵢ = findmin(node.distances)
     end
 
     if !isassigned(node.children, minᵢ)
         node.children[minᵢ] = VTNode{T,Id}(v, id; capacity = capacity(node))
     else # a node already exists, so recurse
-        _insert!(node.children[minᵢ], v, id, metric)
+        _insert!(node.children[minᵢ], v, id, distance)
     end
 
     nothing
@@ -206,29 +211,30 @@ end
      VoronoiTree(
     v::AbstractVector{T},
     id::Id;
-    metric = EuclideanNorm(),
+    distance = EuclideanNorm(),
     capacity = 8,
     triangle_inequality = true
 )
 
 Construct a Voronoi tree data structure for vector `v` of element type `T` and with identifiers
-`Id`. Each node has the given `capacity` and distances are measured by the given `metric`. `triangle_inequality` should be set to `true`, if `metric` satisfies the triangle inequality. Otherwise, it should be set to `false`.
+`Id`. Each node has the given `capacity` and distances are measured by the given `distance`. 
+`triangle_inequality` should be set to `true`, if `distance` satisfies the triangle inequality. Otherwise, it should be set to `false`.
 """
 mutable struct VoronoiTree{T,Id,M}
     root::VTNode{T,Id}
     nentries::Int
-    metric::M
+    distance::M
     triangle_inequality::Bool
 end
 
 function VoronoiTree{T,Id}(
     d::Int;
-    metric = EuclideanNorm(),
+    distance = EuclideanNorm(),
     capacity::Int = 8,
     triangle_inequality::Bool = true,
 ) where {T,Id}
     root = VTNode(T, Id, d; capacity = capacity)
-    VoronoiTree(root, 0, metric, triangle_inequality)
+    VoronoiTree(root, 0, distance, triangle_inequality)
 end
 
 function VoronoiTree(v::AbstractVector{T}, id::Id; kwargs...) where {T,Id}
@@ -254,7 +260,7 @@ function Base.insert!(
     id::Id;
     use_distances::Bool = false,
 ) where {T,Id}
-    _insert!(tree.root, v, id, tree.metric; use_distances = use_distances)
+    _insert!(tree.root, v, id, tree.distance; use_distances = use_distances)
     tree.nentries += 1
     tree
 end
@@ -266,7 +272,7 @@ Search whether the given `tree` contains a point `p` with distances at most `tol
 Returns `nothing` if no point exists, otherwise the identifier of `p` is returned.
 """
 function search_in_radius(tree::VoronoiTree, v, tol::Real)
-    search_in_radius(tree.root, v, tol, tree.metric, tree.triangle_inequality)
+    search_in_radius(tree.root, v, tol, tree.distance, tree.triangle_inequality)
 end
 
 

test/monodromy_test.jl

@@ -69,14 +69,22 @@
             monodromy_solve(F.expressions, [x₀, rand(6)], p₀, parameters = F.parameters)
         @test length(solutions(result)) == 21
 
-        # different distance function
+        # distance function that satisfies triangle inequality
         result = monodromy_solve(F, x₀, p₀, distance = (x, y) -> 0.0)
         @test length(solutions(result)) == 1
 
+        # distance function that does not satisfy triangle inequality
+        result = monodromy_solve(F, x₀, p₀, distance = (x, y) -> norm(x - y, 2)^2)
+        @test length(solutions(result)) == 21
+
         # don't use triangle inequality
         result = monodromy_solve(F, x₀, p₀, triangle_inequality = false)
         @test length(solutions(result)) == 21
 
+        # use triangle inequality
+        result = monodromy_solve(F, x₀, p₀, triangle_inequality = true)
+        @test length(solutions(result)) == 21
+
         # Test stop heuristic with no target solutions count
         result = monodromy_solve(F, x₀, p₀)
         @test is_heuristic_stop(result)
@@ -440,7 +448,7 @@
             target_solutions_count = 305,
         )
 
-        UP = unique_points(solutions(points), metric = dist, rtol = 1e-8, atol = 1e-14)
+        UP = unique_points(solutions(points), distance = dist, rtol = 1e-8, atol = 1e-14)
 
         @test length(solutions(points)) == length(UP)
     end

test/unique_points_test.jl

@@ -338,7 +338,7 @@ end
     M = multiplicities(V)
     @test length(M) == 0
 
-    N = multiplicities(W, metric = InfNorm(), atol = 1e-5)
+    N = multiplicities(W, distance = InfNorm(), atol = 1e-5)
     sort!(N, by = first)
     @test length(N) == 10
     @test unique([length(m) for m in N]) == [2]
@@ -347,7 +347,11 @@ end
     O = multiplicities([U; U])
     @test length(O) == 20
 
-    P = multiplicities(X, metric = (x, y) -> 1 - abs(LinearAlgebra.dot(x, y)), atol = 1e-5)
+    P = multiplicities(
+        X,
+        distance = (x, y) -> 1 - abs(LinearAlgebra.dot(x, y)),
+        atol = 1e-5,
+    )
     @test length(P) == 3
 
     # Test with group action