From 879fbd3b9f7325485e57ccf9d61e6d406251d85b Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 26 Dec 2024 16:13:07 +0530
Subject: [PATCH 01/24] feat: add `NonlinearSolveHomotopyContinuation.jl`

---
 .../LICENSE                                   | 21 ++++++++++++
 .../Project.toml                              | 34 +++++++++++++++++++
 .../src/NonlinearSolveHomotopyContinuation.jl | 27 +++++++++++++++
 .../test/runtests.jl                          | 10 ++++++
 4 files changed, 92 insertions(+)
 create mode 100644 lib/NonlinearSolveHomotopyContinuation/LICENSE
 create mode 100644 lib/NonlinearSolveHomotopyContinuation/Project.toml
 create mode 100644 lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
 create mode 100644 lib/NonlinearSolveHomotopyContinuation/test/runtests.jl

diff --git a/lib/NonlinearSolveHomotopyContinuation/LICENSE b/lib/NonlinearSolveHomotopyContinuation/LICENSE
new file mode 100644
index 000000000..d0d72edb9
--- /dev/null
+++ b/lib/NonlinearSolveHomotopyContinuation/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2024 Aayush Sabharwal <aayush.sabharwal@gmail.com> and contributors
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/lib/NonlinearSolveHomotopyContinuation/Project.toml b/lib/NonlinearSolveHomotopyContinuation/Project.toml
new file mode 100644
index 000000000..6203ce640
--- /dev/null
+++ b/lib/NonlinearSolveHomotopyContinuation/Project.toml
@@ -0,0 +1,34 @@
+name = "NonlinearSolveHomotopyContinuation"
+uuid = "2ac3b008-d579-4536-8c91-a1a5998c2f8b"
+authors = ["Aayush Sabharwal <aayush.sabharwal@gmail.com> and contributors"]
+version = "0.1.0"
+
+[deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
+ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
+
+[compat]
+ADTypes = "1.11.0"
+CommonSolve = "0.2.4"
+ConcreteStructs = "0.2.3"
+DifferentiationInterface = "0.6.27"
+HomotopyContinuation = "2.12.0"
+LinearAlgebra = "1.11.0"
+NonlinearSolveBase = "1.3.3"
+SciMLBase = "2.68.1"
+SymbolicIndexingInterface = "0.3.36"
+julia = "1.10"
+
+[extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Aqua", "Test"]
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
new file mode 100644
index 000000000..0f967d2b0
--- /dev/null
+++ b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
@@ -0,0 +1,27 @@
+module NonlinearSolveHomotopyContinuation
+
+using SciMLBase
+using NonlinearSolveBase
+using SymbolicIndexingInterface
+using LinearAlgebra
+using ADTypes
+import CommonSolve
+using ConcreteStructs: @concrete
+
+export HomotopyContinuationJL, HomotopyNonlinearFunction
+
+@concrete struct HomotopyContinuationJL{AllRoots} <: NonlinearSolveBase.AbstractNonlinearSolveAlgorithm
+    autodiff
+    kwargs
+end
+
+function HomotopyContinuationJL{AllRoots}(; autodiff = true, kwargs...) where {AllRoots}
+    if autodiff isa Bool
+        autodiff = autodiff ? AutoForwardDiff() : AutoFiniteDiff()
+    end
+    HomotopyContinuationJL{AllRoots}(autodiff, kwargs)
+end
+
+HomotopyContinuationJL(; kwargs...) = HomotopyContinuationJL{false}(; kwargs...)
+
+end
diff --git a/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl b/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl
new file mode 100644
index 000000000..aedd9485a
--- /dev/null
+++ b/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl
@@ -0,0 +1,10 @@
+using NonlinearSolveHomotopyContinuation
+using Test
+using Aqua
+
+@testset "NonlinearSolveHomotopyContinuation.jl" begin
+    @testset "Code quality (Aqua.jl)" begin
+        Aqua.test_all(NonlinearSolveHomotopyContinuation)
+    end
+    # Write your tests here.
+end

From 6addac36a55e41041056da21020f1448ce98d291 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 26 Dec 2024 16:13:23 +0530
Subject: [PATCH 02/24] feat: add `solve` implementations

---
 .../src/NonlinearSolveHomotopyContinuation.jl |   7 +
 .../src/interface_types.jl                    |  87 +++++++++++++
 .../src/solve.jl                              | 123 ++++++++++++++++++
 3 files changed, 217 insertions(+)
 create mode 100644 lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
 create mode 100644 lib/NonlinearSolveHomotopyContinuation/src/solve.jl

diff --git a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
index 0f967d2b0..46b537e35 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
@@ -1,11 +1,15 @@
 module NonlinearSolveHomotopyContinuation
 
+using SciMLBase: AbstractNonlinearProblem
 using SciMLBase
 using NonlinearSolveBase
 using SymbolicIndexingInterface
 using LinearAlgebra
 using ADTypes
 import CommonSolve
+import HomotopyContinuation as HC
+import DifferentiationInterface as DI
+
 using ConcreteStructs: @concrete
 
 export HomotopyContinuationJL, HomotopyNonlinearFunction
@@ -24,4 +28,7 @@ end
 
 HomotopyContinuationJL(; kwargs...) = HomotopyContinuationJL{false}(; kwargs...)
 
+include("interface_types.jl")
+include("solve.jl")
+
 end
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
new file mode 100644
index 000000000..4b47d2212
--- /dev/null
+++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
@@ -0,0 +1,87 @@
+abstract type HomotopySystemVariant end
+
+struct Inplace <: HomotopySystemVariant end
+struct OutOfPlace <: HomotopySystemVariant end
+struct Scalar <: HomotopySystemVariant end
+
+@concrete struct HomotopySystemWrapper{variant <: HomotopySystemVariant} <: HC.AbstractSystem
+    prob
+    autodiff
+    prep
+    vars
+end
+
+Base.size(sys::HomotopySystemWrapper) = (length(sys.prob.u0), length(sys.prob.u0))
+HC.ModelKit.variables(sys::HomotopySystemWrapper) = sys.vars
+
+function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{Inplace}, x, p = nothing)
+    sys.prob.f.f(u, x, parameter_values(sys.prob))
+end
+
+function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing)
+    values = sys.prob.f.f(x, parameter_values(sys.prob))
+    copyto!(u, values)
+end
+
+function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{Scalar}, x, p = nothing)
+    u[1] = sys.prob.f.f(only(x), parameter_values(sys.prob))
+end
+
+function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Inplace}, x, p = nothing)
+    f = sys.prob.f
+    p = parameter_values(sys.prob)
+    if SciMLBase.has_jac(f)
+        f.f(u, x, p)
+        f.jac(U, x, p)
+        return
+    end
+
+    DI.value_and_jacobian!(f.f, u, U, sys.prep, sys.autodiff, x, DI.Constant(p))
+end
+
+function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing)
+    f = sys.prob.f
+    p = parameter_values(sys.prob)
+    if SciMLBase.has_jac(f)
+        u_tmp = f.f(x, p)
+        copyto!(u, u_tmp)
+        j_tmp = f.jac(U, x, p)
+        copyto!(U, j_tmp)
+        return
+    end
+    u_tmp, _ = DI.value_and_jacobian(f.f, U, sys.prep, sys.autodiff, x, DI.Constant(p))
+    copyto!(u, u_tmp)
+end
+
+function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Scalar}, x, p = nothing)
+    f = sys.prob.f
+    p = parameter_values(sys.prob)
+    if SciMLBase.has_jac(f)
+        HC.ModelKit.evaluate!(u, sys, x, p)
+        U[1] = f.jac(only(x), p)
+    else
+        u[1], U[1] = DI.value_and_derivative(f.f, sys.prep, sys.autodiff, only(x), DI.Constant(p))
+    end
+    return nothing
+end
+
+@concrete struct GuessHomotopy <: HC.AbstractHomotopy
+    sys <: HC.AbstractSystem
+    fu0
+end
+
+Base.size(h::GuessHomotopy) = size(h.sys)
+
+function HC.ModelKit.evaluate!(u, h::GuessHomotopy, x, t, p = nothing)
+    HC.ModelKit.evaluate!(u, h.sys, x, p)
+    @inbounds for i in eachindex(u)
+        u[i] += (t + 1) * h.fu0[i]
+    end
+end
+
+function HC.ModelKit.evaluate_and_jacobian!(u, U, h::GuessHomotopy, x, t, p = nothing)
+    HC.ModelKit.evaluate_and_jacobian!(u, U, h.sys, x, p)
+    @inbounds for i in eachindex(u)
+        u[i] += (t + 1) * h.fu0[i]
+    end
+end
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
new file mode 100644
index 000000000..2d30874fe
--- /dev/null
+++ b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
@@ -0,0 +1,123 @@
+function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::HomotopyContinuationJL)
+    # cast to a `HomotopyNonlinearFunction`
+    f = if prob.f isa HomotopyNonlinearFunction
+        prob.f
+    else
+        HomotopyNonlinearFunction(prob.f)
+    end
+
+    u0 = state_values(prob)
+    p = parameter_values(prob)
+    isscalar = u0 isa Number
+    iip = SciMLBase.isinplace(prob)
+
+    # jacobian handling
+    if SciMLBase.has_jac(f)
+        # use if present
+        prep = nothing
+    elseif iip
+        # prepare a DI jacobian if not
+        prep = DI.prepare_jacobian(f, copy(state_values(prob)), alg.autodiff, u0, DI.Constant(p))
+    elseif isscalar
+        prep = DI.prepare_derivative(f, alg.autodiff, u0, DI.Constant(p))
+    else
+        prep = DI.prepare_jacobian(f, alg.autodiff, u0, DI.Constant(p))
+    end
+
+    # variables for HC to use
+    vars = if isscalar
+        HC.variables(:x)
+    else
+        HC.variables(:x, axes(u0)...)
+    end
+
+    # HC-compatible system
+    variant = iip ? Inplace : isscalar ? Scalar : OutOfPlace
+    hcsys = HomotopySystemWrapper{variant}(prob, alg.autodiff, prep, vars)
+
+    return f, hcsys
+end
+
+function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{true}; denominator_abstol = 1e-7, kwargs...)
+    f, hcsys = homotopy_continuation_preprocessing(prob, alg)
+
+    u0 = state_values(prob)
+    p = parameter_values(prob)
+    isscalar = u0 isa Number
+
+    orig_sol = HC.solve(hcsys; alg.kwargs..., kwargs...)
+    realsols = HC.results(orig_sol; only_real = true)
+    # no real solutions
+    if isempty(realsols)
+        retcode = SciMLBase.ReturnCode.ConvergenceFailure
+        resid = NonlinearSolveBase.Utils.evaluate_f(prob, u0)
+        nlsol = SciMLBase.build_solution(prob, alg, u0, resid; retcode, original = orig_sol)
+        return SciMLBase.EnsembleSolution([nlsol], 0.0, false, nothing)
+    end
+    T = eltype(u0)
+    validsols = isscalar ? T[] : Vector{T}[]
+    for result in realsols
+        # ignore ones which make the denominator zero
+        real_u = real.(result.solution)
+        if isscalar
+            real_u = only(real_u)
+        end
+        if any(<=(denominator_abstol) ∘ abs, f.denominator(real_u, p))
+            continue
+        end
+        # unpack solutions and make them real
+        u = isscalar ? only(result.solution) : result.solution
+        append!(validsols, f.unpolynomialize(real.(u), p))
+    end
+
+    # if there are no valid solutions
+    if isempty(validsols)
+        retcode = SciMLBase.ReturnCode.Infeasible
+        resid = NonlinearSolveBase.Utils.evaluate_f(prob, u0)
+        nlsol = SciMLBase.build_solution(prob, alg, u0, resid; retcode, original = orig_sol)
+        return SciMLBase.EnsembleSolution([nlsol], 0.0, false, nothing)
+    end
+
+    retcode = SciMLBase.ReturnCode.Success
+    nlsols = map(validsols) do u
+        resid = NonlinearSolveBase.Utils.evaluate_f(prob, u)
+        return SciMLBase.build_solution(prob, alg, u, resid; retcode, original = orig_sol)
+    end
+    return SciMLBase.EnsembleSolution(nlsols, 0.0, true, nothing)
+end
+
+function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{false}; denominator_abstol = 1e-7, kwargs...)
+    f, hcsys = homotopy_continuation_preprocessing(prob, alg)
+
+    u0 = state_values(prob)
+    p = parameter_values(prob)
+
+    u0_p = f.polynomialize(u0, p)
+    fu0 = NonlinearSolveBase.Utils.evaluate_f(prob, u0_p)
+
+    homotopy = GuessHomotopy(hcsys, fu0)
+    orig_sol = HC.solve(homotopy, [fu0]; alg.kwargs..., kwargs...)
+    realsols = HC.results(orig_sol; only_real = true)
+
+    # no real solutions or infeasible solution
+    if isempty(realsols) || any(<=(denominator_abstol), f.denominator(real.(only(realsols).solution), p))
+        retcode = if isempty(realsols)
+            SciMLBase.ReturnCode.ConvergenceFailure
+        else
+            SciMLBase.ReturnCode.Infeasible
+        end
+        resid = NonlinearSolveBase.Utils.evaluate_f(prob, u0)
+        return SciMLBase.build_solution(prob, alg, u0, resid; retcode, original = orig_sol)
+    end
+    realsol = only(realsols)
+    T = eltype(u0)
+    validsols = f.unpolynomialize(realsol.solution, p)
+    sol, idx = findmin(validsols) do sol
+        norm(sol - u0)
+    end
+    resid = NonlinearSolveBase.Utils.evaluate_f(prob, u0)
+
+    retcode = SciMLBase.ReturnCode.Success
+    return SciMLBase.build_solution(prob, alg, u0, resid; retcode, original = orig_sol)
+end
+

From 12ea8840c3a09f6884072820e5e187ad4901b5b5 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 26 Dec 2024 16:13:34 +0530
Subject: [PATCH 03/24] test: add tests for all roots solve

