From 16a32884c622f2c93908ee0b53392b61a13f1b74 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Mon, 15 Jul 2024 18:48:46 +0200 Subject: [PATCH 1/4] Expand observed equations into lhs => rhs defaults, but flip if lhs is given --- src/systems/diffeqs/abstractodesystem.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 6c54cf2200..b03819a633 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -725,10 +725,16 @@ function get_u0( defs = mergedefaults(defs, parammap, ps) end - obs = filter!(x -> !(x[1] isa Number), - map(x -> isparameter(x.rhs) ? x.lhs => x.rhs : x.rhs => x.lhs, observed(sys))) - observedmap = isempty(obs) ? Dict() : todict(obs) - defs = mergedefaults(defs, observedmap, u0map, dvs) + # Convert observed equations "lhs ~ rhs" into defaults. + # Use the order "lhs => rhs" by default, but flip it to "rhs => lhs" + # if "lhs" is known by other means (parameter, another default, ...) + # TODO: Is there a better way to determine which equations to flip? + obs = map(x -> x.lhs => x.rhs, observed(sys)) + obs = map(x -> isparameter(x[1]) || x[1] in keys(defs) ? reverse(x) : x, obs) + obs = filter!(x -> !(x[1] isa Number), obs) # exclude e.g. "0 => x^2 + y^2 - 25" + obsmap = isempty(obs) ? Dict() : todict(obs) + + defs = mergedefaults(defs, obsmap, u0map, dvs) if symbolic_u0 u0 = varmap_to_vars( u0map, dvs; defaults = defs, tofloat = false, use_union = false, toterm) From 6e636b2cdd6305a6d24e94294bce720346431c99 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Mon, 15 Jul 2024 18:56:39 +0200 Subject: [PATCH 2/4] Activate unactivated tests --- test/guess_propagation.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/guess_propagation.jl b/test/guess_propagation.jl index 6d71ce6ca1..738e930adc 100644 --- a/test/guess_propagation.jl +++ b/test/guess_propagation.jl @@ -90,21 +90,21 @@ prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) @variables y(t) = x0 @mtkbuild sys = ODESystem([x ~ x0, D(y) ~ x], t) prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) -prob[x] == 1.0 -prob[y] == 1.0 +@test prob[x] == 1.0 +@test prob[y] == 1.0 @parameters x0 @variables x(t) @variables y(t) = x0 @mtkbuild sys = ODESystem([x ~ y, D(y) ~ x], t) prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) -prob[x] == 1.0 -prob[y] == 1.0 +@test prob[x] == 1.0 +@test prob[y] == 1.0 @parameters x0 @variables x(t) = x0 @variables y(t) = x @mtkbuild sys = ODESystem([x ~ y, D(y) ~ x], t) prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) -prob[x] == 1.0 -prob[y] == 1.0 +@test prob[x] == 1.0 +@test prob[y] == 1.0 From 5e279cee3564f5821e92c273cd53c05eca27ca88 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Mon, 15 Jul 2024 19:04:28 +0200 Subject: [PATCH 3/4] Added test from issue (and a simplified version of it) --- test/odesystem.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/test/odesystem.jl b/test/odesystem.jl index ab28f5cb47..c7da5ba702 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -1194,3 +1194,23 @@ end @test_nowarn obsfn(buffer, [1.0], ps..., 3.0) @test buffer ≈ [2.0, 3.0, 4.0] end + +# https://github.com/SciML/ModelingToolkit.jl/issues/2859 +@testset "Initialization with defaults from observed equations (edge case)" begin + @variables x(t) y(t) z(t) + eqs = [D(x) ~ 0, y ~ x, D(z) ~ 0] + defaults = [x => 1, z => y] + @named sys = ODESystem(eqs, t; defaults) + ssys = structural_simplify(sys) + prob = ODEProblem(ssys, [], (0.0, 1.0), []) + @test prob[x] == prob[y] == prob[z] == 1.0 + + @parameters y0 + @variables x(t) y(t) z(t) + eqs = [D(x) ~ 0, y ~ y0 / x, D(z) ~ y] + defaults = [y0 => 1, x => 1, z => y] + @named sys = ODESystem(eqs, t; defaults) + ssys = structural_simplify(sys) + prob = ODEProblem(ssys, [], (0.0, 1.0), []) + @test prob[x] == prob[y] == prob[z] == 1.0 +end From 8862b97e86b8f677bb641f2c0ab8a2f8d1172e8e Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Mon, 15 Jul 2024 20:15:15 +0200 Subject: [PATCH 4/4] Remove redundant isparameter(lhs) check --- src/systems/diffeqs/abstractodesystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index b03819a633..ab0565eb31 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -730,7 +730,7 @@ function get_u0( # if "lhs" is known by other means (parameter, another default, ...) # TODO: Is there a better way to determine which equations to flip? obs = map(x -> x.lhs => x.rhs, observed(sys)) - obs = map(x -> isparameter(x[1]) || x[1] in keys(defs) ? reverse(x) : x, obs) + obs = map(x -> x[1] in keys(defs) ? reverse(x) : x, obs) obs = filter!(x -> !(x[1] isa Number), obs) # exclude e.g. "0 => x^2 + y^2 - 25" obsmap = isempty(obs) ? Dict() : todict(obs)