Bounds update II (#23)

* Compat bounds update

* Format

Project.toml

@@ -20,11 +20,11 @@ Aqua = "0.7"
 CommonSolve = "0.2"
 ComponentArrays = "0.15"
 DocStringExtensions = "0.9"
-ModelingToolkit = "<8.73.0"
-Optimization = "<3.20.0"
+ModelingToolkit = "~8.72.0"
+Optimization = "~3.19.0"
 OrdinaryDiffEq = "6"
 Reexport = "1.2"
-SciMLBase = "<2.8.0"
+SciMLBase = "~2.7.0"
 SciMLSensitivity = "7"
 julia = "1.6"
 

docs/make.jl

@@ -12,10 +12,10 @@ makedocs(modules = [DynamicOED],
         "Home" => "index.md",
         "Examples" => [
             "Linear System" => "examples/1D.md",
-            "Lotka-Volterra" => "examples/lotka.md",
+            "Lotka-Volterra" => "examples/lotka.md"
         ],
         "Theory" => "theory.md",
-        "API" => "api.md",
+        "API" => "api.md"
     ])
 
 deploydocs(repo = "github.com/mathopt/DynamicOED.jl";

src/augment.jl

@@ -191,12 +191,12 @@ function build_augmented_system(sys::ModelingToolkit.AbstractODESystem,
     # TODO Maybe switch to variables here.
     @variables F(t)[1:np, 1:np]=zeros(T, (np, np)) [
         description = "Fisher Information Matrix",
-        fisher_state = true,
+        fisher_state = true
     ]
     @variables G(t)[1:nx, 1:np]=G_init [description = "Sensitivity State"]
     @variables Q(t)[1:n_obs, 1:np] [
         description = "Unweighted Fisher Information Derivative",
-        unweighted_information_gain = true,
+        unweighted_information_gain = true
     ]
 
     w = Vector{Num}(undef, n_obs)
@@ -233,9 +233,9 @@ function build_augmented_system(sys::ModelingToolkit.AbstractODESystem,
     # We do not need to do anything more than push this into the equations
     # Simplify will figure out the rest
     ODESystem([vec(dynamic_eqs);
-            vec(sens_eqs);
-            vec(new_obs);
-            vec(Q .~ hx * G)],
+               vec(sens_eqs);
+               vec(new_obs);
+               vec(Q .~ hx * G)],
         t,
         vcat(x, vec(G), vec(F[idx]), vec(Q)),
         vcat(union(p, p_tuneable), w),

src/discretize.jl

@@ -53,11 +53,12 @@ function Timegrid(tspan::Tuple, x::Num...; time_tolerance::Real = 1e-6)
 
     completegrid = [tspan_
                     for tspan_ in zip(_completegrid[1:(end - 1)], _completegrid[2:end])
-                        if abs(-(tspan_...)) >= time_tolerance]
+                    if abs(-(tspan_...)) >= time_tolerance]
 
     timegrids = map(_timegrids) do grid
-        [tspan_ for tspan_ in zip(grid[1:(end - 1)], grid[2:end])
-             if abs(-(tspan_...)) >= time_tolerance]
+        [tspan_
+         for tspan_ in zip(grid[1:(end - 1)], grid[2:end])
+         if abs(-(tspan_...)) >= time_tolerance]
     end
 
     indicators = zeros(Int, length(x), size(completegrid, 1))
@@ -296,7 +297,8 @@ function _sequential_solve(remaker::OEDRemake,
         u0::AbstractVector{T},
         p0::AbstractVector{T},
         idxs::TP;
-        kwargs...)::Tuple{AbstractArray{T, 2}, AbstractVector{T}} where {P, A, T, TP <: Tuple}
+        kwargs...)::Tuple{
+        AbstractArray{T, 2}, AbstractVector{T}} where {P, A, T, TP <: Tuple}
     __sequential_solve(remaker, prob, alg, parameters, u0, p0, idxs...)
 end
 

src/problem.jl

@@ -61,7 +61,8 @@ function ModelingToolkit.states(prob::OEDProblem)
 
     (;
         initial_conditions, controls = control_variables,
-        measurements = measurement_variables, regularization = regularization) |> sortkeys |> ComponentVector
+        measurements = measurement_variables, regularization = regularization) |>
+    sortkeys |> ComponentVector
 end
 
 function get_timegrids(prob::OEDProblem)