---
 .../test/allroots.jl                          | 122 ++++++++++++++++++
 .../test/runtests.jl                          |   3 +
 2 files changed, 125 insertions(+)
 create mode 100644 lib/NonlinearSolveHomotopyContinuation/test/allroots.jl

diff --git a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
new file mode 100644
index 000000000..581ea0b22
--- /dev/null
+++ b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
@@ -0,0 +1,122 @@
+using NonlinearSolve
+using NonlinearSolveHomotopyContinuation
+using SciMLBase: NonlinearSolution
+
+alg = HomotopyContinuationJL{true}(; threading = false)
+
+@testset "scalar u" begin
+    @testset "`NonlinearProblem`" begin
+        prob = NonlinearProblem(1.0, 2.0) do u, p
+            return u * u - 2p * u + p * p
+        end
+        sol = solve(prob, alg)
+
+        @test sol isa EnsembleSolution
+        @test length(sol) == 1
+        @test sol.u[1] isa NonlinearSolution
+        @test sol.u[1][1] ≈ 2.0
+
+        @testset "no real solutions" begin
+            prob = NonlinearProblem(1.0, 0.5) do u, p
+                return u * u - 2p * u  + p
+            end
+            sol = solve(prob, alg)
+            @test length(sol) == 1
+            @test sol.u[1].retcode == SciMLBase.ReturnCode.ConvergenceFailure
+            @test !sol.converged
+        end
+    end
+
+    @testset "`HomotopyContinuationFunction`" begin
+        denominator = function (u, p)
+            return [u - 0.7]
+        end
+        polynomialize = function (u, p)
+            return sin(u)
+        end
+        unpolynomialize = function (u, p)
+            return [asin(u)]
+        end
+        fn = HomotopyNonlinearFunction(; denominator, polynomialize, unpolynomialize) do u, p
+            return (u - p[1]) * (u - p[2])
+        end
+        prob = NonlinearProblem(fn, 0.0, [0.5, 0.7])
+
+        sol = solve(prob, alg)
+        @test length(sol) == 1
+        @test sin(sol.u[1][1]) ≈ 0.5
+
+        @testset "no valid solutions" begin
+            prob2 = remake(prob; p = [0.7, 0.7])
+            sol2 = solve(prob2, alg)
+            @test !sol2.converged
+            @test length(sol2) == 1
+            @test sol2.u[1].retcode == SciMLBase.ReturnCode.Infeasible
+        end
+
+        @testset "multiple solutions" begin
+            prob3 = remake(prob; p = [0.5, 0.6])
+            sol3 = solve(prob3, alg)
+            @test length(sol3) == 2
+            @test sort(sin.([sol3.u[1][1], sol3.u[2][1]])) ≈ [0.5, 0.6]
+        end
+    end
+end
+
+function f!(du, u, p)
+    du[1] = u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1
+    du[2] = u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]
+end
+
+function f(u, p)
+    [u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1, u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]]
+end
+
+@testset "vector u - $iip" for (rhs, iip) in [(f, "oop"), (f!, "iip")]
+    sol = nothing
+    @testset "`NonlinearProblem`" begin
+        prob = NonlinearProblem(rhs, [1.0, 2.0], [2.0, 3.0])
+        sol = solve(prob, alg)
+        @test sol isa EnsembleSolution
+        @test sol.converged
+        for nlsol in sol.u
+            @test SciMLBase.successful_retcode(nlsol)
+            @test f(nlsol.u, prob.p) ≈ [0.0, 0.0] atol = 1e-10
+        end
+
+        @testset "no real solutions" begin
+            _prob = remake(prob; p = zeros(2))
+            _sol = solve(_prob, alg)
+            @test !_sol.converged
+            @test length(_sol) == 1
+            @test !SciMLBase.successful_retcode(_sol.u[1])
+        end
+    end
+
+    @testset "`HomotopyNonlinearFunction`" begin
+        denominator = function (u, p)
+            return [u[1] - p[3], u[2] - p[4]]
+        end
+        unpolynomialize = function (u, p)
+            return [[cbrt(u[1]), sin(u[2] / 40)]]
+        end
+        polynomialize = function (u, p)
+            return [u[1] ^ 3, 40asin(u[2])]
+        end
+        fn = HomotopyNonlinearFunction(rhs; denominator, polynomialize, unpolynomialize)
+        prob = NonlinearProblem(fn, [1.0, 2.0], [2.0, 3.0, 4.0, 5.0])
+        sol2 = solve(prob, alg)
+        @test sol2 isa EnsembleSolution
+        @test sol2.converged
+        @test length(sol.u) == length(sol2.u)
+        for nlsol2 in sol2.u
+            @test any(nlsol -> isapprox(polynomialize(nlsol2.u, prob.p), nlsol.u; rtol = 1e-8), sol.u)
+        end
+
+        @testset "some invalid solutions" begin
+            prob3 = remake(prob; p = [2.0, 3.0, polynomialize(sol2.u[1].u, prob.p)...])
+            sol3 = solve(prob3, alg)
+            @test length(sol3.u) == length(sol2.u) - 1
+        end
+    end
+end
diff --git a/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl b/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl
index aedd9485a..50c1777bd 100644
--- a/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl
@@ -6,5 +6,8 @@ using Aqua
     @testset "Code quality (Aqua.jl)" begin
         Aqua.test_all(NonlinearSolveHomotopyContinuation)
     end
+    @testset "AllRoots" begin
+        include("allroots.jl")
+    end
     # Write your tests here.
 end

From c83c95a8a2664272f58a7ffae5c2b8cc4f53d8f3 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 6 Jan 2025 17:39:52 +0530
Subject: [PATCH 04/24] fix: fix jacobian implementation, add taylor
 implementation

---
 .../src/NonlinearSolveHomotopyContinuation.jl |  1 +
 .../src/interface_types.jl                    | 58 ++++++++++++++++++-
 .../src/solve.jl                              | 19 +++++-
 3 files changed, 74 insertions(+), 4 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
index 46b537e35..f9fb54fdf 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
@@ -6,6 +6,7 @@ using NonlinearSolveBase
 using SymbolicIndexingInterface
 using LinearAlgebra
 using ADTypes
+using TaylorSeries
 import CommonSolve
 import HomotopyContinuation as HC
 import DifferentiationInterface as DI
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
index 4b47d2212..4cdff84ad 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
@@ -9,6 +9,8 @@ struct Scalar <: HomotopySystemVariant end
     autodiff
     prep
     vars
+    taylorvars
+    jacobian_buffers
 end
 
 Base.size(sys::HomotopySystemWrapper) = (length(sys.prob.u0), length(sys.prob.u0))
@@ -36,7 +38,11 @@ function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Inp
         return
     end
 
-    DI.value_and_jacobian!(f.f, u, U, sys.prep, sys.autodiff, x, DI.Constant(p))
+    x_tmp, u_tmp, U_tmp = sys.jacobian_buffers
+    copyto!(x_tmp, x)
+    DI.value_and_jacobian!(f.f, u_tmp, U_tmp, sys.prep, sys.autodiff, x_tmp, DI.Constant(p))
+    copyto!(u, u_tmp)
+    copyto!(U, U_tmp)
 end
 
 function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing)
@@ -49,8 +55,11 @@ function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Out
         copyto!(U, j_tmp)
         return
     end
-    u_tmp, _ = DI.value_and_jacobian(f.f, U, sys.prep, sys.autodiff, x, DI.Constant(p))
+    x_tmp, U_tmp = sys.jacobian_buffers
+    copyto!(x_tmp, x)
+    u_tmp, _ = DI.value_and_jacobian!(f.f, U_tmp, sys.prep, sys.autodiff, x_tmp, DI.Constant(p))
     copyto!(u, u_tmp)
+    copyto!(U, U_tmp)
 end
 
 function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Scalar}, x, p = nothing)
@@ -60,11 +69,54 @@ function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Sca
         HC.ModelKit.evaluate!(u, sys, x, p)
         U[1] = f.jac(only(x), p)
     else
-        u[1], U[1] = DI.value_and_derivative(f.f, sys.prep, sys.autodiff, only(x), DI.Constant(p))
+        x = real(first(x))
+        u[1], U[1] = DI.value_and_derivative(f.f, sys.prep, sys.autodiff, x, DI.Constant(p))
     end
     return nothing
 end
 
+function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{Inplace}, x::HC.ModelKit.TaylorVector{M}, p = nothing) where {N, M}
+    f = sys.prob.f
+    p = parameter_values(sys.prob)
+    buffer, vars = sys.taylorvars
+    for i in eachindex(vars)
+        for j in 0:M
+            vars[i][j] = x[i, j + 1]
+        end
+    end
+    f.f(buffer, vars, p)
+    for i in eachindex(vars)
+        u[i] = ntuple(j -> buffer[i][j - 1], Val(N + 1))
+    end
+end
+
+function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{OutOfPlace}, x::HC.ModelKit.TaylorVector{M}, p = nothing) where {N, M}
+    f = sys.prob.f
+    p = parameter_values(sys.prob)
+    vars = sys.taylorvars
+    for i in eachindex(vars)
+        for j in 0:M
+            vars[i][j] = x[i, j + 1]
+        end
+    end
+    buffer = f.f(vars, p)
+    for i in eachindex(vars)
+        u[i] = ntuple(j -> buffer[i][j - 1], Val(N + 1))
+    end
+end
+
+function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{Scalar}, x::HC.ModelKit.TaylorVector{M}, p = nothing) where {N, M}
+    f = sys.prob.f
+    p = parameter_values(sys.prob)
+    var = sys.taylorvars
+    for i in 0:M
+        var[i] = x[1, i + 1]
+    end
+    taylor = f.f(var, p)
+    val = ntuple(i -> taylor[i - 1], Val(N + 1))
+    u[1] = val
+end
+
 @concrete struct GuessHomotopy <: HC.AbstractHomotopy
     sys <: HC.AbstractSystem
     fu0
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
index 2d30874fe..b37fa23b2 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
@@ -31,9 +31,26 @@ function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::Homoto
         HC.variables(:x, axes(u0)...)
     end
 
+    # TODO: Is there an upper bound for the order?
+    taylorvars = if isscalar
+        Taylor1(zeros(ComplexF64, 5), 4)
+    elseif iip
+        ([Taylor1(zeros(ComplexF64, 5), 4) for _ in u0], [Taylor1(zeros(ComplexF64, 5), 4) for _ in u0])
+    else
+        [Taylor1(zeros(ComplexF64, 5), 4) for _ in u0]
+    end
+
+    jacobian_buffers = if isscalar
+        nothing
+    elseif iip
+        (similar(u0), similar(u0), similar(u0, length(u0), length(u0)))
+    else
+        (similar(u0), similar(u0, length(u0), length(u0)))
+    end
+
     # HC-compatible system
     variant = iip ? Inplace : isscalar ? Scalar : OutOfPlace
-    hcsys = HomotopySystemWrapper{variant}(prob, alg.autodiff, prep, vars)
+    hcsys = HomotopySystemWrapper{variant}(prob, alg.autodiff, prep, vars, taylorvars, jacobian_buffers)
 
     return f, hcsys
 end

From fea76097fc603ea9ec853b894ab679aee39ff5ee Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 10 Jan 2025 10:49:41 +0530
Subject: [PATCH 05/24] fix: fix system evaluate and jacobian

---
 .../src/interface_types.jl                                 | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
index 4cdff84ad..86fdfc1cb 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
@@ -18,15 +18,18 @@ HC.ModelKit.variables(sys::HomotopySystemWrapper) = sys.vars
 
 function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{Inplace}, x, p = nothing)
     sys.prob.f.f(u, x, parameter_values(sys.prob))
+    return u
 end
 
 function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing)
     values = sys.prob.f.f(x, parameter_values(sys.prob))
     copyto!(u, values)
+    return u
 end
 
 function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{Scalar}, x, p = nothing)
     u[1] = sys.prob.f.f(only(x), parameter_values(sys.prob))
+    return u
 end
 
 function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Inplace}, x, p = nothing)
@@ -43,6 +46,7 @@ function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Inp
     DI.value_and_jacobian!(f.f, u_tmp, U_tmp, sys.prep, sys.autodiff, x_tmp, DI.Constant(p))
     copyto!(u, u_tmp)
     copyto!(U, U_tmp)
+    return u, U
 end
 
 function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing)
@@ -60,6 +64,7 @@ function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Out
     u_tmp, _ = DI.value_and_jacobian!(f.f, U_tmp, sys.prep, sys.autodiff, x_tmp, DI.Constant(p))
     copyto!(u, u_tmp)
     copyto!(U, U_tmp)
+    return u, U
 end
 
 function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Scalar}, x, p = nothing)
@@ -72,7 +77,7 @@ function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Sca
         x = real(first(x))
         u[1], U[1] = DI.value_and_derivative(f.f, sys.prep, sys.autodiff, x, DI.Constant(p))
     end
-    return nothing
+    return u, U
 end
 
 function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{Inplace}, x::HC.ModelKit.TaylorVector{M}, p = nothing) where {N, M}

From a048ca009458027a46ce579ee8885ac5dc68a44f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 10 Jan 2025 10:50:36 +0530
Subject: [PATCH 06/24] fix: fix inplace system taylor

---
 .../src/interface_types.jl                       | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
index 86fdfc1cb..81ddcbf91 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
@@ -85,14 +85,24 @@ function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWra
     p = parameter_values(sys.prob)
     buffer, vars = sys.taylorvars
     for i in eachindex(vars)
-        for j in 0:M
+        for j in 0:M-1
             vars[i][j] = x[i, j + 1]
         end
+        for j in M:4
+            vars[i][j] = zero(vars[i][j])
+        end
     end
     f.f(buffer, vars, p)
-    for i in eachindex(vars)
-        u[i] = ntuple(j -> buffer[i][j - 1], Val(N + 1))
+    if u isa Vector
+        for i in eachindex(vars)
+            u[i] = buffer[i][N]
+        end
+    else
+        for i in eachindex(vars)
+            u[i] = ntuple(j -> buffer[i][j - 1], Val(N + 1))
+        end
     end
+    return u
 end
 
 function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{OutOfPlace}, x::HC.ModelKit.TaylorVector{M}, p = nothing) where {N, M}

