diff --git a/docs/src/systems.md b/docs/src/systems.md index fdcead3e..87e6e094 100644 --- a/docs/src/systems.md +++ b/docs/src/systems.md @@ -72,3 +72,8 @@ fix_parameters(F::AbstractSystem, p) ```@docs RandomizedSystem ``` + +## Real Systems +```@docs +is_real +``` diff --git a/src/certification.jl b/src/certification.jl index 15e7f726..a9b15789 100644 --- a/src/certification.jl +++ b/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 @@ -236,6 +236,7 @@ function squared_distance_interval( n = length(solution_candidate(cert)) d = zero(IntervalArithmetic.Interval{Float64}) for i = 1:n + @show cert yᵢ = IComplexF64(cert.I[i], a, b) d += IntervalArithmetic.sqr(real(yᵢ) - real(reference_point[i])) + @@ -320,7 +321,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 +752,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 +795,8 @@ function _certify( s, Tp[tid], Tcache[tid], - i; + i, + is_real_system; extended_certificate = extended_certificate, certify_solution_kwargs..., ) @@ -853,7 +856,8 @@ function _certify( s, p, cache, - i; + i, + is_real_system; extended_certificate = extended_certificate, certify_solution_kwargs..., ) @@ -936,7 +940,8 @@ function certify_solution( solution_candidate::AbstractVector, cert_params::Union{Nothing,CertificationParameters}, cert_cache::CertificationCache, - index::Int; + index::Int, + is_real_system::Bool; max_precision::Int = 256, refine_solution::Bool = true, extended_certificate::Bool = false, @@ -966,13 +971,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 +999,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 +1014,8 @@ function certify_solution( x̃₀, cert_params, cert_cache, - index; + index, + is_real_system; max_precision = max_precision, extended_certificate = extended_certificate, ) @@ -1013,7 +1027,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 +1051,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 +1074,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 +1089,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, @@ -1152,6 +1175,8 @@ function ε_inflation_krawczyk(x̃₀, p::Union{Nothing,CertificationParameters} sqr_mul!(x₁, M, δx) x₁ .= (x̃₀ .- Δx₀) .- x₁ certified = all2(isinterior, x₁, x₀) + @show x₀ + @show x₁ is_real = certified ? all2((x₁_i, x₀_i) -> isinterior(conj(x₁_i), x₀_i), x₁, x₀) : false else diff --git a/src/model_kit/abstract_system_homotopy.jl b/src/model_kit/abstract_system_homotopy.jl index a360844a..0e6ba22e 100644 --- a/src/model_kit/abstract_system_homotopy.jl +++ b/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 ## ##############