fix: handle Initial parameters in code generation

Author: AayushSabharwal

src/inputoutput.jl

@@ -211,7 +211,7 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
     sys, _ = io_preprocessing(sys, inputs, []; simplify, kwargs...)
 
     dvs = unknowns(sys)
-    ps = parameters(sys)
+    ps = parameters(sys; initial_parameters = true)
     ps = setdiff(ps, inputs)
     if disturbance_inputs !== nothing
         # remove from inputs since we do not want them as actual inputs to the dynamics
@@ -234,16 +234,14 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
            [eq.rhs for eq in eqs]
 
     # TODO: add an optional check on the ordering of observed equations
-    u = map(x -> time_varying_as_func(value(x), sys), dvs)
-    p = map(x -> time_varying_as_func(value(x), sys), ps)
-    p = reorder_parameters(sys, p)
+    p = reorder_parameters(sys, ps)
     t = get_iv(sys)
 
     # pre = has_difference ? (ex -> ex) : get_postprocess_fbody(sys)
     if disturbance_argument
-        args = (u, inputs, p..., t, disturbance_inputs)
+        args = (dvs, inputs, p..., t, disturbance_inputs)
     else
-        args = (u, inputs, p..., t)
+        args = (dvs, inputs, p..., t)
     end
     if implicit_dae
         ddvs = map(Differential(get_iv(sys)), dvs)

src/systems/abstractsystem.jl