@@ -152,7 +153,7 @@ function Optimization.OptimizationProblem(prob::OEDProblem,
     # Declare the Optimization function
     opt_f = OptimizationFunction(objective, AD;
         syms = syms,
-        cons = cons,)
+        cons = cons)
 
     # Return the optimization problem
     OptimizationProblem(opt_f, u0, p, lb = lb, ub = ub, int = integrality,

test/criteria.jl

@@ -13,10 +13,10 @@ end
 
 @testset "Criteria" begin
     F = [5.68145 -0.148312 2.00382 5.1841 3.57814
-        -0.148312 2.17276 -1.18576 1.57651 -3.95544
-        2.00382 -1.18576 3.8703 -0.14832 2.69545
-        5.1841 1.57651 -0.14832 10.9617 0.0774284
-        3.57814 -3.95544 2.69545 0.0774284 24.8495]
+         -0.148312 2.17276 -1.18576 1.57651 -3.95544
+         2.00382 -1.18576 3.8703 -0.14832 2.69545
+         5.1841 1.57651 -0.14832 10.9617 0.0774284
+         3.57814 -3.95544 2.69545 0.0774284 24.8495]
 
     test_results = [(FisherACriterion(), -47.53571),
         (FisherDCriterion(), -2180.4107236837117),

test/discretize.jl

@@ -7,7 +7,7 @@ using Test
 @variables x(t)=2.0 [
     description = "State with uncertain initial condition",
     tunable = false,
-    bounds = (0.1, 5.0),
+    bounds = (0.1, 5.0)
 ] # The first state is uncertain
 @parameters p[1:1]=-2.0 [description = "Fixed parameter", tunable = true]
 @variables y₁(t) [description = "Observed", measurement_rate = 1.0]
@@ -68,7 +68,7 @@ end
         (0.9, 1.2),
         (1.2, 1.5),
         (1.5, 1.8),
-        (1.8, 2.0),
+        (1.8, 2.0)
     ]
     @test timegrid.timespans != timegrid.timegrids[1] != timegrid.timegrids[2]
     @test size(timegrid.timespans) == (8,) # We have
@@ -80,7 +80,7 @@ end
         (1.0, 1.2),
         (1.2, 1.5),
         (1.5, 1.8),
-        (1.8, 2.0),
+        (1.8, 2.0)
     ]
     # Right assignments
     @test DynamicOED.get_variable_idx(timegrid, Symbol("w₁"), 4) == 1
@@ -119,7 +119,7 @@ end
         (1.6, 1.7),
         (1.7, 1.8),
         (1.8, 1.9),
-        (1.9, 2.0),
+        (1.9, 2.0)
     ]
     # Right assignments
     @test DynamicOED.get_variable_idx(timegrid, Symbol("w₁"), 4) == 1

test/mtk_extensions.jl

@@ -6,7 +6,7 @@ using DynamicOED
 @variables x(t)=2.0 [
     description = "State with uncertain initial condition",
     tunable = false,
-    bounds = (0.1, 5.0),
+    bounds = (0.1, 5.0)
 ] # The first state is uncertain
 @parameters p[1:1]=-2.0 [description = "Fixed parameter", tunable = true]
 @variables y₁(t) [description = "Observed", measurement_rate = 1.0]

test/references/1D.jl

@@ -7,7 +7,7 @@ using Optimization, OptimizationMOI, Ipopt, Juniper
 
 const TEST_CRITERIA = [
     FisherACriterion(), FisherDCriterion(), FisherECriterion(),
-    ACriterion(), DCriterion(), ECriterion(),
+    ACriterion(), DCriterion(), ECriterion()
 ]
 
 ## Define the system
@@ -19,7 +19,7 @@ D = Differential(t)
 
 # Define the eqs
 @named simple_system = ODESystem([
-        D(x) ~ p[1] * x,
+        D(x) ~ p[1] * x
     ], tspan = (0.0, 1.0),
     observed = obs .~ [x])
 
@@ -40,7 +40,7 @@ oed = structural_simplify(oed)
         end
 
         constraints = [
-            0 ≲ 0.2 .- 0.5 * sum(Δts.w₁ .* optimization_variables.measurements.w₁),
+            0 ≲ 0.2 .- 0.5 * sum(Δts.w₁ .* optimization_variables.measurements.w₁)
         ]
 
         # Define an MTK Constraint system
