Move around the tests a bit

Author: avik-pal

.github/workflows/CI.yml

@@ -18,8 +18,8 @@ jobs:
       fail-fast: false
       matrix:
         group:
-          - Core
-          - NLLS
+          - RootFinding
+          - NLLSSolvers
           - 23TestProblems
           - Wrappers
           - Miscellaneous

Project.toml

@@ -97,6 +97,7 @@ SparseDiffTools = "2.14"
 SpeedMapping = "0.3"
 StableRNGs = "1"
 StaticArrays = "1.7"
+Sundials = "4.23.1"
 Symbolics = "5.13"
 Test = "1"
 UnPack = "1.0"
@@ -128,9 +129,10 @@ SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs", "MINPACK", "NLsolve", "OrdinaryDiffEq", "SpeedMapping", "FixedPointAcceleration", "SIAMFANLEquations"]
+test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs", "MINPACK", "NLsolve", "OrdinaryDiffEq", "SpeedMapping", "FixedPointAcceleration", "SIAMFANLEquations", "Sundials"]

docs/src/solvers/FixedPointSolvers.md

@@ -40,3 +40,10 @@ We are only listing the methods that natively solve fixed point problems.
 
   - `FixedPointAccelerationJL()`: accelerates the convergence of a mapping to a fixed point
     by the Anderson acceleration algorithm and a few other methods.
+
+### NLsolve.jl
+
+In our tests, we have found the anderson method implemented here to NOT be the most
+robust.
+
+  - `NLsolveJL(; method = :anderson)`: Anderson acceleration for fixed point problems.

ext/NonlinearSolveFastLevenbergMarquardtExt.jl

@@ -15,6 +15,7 @@ import FiniteDiff, ForwardDiff
     end
 end
 
+# TODO: Implement reinit
 @concrete struct FastLevenbergMarquardtJLCache
     f!
     J!

ext/NonlinearSolveFixedPointAccelerationExt.jl

@@ -9,41 +9,25 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::FixedPointAccelerationJL
     @assert (termination_condition ===
              nothing)||(termination_condition isa AbsNormTerminationMode) "FixedPointAccelerationJL does not support termination conditions!"
 
-    u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
-    u_size = size(u0)
-    T = eltype(u0)
-    iip = isinplace(prob)
-    p = prob.p
-
-    if !iip && prob.u0 isa Number
-        # FixedPointAcceleration makes the scalar problem into a vector problem
-        f = (u) -> [prob.f(u[1], p) .+ u[1]]
-    elseif !iip && prob.u0 isa AbstractVector
-        f = (u) -> (prob.f(u, p) .+ u)
-    elseif !iip && prob.u0 isa AbstractArray
-        f = (u) -> vec(prob.f(reshape(u, u_size), p) .+ u)
-    elseif iip && prob.u0 isa AbstractVector
-        du = similar(u0)
-        f = (u) -> (prob.f(du, u, p); du .+ u)
-    else
-        du = similar(u0)
-        f = (u) -> (prob.f(du, reshape(u, u_size), p); vec(du) .+ u)
-    end
+    f, u0 = NonlinearSolve.__construct_f(prob; alias_u0, make_fixed_point = Val(true),
+        force_oop = Val(true))
 
     tol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u0))
 