From 80ff4a78d3b95e8ae87dcece0dd304b957039ea2 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 10 Jan 2025 10:51:27 +0530
Subject: [PATCH 07/24] fix: fix homotopy and implement taylor

---
 .../src/interface_types.jl                    | 34 +++++++++++++++++--
 1 file changed, 32 insertions(+), 2 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
index 81ddcbf91..0f86982e5 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
@@ -135,20 +135,50 @@ end
 @concrete struct GuessHomotopy <: HC.AbstractHomotopy
     sys <: HC.AbstractSystem
     fu0
+    taylorbuffer::HC.ModelKit.TaylorVector{5, ComplexF64}
+end
+
+function GuessHomotopy(sys, fu0)
+    return GuessHomotopy(sys, fu0, HC.ModelKit.TaylorVector{5}(ComplexF64, length(fu0)))
 end
 
 Base.size(h::GuessHomotopy) = size(h.sys)
 
+# H(x, t) = (1 - t) * F(x) + t * (F(x) - F(x0))
+#         = F(x) - t * F(x0)
 function HC.ModelKit.evaluate!(u, h::GuessHomotopy, x, t, p = nothing)
     HC.ModelKit.evaluate!(u, h.sys, x, p)
     @inbounds for i in eachindex(u)
-        u[i] += (t + 1) * h.fu0[i]
+        u[i] -= t * h.fu0[i]
     end
+    return u
 end
 
 function HC.ModelKit.evaluate_and_jacobian!(u, U, h::GuessHomotopy, x, t, p = nothing)
     HC.ModelKit.evaluate_and_jacobian!(u, U, h.sys, x, p)
     @inbounds for i in eachindex(u)
-        u[i] += (t + 1) * h.fu0[i]
+        u[i] -= t * h.fu0[i]
+    end
+    return u, U
+end
+
+HC.ModelKit.taylor!(u, v::Val{N}, H::GuessHomotopy, tx, t, incremental::Bool) where {N} =
+    HC.ModelKit.taylor!(u, v, H, tx, t)
+
+function HC.ModelKit.taylor!(u, ::Val{1}, h::GuessHomotopy, x::AbstractVector{<:Number}, t, p = nothing)
+    HC.ModelKit.evaluate!(u, h.sys, x, p)
+    @inbounds for i in eachindex(u)
+        u[i] -= h.fu0[i]
+    end
+    return u
+end
+
+function HC.ModelKit.taylor!(u, ::Val{N}, h::GuessHomotopy, x, t, p = nothing) where {N}
+    HC.ModelKit.taylor!(h.taylorbuffer, Val(N), h.sys, x, p)
+    @inbounds for i in eachindex(u)
+        h.taylorbuffer[i, 1] -= t * h.fu0[i]
+        h.taylorbuffer[i, 2] -= h.fu0[i]
+        u[i] = h.taylorbuffer[i, N + 1]
     end
+    return u
 end

From 48d0f7a3927efe8c7e3a2ecfaadc42053f9674de Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 10 Jan 2025 10:51:44 +0530
Subject: [PATCH 08/24] fix: fix jacobian building for systems

---
 lib/NonlinearSolveHomotopyContinuation/src/solve.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
index b37fa23b2..e03957060 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
@@ -17,11 +17,11 @@ function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::Homoto
         prep = nothing
     elseif iip
         # prepare a DI jacobian if not
-        prep = DI.prepare_jacobian(f, copy(state_values(prob)), alg.autodiff, u0, DI.Constant(p))
+        prep = DI.prepare_jacobian(prob.f.f, copy(state_values(prob)), alg.autodiff, u0, DI.Constant(p))
     elseif isscalar
-        prep = DI.prepare_derivative(f, alg.autodiff, u0, DI.Constant(p))
+        prep = DI.prepare_derivative(prob.f.f, alg.autodiff, u0, DI.Constant(p))
     else
-        prep = DI.prepare_jacobian(f, alg.autodiff, u0, DI.Constant(p))
+        prep = DI.prepare_jacobian(prob.f.f, alg.autodiff, u0, DI.Constant(p))
     end
 
     # variables for HC to use

From af30f5e5350b6396006cc214061ad517aeebb30c Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 10 Jan 2025 10:52:08 +0530
Subject: [PATCH 09/24] fix: fix single root solve implementation

---
 lib/NonlinearSolveHomotopyContinuation/src/solve.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
index e03957060..38dea0eda 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
@@ -113,7 +113,10 @@ function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{f
     fu0 = NonlinearSolveBase.Utils.evaluate_f(prob, u0_p)
 
     homotopy = GuessHomotopy(hcsys, fu0)
-    orig_sol = HC.solve(homotopy, [fu0]; alg.kwargs..., kwargs...)
+    if u0_p isa Number
+        u0_p = [u0_p]
+    end
+    orig_sol = HC.solve(homotopy, [u0_p]; alg.kwargs..., kwargs...)
     realsols = HC.results(orig_sol; only_real = true)
 
     # no real solutions or infeasible solution

From 2026034c46f1d00637bc5b4188f433d71570e1f6 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 10 Jan 2025 10:52:47 +0530
Subject: [PATCH 10/24] build: add TaylorSeries as a dependency

---
 lib/NonlinearSolveHomotopyContinuation/Project.toml | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/Project.toml b/lib/NonlinearSolveHomotopyContinuation/Project.toml
index 6203ce640..5306c9ce9 100644
--- a/lib/NonlinearSolveHomotopyContinuation/Project.toml
+++ b/lib/NonlinearSolveHomotopyContinuation/Project.toml
@@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
+TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
 [compat]
 ADTypes = "1.11.0"
@@ -24,11 +25,13 @@ LinearAlgebra = "1.11.0"
 NonlinearSolveBase = "1.3.3"
 SciMLBase = "2.68.1"
 SymbolicIndexingInterface = "0.3.36"
+TaylorSeries = "0.18.2"
 julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test"]
+test = ["Aqua", "Test", "NonlinearSolve"]

From db21cf89bfa41b5ac330406bb84d7158eb4a98cd Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sun, 26 Jan 2025 15:17:01 +0530
Subject: [PATCH 11/24] fix: fix jacobian and taylor implementations

---
 .../src/interface_types.jl                    | 184 ++++++++++--------
 .../src/solve.jl                              |  56 +++---
 2 files changed, 138 insertions(+), 102 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
index 0f86982e5..db8bf22b5 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
@@ -4,8 +4,38 @@ struct Inplace <: HomotopySystemVariant end
 struct OutOfPlace <: HomotopySystemVariant end
 struct Scalar <: HomotopySystemVariant end
 
+@concrete struct ComplexJacobianWrapper{variant <: HomotopySystemVariant}
+    f
+end
+
+function (cjw::ComplexJacobianWrapper{Inplace})(u::AbstractVector{T}, x::AbstractVector{T}, p) where {T}
+    x = reinterpret(Complex{T}, x)
+    u = reinterpret(Complex{T}, u)
+    cjw.f(u, x, p)
+    u = parent(u)
+    return u
+end
+
+function (cjw::ComplexJacobianWrapper{OutOfPlace})(u::AbstractVector{T}, x::AbstractVector{T}, p) where {T}
+    x = reinterpret(Complex{T}, x)
+    u_tmp = cjw.f(x, p)
+    u_tmp = reinterpret(T, u_tmp)
+    copyto!(u, u_tmp)
+    return u
+end
+
+function (cjw::ComplexJacobianWrapper{Scalar})(u::AbstractVector{T}, x::AbstractVector{T}, p) where {T}
+    x = reinterpret(Complex{T}, x)
+    u_tmp = cjw.f(x[1], p)
+    u[1] = real(u_tmp)
+    u[2] = imag(u_tmp)
+    return u
+end
+
 @concrete struct HomotopySystemWrapper{variant <: HomotopySystemVariant} <: HC.AbstractSystem
-    prob
+    f
+    jac
+    p
     autodiff
     prep
     vars
@@ -13,77 +43,70 @@ struct Scalar <: HomotopySystemVariant end
     jacobian_buffers
 end
 
-Base.size(sys::HomotopySystemWrapper) = (length(sys.prob.u0), length(sys.prob.u0))
+Base.size(sys::HomotopySystemWrapper) = (length(sys.vars), length(sys.vars))
 HC.ModelKit.variables(sys::HomotopySystemWrapper) = sys.vars
 
 function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{Inplace}, x, p = nothing)
-    sys.prob.f.f(u, x, parameter_values(sys.prob))
+    sys.f(u, x, sys.p)
     return u
 end
 
 function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing)
-    values = sys.prob.f.f(x, parameter_values(sys.prob))
+    values = sys.f(x, sys.p)
     copyto!(u, values)
     return u
 end
 
 function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{Scalar}, x, p = nothing)
-    u[1] = sys.prob.f.f(only(x), parameter_values(sys.prob))
+    u[1] = sys.f(x[1], sys.p)
     return u
 end
 
 function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Inplace}, x, p = nothing)
-    f = sys.prob.f
-    p = parameter_values(sys.prob)
-    if SciMLBase.has_jac(f)
-        f.f(u, x, p)
-        f.jac(U, x, p)
-        return
-    end
-
-    x_tmp, u_tmp, U_tmp = sys.jacobian_buffers
-    copyto!(x_tmp, x)
-    DI.value_and_jacobian!(f.f, u_tmp, U_tmp, sys.prep, sys.autodiff, x_tmp, DI.Constant(p))
-    copyto!(u, u_tmp)
-    copyto!(U, U_tmp)
+    p = sys.p
+    sys.f(u, x, p)
+    sys.jac(U, x, p)
     return u, U
 end
 
 function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing)
-    f = sys.prob.f
-    p = parameter_values(sys.prob)
-    if SciMLBase.has_jac(f)
-        u_tmp = f.f(x, p)
-        copyto!(u, u_tmp)
-        j_tmp = f.jac(U, x, p)
-        copyto!(U, j_tmp)
-        return
-    end
-    x_tmp, U_tmp = sys.jacobian_buffers
-    copyto!(x_tmp, x)
-    u_tmp, _ = DI.value_and_jacobian!(f.f, U_tmp, sys.prep, sys.autodiff, x_tmp, DI.Constant(p))
+    p = sys.p
+    u_tmp = sys.f(x, p)
     copyto!(u, u_tmp)
-    copyto!(U, U_tmp)
+    j_tmp = sys.jac(x, p)
+    copyto!(U, j_tmp)
     return u, U
 end
 
 function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Scalar}, x, p = nothing)
-    f = sys.prob.f
-    p = parameter_values(sys.prob)
-    if SciMLBase.has_jac(f)
-        HC.ModelKit.evaluate!(u, sys, x, p)
-        U[1] = f.jac(only(x), p)
-    else
-        x = real(first(x))
-        u[1], U[1] = DI.value_and_derivative(f.f, sys.prep, sys.autodiff, x, DI.Constant(p))
-    end
+    p = sys.p
+    u[1] = sys.f(x[1], p)
+    U[1] = sys.jac(x[1], p)
     return u, U
 end
 
-function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{Inplace}, x::HC.ModelKit.TaylorVector{M}, p = nothing) where {N, M}
-    f = sys.prob.f
-    p = parameter_values(sys.prob)
-    buffer, vars = sys.taylorvars
+for V in (Inplace, OutOfPlace, Scalar)
+    @eval function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{$V, F, J}, x, p = nothing) where {F, J <: ComplexJacobianWrapper}
+        p = sys.p
+        U_tmp = sys.jacobian_buffers
+        x = reinterpret(Float64, x)
+        u = reinterpret(Float64, u)
+        DI.value_and_jacobian!(sys.jac, u, U_tmp, sys.prep, sys.autodiff, x, DI.Constant(p))
+        U = reinterpret(Float64, U)
+        @inbounds for j in axes(U, 2)
+            jj = 2j - 1
+            for i in axes(U, 1)
+                U[i, j] = U_tmp[i, jj]
+            end
+        end
+        u = parent(u)
+        U = parent(U)
+
+        return u, U
+    end
+end
+
+function update_taylorvars_from_taylorvector!(vars, x::HC.ModelKit.TaylorVector{M}) where {M}
     for i in eachindex(vars)
         for j in 0:M-1
             vars[i][j] = x[i, j + 1]
@@ -92,44 +115,57 @@ function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWra
             vars[i][j] = zero(vars[i][j])
         end
     end
-    f.f(buffer, vars, p)
-    if u isa Vector
-        for i in eachindex(vars)
-            u[i] = buffer[i][N]
-        end
-    else
-        for i in eachindex(vars)
-            u[i] = ntuple(j -> buffer[i][j - 1], Val(N + 1))
+end
+
+function update_taylorvars_from_taylorvector!(vars, x::AbstractVector)
+    for i in eachindex(vars)
+        vars[i][0] = x[i]
+        for j in 1:4
+            vars[i][j] = zero(vars[i][j])
         end
     end
-    return u
 end
 
-function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{OutOfPlace}, x::HC.ModelKit.TaylorVector{M}, p = nothing) where {N, M}
-    f = sys.prob.f
-    p = parameter_values(sys.prob)
-    vars = sys.taylorvars
+function update_maybe_taylorvector_from_taylorvars!(u::Vector, vars, buffer, ::Val{N}) where {N}
     for i in eachindex(vars)
-        for j in 0:M
-            vars[i][j] = x[i, j + 1]
-        end
+        u[i] = buffer[i][N]
     end
-    buffer = f.f(vars, p)
+end
+
+function update_maybe_taylorvector_from_taylorvars!(u::HC.ModelKit.TaylorVector, vars, buffer, ::Val{N}) where {N}
     for i in eachindex(vars)
         u[i] = ntuple(j -> buffer[i][j - 1], Val(N + 1))
     end
 end
 
