From 343af4048e7bf111fa4be9d17ec25178c342a0c1 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 14:12:00 -0400
Subject: [PATCH 01/10] Cleanup the PolyAlg Code

---
 docs/src/solvers/NonlinearSystemSolvers.md |   8 +-
 src/default.jl                             | 312 ++++++++-------------
 src/trustRegion.jl                         |  20 +-
 test/polyalgs.jl                           |   8 +-
 4 files changed, 130 insertions(+), 218 deletions(-)

diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md
index c4f0a66fb..740654d57 100644
--- a/docs/src/solvers/NonlinearSystemSolvers.md
+++ b/docs/src/solvers/NonlinearSystemSolvers.md
@@ -16,7 +16,7 @@ it will make sure to work.
 
 If one is looking for more robustness then `RobustMultiNewton` is a good choice.
 It attempts a set of the most robust methods in succession and only fails if
-all of the methods fail to converge. Additionally, `DynamicSS` can be a good choice 
+all of the methods fail to converge. Additionally, `DynamicSS` can be a good choice
 for high stability.
 
 As a balance, `NewtonRaphson` is a good choice for most problems that aren't too
@@ -25,13 +25,13 @@ but more stable. If the problem is well-conditioned, `Klement` or `Broyden`
 may be faster, but highly dependent on the eigenvalues of the Jacobian being
 sufficiently small.
 
-`NewtonRaphson` and `TrustRegion` are designed for for large systems. 
+`NewtonRaphson` and `TrustRegion` are designed for for large systems.
 They can make use of sparsity patterns for sparse automatic differentiation
 and sparse linear solving of very large systems. Meanwhile,
 `SimpleNewtonRaphson` and `SimpleTrustRegion` are implementations which is specialized for
 small equations. They are non-allocating on static arrays and thus really well-optimized
 for small systems, thus usually outperforming the other methods when such types are
-used for `u0`. 
+used for `u0`.
 
 ## Full List of Methods
 
@@ -57,7 +57,7 @@ features, but have a bit of overhead on very small problems.
     large-scale and numerically-difficult nonlinear systems.
   - `RobustMultiNewton()`: A polyalgorithm that mixes highly robust methods (line searches and
     trust regions) in order to be as robust as possible for difficult problems. If this method
-    fails to converge, then one can be pretty certain that most (all?) other choices would 
+    fails to converge, then one can be pretty certain that most (all?) other choices would
     likely fail.
   - `FastShortcutNonlinearPolyalg`: The default method. A polyalgorithm that mixes fast methods
     with fallbacks to robust methods to allow for solving easy problems quickly without sacrificing
diff --git a/src/default.jl b/src/default.jl
index 44ac26cad..36ac42225 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -1,6 +1,6 @@
 """
-RobustMultiNewton(; concrete_jac = nothing, linsolve = nothing,
-                    precs = DEFAULT_PRECS, adkwargs...)
+    RobustMultiNewton(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
+                        adkwargs...)
 
 A polyalgorithm focused on robustness. It uses a mixture of Newton methods with different
 globalizing techniques (trust region updates, line searches, etc.) in order to find a
@@ -12,20 +12,20 @@ or more precision / more stable linear solver choice is required).
 
 ### Keyword Arguments
 
-- `autodiff`: determines the backend used for the Jacobian. Note that this argument is
+  - `autodiff`: determines the backend used for the Jacobian. Note that this argument is
     ignored if an analytical Jacobian is passed, as that will be used instead. Defaults to
     `AutoForwardDiff()`. Valid choices are types from ADTypes.jl.
-- `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used,
+  - `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used,
     then the Jacobian will not be constructed and instead direct Jacobian-vector products
     `J*v` are computed using forward-mode automatic differentiation or finite differencing
     tricks (without ever constructing the Jacobian). However, if the Jacobian is still needed,
     for example for a preconditioner, `concrete_jac = true` can be passed in order to force
     the construction of the Jacobian.
-- `linsolve`: the [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) used for the
+  - `linsolve`: the [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) used for the
     linear solves within the Newton method. Defaults to `nothing`, which means it uses the
     LinearSolve.jl default algorithm choice. For more information on available algorithm
     choices, see the [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
-- `precs`: the choice of preconditioners for the linear solver. Defaults to using no
+  - `precs`: the choice of preconditioners for the linear solver. Defaults to using no
     preconditioners. For more information on specifying preconditioners for LinearSolve
     algorithms, consult the
     [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
@@ -41,9 +41,8 @@ end
 # Robusters!
 const Robusters = RobustMultiNewton
 
-function RobustMultiNewton(; concrete_jac = nothing, linsolve = nothing, 
-        precs = DEFAULT_PRECS, adkwargs...)
-
+function RobustMultiNewton(; concrete_jac = nothing, linsolve = nothing,
+    precs = DEFAULT_PRECS, adkwargs...)
     return RobustMultiNewton{_unwrap_val(concrete_jac)}(adkwargs, linsolve, precs)
 end
 
@@ -53,45 +52,46 @@ end
     current::Int
 end
 
-function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::RobustMultiNewton, args...; 
-    kwargs...) where {uType, iip}
-
-    adkwargs = alg.adkwargs
-    linsolve = alg.linsolve
-    precs = alg.precs
-
-    RobustMultiNewtonCache{iip}((
-        SciMLBase.__init(prob, TrustRegion(;linsolve, precs, adkwargs...), args...; kwargs...),
-        SciMLBase.__init(prob, TrustRegion(;linsolve, precs, radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...), args...; kwargs...),
-        SciMLBase.__init(prob, NewtonRaphson(;linsolve, precs, linesearch=BackTracking(), adkwargs...), args...; kwargs...),
-        SciMLBase.__init(prob, TrustRegion(;linsolve, precs, radius_update_scheme = RadiusUpdateSchemes.Fan, adkwargs...), args...; kwargs...),
-        ), alg, 1
-    )
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::RobustMultiNewton,
+    args...; kwargs...) where {uType, iip}
+    @unpack adkwargs, linsolve, precs = alg
+
+    algs = (TrustRegion(; linsolve, precs, adkwargs...),
+        TrustRegion(; linsolve, precs,
+            radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...),
+        NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...),
+        TrustRegion(; linsolve, precs,
+            radius_update_scheme = RadiusUpdateSchemes.NLsolve, adkwargs...),
+        TrustRegion(; linsolve, precs,
+            radius_update_scheme = RadiusUpdateSchemes.Fan, adkwargs...))
+
+    return RobustMultiNewtonCache{iip}(map(solver -> SciMLBase.__init(prob, solver, args...;
+                kwargs...), algs), alg, 1)
 end
 
 """
-FastShortcutNonlinearPolyalg(; concrete_jac = nothing, linsolve = nothing,
-                    precs = DEFAULT_PRECS, adkwargs...)
+    FastShortcutNonlinearPolyalg(; concrete_jac = nothing, linsolve = nothing,
+                                   precs = DEFAULT_PRECS, adkwargs...)
 
 A polyalgorithm focused on balancing speed and robustness. It first tries less robust methods
 for more performance and then tries more robust techniques if the faster ones fail.
 
 ### Keyword Arguments
 
-- `autodiff`: determines the backend used for the Jacobian. Note that this argument is
+  - `autodiff`: determines the backend used for the Jacobian. Note that this argument is
     ignored if an analytical Jacobian is passed, as that will be used instead. Defaults to
     `AutoForwardDiff()`. Valid choices are types from ADTypes.jl.
-- `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used,
+  - `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used,
     then the Jacobian will not be constructed and instead direct Jacobian-vector products
     `J*v` are computed using forward-mode automatic differentiation or finite differencing
     tricks (without ever constructing the Jacobian). However, if the Jacobian is still needed,
     for example for a preconditioner, `concrete_jac = true` can be passed in order to force
     the construction of the Jacobian.
-- `linsolve`: the [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) used for the
+  - `linsolve`: the [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) used for the
     linear solves within the Newton method. Defaults to `nothing`, which means it uses the
     LinearSolve.jl default algorithm choice. For more information on available algorithm
     choices, see the [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
-- `precs`: the choice of preconditioners for the linear solver. Defaults to using no
+  - `precs`: the choice of preconditioners for the linear solver. Defaults to using no
     preconditioners. For more information on specifying preconditioners for LinearSolve
     algorithms, consult the
     [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
@@ -102,95 +102,76 @@ for more performance and then tries more robust techniques if the faster ones fa
     precs
 end
 
-function FastShortcutNonlinearPolyalg(; concrete_jac = nothing, linsolve = nothing, 
+function FastShortcutNonlinearPolyalg(; concrete_jac = nothing, linsolve = nothing,
     precs = DEFAULT_PRECS, adkwargs...)
-
-    return FastShortcutNonlinearPolyalg{_unwrap_val(concrete_jac)}(adkwargs, linsolve, precs)
+    return FastShortcutNonlinearPolyalg{_unwrap_val(concrete_jac)}(adkwargs, linsolve,
+        precs)
 end
 
-@concrete mutable struct FastShortcutNonlinearPolyalgCache{iip} <: AbstractNonlinearSolveCache{iip}
+@concrete mutable struct FastShortcutNonlinearPolyalgCache{iip} <:
+                         AbstractNonlinearSolveCache{iip}
     caches
     alg
     current::Int
 end
 
-function FastShortcutNonlinearPolyalgCache(; concrete_jac = nothing, linsolve = nothing, 
-        precs = DEFAULT_PRECS, adkwargs...)
-
-    return FastShortcutNonlinearPolyalgCache{_unwrap_val(concrete_jac)}(adkwargs, linsolve, precs)
+function FastShortcutNonlinearPolyalgCache(; concrete_jac = nothing, linsolve = nothing,
+    precs = DEFAULT_PRECS, adkwargs...)
+    return FastShortcutNonlinearPolyalgCache{_unwrap_val(concrete_jac)}(adkwargs, linsolve,
+        precs)
 end
 
-function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::FastShortcutNonlinearPolyalg, args...; 
-    kwargs...) where {uType, iip}
-
-    adkwargs = alg.adkwargs
-    linsolve = alg.linsolve
-    precs = alg.precs
-
-    FastShortcutNonlinearPolyalgCache{iip}((
-        #SciMLBase.__init(prob, Klement(), args...; kwargs...),
-        #SciMLBase.__init(prob, Broyden(), args...; kwargs...),
-        SciMLBase.__init(prob, NewtonRaphson(;linsolve, precs, adkwargs...), args...; kwargs...),
-        SciMLBase.__init(prob, NewtonRaphson(;linsolve, precs, linesearch=BackTracking(), adkwargs...), args...; kwargs...),
-        SciMLBase.__init(prob, TrustRegion(;linsolve, precs, adkwargs...), args...; kwargs...),
-        SciMLBase.__init(prob, TrustRegion(;linsolve, precs, radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...), args...; kwargs...),
-        ), alg, 1
-    )
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip},
+    alg::FastShortcutNonlinearPolyalg, args...; kwargs...) where {uType, iip}
+    @unpack adkwargs, linsolve, precs = alg
+
+    algs = (
+        # Klement(),
+        # Broyden(),
+        NewtonRaphson(; linsolve, precs, adkwargs...),
+        NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...),
+        TrustRegion(; linsolve, precs, adkwargs...),
+        TrustRegion(; linsolve, precs,
+            radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...))
+
+    return FastShortcutNonlinearPolyalgCache{iip}(map(solver -> SciMLBase.__init(prob,
+                solver, args...; kwargs...), algs), alg, 1)
 end
 