-    sol = fixed_point(f, NonlinearSolve._vec(u0); Algorithm = alg.algorithm,
+    sol = fixed_point(f, u0; Algorithm = alg.algorithm,
         ConvergenceMetricThreshold = tol, MaxIter = maxiters, MaxM = alg.m,
         ExtrapolationPeriod = alg.extrapolation_period, Dampening = alg.dampening,
         PrintReports, ReplaceInvalids = alg.replace_invalids,
         ConditionNumberThreshold = alg.condition_number_threshold, quiet_errors = true)
 
-    res = prob.u0 isa Number ? first(sol.FixedPoint_) : sol.FixedPoint_
-    if res === missing
+    if sol.FixedPoint_ === missing
+        u0 = prob.u0 isa Number ? u0[1] : u0
         resid = NonlinearSolve.evaluate_f(prob, u0)
         res = u0
         converged = false
     else
+        res = prob.u0 isa Number ? first(sol.FixedPoint_) :
+              reshape(sol.FixedPoint_, size(prob.u0))
         resid = NonlinearSolve.evaluate_f(prob, res)
         converged = maximum(abs, resid) ≤ tol
     end

ext/NonlinearSolveMINPACKExt.jl

@@ -31,7 +31,7 @@ function SciMLBase.__solve(prob::Union{NonlinearProblem{uType, iip},
         original = MINPACK.fsolve(f!, u0, m; tol, show_trace, tracing, method,
             iterations = maxiters)
     else
-        jac! = @closure((J, u) -> (jac!_(J, u); Cint(0)))
+        jac! = @closure((J, u)->(jac!_(J, u); Cint(0)))
         original = MINPACK.fsolve(f!, jac!, u0, m; tol, show_trace, tracing, method,
             iterations = maxiters)
     end

ext/NonlinearSolveSIAMFANLEquationsExt.jl

@@ -58,7 +58,8 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
 
         retcode = __siam_fanl_equations_retcode_mapping(sol)
         stats = __siam_fanl_equations_stats_mapping(method, sol)
-        return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode,
+        resid = NonlinearSolve.evaluate_f(prob, sol.solution)
+        return SciMLBase.build_solution(prob, alg, sol.solution, resid; retcode,
             stats, original = sol)
     end
 
@@ -87,7 +88,8 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
 
         retcode = __siam_fanl_equations_retcode_mapping(sol)
         stats = __siam_fanl_equations_stats_mapping(method, sol)
-        return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode,
+        resid = NonlinearSolve.evaluate_f(prob, sol.solution)
+        return SciMLBase.build_solution(prob, alg, sol.solution, resid; retcode,
             stats, original = sol)
     end
 
@@ -116,7 +118,8 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
 
     retcode = __siam_fanl_equations_retcode_mapping(sol)
     stats = __siam_fanl_equations_stats_mapping(method, sol)
-    return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode, stats,
+    resid = NonlinearSolve.evaluate_f(prob, sol.solution)
+    return SciMLBase.build_solution(prob, alg, sol.solution, resid; retcode, stats,
         original = sol)
 end
 

ext/NonlinearSolveSpeedMappingExt.jl

@@ -17,7 +17,7 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SpeedMappingJL, args...;
     sol = speedmapping(u0; m!, tol, Lp = Inf, maps_limit = maxiters, alg.orders,
         alg.check_obj, store_info, alg.σ_min, alg.stabilize)
     res = prob.u0 isa Number ? first(sol.minimizer) : sol.minimizer
-    resid = NonlinearSolve.evaluate_f(prob, sol.minimizer)
+    resid = NonlinearSolve.evaluate_f(prob, res)
 
     return SciMLBase.build_solution(prob, alg, res, resid;
         retcode = sol.converged ? ReturnCode.Success : ReturnCode.Failure,

src/function_wrappers.jl

@@ -2,11 +2,14 @@
 # downstream packages. Make conversion to those easier.
 function __construct_f(prob; alias_u0::Bool = false, can_handle_oop::Val{OOP} = Val(false),
         can_handle_scalar::Val{SCALAR} = Val(false), make_fixed_point::Val{FP} = Val(false),
-        can_handle_arbitrary_dims::Val{DIMS} = Val(false)) where {SCALAR, OOP, DIMS, FP}
+        can_handle_arbitrary_dims::Val{DIMS} = Val(false),
+        force_oop::Val{FOOP} = Val(false)) where {SCALAR, OOP, DIMS, FP, FOOP}
     if !OOP && SCALAR
         error("Incorrect Specification: OOP not supported but scalar supported.")
     end
 
+    resid = evaluate_f(prob, prob.u0)
+
     if SCALAR || !(prob.u0 isa Number)
         u0 = __maybe_unaliased(prob.u0, alias_u0)
     else
@@ -20,51 +23,69 @@ function __construct_f(prob; alias_u0::Bool = false, can_handle_oop::Val{OOP} =
                 @. du += u
             end
         else
-            @closure (du, u, p) -> prob.f(du, u, p) .+ u
+            @closure (u, p) -> prob.f(u, p) .+ u
         end
     else
         prob.f
     end
 
-    f_final = if isinplace(prob)
+    ff = if isinplace(prob)
+        ninputs = 2
         if DIMS || u0 isa AbstractVector
             @closure (du, u) -> (f(du, u, prob.p); du)
         else
             u0_size = size(u0)
-            du_size = size(evaluate_f(prob, u0))
+            du_size = size(resid)
             @closure (du, u) -> (f(reshape(du, du_size), reshape(u, u0_size), prob.p); du)
         end
     else
         if prob.u0 isa Number
             if SCALAR
-                @closure (u) -> prob.f(u, prob.p)
+                ninputs = 1
+                @closure (u) -> f(u, prob.p)
             elseif OOP
-                @closure (u) -> [prob.f(first(u), prob.p)]
+                ninputs = 1
+                @closure (u) -> [f(first(u), prob.p)]
             else
-                @closure (du, u) -> (du[1] = prob.f(first(u), prob.p); du)
+                ninputs = 2
+                resid = [resid]
+                @closure (du, u) -> (du[1] = f(first(u), prob.p); du)
             end
         else
             if OOP
+                ninputs = 1
                 if DIMS
-                    @closure (u) -> prob.f(u, prob.p)
+                    @closure (u) -> f(u, prob.p)
                 else
                     u0_size = size(u0)
-                    @closure (u) -> _vec(prob.f(reshape(u, u0_size), prob.p))
+                    @closure (u) -> _vec(f(reshape(u, u0_size), prob.p))
                 end
             else
+                ninputs = 2
                 if DIMS
-                    @closure (du, u) -> (copyto!(du, prob.f(u, prob.p)); du)
+                    @closure (du, u) -> (copyto!(du, f(u, prob.p)); du)
                 else
                     u0_size = size(u0)
                     @closure (du, u) -> begin
-                        copyto!(vec(du), vec(prob.f(reshape(u, u0_size), prob.p)))
+                        copyto!(vec(du), vec(f(reshape(u, u0_size), prob.p)))
                         return du
                     end
                 end
             end
         end
     end
 
+    f_final = if FOOP
+        if ninputs == 1
+            ff
+        else
+            du_ = DIMS ? similar(resid) : _vec(similar(resid))
+            @closure (u) -> (ff(du_, u); du_)
+        end
+    else
+        ff
+    end
+
     return f_final, ifelse(DIMS, u0, _vec(u0))
 end