-function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{Scalar}, x::HC.ModelKit.TaylorVector{M}, p = nothing) where {N, M}
-    f = sys.prob.f
-    p = parameter_values(sys.prob)
+function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{Inplace}, x, p = nothing) where {N}
+    f = sys.f
+    p = sys.p
+    buffer, vars = sys.taylorvars
+    update_taylorvars_from_taylorvector!(vars, x)
+    f(buffer, vars, p)
+    update_maybe_taylorvector_from_taylorvars!(u, vars, buffer, Val(N))
+    return u
+end
+
+function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing) where {N}
+    f = sys.f
+    p = sys.p
+    vars = sys.taylorvars
+    update_taylorvars_from_taylorvector!(vars, x)
+    buffer = f(vars, p)
+    update_maybe_taylorvector_from_taylorvars!(u, vars, buffer, Val(N))
+    return u
+end
+
+function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{Scalar}, x, p = nothing) where {N}
+    f = sys.f
+    p = sys.p
     var = sys.taylorvars
-    for i in 0:M
-        var[i] = x[1, i + 1]
-    end
-    taylor = f.f(var, p)
-    val = ntuple(i -> taylor[i - 1], Val(N + 1))
-    u[1] = val
+    update_taylorvars_from_taylorvector!((var,), x)
+    buffer = f(var, p)
+    update_maybe_taylorvector_from_taylorvars!(u, (var,), (buffer,), Val(N))
+    return u
 end
 
 @concrete struct GuessHomotopy <: HC.AbstractHomotopy
@@ -165,14 +201,6 @@ end
 HC.ModelKit.taylor!(u, v::Val{N}, H::GuessHomotopy, tx, t, incremental::Bool) where {N} =
     HC.ModelKit.taylor!(u, v, H, tx, t)
 
-function HC.ModelKit.taylor!(u, ::Val{1}, h::GuessHomotopy, x::AbstractVector{<:Number}, t, p = nothing)
-    HC.ModelKit.evaluate!(u, h.sys, x, p)
-    @inbounds for i in eachindex(u)
-        u[i] -= h.fu0[i]
-    end
-    return u
-end
-
 function HC.ModelKit.taylor!(u, ::Val{N}, h::GuessHomotopy, x, t, p = nothing) where {N}
     HC.ModelKit.taylor!(h.taylorbuffer, Val(N), h.sys, x, p)
     @inbounds for i in eachindex(u)
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
index 38dea0eda..76ed740fd 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
@@ -11,17 +11,22 @@ function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::Homoto
     isscalar = u0 isa Number
     iip = SciMLBase.isinplace(prob)
 
+    variant = iip ? Inplace : isscalar ? Scalar : OutOfPlace
+
     # jacobian handling
-    if SciMLBase.has_jac(f)
+    if SciMLBase.has_jac(f.f)
         # use if present
         prep = nothing
-    elseif iip
-        # prepare a DI jacobian if not
-        prep = DI.prepare_jacobian(prob.f.f, copy(state_values(prob)), alg.autodiff, u0, DI.Constant(p))
-    elseif isscalar
-        prep = DI.prepare_derivative(prob.f.f, alg.autodiff, u0, DI.Constant(p))
+        jac = f.jac
     else
-        prep = DI.prepare_jacobian(prob.f.f, alg.autodiff, u0, DI.Constant(p))
+        # prepare a DI jacobian if not
+        jac = ComplexJacobianWrapper{variant}(f.f.f)
+        tmp = if isscalar
+            Vector{Float64}(undef, 2)
+        else
+            similar(u0, Float64, 2length(u0))
+        end
+        prep = DI.prepare_jacobian(jac, tmp, alg.autodiff, copy(tmp), DI.Constant(p))
     end
 
     # variables for HC to use
@@ -41,16 +46,13 @@ function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::Homoto
     end
 
     jacobian_buffers = if isscalar
-        nothing
-    elseif iip
-        (similar(u0), similar(u0), similar(u0, length(u0), length(u0)))
+        Matrix{Float64}(undef, 2, 2)
     else
-        (similar(u0), similar(u0, length(u0), length(u0)))
+        similar(u0, Float64, 2length(u0), 2length(u0))
     end
 
     # HC-compatible system
-    variant = iip ? Inplace : isscalar ? Scalar : OutOfPlace
-    hcsys = HomotopySystemWrapper{variant}(prob, alg.autodiff, prep, vars, taylorvars, jacobian_buffers)
+    hcsys = HomotopySystemWrapper{variant}(f.f.f, jac, p, alg.autodiff, prep, vars, taylorvars, jacobian_buffers)
 
     return f, hcsys
 end
@@ -113,14 +115,14 @@ function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{f
     fu0 = NonlinearSolveBase.Utils.evaluate_f(prob, u0_p)
 
     homotopy = GuessHomotopy(hcsys, fu0)
-    if u0_p isa Number
-        u0_p = [u0_p]
+    orig_sol = HC.solve(homotopy, u0_p isa Number ? [[u0_p]] : [u0_p]; alg.kwargs..., kwargs...)
+    realsols = map(res -> res.solution, HC.results(orig_sol; only_real = true))
+    if u0 isa Number
+        realsols = map(only, realsols)
     end
-    orig_sol = HC.solve(homotopy, [u0_p]; alg.kwargs..., kwargs...)
-    realsols = HC.results(orig_sol; only_real = true)
 
     # no real solutions or infeasible solution
-    if isempty(realsols) || any(<=(denominator_abstol), f.denominator(real.(only(realsols).solution), p))
+    if isempty(realsols) || any(<=(denominator_abstol), map(abs, f.denominator(real.(only(realsols)), p)))
         retcode = if isempty(realsols)
             SciMLBase.ReturnCode.ConvergenceFailure
         else
@@ -129,15 +131,21 @@ function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{f
         resid = NonlinearSolveBase.Utils.evaluate_f(prob, u0)
         return SciMLBase.build_solution(prob, alg, u0, resid; retcode, original = orig_sol)
     end
-    realsol = only(realsols)
+
+    realsol = real(only(realsols))
     T = eltype(u0)
-    validsols = f.unpolynomialize(realsol.solution, p)
-    sol, idx = findmin(validsols) do sol
-        norm(sol - u0)
+    validsols = f.unpolynomialize(realsol, p)
+    _, idx = findmin(validsols) do sol
+        norm(sol - u0_p)
+    end
+    u = map(real, validsols[idx])
+
+    if u0 isa Number
+        u = only(u)
     end
-    resid = NonlinearSolveBase.Utils.evaluate_f(prob, u0)
+    resid = NonlinearSolveBase.Utils.evaluate_f(prob, u)
 
     retcode = SciMLBase.ReturnCode.Success
-    return SciMLBase.build_solution(prob, alg, u0, resid; retcode, original = orig_sol)
+    return SciMLBase.build_solution(prob, alg, u, resid; retcode, original = orig_sol)
 end
 

From 75aabbb973dfea1719dd905f6a49a1d9540c748f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sun, 26 Jan 2025 15:23:06 +0530
Subject: [PATCH 12/24] test: add tests for single roots, fix allroots tests

---
 .../test/allroots.jl                          |  45 +++++--
 .../test/runtests.jl                          |   4 +-
 .../test/single_root.jl                       | 125 ++++++++++++++++++
 3 files changed, 161 insertions(+), 13 deletions(-)
 create mode 100644 lib/NonlinearSolveHomotopyContinuation/test/single_root.jl

diff --git a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
index 581ea0b22..18a768574 100644
--- a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
@@ -5,16 +5,24 @@ using SciMLBase: NonlinearSolution
 alg = HomotopyContinuationJL{true}(; threading = false)
 
 @testset "scalar u" begin
-    @testset "`NonlinearProblem`" begin
-        prob = NonlinearProblem(1.0, 2.0) do u, p
-            return u * u - 2p * u + p * p
-        end
+    rhs = function (u, p)
+        return u * u - 2p * u + p * p
+    end
+    jac = function (u, p)
+        return 2u - 2p
+    end
+    @testset "`NonlinearProblem` - $name" for (jac, name) in [(nothing, "no jac"), (jac, "jac")]
+        fn = NonlinearFunction(rhs; jac)
+        prob = NonlinearProblem(fn, 1.0, 2.0)
         sol = solve(prob, alg)
 
         @test sol isa EnsembleSolution
-        @test length(sol) == 1
-        @test sol.u[1] isa NonlinearSolution
-        @test sol.u[1][1] ≈ 2.0
+        @test sol.converged
+        for nlsol in sol
+            @test nlsol isa NonlinearSolution
+            @test SciMLBase.successful_retcode(nlsol)
+            @test nlsol.u[1] ≈ 2.0 atol = 1e-7
+        end
 
         @testset "no real solutions" begin
             prob = NonlinearProblem(1.0, 0.5) do u, p
@@ -63,19 +71,31 @@ alg = HomotopyContinuationJL{true}(; threading = false)
     end
 end
 
-function f!(du, u, p)
+f! = function (du, u, p)
     du[1] = u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1
     du[2] = u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]
 end
 
-function f(u, p)
+f = function (u, p)
     [u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1, u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]]
 end
 
-@testset "vector u - $iip" for (rhs, iip) in [(f, "oop"), (f!, "iip")]
+jac! = function (j, u, p)
+    j[1, 1] = 2u[1]
+    j[1, 2] = -p[1] + 3 * u[2]^2
+    j[2, 1] = 2 * p[2] * u[2]
+    j[2, 2] = 3 * u[2]^2 + 2 * p[2] * u[1] + 1
+end
+jac = function (u, p)
+    [2u[1] -p[1] + 3 * u[2]^2;
+     2 * p[2] * u[2] 3 * u[2]^2 + 2 * p[2] * u[1] + 1]
+end
+
+@testset "vector u - $name" for (rhs, jac, name) in [(f, nothing, "oop"), (f, jac, "oop + jac"), (f!, nothing, "iip"), (f!, jac!, "iip + jac")]
     sol = nothing
     @testset "`NonlinearProblem`" begin
-        prob = NonlinearProblem(rhs, [1.0, 2.0], [2.0, 3.0])
+        fn = NonlinearFunction(rhs; jac)
+        prob = NonlinearProblem(fn, [1.0, 2.0], [2.0, 3.0])
         sol = solve(prob, alg)
         @test sol isa EnsembleSolution
         @test sol.converged
@@ -103,7 +123,8 @@ end
         polynomialize = function (u, p)
             return [u[1] ^ 3, 40asin(u[2])]
         end
-        fn = HomotopyNonlinearFunction(rhs; denominator, polynomialize, unpolynomialize)
+        nlfn = NonlinearFunction(rhs; jac)
+        fn = HomotopyNonlinearFunction(nlfn; denominator, polynomialize, unpolynomialize)
         prob = NonlinearProblem(fn, [1.0, 2.0], [2.0, 3.0, 4.0, 5.0])
         sol2 = solve(prob, alg)
         @test sol2 isa EnsembleSolution
diff --git a/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl b/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl
index 50c1777bd..94bc8bddd 100644
--- a/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl
@@ -9,5 +9,7 @@ using Aqua
     @testset "AllRoots" begin
         include("allroots.jl")
     end
-    # Write your tests here.
+    @testset "Single Root" begin
+        include("single_root.jl")
+    end
 end
diff --git a/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl b/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl
new file mode 100644
index 000000000..3cd5b81cf
--- /dev/null
+++ b/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl
@@ -0,0 +1,125 @@
+using NonlinearSolve
+using NonlinearSolveHomotopyContinuation
+using SciMLBase: NonlinearSolution
+
+alg = HomotopyContinuationJL{false}(; threading = false)
+
+@testset "scalar u" begin
+    rhs = function (u, p)
+        return (u - 3.0) * (u - p)
+    end
+    jac = function (u, p)
+        return 2u - (p + 3)
+    end
+    @testset "`NonlinearProblem` - $name" for (jac, name) in [(nothing, "no jac"), (jac, "jac")]
+        fn = NonlinearFunction(rhs; jac)
+        prob = NonlinearProblem(fn, 1.0, 2.0)
+        sol = solve(prob, alg)
+
+        @test sol isa NonlinearSolution
+        @test sol.u ≈ 2.0 atol = 1e-10
+
+        @testset "no real solutions" begin
+            prob = NonlinearProblem(1.0, 1.0) do u, p
+                return u * u - 2p * u + p
+            end
+            sol = solve(prob, alg)
+            @test sol.retcode == SciMLBase.ReturnCode.ConvergenceFailure
+        end
+    end
+
+    @testset "`HomotopyContinuationFunction`" begin
+        denominator = function (u, p)
+            return [u - 0.7]
+        end
+        polynomialize = function (u, p)
+            return sin(u)
+        end
+        unpolynomialize = function (u, p)
+            return [asin(u)]
+        end
+        fn = HomotopyNonlinearFunction(; denominator, polynomialize, unpolynomialize) do u, p
+            return (u - p[1]) * (u - p[2])
+        end
+        prob = NonlinearProblem(fn, 0.0, [0.5, 0.7])
+
+        sol = solve(prob, alg)
+        @test sin(sol.u[1]) ≈ 0.5 atol=1e-10
+
+        @testset "no valid solutions" begin
+            prob2 = remake(prob; p = [0.7, 0.7])
+            sol2 = solve(prob2, alg)
+            @test sol2.retcode == SciMLBase.ReturnCode.Infeasible
+        end
+
+        @testset "closest root" begin
+            prob3 = remake(prob; p = [0.5, 0.6], u0 = asin(0.4))
+            sol3 = solve(prob3, alg)
+            @test sin(sol3.u) ≈ 0.5 atol = 1e-10
+            prob4 = remake(prob3; u0 = asin(0.7))
+            sol4 = solve(prob4, alg)
+            @test sin(sol4.u) ≈ 0.6 atol = 1e-10
+        end
+    end
+end
+
+f! = function (du, u, p)
+    du[1] = u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1
+    du[2] = u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]
+end
+
+f = function (u, p)
+    [u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1, u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]]
+end
+
+jac! = function (j, u, p)
+    j[1, 1] = 2u[1]
+    j[1, 2] = -p[1] + 3 * u[2]^2
+    j[2, 1] = 2 * p[2] * u[2]
+    j[2, 2] = 3 * u[2]^2 + 2 * p[2] * u[1] + 1
+end
+
+jac = function (u, p)
+    [2u[1] -p[1] + 3 * u[2]^2;
+     2 * p[2] * u[2] 3 * u[2]^2 + 2 * p[2] * u[1] + 1]
+end
+
+@testset "vector u - $name" for (rhs, jac, name) in [(f, nothing, "oop"), (f, jac, "oop + jac"), (f!, nothing, "iip"), (f!, jac!, "iip + jac")]
+    @testset "`NonlinearProblem`" begin
+        fn = NonlinearFunction(rhs; jac)
+        prob = NonlinearProblem(fn, [1.0, 2.0], [2.0, 3.0])
+        sol = solve(prob, alg)
+        @test SciMLBase.successful_retcode(sol)
+        @test f(sol.u, prob.p) ≈ [0.0, 0.0] atol = 1e-10
+
+        @testset "no real solutions" begin
+            prob2 = remake(prob; p = zeros(2))
+            sol2 = solve(prob2, alg)
+            @test !SciMLBase.successful_retcode(sol2)
+        end
+    end
+
+    @testset "`HomotopyNonlinearFunction`" begin
+        denominator = function (u, p)
+            return [u[1] - p[3], u[2] - p[4]]
+        end
+        unpolynomialize = function (u, p)
+            return [[cbrt(u[1]), sin(u[2] / 40)]]
+        end
+        polynomialize = function (u, p)
+            return [u[1] ^ 3, 40asin(u[2])]
+        end
+        nlfn = NonlinearFunction(rhs; jac)
+        fn = HomotopyNonlinearFunction(nlfn; denominator, polynomialize, unpolynomialize)
+        prob = NonlinearProblem(fn, [1.0, 1.0], [2.0, 3.0, 4.0, 5.0])
+        sol = solve(prob, alg)
+        @test SciMLBase.successful_retcode(sol)
+        @test f(polynomialize(sol.u, prob.p), prob.p) ≈ zeros(2) atol = 1e-10
+
+        @testset "some invalid solutions" begin
+            prob2 = remake(prob; p = [2.0, 3.0, polynomialize(sol.u, prob.p)...])
+            sol2 = solve(prob2, alg)
+            @test !SciMLBase.successful_retcode(sol2)
+        end
+    end
+end

