From f6aa73f2090a05ca3af7eb3d0d159561576696dc Mon Sep 17 00:00:00 2001 From: DhairyaLGandhi Date: Wed, 11 Jun 2025 17:23:20 +0530 Subject: [PATCH 01/10] test: update for MTK v10 --- test/downstream/observables_autodiff.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/test/downstream/observables_autodiff.jl b/test/downstream/observables_autodiff.jl index ed9ab796a..0e236800a 100644 --- a/test/downstream/observables_autodiff.jl +++ b/test/downstream/observables_autodiff.jl @@ -33,7 +33,7 @@ sol = solve(prob, Tsit5()) gs, = gradient(sol) do sol sum(sol[sys.w]) end - du_ = [0.0, 1.0, 1.0, 1.0] + du_ = [1.0, 1.0, 1.0, 0.0] du = [du_ for _ in sol.u] @test du == gs.u @@ -41,7 +41,7 @@ sol = solve(prob, Tsit5()) gs, = gradient(sol) do sol sum(sum.(sol[[sys.w, sys.x]])) end - du_ = [0.0, 1.0, 1.0, 2.0] + du_ = [1.0, 1.0, 2.0, 0.0] du = [du_ for _ in sol.u] @test du == gs.u end @@ -118,14 +118,19 @@ end end @testset "Adjoints with DAE" begin - gs_mtkp, gs_p_new = gradient(prob.p, prob.p.tunable) do p, new_tunables + model = create_model() + sys = mtkcompile(model) + prob = ODEProblem(sys, [], (0.0, 1.0)) + tunables, _, _ = SS.canonicalize(SS.Tunable(), prob.p) + + gs_mtkp, gs_p_new = gradient(prob.p, tunables) do p, new_tunables new_p = SS.replace(SS.Tunable(), p, new_tunables) new_prob = remake(prob, p = new_p) sol = solve(new_prob, Rodas4()) - mean(abs.(sol[sys.ampermeter.i] .- gt)) sum(sol[sys.ampermeter.i]) end @test isnothing(gs_mtkp) - @test length(gs_p_new) == length(p_new) + @test !isnothing(gs_p_new) + @test length(gs_p_new) == length(tunables) end From f99fd5918ccbc3059d07dff85ab4dbd386209397 Mon Sep 17 00:00:00 2001 From: Dhairya Gandhi Date: Wed, 11 Jun 2025 18:36:31 +0530 Subject: [PATCH 02/10] Update test/downstream/observables_autodiff.jl Co-authored-by: Christopher Rackauckas --- test/downstream/observables_autodiff.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/downstream/observables_autodiff.jl b/test/downstream/observables_autodiff.jl index 0e236800a..d5921bf07 100644 --- a/test/downstream/observables_autodiff.jl +++ b/test/downstream/observables_autodiff.jl @@ -42,7 +42,7 @@ sol = solve(prob, Tsit5()) sum(sum.(sol[[sys.w, sys.x]])) end du_ = [1.0, 1.0, 2.0, 0.0] - du = [du_ for _ in sol.u] + du = [du_ for _ in sol[[D(x), x, y, z]] @test du == gs.u end From cc9a7d53bb85abb0088f0f12ad26ea6e65e87e38 Mon Sep 17 00:00:00 2001 From: DhairyaLGandhi Date: Wed, 11 Jun 2025 18:40:27 +0530 Subject: [PATCH 03/10] test: syntax --- test/downstream/observables_autodiff.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/downstream/observables_autodiff.jl b/test/downstream/observables_autodiff.jl index d5921bf07..54bea8c99 100644 --- a/test/downstream/observables_autodiff.jl +++ b/test/downstream/observables_autodiff.jl @@ -42,7 +42,7 @@ sol = solve(prob, Tsit5()) sum(sum.(sol[[sys.w, sys.x]])) end du_ = [1.0, 1.0, 2.0, 0.0] - du = [du_ for _ in sol[[D(x), x, y, z]] + du = [du_ for _ in sol[[D(x), x, y, z]]] @test du == gs.u end From 8850e1df59fa883ea910f5c3dfaaec800290c136 Mon Sep 17 00:00:00 2001 From: DhairyaLGandhi Date: Wed, 11 Jun 2025 20:52:26 +0530 Subject: [PATCH 04/10] test: order in first case too --- test/downstream/observables_autodiff.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/downstream/observables_autodiff.jl b/test/downstream/observables_autodiff.jl index 54bea8c99..537a81bdb 100644 --- a/test/downstream/observables_autodiff.jl +++ b/test/downstream/observables_autodiff.jl @@ -34,7 +34,7 @@ sol = solve(prob, Tsit5()) sum(sol[sys.w]) end du_ = [1.0, 1.0, 1.0, 0.0] - du = [du_ for _ in sol.u] + du = [du_ for _ in sol[[D(x), x, y, z]]] @test du == gs.u # Observable in a vector From 9c30360045f9b85204566de827e4a72ffee7c8b8 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 17 Jun 2025 13:49:05 +0530 Subject: [PATCH 05/10] build: bump Zygote compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f600c85ea..e8d0e607e 100644 --- a/Project.toml +++ b/Project.toml @@ -93,7 +93,7 @@ StaticArraysCore = "1.4" Statistics = "1.10" SymbolicIndexingInterface = "0.3.36" Tables = "1.11" -Zygote = "0.6.67, 0.7" +Zygote = "0.7.10" julia = "1.10" [extras] From 262019ca2811927c737686ea26ebf1dab9b5b3e9 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 17 Jun 2025 09:00:11 +0000 Subject: [PATCH 06/10] Update solution_interface.jl --- test/downstream/solution_interface.jl | 238 +++++++++++++------------- 1 file changed, 121 insertions(+), 117 deletions(-) diff --git a/test/downstream/solution_interface.jl b/test/downstream/solution_interface.jl index c31957db9..674c1d763 100644 --- a/test/downstream/solution_interface.jl +++ b/test/downstream/solution_interface.jl @@ -6,27 +6,29 @@ using Plots: Plots, plot ### Tests on non-layered model (everything should work). ### -@parameters a b c d -@variables s1(t) s2(t) - -eqs = [D(s1) ~ a * s1 / (1 + s1 + s2) - b * s1, - D(s2) ~ +c * s2 / (1 + s1 + s2) - d * s2] - -@mtkcompile population_model = System(eqs, t) - -# Tests on ODEProblem. -u0 = [s1 => 2.0, s2 => 1.0] -p = [a => 2.0, b => 1.0, c => 1.0, d => 1.0] -tspan = (0.0, 1000000.0) -oprob = ODEProblem(population_model, [u0; p], tspan) -sol = solve(oprob, Rodas4()) - -@test sol[s1] == sol[population_model.s1] == sol[:s1] -@test sol[s2] == sol[population_model.s2] == sol[:s2] -@test sol[s1][end] ≈ 1.0 -@test_throws Exception sol[a] -@test_throws Exception sol[population_model.a] -@test_throws Exception sol[:a] +@testset "Basic indexing" begin + @parameters a b c d + @variables s1(t) s2(t) + + eqs = [D(s1) ~ a * s1 / (1 + s1 + s2) - b * s1, + D(s2) ~ +c * s2 / (1 + s1 + s2) - d * s2] + + @mtkcompile population_model = System(eqs, t) + + # Tests on ODEProblem. + u0 = [s1 => 2.0, s2 => 1.0] + p = [a => 2.0, b => 1.0, c => 1.0, d => 1.0] + tspan = (0.0, 1000000.0) + oprob = ODEProblem(population_model, [u0; p], tspan) + sol = solve(oprob, Rodas4()) + + @test sol[s1] == sol[population_model.s1] == sol[:s1] + @test sol[s2] == sol[population_model.s2] == sol[:s2] + @test sol[s1][end] ≈ 1.0 + @test_throws Exception sol[a] + @test_throws Exception sol[population_model.a] + @test_throws Exception sol[:a] +end @testset "plot ODE solution" begin Plots.unicodeplots() @@ -51,102 +53,104 @@ sol = solve(oprob, Rodas4()) @test_nowarn plot(sol; plot_analytic = true) end -# Tests on SDEProblem -noiseeqs = [0.1 * s1, - 0.1 * s2] -@named noisy_population_model = SDESystem(population_model, noiseeqs) -noisy_population_model = complete(noisy_population_model) -sprob = SDEProblem(noisy_population_model, [u0; p], (0.0, 100.0)) -sol = solve(sprob, ImplicitEM()) - -@test sol[s1] == sol[noisy_population_model.s1] == sol[:s1] -@test sol[s2] == sol[noisy_population_model.s2] == sol[:s2] -@test_throws Exception sol[a] -@test_throws Exception sol[noisy_population_model.a] -@test_throws Exception sol[:a] -@test_nowarn sol(0.5, idxs = noisy_population_model.s1) -### Tests on layered model (some things should not work). ### - -@parameters σ ρ β -@variables x(t) y(t) z(t) - -eqs = [D(x) ~ σ * (y - x), - D(y) ~ x * (ρ - z) - y, - D(z) ~ x * y - β * z] - -@named lorenz1 = System(eqs, t) -@named lorenz2 = System(eqs, t) - -@parameters γ -@variables a(t) α(t) -connections = [0 ~ lorenz1.x + lorenz2.y + a * γ, - α ~ 2lorenz1.x + a * γ] -@mtkcompile sys = System(connections, t, [a, α], [γ], systems = [lorenz1, lorenz2]) - -u0 = [lorenz1.x => 1.0, - lorenz1.y => 0.0, - lorenz1.z => 0.0, - lorenz2.x => 0.0, - lorenz2.y => 1.0, - lorenz2.z => 0.0] - -p = [lorenz1.σ => 10.0, - lorenz1.ρ => 28.0, - lorenz1.β => 8 / 3, - lorenz2.σ => 10.0, - lorenz2.ρ => 28.0, - lorenz2.β => 8 / 3, - γ => 2.0] - -tspan = (0.0, 100.0) -prob = ODEProblem(sys, [u0; p], tspan) -sol = solve(prob, Rodas4()) - -@test_throws ArgumentError sol[x] -@test in(sol[lorenz1.x], [getindex.(sol.u, i) for i in 1:length(unknowns(sol.prob.f.sys))]) -@test_throws KeyError sol[:x] - -### Non-symbolic indexing tests -@test sol[:, 1] isa AbstractVector -@test sol[:, 1:2] isa AbstractDiffEqArray -@test sol[:, [1, 2]] isa AbstractDiffEqArray - -sol1 = sol(0.0:1.0:10.0) -@test sol1.u isa Vector -@test first(sol1.u) isa Vector -@test length(sol1.u) == 11 -@test length(sol1.t) == 11 - -sol2 = sol(0.1) -@test sol2 isa Vector -@test length(sol2) == length(unknowns(sys)) -@test first(sol2) isa Real - -sol3 = sol(0.0:1.0:10.0, idxs = [lorenz1.x, lorenz2.x]) - -sol7 = sol(0.0:1.0:10.0, idxs = [2, 1]) -@test sol7.u isa Vector -@test first(sol7.u) isa Vector -@test length(sol7.u) == 11 -@test length(sol7.t) == 11 -@test collect(sol7[t]) ≈ sol3.t -@test collect(sol7[t, 1:5]) ≈ sol3.t[1:5] - -sol8 = sol(0.1, idxs = [2, 1]) -@test sol8 isa Vector -@test length(sol8) == 2 -@test first(sol8) isa Real - -sol9 = sol(0.0:1.0:10.0, idxs = 2) -@test sol9.u isa Vector -@test first(sol9.u) isa Real -@test length(sol9.u) == 11 -@test length(sol9.t) == 11 -@test collect(sol9[t]) ≈ sol3.t -@test collect(sol9[t, 1:5]) ≈ sol3.t[1:5] - -sol10 = sol(0.1, idxs = 2) -@test sol10 isa Real +@testset "Symbolic Indexing" begin + # Tests on SDEProblem + noiseeqs = [0.1 * s1, + 0.1 * s2] + @named noisy_population_model = SDESystem(population_model, noiseeqs) + noisy_population_model = complete(noisy_population_model) + sprob = SDEProblem(noisy_population_model, [u0; p], (0.0, 100.0)) + sol = solve(sprob, ImplicitEM()) + + @test sol[s1] == sol[noisy_population_model.s1] == sol[:s1] + @test sol[s2] == sol[noisy_population_model.s2] == sol[:s2] + @test_throws Exception sol[a] + @test_throws Exception sol[noisy_population_model.a] + @test_throws Exception sol[:a] + @test_nowarn sol(0.5, idxs = noisy_population_model.s1) + ### Tests on layered model (some things should not work). ### + + @parameters σ ρ β + @variables x(t) y(t) z(t) + + eqs = [D(x) ~ σ * (y - x), + D(y) ~ x * (ρ - z) - y, + D(z) ~ x * y - β * z] + + @named lorenz1 = System(eqs, t) + @named lorenz2 = System(eqs, t) + + @parameters γ + @variables a(t) α(t) + connections = [0 ~ lorenz1.x + lorenz2.y + a * γ, + α ~ 2lorenz1.x + a * γ] + @mtkcompile sys = System(connections, t, [a, α], [γ], systems = [lorenz1, lorenz2]) + + u0 = [lorenz1.x => 1.0, + lorenz1.y => 0.0, + lorenz1.z => 0.0, + lorenz2.x => 0.0, + lorenz2.y => 1.0, + lorenz2.z => 0.0] + + p = [lorenz1.σ => 10.0, + lorenz1.ρ => 28.0, + lorenz1.β => 8 / 3, + lorenz2.σ => 10.0, + lorenz2.ρ => 28.0, + lorenz2.β => 8 / 3, + γ => 2.0] + + tspan = (0.0, 100.0) + prob = ODEProblem(sys, [u0; p], tspan) + sol = solve(prob, Rodas4()) + + @test_throws ArgumentError sol[x] + @test in(sol[lorenz1.x], [getindex.(sol.u, i) for i in 1:length(unknowns(sol.prob.f.sys))]) + @test_throws KeyError sol[:x] + + ### Non-symbolic indexing tests + @test sol[:, 1] isa AbstractVector + @test sol[:, 1:2] isa AbstractDiffEqArray + @test sol[:, [1, 2]] isa AbstractDiffEqArray + + sol1 = sol(0.0:1.0:10.0) + @test sol1.u isa Vector + @test first(sol1.u) isa Vector + @test length(sol1.u) == 11 + @test length(sol1.t) == 11 + + sol2 = sol(0.1) + @test sol2 isa Vector + @test length(sol2) == length(unknowns(sys)) + @test first(sol2) isa Real + + sol3 = sol(0.0:1.0:10.0, idxs = [lorenz1.x, lorenz2.x]) + + sol7 = sol(0.0:1.0:10.0, idxs = [2, 1]) + @test sol7.u isa Vector + @test first(sol7.u) isa Vector + @test length(sol7.u) == 11 + @test length(sol7.t) == 11 + @test collect(sol7[t]) ≈ sol3.t + @test collect(sol7[t, 1:5]) ≈ sol3.t[1:5] + + sol8 = sol(0.1, idxs = [2, 1]) + @test sol8 isa Vector + @test length(sol8) == 2 + @test first(sol8) isa Real + + sol9 = sol(0.0:1.0:10.0, idxs = 2) + @test sol9.u isa Vector + @test first(sol9.u) isa Real + @test length(sol9.u) == 11 + @test length(sol9.t) == 11 + @test collect(sol9[t]) ≈ sol3.t + @test collect(sol9[t, 1:5]) ≈ sol3.t[1:5] + + sol10 = sol(0.1, idxs = 2) + @test sol10 isa Real +end @testset "Plot idxs" begin @variables x(t) y(t) From c85becff358015aa8e937d37f6657aefe2bc53bd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 17 Jun 2025 09:06:06 +0000 Subject: [PATCH 07/10] Update solution_interface.jl --- test/downstream/solution_interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/downstream/solution_interface.jl b/test/downstream/solution_interface.jl index 674c1d763..01e6ffdb3 100644 --- a/test/downstream/solution_interface.jl +++ b/test/downstream/solution_interface.jl @@ -202,7 +202,7 @@ end (0.0, 1.0); build_initializeprob = false) dae_sol = solve(prob, DFBDF(); save_idxs = [x]) - @brownian a b + @brownians a b @mtkcompile sys = System([D(x) ~ x + p * y + x * a, D(y) ~ 2p + x^2 + y * b], t) xidx = variable_index(sys, x) prob = SDEProblem(sys, [x => 1.0, y => 2.0, p => 2.0], (0.0, 1.0)) From db6457052fd265a9d1e933c13d791a6ab6e2528a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 17 Jun 2025 09:44:49 +0000 Subject: [PATCH 08/10] Update solution_interface.jl --- test/downstream/solution_interface.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/downstream/solution_interface.jl b/test/downstream/solution_interface.jl index 01e6ffdb3..2821a79d0 100644 --- a/test/downstream/solution_interface.jl +++ b/test/downstream/solution_interface.jl @@ -28,9 +28,7 @@ using Plots: Plots, plot @test_throws Exception sol[a] @test_throws Exception sol[population_model.a] @test_throws Exception sol[:a] -end -@testset "plot ODE solution" begin Plots.unicodeplots() f = ODEFunction((u, p, t) -> -u, analytic = (u0, p, t) -> u0 * exp(-t)) @@ -51,9 +49,7 @@ end sol = solve(ode, Tsit5()) @test_nowarn plot(sol) @test_nowarn plot(sol; plot_analytic = true) -end -@testset "Symbolic Indexing" begin # Tests on SDEProblem noiseeqs = [0.1 * s1, 0.1 * s2] From 692939cb2a70e6feef902f99fb5c48e8ee2514d8 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 17 Jun 2025 10:20:34 +0000 Subject: [PATCH 09/10] Update modelingtoolkit_remake.jl --- test/downstream/modelingtoolkit_remake.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/downstream/modelingtoolkit_remake.jl b/test/downstream/modelingtoolkit_remake.jl index 5f20ba044..7e29bee56 100644 --- a/test/downstream/modelingtoolkit_remake.jl +++ b/test/downstream/modelingtoolkit_remake.jl @@ -18,7 +18,7 @@ eqs = [D(x) ~ σ * (y - x), D(y) ~ x * (ρ - z) - y, D(z) ~ x * y - β * z] -@named sys = System(eqs, t; parameter_dependencies = [q => 3β]) +@named sys = System(eqs, t; parameter_dependencies = [q ~ 3β]) sys = complete(sys) u0 = [x => 1.0, y => 0.0, From f8e347545e4cb6f993e5a6efadcba2f7582e1e70 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 17 Jun 2025 11:12:05 +0000 Subject: [PATCH 10/10] Update modelingtoolkit_remake.jl --- test/downstream/modelingtoolkit_remake.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/downstream/modelingtoolkit_remake.jl b/test/downstream/modelingtoolkit_remake.jl index 7e29bee56..7bb5c270d 100644 --- a/test/downstream/modelingtoolkit_remake.jl +++ b/test/downstream/modelingtoolkit_remake.jl @@ -70,7 +70,7 @@ push!(probs, OptimizationProblem(optsys, u0, p)) k = ShiftIndex(t) @mtkcompile discsys = System( [x ~ x(k - 1) * ρ + y(k - 2), y ~ y(k - 1) * σ - z(k - 2), z ~ z(k - 1) * β + x(k - 2)], - t; defaults = [x => 1.0, y => 1.0, z => 1.0]) + t; defaults = [x => 1.0, y => 1.0, z => 1.0, x(k-1) => 0.0, y(k-1) => 0.0, z(k-1) => 0.0]) # Roundabout method to avoid having to specify values for previous timestep discprob = DiscreteProblem(discsys, p, (0, 10)) for (var, v) in u0