From 39f8271ea334c3c8659499fd788c3e3c54c60dac Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Sun, 19 Jan 2025 14:00:36 +0530
Subject: [PATCH 1/2] fix: only solve parameter initialization for
 `NonlinearSystem`

---
 src/systems/diffeqs/abstractodesystem.jl  |   2 +-
 src/systems/nonlinear/initializesystem.jl | 106 +++++-----
 src/systems/problem_utils.jl              |  38 +++-
 test/initializationsystem.jl              | 239 +++++++++-------------
 4 files changed, 180 insertions(+), 205 deletions(-)

diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 00b3c11ec3..ac28d41dd4 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -1335,7 +1335,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem,
 
     # TODO: throw on uninitialized arrays
     filter!(x -> !(x isa Symbolics.Arr), uninit)
-    if !isempty(uninit)
+    if is_time_dependent(sys) && !isempty(uninit)
         allow_incomplete || throw(IncompleteInitializationError(uninit))
         # for incomplete initialization, we will add the missing variables as parameters.
         # they will be updated by `update_initializeprob!` and `initializeprobmap` will
diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 0b556eef68..a7f98d79b1 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -42,51 +42,53 @@ function generate_initializesystem(sys::AbstractSystem;
         diffmap = Dict()
     end
 
-    if has_schedule(sys) && (schedule = get_schedule(sys); !isnothing(schedule))
-        # 2) process dummy derivatives and u0map into initialization system
-        # prepare map for dummy derivative substitution
-        for x in filter(x -> !isnothing(x[1]), schedule.dummy_sub)
-            # set dummy derivatives to default_dd_guess unless specified
-            push!(defs, x[1] => get(guesses, x[1], default_dd_guess))
-        end
-        function process_u0map_with_dummysubs(y, x)
-            y = get(schedule.dummy_sub, y, y)
-            y = fixpoint_sub(y, diffmap)
-            if y ∈ vars_set
-                # variables specified in u0 overrides defaults
-                push!(defs, y => x)
-            elseif y isa Symbolics.Arr
-                # TODO: don't scalarize arrays
-                merge!(defs, Dict(scalarize(y .=> x)))
-            elseif y isa Symbolics.BasicSymbolic
-                # y is a derivative expression expanded; add it to the initialization equations
-                push!(eqs_ics, y ~ x)
-            else
-                error("Initialization expression $y is currently not supported. If its a higher order derivative expression, then only the dummy derivative expressions are supported.")
+    if is_time_dependent(sys)
+        if has_schedule(sys) && (schedule = get_schedule(sys); !isnothing(schedule))
+            # 2) process dummy derivatives and u0map into initialization system
+            # prepare map for dummy derivative substitution
+            for x in filter(x -> !isnothing(x[1]), schedule.dummy_sub)
+                # set dummy derivatives to default_dd_guess unless specified
+                push!(defs, x[1] => get(guesses, x[1], default_dd_guess))
             end
-        end
-        for (y, x) in u0map
-            if Symbolics.isarraysymbolic(y)
-                process_u0map_with_dummysubs.(collect(y), collect(x))
-            else
-                process_u0map_with_dummysubs(y, x)
+            function process_u0map_with_dummysubs(y, x)
+                y = get(schedule.dummy_sub, y, y)
+                y = fixpoint_sub(y, diffmap)
+                if y ∈ vars_set
+                    # variables specified in u0 overrides defaults
+                    push!(defs, y => x)
+                elseif y isa Symbolics.Arr
+                    # TODO: don't scalarize arrays
+                    merge!(defs, Dict(scalarize(y .=> x)))
+                elseif y isa Symbolics.BasicSymbolic
+                    # y is a derivative expression expanded; add it to the initialization equations
+                    push!(eqs_ics, y ~ x)
+                else
+                    error("Initialization expression $y is currently not supported. If its a higher order derivative expression, then only the dummy derivative expressions are supported.")
+                end
+            end
+            for (y, x) in u0map
+                if Symbolics.isarraysymbolic(y)
+                    process_u0map_with_dummysubs.(collect(y), collect(x))
+                else
+                    process_u0map_with_dummysubs(y, x)
+                end
+            end
+        else
+            # 2) System doesn't have a schedule, so dummy derivatives don't exist/aren't handled (SDESystem)
+            for (k, v) in u0map
+                defs[k] = v
             end
         end
-    else
-        # 2) System doesn't have a schedule, so dummy derivatives don't exist/aren't handled (SDESystem)
-        for (k, v) in u0map
-            defs[k] = v
-        end
-    end
 
-    # 3) process other variables
-    for var in vars
-        if var ∈ keys(defs)
-            push!(eqs_ics, var ~ defs[var])
-        elseif var ∈ keys(guesses)
-            push!(defs, var => guesses[var])
-        elseif check_defguess
-            error("Invalid setup: variable $(var) has no default value or initial guess")
+        # 3) process other variables
+        for var in vars
+            if var ∈ keys(defs)
+                push!(eqs_ics, var ~ defs[var])
+            elseif var ∈ keys(guesses)
+                push!(defs, var => guesses[var])
+            elseif check_defguess
+                error("Invalid setup: variable $(var) has no default value or initial guess")
+            end
         end
     end
 
@@ -180,16 +182,24 @@ function generate_initializesystem(sys::AbstractSystem;
     pars = Vector{SymbolicParam}(filter(p -> !haskey(paramsubs, p), parameters(sys)))
     is_time_dependent(sys) && push!(pars, get_iv(sys))
 
-    # 8) use observed equations for guesses of observed variables if not provided
-    for eq in trueobs
-        haskey(defs, eq.lhs) && continue
-        any(x -> isequal(default_toterm(x), eq.lhs), keys(defs)) && continue
+    if is_time_dependent(sys)
+        # 8) use observed equations for guesses of observed variables if not provided
+        for eq in trueobs
+            haskey(defs, eq.lhs) && continue
+            any(x -> isequal(default_toterm(x), eq.lhs), keys(defs)) && continue
 
-        defs[eq.lhs] = eq.rhs
+            defs[eq.lhs] = eq.rhs
+        end
+        append!(eqs_ics, trueobs)
+    end
+
+    eqs_ics = Symbolics.substitute.(eqs_ics, (paramsubs,))
+    if is_time_dependent(sys)
+        vars = [vars; collect(values(paramsubs))]
+    else
+        vars = collect(values(paramsubs))
     end
 
-    eqs_ics = Symbolics.substitute.([eqs_ics; trueobs], (paramsubs,))
-    vars = [vars; collect(values(paramsubs))]
     for k in keys(defs)
         defs[k] = substitute(defs[k], paramsubs)
     end
diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index 4cbf495d87..a8d443349b 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -546,30 +546,46 @@ function maybe_build_initialization_problem(
         initializeprob = ModelingToolkit.InitializationProblem(
             sys, t, u0map, pmap; guesses, kwargs...)
 
-        all_init_syms = Set(all_symbols(initializeprob))
-        solved_unknowns = filter(var -> var in all_init_syms, unknowns(sys))
-        initializeprobmap = getu(initializeprob, solved_unknowns)
+        if is_time_dependent(sys)
+            all_init_syms = Set(all_symbols(initializeprob))
+            solved_unknowns = filter(var -> var in all_init_syms, unknowns(sys))
+            initializeprobmap = getu(initializeprob, solved_unknowns)
+        else
+            initializeprobmap = nothing
+        end
 
         punknowns = [p
                      for p in all_variable_symbols(initializeprob)
                      if is_parameter(sys, p)]
-        getpunknowns = getu(initializeprob, punknowns)
-        setpunknowns = setp(sys, punknowns)
-        initializeprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns)
+        if isempty(punknowns)
+            initializeprobpmap = nothing
+        else
+            getpunknowns = getu(initializeprob, punknowns)
+            setpunknowns = setp(sys, punknowns)
+            initializeprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns)
+        end
 
         reqd_syms = parameter_symbols(initializeprob)
-        update_initializeprob! = UpdateInitializeprob(
-            getu(sys, reqd_syms), setu(initializeprob, reqd_syms))
+        # we still want the `initialization_data` because it helps with `remake`
+        if initializeprobmap === nothing && initializeprobpmap === nothing
+            update_initializeprob! = nothing
+        else
+            update_initializeprob! = UpdateInitializeprob(
+                getu(sys, reqd_syms), setu(initializeprob, reqd_syms))
+        end
+
         for p in punknowns
             p = unwrap(p)
             stype = symtype(p)
             op[p] = get_temporary_value(p)
         end
 
-        for v in missing_unknowns
-            op[v] = zero_var(v)
+        if is_time_dependent(sys)
+            for v in missing_unknowns
+                op[v] = zero_var(v)
+            end
+            empty!(missing_unknowns)
         end
-        empty!(missing_unknowns)
         return (;
             initialization_data = SciMLBase.OverrideInitData(
                 initializeprob, update_initializeprob!, initializeprobmap,
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 62b6bba707..0a4ef18038 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -583,14 +583,6 @@ sol = solve(oprob_2nd_order_2, Rosenbrock23()) # retcode: Success
     @test all(sol(1.0, idxs = sys.x) .≈ +exp(1)) && all(sol(1.0, idxs = sys.y) .≈ -exp(1))
 end
 
-NonlinearSystemWrapper(eqs, t; kws...) = NonlinearSystem(eqs; kws...)
-function NonlinearProblemWrapper(sys, u0, tspan, args...; kwargs...)
-    NonlinearProblem(sys, u0, args...; kwargs...)
-end
-function NLLSProblemWrapper(sys, u0, tspan, args...; kwargs...)
-    NonlinearLeastSquaresProblem(sys, u0, args...; kwargs...)
-end
-
 @testset "Initialization of parameters" begin
     @variables _x(..) y(t)
     @parameters p q
@@ -606,31 +598,11 @@ end
         (ModelingToolkit.System, ODEProblem, Tsit5(), zeros(2)),
         (ModelingToolkit.System, SDEProblem, ImplicitEM(), [a, b]),
         (ModelingToolkit.System, DDEProblem, MethodOfSteps(Tsit5()), [_x(t - 0.1), 0.0]),
-        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), [_x(t - 0.1) + a, b]),
-        # polyalg cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper,
-            FastShortcutNonlinearPolyalg(), zeros(2)),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, NewtonRaphson(), zeros(2)),
-        # quasi newton cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, Klement(), zeros(2)),
-        # noinit cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, SimpleNewtonRaphson(), zeros(2)),
-        # DFSane cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, DFSane(), zeros(2)),
-        # Least squares
-        # polyalg cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, FastShortcutNLLSPolyalg(), zeros(2)),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, LevenbergMarquardt(), zeros(2)),
-        # noinit cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, SimpleGaussNewton(), zeros(2))
+        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), [_x(t - 0.1) + a, b])
     ]
