Merge pull request #564 from JuliaHomotopyContinuation/certification_…

…reality_bugfix

Fixing the bug in issue #563

docs/src/systems.md

@@ -72,3 +72,8 @@ fix_parameters(F::AbstractSystem, p)
 ```@docs
 RandomizedSystem
 ```
+
+## Real Systems
+```@docs
+is_real
+```

src/certification.jl

@@ -107,7 +107,7 @@ is_real(C::AbstractSolutionCertificate) = C.real
     is_complex(certificate::AbstractSolutionCertificate)
 
 Returns `true` if `certificate` certifies that the certified solution interval
-contains a true complex zero of the system.
+contains a non-real complex zero of the system.
 """
 is_complex(C::AbstractSolutionCertificate) = C.complex
 
@@ -320,7 +320,6 @@ function add_certificate!(
     cert::AbstractSolutionCertificate,
 )
     d = squared_distance_interval(cert, distinct_sols.reference_point)
-
     for match in IntervalTrees.intersect(distinct_sols.distinct_tree, d)
         certᵢ = IntervalTrees.value(match)
         if Bool(Arblib.overlaps(cert.I, certᵢ.I))
@@ -752,6 +751,8 @@ function _certify(
         throw(ArgumentError("The given system expects parameters but none are given."))
     end
 
+    is_real_system = ModelKit.is_real(F)
+
     distinct_sols =
         DistinctSolutionCertificates(n; extended_certificate = extended_certificate)
     ncertified = Threads.Atomic{Int}(0)
@@ -793,7 +794,8 @@ function _certify(
                 s,
                 Tp[tid],
                 Tcache[tid],
-                i;
+                i,
+                is_real_system;
                 extended_certificate = extended_certificate,
                 certify_solution_kwargs...,
             )
@@ -853,7 +855,8 @@ function _certify(
                 s,
                 p,
                 cache,
-                i;
+                i,
+                is_real_system;
                 extended_certificate = extended_certificate,
                 certify_solution_kwargs...,
             )
@@ -930,7 +933,6 @@ end
     )
 end
 
-
 function certify_solution(
     F::AbstractSystem,
     solution_candidate::AbstractVector,
@@ -940,6 +942,30 @@ function certify_solution(
     max_precision::Int = 256,
     refine_solution::Bool = true,
     extended_certificate::Bool = false,
+)
+    certify_solution(
+        F,
+        solution_candidate,
+        cert_params,
+        cert_cache,
+        index,
+        ModelKit.is_real(F);
+        max_precision = max_precision,
+        refine_solution = refine_solution,
+        extended_certificate = extended_certificate,
+    )
+end
+
+function certify_solution(
+    F::AbstractSystem,
+    solution_candidate::AbstractVector,
+    cert_params::Union{Nothing,CertificationParameters},
+    cert_cache::CertificationCache,
+    index::Int,
+    is_real_system::Bool;
+    max_precision::Int = 256,
+    refine_solution::Bool = true,
+    extended_certificate::Bool = false,
 )
     @unpack C, arb_C, arb_x̃₀ = cert_cache
 
@@ -966,13 +992,21 @@ function certify_solution(
 
     certified, x₁, x₀, is_real = ε_inflation_krawczyk(x̃₀, cert_params, C, cert_cache)
 
+    if !is_real_system
+        is_real = false
+    end
+
     if certified
+        is_complex = any(xi -> !(0.0 in imag(xi)), x₁)
+        if is_complex
+            is_real = false
+        end
         if extended_certificate
             return ExtendedSolutionCertificate(
                 solution_candidate = solution_candidate,
                 certified = true,
                 real = is_real,
-                complex = any(xi -> !(0.0 in imag(xi)), x₁),
+                complex = is_complex,
                 index = index,
                 prec = 53,
                 I = AcbMatrix(x₀; prec = 53),
@@ -986,7 +1020,7 @@ function certify_solution(
                 solution_candidate = solution_candidate,
                 certified = true,
                 real = is_real,
-                complex = any(xi -> !(0.0 in imag(xi)), x₁),
+                complex = is_complex,
                 index = index,
                 prec = 53,
                 I = AcbMatrix(x₁; prec = 53),
@@ -1001,7 +1035,8 @@ function certify_solution(
         x̃₀,
         cert_params,
         cert_cache,
-        index;
+        index,
+        is_real_system;
         max_precision = max_precision,
         extended_certificate = extended_certificate,
     )
@@ -1013,7 +1048,8 @@ function extended_prec_certify_solution(
     x̃₀::AbstractVector,
     cert_params::Union{Nothing,CertificationParameters},
     cert_cache::CertificationCache,
-    index::Int;
+    index::Int,
+    is_real_system::Bool;
     max_precision::Int = 256,
     extended_certificate::Bool = false,
 )
@@ -1036,9 +1072,17 @@ function extended_prec_certify_solution(
         certified, arb_x₁, arb_x₀, is_real =
             arb_ε_inflation_krawczyk(arb_x̃₀, cert_params, arb_C, cert_cache; prec = prec)
 
+        if !is_real_system
+            is_real = false
+        end
+
         if certified
             I = AcbMatrix(n, 1; prec = prec)
             Arblib.set!(I, arb_x₀)
+            is_complex = any(xi -> !Arblib.contains_zero(Arblib.imagref(xi)), arb_x₁)
+            if is_complex
+                is_real = false
+            end
 
             if extended_certificate
                 I′ = AcbMatrix(n, 1; prec = prec)
@@ -1051,7 +1095,7 @@ function extended_prec_certify_solution(
                     solution_candidate = solution_candidate,
                     certified = true,
                     real = is_real,
-                    complex = any(xi -> !Arblib.contains_zero(Arblib.imagref(xi)), arb_x₁),
+                    complex = is_complex,
                     index = index,
                     prec = prec,
                     I = I,
@@ -1066,7 +1110,7 @@ function extended_prec_certify_solution(
                     solution_candidate = solution_candidate,
                     certified = true,
                     real = is_real,
-                    complex = any(xi -> !Arblib.contains_zero(Arblib.imagref(xi)), arb_x₁),
+                    complex = is_complex,
                     index = index,
                     prec = prec,
                     I = I,

src/model_kit/abstract_system_homotopy.jl

@@ -1,4 +1,4 @@
-export AbstractSystem, AbstractHomotopy, jacobian
+export AbstractSystem, AbstractHomotopy, jacobian, is_real
 
 """
     AbstractSystem