From d7972af18a2925bbe09457b6bed175cea29d9303 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sun, 26 Jan 2025 15:23:14 +0530
Subject: [PATCH 13/24] build: add compat entries

---
 lib/NonlinearSolveHomotopyContinuation/Project.toml | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/Project.toml b/lib/NonlinearSolveHomotopyContinuation/Project.toml
index 5306c9ce9..f0d65f647 100644
--- a/lib/NonlinearSolveHomotopyContinuation/Project.toml
+++ b/lib/NonlinearSolveHomotopyContinuation/Project.toml
@@ -16,16 +16,19 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
 [compat]
+Aqua = "0.8"
 ADTypes = "1.11.0"
 CommonSolve = "0.2.4"
 ConcreteStructs = "0.2.3"
 DifferentiationInterface = "0.6.27"
 HomotopyContinuation = "2.12.0"
 LinearAlgebra = "1.11.0"
+NonlinearSolve = "4"
 NonlinearSolveBase = "1.3.3"
-SciMLBase = "2.68.1"
+SciMLBase = "2.71"
 SymbolicIndexingInterface = "0.3.36"
 TaylorSeries = "0.18.2"
+Test = "1.10"
 julia = "1.10"
 
 [extras]

From 4a8def38480bccff2899585bf8511c1b6a370aa4 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sun, 26 Jan 2025 15:59:21 +0530
Subject: [PATCH 14/24] docs: add docstrings to
 NonlinearSolveHomotopyContinuation.jl

---
 .../Project.toml                              |  4 +-
 .../src/NonlinearSolveHomotopyContinuation.jl | 30 ++++++++
 .../src/interface_types.jl                    | 76 ++++++++++++++++++-
 .../src/solve.jl                              |  7 +-
 4 files changed, 113 insertions(+), 4 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/Project.toml b/lib/NonlinearSolveHomotopyContinuation/Project.toml
index f0d65f647..517d7b27d 100644
--- a/lib/NonlinearSolveHomotopyContinuation/Project.toml
+++ b/lib/NonlinearSolveHomotopyContinuation/Project.toml
@@ -8,6 +8,7 @@ ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
@@ -16,11 +17,12 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
 [compat]
-Aqua = "0.8"
 ADTypes = "1.11.0"
+Aqua = "0.8"
 CommonSolve = "0.2.4"
 ConcreteStructs = "0.2.3"
 DifferentiationInterface = "0.6.27"
+DocStringExtensions = "0.9.3"
 HomotopyContinuation = "2.12.0"
 LinearAlgebra = "1.11.0"
 NonlinearSolve = "4"
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
index f9fb54fdf..4513d3349 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
@@ -7,6 +7,7 @@ using SymbolicIndexingInterface
 using LinearAlgebra
 using ADTypes
 using TaylorSeries
+using DocStringExtensions
 import CommonSolve
 import HomotopyContinuation as HC
 import DifferentiationInterface as DI
@@ -15,6 +16,35 @@ using ConcreteStructs: @concrete
 
 export HomotopyContinuationJL, HomotopyNonlinearFunction
 
+"""
+    HomotopyContinuationJL{AllRoots}(; autodiff = true, kwargs...)
+    HomotopyContinuationJL(; kwargs...) = HomotopyContinuationJL{false}(; kwargs...)
+
+This algorithm is an interface to `HomotopyContinuation.jl`. It is only valid for
+fully determined polynomial systems. The `AllRoots` type parameter can be `true` or
+`false` and controls whether the solver will find all roots of the polynomial
+or a single root close to the initial guess provided to the `NonlinearProblem`.
+The polynomial function must allow complex numbers to be provided as the state.
+
+If `AllRoots` is `true`, the initial guess in the `NonlinearProblem` is ignored.
+The function must be traceable using HomotopyContinuation.jl's symbolic variables.
+Note that higher degree polynomials and systems with multiple unknowns can increase
+solve time significantly.
+
+If `AllRoots` is `false`, a single path is traced during the homotopy. The traced path
+depends on the initial guess provided to the `NonlinearProblem` being solved. This method
+does not require that the polynomial function is traceable via HomotopyContinuation.jl's
+symbolic variables.
+
+HomotopyContinuation.jl requires the jacobian of the system. In case a jacobian function
+is provided, it will be used. Otherwise, the `autodiff` keyword argument controls the
+autodiff method used to compute the jacobian. A value of `true` refers to
+`AutoForwardDiff` and `false` refers to `AutoFiniteDiff`. Alternate algorithms can be
+specified using ADTypes.jl.
+
+HomotopyContinuation.jl requires the taylor series of the polynomial system for the single
+root method. This is automatically computed using TaylorSeries.jl.
+"""
 @concrete struct HomotopyContinuationJL{AllRoots} <: NonlinearSolveBase.AbstractNonlinearSolveAlgorithm
     autodiff
     kwargs
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
index db8bf22b5..24f23baea 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
@@ -4,6 +4,16 @@ struct Inplace <: HomotopySystemVariant end
 struct OutOfPlace <: HomotopySystemVariant end
 struct Scalar <: HomotopySystemVariant end
 
+"""
+    $(TYPEDEF)
+
+A simple struct that wraps a polynomial function which takes complex input and returns
+complex output in a form that supports automatic differentiation. If the wrapped
+function if ``f: \\mathbb{C}^n \\rightarrow \\mathbb{C}^n`` then it is assumed
+the input arrays are real-valued and have length ``2n``. They are `reinterpret`ed
+into complex arrays and passed into the function. This struct has an in-place signature
+regardless of the signature of ``f``.
+"""
 @concrete struct ComplexJacobianWrapper{variant <: HomotopySystemVariant}
     f
 end
