Merge pull request #5 from JuliaAlgebra/performance-updates

Performance updates

Author: PBrdng

README.md

@@ -38,20 +38,20 @@ julia> f_1 = x^2 + y^2 - 1;
 julia> f_2 = x^2 + y^2 - 4;
 julia> f = [f_1; f_2]
 julia> regions(f)
-RegionsResult with 3 regions:
-===============================
+egionsResult with 3 regions:
+=============================
 9 complex critical points
 5 real critical points
-╭──────────────┬─────────────────────────────────╮
-│ sign pattern │ regions                        │
-├──────────────┼─────────────────────────────────┤
-│ + +          │ number = 1                      │
-│              │ χ = 0, μ = [1, 1, 0], unbounded │
-│ - -          │ number = 1                      │
-│              │ χ = 1, μ = [1, 0, 0], bounded   │
-│ + -          │ number = 1                      │
-│              │ χ = 0, μ = [1, 1, 0], bounded   │
-╰──────────────┴─────────────────────────────────╯
+╭──────────────┬───────────────────────╮
+│ sign pattern │ regions               │
+├──────────────┼───────────────────────┤
+│ + +          │ number = 1            │
+│              │  χ = 0, μ = [1, 1, 0] │
+│ - -          │ number = 1            │
+│              │  χ = 1, μ = [1, 0, 0] │
+│ + -          │ number = 1            │
+│              │  χ = 0, μ = [1, 1, 0] │
+╰──────────────┴───────────────────────╯
 ```
 
 The output shows that $\mathcal U$ has three regions. The first region has sign pattern $++$. This means that, on this region, both $f_1$ and $f_2$ are positive. On the second region, both $f_1$ and $f_2$ are negative, so it is the contractible region in the middle. The software correctly reports that this region has Euler characteristic 1. The other two regions each have one hole and thus have Euler characteristic 0. 

docs/.documenter-siteinfo.json

@@ -1 +1 @@
-{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-09-09T16:35:10","documenter_version":"1.6.0"}}
+{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-09-12T14:21:35","documenter_version":"1.6.0"}}

docs/index.html

@@ -27,20 +27,20 @@ julia> f_1 = x^2 + y^2 - 1;
 julia> f_2 = x^2 + y^2 - 4;
 julia> f = [f_1; f_2]
 julia> C = regions(f)
-RegionsResult with 3 regions:
-===============================
+egionsResult with 3 regions:
+=============================
 9 complex critical points
 5 real critical points
-╭──────────────┬─────────────────────────────────╮
-│ sign pattern │ regions                        │
-├──────────────┼─────────────────────────────────┤
-│ + +          │ number = 1                      │
-│              │ χ = 0, μ = [1, 1, 0], unbounded │
-│ - -          │ number = 1                      │
-│              │ χ = 1, μ = [1, 0, 0], bounded   │
-│ + -          │ number = 1                      │
-│              │ χ = 0, μ = [1, 1, 0], bounded   │
-╰──────────────┴─────────────────────────────────╯
+╭──────────────┬───────────────────────╮
+│ sign pattern │ regions               │
+├──────────────┼───────────────────────┤
+│ + +          │ number = 1            │
+│              │  χ = 0, μ = [1, 1, 0] │
+│ - -          │ number = 1            │
+│              │  χ = 1, μ = [1, 0, 0] │
+│ + -          │ number = 1            │
+│              │  χ = 0, μ = [1, 1, 0] │
+╰──────────────┴───────────────────────╯
 ```
 
 The output shows that $\mathcal U$ has three regions. The first region has sign pattern $++$. This means that, on this region, both $f_1$ and $f_2$ are positive. On the second region, both $f_1$ and $f_2$ are negative, so it is the contractible region in the middle. The software correctly reports that this region has Euler characteristic 1. The other two regions each have one hole and thus have Euler characteristic 0. 
@@ -101,24 +101,12 @@ Region
 
 ```@docs
 regions
-affine_regions
 projective_regions
 ```
 
 
 ## Documentation: Helper functions
 
-Functions to call on [Region](@ref).
-```@docs
-χ
-μ
-critical_points
-is_bounded
-is_unbounded
-is_undecided
-number
-```
-
 Functions to call on [RegionsResult](@ref).
 ```@docs
 nregions
@@ -134,4 +122,14 @@ ncritical_complex
 ncritical_real
 g
 variables
+```
+
+Functions to call on [Region](@ref).
+```@docs
+χ
+μ
+critical_points
+is_bounded
+is_unbounded
+is_undecided
 ```

docs/objects.inv

