Merge branch 'master' into myb/nonsq

Author: YingboMa

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "9.5.0"
+version = "9.6.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

docs/src/tutorials/SampledData.md

@@ -39,7 +39,7 @@ k = ShiftIndex(clock)
 
 eqs = [
     x(k) ~ 0.5x(k - 1) + u(k - 1),
-    y(k) ~ x(k - 1)
+    y ~ x
 ]
 ```
 
@@ -49,13 +49,13 @@ A few things to note in this basic example:
   - `x` and `u` are automatically inferred to be discrete-time variables, since they appear in an equation with a discrete-time [`ShiftIndex`](@ref) `k`.
   - `y` is also automatically inferred to be a discrete-time-time variable, since it appears in an equation with another discrete-time variable `x`. `x,u,y` all belong to the same discrete-time partition, i.e., they are all updated at the same *instantaneous point in time* at which the clock ticks.
   - The equation `y ~ x` does not use any shift index, this is equivalent to `y(k) ~ x(k)`, i.e., discrete-time variables without shift index are assumed to refer to the variable at the current time step.
-  - The equation `x(k) ~ 0.5x(k-1) + u(k-1)` indicates how `x` is updated, i.e., what the value of `x` will be at the *current* time step in terms of the *past* value. The output `y`, is given by the value of `x` at the *past* time step, i.e., `y(k) ~ x(k-1)`. If this logic was implemented in an imperative programming style, the logic would thus be
+  - The equation `x(k) ~ 0.5x(k-1) + u(k-1)` indicates how `x` is updated, i.e., what the value of `x` will be at the *current* time step in terms of the *past* value. The output `y`, is given by the value of `x` at the *current* time step, i.e., `y(k) ~ x(k)`. If this logic was implemented in an imperative programming style, the logic would thus be
 
 ```julia
 function discrete_step(x, u)
-    y = x # y is assigned the current value of x, y(k) = x(k-1)
     x = 0.5x + u # x is updated to a new value, i.e., x(k) is computed
-    return x, y # The state x now refers to x at the current time step, x(k), while y refers to x at the past time step, y(k) = x(k-1)
+    y = x # y is assigned the current value of x, y(k) = x(k)
+    return x, y # The state x now refers to x at the current time step, x(k), and y equals x, y(k) = x(k)
 end
 ```
 
@@ -68,22 +68,21 @@ eqs = [
 ]
 ```
 
-but the use of positive time shifts is not yet supported.
-Instead, we have *shifted all indices* by `-1`, resulting in exactly the same difference equations. However, the next system is *not equivalent* to the previous one:
+but the use of positive time shifts is not yet supported. Instead, we *shifted all indices* by `-1` above, resulting in exactly the same difference equations. However, the next system is *not equivalent* to the previous one:
 
 ```@example clocks
 eqs = [
-    x(k) ~ 0.5x(k - 1) + u(k - 1),
+    x(k) ~ 0.5x(k - 1) + u(k),
     y ~ x
 ]
 ```
 
-In this last example, `y` refers to the updated `x(k)`, i.e., this system is equivalent to
+In this last example, `u(k)` refers to the input at the new time point `k`., this system is equivalent to
 
 ```
 eqs = [
-    x(k+1) ~ 0.5x(k) + u(k),
-    y(k+1) ~ x(k+1)
+    x(k+1) ~ 0.5x(k) + u(k+1),
+    y(k) ~ x(k)
 ]
 ```
 

docs/src/tutorials/parameter_identifiability.md

@@ -69,7 +69,7 @@ After that, we are ready to check the system for local identifiability:
 ```julia
 # query local identifiability
 # we pass the ode-system
-local_id_all = assess_local_identifiability(de, p = 0.99)
+local_id_all = assess_local_identifiability(de, prob_threshold = 0.99)
 ```
 
 We can see that all unknowns (except $x_7$) and all parameters are locally identifiable with probability 0.99.
@@ -78,7 +78,7 @@ Let's try to check specific parameters and their combinations
 
 ```julia
 to_check = [de.k5, de.k7, de.k10 / de.k9, de.k5 + de.k6]
-local_id_some = assess_local_identifiability(de, funcs_to_check = to_check, p = 0.99)
+local_id_some = assess_local_identifiability(de, funcs_to_check = to_check, prob_threshold = 0.99)
 ```
 
 Notice that in this case, everything (except the unknown variable $x_7$) is locally identifiable, including combinations such as $k_{10}/k_9, k_5+k_6$
@@ -183,7 +183,7 @@ end
 # check only 2 parameters
 to_check = [ode.b, ode.c]
 
