test: add JET and fix issues

.gitignore

@@ -9,3 +9,4 @@ docs/build/
 **/**/tmp_params.txt
 **/**/tmp_sols.txt
 **/**/tmp_cert.txt
+.claude/*

Project.toml

@@ -1,7 +1,7 @@
 name = "HomotopyContinuation"
 uuid = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+version = "2.18.0"
 authors = ["Sascha Timme <sascha@timme.xyz>", "Paul Breiding <pbrdng@gmail.com>"]
-version = "2.18"
 
 [deps]
 Arblib = "fb37089c-8514-4489-9461-98f9c8763369"

benchmarks/judge_system_benchmark.jl

@@ -49,4 +49,3 @@ function judge_results(result1_name, result2_name)
     pretty_table(data; header = ["System", "evaluate", "evaluate_and_jacobian"])
     return r1, r2
 end
-

src/DoubleDouble.jl

@@ -348,8 +348,8 @@ Base.rem(a::DoubleF64, b::Union{Float64,DoubleF64}, r::RoundingMode) =
 end
 
 function Base.mod(x::DoubleF64, y::DoubleF64)
-    n = round(a / b)
-    return (a - b * n)
+    n = round(x / y)
+    return (x - y * n)
 end
 
 #
@@ -1029,7 +1029,7 @@ function Base.asin(a::DoubleF64)
     abs_a = abs(a)
 
     if abs_a > 1.0
-        throw(DomainError())
+        throw(DomainError(a, "asin requires abs(a) <= 1"))
     end
 
     if isone(abs_a)
@@ -1043,7 +1043,7 @@ function Base.acos(a::DoubleF64)
     abs_a = abs(a)
 
     if abs_a > 1.0
-        throw(DomainError())
+        throw(DomainError(a, "acos requires abs(a) <= 1"))
     end
 
     if isone(abs_a)

src/certification.jl

@@ -245,7 +245,7 @@ function squared_distance_interval(
     n = length(solution_candidate(cert))
     d = zero(IntervalArithmetic.Interval{Float64})
     for i = 1:n
-        yᵢ = IComplexF64(cert.I[i], a, b)
+        yᵢ = IComplexF64(something(cert.I)[i], a, b)
         d +=
             IntervalArithmetic.sqr(real(yᵢ) - real(reference_point[i])) +
             IntervalArithmetic.sqr(imag(yᵢ) - imag(reference_point[i]))
@@ -331,7 +331,7 @@ function add_certificate!(
     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))
+        if Bool(Arblib.overlaps(something(cert.I), something(certᵢ.I)))
             return (false, certᵢ)
             break
         end
@@ -362,7 +362,7 @@ function is_solution_candidate_guaranteed_duplicate(
         end
 
         # Check if s is contained in certᵢ.I
-        if Bool(Arblib.contains(certᵢ.I, distinct_sols.acb_solution_candidate))
+        if Bool(Arblib.contains(something(certᵢ.I), distinct_sols.acb_solution_candidate))
             return true
             break
         end
@@ -608,7 +608,10 @@ function CertificationCache(F::AbstractSystem)
         jac_interpreter_F64 = jac_interpreter_F64,
         eval_interpreter_acb = eval_interpreter_acb,
         jac_interpreter_acb = jac_interpreter_acb,
-        newton_cache = NewtonCache(F; optimize_data_structure = false),
+        newton_cache = NewtonCache(
+            F;
+            optimize_data_structure = false,
+        )::NewtonCache{MatrixWorkspace{Matrix{ComplexF64}}},
         C = zeros(ComplexF64, m, m),
         r₀ = zeros(IComplexF64, m),
         Δx₀ = zeros(IComplexF64, m),
@@ -1173,9 +1176,13 @@ function extended_prec_certify_solution(
 end
 
 
-function Base.Vector{ComplexF64}(A::Arblib.AcbMatrixLike)
+function Base.Vector{ComplexF64}(A::Arblib.AcbMatrix)
     @assert size(A, 2) == 1
-    [ComplexF64(Arblib.ref(A, i, 1).acb_ptr) for i = 1:size(A, 1)]
+    [ComplexF64(Arblib.ref(A, i, 1)) for i = 1:size(A, 1)]
+end
+function Base.Vector{ComplexF64}(A::Arblib.AcbRefMatrix)
+    @assert size(A, 2) == 1
+    [ComplexF64(Arblib.ref(A, i, 1)) for i = 1:size(A, 1)]
 end
 
 function ε_inflation_krawczyk(x̃₀, p::Union{Nothing,CertificationParameters}, C, cert_cache)
@@ -1509,7 +1516,7 @@ Base.@kwdef mutable struct DistinctCertifiedSolutions{
     C<:AbstractSolutionCertificate,
 }
     system::Vector{S}
-    parameters::Union{Vector{Nothing},Vector{CertificationParameters}}
+    parameters::Vector{Union{Nothing,CertificationParameters}}
     cache::Vector{CertificationCache}
     access_lock::ReentrantLock
     distinct_solution_certificates::DistinctSolutionCertificates{C}
@@ -1543,7 +1550,10 @@ function DistinctCertifiedSolutions(
     nthreads = thread_safe ? Threads.nthreads() : 1
     return DistinctCertifiedSolutions(
         [F; [deepcopy(F) for _ = 2:nthreads]],
-        [parameters; [deepcopy(parameters) for _ = 2:nthreads]],
+        Union{Nothing,CertificationParameters}[
+            parameters
+            [deepcopy(parameters) for _ = 2:nthreads]
+        ],
         [cache; [deepcopy(cache) for _ = 2:nthreads]],
         access_lock,
         distinct_solution_certificates,

src/endgame_tracker.jl

@@ -595,7 +595,7 @@ function singular_endgame_step!(endgame_tracker::EndgameTracker, debug::Bool = f
     end
 
     @label prediction
-    m = state.winding_number
+    m = something(state.winding_number)
     κ = state.sample_conds[3]
     zero_cond = 1 / (m + 1)
     for i = 1:length(state.prediction)

src/homotopies/mixed_homotopy.jl

@@ -10,7 +10,8 @@ struct MixedHomotopy{ID} <: AbstractHomotopy
     compiled::CompiledHomotopy{ID}
     interpreted::InterpretedHomotopy
 end
-MixedHomotopy(H::Homotopy) = MixedHomotopy(CompiledHomotopy(H), InterpretedHomotopy(H))
+MixedHomotopy(H::Homotopy; kwargs...) =
+    MixedHomotopy(CompiledHomotopy(H), InterpretedHomotopy(H))
 
 Base.size(H::MixedHomotopy) = size(H.compiled)
 ModelKit.variables(H::MixedHomotopy) = variables(H.interpreted)

src/homotopies/subspace_homotopies.jl

@@ -43,6 +43,8 @@ function GrassmannianGeodesic(start, target; embedded_projective::Bool = false)
         else
             Q_cos = target.Y * U
         end
+    else
+        error("start must be an ExtrinsicDescription or IntrinsicDescription")
     end
 
 
@@ -104,8 +106,8 @@ Base.@kwdef mutable struct ExtrinsicSubspaceHomotopy{S<:AbstractSystem} <: Abstr
     b0::Vector{ComplexF64}
 
     # For the offset part (linear interpolation)
-    a_minus_b::Union{Nothing,Vector{ComplexF64}}
-    offset::Union{Nothing,Vector{ComplexF64}}
+    a_minus_b::Vector{ComplexF64}
+    offset::Vector{ComplexF64}
 
     # caches for t
     t_cache::Base.RefValue{ComplexF64}
@@ -177,8 +179,8 @@ function ExtrinsicSubspaceHomotopy(
     # get correct coordinates for A and b in the Stiefel homotopy
     # extrinsic(start).A is replaced by transpose(path.γ1)
     # extrinsic(target).A is replaced by transpose(path.Q_cos)    
-    a0 = path.B_start * extrinsic(start).b
-    b0 = path.B_target * extrinsic(target).b
+    a0 = something(path.B_start) * extrinsic(start).b
+    b0 = something(path.B_target) * extrinsic(target).b
 
     # Prepare offset data for linear interpolation
     a_minus_b = a0 - b0
@@ -237,10 +239,10 @@ Base.@kwdef mutable struct IntrinsicSubspaceHomotopy{S<:AbstractSystem} <: Abstr
     path::GrassmannianGeodesic
 
     # For the offset part (linear interpolation)
-    a_minus_b::Union{Nothing,Vector{ComplexF64}}
-    offset::Union{Nothing,Vector{ComplexF64}}
+    a_minus_b::Vector{ComplexF64}
+    offset::Vector{ComplexF64}
 
-    # cache for t 
+    # cache for t
     t_cache::Base.RefValue{ComplexF64}
     offset_t_cache::Base.RefValue{ComplexF64}
 
@@ -493,8 +495,8 @@ function set_subspaces!(H::SubspaceHomotopy, start::LinearSubspace, target::Line
 
     # Update the offsets
     if isa(H, ExtrinsicSubspaceHomotopy)
-        LA.mul!(H.a0, H.path.B_start, extrinsic(start).b)
-        LA.mul!(H.b0, H.path.B_target, extrinsic(target).b)
+        LA.mul!(H.a0, something(H.path.B_start), extrinsic(start).b)
+        LA.mul!(H.b0, something(H.path.B_target), extrinsic(target).b)
         H.a_minus_b .= H.a0 .- H.b0
         H.offset .= H.b0
     elseif isa(H, IntrinsicSubspaceHomotopy)

src/interval_arithmetic.jl

@@ -47,7 +47,7 @@ is_valid_interval(a::Real, b::Real) = isfinite(a) && isfinite(b) && a ≤ b
     )
 end
 
-Base.hash(x::Interval, h::UInt) = hash(x.hi, hash(x.lo, u))
+Base.hash(x::Interval, h::UInt) = hash(x.hi, hash(x.lo, h))
 Base.:(==)(a::Interval, b::Interval) = a.lo == b.lo && a.hi == b.hi
 Base.eltype(x::Interval{T}) where {T} = T
 
@@ -246,7 +246,7 @@ Base.literal_pow(::typeof(^), a::Interval, ::Val{2}) = sqr(a)
 
 function ^(x::Interval, n::Integer)  # fast integer power
     if n < 0
-        return inv(pow(x, -n))
+        return inv(x^(-n))
     end
     isempty(x) && return x
     if iseven(n) && 0 ∈ x

src/linear_algebra.jl

@@ -27,15 +27,20 @@ function MatrixWorkspace(Â::AbstractMatrix; optimize_data_structure = true)
     m, n = size(Â)
     m ≥ n || throw(ArgumentError("Expected system with more rows than columns."))
 
-    A = Matrix{ComplexF64}(Â)
-    d = ones(m)
-    factorized = Ref(false)
-    qr = LA.qrfactUnblocked!(copy(A))
+    A_mat = Matrix{ComplexF64}(Â)
     # experiments show that for m > 25 the data layout as a
     # struct array is beneficial
     if m > 25 && optimize_data_structure
-        A = StructArrays.StructArray(A)
+        return _make_matrix_workspace(StructArrays.StructArray(A_mat), A_mat, m, n)
+    else
+        return _make_matrix_workspace(A_mat, A_mat, m, n)
     end
+end
+
+function _make_matrix_workspace(A, A_mat::Matrix{ComplexF64}, m, n)
+    d = ones(m)
+    factorized = Ref(false)
+    qr = LA.qrfactUnblocked!(copy(A_mat))
     row_scaling = ones(m)
     scaled = Ref(false)
 
@@ -339,7 +344,7 @@ function ldiv_adj_upper!(
         for i = 1:j-1
             z -= conj(A[i, j]) * x[i]
         end
-        iszero(A[j, j]) && singular_exception && throw(SingularException(j))
+        iszero(A[j, j]) && singular_exception && throw(LA.SingularException(j))
         x[j] = @fastmath conj(A[j, j]) \ z
     end
     x
@@ -752,7 +757,7 @@ function LA.cond(
         if isa(d_l, Nothing) && isa(d_r, Nothing)
             inv(abs(WS.A[1, 1]))
         elseif isa(d_l, Nothing)
-            inv(abs(WS.A[1, 1]) * d_r[1])
+            inv(abs(WS.A[1, 1]) * something(d_r)[1])
         elseif isa(d_r, Nothing)
             inv(d_l[1] * abs(WS.A[1, 1]))
         else

src/model_kit/compiled_system_homotopy.jl

@@ -19,7 +19,7 @@ struct CompiledSystem{HI} <: AbstractSystem
     system::System
 end
 
-function CompiledSystem(F::System)
+function CompiledSystem(F::System; kwargs...)
     n = length(F)
     nvars = nvariables(F)
     nparams = nparameters(F)
@@ -104,7 +104,7 @@ struct CompiledHomotopy{HI} <: AbstractHomotopy
     homotopy::Homotopy
 end
 
-function CompiledHomotopy(H::Homotopy)
+function CompiledHomotopy(H::Homotopy; kwargs...)
     n = length(H)
     nvars = nvariables(H)
     nparams = nparameters(H)

src/model_kit/instruction_interpreter.jl

@@ -275,7 +275,7 @@ for has_parameters in [true, false],
     has_second_output in [true, false]
 
     @eval Base.@propagate_inbounds function execute!(
-        u::Union{Nothing,AbstractArray},
+        $(has_second_output ? :(u::Union{Nothing,AbstractArray}) : :(u::AbstractArray)),
         $((has_second_output ? (:(U::AbstractArray),) : ())...),
         I::Interpreter,
         x::AbstractArray,
@@ -288,9 +288,11 @@ for has_parameters in [true, false],
         checkbounds(x, 1:length(vars_range))
         isnothing(parameters) || checkbounds(parameters, 1:length(params_range))
 
-        @inbounds for (i, k) in enumerate(params_range)
-            I.tape[k] = parameters[i]
-        end
+        $(has_parameters ? quote
+            @inbounds for (i, k) in enumerate(params_range)
+                I.tape[k] = parameters[i]
+            end
+        end : :())
         $(
             has_continuation_parameter ?
             quote
@@ -444,7 +446,7 @@ end
 for has_parameters in [true, false], has_continuation_parameter in [true, false]
 
     @eval Base.@propagate_inbounds function execute_taylor!(
-        u::Union{Nothing,AbstractArray},
+        u::AbstractArray,
         V::Val{K},
         I::Interpreter,
         x::AbstractArray,
@@ -459,9 +461,11 @@ for has_parameters in [true, false], has_continuation_parameter in [true, false]
         checkbounds(x, 1:length(vars_range))
         isnothing(parameters) || checkbounds(parameters, 1:length(params_range))
 
-        @inbounds for (i, k) in enumerate(params_range)
-            I.tape[k] = parameters[i]
-        end
+        $(has_parameters ? quote
+            @inbounds for (i, k) in enumerate(params_range)
+                I.tape[k] = parameters[i]
+            end
+        end : :())
         $(
             has_continuation_parameter ?
             quote

src/model_kit/intermediate_representation.jl

@@ -381,7 +381,7 @@ function split_into_num_denom!(ir, ex, cse, pse)
             xref = expr_to_ir_statements!(ir, x, cse, pse)
             k = to_number(k_ex)
             if k isa Basic
-                throw(ExprError("Cannot handle non-constant exponents"))
+                error("Cannot handle non-constant exponents")
             end
             if k < 0
                 push!(denoms, pow!(ir, xref, -k))

src/model_kit/interpreted_homotopy.jl

@@ -23,7 +23,7 @@ mutable struct InterpretedHomotopy <: AbstractHomotopy
     jac_acb::Union{Nothing,Interpreter{AcbRefVector}}
 end
 
-function InterpretedHomotopy(H::Homotopy)
+function InterpretedHomotopy(H::Homotopy; kwargs...)
     eval_ComplexF64 = interpreter(ComplexF64, H)
     eval_ComplexDF64 = interpreter(ComplexDF64, eval_ComplexF64)
     eval_acb = nothing
@@ -46,6 +46,7 @@ Base.size(H::InterpretedHomotopy) = size(H.homotopy)
 variables(H::InterpretedHomotopy) = variables(H.homotopy)
 parameters(H::InterpretedHomotopy) = parameters(H.homotopy)
 variable_groups(H::InterpretedHomotopy) = variable_groups(H.homotopy)
+Homotopy(H::InterpretedHomotopy) = H.homotopy
 Base.:(==)(H::InterpretedHomotopy, G::InterpretedHomotopy) = H.homotopy == G.homotopy
 
 function Base.show(io::IO, H::InterpretedHomotopy)
@@ -112,8 +113,9 @@ function evaluate!(
     if isnothing(H.eval_acb)
         H.eval_acb = interpreter(AcbRefVector, H.eval_ComplexF64)
     end
-    setprecision!(H.eval_acb, prec)
-    execute!(u, H.eval_acb, x, t, p)
+    I = H.eval_acb::Interpreter{AcbRefVector}
+    setprecision!(I, prec)
+    execute!(u, I, x, t, p)
 end
 
 function evaluate_and_jacobian!(
@@ -128,9 +130,10 @@ function evaluate_and_jacobian!(
     if isnothing(H.jac_acb)
         H.jac_acb = interpreter(AcbRefVector, H.jac_ComplexF64)
     end
-    setprecision!(H.jac_acb, prec)
+    I = H.jac_acb::Interpreter{AcbRefVector}
+    setprecision!(I, prec)
 
-    execute!(u, U, H.jac_acb, x, t, p)
+    execute!(u, U, I, x, t, p)
     nothing
 end