-        is_nlsolve = alg isa SciMLBase.AbstractNonlinearAlgorithm
-
         function test_parameter(prob, sym, val, initialval = zero(val))
             @test prob.ps[sym] ≈ initialval
-            if !is_nlsolve || prob.u0 !== nothing
+            if prob.u0 !== nothing
                 @test init(prob, alg).ps[sym] ≈ val
             end
             @test solve(prob, alg).ps[sym] ≈ val
@@ -641,7 +613,6 @@ end
             @test is_variable(isys, p)
             @test equation in equations(isys) || (0 ~ -equation.rhs) in equations(isys)
         end
-        D = is_nlsolve ? v -> v^3 : Differential(t)
 
         u0map = Dict(x => 1.0, y => 1.0)
         pmap = Dict()
@@ -777,6 +748,90 @@ end
     @test solve(prob, Tsit5()).ps[spring.s_rel0] ≈ -3.905
 end
 
+@testset "NonlinearSystem initialization" begin
+    nl_algs = [FastShortcutNonlinearPolyalg(), NewtonRaphson(),
+        Klement(), SimpleNewtonRaphson(), DFSane()]
+    nlls_algs = [FastShortcutNLLSPolyalg(), LevenbergMarquardt(), SimpleGaussNewton()]
+
+    @testset "No initialization for variables" begin
+        @variables x=1.0 y=0.0 z=0.0
+        @parameters σ=10.0 ρ=26.0 β=8 / 3
+
+        eqs = [0 ~ σ * (y - x),
+            0 ~ x * (ρ - z) - y,
+            0 ~ x * y - β * z]
+        @mtkbuild ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
+
+        prob = NonlinearProblem(ns, [])
+        @test prob.f.initialization_data.update_initializeprob! === nothing
+        @test prob.f.initialization_data.initializeprobmap === nothing
+        @test prob.f.initialization_data.initializeprobpmap === nothing
+        for alg in nl_algs
+            @test SciMLBase.successful_retcode(solve(prob, alg))
+        end
+
+        prob = NonlinearLeastSquaresProblem(ns, [])
+        @test prob.f.initialization_data.update_initializeprob! === nothing
+        @test prob.f.initialization_data.initializeprobmap === nothing
+        @test prob.f.initialization_data.initializeprobpmap === nothing
+        for alg in nlls_algs
+            @test SciMLBase.successful_retcode(solve(prob, alg))
+        end
+    end
+
+    prob_alg_combinations = zip(
+        [NonlinearProblem, NonlinearLeastSquaresProblem], [nl_algs, nlls_algs])
+    @testset "Parameter initialization" begin
+        function test_parameter(prob, alg, param, val)
+            integ = init(prob, alg)
+            @test integ.ps[param]≈val rtol=1e-6
+            sol = solve(prob, alg)
+            @test sol.ps[param]≈val rtol=1e-6
+            @test SciMLBase.successful_retcode(sol)
+        end
+        @variables x=1.0 y=3.0
+        @parameters p q
+
+        @mtkbuild sys = NonlinearSystem(
+            [(x - p)^2 + (y - q)^3 ~ 0, x - q ~ 0]; defaults = [q => missing],
+            guesses = [q => 1.0], initialization_eqs = [p^2 + q^2 + 2p * q ~ 0])
+
+        for (probT, algs) in prob_alg_combinations
+            prob = probT(sys, [], [p => 2.0])
+            @test prob.f.initialization_data !== nothing
+            @test prob.f.initialization_data.initializeprobmap === nothing
+            for alg in algs
+                test_parameter(prob, alg, q, -2.0)
+            end
+
+            # `update_initializeprob!` works
+            prob.ps[p] = -2.0
+            for alg in algs
+                test_parameter(prob, alg, q, 2.0)
+            end
+            prob.ps[p] = 2.0
+
+            # `remake` works
+            prob2 = remake(prob; p = [p => -2.0])
+            @test prob2.f.initialization_data !== nothing
+            @test prob2.f.initialization_data.initializeprobmap === nothing
+            for alg in algs
+                test_parameter(prob2, alg, q, 2.0)
+            end
+
+            # changing types works
+            ps = parameter_values(prob)
+            newps = SciMLStructures.replace(Tunable(), ps, ForwardDiff.Dual.(ps.tunable))
+            prob3 = remake(prob; p = newps)
+            @test prob3.f.initialization_data !== nothing
+            @test eltype(state_values(prob3.f.initialization_data.initializeprob)) <:
+                  ForwardDiff.Dual
+            @test eltype(prob3.f.initialization_data.initializeprob.p.tunable) <:
+                  ForwardDiff.Dual
+        end
+    end
+end
+
 @testset "Update initializeprob parameters" begin
     @variables _x(..) y(t)
     @parameters p q