@@ -32,14 +42,49 @@ function (cjw::ComplexJacobianWrapper{Scalar})(u::AbstractVector{T}, x::Abstract
     return u
 end
 
+"""
+    $(TYPEDEF)
+
+A struct which implements the `HomotopyContinuation.AbstractSystem` interface for
+polynomial systems specified using `NonlinearProblem`.
+
+# Fields
+
+$(FIELDS)
+"""
 @concrete struct HomotopySystemWrapper{variant <: HomotopySystemVariant} <: HC.AbstractSystem
+    """
+    The wrapped polynomial function.
+    """
     f
+    """
+    The jacobian function, if provided to the `NonlinearProblem` being solved. Otherwise,
+    a `ComplexJacobianWrapper` wrapping `f` used for automatic differentiation.
+    """
     jac
+    """
+    The parameter object.
+    """
     p
+    """
+    The ADType for automatic differentiation.
+    """
     autodiff
+    """
+    The result from `DifferentiationInterface.prepare_jacobian`.
+    """
     prep
+    """
+    HomotopyContinuation.jl's symbolic variables for the system.
+    """
     vars
+    """
+    The `TaylorSeries.Taylor1` objects used to compute the taylor series of `f`.
+    """
     taylorvars
+    """
+    Preallocated intermediate buffers used for calculating the jacobian.
+    """
     jacobian_buffers
 end
 
@@ -168,9 +213,38 @@ function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWra
     return u
 end
 
+"""
+    $(TYPEDEF)
+
+A `HomotopyContinuation.AbstractHomotopy` which uses an inital guess ``x_0`` to construct
+the start system for the homotopy. The homotopy is
+
+```math
+\\begin{aligned}
+H(x, t) &= t * (F(x) - F(x_0)) + (1 - t) * F(x)
+        &= F(x) - t * F(x_0)
+\\end{aligned}
+```
+
+Where ``F: \\mathbb{R}^n \\rightarrow \\mathbb{R}^n`` is the polynomial system and
+``x_0 \\in \\mathbb{R}^n`` is the guess provided to the `NonlinearProblem`.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
 @concrete struct GuessHomotopy <: HC.AbstractHomotopy
+    """
+    The `HomotopyContinuation.AbstractSystem` to solve.
+    """
     sys <: HC.AbstractSystem
+    """
+    The residual of the initial guess, i.e. ``F(x_0)``.
+    """
     fu0
+    """
+    A preallocated `TaylorVector` for efficient `taylor!` implementation.
+    """
     taylorbuffer::HC.ModelKit.TaylorVector{5, ComplexF64}
 end
 
@@ -180,8 +254,6 @@ end
 
 Base.size(h::GuessHomotopy) = size(h.sys)
 
-# H(x, t) = (1 - t) * F(x) + t * (F(x) - F(x0))
-#         = F(x) - t * F(x0)
 function HC.ModelKit.evaluate!(u, h::GuessHomotopy, x, t, p = nothing)
     HC.ModelKit.evaluate!(u, h.sys, x, p)
     @inbounds for i in eachindex(u)
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
index 76ed740fd..913f66b44 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
@@ -1,3 +1,9 @@
+"""
+    $(TYPEDSIGNATURES)
+
+Create and return the appropriate `HomotopySystemWrapper` to use for solving the given
+`prob` with `alg`.
+"""
 function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::HomotopyContinuationJL)
     # cast to a `HomotopyNonlinearFunction`
     f = if prob.f isa HomotopyNonlinearFunction
@@ -36,7 +42,6 @@ function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::Homoto
         HC.variables(:x, axes(u0)...)
     end
 
-    # TODO: Is there an upper bound for the order?
     taylorvars = if isscalar
         Taylor1(zeros(ComplexF64, 5), 4)
     elseif iip

From c2aa81239ad384b9970b5e63bf9fad19baec9c9f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sun, 26 Jan 2025 16:50:26 +0530
Subject: [PATCH 15/24] ci: add NonlinearSolveHomotopyContinuation.jl to CI

---
 .../CI_NonlinearSolveHomotopyContinuation.yml | 108 ++++++++++++++++++
 1 file changed, 108 insertions(+)
 create mode 100644 .github/workflows/CI_NonlinearSolveHomotopyContinuation.yml

diff --git a/.github/workflows/CI_NonlinearSolveHomotopyContinuation.yml b/.github/workflows/CI_NonlinearSolveHomotopyContinuation.yml
new file mode 100644
index 000000000..5ccc6fd1f
--- /dev/null
+++ b/.github/workflows/CI_NonlinearSolveHomotopyContinuation.yml
@@ -0,0 +1,108 @@
+name: CI (NonlinearSolveHomotopyContinuation)
+
+on:
+  pull_request:
+    branches:
+      - master
+    paths:
+      - "lib/NonlinearSolveHomotopyContinuation/**"
+      - ".github/workflows/CI_NonlinearSolveHomotopyContinuation.yml"
+      - "lib/NonlinearSolveBase/**"
+  push:
+    branches:
+      - master
+
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - "1.10"
+          - "1"
+        os:
+          - ubuntu-latest
+          - macos-latest
+          - windows-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+      - uses: actions/cache@v4
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - name: "Install Dependencies and Run Tests"
+        run: |
+          import Pkg
+          Pkg.Registry.update()
+          # Install packages present in subdirectories
+          dev_pks = Pkg.PackageSpec[]
+          for path in ("lib/NonlinearSolveBase",)
+              push!(dev_pks, Pkg.PackageSpec(; path))
+          end
+          Pkg.develop(dev_pks)
+          Pkg.instantiate()
+          Pkg.test(; coverage="user")
+        shell: julia --color=yes --code-coverage=user --depwarn=yes --project=lib/NonlinearSolveHomotopyContinuation {0}
+      - uses: julia-actions/julia-processcoverage@v1
+        with:
+          directories: lib/NonlinearSolveBase/src,lib/NonlinearSolveBase/ext,lib/NonlinearSolveHomotopyContinuation/src
+      - uses: codecov/codecov-action@v5
+        with:
+          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          verbose: true
+          fail_ci_if_error: false
+
+  downgrade:
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - "1.10"
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+      - uses: julia-actions/julia-downgrade-compat@v1
+        with:
+          skip: NonlinearSolveBase, SciMLJacobianOperators
+      - name: "Install Dependencies and Run Tests"
+        run: |
+          import Pkg
+          Pkg.Registry.update()
+          # Install packages present in subdirectories
+          dev_pks = Pkg.PackageSpec[]
+          for path in ("lib/NonlinearSolveBase",)
+              push!(dev_pks, Pkg.PackageSpec(; path))
+          end
+          Pkg.develop(dev_pks)
+          Pkg.instantiate()
+          Pkg.test(; coverage="user")
+        shell: julia --color=yes --code-coverage=user --depwarn=yes --project=lib/NonlinearSolveHomotopyContinuation {0}
+      - uses: julia-actions/julia-processcoverage@v1
+        with:
+          directories: lib/NonlinearSolveBase/src,lib/NonlinearSolveBase/ext,lib/NonlinearSolveHomotopyContinuation/src
+      - uses: codecov/codecov-action@v5
+        with:
+          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          verbose: true
+          fail_ci_if_error: false

From a6f6e7bbbe6a0fa977fa06435dfdb91f4c10e93b Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sun, 26 Jan 2025 16:53:15 +0530
Subject: [PATCH 16/24] docs: add NonlinearSolveHomotopyContinuation to docs

---
 .github/workflows/Documentation.yml  |  2 +-
 docs/Project.toml                    |  1 +
 docs/pages.jl                        |  3 ++-
 docs/src/api/homotopycontinuation.md | 18 ++++++++++++++++++
 4 files changed, 22 insertions(+), 2 deletions(-)
 create mode 100644 docs/src/api/homotopycontinuation.md

diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index 2c81172a3..63942ca52 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -22,7 +22,7 @@ jobs:
           Pkg.Registry.update()
           # Install packages present in subdirectories
           dev_pks = Pkg.PackageSpec[]
-          for path in ("lib/SciMLJacobianOperators", ".", "lib/SimpleNonlinearSolve", "lib/NonlinearSolveBase", "lib/BracketingNonlinearSolve", "lib/NonlinearSolveFirstOrder", "lib/NonlinearSolveQuasiNewton", "lib/NonlinearSolveSpectralMethods")
+          for path in ("lib/SciMLJacobianOperators", ".", "lib/SimpleNonlinearSolve", "lib/NonlinearSolveBase", "lib/BracketingNonlinearSolve", "lib/NonlinearSolveFirstOrder", "lib/NonlinearSolveQuasiNewton", "lib/NonlinearSolveSpectralMethods", "lib/NonlinearSolveHomotopyContinuation")
               push!(dev_pks, Pkg.PackageSpec(; path))
           end
           Pkg.develop(dev_pks)
diff --git a/docs/Project.toml b/docs/Project.toml
index c6c28820c..501497ceb 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -15,6 +15,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
+NonlinearSolveHomotopyContinuation = "2ac3b008-d579-4536-8c91-a1a5998c2f8b"
 NonlinearSolveQuasiNewton = "9a2c21bd-3a47-402d-9113-8faf9a0ee114"
 NonlinearSolveSpectralMethods = "26075421-4e9a-44e1-8bd1-420ed7ad02b2"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
diff --git a/docs/pages.jl b/docs/pages.jl
index 2834615d0..80dd4290a 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -49,7 +49,8 @@ pages = [
         "api/petsc.md",
         "api/siamfanlequations.md",
         "api/speedmapping.md",
-        "api/sundials.md"
+        "api/sundials.md",
+        "api/homotopycontinuation.md"
     ],
     "Sub-Packages" => Any[
         "api/SciMLJacobianOperators.md",
diff --git a/docs/src/api/homotopycontinuation.md b/docs/src/api/homotopycontinuation.md
new file mode 100644
index 000000000..11f7a6aee
--- /dev/null
+++ b/docs/src/api/homotopycontinuation.md
@@ -0,0 +1,18 @@
+# HomotopyContinuation.jl
+
+NonlinearSolve wraps the homotopy continuation algorithm implemented in
+HomotopyContinuation.jl. This solver is not included by default and needs
+to be installed separately:
+
+```julia
+using Pkg
+Pkg.add("NonlinearSolveHomotopyContinuation")
+using NonlinearSolveHomotopyContinuation, NonlinearSolve
+```
+
+# Solver API
+
+```@docs
+HomotopyContinuationJL
+SciMLBase.HomotopyContinuationFunction
+```

From 9c9bc79e46590f8ea8fa6e3fb8fe2eb96a9bce10 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sun, 26 Jan 2025 17:39:37 +0530
Subject: [PATCH 17/24] build: fix LinearAlgebra compat

---
 lib/NonlinearSolveHomotopyContinuation/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/Project.toml b/lib/NonlinearSolveHomotopyContinuation/Project.toml
index 517d7b27d..a9d8e5734 100644
--- a/lib/NonlinearSolveHomotopyContinuation/Project.toml
+++ b/lib/NonlinearSolveHomotopyContinuation/Project.toml
@@ -24,7 +24,7 @@ ConcreteStructs = "0.2.3"
 DifferentiationInterface = "0.6.27"
 DocStringExtensions = "0.9.3"
 HomotopyContinuation = "2.12.0"
-LinearAlgebra = "1.11.0"
+LinearAlgebra = "1.10"
 NonlinearSolve = "4"
 NonlinearSolveBase = "1.3.3"
 SciMLBase = "2.71"

From d9b337870954350246f0279c68ef2b65c111596c Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 27 Jan 2025 11:39:34 +0530
Subject: [PATCH 18/24] test: do not rely on singular roots for tests

---
 .../test/allroots.jl                            | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
index 18a768574..854f8d1bf 100644
--- a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
@@ -6,23 +6,24 @@ alg = HomotopyContinuationJL{true}(; threading = false)
 
 @testset "scalar u" begin
     rhs = function (u, p)
-        return u * u - 2p * u + p * p
+        return u * u - p[1] * u + p[2]
     end
     jac = function (u, p)
-        return 2u - 2p
+        return 2u - p[1]
     end
     @testset "`NonlinearProblem` - $name" for (jac, name) in [(nothing, "no jac"), (jac, "jac")]
         fn = NonlinearFunction(rhs; jac)
-        prob = NonlinearProblem(fn, 1.0, 2.0)
+        prob = NonlinearProblem(fn, 1.0, [5.0, 6.0])
         sol = solve(prob, alg)
 
         @test sol isa EnsembleSolution
         @test sol.converged
-        for nlsol in sol
-            @test nlsol isa NonlinearSolution
-            @test SciMLBase.successful_retcode(nlsol)
-            @test nlsol.u[1] ≈ 2.0 atol = 1e-7
-        end
+        @test sol.u[1] isa NonlinearSolution
+        @test SciMLBase.successful_retcode(sol.u[1])
+        @test sol.u[1].u ≈ 2.0 atol = 1e-10
+        @test sol.u[2] isa NonlinearSolution
+        @test SciMLBase.successful_retcode(sol.u[2])
+        @test sol.u[2].u ≈ 3.0 atol = 1e-10
 
         @testset "no real solutions" begin
             prob = NonlinearProblem(1.0, 0.5) do u, p

From f66af4e121d04318a7a7c9e119fdcfa1a8cbe178 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 27 Jan 2025 12:52:42 +0530
Subject: [PATCH 19/24] test: sort roots in `allroots` test

---
 lib/NonlinearSolveHomotopyContinuation/test/allroots.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
index 854f8d1bf..caef4d932 100644
--- a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
@@ -18,6 +18,7 @@ alg = HomotopyContinuationJL{true}(; threading = false)
 
         @test sol isa EnsembleSolution
         @test sol.converged
+        sort!(sol.u; by = x -> x.u)
         @test sol.u[1] isa NonlinearSolution
         @test SciMLBase.successful_retcode(sol.u[1])
         @test sol.u[1].u ≈ 2.0 atol = 1e-10

From ac6c3dc494cfa3191d65245ff7f10f9f6b7ca349 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 27 Jan 2025 13:14:44 +0530
Subject: [PATCH 20/24] refactor: format

---
 .../src/NonlinearSolveHomotopyContinuation.jl |  3 +-
 .../src/interface_types.jl                    | 49 +++++++++++++------
 .../src/solve.jl                              | 22 ++++++---
 .../test/allroots.jl                          | 34 +++++++------
 .../test/single_root.jl                       | 34 +++++++------
 5 files changed, 88 insertions(+), 54 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
index 4513d3349..0dbe08edf 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
@@ -45,7 +45,8 @@ specified using ADTypes.jl.
 HomotopyContinuation.jl requires the taylor series of the polynomial system for the single
 root method. This is automatically computed using TaylorSeries.jl.
 """
-@concrete struct HomotopyContinuationJL{AllRoots} <: NonlinearSolveBase.AbstractNonlinearSolveAlgorithm
+@concrete struct HomotopyContinuationJL{AllRoots} <:
+                 NonlinearSolveBase.AbstractNonlinearSolveAlgorithm
     autodiff
     kwargs
 end
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
index 24f23baea..6108d7676 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
@@ -18,7 +18,8 @@ regardless of the signature of ``f``.
     f
 end
 
-function (cjw::ComplexJacobianWrapper{Inplace})(u::AbstractVector{T}, x::AbstractVector{T}, p) where {T}
+function (cjw::ComplexJacobianWrapper{Inplace})(
+        u::AbstractVector{T}, x::AbstractVector{T}, p) where {T}
     x = reinterpret(Complex{T}, x)
     u = reinterpret(Complex{T}, u)
     cjw.f(u, x, p)
@@ -26,7 +27,8 @@ function (cjw::ComplexJacobianWrapper{Inplace})(u::AbstractVector{T}, x::Abstrac
     return u
 end
 
-function (cjw::ComplexJacobianWrapper{OutOfPlace})(u::AbstractVector{T}, x::AbstractVector{T}, p) where {T}
+function (cjw::ComplexJacobianWrapper{OutOfPlace})(
+        u::AbstractVector{T}, x::AbstractVector{T}, p) where {T}
     x = reinterpret(Complex{T}, x)
     u_tmp = cjw.f(x, p)
     u_tmp = reinterpret(T, u_tmp)
@@ -34,7 +36,8 @@ function (cjw::ComplexJacobianWrapper{OutOfPlace})(u::AbstractVector{T}, x::Abst
     return u
 end
 
-function (cjw::ComplexJacobianWrapper{Scalar})(u::AbstractVector{T}, x::AbstractVector{T}, p) where {T}
+function (cjw::ComplexJacobianWrapper{Scalar})(
+        u::AbstractVector{T}, x::AbstractVector{T}, p) where {T}
     x = reinterpret(Complex{T}, x)
     u_tmp = cjw.f(x[1], p)
     u[1] = real(u_tmp)
@@ -52,7 +55,8 @@ polynomial systems specified using `NonlinearProblem`.
 
 $(FIELDS)
 """
-@concrete struct HomotopySystemWrapper{variant <: HomotopySystemVariant} <: HC.AbstractSystem
+@concrete struct HomotopySystemWrapper{variant <: HomotopySystemVariant} <:
+                 HC.AbstractSystem
     """
     The wrapped polynomial function.
     """
@@ -107,14 +111,16 @@ function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{Scalar}, x, p = not
     return u
 end
 
-function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Inplace}, x, p = nothing)
+function HC.ModelKit.evaluate_and_jacobian!(
+        u, U, sys::HomotopySystemWrapper{Inplace}, x, p = nothing)
     p = sys.p
     sys.f(u, x, p)
     sys.jac(U, x, p)
     return u, U
 end
 
-function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing)
+function HC.ModelKit.evaluate_and_jacobian!(
+        u, U, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing)
     p = sys.p
     u_tmp = sys.f(x, p)
     copyto!(u, u_tmp)
@@ -123,7 +129,8 @@ function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Out
     return u, U
 end
 
-function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Scalar}, x, p = nothing)
+function HC.ModelKit.evaluate_and_jacobian!(
+        u, U, sys::HomotopySystemWrapper{Scalar}, x, p = nothing)
     p = sys.p
     u[1] = sys.f(x[1], p)
     U[1] = sys.jac(x[1], p)
@@ -131,7 +138,9 @@ function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{Sca
 end
 
 for V in (Inplace, OutOfPlace, Scalar)
-    @eval function HC.ModelKit.evaluate_and_jacobian!(u, U, sys::HomotopySystemWrapper{$V, F, J}, x, p = nothing) where {F, J <: ComplexJacobianWrapper}
+    @eval function HC.ModelKit.evaluate_and_jacobian!(
+            u, U, sys::HomotopySystemWrapper{$V, F, J}, x,
+            p = nothing) where {F, J <: ComplexJacobianWrapper}
         p = sys.p
         U_tmp = sys.jacobian_buffers
         x = reinterpret(Float64, x)
@@ -151,9 +160,10 @@ for V in (Inplace, OutOfPlace, Scalar)
     end
 end
 
-function update_taylorvars_from_taylorvector!(vars, x::HC.ModelKit.TaylorVector{M}) where {M}
+function update_taylorvars_from_taylorvector!(
+        vars, x::HC.ModelKit.TaylorVector{M}) where {M}
     for i in eachindex(vars)
-        for j in 0:M-1
+        for j in 0:(M - 1)
             vars[i][j] = x[i, j + 1]
         end
         for j in M:4
@@ -171,19 +181,22 @@ function update_taylorvars_from_taylorvector!(vars, x::AbstractVector)
     end
 end
 
-function update_maybe_taylorvector_from_taylorvars!(u::Vector, vars, buffer, ::Val{N}) where {N}
+function update_maybe_taylorvector_from_taylorvars!(
+        u::Vector, vars, buffer, ::Val{N}) where {N}
     for i in eachindex(vars)
         u[i] = buffer[i][N]
     end
 end
 
-function update_maybe_taylorvector_from_taylorvars!(u::HC.ModelKit.TaylorVector, vars, buffer, ::Val{N}) where {N}
+function update_maybe_taylorvector_from_taylorvars!(
+        u::HC.ModelKit.TaylorVector, vars, buffer, ::Val{N}) where {N}
     for i in eachindex(vars)
         u[i] = ntuple(j -> buffer[i][j - 1], Val(N + 1))
     end
 end
 
-function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{Inplace}, x, p = nothing) where {N}
+function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N},
+        sys::HomotopySystemWrapper{Inplace}, x, p = nothing) where {N}
     f = sys.f
     p = sys.p
     buffer, vars = sys.taylorvars
@@ -193,7 +206,8 @@ function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWra
     return u
 end
 
-function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing) where {N}
+function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N},
+        sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing) where {N}
     f = sys.f
     p = sys.p
     vars = sys.taylorvars
@@ -203,7 +217,8 @@ function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWra
     return u
 end
 
-function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, sys::HomotopySystemWrapper{Scalar}, x, p = nothing) where {N}
+function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N},
+        sys::HomotopySystemWrapper{Scalar}, x, p = nothing) where {N}
     f = sys.f
     p = sys.p
     var = sys.taylorvars
@@ -270,8 +285,10 @@ function HC.ModelKit.evaluate_and_jacobian!(u, U, h::GuessHomotopy, x, t, p = no
     return u, U
 end
 
-HC.ModelKit.taylor!(u, v::Val{N}, H::GuessHomotopy, tx, t, incremental::Bool) where {N} =
+function HC.ModelKit.taylor!(
+        u, v::Val{N}, H::GuessHomotopy, tx, t, incremental::Bool) where {N}
     HC.ModelKit.taylor!(u, v, H, tx, t)
+end
 
 function HC.ModelKit.taylor!(u, ::Val{N}, h::GuessHomotopy, x, t, p = nothing) where {N}
     HC.ModelKit.taylor!(h.taylorbuffer, Val(N), h.sys, x, p)
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
index 913f66b44..f14020152 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
@@ -4,7 +4,8 @@
 Create and return the appropriate `HomotopySystemWrapper` to use for solving the given
 `prob` with `alg`.
 """
-function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::HomotopyContinuationJL)
+function homotopy_continuation_preprocessing(
+        prob::NonlinearProblem, alg::HomotopyContinuationJL)
     # cast to a `HomotopyNonlinearFunction`
     f = if prob.f isa HomotopyNonlinearFunction
         prob.f
@@ -45,7 +46,8 @@ function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::Homoto
     taylorvars = if isscalar
         Taylor1(zeros(ComplexF64, 5), 4)
     elseif iip
-        ([Taylor1(zeros(ComplexF64, 5), 4) for _ in u0], [Taylor1(zeros(ComplexF64, 5), 4) for _ in u0])
+        ([Taylor1(zeros(ComplexF64, 5), 4) for _ in u0],
+            [Taylor1(zeros(ComplexF64, 5), 4) for _ in u0])
     else
         [Taylor1(zeros(ComplexF64, 5), 4) for _ in u0]
     end
@@ -57,12 +59,14 @@ function homotopy_continuation_preprocessing(prob::NonlinearProblem, alg::Homoto
     end
 
     # HC-compatible system
-    hcsys = HomotopySystemWrapper{variant}(f.f.f, jac, p, alg.autodiff, prep, vars, taylorvars, jacobian_buffers)
+    hcsys = HomotopySystemWrapper{variant}(
+        f.f.f, jac, p, alg.autodiff, prep, vars, taylorvars, jacobian_buffers)
 
     return f, hcsys
 end
 
-function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{true}; denominator_abstol = 1e-7, kwargs...)
+function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{true};
+        denominator_abstol = 1e-7, kwargs...)
     f, hcsys = homotopy_continuation_preprocessing(prob, alg)
 
     u0 = state_values(prob)
@@ -110,7 +114,8 @@ function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{t
     return SciMLBase.EnsembleSolution(nlsols, 0.0, true, nothing)
 end
 
-function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{false}; denominator_abstol = 1e-7, kwargs...)
+function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{false};
+        denominator_abstol = 1e-7, kwargs...)
     f, hcsys = homotopy_continuation_preprocessing(prob, alg)
 
     u0 = state_values(prob)