@@ -81,7 +81,7 @@ end
         end
 
         constraints = [
-            0 ≲ 0.2 .- 0.5 * sum(Δts.w₁ .* optimization_variables.measurements.w₁),
+            0 ≲ 0.2 .- 0.5 * sum(Δts.w₁ .* optimization_variables.measurements.w₁)
         ]
 
         # Define an MTK Constraint system

test/references/LotkaVolterra.jl

@@ -8,29 +8,31 @@ using Optimization, OptimizationMOI, Ipopt, Juniper
 @testset "Relaxed" begin
     @variables t [description = "Time"]
     @variables x(t)=0.5 [description = "Biomass Prey"] y(t)=0.7 [
-        description = "Biomass Predator",
+        description = "Biomass Predator"
     ]
     @parameters p[1:4]=[1.0; 1.0; 0.4; 0.6] [
         description = "Fixed parameters",
-        tunable = false,
+        tunable = false
     ]
     @parameters c[1:2]=[1.0; 1.0] [description = "Uncertain parameters", tunable = true]
     @parameters u=0.0 [
         description = "Binary control variable, measured 10 times over the provided time span, relaxed",
         measurement_rate = 0.3,
         input = true,
-        bounds = (0.0, 1.0),
+        bounds = (0.0, 1.0)
     ]
     @variables obs(t)[1:2] [
         description = "Observed variable measured 10 times over the provided time span",
-        measurement_rate = 0.1,
+        measurement_rate = 0.1
     ]
     obs = collect(obs)
     D = Differential(t)
 
     # Define the ODE System
-    @named lotka_volterra = ODESystem([D(x) ~ p[1] * x - c[1] * x * y - p[3] * x * u;
-            D(y) ~ -p[2] * y + c[2] * x * y - p[4] * y * u], tspan = (0.0, 3.0),
+    @named lotka_volterra = ODESystem(
+        [D(x) ~ p[1] * x - c[1] * x * y - p[3] * x * u;
+         D(y) ~ -p[2] * y + c[2] * x * y - p[4] * y * u],
+        tspan = (0.0, 3.0),
         observed = obs .~ [x; y])
 
     @named oed_lotka = OEDSystem(lotka_volterra)
@@ -51,7 +53,7 @@ using Optimization, OptimizationMOI, Ipopt, Juniper
     constraints = [
         0 ≲ 0.2 .- 0.5 * sum(Δts.w₁ .* optimization_variables.measurements.w₁),
         0 ≲ 0.2 .- 0.5 * sum(Δts.w₁ .* optimization_variables.measurements.w₂),
-        1.0 ≳ sum(Δts.u .* optimization_variables.controls.u),
+        1.0 ≳ sum(Δts.u .* optimization_variables.controls.u)
     ]
 
     # Define an MTK Constraint system
@@ -77,7 +79,7 @@ using Optimization, OptimizationMOI, Ipopt, Juniper
             -8.770673425035107e-9,
             -8.822579849602493e-9,
             -8.756712678420366e-9,
-            -7.15514843922874e-9,
+            -7.15514843922874e-9
         ], atol = 1e-2, rtol = 1e-5)
 
     @test isapprox(u_opt.measurements.w₁,
@@ -111,7 +113,7 @@ using Optimization, OptimizationMOI, Ipopt, Juniper
             0.9999999700398273,
             0.9999999928632054,
             0.9999999994598789,
-            1.000000002696486,
+            1.000000002696486
         ], atol = 1e-2, rtol = 1e-5)
 
     @test isapprox(u_opt.measurements.w₂,
@@ -145,7 +147,7 @@ using Optimization, OptimizationMOI, Ipopt, Juniper
             0.9999999865101933,
             0.9999999999387468,
             1.0000000041944848,
-            1.0000000061707939,
+            1.0000000061707939
         ], atol = 1e-2, rtol = 1e-5)
 
     @info res.objective
@@ -155,29 +157,31 @@ end
 @testset "Integer" begin
     @variables t [description = "Time"]
     @variables x(t)=0.5 [description = "Biomass Prey"] y(t)=0.7 [
-        description = "Biomass Predator",
+        description = "Biomass Predator"
     ]
     @parameters p[1:4]=[1.0; 1.0; 0.4; 0.6] [
         description = "Fixed parameters",
-        tunable = false,
+        tunable = false
     ]
     @parameters c[1:2]=[1.0; 1.0] [description = "Uncertain parameters", tunable = true]
     @parameters u::Int=0 [
         description = "Binary control variable, measured 10 times over the provided time span, discrete",
         measurement_rate = 0.3,
         input = true,
-        bounds = (0.0, 1.0),
+        bounds = (0.0, 1.0)
     ]
     @variables obs(t)[1:2] [
         description = "Observed variable measured 10 times over the provided time span",
-        measurement_rate = 0.1,
+        measurement_rate = 0.1
     ]
     obs = collect(obs)
     D = Differential(t)
 
     # Define the ODE System