-global_id = assess_identifiability(ode, funcs_to_check = to_check, p = 0.9)
+global_id = assess_identifiability(ode, funcs_to_check = to_check, prob_threshold = 0.9)
 ```
 
 Both parameters `b, c` are globally identifiable with probability `0.9` in this case.

src/systems/diffeqs/abstractodesystem.jl

@@ -870,12 +870,15 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
         if eltype(u0map) <: Number
             u0map = unknowns(sys) .=> u0map
         end
+        if isempty(u0map)
+            u0map = Dict()
+        end
         initializeprob = ModelingToolkit.InitializationProblem(
             sys, t, u0map, parammap; guesses, warn_initialize_determined)
         initializeprobmap = getu(initializeprob, unknowns(sys))
 
-        zerovars = setdiff(unknowns(sys), keys(defaults(sys))) .=> 0.0
-        trueinit = identity.([zerovars; u0map])
+        zerovars = Dict(setdiff(unknowns(sys), keys(defaults(sys))) .=> 0.0)
+        trueinit = collect(merge(zerovars, eltype(u0map) <: Pair ? todict(u0map) : u0map))
         u0map isa StaticArraysCore.StaticArray &&
             (trueinit = SVector{length(trueinit)}(trueinit))
     else
@@ -913,7 +916,6 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
         du0 = nothing
         ddvs = nothing
     end
-
     check_eqs_u0(eqs, dvs, u0; kwargs...)
 
     f = constructor(sys, dvs, ps, u0; ddvs = ddvs, tgrad = tgrad, jac = jac,

src/systems/systemstructure.jl

@@ -630,7 +630,7 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     if has_io
         ModelingToolkit.markio!(state, orig_inputs, io...)
     end
-    if io !== nothing || any(isinput, state.fullvars)
+    if io !== nothing
         state, input_idxs = ModelingToolkit.inputs_to_parameters!(state, io)
     else
         input_idxs = 0:-1 # Empty range

test/input_output_handling.jl

@@ -50,7 +50,7 @@ end
 @test !is_bound(sys31, sys1.v[2])
 
 # simplification turns input variables into parameters
-ssys = structural_simplify(sys)
+ssys, _ = structural_simplify(sys, ([u], []))
 @test ModelingToolkit.isparameter(unbound_inputs(ssys)[])
 @test !is_bound(ssys, u)
 @test u ∈ Set(unbound_inputs(ssys))
@@ -236,7 +236,7 @@ i = findfirst(isequal(u[1]), out)
 @variables x(t) u(t) [input = true]
 eqs = [D(x) ~ u]
 @named sys = ODESystem(eqs, t)
-@test_nowarn structural_simplify(sys)
+@test_nowarn structural_simplify(sys, ([u], []))
 
 #=
 ## Disturbance input handling

test/reduction.jl

@@ -234,7 +234,7 @@ eqs = [D(x) ~ σ * (y - x)
        u ~ z + a]
 
 lorenz1 = ODESystem(eqs, t, name = :lorenz1)
-lorenz1_reduced = structural_simplify(lorenz1)
+lorenz1_reduced, _ = structural_simplify(lorenz1, ([z], []))
 @test z in Set(parameters(lorenz1_reduced))
 
 # #2064

test/serialization.jl

@@ -65,4 +65,4 @@ probexpr = ODEProblemExpr{true}(ss, [capacitor.v => 0.0], (0, 0.1); observedfun_
 prob_obs = eval(probexpr)
 sol_obs = solve(prob_obs, ImplicitEuler())
 @show all_obs
-@test_broken sol_obs[all_obs] == sol[all_obs]
+@test sol_obs[all_obs] == sol[all_obs]

test/structural_transformation/tearing.jl

@@ -164,9 +164,10 @@ prob.f(du, u, pr, tt)
 @test du≈[u[2], u[1] + sin(u[2]) - pr * tt] atol=1e-5
 
 # test the initial guess is respected
-@named sys = ODESystem(eqs, t, defaults = Dict(z => Inf))
+@named sys = ODESystem(eqs, t, defaults = Dict(z => NaN))
 infprob = ODEProblem(structural_simplify(sys), [x => 1.0], (0, 1.0), [p => 0.2])
-@test_throws Any infprob.f(du, infprob.u0, pr, tt)
+infprob.f(du, infprob.u0, pr, tt)
+@test any(isnan, du)
 
 sol1 = solve(prob, RosShamp4(), reltol = 8e-7)
 sol2 = solve(ODEProblem{false}((u, p, t) -> [-asin(u[1] - pr * t)],