@@ -120,14 +125,16 @@ function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{f
     fu0 = NonlinearSolveBase.Utils.evaluate_f(prob, u0_p)
 
     homotopy = GuessHomotopy(hcsys, fu0)
-    orig_sol = HC.solve(homotopy, u0_p isa Number ? [[u0_p]] : [u0_p]; alg.kwargs..., kwargs...)
+    orig_sol = HC.solve(
+        homotopy, u0_p isa Number ? [[u0_p]] : [u0_p]; alg.kwargs..., kwargs...)
     realsols = map(res -> res.solution, HC.results(orig_sol; only_real = true))
     if u0 isa Number
         realsols = map(only, realsols)
     end
 
     # no real solutions or infeasible solution
-    if isempty(realsols) || any(<=(denominator_abstol), map(abs, f.denominator(real.(only(realsols)), p)))
+    if isempty(realsols) ||
+       any(<=(denominator_abstol), map(abs, f.denominator(real.(only(realsols)), p)))
         retcode = if isempty(realsols)
             SciMLBase.ReturnCode.ConvergenceFailure
         else
@@ -153,4 +160,3 @@ function CommonSolve.solve(prob::NonlinearProblem, alg::HomotopyContinuationJL{f
     retcode = SciMLBase.ReturnCode.Success
     return SciMLBase.build_solution(prob, alg, u, resid; retcode, original = orig_sol)
 end
-
diff --git a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
index caef4d932..a715e0189 100644
--- a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
@@ -11,7 +11,8 @@ alg = HomotopyContinuationJL{true}(; threading = false)
     jac = function (u, p)
         return 2u - p[1]
     end
-    @testset "`NonlinearProblem` - $name" for (jac, name) in [(nothing, "no jac"), (jac, "jac")]
+    @testset "`NonlinearProblem` - $name" for (jac, name) in [
+        (nothing, "no jac"), (jac, "jac")]
         fn = NonlinearFunction(rhs; jac)
         prob = NonlinearProblem(fn, 1.0, [5.0, 6.0])
         sol = solve(prob, alg)
@@ -21,14 +22,14 @@ alg = HomotopyContinuationJL{true}(; threading = false)
         sort!(sol.u; by = x -> x.u)
         @test sol.u[1] isa NonlinearSolution
         @test SciMLBase.successful_retcode(sol.u[1])
-        @test sol.u[1].u ≈ 2.0 atol = 1e-10
+        @test sol.u[1].u≈2.0 atol=1e-10
         @test sol.u[2] isa NonlinearSolution
         @test SciMLBase.successful_retcode(sol.u[2])
-        @test sol.u[2].u ≈ 3.0 atol = 1e-10
+        @test sol.u[2].u≈3.0 atol=1e-10
 
         @testset "no real solutions" begin
             prob = NonlinearProblem(1.0, 0.5) do u, p
-                return u * u - 2p * u  + p
+                return u * u - 2p * u + p
             end
             sol = solve(prob, alg)
             @test length(sol) == 1
@@ -47,7 +48,8 @@ alg = HomotopyContinuationJL{true}(; threading = false)
         unpolynomialize = function (u, p)
             return [asin(u)]
         end
-        fn = HomotopyNonlinearFunction(; denominator, polynomialize, unpolynomialize) do u, p
+        fn = HomotopyNonlinearFunction(;
+            denominator, polynomialize, unpolynomialize) do u, p
             return (u - p[1]) * (u - p[2])
         end
         prob = NonlinearProblem(fn, 0.0, [0.5, 0.7])
@@ -74,12 +76,12 @@ alg = HomotopyContinuationJL{true}(; threading = false)
 end
 
 f! = function (du, u, p)
-    du[1] = u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1
-    du[2] = u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]
+    du[1] = u[1] * u[1] - p[1] * u[2] + u[2]^3 + 1
+    du[2] = u[2]^3 + 2 * p[2] * u[1] * u[2] + u[2]
 end
 
 f = function (u, p)
-    [u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1, u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]]
+    [u[1] * u[1] - p[1] * u[2] + u[2]^3 + 1, u[2]^3 + 2 * p[2] * u[1] * u[2] + u[2]]
 end
 
 jac! = function (j, u, p)
@@ -89,11 +91,13 @@ jac! = function (j, u, p)
     j[2, 2] = 3 * u[2]^2 + 2 * p[2] * u[1] + 1
 end
 jac = function (u, p)
-    [2u[1] -p[1] + 3 * u[2]^2;
-     2 * p[2] * u[2] 3 * u[2]^2 + 2 * p[2] * u[1] + 1]
+    [2u[1] -p[1]+3 * u[2]^2;
+     2*p[2]*u[2] 3*u[2]^2+2*p[2]*u[1]+1]
 end
 
-@testset "vector u - $name" for (rhs, jac, name) in [(f, nothing, "oop"), (f, jac, "oop + jac"), (f!, nothing, "iip"), (f!, jac!, "iip + jac")]
+@testset "vector u - $name" for (rhs, jac, name) in [
+    (f, nothing, "oop"), (f, jac, "oop + jac"),
+    (f!, nothing, "iip"), (f!, jac!, "iip + jac")]
     sol = nothing
     @testset "`NonlinearProblem`" begin
         fn = NonlinearFunction(rhs; jac)
@@ -103,7 +107,7 @@ end
         @test sol.converged
         for nlsol in sol.u
             @test SciMLBase.successful_retcode(nlsol)
-            @test f(nlsol.u, prob.p) ≈ [0.0, 0.0] atol = 1e-10
+            @test f(nlsol.u, prob.p)≈[0.0, 0.0] atol=1e-10
         end
 
         @testset "no real solutions" begin
@@ -123,7 +127,7 @@ end
             return [[cbrt(u[1]), sin(u[2] / 40)]]
         end
         polynomialize = function (u, p)
-            return [u[1] ^ 3, 40asin(u[2])]
+            return [u[1]^3, 40asin(u[2])]
         end
         nlfn = NonlinearFunction(rhs; jac)
         fn = HomotopyNonlinearFunction(nlfn; denominator, polynomialize, unpolynomialize)
@@ -133,7 +137,9 @@ end
         @test sol2.converged
         @test length(sol.u) == length(sol2.u)
         for nlsol2 in sol2.u
-            @test any(nlsol -> isapprox(polynomialize(nlsol2.u, prob.p), nlsol.u; rtol = 1e-8), sol.u)
+            @test any(
+                nlsol -> isapprox(polynomialize(nlsol2.u, prob.p), nlsol.u; rtol = 1e-8),
+                sol.u)
         end
 
         @testset "some invalid solutions" begin
diff --git a/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl b/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl
index 3cd5b81cf..9e41ee6fa 100644
--- a/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl
@@ -11,13 +11,14 @@ alg = HomotopyContinuationJL{false}(; threading = false)
     jac = function (u, p)
         return 2u - (p + 3)
     end
-    @testset "`NonlinearProblem` - $name" for (jac, name) in [(nothing, "no jac"), (jac, "jac")]
+    @testset "`NonlinearProblem` - $name" for (jac, name) in [
+        (nothing, "no jac"), (jac, "jac")]
         fn = NonlinearFunction(rhs; jac)
         prob = NonlinearProblem(fn, 1.0, 2.0)
         sol = solve(prob, alg)
 
         @test sol isa NonlinearSolution
-        @test sol.u ≈ 2.0 atol = 1e-10
+        @test sol.u≈2.0 atol=1e-10
 
         @testset "no real solutions" begin
             prob = NonlinearProblem(1.0, 1.0) do u, p
@@ -38,13 +39,14 @@ alg = HomotopyContinuationJL{false}(; threading = false)
         unpolynomialize = function (u, p)
             return [asin(u)]
         end
-        fn = HomotopyNonlinearFunction(; denominator, polynomialize, unpolynomialize) do u, p
+        fn = HomotopyNonlinearFunction(;
+            denominator, polynomialize, unpolynomialize) do u, p
             return (u - p[1]) * (u - p[2])
         end
         prob = NonlinearProblem(fn, 0.0, [0.5, 0.7])
 
         sol = solve(prob, alg)
-        @test sin(sol.u[1]) ≈ 0.5 atol=1e-10
+        @test sin(sol.u[1])≈0.5 atol=1e-10
 
         @testset "no valid solutions" begin
             prob2 = remake(prob; p = [0.7, 0.7])
@@ -55,21 +57,21 @@ alg = HomotopyContinuationJL{false}(; threading = false)
         @testset "closest root" begin
             prob3 = remake(prob; p = [0.5, 0.6], u0 = asin(0.4))
             sol3 = solve(prob3, alg)
-            @test sin(sol3.u) ≈ 0.5 atol = 1e-10
+            @test sin(sol3.u)≈0.5 atol=1e-10
             prob4 = remake(prob3; u0 = asin(0.7))
             sol4 = solve(prob4, alg)
-            @test sin(sol4.u) ≈ 0.6 atol = 1e-10
+            @test sin(sol4.u)≈0.6 atol=1e-10
         end
     end
 end
 
 f! = function (du, u, p)
-    du[1] = u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1
-    du[2] = u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]
+    du[1] = u[1] * u[1] - p[1] * u[2] + u[2]^3 + 1
+    du[2] = u[2]^3 + 2 * p[2] * u[1] * u[2] + u[2]
 end
 
 f = function (u, p)
-    [u[1] * u[1] - p[1] * u[2] + u[2] ^ 3 + 1, u[2] ^ 3 + 2 * p[2] * u[1] * u[2] + u[2]]
+    [u[1] * u[1] - p[1] * u[2] + u[2]^3 + 1, u[2]^3 + 2 * p[2] * u[1] * u[2] + u[2]]
 end
 
 jac! = function (j, u, p)
@@ -80,17 +82,19 @@ jac! = function (j, u, p)
 end
 
 jac = function (u, p)
-    [2u[1] -p[1] + 3 * u[2]^2;
-     2 * p[2] * u[2] 3 * u[2]^2 + 2 * p[2] * u[1] + 1]
+    [2u[1] -p[1]+3 * u[2]^2;
+     2*p[2]*u[2] 3*u[2]^2+2*p[2]*u[1]+1]
 end
 
-@testset "vector u - $name" for (rhs, jac, name) in [(f, nothing, "oop"), (f, jac, "oop + jac"), (f!, nothing, "iip"), (f!, jac!, "iip + jac")]
+@testset "vector u - $name" for (rhs, jac, name) in [
+    (f, nothing, "oop"), (f, jac, "oop + jac"),
+    (f!, nothing, "iip"), (f!, jac!, "iip + jac")]
     @testset "`NonlinearProblem`" begin
         fn = NonlinearFunction(rhs; jac)
         prob = NonlinearProblem(fn, [1.0, 2.0], [2.0, 3.0])
         sol = solve(prob, alg)
         @test SciMLBase.successful_retcode(sol)