-    @named lotka_volterra = ODESystem([D(x) ~ p[1] * x - c[1] * x * y - p[3] * x * u;
-            D(y) ~ -p[2] * y + c[2] * x * y - p[4] * y * u], tspan = (0.0, 3.0),
+    @named lotka_volterra = ODESystem(
+        [D(x) ~ p[1] * x - c[1] * x * y - p[3] * x * u;
+         D(y) ~ -p[2] * y + c[2] * x * y - p[4] * y * u],
+        tspan = (0.0, 3.0),
         observed = obs .~ [x; y])
 
     @named oed_lotka = OEDSystem(lotka_volterra)
@@ -201,7 +205,7 @@ end
         0 ≲ 0.2 .- 0.5 * sum(Δts.w₁ .* optimization_variables.measurements.w₁),
         0 ≲ 0.2 .- 0.5 * sum(Δts.w₁ .* optimization_variables.measurements.w₂),
         2.0 ≳ sum(optimization_variables.controls.u),
-        sum(optimization_variables.controls.u) ≳ 1.0,
+        sum(optimization_variables.controls.u) ≳ 1.0
     ]
 
     # Define an MTK Constraint system
@@ -226,7 +230,7 @@ end
             -8.072531050877907e-9,
             -8.10418072743497e-9,
             -7.828539249999535e-9,
-            2.2919106922838742e-8,
+            2.2919106922838742e-8
         ], atol = 1e-2, rtol = 1e-5)
     @test isapprox(u_opt.measurements.w₁,
         [
@@ -259,7 +263,7 @@ end
             0.9999999785713451,
             0.9999999930680564,
             1.0000000001481675,
-            1.0000000033846281,
+            1.0000000033846281
         ], atol = 1e-2, rtol = 1e-5)
 
     @test isapprox(u_opt.measurements.w₂,
@@ -293,7 +297,7 @@ end
             0.999999969264283,
             1.000000000489115,
             1.000000004387359,
-            1.0000000064340413,
+            1.0000000064340413
         ], atol = 1e-2, rtol = 1e-5)
     @test isapprox(res.objective, 0.34, atol = 1e-2)
 end
@@ -303,20 +307,22 @@ end
     @variables x(t)=0.49 [
         description = "Biomass Prey",
         tunable = true,
-        bounds = (0.1, 1.0),
+        bounds = (0.1, 1.0)
     ] y(t)=0.7 [description = "Biomass Predator", tunable = false]
     @parameters p[1:3]=[1.0; 1.0; 1.0] [description = "Fixed parameters", tunable = false]
     @parameters c[1:1]=[1.0;] [description = "Uncertain parameters", tunable = true]
     @variables obs(t)[1:1] [
         description = "Observed variable measured 10 times over the provided time span",
-        measurement_rate = 10,
+        measurement_rate = 10
     ]
     obs = collect(obs)
     D = Differential(t)
 
     # Define the ODE System
-    @named lotka_volterra = ODESystem([D(x) ~ p[1] * x - c[1] * x * y;
-            D(y) ~ -p[2] * y + p[3] * x * y], tspan = (0.0, 5.0),
+    @named lotka_volterra = ODESystem(
+        [D(x) ~ p[1] * x - c[1] * x * y;
+         D(y) ~ -p[2] * y + p[3] * x * y],
+        tspan = (0.0, 5.0),
         observed = obs .~ [y + x])
 
     @named oed_lotka = OEDSystem(lotka_volterra)
@@ -334,7 +340,7 @@ end
     end
 
     constraints = [
-        0 ≲ 1 .- sum(Δts.w₁ .* optimization_variables.measurements.w₁),
+        0 ≲ 1 .- sum(Δts.w₁ .* optimization_variables.measurements.w₁)
     ]
 
     # Define an MTK Constraint system
@@ -360,7 +366,7 @@ end
             4.771191106750985e-6,
             7.588080619413368e-5,
             0.9999093339443416,
-            0.9999941524754293,
+            0.9999941524754293
         ],
         atol = 1e-2,
         rtol = 1e-5)