@@ -787,29 +842,8 @@ end
         (ModelingToolkit.System, ODEProblem, Tsit5(), zeros(2)),
         (ModelingToolkit.System, SDEProblem, ImplicitEM(), [a, b]),
         (ModelingToolkit.System, DDEProblem, MethodOfSteps(Tsit5()), [_x(t - 0.1), 0.0]),
-        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), [_x(t - 0.1) + a, b]),
-        # polyalg cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper,
-            FastShortcutNonlinearPolyalg(), zeros(2)),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, NewtonRaphson(), zeros(2)),
-        # quasi newton cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, Klement(), zeros(2)),
-        # noinit cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, SimpleNewtonRaphson(), zeros(2)),
-        # DFSane cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, DFSane(), zeros(2)),
-        # Least squares
-        # polyalg cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, FastShortcutNLLSPolyalg(), zeros(2)),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, LevenbergMarquardt(), zeros(2)),
-        # noinit cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, SimpleGaussNewton(), zeros(2))
+        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), [_x(t - 0.1) + a, b])
     ]
-        is_nlsolve = alg isa SciMLBase.AbstractNonlinearAlgorithm
-        D = is_nlsolve ? v -> v^3 : Differential(t)
-
         @mtkbuild sys = System(
             [D(x) ~ x + rhss[1], p ~ x + y + rhss[2]], t; guesses = [x => 0.0, p => 0.0])
         prob = Problem(sys, [y => 1.0], (0.0, 1.0), [p => 3.0])
