test(NonlinearSolveSpectralMethods): add tests and ci scripts

Author: avik-pal

.github/workflows/CI_NonlinearSolveSpectralMethods.yml

@@ -0,0 +1,71 @@
+name: CI (NonlinearSolveSpectralMethods)
+
+on:
+  pull_request:
+    branches:
+      - master
+    paths:
+      - "lib/NonlinearSolveSpectralMethods/**"
+      - ".github/workflows/CI_NonlinearSolveSpectralMethods.yml"
+      - "lib/NonlinearSolveBase/**"
+      - "lib/SciMLJacobianOperators/**"
+  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:
+          - "lts"
+          - "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/SciMLJacobianOperators", "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/NonlinearSolveSpectralMethods {0}
+      - uses: julia-actions/julia-processcoverage@v1
+        with:
+          directories: lib/NonlinearSolveSpectralMethods/src,lib/NonlinearSolveBase/src,lib/NonlinearSolveBase/ext,lib/SciMLJacobianOperators/src
+      - uses: codecov/codecov-action@v4
+        with:
+          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          verbose: true
+          fail_ci_if_error: true

common/common_core_testing.jl

@@ -0,0 +1,54 @@
+using NonlinearSolveBase, SciMLBase
+
+const TERMINATION_CONDITIONS = [
+    NormTerminationMode(Base.Fix1(maximum, abs)),
+    RelTerminationMode(),
+    RelNormTerminationMode(Base.Fix1(maximum, abs)),
+    RelNormSafeTerminationMode(Base.Fix1(maximum, abs)),
+    RelNormSafeBestTerminationMode(Base.Fix1(maximum, abs)),
+    AbsTerminationMode(),
+    AbsNormTerminationMode(Base.Fix1(maximum, abs)),
+    AbsNormSafeTerminationMode(Base.Fix1(maximum, abs)),
+    AbsNormSafeBestTerminationMode(Base.Fix1(maximum, abs))
+]
+
+quadratic_f(u, p) = u .* u .- p
+quadratic_f!(du, u, p) = (du .= u .* u .- p)
+quadratic_f2(u, p) = @. p[1] * u * u - p[2]
+
+function newton_fails(u, p)
+    return 0.010000000000000002 .+
+           10.000000000000002 ./ (1 .+
+            (0.21640425613334457 .+
+             216.40425613334457 ./ (1 .+
+              (0.21640425613334457 .+
+               216.40425613334457 ./ (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^
+            2.0) .- 0.0011552453009332421u .- p
+end
+
+function solve_oop(f, u0, p = 2.0; solver, kwargs...)
+    prob = NonlinearProblem{false}(f, u0, p)
+    return solve(prob, solver; abstol = 1e-9, kwargs...)
+end
+
+function solve_iip(f, u0, p = 2.0; solver, kwargs...)
+    prob = NonlinearProblem{true}(f, u0, p)
+    return solve(prob, solver; abstol = 1e-9, kwargs...)
+end
+
+function nlprob_iterator_interface(f, p_range, isinplace, solver)
+    probN = NonlinearProblem{isinplace}(f, isinplace ? [0.5] : 0.5, p_range[begin])
+    cache = init(probN, solver; maxiters = 100, abstol = 1e-10)
+    sols = zeros(length(p_range))
+    for (i, p) in enumerate(p_range)
+        reinit!(cache, isinplace ? [cache.u[1]] : cache.u; p = p)
+        sol = solve!(cache)
+        sols[i] = isinplace ? sol.u[1] : sol.u
+    end
+    return sols
+end
+
+export TERMINATION_CONDITIONS
+export quadratic_f, quadratic_f!, quadratic_f2, newton_fails
+export solve_oop, solve_iip
+export nlprob_iterator_interface

lib/NonlinearSolveBase/src/utils.jl

@@ -4,7 +4,7 @@ using ArrayInterface: ArrayInterface
 using FastClosures: @closure
 using LinearAlgebra: LinearAlgebra, Diagonal, Symmetric, norm, dot, cond, diagind, pinv
 using MaybeInplace: @bb
-using RecursiveArrayTools: AbstractVectorOfArray, ArrayPartition
+using RecursiveArrayTools: AbstractVectorOfArray, ArrayPartition, recursivecopy!
 using SciMLOperators: AbstractSciMLOperator
 using SciMLBase: SciMLBase, AbstractNonlinearProblem, NonlinearFunction
 using StaticArraysCore: StaticArray, SArray, SMatrix
@@ -242,6 +242,17 @@ function make_identity!!(A::AbstractMatrix{T}, α) where {T}
     return A
 end
 
+function reinit_common!(cache, u0, p, alias_u0::Bool)
+    if SciMLBase.isinplace(cache)
+        recursivecopy!(cache.u, u0)
+        cache.prob.f(cache.fu, cache.u, p)
+    else
+        cache.u = maybe_unaliased(u0, alias_u0)
+        NonlinearSolveBase.set_fu!(cache, cache.prob.f(u0, p))
+    end
+    cache.p = p
+end
+
 function clean_sprint_struct(x)
     x isa Symbol && return "$(Meta.quot(x))"
     x isa Number && return string(x)

lib/NonlinearSolveSpectralMethods/Project.toml

@@ -8,7 +8,6 @@ CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MaybeInplace = "bb5d69b7-63fc-4a16-80bd-7e42200c7bdb"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -17,14 +16,14 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
 [compat]
 Aqua = "0.8"
+BenchmarkTools = "1.5.0"
 CommonSolve = "0.2.4"
 ConcreteStructs = "0.2.3"
 DiffEqBase = "6.155.3"
 ExplicitImports = "1.5"
 Hwloc = "3"
 InteractiveUtils = "<0.0.1, 1"
 LineSearch = "0.1.4"
-LinearAlgebra = "1.11.0"
 MaybeInplace = "0.1.4"
 NonlinearProblemLibrary = "0.1.2"
 NonlinearSolveBase = "1.1"
@@ -34,19 +33,22 @@ ReTestItems = "1.24"
 Reexport = "1"
 SciMLBase = "2.54"
 StableRNGs = "1"
+StaticArrays = "1.9.8"
 Test = "1.10"
 julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 Hwloc = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ExplicitImports", "Hwloc", "InteractiveUtils", "NonlinearProblemLibrary", "Pkg", "ReTestItems", "StableRNGs", "Test"]
+test = ["Aqua", "BenchmarkTools", "ExplicitImports", "Hwloc", "InteractiveUtils", "NonlinearProblemLibrary", "Pkg", "ReTestItems", "StableRNGs", "StaticArrays", "Test"]

lib/NonlinearSolveSpectralMethods/src/NonlinearSolveSpectralMethods.jl

@@ -1,12 +1,11 @@
 module NonlinearSolveSpectralMethods
 
+using ConcreteStructs: @concrete
 using Reexport: @reexport
 using PrecompileTools: @compile_workload, @setup_workload
 
 using CommonSolve: CommonSolve
-using ConcreteStructs: @concrete
 using DiffEqBase: DiffEqBase # Needed for `init` / `solve` dispatches
-using LinearAlgebra: dot
 using LineSearch: RobustNonMonotoneLineSearch
 using MaybeInplace: @bb
 using NonlinearSolveBase: NonlinearSolveBase, AbstractNonlinearSolveAlgorithm,
@@ -19,9 +18,7 @@ include("dfsane.jl")
 include("solve.jl")
 
 @setup_workload begin
-    include(joinpath(
-        @__DIR__, "..", "..", "..", "common", "nonlinear_problem_workloads.jl"
-    ))
+    include("../../../common/nonlinear_problem_workloads.jl")
 
     algs = [DFSane()]
 

lib/NonlinearSolveSpectralMethods/src/solve.jl

@@ -70,43 +70,41 @@ end
     kwargs
 end
 
-# XXX: Implement
-# function __reinit_internal!(
-#         cache::GeneralizedDFSaneCache{iip}, args...; p = cache.p, u0 = cache.u,
-#         alias_u0::Bool = false, maxiters = 1000, maxtime = nothing, kwargs...) where {iip}
-#     if iip
-#         recursivecopy!(cache.u, u0)
-#         cache.prob.f(cache.fu, cache.u, p)
-#     else
-#         cache.u = __maybe_unaliased(u0, alias_u0)
-#         set_fu!(cache, cache.prob.f(cache.u, p))
-#     end
-#     cache.p = p
-
-#     if cache.alg.σ_1 === nothing
-#         σ_n = dot(cache.u, cache.u) / dot(cache.u, cache.fu)
-#         # Spectral parameter bounds check
-#         if !(cache.alg.σ_min ≤ abs(σ_n) ≤ cache.alg.σ_max)
-#             test_norm = dot(cache.fu, cache.fu)
-#             σ_n = clamp(inv(test_norm), T(1), T(1e5))
-#         end
-#     else
-#         σ_n = T(cache.alg.σ_1)
-#     end
-#     cache.σ_n = σ_n
-
-#     reset_timer!(cache.timer)
-#     cache.total_time = 0.0
-
-#     reset!(cache.trace)
-#     reinit!(cache.termination_cache, get_fu(cache), get_u(cache); kwargs...)
-#     __reinit_internal!(cache.stats)
-#     cache.nsteps = 0
-#     cache.maxiters = maxiters
-#     cache.maxtime = maxtime
-#     cache.force_stop = false
-#     cache.retcode = ReturnCode.Default
-# end
+function InternalAPI.reinit_self!(
+        cache::GeneralizedDFSaneCache, args...; p = cache.p, u0 = cache.u,
+        alias_u0::Bool = false, maxiters = 1000, maxtime = nothing, kwargs...
+)
+    Utils.reinit_common!(cache, u0, p, alias_u0)
+
+    if cache.alg.σ_1 === nothing
+        σ_n = Utils.safe_dot(cache.u, cache.u) / Utils.safe_dot(cache.u, cache.fu)
+        # Spectral parameter bounds check
+        if !(cache.alg.σ_min ≤ abs(σ_n) ≤ cache.alg.σ_max)
+            test_norm = NonlinearSolveBase.L2_NORM(cache.fu)
+            σ_n = clamp(inv(test_norm), T(1), T(1e5))
+        end
+    else
+        σ_n = T(cache.alg.σ_1)
+    end
+    cache.σ_n = σ_n
+
+    NonlinearSolveBase.reset_timer!(cache.timer)
+    cache.total_time = 0.0
+
+    NonlinearSolveBase.reset!(cache.trace)
+    SciMLBase.reinit!(
+        cache.termination_cache, NonlinearSolveBase.get_fu(cache),
+        NonlinearSolveBase.get_u(cache); kwargs...
+    )
+
+    InternalAPI.reinit!(cache.stats)
+    cache.nsteps = 0
+    cache.maxiters = maxiters
+    cache.maxtime = maxtime
+    cache.force_stop = false
+    cache.retcode = ReturnCode.Default
+    return
+end
 
 NonlinearSolveBase.@internal_caches GeneralizedDFSaneCache :linesearch_cache
 
@@ -137,7 +135,7 @@ function SciMLBase.__init(
         )
 
         if alg.σ_1 === nothing
-            σ_n = dot(u, u) / dot(u, fu)
+            σ_n = Utils.safe_dot(u, u) / Utils.safe_dot(u, fu)
             # Spectral parameter bounds check
             if !(alg.σ_min ≤ abs(σ_n) ≤ alg.σ_max)
                 test_norm = NonlinearSolveBase.L2_NORM(fu)

lib/NonlinearSolveSpectralMethods/test/core_tests.jl

@@ -0,0 +1,83 @@
+@testsetup module CoreRootfindTesting
+
+include("../../../common/common_core_testing.jl")
+
+end
+
+@testitem "DFSane" setup=[CoreRootfindTesting] tags=[:core] begin
+    using BenchmarkTools: @ballocated
+    using StaticArrays: @SVector
+
+    u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
+
+    @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
+        sol = solve_oop(quadratic_f, u0; solver = DFSane())
+        @test SciMLBase.successful_retcode(sol)
+        @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+
+        cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), DFSane(), abstol = 1e-9)
+        @test (@ballocated solve!($cache)) < 200
+    end
+
+    @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
+        sol = solve_iip(quadratic_f!, u0; solver = DFSane())
+        @test SciMLBase.successful_retcode(sol)
+        @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+
+        cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), DFSane(), abstol = 1e-9)
+        @test (@ballocated solve!($cache)) ≤ 64
+    end
+end
+
+@testitem "DFSane Iterator Interface" setup=[CoreRootfindTesting] tags=[:core] begin
+    p = range(0.01, 2, length = 200)
+    @test nlprob_iterator_interface(quadratic_f, p, false, DFSane()) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f!, p, true, DFSane()) ≈ sqrt.(p)
+end
+
+@testitem "DFSane NewtonRaphson Fails" setup=[CoreRootfindTesting] tags=[:core] begin
+    u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
+    p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    sol = solve_oop(newton_fails, u0, p; solver = DFSane())
+    @test SciMLBase.successful_retcode(sol)
+    @test all(abs.(newton_fails(sol.u, p)) .< 1e-9)
+end
+
+@testitem "DFSane: Kwargs" setup=[CoreRootfindTesting] tags=[:core] begin
+    σ_min = [1e-10, 1e-5, 1e-4]
+    σ_max = [1e10, 1e5, 1e4]
+    σ_1 = [1.0, 0.5, 2.0]
+    M = [10, 1, 100]
+    γ = [1e-4, 1e-3, 1e-5]
+    τ_min = [0.1, 0.2, 0.3]
+    τ_max = [0.5, 0.8, 0.9]
+    nexp = [2, 1, 2]
+    η_strategy = [
+        (f_1, k, x, F) -> f_1 / k^2, (f_1, k, x, F) -> f_1 / k^3,
+        (f_1, k, x, F) -> f_1 / k^4
+    ]
+
+    list_of_options = zip(σ_min, σ_max, σ_1, M, γ, τ_min, τ_max, nexp, η_strategy)
+    for options in list_of_options
+        local probN, sol, alg
+        alg = DFSane(;
+            sigma_min = options[1], sigma_max = options[2], sigma_1 = options[3],
+            M = options[4], gamma = options[5], tau_min = options[6],
+            tau_max = options[7], n_exp = options[8], eta_strategy = options[9]
+        )
+
+        probN = NonlinearProblem{false}(quadratic_f, [1.0, 1.0], 2.0)
+        sol = solve(probN, alg, abstol = 1e-11)
+        @test all(abs.(quadratic_f(sol.u, 2.0)) .< 1e-6)
+    end
+end
+
+@testitem "DFSane Termination Conditions" setup=[CoreRootfindTesting] tags=[:core] begin
+    @testset "TC: $(nameof(typeof(termination_condition)))" for termination_condition in TERMINATION_CONDITIONS
+        @testset "u0: $(typeof(u0))" for u0 in ([1.0, 1.0], 1.0)
+            probN = NonlinearProblem(quadratic_f, u0, 2.0)
+            sol = solve(probN, DFSane(); termination_condition)
+            @test all(abs.(quadratic_f(sol.u, 2.0)) .< 1e-10)
+        end
+    end
+end