From b2dbf919fa7a38bda16b45b0e1b98345e3596c49 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 6 Aug 2024 10:57:25 +0530
Subject: [PATCH 01/19] feat: allow parameters to be unknowns in the
 initialization system

---
 src/systems/abstractsystem.jl             |   9 ++
 src/systems/diffeqs/abstractodesystem.jl  | 102 +++++++++++++++++++--
 src/systems/nonlinear/initializesystem.jl | 106 ++++++++++++++++++++--
 src/systems/parameter_buffer.jl           |   3 +
 src/utils.jl                              |   4 +-
 5 files changed, 207 insertions(+), 17 deletions(-)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index aa02a5fb73..661a5fd4b4 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -736,6 +736,15 @@ function has_observed_with_lhs(sys, sym)
     end
 end
 
+function has_parameter_dependency_with_lhs(sys, sym)
+    has_parameter_dependencies(sys) || return false
+    if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
+        return any(isequal(sym), ic.dependent_pars)
+    else
+        return any(isequal(sym), [eq.lhs for eq in parameter_dependencies(sys)])
+    end
+end
+
 function _all_ts_idxs!(ts_idxs, ::NotSymbolic, sys, sym)
     if is_variable(sys, sym) || is_independent_variable(sys, sym)
         push!(ts_idxs, ContinuousTimeseries())
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index a7ae411def..623b2b47af 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -357,7 +357,9 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem,
         analytic = nothing,
         split_idxs = nothing,
         initializeprob = nothing,
+        update_initializeprob! = nothing,
         initializeprobmap = nothing,