-function SciMLBase.__solve(prob::NonlinearProblem{uType, false}, alg::FastShortcutNonlinearPolyalg, args...; 
-    kwargs...) where {uType}
-
-    adkwargs = alg.adkwargs
-    linsolve = alg.linsolve
-    precs = alg.precs
-
-    sol1 = SciMLBase.__solve(prob, Klement(), args...; kwargs...)
-    if SciMLBase.successful_retcode(sol1)
-        return SciMLBase.build_solution(prob, alg, sol1.u, sol1.resid;
-                                        sol1.retcode, sol1.stats)
-    end
-
-    sol2 = SciMLBase.__solve(prob, Broyden(), args...; kwargs...)
-    if SciMLBase.successful_retcode(sol2)
-        return SciMLBase.build_solution(prob, alg, sol2.u, sol2.resid;
-                                        sol2.retcode, sol2.stats)
-    end
-
-    sol3 = SciMLBase.__solve(prob, NewtonRaphson(;linsolve, precs, adkwargs...), args...; kwargs...)
-    if SciMLBase.successful_retcode(sol3)
-        return SciMLBase.build_solution(prob, alg, sol3.u, sol3.resid;
-                                        sol3.retcode, sol3.stats)
-    end
-
-    sol4 = SciMLBase.__solve(prob,  TrustRegion(;linsolve, precs, adkwargs...), args...; kwargs...)
-    if SciMLBase.successful_retcode(sol4)
-        return SciMLBase.build_solution(prob, alg, sol4.u, sol4.resid;
-                                        sol4.retcode, sol4.stats)
+# This version doesn't allocate all the caches!
+function SciMLBase.__solve(prob::NonlinearProblem{uType, iip},
+    alg::FastShortcutNonlinearPolyalg, args...; kwargs...) where {uType, iip}
+    @unpack adkwargs, linsolve, precs = alg
+
+    algs = [
+        iip ? Klement() : nothing, # Klement not yet implemented for IIP
+        iip ? Broyden() : nothing, # Broyden not yet implemented for IIP
+        NewtonRaphson(; linsolve, precs, adkwargs...),
+        NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...),
+        TrustRegion(; linsolve, precs, adkwargs...),
+        TrustRegion(; linsolve, precs,
+            radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...),
+    ]
+    filter!(!isnothing, algs)
+
+    sols = Vector{SciMLBase.NonlinearSolution}(undef, length(algs))
+
+    for (i, solver) in enumerate(algs)
+        sols[i] = SciMLBase.__solve(prob, solver, args...; kwargs...)
+        if SciMLBase.successful_retcode(sols[i])
+            return SciMLBase.build_solution(prob, alg, sols[i].u, sols[i].resid;
+                sols[i].retcode, sols[i].stats, original = sols[i])
+        end
     end
 
-    resids = (sol1.resid, sol2.resid, sol3.resid, sol4.resid)
+    resids = map(Base.Fix2(getproperty, resid), sols)
     minfu, idx = findmin(DEFAULT_NORM, resids)
 
-    if idx == 1
-        SciMLBase.build_solution(prob, alg, sol1.u, sol1.resid;
-                                        sol1.retcode, sol1.stats)
-    elseif idx == 2
-        SciMLBase.build_solution(prob, alg, sol2.u, sol2.resid;
-                                        sol2.retcode, sol2.stats)
-    elseif idx == 3
-        SciMLBase.build_solution(prob, alg, sol3.u, sol3.resid;
-                                sol3.retcode, sol3.stats)
-    elseif idx == 4
-        SciMLBase.build_solution(prob, alg, sol4.u, sol4.resid;
-                                sol4.retcode, sol4.stats)
-    else
-        error("Unreachable reached, 박정석")
-    end
-
+    return SciMLBase.build_solution(prob, alg, sols[idx].u, sols[idx].resid;
+        sols[idx].retcode, sols[idx].stats, original = sols[idx])
 end
 
-function SciMLBase.__solve(prob::NonlinearProblem{uType, true}, alg::FastShortcutNonlinearPolyalg, args...; 
+function SciMLBase.__solve(prob::NonlinearProblem{uType, true}, alg::FastShortcutNonlinearPolyalg, args...;
     kwargs...) where {uType}
 
     adkwargs = alg.adkwargs
@@ -244,122 +225,57 @@ end
 
 ## General shared polyalg functions
 
-function perform_step!(cache::Union{RobustMultiNewtonCache, FastShortcutNonlinearPolyalgCache})
+function perform_step!(cache::Union{RobustMultiNewtonCache,
+    FastShortcutNonlinearPolyalgCache})
     current = cache.current
+    1 ≤ current ≤ length(cache.caches) || error("Current choices shouldn't get here!")
 
-    while true
-        if current == 1
-            perform_step!(cache.caches[1])
-        elseif current == 2
-            perform_step!(cache.caches[2])
-        elseif current == 3
-            perform_step!(cache.caches[3])
-        elseif current == 4
-            perform_step!(cache.caches[4])
-        else
-            error("Current choices shouldn't get here!")
-        end
+    current_cache = cache.caches[current]
+    while not_terminated(current_cache)
+        perform_step!(current_cache)
     end
 
     return nothing
 end
 
-function SciMLBase.solve!(cache::Union{RobustMultiNewtonCache, FastShortcutNonlinearPolyalgCache})
+function SciMLBase.solve!(cache::Union{RobustMultiNewtonCache,
+    FastShortcutNonlinearPolyalgCache})
     current = cache.current
-
-    while current < 5 && all(not_terminated, cache.caches)
-        if current == 1
-            perform_step!(cache.caches[1])
-            !not_terminated(cache.caches[1]) && (cache.current += 1)
-        elseif current == 2
-            perform_step!(cache.caches[2])
-            !not_terminated(cache.caches[2]) && (cache.current += 1)
-        elseif current == 3
-            perform_step!(cache.caches[3])
-            !not_terminated(cache.caches[3]) && (cache.current += 1)
-        elseif current == 4
-            perform_step!(cache.caches[4])
-            !not_terminated(cache.caches[4]) && (cache.current += 1)
-        else
-            error("Current choices shouldn't get here!")
-        end
-
-        
-
-        #cache.stats.nsteps += 1
+    1 ≤ current ≤ length(cache.caches) || error("Current choices shouldn't get here!")
+
+    current_cache = cache.caches[current]
+    while current ≤ length(cache.caches) #  && !all(terminated[current:end])
+        sol_tmp = solve!(current_cache)
+        SciMLBase.successful_retcode(sol_tmp) && break
+        current += 1
+        cache.current = current
+        current_cache = cache.caches[current]
     end
 
-    if current < 5
-        stats = if current == 1
-            cache.caches[1].stats
-        elseif current == 2
-            cache.caches[2].stats
-        elseif current == 3
-            cache.caches[3].stats
-        elseif current == 4
-            cache.caches[4].stats
-        end
+    if current ≤ length(cache.caches)
+        retcode = ReturnCode.Success
 
-        u = if current == 1
-            cache.caches[1].u
-        elseif current == 2
-            cache.caches[2].u
-        elseif current == 3
-            cache.caches[3].u
-        elseif current == 4
-            cache.caches[4].u
-        end
-
-        fu = if current == 1
-            get_fu(cache.caches[1])
-        elseif current == 2
-            get_fu(cache.caches[2])
-        elseif current == 3
-            get_fu(cache.caches[3])
-        elseif current == 4
-            get_fu(cache.caches[4])
-        end
-
-        retcode = if stats.nsteps == cache.caches[1].maxiters
-            ReturnCode.MaxIters
-        else
-            ReturnCode.Success
-        end
+        stats = cache.caches[current].stats
+        u = cache.caches[current].u
+        fu = get_fu(cache.caches[current])
 
         return SciMLBase.build_solution(cache.caches[1].prob, cache.alg, u, fu;
             retcode, stats)
     else
         retcode = ReturnCode.MaxIters
 
-        fus = (get_fu(cache.caches[1]), get_fu(cache.caches[2]), get_fu(cache.caches[3]), get_fu(cache.caches[4]))
+        fus = get_fu.(cache.caches)
         minfu, idx = findmin(cache.caches[1].internalnorm, fus)
+        stats = cache.caches[idx].stats
+        u = cache.caches[idx].u
 
-        stats = if idx == 1
-            cache.caches[1].stats
-        elseif idx == 2
-            cache.caches[2].stats
-        elseif idx == 3
-            cache.caches[3].stats
-        elseif idx == 4
-            cache.caches[4].stats
-        end
-
-        u = if idx == 1
-            cache.caches[1].u
-        elseif idx == 2
-            cache.caches[2].u
-        elseif idx == 3
-            cache.caches[3].u
-        elseif idx == 4
-            cache.caches[4].u
-        end
-        
-        return SciMLBase.build_solution(cache.caches[1].prob, cache.alg, u, fu;
-                                        retcode, stats)
+        return SciMLBase.build_solution(cache.caches[idx].prob, cache.alg, u, fus[idx];
+            retcode, stats)
     end
 end
 
-function SciMLBase.reinit!(cache::Union{RobustMultiNewtonCache, FastShortcutNonlinearPolyalgCache}, args...; kwargs...)
+function SciMLBase.reinit!(cache::Union{RobustMultiNewtonCache,
+        FastShortcutNonlinearPolyalgCache}, args...; kwargs...)
     for c in cache.caches
         SciMLBase.reinit!(c, args...; kwargs...)
     end
@@ -367,14 +283,12 @@ end
 
 ## Defaults
 
-function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::Nothing, args...; 
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::Nothing, args...;
     kwargs...) where {uType, iip}
-
     SciMLBase.__init(prob, FastShortcutNonlinearPolyalg(), args...; kwargs...)
 end
 
