Merge pull request #593 from JuliaHomotopyContinuation/nid_bugfixes

NID bug fixes

docs/src/witness_sets.md

@@ -42,7 +42,7 @@ NumericalIrreducibleDecomposition
 ```
 
 ```@docs
-n_components
+ncomponents
 witness_sets
 degrees
 ```

src/numerical_irreducible_decomposition.jl

@@ -4,6 +4,7 @@ export NumericalIrreducibleDecomposition,
     regeneration,
     witness_sets,
     decompose,
+    ncomponents,
     n_components,
     seed
 
@@ -104,6 +105,7 @@ function update_progress!(progress::WitnessSetsProgress, W::Union{WitnessPoints,
         showvalues = showvalues(progress),
     )
 end
+update_progress!(progress::Union{Nothing,WitnessSetsProgress}, W::Nothing) = nothing
 update_progress!(progress::Nothing) = nothing
 finish_progress!(progress::Nothing) = nothing
 function finish_progress!(progress::WitnessSetsProgress)
@@ -147,7 +149,7 @@ This is the core routine of the regeneration algorithm. It intersects a set of [
 """
 function intersect_with_hypersurface!(
     W::WitnessPoints{T1,T2,Vector{ComplexF64}},
-    X::WitnessPoints{T3,T4,Vector{ComplexF64}},
+    X::Union{WitnessPoints{T3,T4,Vector{ComplexF64}},Nothing},
     F::AS,
     H::WitnessSet{T5,T6,Vector{ComplexF64}},
     u::Variable,
@@ -174,6 +176,11 @@ function intersect_with_hypersurface!(
     deleteat!(P, m)
     update_progress!(progress, W)
 
+    if isnothing(X)
+        return nothing
+    end
+
+
     # Step 2:
     # the points in P_next are used as starting points for a homotopy.
     # where u^d-1 (u is the extra variable in u-regeneration) is deformed into g 
@@ -203,10 +210,12 @@ function intersect_with_hypersurface!(
 
     # here comes the loop for tracking
     l_start = length(start)
+
     for (i, s) in enumerate(start)
         p = s[1]
         p[end] = s[2] # the last entry of s[1] is zero. we replace it with a d-th root of unity.
 
+
         res = track(tracker, p, 1)
         if is_success(res) && is_finite(res) && is_nonsingular(res)
             new = copy(tracker.state.solution)
@@ -217,6 +226,7 @@ function intersect_with_hypersurface!(
     end
 
 
+
     nothing
 end
 
@@ -330,6 +340,7 @@ function is_contained!(
             end
             # check if x is among the points in U
             _, added = add!(U, x, 0)
+
             if added
                 return false
             else
@@ -537,11 +548,14 @@ function regeneration!(
                     # we have already intersected X ∩ Hᵢ.
                     update_progress!(progress, true)
                     for (k, W) in reverse(E) # k = codim(W) for W in Ws
-                        if k < min(i, n)
-                            X = out[k+1]
+                        if k < i
+                            if k < n
+                                X = out[k+1]
+                            else
+                                X = nothing
+                            end
                             # here is the intersection step
-                            # the next witness superset X is also passed to this function,
-                            # because we add points that do not belong to W∩Hᵢ to X.
+                            # if k < min(i,n), the next witness superset X is also passed to this function, because we add points that do not belong to W∩Hᵢ to X.
                             # at this point the equation for W is f[1:(i-1)]
                             intersect_with_hypersurface!(
                                 W,
@@ -559,7 +573,10 @@ function regeneration!(
                         end
                     end
 
+
+
                     Fᵢ = fixed(System(f[1:i], variables = vars), compile = false)
+
                     # after the first loop that takes care of intersecting with Hᵢ
                     # we now check if we have added points that are already contained in 
                     # witness sets of higher dimension.
@@ -1110,13 +1127,13 @@ witness_sets(N::NumericalIrreducibleDecomposition, dim::Int) = witness_sets(N, [
 seed(N::NumericalIrreducibleDecomposition) = N.seed
 
 """
-    n_components(N::NumericalIrreducibleDecomposition;
+    ncomponents(N::NumericalIrreducibleDecomposition;
         dims::Union{Vector{Int},Nothing} = nothing)
 
 Returns the total number of components in `N`. 
 `dims` specifies the dimensions that should be considered.
 """
-function n_components(
+function ncomponents(
     N::NumericalIrreducibleDecomposition;
     dims::Union{Vector{Int},Nothing} = nothing,
 )
@@ -1134,7 +1151,9 @@ function n_components(
 
     out
 end
-n_components(N::NumericalIrreducibleDecomposition, dim::Int) = n_components(N, dims = [dim])
+ncomponents(N::NumericalIrreducibleDecomposition, dim::Int) = ncomponents(N; dims = [dim])
+n_components(N; dims = nothing) = ncomponents(N; dims = dims)
+n_components(N, dim) = ncomponents(N; dims = [dim])
 
 """
 

test/nid_test.jl

@@ -55,6 +55,9 @@
         @test isa(N3, NumericalIrreducibleDecomposition)
 
         # number of components
+        @test ncomponents(N3) == 4
+        @test ncomponents(N3, dims = [1, 2]) == 3
+        @test ncomponents(N3, 1) == 2
         @test n_components(N3) == 4
         @test n_components(N3, dims = [1, 2]) == 3
         @test n_components(N3, 1) == 2
@@ -67,7 +70,7 @@
 
         N_Hyp = numerical_irreducible_decomposition(Hyp)
         @test degrees(N_Hyp) == Dict(3 => [5])
-        @test n_components(N_Hyp) == 1
+        @test ncomponents(N_Hyp) == 1
     end
 
     @testset "Curve of degree 6" begin
@@ -78,7 +81,7 @@
 
         N_Curve = nid(Curve)
         @test degrees(N_Curve) == Dict(1 => [6])
-        @test n_components(N_Curve) == 1
+        @test ncomponents(N_Curve) == 1
     end
 
     @testset "Overdetermined Test" begin
@@ -88,7 +91,7 @@
 
         N_TwistedCubicSphere = nid(TwistedCubicSphere)
         @test degrees(N_TwistedCubicSphere) == Dict(0 => [3])
-        @test n_components(N_TwistedCubicSphere) == 1
+        @test ncomponents(N_TwistedCubicSphere) == 1
     end
 
     @testset "Three Lines" begin
@@ -101,7 +104,7 @@
 
         N_ThreeLines = nid(ThreeLines)
         @test degrees(N_ThreeLines) == Dict(1 => [1; 1; 1])
-        @test n_components(N_ThreeLines) == 3
+        @test ncomponents(N_ThreeLines) == 3
     end
 
     @testset "Bricard6R" begin
@@ -198,7 +201,7 @@
 
         N_Bricard6R = nid(Bricard6R)
         @test degrees(N_Bricard6R) == Dict(1 => [8])
-        @test n_components(N_Bricard6R) == 1
+        @test ncomponents(N_Bricard6R) == 1
     end
 
     @testset "ACR" begin
@@ -242,6 +245,22 @@
     end
 
 
+    @testset "Union of a sphere, a line, and a point" begin
+        @var x, y, z
+
+        S = [x^2 + y^2 + z^2 - 1] #The sphere
+        L = [2 * x - z, 2 * y - z] #The line (t,t,2t)
+        P = [x + y + 2 * z - 4, y - z, x - z] #The point (1,1,1)
+
+        F = System([s * l * p for s in S for l in L for p in P])
+
+        seed = rand(UInt32)
+        NID = numerical_irreducible_decomposition(F; seed = 0x7a4845b9)
+
+        @test ncomponents(NID, 0) == 1
+        @test degrees(NID)[0] == [1]
+    end
+
     #     @testset "426" begin
     #         ### Example thanks to Julian Vill
     #         @var a, b, c, d, e
@@ -304,7 +323,7 @@
     #             endgame_options = EndgameOptions(; sing_accuracy = 1e-10),
     #         )
 
-    #         @test n_components(N) == 1
+    #         @test ncomponents(N) == 1
     #         @test degrees(N) == Dict(2 => [426])
     #     end
 end