@@ -47,34 +47,39 @@ function _affine_regions(
     f0::System,
     progress::Union{Nothing,RegionsProgress};
     s::Union{Nothing,Vector{T}} = nothing,
-    target_parameters::Union{Nothing, Vector{T1}} = nothing,
+    target_parameters::Union{Nothing,Vector{T1}} = nothing,
     epsilon::Float64 = 1e-6,
     reltol::Float64 = 1e-6,
     abstol::Float64 = 1e-9,
     monodromy_options = HC.MonodromyOptions(max_loops_no_progress = 10),
     start_pair_using_newton::Bool = false,
     seed = nothing,
     kwargs...,
-) where {T<:Real, T1<:Number}
+) where {T<:Real,T1<:Number}
 
     if isnothing(f0)
         f = f0
     else
-        variable_list = 
-        f = System(f0(HC.variables(f0), target_parameters), variables = HC.variables(f0))
+        variable_list =
+            f = System(
+                f0(HC.variables(f0), target_parameters),
+                variables = HC.variables(f0),
+            )
     end
 
     if !isnothing(seed)
         Random.seed!(seed)
     end
-    
-    s = set_up_s(s,f)
-    M_1, f_list = compute_critical_points(f,
-                        s,
-                        monodromy_options,
-                        progress,
-                        start_pair_using_newton;
-                        kwargs...)
+
+    s = set_up_s(s, f)
+    M_1, f_list = compute_critical_points(
+        f,
+        s,
+        monodromy_options,
+        progress,
+        start_pair_using_newton;
+        kwargs...,
+    )
     if isnothing(M_1)
         println("Couldn't solve the critical equations")
         return nothing
@@ -168,25 +173,25 @@ function compute_critical_points(
     monodromy_options::MonodromyOptions,
     progress::Union{Nothing,RegionsProgress},
     start_pair_using_newton::Bool;
-    target_parameters::Union{Nothing, Vector{T1}, Vector{Vector{T2}}} = nothing, 
+    target_parameters::Union{Nothing,Vector{T1},Vector{Vector{T2}}} = nothing,
     kwargs...,
-) where {T<:Real, T1,T2 <: Number}
+) where {T<:Real,T1,T2<:Number}
 
     variable_list = HC.variables(f)
     parameter_list = HC.parameters(f)
-    if isnothing(target_parameters) 
+    if isnothing(target_parameters)
         if length(parameter_list) > 0
             error("Number of parameters don't match.")
         end
-    elseif typeof(target_parameters) <: Vector{Vector{T}} where T<:Number
+    elseif typeof(target_parameters) <: Vector{Vector{T}} where {T<:Number}
         if length(target_parameters[1]) != length(parameter_list)
             error("Number of parameters don't match.")
         end
     elseif length(target_parameters) != length(parameter_list)
         error("Number of parameters don't match.")
     end
 
-    if typeof(target_parameters) <: Vector{T} where T<:Number
+    if typeof(target_parameters) <: Vector{T} where {T<:Number}
         target_parameters = [target_parameters]
     end
 
@@ -221,7 +226,7 @@ function compute_critical_points(
             options = monodromy_options,
             show_progress = show_progress,
         )
-        if isnothing(target_parameters) 
+        if isnothing(target_parameters)
             all_target_parameters = s
         else
             all_target_parameters = [[s; p] for p in target_parameters]