+        initializeprobpmap = nothing,
         kwargs...) where {iip, specialize}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEFunction`")
@@ -459,7 +461,9 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem,
         sparsity = sparsity ? jacobian_sparsity(sys) : nothing,
         analytic = analytic,
         initializeprob = initializeprob,
-        initializeprobmap = initializeprobmap)
+        update_initializeprob! = update_initializeprob!,
+        initializeprobmap = initializeprobmap,
+        initializeprobpmap = initializeprobpmap)
 end
 
 """
@@ -789,6 +793,45 @@ function get_u0(
     return u0, defs
 end
 
+struct GetUpdatedMTKParameters{G, S}
+    # `getu` functor which gets parameters that are unknowns during initialization
+    getpunknowns::G
+    # `setu` functor which returns a modified MTKParameters using those parameters
+    setpunknowns::S
+end
+
+function (f::GetUpdatedMTKParameters)(prob, initializesol)
+    mtkp = copy(parameter_values(prob))
+    f.setpunknowns(mtkp, f.getpunknowns(initializesol))
+    mtkp
+end
+
+struct UpdateInitializeprob{G, S}
+    # `getu` functor which gets all values from prob
+    getvals::G
+    # `setu` functor which updates initializeprob with values
+    setvals::S
+end
+
+function (f::UpdateInitializeprob)(initializeprob, prob)
+    f.setvals(initializeprob, f.getvals(prob))
+end
+
+function get_temporary_value(p)
+    stype = symtype(unwrap(p))
+    return if stype == Real
+        zero(Float64)
+    elseif stype <: AbstractArray{Real}
+        zeros(Float64, size(p))
+    elseif stype <: Real
+        zero(stype)
+    elseif stype <: AbstractArray
+        zeros(eltype(stype), size(p))
+    else
+        error("Nonnumeric parameter $p with symtype $stype cannot be solved for during initialization")
+    end
+end
+
 function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
         implicit_dae = false, du0map = nothing,
         version = nothing, tgrad = false,
@@ -829,18 +872,38 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
     end
 
     if eltype(parammap) <: Pair
-        parammap = Dict(unwrap(k) => v for (k, v) in todict(parammap))
+        parammap = Dict{Any, Any}(unwrap(k) => v for (k, v) in parammap)
     elseif parammap isa AbstractArray
         if isempty(parammap)
             parammap = SciMLBase.NullParameters()
         else
-            parammap = Dict(unwrap.(parameters(sys)) .=> parammap)
+            parammap = Dict{Any, Any}(unwrap.(parameters(sys)) .=> parammap)
         end
     end
-
+    defs = defaults(sys)
+    if has_guesses(sys)
+        guesses = merge(
+            ModelingToolkit.guesses(sys), isempty(guesses) ? Dict() : todict(guesses))
+        solvablepars = [p
+                        for p in parameters(sys)
+                        if is_parameter_solvable(p, parammap, defs, guesses)]
+
+        pvarmap = if parammap === nothing || parammap == SciMLBase.NullParameters() || !(eltype(parammap) <: Pair) && isempty(parammap)
+            defs
+        else
+            merge(defs, todict(parammap))
+        end
+        setparobserved = filter(keys(pvarmap)) do var
+            has_parameter_dependency_with_lhs(sys, var)
+        end
+    else
+        solvablepars = ()
+        setparobserved = ()
+    end
     # ModelingToolkit.get_tearing_state(sys) !== nothing => Requires structural_simplify first
     if sys isa ODESystem && build_initializeprob &&
-       (((implicit_dae || !isempty(missingvars) || !isempty(setobserved)) &&
+       (((implicit_dae || !isempty(missingvars) || !isempty(solvablepars) ||
+          !isempty(setobserved) || !isempty(setparobserved)) &&
          ModelingToolkit.get_tearing_state(sys) !== nothing) ||
         !isempty(initialization_equations(sys))) && t !== nothing
         if eltype(u0map) <: Number
@@ -854,14 +917,32 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
             sys, t, u0map, parammap; guesses, warn_initialize_determined,
             initialization_eqs, eval_expression, eval_module, fully_determined, check_units)
         initializeprobmap = getu(initializeprob, unknowns(sys))
+        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)
+        reqd_syms = parameter_symbols(initializeprob)
+        update_initializeprob! = UpdateInitializeprob(
+            getu(sys, reqd_syms), setu(initializeprob, reqd_syms))
 
         zerovars = Dict(setdiff(unknowns(sys), keys(defaults(sys))) .=> 0.0)
+        if parammap isa SciMLBase.NullParameters
+            parammap = Dict()
+        end
+        for p in punknowns
+            p = unwrap(p)
+            stype = symtype(p)
+            parammap[p] = get_temporary_value(p)
+        end
         trueinit = collect(merge(zerovars, eltype(u0map) <: Pair ? todict(u0map) : u0map))
         u0map isa StaticArraysCore.StaticArray &&
             (trueinit = SVector{length(trueinit)}(trueinit))
     else
         initializeprob = nothing
+        update_initializeprob! = nothing
         initializeprobmap = nothing
+        initializeprobpmap = nothing
         trueinit = u0map
     end
 
@@ -909,7 +990,9 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
         sparse = sparse, eval_expression = eval_expression,
         eval_module = eval_module,
         initializeprob = initializeprob,
+        update_initializeprob! = update_initializeprob!,
         initializeprobmap = initializeprobmap,
+        initializeprobpmap = initializeprobpmap,
         kwargs...)
     implicit_dae ? (f, du0, u0, p) : (f, u0, p)
 end
@@ -1471,10 +1554,12 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
         isys = get_initializesystem(sys; initialization_eqs, check_units)
     elseif isempty(u0map) && get_initializesystem(sys) === nothing
         isys = structural_simplify(
-            generate_initializesystem(sys; initialization_eqs, check_units); fully_determined)
+            generate_initializesystem(
+                sys; initialization_eqs, check_units, pmap = parammap); fully_determined)
     else
         isys = structural_simplify(
-            generate_initializesystem(sys; u0map, initialization_eqs, check_units); fully_determined)
+            generate_initializesystem(
+                sys; u0map, initialization_eqs, check_units, pmap = parammap); fully_determined)
     end
 
     uninit = setdiff(unknowns(sys), [unknowns(isys); getfield.(observed(isys), :lhs)])
@@ -1498,6 +1583,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
     parammap = parammap isa DiffEqBase.NullParameters || isempty(parammap) ?
                [get_iv(sys) => t] :
                merge(todict(parammap), Dict(get_iv(sys) => t))
+    parammap = Dict(k => v for (k, v) in parammap if v !== missing)
     if isempty(u0map)
         u0map = Dict()
     end
@@ -1505,7 +1591,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
         guesses = Dict()
     end
 
-    u0map = merge(todict(guesses), todict(u0map))
+    u0map = merge(ModelingToolkit.guesses(sys), todict(guesses), todict(u0map))
     if neqs == nunknown
         NonlinearProblem(isys, u0map, parammap; kwargs...)
     else
diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 2097e2955e..d4253afda5 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -5,6 +5,7 @@ Generate `NonlinearSystem` which initializes an ODE problem from specified initi
 """
 function generate_initializesystem(sys::ODESystem;
         u0map = Dict(),
+        pmap = Dict(),
         initialization_eqs = [],
         guesses = Dict(),
         default_dd_guess = 0.0,
@@ -74,12 +75,103 @@ function generate_initializesystem(sys::ODESystem;
         end
     end
 
-    pars = [parameters(sys); get_iv(sys)] # include independent variable as pseudo-parameter
-    eqs_ics = [eqs_ics; observed(sys)]
-    return NonlinearSystem(
-        eqs_ics, vars, pars;
-        defaults = defs, parameter_dependencies = parameter_dependencies(sys),
-        checks = check_units,
-        name, kwargs...
+    # 4) process parameters as initialization unknowns
+    paramsubs = Dict()
+    if pmap isa SciMLBase.NullParameters
+        pmap = Dict()
+    end
+    pmap = todict(pmap)
+    for p in parameters(sys)
+        if is_parameter_solvable(p, pmap, defs, guesses)
+            # If either of them are `missing` the parameter is an unknown
+            # But if the parameter is passed a value, use that as an additional
+            # equation in the system
+            _val1 = get(pmap, p, nothing)
+            _val2 = get(defs, p, nothing)
+            _val3 = get(guesses, p, nothing)
+            varp = tovar(p)
+            paramsubs[p] = varp
+            # Has a default of `missing`, and (either an equation using the value passed to `ODEProblem` or a guess)
+            if _val2 === missing
+                if _val1 !== nothing && _val1 !== missing
+                    push!(eqs_ics, varp ~ _val1)
+                    push!(defs, varp => _val1)
+                elseif _val3 !== nothing
+                    # assuming an equation exists (either via algebraic equations or initialization_eqs)
+                    push!(defs, varp => _val3)
+                elseif check_defguess
+                    error("Invalid setup: parameter $(p) has no default value, initial value, or guess")
+                end
+                # `missing` passed to `ODEProblem`, and (either an equation using default or a guess)
+            elseif _val1 === missing
+                if _val2 !== nothing && _val2 !== missing
+                    push!(eqs_ics, varp ~ _val2)
+                    push!(defs, varp => _val2)
+                elseif _val3 !== nothing
+                    push!(defs, varp => _val3)
+                elseif check_defguess
+                    error("Invalid setup: parameter $(p) has no default value, initial value, or guess")
+                end
+                # given a symbolic value to ODEProblem
+            elseif symbolic_type(_val1) != NotSymbolic()
+                push!(eqs_ics, varp ~ _val1)
+                push!(defs, varp => _val3)
+                # No value passed to `ODEProblem`, but a default and a guess are present
+                # _val2 !== missing is implied by it falling this far in the elseif chain
+            elseif _val1 === nothing && _val2 !== nothing
+                push!(eqs_ics, varp ~ _val2)
+                push!(defs, varp => _val3)
+            else
+                # _val1 !== missing and _val1 !== nothing, so a value was provided to ODEProblem
+                # This would mean `is_parameter_solvable` returned `false`, so we never end up
+                # here
+                error("This should never be reached")
+            end
+        end
+    end
+
+    # 5) parameter dependencies become equations, their LHS become unknowns
+    for eq in parameter_dependencies(sys)
+        varp = tovar(eq.lhs)
+        paramsubs[eq.lhs] = varp
+        push!(eqs_ics, eq)
+        guessval = get(guesses, eq.lhs, eq.rhs)
+        push!(defs, varp => guessval)
+    end
+
+    # 6) handle values provided for dependent parameters similar to values for observed variables
+    for (k, v) in merge(defaults(sys), pmap)
+        if is_variable_floatingpoint(k) && has_parameter_dependency_with_lhs(sys, k)
+            push!(eqs_ics, paramsubs[k] ~ v)
+        end
+    end
+
+    # parameters do not include ones that became initialization unknowns
+    pars = vcat(
+        [get_iv(sys)], # include independent variable as pseudo-parameter
+        [p for p in parameters(sys) if !haskey(paramsubs, p)]
     )
+
+    eqs_ics = Symbolics.substitute.([eqs_ics; observed(sys)], (paramsubs,))
+    vars = [vars; collect(values(paramsubs))]
+    for k in keys(defs)
+        defs[k] = substitute(defs[k], paramsubs)
+    end
+    return NonlinearSystem(eqs_ics,
+        vars,
+        pars;
+        defaults = defs,
+        checks = check_units,
+        name,
+        kwargs...)
+end
+
+function is_parameter_solvable(p, pmap, defs, guesses)
+    _val1 = pmap isa AbstractDict ? get(pmap, p, nothing) : nothing
+    _val2 = get(defs, p, nothing)
+    _val3 = get(guesses, p, nothing)
+    # either (missing is a default or was passed to the ODEProblem) or (nothing was passed to
+    # the ODEProblem and it has a default and a guess)
+    return ((_val1 === missing || _val2 === missing) ||
+           (_val1 === nothing && _val2 !== nothing)) && _val3 !== nothing
 end
diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl
index 20ba3d731a..67b26ea50c 100644
--- a/src/systems/parameter_buffer.jl
+++ b/src/systems/parameter_buffer.jl
@@ -522,6 +522,9 @@ function _remake_buffer(indp, oldbuf::MTKParameters, idxs, vals; validate = true
     ic = get_index_cache(indp_to_system(indp))
     for (idx, val) in zip(idxs, vals)
         sym = nothing
+        if val === missing
+            val = get_temporary_value(idx)
+        end
         if symbolic_type(idx) == ScalarSymbolic()
             sym = idx
             idx = parameter_index(ic, sym)
diff --git a/src/utils.jl b/src/utils.jl
index 25731f495a..a3f99731af 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -494,8 +494,8 @@ function collect_var!(unknowns, parameters, var, iv)
         push!(unknowns, var)
     end
     # Add also any parameters that appear only as defaults in the var
-    if hasdefault(var)
-        collect_vars!(unknowns, parameters, getdefault(var), iv)
+    if hasdefault(var) && (def = getdefault(var)) !== missing
+        collect_vars!(unknowns, parameters, def, iv)
     end
     return nothing
 end

From 8ec6d3e467dede8db4c5c97fefb994b38c6a3af0 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 21 Aug 2024 12:36:44 +0530
Subject: [PATCH 02/19] feat: collect guesses from parameter metadata in
 ODESystem

---
 src/systems/diffeqs/odesystem.jl | 15 ++++++++++-----
 src/variables.jl                 | 10 ++++++++++
 2 files changed, 20 insertions(+), 5 deletions(-)

diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 784b071ef1..48c00a933b 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -249,11 +249,16 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     process_variables!(var_to_name, defaults, dvs′)
     process_variables!(var_to_name, defaults, ps′)
 
-    sysguesses = [ModelingToolkit.getguess(st) for st in dvs′]
-    hasaguess = findall(!isnothing, sysguesses)
-    var_guesses = dvs′[hasaguess] .=> sysguesses[hasaguess]
-    sysguesses = isempty(var_guesses) ? Dict() : todict(var_guesses)
-    guesses = merge(sysguesses, todict(guesses))
+    sysdvsguesses = [ModelingToolkit.getguess(st) for st in dvs′]
+    hasaguess = findall(!isnothing, sysdvsguesses)
+    var_guesses = dvs′[hasaguess] .=> sysdvsguesses[hasaguess]
+    sysdvsguesses = isempty(var_guesses) ? Dict() : todict(var_guesses)
+    syspsguesses = [ModelingToolkit.getguess(st) for st in ps′]
+    hasaguess = findall(!isnothing, syspsguesses)
+    ps_guesses = ps′[hasaguess] .=> syspsguesses[hasaguess]
+    syspsguesses = isempty(ps_guesses) ? Dict() : todict(ps_guesses)
+
+    guesses = merge(sysdvsguesses, syspsguesses, todict(guesses))
     guesses = Dict{Any, Any}(value(k) => value(v) for (k, v) in pairs(guesses))
 
     isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed))
diff --git a/src/variables.jl b/src/variables.jl
index 4e5f860513..1c22540c63 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -491,6 +491,16 @@ function getguess(x)
     Symbolics.getmetadata(x, VariableGuess, nothing)
 end
 
+"""
+    setguess(x, v)
+
+Set the guess for the initial value associated with symbolic variable `x` to `v`.
+See also [`hasguess`](@ref).
+"""
+function setguess(x, v)
+    Symbolics.setmetadata(x, VariableGuess, v)
+end
+
 """
     hasguess(x)
 

From 6d63ffdfef27fef051390106f3bdabf09a002bb2 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 6 Aug 2024 16:19:26 +0530
Subject: [PATCH 03/19] test: test parameters as initialization unknowns

---
 test/initializationsystem.jl   | 130 +++++++++++++++++++++++++++++++++
 test/parameter_dependencies.jl |   4 +-
 2 files changed, 132 insertions(+), 2 deletions(-)

diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index e0caf22ba4..e6192144c8 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -1,4 +1,5 @@
 using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve, Test
+using SymbolicIndexingInterface
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @parameters g
@@ -576,3 +577,132 @@ sol = solve(oprob_2nd_order_2, Rosenbrock23()) # retcode: Success
     sol = solve(prob, Tsit5(), reltol = 1e-4)
     @test all(sol(1.0, idxs = sys.x) .≈ +exp(1)) && all(sol(1.0, idxs = sys.y) .≈ -exp(1))
 end
+
+@testset "Initialization of parameters" begin
+    function test_parameter(prob, sym, val, initialval = zero(val))
+        @test prob.ps[sym] ≈ initialval
+        @test init(prob, Tsit5()).ps[sym] ≈ val
+        @test solve(prob, Tsit5()).ps[sym] ≈ val
+    end
+    function test_initializesystem(sys, u0map, pmap, p, equation)
+        isys = ModelingToolkit.generate_initializesystem(
+            sys; u0map, pmap, guesses = ModelingToolkit.guesses(sys))
+        @test is_variable(isys, p)
+        @test equation in equations(isys) || (0 ~ -equation.rhs) in equations(isys)
+    end
+    @variables x(t) y(t)
+    @parameters p q
+    u0map = Dict(x => 1.0, y => 1.0)
+    pmap = Dict()
+    pmap[q] = 1.0
+    # `missing` default, equation from ODEProblem
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x * q, D(y) ~ y * p], t; defaults = [p => missing], guesses = [p => 1.0])
+    pmap[p] = 2q
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    # `missing` default, provided guess
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x, p ~ x + y], t; defaults = [p => missing], guesses = [p => 0.0])
+    prob = ODEProblem(sys, u0map, (0.0, 1.0))
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ p - x - y)
+
+    # `missing` to ODEProblem, equation from default
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x * q, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 1.0])
+    pmap[p] = missing
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ 2q - p)
+    test_parameter(prob2, p, 2.0)
+    # `missing` to ODEProblem, provided guess
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x, p ~ x + y], t; guesses = [p => 0.0])
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ x + y - p)
+
+    # No `missing`, default and guess
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x * q, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 0.0])
+    delete!(pmap, p)
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ 2q - p)
+
+    # Should not be solved for:
+
+    # ODEProblem value with guess, no `missing`
+    @mtkbuild sys = ODESystem([D(x) ~ x * q, D(y) ~ y * p], t; guesses = [p => 0.0])
+    _pmap = merge(pmap, Dict(p => 3q))
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), _pmap)
+    @test prob.ps[p] ≈ 3.0
+    @test prob.f.initializeprob === nothing
+    # Default overridden by ODEProblem, guess provided
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ q * x, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 1.0])
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), _pmap)
+    @test prob.ps[p] ≈ 3.0
+    @test prob.f.initializeprob === nothing
+
+    @mtkbuild sys = ODESystem([D(x) ~ x, p ~ x + y], t; guesses = [p => 0.0])
+    @test_throws ModelingToolkit.MissingParametersError ODEProblem(
+        sys, [x => 1.0, y => 1.0], (0.0, 1.0))
+
+    @testset "Null system" begin
+        @variables x(t) y(t) s(t)
+        @parameters x0 y0
+        @mtkbuild sys = ODESystem([x ~ x0, y ~ y0, s ~ x + y], t; guesses = [y0 => 0.0])
+        prob = ODEProblem(sys, [s => 1.0], (0.0, 1.0), [x0 => 0.3, y0 => missing])
+        test_parameter(prob, y0, 0.7)
+    end
+
+    using ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica: Fixed, Mass,
+                                                                           Spring, Force,
+                                                                           Damper
+    using ModelingToolkitStandardLibrary.Mechanical: TranslationalModelica as TM
+    using ModelingToolkitStandardLibrary.Blocks: Constant
+
+    @named mass = TM.Mass(; m = 1.0, s = 1.0, v = 0.0, a = 0.0)
+    @named fixed = Fixed(; s0 = 0.0)
+    @named spring = Spring(; c = 2.0, s_rel0 = nothing)
+    @named gravity = Force()
+    @named constant = Constant(; k = 9.81)
+    @named damper = TM.Damper(; d = 0.1)
+    @mtkbuild sys = ODESystem(
+        [connect(fixed.flange, spring.flange_a), connect(spring.flange_b, mass.flange_a),
+            connect(mass.flange_a, gravity.flange), connect(constant.output, gravity.f),
+            connect(fixed.flange, damper.flange_a), connect(damper.flange_b, mass.flange_a)],
+        t;
+        systems = [fixed, spring, mass, gravity, constant, damper],
+        guesses = [spring.s_rel0 => 1.0])
+    prob = ODEProblem(sys, [], (0.0, 1.0), [spring.s_rel0 => missing])
+    test_parameter(prob, spring.s_rel0, -3.905)
+end
+
+@testset "Update initializeprob parameters" begin
+    @variables x(t) y(t)
+    @parameters p q
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x, p ~ x + y], t; guesses = [x => 0.0, p => 0.0])
+    prob = ODEProblem(sys, [y => 1.0], (0.0, 1.0), [p => 3.0])
+    @test prob.f.initializeprob.ps[p] ≈ 3.0
+    @test init(prob, Tsit5())[x] ≈ 2.0
+    prob.ps[p] = 2.0
+    @test prob.f.initializeprob.ps[p] ≈ 3.0
+    @test init(prob, Tsit5())[x] ≈ 1.0
+    ModelingToolkit.defaults(prob.f.sys)[p] = missing
+end
+
+@testset "Equations for dependent parameters" begin
+    @variables x(t)
+    @parameters p q=5 r
+    @mtkbuild sys = ODESystem(
+        D(x) ~ 2x + r, t; parameter_dependencies = [r ~ p + 2q, q ~ p + 3],
+        guesses = [p => 1.0])
+    prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p => missing])
+    @test length(equations(ModelingToolkit.get_parent(prob.f.initializeprob.f.sys))) == 4
+    integ = init(prob, Tsit5())
+    @test integ.ps[p] ≈ 2
+end
diff --git a/test/parameter_dependencies.jl b/test/parameter_dependencies.jl
index 9cfd4ca5c5..1f1aab9953 100644
--- a/test/parameter_dependencies.jl
+++ b/test/parameter_dependencies.jl
@@ -10,7 +10,7 @@ using SymbolicIndexingInterface
 using NonlinearSolve
 
 @testset "ODESystem with callbacks" begin
-    @parameters p1=1.0 p2=1.0
+    @parameters p1=1.0 p2
     @variables x(t)
     cb1 = [x ~ 2.0] => [p1 ~ 2.0] # triggers at t=-2+√6
     function affect1!(integ, u, p, ctx)
@@ -31,7 +31,7 @@ using NonlinearSolve
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.5), jac = true)
     @test prob.ps[p1] == 1.0
     @test prob.ps[p2] == 2.0
-    @test_nowarn solve(prob, Tsit5())
+    @test SciMLBase.successful_retcode(solve(prob, Tsit5()))
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.5), [p1 => 1.0], jac = true)
     @test prob.ps[p1] == 1.0
     @test prob.ps[p2] == 2.0

From 4682f8d2d968d2ed7c44a1fc87d3a864783ce5a7 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 6 Aug 2024 17:53:15 +0530
Subject: [PATCH 04/19] fix: check symbolic_type in `namespace_expr`

---
 src/systems/abstractsystem.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 661a5fd4b4..bad83133a1 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -1356,6 +1356,7 @@ end
 function namespace_expr(
         O, sys, n = nameof(sys); ivs = independent_variables(sys))
     O = unwrap(O)
+    symbolic_type(O) == NotSymbolic() && return O
     if any(isequal(O), ivs)
         return O
     elseif iscall(O)

From 8425bb83060d9d57dfe61fc5860438c7d233b6bc Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 6 Aug 2024 17:53:31 +0530
Subject: [PATCH 05/19] refactor: allow parameter default to be `missing` in
 `@mtkmodel`

---
 src/systems/model_parsing.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index fce2add6b7..70fab5a7cc 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -154,9 +154,9 @@ function update_kwargs_and_metadata!(dict, kwargs, a, def, indices, type, var,
         if !isnothing(meta) && haskey(meta, VariableUnit)
             uvar = gensym()
             push!(where_types, uvar)
-            push!(kwargs, Expr(:kw, :($a::Union{Nothing, $uvar}), nothing))
+            push!(kwargs, Expr(:kw, :($a::Union{Nothing, Missing, $uvar}), nothing))
         else
-            push!(kwargs, Expr(:kw, :($a::Union{Nothing, $type}), nothing))
+            push!(kwargs, Expr(:kw, :($a::Union{Nothing, Missing, $type}), nothing))
         end
         dict[:kwargs][getname(var)] = Dict(:value => def, :type => type)
     else

From 5a538c3a7807f1744abeebcd60d98c1b6cd354e3 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 5 Jul 2024 19:56:35 +0530
Subject: [PATCH 06/19] feat: support `SciMLBase.remake_initializeprob`

---
 src/systems/nonlinear/initializesystem.jl | 53 +++++++++++++++++++++++
 test/initializationsystem.jl              |  5 +++
 2 files changed, 58 insertions(+)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index d4253afda5..a4a29f223b 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -175,3 +175,56 @@ function is_parameter_solvable(p, pmap, defs, guesses)
     return ((_val1 === missing || _val2 === missing) ||
            (_val1 === nothing && _val2 !== nothing)) && _val3 !== nothing
 end
+
+function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p)
+    if (u0 === missing || !(eltype(u0) <: Pair) || isempty(u0)) &&
+       (p === missing || !(eltype(p) <: Pair) || isempty(p))
+        return odefn.initializeprob, odefn.update_initializeprob!, odefn.initializeprobmap,
+        odefn.initializeprobpmap
+    end
+    if u0 === missing || isempty(u0)
+        u0 = Dict()
+    elseif !(eltype(u0) <: Pair)
+        u0 = Dict(unknowns(sys) .=> u0)
+    end
+    if p === missing
+        p = Dict()
+    end
+    if t0 === nothing
+        t0 = 0.0
+    end
+    u0 = todict(u0)
+    defs = defaults(sys)
+    varmap = merge(defs, u0)
+    varmap = canonicalize_varmap(varmap)
+    missingvars = setdiff(unknowns(sys), collect(keys(varmap)))
+    setobserved = filter(keys(varmap)) do var
+        has_observed_with_lhs(sys, var) || has_observed_with_lhs(sys, default_toterm(var))
+    end
+    p = todict(p)
+    guesses = ModelingToolkit.guesses(sys)
+    solvablepars = [par
+                    for par in parameters(sys)
+                    if is_parameter_solvable(par, p, defs, guesses)]
+    pvarmap = merge(defs, p)
+    setparobserved = filter(keys(pvarmap)) do var
+        has_parameter_dependency_with_lhs(sys, var)
+    end
+    if (((!isempty(missingvars) || !isempty(solvablepars) ||
+          !isempty(setobserved) || !isempty(setparobserved)) &&
+         ModelingToolkit.get_tearing_state(sys) !== nothing) ||
+        !isempty(initialization_equations(sys)))
+        initprob = InitializationProblem(sys, t0, u0, p)
+        initprobmap = getu(initprob, unknowns(sys))
+        punknowns = [p for p in all_variable_symbols(initprob) if is_parameter(sys, p)]
+        getpunknowns = getu(initprob, punknowns)
+        setpunknowns = setp(sys, punknowns)
+        initprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns)
+        reqd_syms = parameter_symbols(initprob)
+        update_initializeprob! = UpdateInitializeprob(
+            getu(sys, reqd_syms), setu(initprob, reqd_syms))
+        return initprob, update_initializeprob!, initprobmap, initprobpmap
+    else
+        return nothing, nothing, nothing, nothing
+    end
+end
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index e6192144c8..353eaa24e2 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -693,6 +693,11 @@ end
     @test prob.f.initializeprob.ps[p] ≈ 3.0
     @test init(prob, Tsit5())[x] ≈ 1.0
     ModelingToolkit.defaults(prob.f.sys)[p] = missing
+    prob2 = remake(prob; u0 = [y => 1.0], p = [p => 3x])
+    @test !is_variable(prob2.f.initializeprob, p) &&
+          !is_parameter(prob2.f.initializeprob, p)
+    @test init(prob2, Tsit5())[x] ≈ 0.5
+    @test_nowarn solve(prob2, Tsit5())
 end
 
 @testset "Equations for dependent parameters" begin

From f783cb105f9e88117bf637151551f004e42a6c5a Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 21 Aug 2024 14:19:08 +0530
Subject: [PATCH 07/19] test: test remaking of initialization problem

---
 test/initializationsystem.jl | 35 +++++++++++++++++++++++++++++++++++
 test/mtkparameters.jl        |  2 +-
 test/split_parameters.jl     |  5 +++--
 3 files changed, 39 insertions(+), 3 deletions(-)

diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 353eaa24e2..d9757ffc5b 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -601,12 +601,18 @@ end
     pmap[p] = 2q
     prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
     test_parameter(prob, p, 2.0)
+    prob2 = remake(prob; u0 = u0map, p = pmap)
+    prob2.ps[p] = 0.0
+    test_parameter(prob2, p, 2.0)
     # `missing` default, provided guess
     @mtkbuild sys = ODESystem(
         [D(x) ~ x, p ~ x + y], t; defaults = [p => missing], guesses = [p => 0.0])
     prob = ODEProblem(sys, u0map, (0.0, 1.0))
     test_parameter(prob, p, 2.0)
     test_initializesystem(sys, u0map, pmap, p, 0 ~ p - x - y)
+    prob2 = remake(prob; u0 = u0map)
+    prob2.ps[p] = 0.0
+    test_parameter(prob2, p, 2.0)
 
     # `missing` to ODEProblem, equation from default
     @mtkbuild sys = ODESystem(
@@ -615,6 +621,8 @@ end
     prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
     test_parameter(prob, p, 2.0)
     test_initializesystem(sys, u0map, pmap, p, 0 ~ 2q - p)
+    prob2 = remake(prob; u0 = u0map, p = pmap)
+    prob2.ps[p] = 0.0
     test_parameter(prob2, p, 2.0)
     # `missing` to ODEProblem, provided guess
     @mtkbuild sys = ODESystem(
@@ -622,6 +630,9 @@ end
     prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
     test_parameter(prob, p, 2.0)
     test_initializesystem(sys, u0map, pmap, p, 0 ~ x + y - p)
+    prob2 = remake(prob; u0 = u0map, p = pmap)
+    prob2.ps[p] = 0.0
+    test_parameter(prob2, p, 2.0)
 
     # No `missing`, default and guess
     @mtkbuild sys = ODESystem(
@@ -630,6 +641,9 @@ end
     prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
     test_parameter(prob, p, 2.0)
     test_initializesystem(sys, u0map, pmap, p, 0 ~ 2q - p)
+    prob2 = remake(prob; u0 = u0map, p = pmap)
+    prob2.ps[p] = 0.0
+    test_parameter(prob2, p, 2.0)
 
     # Should not be solved for:
 
@@ -711,3 +725,24 @@ end
     integ = init(prob, Tsit5())
     @test integ.ps[p] ≈ 2
 end
+
+@testset "Re-creating initialization problem on remake" begin
+    @variables x(t) y(t)
+    @parameters p q
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x, p ~ x + y], t; defaults = [p => missing], guesses = [x => 0.0, p => 0.0])
+    prob = ODEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0))
+    @test init(prob, Tsit5()).ps[p] ≈ 2.0
+    # nonsensical value for y just to test that equations work
+    prob2 = remake(prob; u0 = [x => 1.0, y => 2x + exp(t)])
+    @test init(prob2, Tsit5()).ps[p] ≈ 4.0
+    # solve for `x` given `p` and `y`
+    prob3 = remake(prob; u0 = [x => nothing, y => 1.0], p = [p => 2x + exp(t)])
+    @test init(prob3, Tsit5())[x] ≈ 0.0
+    @test_logs (:warn, r"overdetermined") remake(
+        prob; u0 = [x => 1.0, y => 2.0], p = [p => 4.0])
+    prob4 = remake(prob; u0 = [x => 1.0, y => 2.0], p = [p => 4.0])
+    @test solve(prob4, Tsit5()).retcode == ReturnCode.InitialFailure
+    prob5 = remake(prob)
+    @test init(prob, Tsit5()).ps[p] ≈ 2.0
+end
diff --git a/test/mtkparameters.jl b/test/mtkparameters.jl
index e0aee8c289..b4fc1a9677 100644
--- a/test/mtkparameters.jl
+++ b/test/mtkparameters.jl
@@ -293,7 +293,7 @@ end
     ps_vec = [k => [2.0, 3.0, 4.0, 5.0]]
     ps_scal = [k[1] => 1.0, k[2] => 2.0, k[3] => 3.0, k[4] => 4.0]
     oprob_scal_scal = ODEProblem(osys_scal, u0, 1.0, ps_scal)
-    newoprob = remake(oprob_scal_scal; p = ps_vec)
+    newoprob = remake(oprob_scal_scal; p = ps_vec, build_initializeprob = false)
     @test newoprob.ps[k] == [2.0, 3.0, 4.0, 5.0]
 end
 
diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index 9f124ef45d..b8651238ea 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -114,7 +114,7 @@ eqs = [D(y) ~ dy * a
 sys = structural_simplify(model; split = false)
 
 tspan = (0.0, t_end)
-prob = ODEProblem(sys, [], tspan, [])
+prob = ODEProblem(sys, [], tspan, []; build_initializeprob = false)
 
 @test prob.p isa Vector{Float64}
 sol = solve(prob, ImplicitEuler());
@@ -122,7 +122,8 @@ sol = solve(prob, ImplicitEuler());
 
 # ------------------------ Mixed Type Conserved
 
-prob = ODEProblem(sys, [], tspan, []; tofloat = false, use_union = true)
+prob = ODEProblem(
+    sys, [], tspan, []; tofloat = false, use_union = true, build_initializeprob = false)
 
 @test prob.p isa Tuple{Vector{Float64}, Vector{Int64}}
 sol = solve(prob, ImplicitEuler());

From 26c3f3acc61516be334109f4e882dc0431f97837 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 4 Sep 2024 19:27:52 +0530
Subject: [PATCH 08/19] fix: handle edge cases in namespacing

---
 src/systems/abstractsystem.jl | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index bad83133a1..9908d6f1e3 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -1353,10 +1353,21 @@ function namespace_assignment(eq::Assignment, sys)
     Assignment(_lhs, _rhs)
 end
 
+function is_array_of_symbolics(x)
+    symbolic_type(x) == ArraySymbolic() && return true
+    symbolic_type(x) == ScalarSymbolic() && return false
+    x isa AbstractArray &&
+        any(y -> symbolic_type(y) != NotSymbolic() || is_array_of_symbolics(y), x)
+end
+
 function namespace_expr(
         O, sys, n = nameof(sys); ivs = independent_variables(sys))
     O = unwrap(O)
-    symbolic_type(O) == NotSymbolic() && return O
+    # Exceptions for arrays of symbolic and Ref of a symbolic, the latter
+    # of which shows up in broadcasts
+    if symbolic_type(O) == NotSymbolic() && !(O isa AbstractArray) && !(O isa Ref)
+        return O
+    end
     if any(isequal(O), ivs)
         return O
     elseif iscall(O)
@@ -1378,7 +1389,7 @@ function namespace_expr(
         end
     elseif isvariable(O)
         renamespace(n, O)
-    elseif O isa Array
+    elseif O isa AbstractArray && is_array_of_symbolics(O)
         let sys = sys, n = n
             map(o -> namespace_expr(o, sys, n; ivs), O)
         end

From de94dc1a724225c53fba2073f7874d6b9a484f1d Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 6 Sep 2024 17:01:45 +0530
Subject: [PATCH 09/19] fix: allow passing `nothing` to override default for
 `@mtkmodel` models

---
 src/systems/model_parsing.jl | 18 ++++++++++++------
 test/model_parsing.jl        |  2 +-
 2 files changed, 13 insertions(+), 7 deletions(-)

diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index 70fab5a7cc..7b6705bbd5 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -148,15 +148,20 @@ end
 
 pop_structure_dict!(dict, key) = length(dict[key]) == 0 && pop!(dict, key)
 
+struct NoValue end
+const NO_VALUE = NoValue()
+
 function update_kwargs_and_metadata!(dict, kwargs, a, def, indices, type, var,
         varclass, where_types, meta)
     if indices isa Nothing
         if !isnothing(meta) && haskey(meta, VariableUnit)
             uvar = gensym()
             push!(where_types, uvar)
-            push!(kwargs, Expr(:kw, :($a::Union{Nothing, Missing, $uvar}), nothing))
+            push!(kwargs,
+                Expr(:kw, :($a::Union{Nothing, Missing, $NoValue, $uvar}), NO_VALUE))
         else
-            push!(kwargs, Expr(:kw, :($a::Union{Nothing, Missing, $type}), nothing))
+            push!(kwargs,
+                Expr(:kw, :($a::Union{Nothing, Missing, $NoValue, $type}), NO_VALUE))
         end
         dict[:kwargs][getname(var)] = Dict(:value => def, :type => type)
     else
@@ -164,8 +169,9 @@ function update_kwargs_and_metadata!(dict, kwargs, a, def, indices, type, var,
         push!(kwargs,
             Expr(:kw,
                 Expr(:(::), a,
-                    Expr(:curly, :Union, :Nothing, Expr(:curly, :AbstractArray, vartype))),
-                nothing))
+                    Expr(:curly, :Union, :Nothing, :Missing, NoValue,
+                        Expr(:curly, :AbstractArray, vartype))),
+                NO_VALUE))
         if !isnothing(meta) && haskey(meta, VariableUnit)
             push!(where_types, vartype)
         else
@@ -679,7 +685,7 @@ function parse_variable_arg(dict, mod, arg, varclass, kwargs, where_types)
     varexpr = if haskey(metadata_with_exprs, VariableUnit)
         unit = metadata_with_exprs[VariableUnit]
         quote
-            $name = if $name === nothing
+            $name = if $name === $NO_VALUE
                 $setdefault($vv, $def)
             else
                 try
@@ -699,7 +705,7 @@ function parse_variable_arg(dict, mod, arg, varclass, kwargs, where_types)
         end
     else
         quote
-            $name = if $name === nothing
+            $name = if $name === $NO_VALUE
                 $setdefault($vv, $def)
             else
                 $setdefault($vv, $name)
diff --git a/test/model_parsing.jl b/test/model_parsing.jl
index 6332dc12a5..f80dc77dc4 100644
--- a/test/model_parsing.jl
+++ b/test/model_parsing.jl
@@ -840,7 +840,7 @@ end
             n[1:3] = if flag
                 [2, 2, 2]
             else
-                1
+                [1, 1, 1]
             end
         end
     end

From 5874f97ed519a8bb18c83d9ec9b8bb8f794e4f74 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 6 Sep 2024 17:02:02 +0530
Subject: [PATCH 10/19] fix: fix `nothing` default overrides being ignored

---
 src/systems/diffeqs/odesystem.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 48c00a933b..539fd1a975 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -242,12 +242,12 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
             "`default_u0` and `default_p` are deprecated. Use `defaults` instead.",
             :ODESystem, force = true)
     end
-    defaults = todict(defaults)
-    defaults = Dict{Any, Any}(value(k) => value(v)
-    for (k, v) in pairs(defaults) if value(v) !== nothing)
+    defaults = Dict{Any, Any}(todict(defaults))
     var_to_name = Dict()
     process_variables!(var_to_name, defaults, dvs′)
     process_variables!(var_to_name, defaults, ps′)
+    defaults = Dict{Any, Any}(value(k) => value(v)
+    for (k, v) in pairs(defaults) if v !== nothing)
 
     sysdvsguesses = [ModelingToolkit.getguess(st) for st in dvs′]
     hasaguess = findall(!isnothing, sysdvsguesses)

From d7755bfb9464469cd983f437fc9ab045c4398b08 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 26 Sep 2024 19:24:33 +0530
Subject: [PATCH 11/19] feat: track parameter dependency defaults and guesses

---
 src/systems/diffeqs/odesystem.jl | 17 +++++++++++++----
 1 file changed, 13 insertions(+), 4 deletions(-)

diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 539fd1a975..bda5690a67 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -237,6 +237,8 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     ctrl′ = value.(controls)
     dvs′ = value.(dvs)
     dvs′ = filter(x -> !isdelay(x, iv), dvs′)
+    parameter_dependencies, ps′ = process_parameter_dependencies(
+        parameter_dependencies, ps′)
     if !(isempty(default_u0) && isempty(default_p))
         Base.depwarn(
             "`default_u0` and `default_p` are deprecated. Use `defaults` instead.",
@@ -246,6 +248,8 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     var_to_name = Dict()
     process_variables!(var_to_name, defaults, dvs′)
     process_variables!(var_to_name, defaults, ps′)
+    process_variables!(var_to_name, defaults, [eq.lhs for eq in parameter_dependencies])
+    process_variables!(var_to_name, defaults, [eq.rhs for eq in parameter_dependencies])
     defaults = Dict{Any, Any}(value(k) => value(v)
     for (k, v) in pairs(defaults) if v !== nothing)
 
@@ -257,9 +261,15 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     hasaguess = findall(!isnothing, syspsguesses)
     ps_guesses = ps′[hasaguess] .=> syspsguesses[hasaguess]
     syspsguesses = isempty(ps_guesses) ? Dict() : todict(ps_guesses)
+    syspdepguesses = [ModelingToolkit.getguess(eq.lhs) for eq in parameter_dependencies]
+    hasaguess = findall(!isnothing, syspdepguesses)
+    pdep_guesses = [eq.lhs for eq in parameter_dependencies][hasaguess] .=>
+        syspdepguesses[hasaguess]
+    syspdepguesses = isempty(pdep_guesses) ? Dict() : todict(pdep_guesses)
 
-    guesses = merge(sysdvsguesses, syspsguesses, todict(guesses))
-    guesses = Dict{Any, Any}(value(k) => value(v) for (k, v) in pairs(guesses))
+    guesses = merge(sysdvsguesses, syspsguesses, syspdepguesses, todict(guesses))
+    guesses = Dict{Any, Any}(value(k) => value(v)
+    for (k, v) in pairs(guesses) if v !== nothing)
 
     isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed))
 
@@ -274,8 +284,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     end
     cont_callbacks = SymbolicContinuousCallbacks(continuous_events)
     disc_callbacks = SymbolicDiscreteCallbacks(discrete_events)
-    parameter_dependencies, ps′ = process_parameter_dependencies(
-        parameter_dependencies, ps′)
+
     if is_dde === nothing
         is_dde = _check_if_dde(deqs, iv′, systems)
     end

From 96faa983161c677c33e351361238163b65a75dd8 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 1 Oct 2024 15:08:21 +0530
Subject: [PATCH 12/19] fix: handle type promotion in `InitializationProblem`
 with observed

---
 src/systems/diffeqs/abstractodesystem.jl | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 623b2b47af..8ebe6d93d0 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -888,7 +888,8 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
                         for p in parameters(sys)
                         if is_parameter_solvable(p, parammap, defs, guesses)]
 
-        pvarmap = if parammap === nothing || parammap == SciMLBase.NullParameters() || !(eltype(parammap) <: Pair) && isempty(parammap)
+        pvarmap = if parammap === nothing || parammap == SciMLBase.NullParameters() ||
+                     !(eltype(parammap) <: Pair) && isempty(parammap)
             defs
         else
             merge(defs, todict(parammap))
@@ -1592,6 +1593,22 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
     end
 
     u0map = merge(ModelingToolkit.guesses(sys), todict(guesses), todict(u0map))
+    fullmap = merge(u0map, parammap)
+    u0T = Union{}
+    for sym in unknowns(isys)
+        haskey(fullmap, sym) || continue
+        symbolic_type(fullmap[sym]) == NotSymbolic() || continue
+        u0T = promote_type(u0T, typeof(fullmap[sym]))
+    end
+    for eq in observed(isys)
+        haskey(fullmap, eq.lhs) || continue
+        symbolic_type(fullmap[eq.lhs]) == NotSymbolic() || continue
+        u0T = promote_type(u0T, typeof(fullmap[eq.lhs]))
+    end
+    if u0T != Union{}
+        u0map = Dict(k => symbolic_type(v) == NotSymbolic() ? u0T(v) : v
+        for (k, v) in u0map)
+    end
     if neqs == nunknown
         NonlinearProblem(isys, u0map, parammap; kwargs...)
     else

From 8836abec485bd0d2d7504c9960b37bfa73fee3d5 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 1 Oct 2024 15:08:44 +0530
Subject: [PATCH 13/19] fix: improve type promotion in `remake_initializeprob`

---
 src/systems/nonlinear/initializesystem.jl | 45 +++++++++++++++++++++--
 test/initializationsystem.jl              | 36 +++++++++++++++++-
 2 files changed, 77 insertions(+), 4 deletions(-)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index a4a29f223b..5b465a677c 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -173,15 +173,54 @@ function is_parameter_solvable(p, pmap, defs, guesses)
     # either (missing is a default or was passed to the ODEProblem) or (nothing was passed to
     # the ODEProblem and it has a default and a guess)
     return ((_val1 === missing || _val2 === missing) ||
-           (_val1 === nothing && _val2 !== nothing)) && _val3 !== nothing
+            (_val1 === nothing && _val2 !== nothing)) && _val3 !== nothing
 end
 
 function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p)
-    if (u0 === missing || !(eltype(u0) <: Pair) || isempty(u0)) &&
-       (p === missing || !(eltype(p) <: Pair) || isempty(p))
+    if u0 === missing && p === missing
         return odefn.initializeprob, odefn.update_initializeprob!, odefn.initializeprobmap,
         odefn.initializeprobpmap
     end
+    if !(eltype(u0) <: Pair) && !(eltype(p) <: Pair)
+        oldinitprob = odefn.initializeprob
+        if oldinitprob === nothing || !SciMLBase.has_sys(oldinitprob.f) ||
+           !(oldinitprob.f.sys isa NonlinearSystem)
+            return oldinitprob, odefn.update_initializeprob!, odefn.initializeprobmap,
+            odefn.initializeprobpmap
+        end
+        pidxs = ParameterIndex[]
+        pvals = []
+        u0idxs = Int[]
+        u0vals = []
+        for sym in variable_symbols(oldinitprob)
+            if is_variable(sys, sym)
+                u0 !== missing || continue
+                idx = variable_index(oldinitprob, sym)
+                push!(u0idxs, idx)
+                push!(u0vals, eltype(u0)(state_values(oldinitprob, idx)))
+            else
+                p !== missing || continue
+                idx = variable_index(oldinitprob, sym)
+                push!(u0idxs, idx)
+                push!(u0vals, typeof(getp(sys, sym)(p))(state_values(oldinitprob, idx)))
+            end
+        end
+        if p !== missing
+            for sym in parameter_symbols(oldinitprob)
+                push!(pidxs, parameter_index(oldinitprob, sym))
+                if isequal(sym, get_iv(sys))
+                    push!(pvals, t0)
+                else
+                    push!(pvals, getp(sys, sym)(p))
+                end
+            end
+        end
+        newu0 = remake_buffer(oldinitprob.f.sys, state_values(oldinitprob), u0idxs, u0vals)
+        newp = remake_buffer(oldinitprob.f.sys, parameter_values(oldinitprob), pidxs, pvals)
+        initprob = remake(oldinitprob; u0 = newu0, p = newp)
+        return initprob, odefn.update_initializeprob!, odefn.initializeprobmap,
+        odefn.initializeprobpmap
+    end
     if u0 === missing || isempty(u0)
         u0 = Dict()
     elseif !(eltype(u0) <: Pair)
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index d9757ffc5b..f30dc5fa3a 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -1,5 +1,7 @@
 using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve, Test
-using SymbolicIndexingInterface
+using ForwardDiff
+using SymbolicIndexingInterface, SciMLStructures
+using SciMLStructures: Tunable
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @parameters g
@@ -746,3 +748,35 @@ end
     prob5 = remake(prob)
     @test init(prob, Tsit5()).ps[p] ≈ 2.0
 end
+
+@testset "`remake` changes initialization problem types" begin
+    @variables x(t) y(t) z(t)
+    @parameters p q
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x * p + y * q, y^2 * q + q^2 * x ~ 0, z * p - p^2 * x * z ~ 0],
+        t; guesses = [x => 0.0, y => 0.0, z => 0.0, p => 0.0, q => 0.0])
+    prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0, q => missing])
+    @test is_variable(prob.f.initializeprob, q)
+    ps = prob.p
+    newps = SciMLStructures.replace(Tunable(), ps, ForwardDiff.Dual.(ps.tunable))
+    prob2 = remake(prob; p = newps)
+    @test eltype(prob2.f.initializeprob.u0) <: ForwardDiff.Dual
+    @test eltype(prob2.f.initializeprob.p.tunable) <: ForwardDiff.Dual
+    @test prob2.f.initializeprob.u0 ≈ prob.f.initializeprob.u0
+
+    prob2 = remake(prob; u0 = ForwardDiff.Dual.(prob.u0))
+    @test eltype(prob2.f.initializeprob.u0) <: ForwardDiff.Dual
+    @test eltype(prob2.f.initializeprob.p.tunable) <: Float64
+    @test prob2.f.initializeprob.u0 ≈ prob.f.initializeprob.u0
+
+    prob2 = remake(prob; u0 = ForwardDiff.Dual.(prob.u0), p = newps)
+    @test eltype(prob2.f.initializeprob.u0) <: ForwardDiff.Dual
+    @test eltype(prob2.f.initializeprob.p.tunable) <: ForwardDiff.Dual
+    @test prob2.f.initializeprob.u0 ≈ prob.f.initializeprob.u0
+
+    prob2 = remake(prob; u0 = [x => ForwardDiff.Dual(1.0)],
+        p = [p => ForwardDiff.Dual(1.0), q => missing])
+    @test eltype(prob2.f.initializeprob.u0) <: ForwardDiff.Dual
+    @test eltype(prob2.f.initializeprob.p.tunable) <: ForwardDiff.Dual
+    @test prob2.f.initializeprob.u0 ≈ prob.f.initializeprob.u0
+end

From 015ce60b2d528f69dad83317100647985a817b02 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 1 Oct 2024 17:30:51 +0530
Subject: [PATCH 14/19] feat: preserve the `u0map` passed to `ODEProblem` and
 use it in `remake_initializeprob`

---
 src/systems/nonlinear/initializesystem.jl | 16 ++++++++++++++++
 test/initializationsystem.jl              | 16 ++++++++++++++++
 2 files changed, 32 insertions(+)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 5b465a677c..f0a07bfc25 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -157,15 +157,22 @@ function generate_initializesystem(sys::ODESystem;
     for k in keys(defs)
         defs[k] = substitute(defs[k], paramsubs)
     end
+    meta = InitializationSystemMetadata(Dict{Any, Any}(u0map), Dict{Any, Any}(pmap))
     return NonlinearSystem(eqs_ics,
         vars,
         pars;
         defaults = defs,
         checks = check_units,
         name,
+        metadata = meta,
         kwargs...)
 end
 
+struct InitializationSystemMetadata
+    u0map::Dict{Any, Any}
+    pmap::Dict{Any, Any}
+end
+
 function is_parameter_solvable(p, pmap, defs, guesses)
     _val1 = pmap isa AbstractDict ? get(pmap, p, nothing) : nothing
     _val2 = get(defs, p, nothing)
@@ -253,6 +260,15 @@ function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p)
           !isempty(setobserved) || !isempty(setparobserved)) &&
          ModelingToolkit.get_tearing_state(sys) !== nothing) ||
         !isempty(initialization_equations(sys)))
+        if SciMLBase.has_initializeprob(odefn)
+            oldsys = odefn.initializeprob.f.sys
+            meta = get_metadata(oldsys)
+            if meta isa InitializationSystemMetadata
+                u0 = merge(meta.u0map, u0)
+                p = merge(meta.pmap, p)
+            end
+        end
+
         initprob = InitializationProblem(sys, t0, u0, p)
         initprobmap = getu(initprob, unknowns(sys))
         punknowns = [p for p in all_variable_symbols(initprob) if is_parameter(sys, p)]
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index f30dc5fa3a..77c015144e 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -780,3 +780,19 @@ end
     @test eltype(prob2.f.initializeprob.p.tunable) <: ForwardDiff.Dual
     @test prob2.f.initializeprob.u0 ≈ prob.f.initializeprob.u0
 end
+
+@testset "`remake` preserves old u0map and pmap" begin
+    @variables x(t) y(t)
+    @parameters p
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x + p * y, y^2 + 4y * p^2 ~ x], t; guesses = [y => 1.0, p => 1.0])
+    prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0])
+    @test is_variable(prob.f.initializeprob, y)
+    prob2 = @test_nowarn remake(prob; p = [p => 3.0]) # ensure no over/under-determined warning
+    @test is_variable(prob.f.initializeprob, y)
+
+    prob = ODEProblem(sys, [y => 1.0, x => 2.0], (0.0, 1.0), [p => missing])
+    @test is_variable(prob.f.initializeprob, p)
+    prob2 = @test_nowarn remake(prob; u0 = [y => 0.5])
+    @test is_variable(prob.f.initializeprob, p)
+end

From cab11eb594002aea5ece2f5948c21a103878ccf1 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 7 Oct 2024 18:42:32 +0530
Subject: [PATCH 15/19] feat: validate parameter type and allow dependent
 initial values in param init

---
 src/systems/nonlinear/initializesystem.jl |  5 +++-
 src/utils.jl                              |  7 ++++++
 test/initializationsystem.jl              | 30 ++++++++++++++++++-----
 3 files changed, 35 insertions(+), 7 deletions(-)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index f0a07bfc25..6a21e9a6f0 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -174,13 +174,16 @@ struct InitializationSystemMetadata
 end
 
 function is_parameter_solvable(p, pmap, defs, guesses)
+    p = unwrap(p)
+    is_variable_floatingpoint(p) || return false
     _val1 = pmap isa AbstractDict ? get(pmap, p, nothing) : nothing
     _val2 = get(defs, p, nothing)
     _val3 = get(guesses, p, nothing)
     # either (missing is a default or was passed to the ODEProblem) or (nothing was passed to
     # the ODEProblem and it has a default and a guess)
     return ((_val1 === missing || _val2 === missing) ||
-            (_val1 === nothing && _val2 !== nothing)) && _val3 !== nothing
+            (symbolic_type(_val1) != NotSymbolic() ||
+             _val1 === nothing && _val2 !== nothing)) && _val3 !== nothing
 end
 
 function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p)
diff --git a/src/utils.jl b/src/utils.jl
index a3f99731af..29e61ede69 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -885,3 +885,10 @@ end
 
 diff2term_with_unit(x, t) = _with_unit(diff2term, x, t)
 lower_varname_with_unit(var, iv, order) = _with_unit(lower_varname, var, iv, iv, order)
+
+function is_variable_floatingpoint(sym)
+    sym = unwrap(sym)
+    T = symtype(sym)
+    return T == Real || T <: AbstractFloat || T <: AbstractArray{Real} ||
+           T <: AbstractArray{<:AbstractFloat}
+end
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 77c015144e..08b653c22e 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -647,19 +647,37 @@ end
     prob2.ps[p] = 0.0
     test_parameter(prob2, p, 2.0)
 
-    # Should not be solved for:
+    # Default overridden by ODEProblem, guess provided
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ q * x, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 1.0])
+    _pmap = merge(pmap, Dict(p => q))
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), _pmap)
+    test_parameter(prob, p, _pmap[q])
+    test_initializesystem(sys, u0map, _pmap, p, 0 ~ q - p)
 
-    # ODEProblem value with guess, no `missing`
+    # ODEProblem dependent value with guess, no `missing`
     @mtkbuild sys = ODESystem([D(x) ~ x * q, D(y) ~ y * p], t; guesses = [p => 0.0])
     _pmap = merge(pmap, Dict(p => 3q))
     prob = ODEProblem(sys, u0map, (0.0, 1.0), _pmap)
-    @test prob.ps[p] ≈ 3.0
-    @test prob.f.initializeprob === nothing
-    # Default overridden by ODEProblem, guess provided
+    test_parameter(prob, p, 3pmap[q])
+
+    # Should not be solved for:
+
+    # Override dependent default with direct value
     @mtkbuild sys = ODESystem(
         [D(x) ~ q * x, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 1.0])
+    _pmap = merge(pmap, Dict(p => 1.0))
     prob = ODEProblem(sys, u0map, (0.0, 1.0), _pmap)
-    @test prob.ps[p] ≈ 3.0
+    @test prob.ps[p] ≈ 1.0
+    @test prob.f.initializeprob === nothing
+
+    # Non-floating point
+    @parameters r::Int s::Int
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ s * x, D(y) ~ y * r], t; defaults = [s => 2r], guesses = [s => 1.0])
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), [r => 1])
+    @test prob.ps[r] == 1
+    @test prob.ps[s] == 2
     @test prob.f.initializeprob === nothing
 
     @mtkbuild sys = ODESystem([D(x) ~ x, p ~ x + y], t; guesses = [p => 0.0])

From 62c2f09d11c72dfeaf212270640fc7ff84000455 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 7 Oct 2024 22:14:57 +0530
Subject: [PATCH 16/19] fix: only allow numeric pdeps to be initialization
 equations

---
 src/systems/nonlinear/initializesystem.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 6a21e9a6f0..1641d66dc0 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -132,6 +132,7 @@ function generate_initializesystem(sys::ODESystem;
 
     # 5) parameter dependencies become equations, their LHS become unknowns
     for eq in parameter_dependencies(sys)
+        is_variable_floatingpoint(eq.lhs) || continue
         varp = tovar(eq.lhs)
         paramsubs[eq.lhs] = varp
         push!(eqs_ics, eq)

From ba82445d0de255f6340348c019c675a47a65609f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 7 Oct 2024 22:15:11 +0530
Subject: [PATCH 17/19] feat: allow `nothing` to override retained values in
 `remake_initializeprob`

---
 src/systems/nonlinear/initializesystem.jl | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 1641d66dc0..38a3095b5c 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -246,6 +246,11 @@ function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p)
     u0 = todict(u0)
     defs = defaults(sys)
     varmap = merge(defs, u0)
+    for k in collect(keys(varmap))
+        if varmap[k] === nothing
+            delete!(varmap, k)
+        end
+    end
     varmap = canonicalize_varmap(varmap)
     missingvars = setdiff(unknowns(sys), collect(keys(varmap)))
     setobserved = filter(keys(varmap)) do var
@@ -272,6 +277,16 @@ function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p)
                 p = merge(meta.pmap, p)
             end
         end
+        for k in collect(keys(u0))
+            if u0[k] === nothing
+                delete!(u0, k)
+            end
+        end
+        for k in collect(keys(p))
+            if p[k] === nothing
+                delete!(p, k)
+            end
+        end
 
         initprob = InitializationProblem(sys, t0, u0, p)
         initprobmap = getu(initprob, unknowns(sys))

From ec75ebee36aff1aefcd8dcc2bfc0ac6d3e196649 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 6 Aug 2024 18:45:41 +0530
Subject: [PATCH 18/19] docs: document parameter initialization in the tutorial

---
 docs/src/tutorials/initialization.md | 134 +++++++++++++++++++++++++++
 1 file changed, 134 insertions(+)

diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md
index f40c28991f..6451832dbc 100644
--- a/docs/src/tutorials/initialization.md
+++ b/docs/src/tutorials/initialization.md
@@ -201,6 +201,87 @@ long enough you will see that `λ = 0` is required for this equation, but since
     problem constructor. Additionally, any warning about not being fully determined can
     be suppressed via passing `warn_initialize_determined = false`.
 
+## Initialization of parameters
+
+Parameters may also be treated as unknowns in the initialization system. Doing so works
+almost identically to the standard case. For a parameter to be an initialization unknown
+(henceforth referred to as "solved parameter") it must represent a floating point number
+(have a `symtype` of `Real` or `<:AbstractFloat`) or an array of such numbers. Additionally,
+it must have a guess and one of the following conditions must be satisfied:
+
+ 1. The value of the parameter as passed to `ODEProblem` is an expression involving other
+    variables/parameters. For example, if `[p => 2q + x]` is passed to `ODEProblem`. In
+    this case, `p ~ 2q + x` is used as an equation during initialization.
+ 2. The parameter has a default (and no value for it is given to `ODEProblem`, since
+    that is condition 1). The default will be used as an equation during initialization.
+ 3. The parameter has a default of `missing`. If `ODEProblem` is given a value for this
+    parameter, it is used as an equation during initialization (whether the value is an
+    expression or not).
+ 4. `ODEProblem` is given a value of `missing` for the parameter. If the parameter has a
+    default, it will be used as an equation during initialization.
+
+All parameter dependencies (where the dependent parameter is a floating point number or
+array thereof) also become equations during initialization, and the dependent parameters
+become unknowns.
+
+`remake` will reconstruct the initialization system and problem, given the new
+constraints provided to it. The new values will be combined with the original
+variable-value mapping provided to `ODEProblem` and used to construct the initialization
+problem.
+
+### Parameter initialization by example
+
+Consider the following system, where the sum of two unknowns is a constant parameter
+`total`.
+
+```@example paraminit
+using ModelingToolkit, OrdinaryDiffEq # hidden
+using ModelingToolkit: t_nounits as t, D_nounits as D # hidden
+
+@variables x(t) y(t)
+@parameters total
+@mtkbuild sys = ODESystem([D(x) ~ -x, total ~ x + y], t;
+    defaults = [total => missing], guesses = [total => 1.0])
+```
+
+Given any two of `x`, `y` and `total` we can determine the remaining variable.
+
+```@example paraminit
+prob = ODEProblem(sys, [x => 1.0, y => 2.0], (0.0, 1.0))
+integ = init(prob, Tsit5())
+@assert integ.ps[total] ≈ 3.0 # hide
+integ.ps[total]
+```
+
+Suppose we want to re-create this problem, but now solve for `x` given `total` and `y`:
+
+```@example paraminit
+prob2 = remake(prob; u0 = [y => 1.0], p = [total => 4.0])
+initsys = prob2.f.initializeprob.f.sys
+```
+
+The system is now overdetermined. In fact:
+
+```@example paraminit
+[equations(initsys); observed(initsys)]
+```
+
+The system can never be satisfied and will always lead to an `InitialFailure`. This is
+due to the aforementioned behavior of retaining the original variable-value mapping
+provided to `ODEProblem`. To fix this, we pass `x => nothing` to `remake` to remove its
+retained value.
+
+```@example paraminit
+prob2 = remake(prob; u0 = [y => 1.0, x => nothing], p = [total => 4.0])
+initsys = prob2.f.initializeprob.f.sys
+```
+
+The system is fully determined, and the equations are solvable.
+
+```@example
+[equations(initsys); observed(initsys)]
+```
+
 ## Diving Deeper: Constructing the Initialization System
 
 To get a better sense of the initialization system and to help debug it, you can construct
@@ -383,3 +464,56 @@ sol[α * x - β * x * y]
 ```@example init
 plot(sol)
 ```
+
+## Solving for parameters during initialization
+
+Sometimes, it is necessary to solve for a parameter during initialization. For example,
+given a spring-mass system we want to find the un-stretched length of the spring given
+that the initial condition of the system is its steady state.
+
+```@example init
+using ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica: Fixed, Mass, Spring,
+                                                                       Force, Damper
+using ModelingToolkitStandardLibrary.Blocks: Constant
+
+@named mass = Mass(; m = 1.0, s = 1.0, v = 0.0, a = 0.0)
+@named fixed = Fixed(; s0 = 0.0)
+@named spring = Spring(; c = 2.0, s_rel0 = missing)
+@named gravity = Force()
+@named constant = Constant(; k = 9.81)
+@named damper = Damper(; d = 0.1)
+@mtkbuild sys = ODESystem(
+    [connect(fixed.flange, spring.flange_a), connect(spring.flange_b, mass.flange_a),
+        connect(mass.flange_a, gravity.flange), connect(constant.output, gravity.f),
+        connect(fixed.flange, damper.flange_a), connect(damper.flange_b, mass.flange_a)],
+    t;
+    systems = [fixed, spring, mass, gravity, constant, damper],
+    guesses = [spring.s_rel0 => 1.0])
+```
+
+Note that we explicitly provide `s_rel0 = missing` to the spring. Parameters are only
+solved for during initialization if their value (either default, or explicitly passed
+to the `ODEProblem` constructor) is `missing`. We also need to provide a guess for the
+parameter.
+
+If a parameter is not given a value of `missing`, and does not have a default or initial
+value, the `ODEProblem` constructor will throw an error. If the parameter _does_ have a
+value of `missing`, it must be given a guess.
+
+```@example init
+prob = ODEProblem(sys, [], (0.0, 1.0))
+prob.ps[spring.s_rel0]
+```
+
+Note that the value of the parameter in the problem is zero, similar to unknowns that
+are solved for during initialization.
+
+```@example init
+integ = init(prob)
+integ.ps[spring.s_rel0]
+```
+
+The un-stretched length of the spring is now correctly calculated. The same result can be
+achieved if `s_rel0 = missing` is omitted when constructing `spring`, and instead
+`spring.s_rel0 => missing` is passed to the `ODEProblem` constructor along with values
+of other parameters.

From 84a1f2e233624a4fd2a60420d03f6349e3c74aa6 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Oct 2024 17:54:17 +0530
Subject: [PATCH 19/19] build: bump compats for SciMLBase, DiffEqBase,
 OrdinaryDiffEqCore

---
 Project.toml | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 1947ea7ce6..cfd424ef8b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -81,7 +81,7 @@ ConstructionBase = "1"
 DataInterpolations = "6.4"
 DataStructures = "0.17, 0.18"
 DeepDiffs = "1"
-DiffEqBase = "6.103.0"
+DiffEqBase = "6.157"
 DiffEqCallbacks = "2.16, 3, 4"
 DiffEqNoiseProcess = "5"
 DiffRules = "0.1, 1.0"
@@ -110,12 +110,13 @@ NonlinearSolve = "3.14"
 OffsetArrays = "1"
 OrderedCollections = "1"
 OrdinaryDiffEq = "6.82.0"
+OrdinaryDiffEqCore = "1.7.0"
 PrecompileTools = "1"
 REPL = "1"
 RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
-SciMLBase = "2.55"
+SciMLBase = "2.56.1"
 SciMLStructures = "1.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"
@@ -148,6 +149,7 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -162,4 +164,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"]
+test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"]