From 7bb635db146cf2562b117bf30e2966f39a072ebd Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 31 Jan 2025 07:04:48 +0530 Subject: [PATCH] feat: add `NonlinearSolveHomotopyContinuation.jl` (#523) * feat: add `NonlinearSolveHomotopyContinuation.jl` * feat: add `solve` implementations * test: add tests for all roots solve * fix: fix jacobian implementation, add taylor implementation * fix: fix system evaluate and jacobian * fix: fix inplace system taylor * fix: fix homotopy and implement taylor * fix: fix jacobian building for systems * fix: fix single root solve implementation * build: add TaylorSeries as a dependency * fix: fix jacobian and taylor implementations * test: add tests for single roots, fix allroots tests * build: add compat entries * docs: add docstrings to NonlinearSolveHomotopyContinuation.jl * ci: add NonlinearSolveHomotopyContinuation.jl to CI * docs: add NonlinearSolveHomotopyContinuation to docs * build: fix LinearAlgebra compat * test: do not rely on singular roots for tests * test: sort roots in `allroots` test * refactor: format * test: do not rely on singular roots * Update homotopycontinuation.md * docs: add `NonlinearSolveHomotopyContinuation` to doc build * refactor: use TaylorDiff.jl instead of TaylorSeries.jl --------- Co-authored-by: Christopher Rackauckas --- .../CI_NonlinearSolveHomotopyContinuation.yml | 108 ++++++ .github/workflows/Documentation.yml | 2 +- docs/Project.toml | 1 + docs/make.jl | 2 + docs/pages.jl | 3 +- docs/src/api/homotopycontinuation.md | 18 + .../LICENSE | 21 ++ .../Project.toml | 42 +++ .../src/NonlinearSolveHomotopyContinuation.jl | 66 ++++ .../src/interface_types.jl | 340 ++++++++++++++++++ .../src/solve.jl | 164 +++++++++ .../test/allroots.jl | 151 ++++++++ .../test/runtests.jl | 15 + .../test/single_root.jl | 129 +++++++ 14 files changed, 1060 insertions(+), 2 deletions(-) create mode 100644 .github/workflows/CI_NonlinearSolveHomotopyContinuation.yml create mode 100644 docs/src/api/homotopycontinuation.md 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/src/interface_types.jl create mode 100644 lib/NonlinearSolveHomotopyContinuation/src/solve.jl create mode 100644 lib/NonlinearSolveHomotopyContinuation/test/allroots.jl create mode 100644 lib/NonlinearSolveHomotopyContinuation/test/runtests.jl create mode 100644 lib/NonlinearSolveHomotopyContinuation/test/single_root.jl 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 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/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 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..1e46c6378 --- /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 +NonlinearSolveHomotopyContinuation.HomotopyContinuationJL +SciMLBase.HomotopyContinuationFunction +``` 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 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..4b4d3a512 --- /dev/null +++ b/lib/NonlinearSolveHomotopyContinuation/Project.toml @@ -0,0 +1,42 @@ +name = "NonlinearSolveHomotopyContinuation" +uuid = "2ac3b008-d579-4536-8c91-a1a5998c2f8b" +authors = ["Aayush Sabharwal 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" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +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" +TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c" + +[compat] +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.10" +NonlinearSolve = "4" +NonlinearSolveBase = "1.3.3" +SciMLBase = "2.71" +SymbolicIndexingInterface = "0.3.36" +TaylorDiff = "0.3.1" +Test = "1.10" +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", "NonlinearSolve"] diff --git a/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl new file mode 100644 index 000000000..5f68e55a8 --- /dev/null +++ b/lib/NonlinearSolveHomotopyContinuation/src/NonlinearSolveHomotopyContinuation.jl @@ -0,0 +1,66 @@ +module NonlinearSolveHomotopyContinuation + +using SciMLBase: AbstractNonlinearProblem +using SciMLBase +using NonlinearSolveBase +using SymbolicIndexingInterface +using LinearAlgebra +using ADTypes +using TaylorDiff +using DocStringExtensions +import CommonSolve +import HomotopyContinuation as HC +import DifferentiationInterface as DI + +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 +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...) + +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..63a64ec36 --- /dev/null +++ b/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl @@ -0,0 +1,340 @@ +abstract type HomotopySystemVariant end + +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 + +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 + +""" + $(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 `TaylorDiff.TaylorScalar` objects used to compute the taylor series of `f`. + """ + taylorvars + """ + Preallocated intermediate buffers used for calculating the jacobian. + """ + jacobian_buffers +end + +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.f(u, x, sys.p) + return u +end + +function HC.ModelKit.evaluate!(u, sys::HomotopySystemWrapper{OutOfPlace}, x, p = nothing) + 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.f(x[1], sys.p) + return u +end + +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) + p = sys.p + u_tmp = sys.f(x, p) + 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) + p = sys.p + u[1] = sys.f(x[1], p) + U[1] = sys.jac(x[1], p) + return u, U +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} + 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) + for i in eachindex(x) + xvar = x[i] + realx = ntuple(Val(4)) do j + j <= length(xvar) ? real(xvar[j - 1]) : 0.0 + end + 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(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) + 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{M}, vars, buffer, ::Val{N}) where {M, N} + for i in eachindex(vars) + 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 + +function HC.ModelKit.taylor!(u::AbstractVector, ::Val{N}, + sys::HomotopySystemWrapper{Inplace}, x, p = nothing) where {N} + f = sys.f + p = sys.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 + +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 + 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 + +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 + 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 + +""" + $(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 + +function GuessHomotopy(sys, fu0) + return GuessHomotopy(sys, fu0, HC.ModelKit.TaylorVector{5}(ComplexF64, length(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 * 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 * h.fu0[i] + end + return u, U +end + +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) + @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 diff --git a/lib/NonlinearSolveHomotopyContinuation/src/solve.jl b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl new file mode 100644 index 000000000..2db789df9 --- /dev/null +++ b/lib/NonlinearSolveHomotopyContinuation/src/solve.jl @@ -0,0 +1,164 @@ +""" + $(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 + prob.f + else + HomotopyNonlinearFunction(prob.f) + end + + u0 = state_values(prob) + p = parameter_values(prob) + isscalar = u0 isa Number + iip = SciMLBase.isinplace(prob) + + variant = iip ? Inplace : isscalar ? Scalar : OutOfPlace + + # jacobian handling + if SciMLBase.has_jac(f.f) + # use if present + prep = nothing + jac = f.jac + else + # 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 + vars = if isscalar + HC.variables(:x) + else + HC.variables(:x, axes(u0)...) + end + + taylorvars = if isscalar + [TaylorScalar(ntuple(Returns(0.0), 4)), TaylorScalar(ntuple(Returns(0.0), 4))] + elseif iip + ( + [TaylorScalar(ntuple(Returns(0.0), 4)) for _ in 1:2length(u0)], + [TaylorScalar(ntuple(Returns(0.0), 4)) for _ in 1:2length(u0)] + ) + else + [TaylorScalar(ntuple(Returns(0.0), 4)) for _ in 1:2length(u0)] + end + + jacobian_buffers = if isscalar + Matrix{Float64}(undef, 2, 2) + else + similar(u0, Float64, 2length(u0), 2length(u0)) + end + + # HC-compatible system + 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...) + 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, 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))) + 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 = real(only(realsols)) + T = eltype(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, u) + + 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 new file mode 100644 index 000000000..c0d025e49 --- /dev/null +++ b/lib/NonlinearSolveHomotopyContinuation/test/allroots.jl @@ -0,0 +1,151 @@ +using NonlinearSolve +using NonlinearSolveHomotopyContinuation +using SciMLBase: NonlinearSolution + +alg = HomotopyContinuationJL{true}(; threading = false) + +@testset "scalar u" begin + rhs = function (u, p) + return u * u - p[1] * u + p[2] + end + jac = function (u, p) + 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, [5.0, 6.0]) + sol = solve(prob, alg) + + @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 + @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 + 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, u - 0.9] + 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.9]) + 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 + +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")] + sol = nothing + @testset "`NonlinearProblem`" begin + 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 + 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 + 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 + @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 new file mode 100644 index 000000000..94bc8bddd --- /dev/null +++ b/lib/NonlinearSolveHomotopyContinuation/test/runtests.jl @@ -0,0 +1,15 @@ +using NonlinearSolveHomotopyContinuation +using Test +using Aqua + +@testset "NonlinearSolveHomotopyContinuation.jl" begin + @testset "Code quality (Aqua.jl)" begin + Aqua.test_all(NonlinearSolveHomotopyContinuation) + end + @testset "AllRoots" begin + include("allroots.jl") + end + @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..049578b37 --- /dev/null +++ b/lib/NonlinearSolveHomotopyContinuation/test/single_root.jl @@ -0,0 +1,129 @@ +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, u - 0.9] + 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.9]) + 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