@@ -75,6 +75,20 @@ function jacobian(F::AbstractSystem, x, p = nothing)
     ModelKit.to_smallest_eltype(U)
 end
 
+"""
+    is_real(F::AbstractSystem)
+
+Checks, if `F(ℝⁿ) ⊂ ℝᵐ`, where `(m,n) = size(F)`. 
+This is an algorithm that works with probability one:
+It randomly samples `x ∈ ℝⁿ` and checks if `F(x) ∈ ℝᵐ`.
+"""
+function is_real(F::T) where {T<:AbstractSystem}
+    x = randn(Float64, size(F, 2))
+    u = F(x)
+    all(ui -> imag(ui) == 0, u)
+end
+
+
 ##############
 ## Homotopy ##
 ##############

test/certification_test.jl

@@ -75,6 +75,20 @@
         @test ndistinct_real_certified(cert) == 4
     end
 
+    @testset "Reality Check" begin
+
+        @var x
+        x0 = [1.0]
+
+        F = System([x - 1 + 1e-16 * im])
+        certF = certify(F, [x0])
+        @test is_real(certificates(certF)[1]) == false
+
+        G = System([x - 1])
+        certG = certify(G, [x0])
+        @test is_real(certificates(certG)[1]) == true
+    end
+
     @testset "Parameters: Input = Vector{Expression}" begin
         @var x y
         f = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
@@ -337,7 +351,7 @@
         f = fixed(F)
         c = HomotopyContinuation.CertificationCache(f)
         q = HomotopyContinuation.certification_parameters(p)
-        cert = HomotopyContinuation.certify_solution(f, sol, q, c, 1)
+        cert = HomotopyContinuation.certify_solution(f, sol, q, c, 1, false)
         @test is_certified(cert)
     end
 end

test/systems_test.jl

@@ -135,4 +135,12 @@ end
         G = SlicedSystem(g, L)
         test_system(G, System(G))
     end
+
+    @testset "is_real" begin
+        @var x
+        f = System([x - 1])
+        g = System([x - 1 + 1e-16 * im])
+        @test ModelKit.is_real(fixed(f)) == true
+        @test ModelKit.is_real(fixed(g)) == false
+    end
 end