@@ -837,29 +871,8 @@ end
         (ModelingToolkit.System, ODEProblem, Tsit5(), 0),
         (ModelingToolkit.System, SDEProblem, ImplicitEM(), a),
         (ModelingToolkit.System, DDEProblem, MethodOfSteps(Tsit5()), _x(t - 0.1)),
-        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), _x(t - 0.1) + a),
-        # polyalg cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper,
-            FastShortcutNonlinearPolyalg(), 0),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, NewtonRaphson(), 0),
-        # quasi newton cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, Klement(), 0),
-        # noinit cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, SimpleNewtonRaphson(), 0),
-        # DFSane cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, DFSane(), 0),
-        # Least squares
-        # polyalg cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, FastShortcutNLLSPolyalg(), 0),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, LevenbergMarquardt(), 0),
-        # noinit cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, SimpleGaussNewton(), 0)
+        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), _x(t - 0.1) + a)
     ]
-        is_nlsolve = alg isa SciMLBase.AbstractNonlinearAlgorithm
-        D = is_nlsolve ? v -> v^3 : Differential(t)
-
         @mtkbuild sys = System(
             [D(x) ~ 2x + r + rhss], t; parameter_dependencies = [r ~ p + 2q, q ~ p + 3],
             guesses = [p => 1.0])