-function SciMLBase.__solve(prob::NonlinearProblem{uType, iip}, alg::Nothing, args...; 
+function SciMLBase.__solve(prob::NonlinearProblem{uType, iip}, alg::Nothing, args...;
     kwargs...) where {uType, iip}
-
     SciMLBase.__solve(prob, FastShortcutNonlinearPolyalg(), args...; kwargs...)
-end
\ No newline at end of file
+end
diff --git a/src/trustRegion.jl b/src/trustRegion.jl
index 50d1078f3..583b1b675 100644
--- a/src/trustRegion.jl
+++ b/src/trustRegion.jl
@@ -28,14 +28,14 @@ states as `RadiusUpdateSchemes.T`. Simply put the desired scheme as follows:
     """
     `RadiusUpdateSchemes.NLsolve`
 
-    The same updating scheme as in NLsolve's (https://github.com/JuliaNLSolvers/NLsolve.jl) trust region dogleg implementation. 
+    The same updating scheme as in NLsolve's (https://github.com/JuliaNLSolvers/NLsolve.jl) trust region dogleg implementation.
     """
     NLsolve
 
     """
     `RadiusUpdateSchemes.NocedalWright`
 
-    Trust region updating scheme as in Nocedal and Wright [see Alg 11.5, page 291]. 
+    Trust region updating scheme as in Nocedal and Wright [see Alg 11.5, page 291].
     """
     NocedalWright
 
@@ -239,7 +239,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion,
     u_gauss_newton = zero(u)
 
     loss_new = loss
-    H = zero(J)
+    H = zero(J' * J)
     g = _mutable_zero(fu1)
     shrink_counter = 0
     fu_new = zero(fu1)
@@ -341,7 +341,7 @@ function perform_step!(cache::TrustRegionCache{true})
         mul!(cache.g, J', fu)
         cache.stats.njacs += 1
 
-        # do not use A = cache.H, b = _vec(cache.g) since it is equivalent 
+        # do not use A = cache.H, b = _vec(cache.g) since it is equivalent
         # to  A = cache.J, b = _vec(fu) as long as the Jacobian is non-singular
         linres = dolinsolve(alg.precs, linsolve, A = J, b = _vec(fu),
             linu = _vec(u_gauss_newton),
@@ -450,12 +450,12 @@ function trust_region_step!(cache::TrustRegionCache)
             cache.make_new_J = false
         end
 
-        # trust region update 
-        if r < 1 // 10 # cache.shrink_threshold 
-            cache.trust_r *= 1 // 2 # cache.shrink_factor 
-        elseif r >= 9 // 10 # cache.expand_threshold 
-            cache.trust_r = 2 * norm(cache.du) # cache.expand_factor * norm(cache.du) 
-        elseif r >= 1 // 2 # cache.p1 
+        # trust region update
+        if r < 1 // 10 # cache.shrink_threshold
+            cache.trust_r *= 1 // 2 # cache.shrink_factor
+        elseif r >= 9 // 10 # cache.expand_threshold
+            cache.trust_r = 2 * norm(cache.du) # cache.expand_factor * norm(cache.du)
+        elseif r >= 1 // 2 # cache.p1
             cache.trust_r = max(cache.trust_r, 2 * norm(cache.du)) # cache.expand_factor * norm(cache.du))
         end
 
diff --git a/test/polyalgs.jl b/test/polyalgs.jl
index 0bdfe42f8..3aa552b24 100644
--- a/test/polyalgs.jl
+++ b/test/polyalgs.jl
@@ -8,7 +8,6 @@ probN = NonlinearProblem(f, u0)
 @time solver = solve(probN, FastShortcutNonlinearPolyalg(), abstol = 1e-9)
 
 # https://github.com/SciML/NonlinearSolve.jl/issues/153
-
 function f(du, u, p)
     s1, s1s2, s2 = u
     k1, c1, Δt = p
@@ -18,13 +17,12 @@ function f(du, u, p)
     du[3] = -0.25 * c1 * k1 * s1 * s2
 end
 
-prob = NonlinearProblem(f, [2.0,2.0,2.0], [1.0, 2.0, 2.5])
+prob = NonlinearProblem(f, [2.0, 2.0, 2.0], [1.0, 2.0, 2.5])
 sol = solve(prob)
 @test SciMLBase.successful_retcode(sol)
 
 # https://github.com/SciML/NonlinearSolve.jl/issues/187
-
-ff(u, p) = 0.5/1.5*log.(u./(1.0.-u)) .- 2.0*u .+1.0
+ff(u, p) = 0.5 / 1.5 * log.(u ./ (1.0 .- u)) .- 2.0 * u .+ 1.0
 
 uspan = (0.02, 0.1)
 prob = IntervalNonlinearProblem(ff, uspan)
@@ -35,4 +33,4 @@ u0 = 0.06
 p = 2.0
 prob = NonlinearProblem(ff, u0, p)
 solver = solve(prob)
-@test SciMLBase.successful_retcode(sol)
\ No newline at end of file
+@test SciMLBase.successful_retcode(sol)

From 8cbd6631e1e61c76422090cce9be78b43ba7a48c Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 14:19:39 -0400
Subject: [PATCH 02/10] Mistake in rebase

---
 Project.toml   |  2 +-
 src/default.jl | 52 --------------------------------------------------
 2 files changed, 1 insertion(+), 53 deletions(-)

diff --git a/Project.toml b/Project.toml
index 1205910b7..7ee0625b0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -50,7 +50,7 @@ PrecompileTools = "1"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1"
 SciMLBase = "2.4"
-SimpleNonlinearSolve = "0.1"
+SimpleNonlinearSolve = "0.1.22"
 SparseDiffTools = "2.6"
 StaticArraysCore = "1.4"
 UnPack = "1.0"
diff --git a/src/default.jl b/src/default.jl
index 36ac42225..c7092341b 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -171,58 +171,6 @@ function SciMLBase.__solve(prob::NonlinearProblem{uType, iip},
         sols[idx].retcode, sols[idx].stats, original = sols[idx])
 end
 
-function SciMLBase.__solve(prob::NonlinearProblem{uType, true}, alg::FastShortcutNonlinearPolyalg, args...;
-    kwargs...) where {uType}
-
-    adkwargs = alg.adkwargs
-    linsolve = alg.linsolve
-    precs = alg.precs
-
-    sol1 = SciMLBase.__solve(prob, NewtonRaphson(;linsolve, precs, adkwargs...), args...; kwargs...)
-    if SciMLBase.successful_retcode(sol1)
-        return SciMLBase.build_solution(prob, alg, sol1.u, sol1.resid;
-                                        sol1.retcode, sol1.stats)
-    end
-
-    sol2 = SciMLBase.__solve(prob, NewtonRaphson(;linsolve, precs, linesearch=BackTracking(), adkwargs...), args...; kwargs...)
-    if SciMLBase.successful_retcode(sol2)
-        return SciMLBase.build_solution(prob, alg, sol2.u, sol2.resid;
-                                        sol2.retcode, sol2.stats)
-    end
-
-    sol3 = SciMLBase.__solve(prob, TrustRegion(;linsolve, precs, adkwargs...), args...; kwargs...)
-    if SciMLBase.successful_retcode(sol3)
-        return SciMLBase.build_solution(prob, alg, sol3.u, sol3.resid;
-                                        sol3.retcode, sol3.stats)
-    end
-
-    sol4 = SciMLBase.__solve(prob,  TrustRegion(;linsolve, precs, radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...), args...; kwargs...)
-    if SciMLBase.successful_retcode(sol4)
-        return SciMLBase.build_solution(prob, alg, sol4.u, sol4.resid;
-                                        sol4.retcode, sol4.stats)
-    end
-
-    resids = (sol1.resid, sol2.resid, sol3.resid, sol4.resid)
-    minfu, idx = findmin(DEFAULT_NORM, resids)
-
-    if idx == 1
-        SciMLBase.build_solution(prob, alg, sol1.u, sol1.resid;
-                                        sol1.retcode, sol1.stats)
-    elseif idx == 2
-        SciMLBase.build_solution(prob, alg, sol2.u, sol2.resid;
-                                        sol2.retcode, sol2.stats)
-    elseif idx == 3
-        SciMLBase.build_solution(prob, alg, sol3.u, sol3.resid;
-                                sol3.retcode, sol3.stats)
-    elseif idx == 4
-        SciMLBase.build_solution(prob, alg, sol4.u, sol4.resid;
-                                sol4.retcode, sol4.stats)
-    else
-        error("Unreachable reached, 박정석")
-    end
-
-end
-
 ## General shared polyalg functions
 
 function perform_step!(cache::Union{RobustMultiNewtonCache,

From 93bc347df448a7f9fe6e7a724d1ee7e1b7920b83 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 14:49:14 -0400
Subject: [PATCH 03/10] Make solve for FastShortcutNonlinearPolyalg type stable

---
 src/default.jl | 66 +++++++++++++++++++++++++++++++++-----------------
 1 file changed, 44 insertions(+), 22 deletions(-)

diff --git a/src/default.jl b/src/default.jl
index c7092341b..1c5dabe9d 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -65,6 +65,8 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::RobustMultiNe
         TrustRegion(; linsolve, precs,
             radius_update_scheme = RadiusUpdateSchemes.Fan, adkwargs...))
 
+    # Partially Type Unstable but can't do much since some upstream caches -- LineSearches
+    # and SparseDiffTools cause the instability
     return RobustMultiNewtonCache{iip}(map(solver -> SciMLBase.__init(prob, solver, args...;
                 kwargs...), algs), alg, 1)
 end
@@ -139,36 +141,56 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip},
 end
 
 # This version doesn't allocate all the caches!
-function SciMLBase.__solve(prob::NonlinearProblem{uType, iip},
+@generated function SciMLBase.__solve(prob::NonlinearProblem{uType, iip},
     alg::FastShortcutNonlinearPolyalg, args...; kwargs...) where {uType, iip}
-    @unpack adkwargs, linsolve, precs = alg
+    calls = [:(@unpack adkwargs, linsolve, precs = alg)]
 
     algs = [
-        iip ? Klement() : nothing, # Klement not yet implemented for IIP
-        iip ? Broyden() : nothing, # Broyden not yet implemented for IIP
-        NewtonRaphson(; linsolve, precs, adkwargs...),
-        NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...),
-        TrustRegion(; linsolve, precs, adkwargs...),
-        TrustRegion(; linsolve, precs,
-            radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...),
+        !iip ? :(Klement()) : nothing, # Klement not yet implemented for IIP
+        !iip ? :(Broyden()) : nothing, # Broyden not yet implemented for IIP
+        :(NewtonRaphson(; linsolve, precs, adkwargs...)),
+        :(NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...)),
+        :(TrustRegion(; linsolve, precs, adkwargs...)),
+        :(TrustRegion(; linsolve, precs,
+            radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...)),
     ]
     filter!(!isnothing, algs)
-
-    sols = Vector{SciMLBase.NonlinearSolution}(undef, length(algs))
-
-    for (i, solver) in enumerate(algs)
-        sols[i] = SciMLBase.__solve(prob, solver, args...; kwargs...)
-        if SciMLBase.successful_retcode(sols[i])
-            return SciMLBase.build_solution(prob, alg, sols[i].u, sols[i].resid;
-                sols[i].retcode, sols[i].stats, original = sols[i])
-        end
+    counter = 1
+    sol_syms = [gensym("sol") for i in 1:length(algs)]
+    for i in 1:length(algs)
+        cur_sol = sol_syms[i]
+        push!(calls,
+            quote
+                $(cur_sol) = SciMLBase.__solve(prob, $(algs[i]), args...; kwargs...)
+                if SciMLBase.successful_retcode($(cur_sol))
+                    return SciMLBase.build_solution(prob, alg, $(cur_sol).u,
+                        $(cur_sol).resid; $(cur_sol).retcode, $(cur_sol).stats,
+                        original = $(cur_sol))
+                end
+            end)
     end
 
-    resids = map(Base.Fix2(getproperty, resid), sols)
-    minfu, idx = findmin(DEFAULT_NORM, resids)
+    resids = map(x -> "$x.resid", sol_syms)
+
+    push!(calls,
+        quote
+            resids = $(Tuple(resids))
+            minfu, idx = findmin(DEFAULT_NORM, resids)
+        end)
+
+    for i in 1:length(algs)
+        push!(calls,
+            quote
+                if idx == $i
+                    return SciMLBase.build_solution(prob, alg, $(sol_syms[i]).u,
+                        $(sol_syms[i]).resid; $(sol_syms[i]).retcode, $(sol_syms[i]).stats,
+                        original = $(sol_syms[i]))
+                end
+            end)
+    end
+    push!(calls, :(error("Current choices shouldn't get here!")))
 
-    return SciMLBase.build_solution(prob, alg, sols[idx].u, sols[idx].resid;
-        sols[idx].retcode, sols[idx].stats, original = sols[idx])
+    return Expr(:block, calls...)
 end
 
 ## General shared polyalg functions

From aa03930aeaaf51aae0fd06dbf853124b2570222c Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 15:04:10 -0400
Subject: [PATCH 04/10] Fast Path for Robust MultiNewton solve

---
 src/default.jl   | 35 ++++++++++++++++++++++++-----------
 test/polyalgs.jl |  6 +++---
 2 files changed, 27 insertions(+), 14 deletions(-)

diff --git a/src/default.jl b/src/default.jl
index 1c5dabe9d..803c09ac2 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -142,20 +142,33 @@ end
 
 # This version doesn't allocate all the caches!
 @generated function SciMLBase.__solve(prob::NonlinearProblem{uType, iip},
-    alg::FastShortcutNonlinearPolyalg, args...; kwargs...) where {uType, iip}
+    alg::Union{FastShortcutNonlinearPolyalg, RobustMultiNewton}, args...;
+    kwargs...) where {uType, iip}
     calls = [:(@unpack adkwargs, linsolve, precs = alg)]
 
-    algs = [
-        !iip ? :(Klement()) : nothing, # Klement not yet implemented for IIP
-        !iip ? :(Broyden()) : nothing, # Broyden not yet implemented for IIP
-        :(NewtonRaphson(; linsolve, precs, adkwargs...)),
-        :(NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...)),
-        :(TrustRegion(; linsolve, precs, adkwargs...)),
-        :(TrustRegion(; linsolve, precs,
-            radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...)),
-    ]
+    algs = if parameterless_type(alg) == RobustMultiNewton
+        [
+            :(TrustRegion(; linsolve, precs, adkwargs...)),
+            :(TrustRegion(; linsolve, precs,
+                radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...)),
+            :(NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...)),
+            :(TrustRegion(; linsolve, precs,
+                radius_update_scheme = RadiusUpdateSchemes.NLsolve, adkwargs...)),
+            :(TrustRegion(; linsolve, precs,
+                radius_update_scheme = RadiusUpdateSchemes.Fan, adkwargs...)),
+        ]
+    else
+        [
+            !iip ? :(Klement()) : nothing, # Klement not yet implemented for IIP
+            !iip ? :(Broyden()) : nothing, # Broyden not yet implemented for IIP
+            :(NewtonRaphson(; linsolve, precs, adkwargs...)),
+            :(NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...)),
+            :(TrustRegion(; linsolve, precs, adkwargs...)),
+            :(TrustRegion(; linsolve, precs,
+                radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...)),
+        ]
+    end
     filter!(!isnothing, algs)