-        @test f(sol.u, prob.p) ≈ [0.0, 0.0] atol = 1e-10
+        @test f(sol.u, prob.p)≈[0.0, 0.0] atol=1e-10
 
         @testset "no real solutions" begin
             prob2 = remake(prob; p = zeros(2))
@@ -107,14 +111,14 @@ end
             return [[cbrt(u[1]), sin(u[2] / 40)]]
         end
         polynomialize = function (u, p)
-            return [u[1] ^ 3, 40asin(u[2])]
+            return [u[1]^3, 40asin(u[2])]
         end
         nlfn = NonlinearFunction(rhs; jac)
         fn = HomotopyNonlinearFunction(nlfn; denominator, polynomialize, unpolynomialize)
         prob = NonlinearProblem(fn, [1.0, 1.0], [2.0, 3.0, 4.0, 5.0])
         sol = solve(prob, alg)
         @test SciMLBase.successful_retcode(sol)
-        @test f(polynomialize(sol.u, prob.p), prob.p) ≈ zeros(2) atol = 1e-10
+        @test f(polynomialize(sol.u, prob.p), prob.p)≈zeros(2) atol=1e-10
 
         @testset "some invalid solutions" begin
             prob2 = remake(prob; p = [2.0, 3.0, polynomialize(sol.u, prob.p)...])

From 8fa843c6808fe0c083725d4132257b7e3f859394 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 27 Jan 2025 15:23:45 +0530
Subject: [PATCH 21/24] test: do not rely on singular roots

---
 lib/NonlinearSolveHomotopyContinuation/test/allroots.jl    | 4 ++--
 lib/NonlinearSolveHomotopyContinuation/test/single_root.jl | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
index a715e0189..c0d025e49 100644
--- a/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl
@@ -40,7 +40,7 @@ alg = HomotopyContinuationJL{true}(; threading = false)
 
     @testset "`HomotopyContinuationFunction`" begin
         denominator = function (u, p)
-            return [u - 0.7]
+            return [u - 0.7, u - 0.9]
         end
         polynomialize = function (u, p)
             return sin(u)
@@ -59,7 +59,7 @@ alg = HomotopyContinuationJL{true}(; threading = false)
         @test sin(sol.u[1][1]) ≈ 0.5
 
         @testset "no valid solutions" begin
-            prob2 = remake(prob; p = [0.7, 0.7])
+            prob2 = remake(prob; p = [0.7, 0.9])
             sol2 = solve(prob2, alg)
             @test !sol2.converged
             @test length(sol2) == 1
diff --git a/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl b/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl
index 9e41ee6fa..049578b37 100644
--- a/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl
@@ -31,7 +31,7 @@ alg = HomotopyContinuationJL{false}(; threading = false)
 
     @testset "`HomotopyContinuationFunction`" begin
         denominator = function (u, p)
-            return [u - 0.7]
+            return [u - 0.7, u - 0.9]
         end
         polynomialize = function (u, p)
             return sin(u)
@@ -49,7 +49,7 @@ alg = HomotopyContinuationJL{false}(; threading = false)
         @test sin(sol.u[1])≈0.5 atol=1e-10
 
         @testset "no valid solutions" begin
-            prob2 = remake(prob; p = [0.7, 0.7])
+            prob2 = remake(prob; p = [0.7, 0.9])
             sol2 = solve(prob2, alg)
             @test sol2.retcode == SciMLBase.ReturnCode.Infeasible
         end

From 1decea630d7cdef80a493f53c597bfda137b4658 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Tue, 28 Jan 2025 09:40:42 -0500
Subject: [PATCH 22/24] Update homotopycontinuation.md

---
 docs/src/api/homotopycontinuation.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/api/homotopycontinuation.md b/docs/src/api/homotopycontinuation.md
index 11f7a6aee..1e46c6378 100644
--- a/docs/src/api/homotopycontinuation.md
+++ b/docs/src/api/homotopycontinuation.md
@@ -13,6 +13,6 @@ using NonlinearSolveHomotopyContinuation, NonlinearSolve
 # Solver API
 
 ```@docs
-HomotopyContinuationJL
+NonlinearSolveHomotopyContinuation.HomotopyContinuationJL
 SciMLBase.HomotopyContinuationFunction
 ```

From 5426e992af2ca37ebcb9dec7d3d6f77354fae36e Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 29 Jan 2025 12:37:10 +0530
Subject: [PATCH 23/24] docs: add `NonlinearSolveHomotopyContinuation` to doc
 build

---
 docs/make.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/make.jl b/docs/make.jl
index 29b1535dd..06a481285 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -5,6 +5,7 @@ using Sundials
 using NonlinearSolveBase, SciMLBase, DiffEqBase
 using SimpleNonlinearSolve, BracketingNonlinearSolve
 using NonlinearSolveFirstOrder, NonlinearSolveQuasiNewton, NonlinearSolveSpectralMethods
+using NonlinearSolveHomotopyContinuation
 using SciMLJacobianOperators
 using NonlinearSolve, SteadyStateDiffEq
 
@@ -35,6 +36,7 @@ makedocs(;
         NonlinearSolveBase, SciMLBase, DiffEqBase,
         SimpleNonlinearSolve, BracketingNonlinearSolve,
         NonlinearSolveFirstOrder, NonlinearSolveQuasiNewton, NonlinearSolveSpectralMethods,
+        NonlinearSolveHomotopyContinuation,
         Sundials,
         SciMLJacobianOperators,
         NonlinearSolve, SteadyStateDiffEq

From d2e8ec6eb29e5a2e0ebdc758ef20f839bd715722 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 29 Jan 2025 18:32:56 +0530
Subject: [PATCH 24/24] refactor: use TaylorDiff.jl instead of TaylorSeries.jl

---
 .../Project.toml                              |  4 +-
 .../src/NonlinearSolveHomotopyContinuation.jl |  2 +-
 .../src/interface_types.jl                    | 85 ++++++++++++++-----
 .../src/solve.jl                              | 10 ++-
 4 files changed, 71 insertions(+), 30 deletions(-)

diff --git a/lib/NonlinearSolveHomotopyContinuation/Project.toml b/lib/NonlinearSolveHomotopyContinuation/Project.toml
index a9d8e5734..4b4d3a512 100644
--- a/lib/NonlinearSolveHomotopyContinuation/Project.toml
+++ b/lib/NonlinearSolveHomotopyContinuation/Project.toml
@@ -14,7 +14,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
+TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c"
 
 [compat]
 ADTypes = "1.11.0"
@@ -29,7 +29,7 @@ NonlinearSolve = "4"
 NonlinearSolveBase = "1.3.3"
 SciMLBase = "2.71"
 SymbolicIndexingInterface = "0.3.36"
-TaylorSeries = "0.18.2"
+TaylorDiff = "0.3.1"
 Test = "1.10"
 julia = "1.10"
 
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
index 0dbe08edf..5f68e55a8 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl
@@ -6,7 +6,7 @@ using NonlinearSolveBase
 using SymbolicIndexingInterface
 using LinearAlgebra
 using ADTypes
-using TaylorSeries
+using TaylorDiff
 using DocStringExtensions
 import CommonSolve
 import HomotopyContinuation as HC
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
index 6108d7676..63a64ec36 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl
@@ -83,7 +83,7 @@ $(FIELDS)
     """
     vars
     """
-    The `TaylorSeries.Taylor1` objects used to compute the taylor series of `f`.
+    The `TaylorDiff.TaylorScalar` objects used to compute the taylor series of `f`.
     """
     taylorvars
     """
@@ -161,37 +161,58 @@ for V in (Inplace, OutOfPlace, Scalar)
 end
 
 function update_taylorvars_from_taylorvector!(
-        vars, x::HC.ModelKit.TaylorVector{M}) where {M}
-    for i in eachindex(vars)
-        for j in 0:(M - 1)
-            vars[i][j] = x[i, j + 1]
+        vars, x::HC.ModelKit.TaylorVector)
+    for i in eachindex(x)
+        xvar = x[i]
+        realx = ntuple(Val(4)) do j
+            j <= length(xvar) ? real(xvar[j - 1]) : 0.0
         end
-        for j in M:4
-            vars[i][j] = zero(vars[i][j])
+        imagx = ntuple(Val(4)) do j
+            j <= length(xvar) ? imag(xvar[j - 1]) : 0.0
         end
+
+        vars[2i - 1] = TaylorScalar(realx)
+        vars[2i] = TaylorScalar(imagx)
     end
 end
 
 function update_taylorvars_from_taylorvector!(vars, x::AbstractVector)
-    for i in eachindex(vars)
-        vars[i][0] = x[i]
-        for j in 1:4
-            vars[i][j] = zero(vars[i][j])
-        end
+    for i in eachindex(x)
+        vars[2i - 1] = TaylorScalar(real(x[i]), ntuple(Returns(0.0), Val(3)))
+        vars[2i] = TaylorScalar(imag(x[i]), ntuple(Returns(0.0), Val(3)))
+    end
+end
+
+function check_taylor_equality(vars, x::HC.ModelKit.TaylorVector)
+    for i in eachindex(x)
+        TaylorDiff.flatten(vars[2i-1]) == map(real, x[i]) || return false
+        TaylorDiff.flatten(vars[2i]) == map(imag, x[i]) || return false
+    end
+    return true
+end
+function check_taylor_equality(vars, x::AbstractVector)
+    for i in eachindex(x)
+        TaylorDiff.value(vars[2i-1]) != real(x[i]) && return false
+        TaylorDiff.value(vars[2i]) != imag(x[i]) && return false
     end
+    return true
 end
 
 function update_maybe_taylorvector_from_taylorvars!(
         u::Vector, vars, buffer, ::Val{N}) where {N}
     for i in eachindex(vars)
-        u[i] = buffer[i][N]
+        rval = TaylorDiff.flatten(real(buffer[i]))
+        ival = TaylorDiff.flatten(imag(buffer[i]))
+        u[i] = rval[N] + im * ival[N]
     end
 end
 
 function update_maybe_taylorvector_from_taylorvars!(
-        u::HC.ModelKit.TaylorVector, vars, buffer, ::Val{N}) where {N}
+        u::HC.ModelKit.TaylorVector{M}, vars, buffer, ::Val{N}) where {M, N}
     for i in eachindex(vars)
-        u[i] = ntuple(j -> buffer[i][j - 1], Val(N + 1))
+        rval = TaylorDiff.flatten(real(buffer[i]))
+        ival = TaylorDiff.flatten(imag(buffer[i]))
+        u[i] = ntuple(i -> rval[i]  + im * ival[i], Val(length(rval)))
     end
 end
 
@@ -199,9 +220,15 @@ function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N},
         sys::HomotopySystemWrapper{Inplace}, x, p = nothing) where {N}
     f = sys.f
     p = sys.p
-    buffer, vars = sys.taylorvars
-    update_taylorvars_from_taylorvector!(vars, x)
-    f(buffer, vars, p)
+    vars, buffer = sys.taylorvars
+    if !check_taylor_equality(vars, x)
+        update_taylorvars_from_taylorvector!(vars, x)
+        vars = reinterpret(Complex{eltype(vars)}, vars)
+        buffer = reinterpret(Complex{eltype(buffer)}, buffer)
+        f(buffer, vars, p)
+    else
+        vars = reinterpret(Complex{eltype(vars)}, vars)
+    end
     update_maybe_taylorvector_from_taylorvars!(u, vars, buffer, Val(N))
     return u
 end
@@ -211,8 +238,14 @@ function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N},
     f = sys.f
     p = sys.p
     vars = sys.taylorvars
-    update_taylorvars_from_taylorvector!(vars, x)
-    buffer = f(vars, p)
+    if !check_taylor_equality(vars, x)
+        update_taylorvars_from_taylorvector!(vars, x)
+        vars = reinterpret(Complex{eltype(vars)}, vars)
+        buffer = f(vars, p)
+        copyto!(vars, buffer)
+    else
+        vars = buffer = reinterpret(Complex{eltype(vars)}, vars)
+    end
     update_maybe_taylorvector_from_taylorvars!(u, vars, buffer, Val(N))
     return u
 end
@@ -222,9 +255,15 @@ function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N},
     f = sys.f
     p = sys.p
     var = sys.taylorvars
-    update_taylorvars_from_taylorvector!((var,), x)
-    buffer = f(var, p)
-    update_maybe_taylorvector_from_taylorvars!(u, (var,), (buffer,), Val(N))
+    if !check_taylor_equality(var, x)
+        update_taylorvars_from_taylorvector!(var, x)
+        var = reinterpret(Complex{eltype(var)}, var)
+        buffer = f(var[1], p)
+        var[1] = buffer
+    else
+        var = buffer = reinterpret(Complex{eltype(var)}, var)
+    end
+    update_maybe_taylorvector_from_taylorvars!(u, var, buffer, Val(N))
     return u
 end
 
diff --git a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
index f14020152..2db789df9 100644
--- a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
+++ b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl
@@ -44,12 +44,14 @@ function homotopy_continuation_preprocessing(
     end
 
     taylorvars = if isscalar
-        Taylor1(zeros(ComplexF64, 5), 4)
+        [TaylorScalar(ntuple(Returns(0.0), 4)), TaylorScalar(ntuple(Returns(0.0), 4))]
     elseif iip
-        ([Taylor1(zeros(ComplexF64, 5), 4) for _ in u0],
-            [Taylor1(zeros(ComplexF64, 5), 4) for _ in u0])
+        (
+            [TaylorScalar(ntuple(Returns(0.0), 4)) for _ in 1:2length(u0)],
+            [TaylorScalar(ntuple(Returns(0.0), 4)) for _ in 1:2length(u0)]
+        )
     else
-        [Taylor1(zeros(ComplexF64, 5), 4) for _ in u0]
+        [TaylorScalar(ntuple(Returns(0.0), 4)) for _ in 1:2length(u0)]
     end
 
     jacobian_buffers = if isscalar