@@ -881,29 +894,8 @@ end
         (ModelingToolkit.System, ODEProblem, Tsit5(), zeros(2)),
         (ModelingToolkit.System, SDEProblem, ImplicitEM(), [a, b]),
         (ModelingToolkit.System, DDEProblem, MethodOfSteps(Tsit5()), [_x(t - 0.1), 0.0]),
-        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), [_x(t - 0.1) + a, b]),
-        # polyalg cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper,
-            FastShortcutNonlinearPolyalg(), zeros(2)),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, NewtonRaphson(), zeros(2)),
-        # quasi newton cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, Klement(), zeros(2)),
-        # noinit cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, SimpleNewtonRaphson(), zeros(2)),
-        # DFSane cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, DFSane(), zeros(2)),
-        # Least squares
-        # polyalg cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, FastShortcutNLLSPolyalg(), zeros(2)),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, LevenbergMarquardt(), zeros(2)),
-        # noinit cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, SimpleGaussNewton(), zeros(2))
+        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), [_x(t - 0.1) + a, b])
     ]
-        is_nlsolve = alg isa SciMLBase.AbstractNonlinearAlgorithm
-        D = is_nlsolve ? v -> v^3 : Differential(t)
-
         @mtkbuild sys = System(
             [D(x) ~ x + rhss[1], p ~ x + y + rhss[2]], t; defaults = [p => missing], guesses = [
                 x => 0.0, p => 0.0])
@@ -934,35 +926,14 @@ end
         (ModelingToolkit.System, ODEProblem, Tsit5(), 0),
         (ModelingToolkit.System, SDEProblem, ImplicitEM(), a),
         (ModelingToolkit.System, DDEProblem, MethodOfSteps(Tsit5()), _x(t - 0.1)),
-        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), _x(t - 0.1) + a),
-        # polyalg cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper,
-            FastShortcutNonlinearPolyalg(), 0),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, NewtonRaphson(), 0),
-        # quasi newton cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, Klement(), 0),
-        # noinit cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, SimpleNewtonRaphson(), 0),
-        # DFSane cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, DFSane(), 0),
-        # Least squares
-        # polyalg cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, FastShortcutNLLSPolyalg(), 0),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, LevenbergMarquardt(), 0),
-        # noinit cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, SimpleGaussNewton(), 0)
+        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), _x(t - 0.1) + a)
     ]
