Try updating compats

.github/workflows/CI.yml

@@ -20,7 +20,6 @@ jobs:
           - Datafit
         version:
           - '1'
-          - '1.6'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1

.github/workflows/Documentation.yml

@@ -15,7 +15,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '1.9'
+          version: '1'
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

Project.toml

@@ -25,19 +25,19 @@ DifferentialEquations = "7"
 Distributions = "0.25"
 GlobalSensitivity = "2"
 LinearAlgebra = "1"
-ModelingToolkit = "8"
+ModelingToolkit = "9.12"
 NLopt = "0.6, 1"
 Optimization = "3"
-OptimizationBBO = "0.1, 0.2"
-OptimizationNLopt = "0.1, 0.2"
+OptimizationBBO = "0.2"
+OptimizationNLopt = "0.2"
 Plots = "1"
 Reexport = "1"
 SafeTestsets = "0.1"
-SciMLBase = "1.93.1, 2"
+SciMLBase = "2.32.1"
 SciMLExpectations = "2"
 Test = "1"
-Turing = "0.22, 0.23, 0.24"
-julia = "1.6"
+Turing = "0.30"
+julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

docs/Project.toml

@@ -22,10 +22,10 @@ DifferentialEquations = "7"
 Documenter = "1"
 EasyModelAnalysis = "0.1"
 Ipopt = "1"
-ModelingToolkit = "8.43"
+ModelingToolkit = "9"
 NLopt = "0.6, 1"
 Optimization = "3"
-OptimizationMOI = "0.1, 0.2, 0.3, 0.4"
+OptimizationMOI = "0.4"
 Plots = "1"
 SBMLToolkit = "0.1"
 UnPack = "1"

src/intervention.jl