@@ -161,7 +161,7 @@ time-independent systems. If `split=true` (the default) was passed to [`complete
 object.
 """
 function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys),
-        ps = parameters(sys);
+        ps = parameters(sys; initial_parameters = true);
         expression = Val{true}, eval_expression = false, eval_module = @__MODULE__,
         cachesyms::Tuple = (), kwargs...)
     if !iscomplete(sys)
@@ -533,7 +533,7 @@ function SymbolicIndexingInterface.get_all_timeseries_indexes(sys::AbstractSyste
 end
 
 function SymbolicIndexingInterface.parameter_symbols(sys::AbstractSystem)
-    return parameters(sys)
+    return parameters(sys; initial_parameters = true)
 end
 
 function SymbolicIndexingInterface.is_independent_variable(sys::AbstractSystem, sym)
@@ -2391,7 +2391,7 @@ function linearize_symbolic(sys::AbstractSystem, inputs,
         kwargs...)
     sts = unknowns(sys)
     t = get_iv(sys)
-    ps = parameters(sys)
+    ps = parameters(sys; initial_parameters = true)
     p = reorder_parameters(sys, ps)
 
     fun_expr = generate_function(sys, sts, ps; expression = Val{true})[1]

src/systems/callbacks.jl

@@ -682,7 +682,7 @@ function compile_affect(eqs::Vector{Equation}, cb, sys, dvs, ps; outputidxs = no
 end
 
 function generate_rootfinding_callback(sys::AbstractTimeDependentSystem,
-        dvs = unknowns(sys), ps = parameters(sys); kwargs...)
+        dvs = unknowns(sys), ps = parameters(sys; initial_parameters = true); kwargs...)
     cbs = continuous_events(sys)
     isempty(cbs) && return nothing
     generate_rootfinding_callback(cbs, sys, dvs, ps; kwargs...)
@@ -693,7 +693,7 @@ generate_rootfinding_callback and thus we can produce a ContinuousCallback inste
 """
 function generate_single_rootfinding_callback(
         eq, cb, sys::AbstractTimeDependentSystem, dvs = unknowns(sys),
-        ps = parameters(sys); kwargs...)
+        ps = parameters(sys; initial_parameters = true); kwargs...)
     if !isequal(eq.lhs, 0)
         eq = 0 ~ eq.lhs - eq.rhs
     end
@@ -736,7 +736,7 @@ end
 
 function generate_vector_rootfinding_callback(
         cbs, sys::AbstractTimeDependentSystem, dvs = unknowns(sys),
-        ps = parameters(sys); rootfind = SciMLBase.RightRootFind,
+        ps = parameters(sys; initial_parameters = true); rootfind = SciMLBase.RightRootFind,
         reinitialization = SciMLBase.CheckInit(), kwargs...)
     eqs = map(cb -> flatten_equations(cb.eqs), cbs)
     num_eqs = length.(eqs)
@@ -861,7 +861,7 @@ function compile_affect_fn(cb, sys::AbstractTimeDependentSystem, dvs, ps, kwargs
 end
 
 function generate_rootfinding_callback(cbs, sys::AbstractTimeDependentSystem,
-        dvs = unknowns(sys), ps = parameters(sys); kwargs...)
+        dvs = unknowns(sys), ps = parameters(sys; initial_parameters = true); kwargs...)
     eqs = map(cb -> flatten_equations(cb.eqs), cbs)
     num_eqs = length.(eqs)
     total_eqs = sum(num_eqs)
@@ -949,7 +949,9 @@ function invalid_variables(sys, expr)
 end
 function unassignable_variables(sys, expr)
     assignable_syms = reduce(
-        vcat, Symbolics.scalarize.(vcat(unknowns(sys), parameters(sys))); init = [])
+        vcat, Symbolics.scalarize.(vcat(
+            unknowns(sys), parameters(sys; initial_parameters = true)));
+        init = [])
     written = reduce(vcat, Symbolics.scalarize.(vars(expr)); init = [])
     return filter(
         x -> !any(isequal(x), assignable_syms), written)
@@ -1075,7 +1077,7 @@ function generate_discrete_callback(cb, sys, dvs, ps; postprocess_affect_expr! =
 end
 
 function generate_discrete_callbacks(sys::AbstractSystem, dvs = unknowns(sys),
-        ps = parameters(sys); kwargs...)
+        ps = parameters(sys; initial_parameters = true); kwargs...)
     has_discrete_events(sys) || return nothing
     symcbs = discrete_events(sys)
     isempty(symcbs) && return nothing

src/systems/diffeqs/abstractodesystem.jl

@@ -105,7 +105,8 @@ function calculate_control_jacobian(sys::AbstractODESystem;
 end
 
 function generate_tgrad(
-        sys::AbstractODESystem, dvs = unknowns(sys), ps = parameters(sys);
+        sys::AbstractODESystem, dvs = unknowns(sys), ps = parameters(
+            sys; initial_parameters = true);
         simplify = false, kwargs...)
     tgrad = calculate_tgrad(sys, simplify = simplify)
     p = reorder_parameters(sys, ps)
@@ -117,7 +118,7 @@ function generate_tgrad(
 end
 
 function generate_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = parameters(sys);
+        ps = parameters(sys; initial_parameters = true);
         simplify = false, sparse = false, kwargs...)
     jac = calculate_jacobian(sys; simplify = simplify, sparse = sparse)
     p = reorder_parameters(sys, ps)
@@ -129,15 +130,15 @@ function generate_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
 end
 
 function generate_control_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = parameters(sys);
+        ps = parameters(sys; initial_parameters = true);
         simplify = false, sparse = false, kwargs...)
     jac = calculate_control_jacobian(sys; simplify = simplify, sparse = sparse)
     p = reorder_parameters(sys, ps)
     return build_function_wrapper(sys, jac, dvs, p..., get_iv(sys); kwargs...)
 end
 
 function generate_dae_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = parameters(sys); simplify = false, sparse = false,
+        ps = parameters(sys; initial_parameters = true); simplify = false, sparse = false,
         kwargs...)
     jac_u = calculate_jacobian(sys; simplify = simplify, sparse = sparse)
     derivatives = Differential(get_iv(sys)).(unknowns(sys))
@@ -153,7 +154,7 @@ function generate_dae_jacobian(sys::AbstractODESystem, dvs = unknowns(sys),
 end
 
 function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
-        ps = parameters(sys);
+        ps = parameters(sys; initial_parameters = true);
         implicit_dae = false,
         ddvs = implicit_dae ? map(Differential(get_iv(sys)), dvs) :
                nothing,
@@ -691,7 +692,7 @@ function SymbolicTstops(
             term(:, t0, unwrap(val), t1; type = AbstractArray{Real})
         end
     end
-    rps = reorder_parameters(sys, parameters(sys))
+    rps = reorder_parameters(sys)
     tstops, _ = build_function_wrapper(sys, tstops,
         rps...,
         t0,
@@ -817,7 +818,7 @@ function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan
 end
 
 function generate_history(sys::AbstractODESystem, u0; expression = Val{false}, kwargs...)
-    p = reorder_parameters(sys, parameters(sys))
+    p = reorder_parameters(sys)
     build_function_wrapper(
         sys, u0, p..., get_iv(sys); expression, p_start = 1, p_end = length(p),
         similarto = typeof(u0), wrap_delays = false, kwargs...)

src/systems/diffeqs/odesystem.jl

@@ -424,7 +424,7 @@ function build_explicit_observed_function(sys, ts;
         eval_module = @__MODULE__,
         output_type = Array,
         checkbounds = true,
-        ps = parameters(sys),
+        ps = parameters(sys; initial_parameters = true),
         return_inplace = false,
         param_only = false,
         op = Operator,

src/systems/diffeqs/sdesystem.jl

@@ -409,7 +409,7 @@ function __get_num_diag_noise(mat)
 end
 
 function generate_diffusion_function(sys::SDESystem, dvs = unknowns(sys),
-        ps = parameters(sys); isdde = false, kwargs...)
+        ps = parameters(sys; initial_parameters = true); isdde = false, kwargs...)
     eqs = get_noiseeqs(sys)
     p = reorder_parameters(sys, ps)
     return build_function_wrapper(sys, eqs, dvs, p..., get_iv(sys); kwargs...)

src/systems/jumps/jumpsystem.jl

@@ -282,7 +282,7 @@ function generate_rate_function(js::JumpSystem, rate)
         csubs = Dict(c => getdefault(c) for c in consts)
         rate = substitute(rate, csubs)
     end
-    p = reorder_parameters(js, parameters(js))
+    p = reorder_parameters(js)
     rf = build_function_wrapper(js, rate, unknowns(js), p...,
         get_iv(js),
         expression = Val{true})
@@ -634,7 +634,7 @@ end
 function JumpSysMajParamMapper(js::JumpSystem, p; jseqs = nothing, rateconsttype = Float64)
     eqs = (jseqs === nothing) ? equations(js) : jseqs
     paramexprs = [maj.scaled_rates for maj in eqs.x[1]]
-    psyms = reduce(vcat, reorder_parameters(js, parameters(js)); init = [])
+    psyms = reduce(vcat, reorder_parameters(js); init = [])
     paramdict = Dict(value(k) => value(v) for (k, v) in zip(psyms, vcat(p...)))
     JumpSysMajParamMapper{typeof(paramexprs), typeof(psyms), rateconsttype}(paramexprs,
         psyms,

src/systems/nonlinear/initializesystem.jl

@@ -339,7 +339,9 @@ end
 
 function ReconstructInitializeprob(
         srcsys::AbstractSystem, dstsys::AbstractSystem; remap = Dict())
-    syms = reduce(vcat, reorder_parameters(dstsys, parameters(dstsys)); init = [])
+    syms = reduce(
+        vcat, reorder_parameters(dstsys, parameters(dstsys; initial_parameters = true));
+        init = [])
     getter = getu(srcsys, map(x -> get(remap, x, x), syms))
     setter = setp_oop(dstsys, syms)
     return ReconstructInitializeprob(getter, setter)

src/systems/nonlinear/nonlinearsystem.jl

@@ -248,7 +248,8 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal
 end
 
 function generate_jacobian(
-        sys::NonlinearSystem, vs = unknowns(sys), ps = parameters(sys);
+        sys::NonlinearSystem, vs = unknowns(sys), ps = parameters(
+            sys; initial_parameters = true);
         sparse = false, simplify = false, kwargs...)
     jac = calculate_jacobian(sys, sparse = sparse, simplify = simplify)
     p = reorder_parameters(sys, ps)
@@ -268,15 +269,17 @@ function calculate_hessian(sys::NonlinearSystem; sparse = false, simplify = fals
 end
 
 function generate_hessian(
-        sys::NonlinearSystem, vs = unknowns(sys), ps = parameters(sys);
+        sys::NonlinearSystem, vs = unknowns(sys), ps = parameters(
+            sys; initial_parameters = true);
         sparse = false, simplify = false, kwargs...)
     hess = calculate_hessian(sys, sparse = sparse, simplify = simplify)
     p = reorder_parameters(sys, ps)
     return build_function_wrapper(sys, hess, vs, p...; kwargs...)
 end
 
 function generate_function(
-        sys::NonlinearSystem, dvs = unknowns(sys), ps = parameters(sys);
+        sys::NonlinearSystem, dvs = unknowns(sys), ps = parameters(
+            sys; initial_parameters = true);
         scalar = false, kwargs...)
     rhss = [deq.rhs for deq in equations(sys)]
     dvs′ = value.(dvs)
@@ -569,7 +572,7 @@ end
 function CacheWriter(sys::AbstractSystem, buffer_types::Vector{TypeT},
         exprs::Dict{TypeT, Vector{Any}}, solsyms, obseqs::Vector{Equation};
         eval_expression = false, eval_module = @__MODULE__)
-    ps = parameters(sys)
+    ps = parameters(sys; initial_parameters = true)
     rps = reorder_parameters(sys, ps)
     obs_assigns = [eq.lhs ← eq.rhs for eq in obseqs]
     body = map(eachindex(buffer_types), buffer_types) do i, T
@@ -589,7 +592,7 @@ struct SCCNonlinearFunction{iip} end
 function SCCNonlinearFunction{iip}(
         sys::NonlinearSystem, _eqs, _dvs, _obs, cachesyms; eval_expression = false,
         eval_module = @__MODULE__, kwargs...) where {iip}
-    ps = parameters(sys)
+    ps = parameters(sys; initial_parameters = true)
     rps = reorder_parameters(sys, ps)
 
     obs_assignments = [eq.lhs ← eq.rhs for eq in _obs]

src/systems/optimization/constraints_system.jl

@@ -171,7 +171,8 @@ function calculate_jacobian(sys::ConstraintsSystem; sparse = false, simplify = f
 end
 
 function generate_jacobian(
-        sys::ConstraintsSystem, vs = unknowns(sys), ps = parameters(sys);
+        sys::ConstraintsSystem, vs = unknowns(sys), ps = parameters(
+            sys; initial_parameters = true);
         sparse = false, simplify = false, kwargs...)
     jac = calculate_jacobian(sys, sparse = sparse, simplify = simplify)
     p = reorder_parameters(sys, ps)
@@ -190,15 +191,16 @@ function calculate_hessian(sys::ConstraintsSystem; sparse = false, simplify = fa
 end
 
 function generate_hessian(
-        sys::ConstraintsSystem, vs = unknowns(sys), ps = parameters(sys);
+        sys::ConstraintsSystem, vs = unknowns(sys), ps = parameters(
+            sys; initial_parameters = true);
         sparse = false, simplify = false, kwargs...)
     hess = calculate_hessian(sys, sparse = sparse, simplify = simplify)
     p = reorder_parameters(sys, ps)
     return build_function_wrapper(sys, hess, vs, p...; kwargs...)
 end
 
 function generate_function(sys::ConstraintsSystem, dvs = unknowns(sys),
-        ps = parameters(sys);
+        ps = parameters(sys; initial_parameters = true);
         kwargs...)
     lhss = generate_canonical_form_lhss(sys)
     p = reorder_parameters(sys, value.(ps))