-    counter = 1
     sol_syms = [gensym("sol") for i in 1:length(algs)]
     for i in 1:length(algs)
         cur_sol = sol_syms[i]
diff --git a/test/polyalgs.jl b/test/polyalgs.jl
index 3aa552b24..cc3f7e7e4 100644
--- a/test/polyalgs.jl
+++ b/test/polyalgs.jl
@@ -3,9 +3,9 @@ using NonlinearSolve, Test
 f(u, p) = u .* u .- 2
 u0 = [1.0, 1.0]
 probN = NonlinearProblem(f, u0)
-@time solver = solve(probN, abstol = 1e-9)
-@time solver = solve(probN, RobustMultiNewton(), abstol = 1e-9)
-@time solver = solve(probN, FastShortcutNonlinearPolyalg(), abstol = 1e-9)
+@time solver = solve(probN; abstol = 1e-9)
+@time solver = solve(probN, RobustMultiNewton(); abstol = 1e-9)
+@time solver = solve(probN, FastShortcutNonlinearPolyalg(); abstol = 1e-9)
 
 # https://github.com/SciML/NonlinearSolve.jl/issues/153
 function f(du, u, p)

From 9319332fd813d247be32624e3435faa0dd2a7efa Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 15:17:01 -0400
Subject: [PATCH 05/10] Add tests for the caching interface

---
 src/default.jl   | 12 +++++++-----
 test/polyalgs.jl | 16 ++++++++++++++++
 2 files changed, 23 insertions(+), 5 deletions(-)

diff --git a/src/default.jl b/src/default.jl
index 803c09ac2..b7c16d735 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -46,7 +46,7 @@ function RobustMultiNewton(; concrete_jac = nothing, linsolve = nothing,
     return RobustMultiNewton{_unwrap_val(concrete_jac)}(adkwargs, linsolve, precs)
 end
 
-@concrete mutable struct RobustMultiNewtonCache{iip} <: AbstractNonlinearSolveCache{iip}
+@concrete mutable struct RobustMultiNewtonCache{iip, N} <: AbstractNonlinearSolveCache{iip}
     caches
     alg
     current::Int
@@ -67,8 +67,8 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::RobustMultiNe
 
     # Partially Type Unstable but can't do much since some upstream caches -- LineSearches
     # and SparseDiffTools cause the instability
-    return RobustMultiNewtonCache{iip}(map(solver -> SciMLBase.__init(prob, solver, args...;
-                kwargs...), algs), alg, 1)
+    return RobustMultiNewtonCache{iip, length(algs)}(map(solver -> SciMLBase.__init(prob,
+                solver, args...; kwargs...), algs), alg, 1)
 end
 
 """
@@ -110,7 +110,7 @@ function FastShortcutNonlinearPolyalg(; concrete_jac = nothing, linsolve = nothi
         precs)
 end
 
-@concrete mutable struct FastShortcutNonlinearPolyalgCache{iip} <:
+@concrete mutable struct FastShortcutNonlinearPolyalgCache{iip, N} <:
                          AbstractNonlinearSolveCache{iip}
     caches
     alg
@@ -136,7 +136,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip},
         TrustRegion(; linsolve, precs,
             radius_update_scheme = RadiusUpdateSchemes.Bastin, adkwargs...))
 
-    return FastShortcutNonlinearPolyalgCache{iip}(map(solver -> SciMLBase.__init(prob,
+    return FastShortcutNonlinearPolyalgCache{iip, length(algs)}(map(solver -> SciMLBase.__init(prob,
                 solver, args...; kwargs...), algs), alg, 1)
 end
 
@@ -159,6 +159,8 @@ end
         ]
     else
         [
+            # FIXME: Broyden and Klement are type unstable
+            #        (upstream SimpleNonlinearSolve.jl issue)
             !iip ? :(Klement()) : nothing, # Klement not yet implemented for IIP
             !iip ? :(Broyden()) : nothing, # Broyden not yet implemented for IIP
             :(NewtonRaphson(; linsolve, precs, adkwargs...)),
diff --git a/test/polyalgs.jl b/test/polyalgs.jl
index cc3f7e7e4..00b8ee82b 100644
--- a/test/polyalgs.jl
+++ b/test/polyalgs.jl
@@ -3,9 +3,25 @@ using NonlinearSolve, Test
 f(u, p) = u .* u .- 2
 u0 = [1.0, 1.0]
 probN = NonlinearProblem(f, u0)
+
+# Uses the `__solve` function
 @time solver = solve(probN; abstol = 1e-9)
+@test SciMLBase.successful_retcode(solver)
 @time solver = solve(probN, RobustMultiNewton(); abstol = 1e-9)
+@test SciMLBase.successful_retcode(solver)
 @time solver = solve(probN, FastShortcutNonlinearPolyalg(); abstol = 1e-9)
+@test SciMLBase.successful_retcode(solver)
+
+# Test the caching interface
+cache = init(probN; abstol = 1e-9)
+@time solver = solve!(cache)
+@test SciMLBase.successful_retcode(solver)
+cache = init(probN, RobustMultiNewton(); abstol = 1e-9)
+@time solver = solve!(cache)
+@test SciMLBase.successful_retcode(solver)
+cache = init(probN, FastShortcutNonlinearPolyalg(); abstol = 1e-9)
+@time solver = solve!(cache)
+@test SciMLBase.successful_retcode(solver)
 
 # https://github.com/SciML/NonlinearSolve.jl/issues/153
 function f(du, u, p)

From bb26efeb943edfd845b1c3c935c8f6e967f91a0f Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 15:56:10 -0400
Subject: [PATCH 06/10] Select the correct AD Type based on problem
 specification

---
 src/gaussnewton.jl | 10 +++++--
 src/levenberg.jl   | 12 +++++++--
 src/raphson.jl     | 10 +++++--
 src/trustRegion.jl | 13 ++++++++--
 src/utils.jl       | 65 ++++++++++++++++++++++++++++++++++++----------
 test/polyalgs.jl   |  8 +++---
 6 files changed, 92 insertions(+), 26 deletions(-)

diff --git a/src/gaussnewton.jl b/src/gaussnewton.jl
index 42521189e..6dcb0f835 100644
--- a/src/gaussnewton.jl
+++ b/src/gaussnewton.jl
@@ -14,7 +14,8 @@ for large-scale and numerically-difficult nonlinear least squares problems.
 
   - `autodiff`: determines the backend used for the Jacobian. Note that this argument is
     ignored if an analytical Jacobian is passed, as that will be used instead. Defaults to
-    `AutoForwardDiff()`. Valid choices are types from ADTypes.jl.
+    `nothing` which means that a default is selected according to the problem specification!
+    Valid choices are types from ADTypes.jl.
   - `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used,
     then the Jacobian will not be constructed and instead direct Jacobian-vector products
     `J*v` are computed using forward-mode automatic differentiation or finite differencing
@@ -41,6 +42,10 @@ for large-scale and numerically-difficult nonlinear least squares problems.
     precs
 end
 
+function set_ad(alg::GaussNewton{CJ}, ad) where {CJ}
+    return GaussNewton{CJ}(ad, alg.linsolve, alg.precs)
+end
+
 function GaussNewton(; concrete_jac = nothing, linsolve = CholeskyFactorization(),
     precs = DEFAULT_PRECS, adkwargs...)
     ad = default_adargs_to_adtype(; adkwargs...)
@@ -71,9 +76,10 @@ end
     stats::NLStats
 end
 
-function SciMLBase.__init(prob::NonlinearLeastSquaresProblem{uType, iip}, alg::GaussNewton,
+function SciMLBase.__init(prob::NonlinearLeastSquaresProblem{uType, iip}, alg_::GaussNewton,
     args...; alias_u0 = false, maxiters = 1000, abstol = 1e-6, internalnorm = DEFAULT_NORM,
     kwargs...) where {uType, iip}
+    alg = get_concrete_algorithm(alg_, prob)
     @unpack f, u0, p = prob
     u = alias_u0 ? u0 : deepcopy(u0)
     if iip
diff --git a/src/levenberg.jl b/src/levenberg.jl
index 07a8b8251..f35f35cb2 100644
--- a/src/levenberg.jl
+++ b/src/levenberg.jl
@@ -14,7 +14,8 @@ numerically-difficult nonlinear systems.
 
   - `autodiff`: determines the backend used for the Jacobian. Note that this argument is
     ignored if an analytical Jacobian is passed, as that will be used instead. Defaults to
-    `AutoForwardDiff()`. Valid choices are types from ADTypes.jl.
+    `nothing` which means that a default is selected according to the problem specification!
+    Valid choices are types from ADTypes.jl.
   - `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used,
     then the Jacobian will not be constructed and instead direct Jacobian-vector products
     `J*v` are computed using forward-mode automatic differentiation or finite differencing
@@ -86,6 +87,12 @@ numerically-difficult nonlinear systems.
     min_damping_D::T
 end
 
+function set_ad(alg::LevenbergMarquardt{CJ}, ad) where {CJ}
+    return LevenbergMarquardt{CJ}(ad, alg.linsolve, alg.precs, alg.damping_initial,
+        alg.damping_increase_factor, alg.damping_decrease_factor,
+        alg.finite_diff_step_geodesic, alg.α_geodesic, alg.b_uphill, alg.min_damping_D)
+end
+
 function LevenbergMarquardt(; concrete_jac = nothing, linsolve = nothing,
     precs = DEFAULT_PRECS, damping_initial::Real = 1.0, damping_increase_factor::Real = 2.0,
     damping_decrease_factor::Real = 3.0, finite_diff_step_geodesic::Real = 0.1,
@@ -141,9 +148,10 @@ end
 end
 
 function SciMLBase.__init(prob::Union{NonlinearProblem{uType, iip},
-        NonlinearLeastSquaresProblem{uType, iip}}, alg::LevenbergMarquardt,
+        NonlinearLeastSquaresProblem{uType, iip}}, alg_::LevenbergMarquardt,
     args...; alias_u0 = false, maxiters = 1000, abstol = 1e-6, internalnorm = DEFAULT_NORM,
     linsolve_kwargs = (;), kwargs...) where {uType, iip}
+    alg = get_concrete_algorithm(alg_, prob)
     @unpack f, u0, p = prob
     u = alias_u0 ? u0 : deepcopy(u0)
     fu1 = evaluate_f(prob, u)
diff --git a/src/raphson.jl b/src/raphson.jl
index a1fa1a1c5..642c7ee36 100644
--- a/src/raphson.jl
+++ b/src/raphson.jl
@@ -10,7 +10,8 @@ for large-scale and numerically-difficult nonlinear systems.
 
   - `autodiff`: determines the backend used for the Jacobian. Note that this argument is
     ignored if an analytical Jacobian is passed, as that will be used instead. Defaults to
-    `AutoForwardDiff()`. Valid choices are types from ADTypes.jl.
+    `nothing` which means that a default is selected according to the problem specification!
+    Valid choices are types from ADTypes.jl.
   - `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used,
     then the Jacobian will not be constructed and instead direct Jacobian-vector products
     `J*v` are computed using forward-mode automatic differentiation or finite differencing
@@ -36,6 +37,10 @@ for large-scale and numerically-difficult nonlinear systems.
     linesearch
 end
 
+function set_ad(alg::NewtonRaphson{CJ}, ad) where {CJ}
+    return NewtonRaphson{CJ}(ad, alg.linsolve, alg.precs, alg.linesearch)
+end
+
 function NewtonRaphson(; concrete_jac = nothing, linsolve = nothing,
     linesearch = LineSearch(), precs = DEFAULT_PRECS, adkwargs...)
     ad = default_adargs_to_adtype(; adkwargs...)
@@ -65,9 +70,10 @@ end
     lscache
 end
 
-function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::NewtonRaphson, args...;
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::NewtonRaphson, args...;
     alias_u0 = false, maxiters = 1000, abstol = 1e-6, internalnorm = DEFAULT_NORM,
     linsolve_kwargs = (;), kwargs...) where {uType, iip}
+    alg = get_concrete_algorithm(alg_, prob)
     @unpack f, u0, p = prob
     u = alias_u0 ? u0 : deepcopy(u0)
     fu1 = evaluate_f(prob, u)
diff --git a/src/trustRegion.jl b/src/trustRegion.jl
index 583b1b675..769f5e75c 100644
--- a/src/trustRegion.jl
+++ b/src/trustRegion.jl
@@ -97,7 +97,8 @@ for large-scale and numerically-difficult nonlinear systems.
 
   - `autodiff`: determines the backend used for the Jacobian. Note that this argument is
     ignored if an analytical Jacobian is passed, as that will be used instead. Defaults to
-    `AutoForwardDiff()`. Valid choices are types from ADTypes.jl.
+    `nothing` which means that a default is selected according to the problem specification!.
+    Valid choices are types from ADTypes.jl.
   - `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used,
     then the Jacobian will not be constructed and instead direct Jacobian-vector products
     `J*v` are computed using forward-mode automatic differentiation or finite differencing
@@ -162,6 +163,13 @@ for large-scale and numerically-difficult nonlinear systems.
     max_shrink_times::Int
 end
 
+function set_ad(alg::TrustRegion{CJ}, ad) where {CJ}
+    return TrustRegion{CJ}(ad, alg.linsolve, alg.precs, alg.radius_update_scheme,
+        alg.max_trust_radius, alg.initial_trust_radius, alg.step_threshold,
+        alg.shrink_threshold, alg.expand_threshold, alg.shrink_factor, alg.expand_factor,
+        alg.max_shrink_times)
+end
+
 function TrustRegion(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
     radius_update_scheme::RadiusUpdateSchemes.T = RadiusUpdateSchemes.Simple, #defaults to conventional radius update
     max_trust_radius::Real = 0 // 1, initial_trust_radius::Real = 0 // 1,
@@ -222,9 +230,10 @@ end
     stats::NLStats
 end
 
-function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, args...;
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::TrustRegion, args...;
     alias_u0 = false, maxiters = 1000, abstol = 1e-8, internalnorm = DEFAULT_NORM,
     linsolve_kwargs = (;), kwargs...) where {uType, iip}
+    alg = get_concrete_algorithm(alg_, prob)
     @unpack f, u0, p = prob
     u = alias_u0 ? u0 : deepcopy(u0)
     u_prev = zero(u)
diff --git a/src/utils.jl b/src/utils.jl
index 6dd1eeb75..173c279be 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -13,28 +13,41 @@ end
 @inline DEFAULT_NORM(u::AbstractArray) = sqrt(real(sum(UNITLESS_ABS2, u)) / length(u))
 @inline DEFAULT_NORM(u) = norm(u)
 
-alg_autodiff(alg::AbstractNewtonAlgorithm{<:AbstractFiniteDifferencesMode}) = false
-alg_autodiff(alg::AbstractNewtonAlgorithm) = true
-alg_autodiff(alg) = false
-
 """
     default_adargs_to_adtype(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), diff_type = Val{:forward})
 
 Construct the AD type from the arguments. This is mostly needed for compatibility with older
 code.
