Merge pull request #2575 from AayushSabharwal/as/initialization-fix

fix: fix initialization and observed bugs

Author: ChrisRackauckas

src/systems/diffeqs/odesystem.jl

@@ -460,9 +460,6 @@ function build_explicit_observed_function(sys, ts;
     else
         ps = (DestructuredArgs(ps, inbounds = !checkbounds),)
     end
-    if isempty(ps)
-        ps = (DestructuredArgs([]),)
-    end
     dvs = DestructuredArgs(unknowns(sys), inbounds = !checkbounds)
     if inputs === nothing
         args = [dvs, ps..., ivs...]

src/systems/parameter_buffer.jl

@@ -39,7 +39,7 @@ function MTKParameters(
         u0 = Dict()
     end
     defs = merge(defs, u0)
-    defs = merge(defs, Dict(eq.lhs => eq.rhs for eq in observed(sys)))
+    defs = merge(Dict(eq.lhs => eq.rhs for eq in observed(sys)), defs)
     p = merge(defs, p)
     p = merge(Dict(unwrap(k) => v for (k, v) in p),
         Dict(default_toterm(unwrap(k)) => v for (k, v) in p))

src/variables.jl

@@ -192,11 +192,11 @@ function _varmap_to_vars(varmap::Dict, varlist; defaults = Dict(), check = false
         toterm = Symbolics.diff2term, initialization_phase = false)
     varmap = canonicalize_varmap(varmap; toterm)
     defaults = canonicalize_varmap(defaults; toterm)
+    varmap = merge(defaults, varmap)
     values = Dict()
     for var in varlist
         var = unwrap(var)
-        val = unwrap(fixpoint_sub(fixpoint_sub(var, varmap; operator = Symbolics.Operator),
-            defaults; operator = Symbolics.Operator))
+        val = unwrap(fixpoint_sub(var, varmap; operator = Symbolics.Operator))
         if symbolic_type(val) === NotSymbolic()
             values[var] = val
         end
@@ -211,6 +211,11 @@ function canonicalize_varmap(varmap; toterm = Symbolics.diff2term)
     for (k, v) in varmap
         new_varmap[unwrap(k)] = unwrap(v)
         new_varmap[toterm(unwrap(k))] = unwrap(v)
+        if Symbolics.isarraysymbolic(k) && Symbolics.shape(k) !== Symbolics.Unknown()
+            for i in eachindex(k)
+                new_varmap[k[i]] = v[i]
+            end
+        end
     end
     return new_varmap
 end

test/initial_values.jl

@@ -33,3 +33,30 @@ getter = getu(sys, [x..., y, z...])
 @test getter(get_u0(
     sys, [y => 2w, w => 3.0, z[1] => 2p1, z[2] => 3p2], [p1 => 3.0, p2 => 4.0])[1]) ==
       [1.0, 2.0, 3.0, 6.0, 6.0, 12.0]
+
+# Issue#2566
+@variables X(t)
+@parameters p1 p2 p3
+
+p_vals = [p1 => 1.0, p2 => 2.0]
+u_vals = [X => 3.0]
+
+var_vals = [p1 => 1.0, p2 => 2.0, X => 3.0]
+desired_values = [p1, p2, p3]
+defaults = Dict([p3 => X])
+vals = ModelingToolkit.varmap_to_vars(var_vals, desired_values; defaults = defaults)
+@test vals == [1.0, 2.0, 3.0]
+
+# Issue#2565
+# Create ODESystem.
+@variables X1(t) X2(t)
+@parameters k1 k2 Γ[1:1]=X1 + X2
+eq = D(X1) ~ -k1 * X1 + k2 * (-X1 + Γ[1])
+obs = X2 ~ Γ[1] - X1
+@mtkbuild osys_m = ODESystem([eq], t, [X1], [k1, k2, Γ[1]]; observed = [X2 ~ Γ[1] - X1])
+
+# Creates ODEProblem.
+u0 = [X1 => 1.0, X2 => 2.0]
+tspan = (0.0, 1.0)
+ps = [k1 => 1.0, k2 => 5.0]
+@test_nowarn oprob = ODEProblem(osys_m, u0, tspan, ps)

test/nonlinearsystem.jl

@@ -242,3 +242,18 @@ testdict = Dict([:test => 1])
     @test prob_.u0 == [1.0, 2.0, 1.0]
     @test prob_.p == MTKParameters(sys, [a => 2.0, b => 1.0, c => 1.0])
 end
+
+@testset "Observed function generation without parameters" begin
+    @variables x y z
+
+    eqs = [0 ~ x + sin(y),
+        0 ~ z - cos(x),
+        0 ~ x * y]
+    @named ns = NonlinearSystem(eqs, [x, y, z], [])
+    ns = complete(ns)
+    vs = [unknowns(ns); parameters(ns)]
+    ss_mtk = structural_simplify(ns)
+    prob = NonlinearProblem(ss_mtk, vs .=> 1.0)
+    sol = solve(prob)
+    @test_nowarn sol[unknowns(ns)]
+end