-        is_nlsolve = alg isa SciMLBase.AbstractNonlinearAlgorithm
-        D = is_nlsolve ? v -> v^3 : Differential(t)
         alge_eqs = [y^2 * q + q^2 * x ~ 0, z * p - p^2 * x * z ~ 0]
 
         @mtkbuild sys = System(
             [D(x) ~ x * p + y^2 * q + rhss; alge_eqs],
             t; guesses = [x => 0.0, y => 0.0, z => 0.0, p => 0.0, q => 0.0])
-        prob = Problem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0, q => missing];
-            initialization_eqs = is_nlsolve ? alge_eqs : [])
+        prob = Problem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0, q => missing])
         @test is_variable(prob.f.initialization_data.initializeprob, q)
         ps = prob.p
         newps = SciMLStructures.replace(Tunable(), ps, ForwardDiff.Dual.(ps.tunable))
@@ -1010,40 +981,18 @@ end
         (ModelingToolkit.System, ODEProblem, Tsit5(), 0),
         (ModelingToolkit.System, SDEProblem, ImplicitEM(), a),
         (ModelingToolkit.System, DDEProblem, MethodOfSteps(Tsit5()), _x(t - 0.1)),
-        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), _x(t - 0.1) + a),
-        # polyalg cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper,
-            FastShortcutNonlinearPolyalg(), 0),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, NewtonRaphson(), 0),
-        # quasi newton cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, Klement(), 0),
-        # noinit cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, SimpleNewtonRaphson(), 0),
-        # DFSane cache
-        (NonlinearSystemWrapper, NonlinearProblemWrapper, DFSane(), 0),
-        # Least squares
-        # polyalg cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, FastShortcutNLLSPolyalg(), 0),
-        # generalized first order cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, LevenbergMarquardt(), 0),
-        # noinit cache
-        (NonlinearSystemWrapper, NLLSProblemWrapper, SimpleGaussNewton(), 0)
+        (ModelingToolkit.System, SDDEProblem, ImplicitEM(), _x(t - 0.1) + a)
     ]
-        is_nlsolve = alg isa SciMLBase.AbstractNonlinearAlgorithm
-        D = is_nlsolve ? v -> v^3 : Differential(t)
         alge_eqs = [y^2 + 4y * p^2 ~ x^3]
         @mtkbuild sys = System(
             [D(x) ~ x + p * y^2 + rhss; alge_eqs], t; guesses = [
                 y => 1.0, p => 1.0])
-        prob = Problem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0];
-            initialization_eqs = is_nlsolve ? alge_eqs : [])
+        prob = Problem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0])
         @test is_variable(prob.f.initialization_data.initializeprob, y)
         prob2 = @test_nowarn remake(prob; p = [p => 3.0]) # ensure no over/under-determined warning
         @test is_variable(prob.f.initialization_data.initializeprob, y)
 
-        prob = Problem(sys, [y => 1.0, x => 2.0], (0.0, 1.0), [p => missing];
-            initialization_eqs = is_nlsolve ? alge_eqs : [])
+        prob = Problem(sys, [y => 1.0, x => 2.0], (0.0, 1.0), [p => missing])
         @test is_variable(prob.f.initialization_data.initializeprob, p)
         prob2 = @test_nowarn remake(prob; u0 = [y => 0.5])
         @test is_variable(prob.f.initialization_data.initializeprob, p)

From ccfc13b12ea766e5e99e0e12d59a97f07dce221c Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 19 Jan 2025 11:25:21 +0100
Subject: [PATCH 2/2] Update Project.toml

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

diff --git a/Project.toml b/Project.toml
index 7ccb4d2110..1f46a4771c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -132,7 +132,7 @@ RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
 SCCNonlinearSolve = "1.0.0"
-SciMLBase = "2.71"
+SciMLBase = "2.71.1"
 SciMLStructures = "1.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"