+
+!!! warning
+    `chunk_size`, `standardtag`, `diff_type`, and `autodiff::Union{Val, Bool}` are
+    deprecated and will be removed in v3. Update your code to directly specify
+    `autodiff=<ADTypes>`.
 """
-function default_adargs_to_adtype(; chunk_size = Val{0}(), autodiff = Val{true}(),
-    standardtag = Val{true}(), diff_type = Val{:forward}())
-    ad = _unwrap_val(autodiff)
-    # Old API
-    if ad isa Bool
-        # FIXME: standardtag is not the Tag
-        ad && return AutoForwardDiff(; chunksize = _unwrap_val(chunk_size),
-            tag = _unwrap_val(standardtag))
-        return AutoFiniteDiff(; fdtype = diff_type)
+function default_adargs_to_adtype(; chunk_size = missing, autodiff = nothing,
+    standardtag = missing, diff_type = missing)
+    # If using the new API short circuit
+    autodiff === nothing && return nothing
+    autodiff isa ADTypes.AbstractADType && return autodiff
+
+    # Deprecate the old version
+    if chunk_size !== missing || standardtag !== missing || diff_type !== missing ||
+       autodiff !== missing
+        Base.depwarn("`chunk_size`, `standardtag`, `diff_type`, \
+            `autodiff::Union{Val, Bool}` kwargs have been deprecated and will be removed in\
+             v3. Update your code to directly specify autodiff=<ADTypes>",
+            :default_adargs_to_adtype)
     end
-    return ad
+    chunk_size === missing && (chunk_size = Val{0}())
+    standardtag === missing && (standardtag = Val{true}())
+    diff_type === missing && (diff_type = Val{:forward}())
+    autodiff === missing && (autodiff = Val{true}())
+
+    ad = _unwrap_val(autodiff)
+    tag = _unwrap_val(standardtag)
+    ad && return AutoForwardDiff{_unwrap_val(chunk_size), typeof(tag)}(tag)
+    return AutoFiniteDiff(; fdtype = diff_type)
 end
 
 """