@@ -112,7 +112,8 @@ function optimal_parameter_intervention_for_threshold(prob, obs, threshold,
     symbolic_cost = Symbolics.unwrap(symbolic_cost)
     #ps = collect(ModelingToolkit.vars(symbolic_cost))
     _cost = Symbolics.build_function(symbolic_cost, ps, expression = Val{false})
-    _cost(prob.p) # just throw when something is wrong during the setup.
+
+    _cost(prob.p[1]) # just throw when something is wrong during the setup.
 
     cost = let _cost = _cost
         (x, grad) -> _cost(x)
@@ -212,7 +213,8 @@ function optimal_parameter_intervention_for_reach(prob, obs, reach,
     symbolic_cost = Symbolics.unwrap(symbolic_cost)
     #ps = collect(ModelingToolkit.vars(symbolic_cost))
     _cost = Symbolics.build_function(symbolic_cost, ps, expression = Val{false})
-    _cost(prob.p) # just throw when something is wrong during the setup.
+
+    _cost(prob.p[1]) # just throw when something is wrong during the setup.
 
     cost = let _cost = _cost, cost_sol = cost_sol
         (x, grad) -> begin
@@ -229,7 +231,7 @@ function optimal_parameter_intervention_for_reach(prob, obs, reach,
                 sol_intervention = solve(prob_intervention; kw...)
                 sol_cost = cost_sol(sol_intervention)
             end
-            p_cost + sol_cost
+            only(p_cost .+ sol_cost)
         end
     end
 

src/threshold.jl

@@ -60,8 +60,8 @@ function prob_violating_threshold(prob, p, thresholds)
         end
         return 0.0
     end
-    exprob = ExpectationProblem(sm, g, h, gd; nout = 1)
-    exsol = solve(exprob, Koopman(), batch = 0, quadalg = HCubatureJL())
+    exprob = ExpectationProblem(sm, g, h, gd)
+    exsol = solve(exprob, Koopman(), quadalg = HCubatureJL())
     exsol.u
 end
 

test/basics.jl

@@ -8,8 +8,7 @@ eqs = [D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys = ODESystem(eqs)
-sys = structural_simplify(sys)
+@mtkbuild sys = ODESystem(eqs, t)
 
 u0 = [D(x) => 2.0,
     x => 1.0,

test/datafit.jl

@@ -8,8 +8,7 @@ eqs = [D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys = ODESystem(eqs)
-sys = structural_simplify(sys)
+@mtkbuild sys = ODESystem(eqs, t)
 
 u0 = [D(x) => 2.0,
     x => 1.0,
@@ -68,8 +67,7 @@ eqs_obs = [D(D(x)) ~ σ * (y - x),
     D(z) ~ x * y - β * z,
     x_2 ~ 2 * x]
 
-@named sys_obs = ODESystem(eqs_obs)
-sys_obs = structural_simplify(sys_obs)
+@mtkbuild sys_obs = ODESystem(eqs_obs, t)
 
 u0_obs = [D(x) => 2.0,
     x => 1.0,
@@ -88,7 +86,7 @@ tsave = collect(10.0:10.0:100.0)
 sol_data = solve(prob, saveat = tsave)
 data = [x => sol_data[x], z => sol_data[z]]
 p_prior = [σ => Normal(26.8, 0.1), β => Normal(2.7, 0.1)]
-p_posterior = @time bayesian_datafit(prob, p_prior, tsave, data)
+p_posterior = @time bayesian_datafit(prob, p_prior, tsave, data, niter = 3000)
 @test var.(getfield.(p_prior, :second)) >= var.(getfield.(p_posterior, :second))
 
 tsave1 = collect(10.0:10.0:100.0)
@@ -97,5 +95,5 @@ tsave2 = collect(10.0:13.5:100.0)
 sol_data2 = solve(prob, saveat = tsave2)
 data_with_t = [x => (tsave1, sol_data1[x]), z => (tsave2, sol_data2[z])]
 
-p_posterior = @time bayesian_datafit(prob, p_prior, data_with_t)
+p_posterior = @time bayesian_datafit(prob, p_prior, data_with_t, niter = 3000)
 @test var.(getfield.(p_prior, :second)) >= var.(getfield.(p_posterior, :second))

test/ensemble.jl

@@ -8,7 +8,7 @@ eqs = [∂(S) ~ -β * c * I / N * S,
     ∂(I) ~ β * c * I / N * S - γ * I,
     ∂(R) ~ γ * I];
 
-@named sys = ODESystem(eqs);
+@mtkbuild sys = ODESystem(eqs, t);
 tspan = (0, 30)
 prob = ODEProblem(sys, [], tspan);
 
@@ -22,7 +22,7 @@ eqs = [∂(S) ~ -β * c * I / N * S,
     ∂(H) ~ h * I - r * H - d * H,
     ∂(D) ~ ρ * I + d * H];
 
-@named sys2 = ODESystem(eqs);
+@mtkbuild sys2 = ODESystem(eqs, t);
 
 prob2 = ODEProblem(sys2, [], tspan);
 
@@ -47,8 +47,7 @@ eqs = [∂(S) ~ -β * c * I_total / N * S - v * Sv,
     I_total ~ I + Iv
 ];
 
-@named sys3 = ODESystem(eqs)
-sys3 = structural_simplify(sys3)
+@mtkbuild sys3 = ODESystem(eqs, t)
 prob3 = ODEProblem(sys3, [], tspan);
 enprob = EnsembleProblem([prob, prob2, prob3])
 

test/examples.jl

@@ -10,7 +10,7 @@ eqs = [Dₜ(S) ~ μₙ * S - μₘ * S - θ * α * S * Iₛ - (1 - θ) * α * S
        Dₜ(Rₐ) ~ βₐ * Iₐ - ρ * Rₐ
        Dₜ(Rₛ) ~ (1 - ωₛ) * βₛ * Iₛ - ρ * Rₛ
        Dₜ(D) ~ ωₛ * βₛ * Iₛ]
-@named asir = ODESystem(eqs)
+@mtkbuild asir = ODESystem(eqs, t)
 prob = ODEProblem(asir, [], (0, 110.0))
 sol = solve(prob)
 # plot(sol)
@@ -33,8 +33,7 @@ eqs = [T ~ S + E + I + R + H + D
        Dₜ(R) ~ γ₁ * I + γ₂ * H
        Dₜ(H) ~ δ * I - (μ + γ₂) * H
        Dₜ(D) ~ μ * H];
-@named seirhd = ODESystem(eqs)
-seirhd = structural_simplify(seirhd)
+@mtkbuild seirhd = ODESystem(eqs, t)
 prob = ODEProblem(seirhd, [], (0, 110.0))
 sol = solve(prob)
 # plot(sol)

test/runtests.jl

@@ -1,15 +1,17 @@
-using SafeTestsets
+using SafeTestsets, Test
 
 const GROUP = get(ENV, "GROUP", "All")
 
 @time begin
     if GROUP == "All" || GROUP == "Core"
-        @time @safetestset "Quality Assurance" include("qa.jl")
-        @time @safetestset "Basic Tests" include("basics.jl")
-        @time @safetestset "Sensitivity Tests" include("sensitivity.jl")
-        VERSION >= v"1.9" && @time @safetestset "Ensemble Tests" include("ensemble.jl")
-        @time @safetestset "Threshold Tests" include("threshold.jl")
-        @time @safetestset "Example Tests" include("examples.jl")
+        @testset "Core" begin
+            @time @safetestset "Quality Assurance" include("qa.jl")
+            @time @safetestset "Basic Tests" include("basics.jl")
+            @time @safetestset "Sensitivity Tests" include("sensitivity.jl")
+            @time @safetestset "Ensemble Tests" include("ensemble.jl")
+            @time @safetestset "Threshold Tests" include("threshold.jl")
+            @time @safetestset "Example Tests" include("examples.jl")
+        end
     end
 
     if GROUP == "All" || GROUP == "Datafit"

test/sensitivity.jl

@@ -8,8 +8,7 @@ eqs = [D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys = ODESystem(eqs)
-sys = structural_simplify(sys)
+@mtkbuild sys = ODESystem(eqs, t)
 
 u0 = [D(x) => 2.0,
     x => 1.0,

test/threshold.jl

@@ -8,8 +8,7 @@ eqs = [D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys = ODESystem(eqs)
-sys = structural_simplify(sys)
+@mtkbuild sys = ODESystem(eqs, t)
 
 u0 = [D(x) => 2.0,
     x => 1.0,
@@ -28,7 +27,7 @@ sol = solve(prob)
 @variables t x(t)
 D = Differential(t)
 eqs = [D(x) ~ x]
-@named sys = ODESystem(eqs)
+@mtkbuild sys = ODESystem(eqs, t)
 prob = ODEProblem(sys, [x => 0.01], (0.0, Inf))
 sol = stop_at_threshold(prob, x^2, 0.1)
 @test sol.u[end][1]^2≈0.1 atol=1e-5
@@ -38,7 +37,7 @@ sol = stop_at_threshold(prob, x^2, 0.1)
 @parameters p
 D = Differential(t)
 eqs = [D(x) ~ p * x]
-@named sys = ODESystem(eqs)
+@mtkbuild sys = ODESystem(eqs, t)
 prob = ODEProblem(sys, [x => 0.01], (0.0, 50), [p => 1.0])
 opt_tspan, (s1, s2, s3), ret = optimal_threshold_intervention(prob, [p => -1.0], x, 3, 50);
 @test -(-(opt_tspan...)) < 25
@@ -85,7 +84,7 @@ opt_ps, (s1, s2, s3), ret = optimal_parameter_intervention_for_reach(
 D = Differential(t)
 eqs = [D(x) ~ p1 * abs(x) + p2 * y
        D(y) ~ p1 * abs(x) + p2 * y]
-@named sys = ODESystem(eqs)
+@mtkbuild sys = ODESystem(eqs, t)
 prob = ODEProblem(sys, [x => 0.01, y => 1], (0.0, 50), [p1 => 0.5, p2 => 0.2])
 opt_ps, s2, ret = optimal_parameter_threshold(prob, x, 2, p1 - p2, [p1, p2],
     [-2.0, -2.0], [2.0, 2],
@@ -102,7 +101,7 @@ eqs = [D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys = ODESystem(eqs)
+@mtkbuild sys = ODESystem(eqs, t)
 sys = structural_simplify(sys)
 
 u0 = [D(x) => 2.0,