@@ -265,7 +270,8 @@ function compute_critical_points(
         if isnothing(target_parameters)
             all_target_parameters = [s; randn(ComplexF64, K - k); 0.0]
         else
-            all_target_parameters = [[s; randn(ComplexF64, K - k); 0.0; p] for p in target_parameters]
+            all_target_parameters =
+                [[s; randn(ComplexF64, K - k); 0.0; p] for p in target_parameters]
         end
     end
 
@@ -297,6 +303,7 @@ function compute_critical_points(
         start_parameters = parameters(M),
         target_parameters = all_target_parameters,
         show_progress = show_progress,
+        endgame_options = HC.EndgameOptions(; at_infinity_check = false),
         kwargs...,
     )
 
@@ -332,7 +339,8 @@ end
 function compute_parameter(J, variable_list, parameter_list, x0, q0)
     M =
         [
-            isa(entry, Expression) ? evaluate(entry, variable_list => x0, parameter_list => q0) : entry for
+            isa(entry, Expression) ?
+            evaluate(entry, variable_list => x0, parameter_list => q0) : entry for
             entry in J
         ] |> Matrix{ComplexF64}
     N = LinearAlgebra.nullspace(M)

docs/search_index.js

@@ -29,7 +29,7 @@ function membership(
     ∇logg::AS;
     kwargs...,
 ) where {AS<:AbstractSystem}
-  
+
     C = _membership(R, p, ∇logg; kwargs...)
 
     return C
@@ -59,6 +59,7 @@ function _membership(
     critical_point_index, failed_info =
         limit_critical_point(set_up_ode(∇logg), p, critical_points, reltol, abstol)
     if critical_point_index == -1
+        # println(failed_info)
         return nothing
     end
 

docs/src/index.md

@@ -71,14 +71,14 @@ critical_points(C::Region) = C.critical_points
 """
     is_unbounded(C::Region)
 
-Returns a boolean that is `true`, if `C` is undecided. 
+Returns a boolean that is `true`, if `C` was labelled unbounded. 
 """
 is_unbounded(C::Region) = C.is_bounded == 0
 
 """
     is_bounded(C::Region)
 
-Returns a boolean that is `true`, if `C` is bounded. 
+Returns a boolean that is `true`, if `C`  was labelled bounded. 
 """
 is_bounded(C::Region) = C.is_bounded == 1
 
@@ -139,9 +139,12 @@ C = regions(f)
 p = projective_regions(C)
 ```
 """
-function projective_regions(C::RegionsResult) 
+function projective_regions(C::RegionsResult)
     c = regions(C)
     p = C.projective_regions
+    if isnothing(p)
+        return p
+    end
     P = map(p) do pᵢ
         filter(ci -> in(number(ci), pᵢ), c)
     end
@@ -159,42 +162,42 @@ nregions(C::RegionsResult) = C.nregions
 """
     nbounded(C::RegionsResult)
 
-Returns the number of (weakly) bounded regions in `C`.
+Returns the number of regions in `C` that were labelled bounded (if that information was computed).
 """
 nbounded(C::RegionsResult) = count(is_bounded, regions(C))
 
 """
     nunbounded(C::RegionsResult)
 
-Returns the number of unbounded regions in `C`.
+Returns the number of regions in `C` that were labelled unbounded (if that information was computed).
 """
 nunbounded(C::RegionsResult) = count(is_unbounded, regions(C))
 
 """
     nundecided(C::RegionsResult)
 
-Returns the number of regions in `C`, where bounded or unbounded could not be decided.
+Returns the number of regions in `C`, where bounded or unbounded could not be decided (if that information was computed).
 """
 nundecided(C::RegionsResult) = count(is_undecided, regions(C))
 
 """
     bounded(C::RegionsResult)
 
-Returns the (weakly) bounded regions in `C`.
+Returns the regions in `C` that were labelled bounded.
 """
 bounded(C::RegionsResult) = filter(is_bounded, regions(C))
 
 """
     unbounded(C::RegionsResult)
 
-Returns the unbounded regions in `C`.
+Returns the regions in `C` that were labelled unbounded.
 """
 unbounded(C::RegionsResult) = filter(is_unbounded, regions(C))
 
 """
     undecided(C::RegionsResult)
 
-Returns the number of regions in `C`, where bounded or unbounded could not be decided.
+Returns the regions in `C` that were labelled undecided.
 """
 undecided(C::RegionsResult) = filter(is_undecided, regions(C))
 
@@ -243,13 +246,13 @@ Returns the order of variables.
 """
 function variables(C::RegionsResult)
     f, _, _ = g(C)
-    HC.HC.variables(f)
+    HC.variables(f)
 end
 
 ###############
 ### Show ###
 ###############
-Base.show(C::RegionsResult; crop = true) = Base.show(stdout, R, crop = crop)
+Base.show(C::RegionsResult; crop = true) = Base.show(stdout, C, crop = crop)
 function Base.show(io::IO, C::RegionsResult; crop = true)
 
     header = "RegionsResult with $(C.nregions) regions:"
@@ -293,7 +296,7 @@ function Base.show(io::IO, C::RegionsResult; crop = true)
             end
             i += 1
         end
-        
+
     end
 
     ds = displaysize()
@@ -340,7 +343,7 @@ function Base.show(io::IO, C::Region)
     if isnothing(b)
         println(io, " χ = $χ, μ = $μ")
     elseif b == 0
-       println(io, "χ = $χ, μ = $μ, unbounded")
+        println(io, "χ = $χ, μ = $μ, unbounded")
     elseif b == 1
         println(io, "χ = $χ, μ = $μ, bounded")
     elseif b == 2

src/affine_regions.jl

@@ -19,8 +19,9 @@ end
 # find the nearest point in a list to a given point
 function finding_nearest_point(v, list_of_vectors)
     v = reshape(v, length(v))
-    distances = [norm(v - x) for x in list_of_vectors]
-    if minimum(distances) < 1e-4
+    λ = LA.norm(v)
+    distances = [LA.norm(v - x) for x in list_of_vectors]
+    if minimum(distances) < 1e-4 * λ
         return argmin(distances)
     else
         return -1
@@ -192,4 +193,4 @@ function limit_critical_point_from_critical_point(
 
         return [starting_point_index, critical_point_index], failed_info
     end
-end
+end