@@ -171,3 +184,27 @@ Defaults to `mul!(C, A, B)`. However, for sparse matrices uses `C .= A * B`.
 """
 __matmul!(C, A, B) = mul!(C, A, B)
 __matmul!(C::AbstractSparseMatrix, A, B) = C .= A * B
+
+# Concretize Algorithms
+function get_concrete_algorithm(alg, prob)
+    !hasfield(typeof(alg), :ad) && return alg
+    alg.ad isa ADTypes.AbstractADType && return alg
+
+    # Figure out the default AD
+    # Now that we have handed trivial cases, we can allow extending this function
+    # for specific algorithms
+    return __get_concrete_algorithm(alg, prob)
+end
+
+function __get_concrete_algorithm(alg, prob)
+    @unpack sparsity, jac_prototype = prob.f
+    use_sparse_ad = sparsity !== nothing || jac_prototype !== nothing
+    ad = if eltype(prob.u0) <: Complex
+        # Use Finite Differencing
+        use_sparse_ad ? AutoSparseFiniteDiff() : AutoFiniteDiff()
+    else
+        use_sparse_ad ? AutoSparseForwardDiff() :
+        AutoForwardDiff{ForwardDiff.pickchunksize(length(prob.u0)), Nothing}(nothing)
+    end
+    return set_ad(alg, ad)
+end
diff --git a/test/polyalgs.jl b/test/polyalgs.jl
index 00b8ee82b..4497eae97 100644
--- a/test/polyalgs.jl
+++ b/test/polyalgs.jl
@@ -2,7 +2,7 @@ using NonlinearSolve, Test
 
 f(u, p) = u .* u .- 2
 u0 = [1.0, 1.0]
-probN = NonlinearProblem(f, u0)
+probN = NonlinearProblem{false}(f, u0)
 
 # Uses the `__solve` function
 @time solver = solve(probN; abstol = 1e-9)
@@ -13,13 +13,13 @@ probN = NonlinearProblem(f, u0)
 @test SciMLBase.successful_retcode(solver)
 
 # Test the caching interface
-cache = init(probN; abstol = 1e-9)
+cache = init(probN; abstol = 1e-9);
 @time solver = solve!(cache)
 @test SciMLBase.successful_retcode(solver)
-cache = init(probN, RobustMultiNewton(); abstol = 1e-9)
+cache = init(probN, RobustMultiNewton(); abstol = 1e-9);
 @time solver = solve!(cache)
 @test SciMLBase.successful_retcode(solver)
-cache = init(probN, FastShortcutNonlinearPolyalg(); abstol = 1e-9)
+cache = init(probN, FastShortcutNonlinearPolyalg(); abstol = 1e-9);
 @time solver = solve!(cache)
 @test SciMLBase.successful_retcode(solver)
 

From cb5448a1a47f2edd22c25ba3e35c9a7e733e2d46 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 16:00:06 -0400
Subject: [PATCH 07/10] Update the tutorial

---
 docs/src/tutorials/large_systems.md | 28 +---------------------------
 1 file changed, 1 insertion(+), 27 deletions(-)

diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md
index 43c7f6cd6..5cbb99167 100644
--- a/docs/src/tutorials/large_systems.md
+++ b/docs/src/tutorials/large_systems.md
@@ -209,8 +209,7 @@ function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
 end
 
 @btime solve(prob_brusselator_2d_sparse,
-    NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = incompletelu,
-        concrete_jac = true));
+    NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true));
 nothing # hide
 ```
 
@@ -275,28 +274,3 @@ nothing # hide
 
 For more information on the preconditioner interface, see the
 [linear solver documentation](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/).
-
-## Speed up Jacobian computation with sparsity exploitation and matrix coloring
-
-To cut down the of Jacobian building overhead, we can choose to exploit the sparsity pattern and deploy matrix coloring during Jacobian construction. With NonlinearSolve.jl, we can simply use `autodiff=AutoSparseForwardDiff()` to automatically exploit the sparsity pattern of Jacobian matrices:
-
-```@example ill_conditioned_nlprob
-@btime solve(prob_brusselator_2d,
-    NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true,
-        autodiff = AutoSparseForwardDiff()));
-nothing # hide
-```
-
-To setup matrix coloring for the jacobian sparsity pattern, we can simply get the coloring vector by using [ArrayInterface.jl](https://github.com/JuliaArrays/ArrayInterface.jl) for the sparsity pattern of `jac_prototype`:
-
-```@example ill_conditioned_nlprob
-using ArrayInterface
-colorvec = ArrayInterface.matrix_colors(jac_sparsity)
-ff = NonlinearFunction(brusselator_2d_loop; jac_prototype = float.(jac_sparsity), colorvec)
-prob_brusselator_2d_sparse = NonlinearProblem(ff, u0, p)
-
-@btime solve(prob_brusselator_2d_sparse,
-    NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true,
-        autodiff = AutoSparseForwardDiff()));
-nothing # hide
-```

From d40b90173039209ed2afcc6fe7e27074a89c7bc6 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 16:05:06 -0400
Subject: [PATCH 08/10] make DFSane conform to the current code style

---
 src/dfsane.jl      | 175 ++++++++++++++++-----------------------------
 test/basictests.jl |  63 ++++++++--------
 2 files changed, 97 insertions(+), 141 deletions(-)

diff --git a/src/dfsane.jl b/src/dfsane.jl
index a77703906..a4bf30c8a 100644
--- a/src/dfsane.jl
+++ b/src/dfsane.jl
@@ -57,97 +57,51 @@ struct DFSane{T, F} <: AbstractNonlinearSolveAlgorithm
     max_inner_iterations::Int
 end
 
-function DFSane(; σ_min = 1e-10,
-                σ_max = 1e+10,
-                σ_1 = 1.0,
-                M = 10,
-                γ = 1e-4,
-                τ_min = 0.1,
-                τ_max = 0.5,
-                n_exp = 2,
-                η_strategy = (fn_1, n, x_n, f_n) -> fn_1 / n^2,
-                max_inner_iterations = 1000)
-    return DFSane{typeof(σ_min), typeof(η_strategy)}(σ_min,
-                                                     σ_max,
-                                                     σ_1,
-                                                     M,
-                                                     γ,
-                                                     τ_min,
-                                                     τ_max,
-                                                     n_exp,
-                                                     η_strategy,
-                                                     max_inner_iterations)
+function DFSane(; σ_min = 1e-10, σ_max = 1e+10, σ_1 = 1.0, M = 10, γ = 1e-4, τ_min = 0.1,
+    τ_max = 0.5, n_exp = 2, η_strategy = (fn_1, n, x_n, f_n) -> fn_1 / n^2,
+    max_inner_iterations = 1000)
+    return DFSane{typeof(σ_min), typeof(η_strategy)}(σ_min, σ_max, σ_1, M, γ, τ_min, τ_max,
+        n_exp, η_strategy, max_inner_iterations)
 end
-mutable struct DFSaneCache{iip, algType, uType, resType, T, pType,
-                           INType,
-                           tolType,
-                           probType}
-    f::Function
-    alg::algType
-    uₙ::uType
-    uₙ₋₁::uType
-    fuₙ::resType
-    fuₙ₋₁::resType
-    𝒹::uType
-    ℋ::Vector{T}
-    f₍ₙₒᵣₘ₎ₙ₋₁::T
-    f₍ₙₒᵣₘ₎₀::T
-    M::Int
-    σₙ::T
-    σₘᵢₙ::T
-    σₘₐₓ::T
-    α₁::T
-    γ::T
-    τₘᵢₙ::T
-    τₘₐₓ::T
+
+@concrete mutable struct DFSaneCache{iip}
+    alg
+    uₙ
+    uₙ₋₁
+    fuₙ
+    fuₙ₋₁
+    𝒹
+    ℋ
+    f₍ₙₒᵣₘ₎ₙ₋₁
+    f₍ₙₒᵣₘ₎₀
+    M
+    σₙ
+    σₘᵢₙ
+    σₘₐₓ
+    α₁
+    γ
+    τₘᵢₙ
+    τₘₐₓ
     nₑₓₚ::Int
-    p::pType
+    p
     force_stop::Bool
     maxiters::Int
-    internalnorm::INType
+    internalnorm
     retcode::SciMLBase.ReturnCode.T
-    abstol::tolType
-    prob::probType
+    abstol
+    prob
     stats::NLStats
-    function DFSaneCache{iip}(f::Function, alg::algType, uₙ::uType, uₙ₋₁::uType,
-                              fuₙ::resType, fuₙ₋₁::resType, 𝒹::uType, ℋ::Vector{T},
-                              f₍ₙₒᵣₘ₎ₙ₋₁::T, f₍ₙₒᵣₘ₎₀::T, M::Int, σₙ::T, σₘᵢₙ::T, σₘₐₓ::T,
-                              α₁::T, γ::T, τₘᵢₙ::T, τₘₐₓ::T, nₑₓₚ::Int, p::pType,
-                              force_stop::Bool, maxiters::Int, internalnorm::INType,
-                              retcode::SciMLBase.ReturnCode.T, abstol::tolType,
-                              prob::probType,
-                              stats::NLStats) where {iip, algType, uType,
-                                                     resType, T, pType, INType,
-                                                     tolType,
-                                                     probType
-                                                     }
-        new{iip, algType, uType, resType, T, pType, INType, tolType,
-            probType
-            }(f, alg, uₙ, uₙ₋₁, fuₙ, fuₙ₋₁, 𝒹, ℋ, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀, M, σₙ,
-              σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ,
-              τₘₐₓ, nₑₓₚ, p, force_stop, maxiters, internalnorm,
-              retcode,
-              abstol, prob, stats)
-    end
 end
 
-function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::DFSane,
-                          args...;
-                          alias_u0 = false,
-                          maxiters = 1000,
-                          abstol = 1e-6,
-                          internalnorm = DEFAULT_NORM,
-                          kwargs...) where {uType, iip}
-    if alias_u0
-        uₙ = prob.u0
-    else
-        uₙ = deepcopy(prob.u0)
-    end
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::DFSane, args...;
+    alias_u0 = false, maxiters = 1000, abstol = 1e-6, internalnorm = DEFAULT_NORM,
+    kwargs...) where {uType, iip}
+    uₙ = alias_u0 ? prob.u0 : deepcopy(prob.u0)
 
     p = prob.p
     T = eltype(uₙ)
     σₘᵢₙ, σₘₐₓ, γ, τₘᵢₙ, τₘₐₓ = T(alg.σ_min), T(alg.σ_max), T(alg.γ), T(alg.τ_min),
-                                T(alg.τ_max)
+    T(alg.τ_max)
     α₁ = one(T)
     γ = T(alg.γ)
     f₍ₙₒᵣₘ₎ₙ₋₁ = α₁
@@ -157,27 +111,27 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::DFSane,
     𝒹, uₙ₋₁, fuₙ, fuₙ₋₁ = copy(uₙ), copy(uₙ), copy(uₙ), copy(uₙ)
 
     if iip
-        f = (dx, x) -> prob.f(dx, x, p)
-        f(fuₙ₋₁, uₙ₋₁)
+        # f = (dx, x) -> prob.f(dx, x, p)
+        # f(fuₙ₋₁, uₙ₋₁)
+        prob.f(fuₙ₋₁, uₙ₋₁, p)
     else
-        f = (x) -> prob.f(x, p)
-        fuₙ₋₁ = f(uₙ₋₁)
+        # f = (x) -> prob.f(x, p)
+        fuₙ₋₁ = prob.f(uₙ₋₁, p) # f(uₙ₋₁)
     end
 
     f₍ₙₒᵣₘ₎ₙ₋₁ = norm(fuₙ₋₁)^nₑₓₚ
     f₍ₙₒᵣₘ₎₀ = f₍ₙₒᵣₘ₎ₙ₋₁
 
     ℋ = fill(f₍ₙₒᵣₘ₎ₙ₋₁, M)
-    return DFSaneCache{iip}(f, alg, uₙ, uₙ₋₁, fuₙ, fuₙ₋₁, 𝒹, ℋ, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀,
-                            M, σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ,
-                            τₘₐₓ, nₑₓₚ, p, false, maxiters,
-                            internalnorm, ReturnCode.Default, abstol, prob,
-                            NLStats(1, 0, 0, 0, 0))
+    return DFSaneCache{iip}(alg, uₙ, uₙ₋₁, fuₙ, fuₙ₋₁, 𝒹, ℋ, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀,
+        M, σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, p, false, maxiters,
+        internalnorm, ReturnCode.Default, abstol, prob, NLStats(1, 0, 0, 0, 0))
 end
 
 function perform_step!(cache::DFSaneCache{true})
-    @unpack f, alg, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀,
-    σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, M = cache
+    @unpack alg, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀, σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, M = cache
+
+    f = iip ? (dx, x) -> cache.prob.f(dx, x, cache.p) : (x) -> cache.prob.f(x, cache.p)
 
     T = eltype(cache.uₙ)
     n = cache.stats.nsteps
@@ -202,10 +156,9 @@ function perform_step!(cache::DFSaneCache{true})
 
         f₍ₙₒᵣₘ₎ₙ ≤ 𝒸 && break
 
-        α₊ = clamp(α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ /
-                   (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₊ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
-                   τₘᵢₙ * α₊,
-                   τₘₐₓ * α₊)
+        α₊ = clamp(α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₊ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
+            τₘᵢₙ * α₊,
+            τₘₐₓ * α₊)
         @. cache.uₙ = cache.uₙ₋₁ - α₋ * cache.𝒹
 
         f(cache.fuₙ, cache.uₙ)
@@ -214,8 +167,8 @@ function perform_step!(cache::DFSaneCache{true})
         f₍ₙₒᵣₘ₎ₙ .≤ 𝒸 && break
 
         α₋ = clamp(α₋^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₋ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
-                   τₘᵢₙ * α₋,
-                   τₘₐₓ * α₋)
+            τₘᵢₙ * α₋,
+            τₘₐₓ * α₋)
 
         @. cache.uₙ = cache.uₙ₋₁ + α₊ * cache.𝒹
         f(cache.fuₙ, cache.uₙ)
@@ -253,8 +206,9 @@ function perform_step!(cache::DFSaneCache{true})
 end
 
 function perform_step!(cache::DFSaneCache{false})
-    @unpack f, alg, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀,
-    σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, M = cache
+    @unpack alg, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀, σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, M = cache
+
+    f = iip ? (dx, x) -> cache.prob.f(dx, x, cache.p) : (x) -> cache.prob.f(x, cache.p)
 
     T = eltype(cache.uₙ)
     n = cache.stats.nsteps
@@ -279,10 +233,8 @@ function perform_step!(cache::DFSaneCache{false})
 
         f₍ₙₒᵣₘ₎ₙ ≤ 𝒸 && break
 
-        α₊ = clamp(α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ /
-                   (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₊ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
-                   τₘᵢₙ * α₊,
-                   τₘₐₓ * α₊)
+        α₊ = clamp(α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₊ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
+            τₘᵢₙ * α₊, τₘₐₓ * α₊)
         cache.uₙ = @. cache.uₙ₋₁ - α₋ * cache.𝒹
 
         cache.fuₙ = f(cache.uₙ)
@@ -291,8 +243,7 @@ function perform_step!(cache::DFSaneCache{false})
         f₍ₙₒᵣₘ₎ₙ .≤ 𝒸 && break
 
         α₋ = clamp(α₋^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₋ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
-                   τₘᵢₙ * α₋,
-                   τₘₐₓ * α₋)
+            τₘᵢₙ * α₋, τₘₐₓ * α₋)
 
         cache.uₙ = @. cache.uₙ₋₁ + α₊ * cache.𝒹
         cache.fuₙ = f(cache.uₙ)
@@ -341,25 +292,23 @@ function SciMLBase.solve!(cache::DFSaneCache)
         cache.retcode = ReturnCode.Success
     end
 
-    SciMLBase.build_solution(cache.prob, cache.alg, cache.uₙ, cache.fuₙ;
-                             retcode = cache.retcode, stats = cache.stats)
+    return SciMLBase.build_solution(cache.prob, cache.alg, cache.uₙ, cache.fuₙ;
+        retcode = cache.retcode, stats = cache.stats)
 end
 
 function SciMLBase.reinit!(cache::DFSaneCache{iip}, u0 = cache.uₙ; p = cache.p,
-                           abstol = cache.abstol, maxiters = cache.maxiters) where {iip}
+    abstol = cache.abstol, maxiters = cache.maxiters) where {iip}
     cache.p = p
     if iip
         recursivecopy!(cache.uₙ, u0)
         recursivecopy!(cache.uₙ₋₁, u0)
-        cache.f = (dx, x) -> cache.prob.f(dx, x, p)
-        cache.f(cache.fuₙ, cache.uₙ)
-        cache.f(cache.fuₙ₋₁, cache.uₙ)
+        cache.prob.f(cache.fuₙ, cache.uₙ, p)
+        cache.prob.f(cache.fuₙ₋₁, cache.uₙ, p)
     else
         cache.uₙ = u0
         cache.uₙ₋₁ = u0
-        cache.f = (x) -> cache.prob.f(x, p)
-        cache.fuₙ = cache.f(cache.uₙ)
-        cache.fuₙ₋₁ = cache.f(cache.uₙ)
+        cache.fuₙ = cache.prob.f(cache.uₙ, p)
+        cache.fuₙ₋₁ = cache.prob.f(cache.uₙ, p)
     end
 
     cache.f₍ₙₒᵣₘ₎ₙ₋₁ = norm(cache.fuₙ₋₁)^cache.nₑₓₚ
diff --git a/test/basictests.jl b/test/basictests.jl
index 51756cf97..efa148bb0 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -391,18 +391,17 @@ end
     end
 end
 
-
 # --- DFSane tests ---
 
 @testset "DFSane" begin
-    function benchmark_nlsolve_oop(f, u0, p=2.0)
+    function benchmark_nlsolve_oop(f, u0, p = 2.0)
         prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, DFSane(), abstol=1e-9)
+        return solve(prob, DFSane(), abstol = 1e-9)
     end
 
-    function benchmark_nlsolve_iip(f, u0, p=2.0)
+    function benchmark_nlsolve_iip(f, u0, p = 2.0)
         prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, DFSane(), abstol=1e-9)
+        return solve(prob, DFSane(), abstol = 1e-9)
     end
 
     u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
@@ -413,7 +412,7 @@ end
         @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
         cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), DFSane(),
-            abstol=1e-9)
+            abstol = 1e-9)
         @test (@ballocated solve!($cache)) < 200
     end
 
@@ -424,11 +423,10 @@ end
         @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
         cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0),
-            DFSane(), abstol=1e-9)
+            DFSane(), abstol = 1e-9)
         @test (@ballocated solve!($cache)) ≤ 64
     end
 
-
     @testset "[OOP] [Immutable AD]" begin
         broken_forwarddiff = [1.6, 2.9, 3.0, 3.5, 4.0, 81.0]
         for p in 1.1:0.1:100.0
@@ -437,14 +435,14 @@ end
             if any(x -> isnan(x) || x <= 1e-5 || x >= 1e5, res)
                 @test_broken all(res .≈ sqrt(p))
                 @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
-                        @SVector[1.0, 1.0], p).u[end], p)) ≈ 1 / (2 * sqrt(p))
+                    @SVector[1.0, 1.0], p).u[end], p)) ≈ 1 / (2 * sqrt(p))
             elseif p in broken_forwarddiff
                 @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
-                        @SVector[1.0, 1.0], p).u[end], p)) ≈ 1 / (2 * sqrt(p))
+                    @SVector[1.0, 1.0], p).u[end], p)) ≈ 1 / (2 * sqrt(p))
             else
                 @test all(res .≈ sqrt(p))
                 @test isapprox(abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
-                            @SVector[1.0, 1.0], p).u[end], p)), 1 / (2 * sqrt(p)))
+                        @SVector[1.0, 1.0], p).u[end], p)), 1 / (2 * sqrt(p)))
             end
         end
     end
@@ -456,12 +454,22 @@ end
 
             if any(x -> isnan(x) || x <= 1e-5 || x >= 1e5, res)
                 @test_broken res ≈ sqrt(p)
-                @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, 1.0, p).u, p)) ≈ 1 / (2 * sqrt(p))
+                @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                        1.0,
+                        p).u,
+                    p)) ≈ 1 / (2 * sqrt(p))
             elseif p in broken_forwarddiff
-                @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, 1.0, p).u, p)) ≈ 1 / (2 * sqrt(p))
+                @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                        1.0,
+                        p).u,
+                    p)) ≈ 1 / (2 * sqrt(p))
             else
                 @test res ≈ sqrt(p)
-                @test isapprox(abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, 1.0, p).u, p)), 1 / (2 * sqrt(p)))
+                @test isapprox(abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                            1.0,
+                            p).u,
+                        p)),
+                    1 / (2 * sqrt(p)))
             end
         end
     end
@@ -475,20 +483,19 @@ end
     # Iterator interface
     function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
         probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
-        cache = init(probN, DFSane(); maxiters=100, abstol=1e-10)
+        cache = init(probN, DFSane(); maxiters = 100, abstol = 1e-10)
         sols = zeros(length(p_range))
         for (i, p) in enumerate(p_range)
-            reinit!(cache, iip ? [cache.uₙ[1]] : cache.uₙ; p=p)
+            reinit!(cache, iip ? [cache.uₙ[1]] : cache.uₙ; p = p)
             sol = solve!(cache)
             sols[i] = iip ? sol.u[1] : sol.u
         end
         return sols
     end
-    p = range(0.01, 2, length=200)
+    p = range(0.01, 2, length = 200)
     @test abs.(nlprob_iterator_interface(quadratic_f, p, Val(false))) ≈ sqrt.(p)
     @test abs.(nlprob_iterator_interface(quadratic_f!, p, Val(true))) ≈ sqrt.(p)
 
-
     # Test that `DFSane` passes a test that `NewtonRaphson` fails on.
     @testset "Newton Raphson Fails" begin
         u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
@@ -518,18 +525,18 @@ end
             η_strategy)
         for options in list_of_options
             local probN, sol, alg
-            alg = DFSane(σ_min=options[1],
-                σ_max=options[2],
-                σ_1=options[3],
-                M=options[4],
-                γ=options[5],
-                τ_min=options[6],
-                τ_max=options[7],
-                n_exp=options[8],
-                η_strategy=options[9])
+            alg = DFSane(σ_min = options[1],
+                σ_max = options[2],
+                σ_1 = options[3],
+                M = options[4],
+                γ = options[5],
+                τ_min = options[6],
+                τ_max = options[7],
+                n_exp = options[8],
+                η_strategy = options[9])
 
             probN = NonlinearProblem{false}(quadratic_f, [1.0, 1.0], 2.0)
-            sol = solve(probN, alg, abstol=1e-11)
+            sol = solve(probN, alg, abstol = 1e-11)
             println(abs.(quadratic_f(sol.u, 2.0)))
             @test all(abs.(quadratic_f(sol.u, 2.0)) .< 1e-10)
         end

From 3724fbc9779a9411a4b0a8fcb39ed2ebf2f59b51 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 18:21:13 -0400
Subject: [PATCH 09/10] Make it generated

---
 src/default.jl | 93 ++++++++++++++++++++++++++------------------------
 src/dfsane.jl  |  4 +--
 2 files changed, 50 insertions(+), 47 deletions(-)

diff --git a/src/default.jl b/src/default.jl
index b7c16d735..913f994fe 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -185,11 +185,14 @@ end
             end)
     end
 
-    resids = map(x -> "$x.resid", sol_syms)
+    resids = map(x -> Symbol("$(x)_resid"), sol_syms)
+    for (sym, resid) in zip(sol_syms, resids)
+        push!(calls, :($(resid) = $(sym).resid))
+    end
 
     push!(calls,
         quote
-            resids = $(Tuple(resids))
+            resids = tuple($(Tuple(resids)...))
             minfu, idx = findmin(DEFAULT_NORM, resids)
         end)
 
@@ -198,8 +201,7 @@ end
             quote
                 if idx == $i
                     return SciMLBase.build_solution(prob, alg, $(sol_syms[i]).u,
-                        $(sol_syms[i]).resid; $(sol_syms[i]).retcode, $(sol_syms[i]).stats,
-                        original = $(sol_syms[i]))
+                        $(sol_syms[i]).resid; $(sol_syms[i]).retcode, $(sol_syms[i]).stats)
                 end
             end)
     end
@@ -210,53 +212,54 @@ end
 
 ## General shared polyalg functions
 
-function perform_step!(cache::Union{RobustMultiNewtonCache,
-    FastShortcutNonlinearPolyalgCache})
-    current = cache.current
-    1 ≤ current ≤ length(cache.caches) || error("Current choices shouldn't get here!")
-
-    current_cache = cache.caches[current]
-    while not_terminated(current_cache)
-        perform_step!(current_cache)
+@generated function SciMLBase.solve!(cache::Union{RobustMultiNewtonCache{iip, N},
+    FastShortcutNonlinearPolyalgCache{iip, N}}) where {iip, N}
+    calls = [
+        quote
+            1 ≤ cache.current ≤ length(cache.caches) ||
+                error("Current choices shouldn't get here!")
+        end,
+    ]
+
+    cache_syms = [gensym("cache") for i in 1:N]
+    sol_syms = [gensym("sol") for i in 1:N]
+    for i in 1:N
+        push!(calls,
+            quote
+                $(cache_syms[i]) = cache.caches[$(i)]
+                if $(i) == cache.current
+                    $(sol_syms[i]) = SciMLBase.solve!($(cache_syms[i]))
+                    if SciMLBase.successful_retcode($(sol_syms[i]))
+                        stats = $(sol_syms[i]).stats
+                        u = $(sol_syms[i]).u
+                        fu = get_fu($(cache_syms[i]))
+                        return SciMLBase.build_solution($(sol_syms[i]).prob, cache.alg, u,
+                            fu; retcode = ReturnCode.Success, stats,
+                            original = $(sol_syms[i]))
+                    end
+                    cache.current = $(i + 1)
+                end
+            end)
     end
 
-    return nothing
-end
-
-function SciMLBase.solve!(cache::Union{RobustMultiNewtonCache,
-    FastShortcutNonlinearPolyalgCache})
-    current = cache.current
-    1 ≤ current ≤ length(cache.caches) || error("Current choices shouldn't get here!")
-
-    current_cache = cache.caches[current]
-    while current ≤ length(cache.caches) #  && !all(terminated[current:end])
-        sol_tmp = solve!(current_cache)
-        SciMLBase.successful_retcode(sol_tmp) && break
-        current += 1
-        cache.current = current
-        current_cache = cache.caches[current]
+    resids = map(x -> Symbol("$(x)_resid"), cache_syms)
+    for (sym, resid) in zip(cache_syms, resids)
+        push!(calls, :($(resid) = get_fu($(sym))))
     end
+    push!(calls,
+        quote
+            retcode = ReturnCode.MaxIters
 
-    if current ≤ length(cache.caches)
-        retcode = ReturnCode.Success
-
-        stats = cache.caches[current].stats
-        u = cache.caches[current].u
-        fu = get_fu(cache.caches[current])
-
-        return SciMLBase.build_solution(cache.caches[1].prob, cache.alg, u, fu;
-            retcode, stats)
-    else
-        retcode = ReturnCode.MaxIters
+            fus = tuple($(Tuple(resids)...))
+            minfu, idx = findmin(cache.caches[1].internalnorm, fus)
+            stats = cache.caches[idx].stats
+            u = cache.caches[idx].u
 
-        fus = get_fu.(cache.caches)
-        minfu, idx = findmin(cache.caches[1].internalnorm, fus)
-        stats = cache.caches[idx].stats
-        u = cache.caches[idx].u
+            return SciMLBase.build_solution(cache.caches[idx].prob, cache.alg, u,
+                fus[idx]; retcode, stats)
+        end)
 
-        return SciMLBase.build_solution(cache.caches[idx].prob, cache.alg, u, fus[idx];
-            retcode, stats)
-    end
+    return Expr(:block, calls...)
 end
 
 function SciMLBase.reinit!(cache::Union{RobustMultiNewtonCache,
diff --git a/src/dfsane.jl b/src/dfsane.jl
index a4bf30c8a..13de5ff6a 100644
--- a/src/dfsane.jl
+++ b/src/dfsane.jl
@@ -131,7 +131,7 @@ end
 function perform_step!(cache::DFSaneCache{true})
     @unpack alg, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀, σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, M = cache
 
-    f = iip ? (dx, x) -> cache.prob.f(dx, x, cache.p) : (x) -> cache.prob.f(x, cache.p)
+    f = (dx, x) -> cache.prob.f(dx, x, cache.p)
 
     T = eltype(cache.uₙ)
     n = cache.stats.nsteps
@@ -208,7 +208,7 @@ end
 function perform_step!(cache::DFSaneCache{false})
     @unpack alg, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀, σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, M = cache
 
-    f = iip ? (dx, x) -> cache.prob.f(dx, x, cache.p) : (x) -> cache.prob.f(x, cache.p)
+    f = x -> cache.prob.f(x, cache.p)
 
     T = eltype(cache.uₙ)
     n = cache.stats.nsteps

From 91c0ca09c7c78672b5fbca905f2ef5491c057305 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 17 Oct 2023 19:57:33 -0400
Subject: [PATCH 10/10] Rebase

---
 docs/pages.jl                                |  3 +-
 docs/src/tutorials/code_optimization.md      |  2 +-
 docs/src/tutorials/getting_started.md        | 40 ++++++++++----------
 docs/src/tutorials/iterator_interface.md     |  2 +-
 docs/src/tutorials/small_compile.md          |  2 +-
 docs/src/tutorials/termination_conditions.md |  2 +-
 test/basictests.jl                           |  3 +-
 7 files changed, 27 insertions(+), 27 deletions(-)

diff --git a/docs/pages.jl b/docs/pages.jl
index 35f924cdd..c29d2874b 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -2,8 +2,7 @@
 
 pages = ["index.md",
     "tutorials/getting_started.md",
-    "Tutorials" => Any[
-        "tutorials/code_optimization.md",
+    "Tutorials" => Any["tutorials/code_optimization.md",
         "tutorials/large_systems.md",
         "tutorials/small_compile.md",
         "tutorials/termination_conditions.md",
diff --git a/docs/src/tutorials/code_optimization.md b/docs/src/tutorials/code_optimization.md
index bb6fbdbd4..24d003206 100644
--- a/docs/src/tutorials/code_optimization.md
+++ b/docs/src/tutorials/code_optimization.md
@@ -55,4 +55,4 @@ prob = NonlinearProblem(f!, u0, p)
 
 linsolve = LinearSolve.KrylovJL_GMRES()
 sol = solve(prob, NewtonRaphson(; linsolve), reltol = 1e-9)
-```
\ No newline at end of file
+```
diff --git a/docs/src/tutorials/getting_started.md b/docs/src/tutorials/getting_started.md
index edcbb6921..9f1685557 100644
--- a/docs/src/tutorials/getting_started.md
+++ b/docs/src/tutorials/getting_started.md
@@ -9,19 +9,19 @@ to understanding the deeper parts of the documentation.
 
 There are three types of nonlinear systems:
 
-1. The "standard nonlinear system", i.e. the `NonlinearProblem`. This is a
-   system of equations with an initial condition where you want to satisfy
-   all equations simultaniously.
-2. The "interval rootfinding problem", i.e. the `IntervalNonlinearProblem`.
-   This is the case where you're given an interval `[a,b]` and need to find
-   where `f(u) = 0` for `u` inside the bounds.
-3. The "steady state problem", i.e. find the `u` such that `u' = f(u) = 0`.
-   While related to (1), it's not entirely the same because there's a uniquely
-   defined privledged root.
-4. The nonlinear least squares problem, which is an overconstrained nonlinear
-   system (i.e. more equations than states) which might not be satisfiable, i.e.
-   there may be no `u` such that `f(u) = 0`, and thus we find the `u` which
-   minimizes `||f(u)||` in the least squares sense.
+ 1. The "standard nonlinear system", i.e. the `NonlinearProblem`. This is a
+    system of equations with an initial condition where you want to satisfy
+    all equations simultaniously.
+ 2. The "interval rootfinding problem", i.e. the `IntervalNonlinearProblem`.
+    This is the case where you're given an interval `[a,b]` and need to find
+    where `f(u) = 0` for `u` inside the bounds.
+ 3. The "steady state problem", i.e. find the `u` such that `u' = f(u) = 0`.
+    While related to (1), it's not entirely the same because there's a uniquely
+    defined privledged root.
+ 4. The nonlinear least squares problem, which is an overconstrained nonlinear
+    system (i.e. more equations than states) which might not be satisfiable, i.e.
+    there may be no `u` such that `f(u) = 0`, and thus we find the `u` which
+    minimizes `||f(u)||` in the least squares sense.
 
 For now let's focus on the first two. The other two are covered in later tutorials,
 but from the first two we can show the general flow of the NonlinearSolve.jl package.
@@ -105,7 +105,7 @@ For a complete list of solver choices, see [the nonlinear system solvers page](@
 Next we can modify the tolerances. Here let's set some really low tolerances to force a tight solution:
 
 ```@example 1
-solve(prob, TrustRegion(), reltol=1e-12, abstol=1e-12)
+solve(prob, TrustRegion(), reltol = 1e-12, abstol = 1e-12)
 ```
 
 There are many more options for doing this configuring. Specifically for handling termination conditions,
@@ -139,10 +139,10 @@ sol = solve(prob_int, ITP(), abstol = 0.01)
 Congrats, you now know how to use the basics of NonlinearSolve.jl! However, there is so much more to
 see. Next check out:
 
-- [Some code optimization tricks to know about with NonlinearSolve.jl](@ref code_optimization)
-- [An iterator interface which lets you step through the solving process step by step](@ref iterator)
-- [How to handle large systems of equations efficiently](@ref large_systems)
-- [Ways to use NonlinearSolve.jl that is faster to startup and can statically compile](@ref fast_startup)
-- [More detailed termination conditions](@ref termination_conditions_tutorial)
+  - [Some code optimization tricks to know about with NonlinearSolve.jl](@ref code_optimization)
+  - [An iterator interface which lets you step through the solving process step by step](@ref iterator)
+  - [How to handle large systems of equations efficiently](@ref large_systems)
+  - [Ways to use NonlinearSolve.jl that is faster to startup and can statically compile](@ref fast_startup)
+  - [More detailed termination conditions](@ref termination_conditions_tutorial)
 
-And also check out the rest of the manual.
\ No newline at end of file
+And also check out the rest of the manual.
diff --git a/docs/src/tutorials/iterator_interface.md b/docs/src/tutorials/iterator_interface.md
index 640fb5f02..c0fb914f4 100644
--- a/docs/src/tutorials/iterator_interface.md
+++ b/docs/src/tutorials/iterator_interface.md
@@ -1,7 +1,7 @@
 # [Nonlinear Solver Iterator Interface](@id iterator)
 
 !!! warn
-
+    
     This iterator interface will be expanded with a `step!` function soon!
 
 There is an iterator form of the nonlinear solver which mirrors the DiffEq integrator interface:
diff --git a/docs/src/tutorials/small_compile.md b/docs/src/tutorials/small_compile.md
index e6aa07312..aab4531b8 100644
--- a/docs/src/tutorials/small_compile.md
+++ b/docs/src/tutorials/small_compile.md
@@ -1,3 +1,3 @@
 # Faster Startup and and Static Compilation
 
-This is a stub article to be completed soon.
\ No newline at end of file
+This is a stub article to be completed soon.
diff --git a/docs/src/tutorials/termination_conditions.md b/docs/src/tutorials/termination_conditions.md
index 152e0abc4..6a7e8a3cd 100644
--- a/docs/src/tutorials/termination_conditions.md
+++ b/docs/src/tutorials/termination_conditions.md
@@ -1,3 +1,3 @@
 # [More Detailed Termination Conditions](@id termination_conditions_tutorial)
 
-This is a stub article to be completed soon.
\ No newline at end of file
+This is a stub article to be completed soon.
diff --git a/test/basictests.jl b/test/basictests.jl
index efa148bb0..06cfc103d 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -225,7 +225,8 @@ end
         radius_update_scheme in radius_update_schemes
 
         probN = NonlinearProblem(quadratic_f, u0, 2.0)
-        @test all(solve(probN, TrustRegion(; autodiff, radius_update_scheme)).u .≈ sqrt(2.0))
+        @test all(solve(probN, TrustRegion(; autodiff, radius_update_scheme)).u .≈
+                  sqrt(2.0))
     end
 
     # Test that `TrustRegion` passes a test that `NewtonRaphson` fails on.