Merge branch 'SciML:master' into iss3707

Author: fchen121

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 = "10.6.0"
+version = "10.7.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -106,7 +106,7 @@ DiffEqBase = "6.170.1"
 DiffEqCallbacks = "2.16, 3, 4"
 DiffEqNoiseProcess = "5"
 DiffRules = "0.1, 1.0"
-DifferentiationInterface = "0.6.47"
+DifferentiationInterface = "0.6.47, 0.7"
 Distributed = "1"
 Distributions = "0.23, 0.24, 0.25"
 DocStringExtensions = "0.7, 0.8, 0.9"

docs/src/tutorials/fmi.md

@@ -85,7 +85,7 @@ We can interpolate the solution object to obtain values at arbitrary time points
 just like a normal solution.
 
 ```@repl fmi
-sol(0.0:0.1:1.0; idxs = sys.mass_a)
+sol(0.0:0.1:1.0; idxs = sys.mass__a)
 ```
 
 FMUs following version 3 of the specification can be simulated with almost the same process. This time,

src/linearization.jl

@@ -174,13 +174,13 @@ struct PreparedJacobian{iip, P, F, B, A}
 end
 
 function PreparedJacobian{true}(f, buf, autodiff, args...)
-    prep = DI.prepare_jacobian(f, buf, autodiff, args...)
+    prep = DI.prepare_jacobian(f, buf, autodiff, args...; strict = Val(false))
     return PreparedJacobian{true, typeof(prep), typeof(f), typeof(buf), typeof(autodiff)}(
         prep, f, buf, autodiff)
 end
 
 function PreparedJacobian{false}(f, autodiff, args...)
-    prep = DI.prepare_jacobian(f, autodiff, args...)
+    prep = DI.prepare_jacobian(f, autodiff, args...; strict = Val(false))
     return PreparedJacobian{true, typeof(prep), typeof(f), Nothing, typeof(autodiff)}(
         prep, f, nothing)
 end

src/systems/abstractsystem.jl

@@ -820,6 +820,9 @@ function refreshed_metadata(meta::Base.ImmutableDict)
         end
         newmeta = Base.ImmutableDict(newmeta, k => v)
     end
+    if !haskey(newmeta, MutableCacheKey)
+        newmeta = Base.ImmutableDict(newmeta, MutableCacheKey => MutableCacheT())
+    end
     return newmeta
 end
 
@@ -2683,22 +2686,27 @@ function Symbolics.substitute(sys::AbstractSystem, rules::Union{Vector{<:Pair},
     elseif sys isa System
         rules = todict(map(r -> Symbolics.unwrap(r[1]) => Symbolics.unwrap(r[2]),
             collect(rules)))
-        eqs = fast_substitute(get_eqs(sys), rules)
-        pdeps = fast_substitute(get_parameter_dependencies(sys), rules)
-        defs = Dict(fast_substitute(k, rules) => fast_substitute(v, rules)
+        newsys = @set sys.eqs = fast_substitute(get_eqs(sys), rules)
+        @set! newsys.unknowns = map(get_unknowns(sys)) do var
+            get(rules, var, var)
+        end
+        @set! newsys.ps = map(get_ps(sys)) do var
+            get(rules, var, var)
+        end
+        @set! newsys.parameter_dependencies = fast_substitute(
+            get_parameter_dependencies(sys), rules)
+        @set! newsys.defaults = Dict(fast_substitute(k, rules) => fast_substitute(v, rules)
         for (k, v) in get_defaults(sys))
-        guess = Dict(fast_substitute(k, rules) => fast_substitute(v, rules)
+        @set! newsys.guesses = Dict(fast_substitute(k, rules) => fast_substitute(v, rules)
         for (k, v) in get_guesses(sys))
-        noise_eqs = fast_substitute(get_noise_eqs(sys), rules)
-        costs = fast_substitute(get_costs(sys), rules)
-        observed = fast_substitute(get_observed(sys), rules)
-        initialization_eqs = fast_substitute(get_initialization_eqs(sys), rules)
-        cstrs = fast_substitute(get_constraints(sys), rules)
-        subsys = map(s -> substitute(s, rules), get_systems(sys))
-        newsys = System(eqs, get_iv(sys); name = nameof(sys), defaults = defs,
-            guesses = guess, systems = subsys, noise_eqs,
-            observed, initialization_eqs, constraints = cstrs)
-        @set! newsys.parameter_dependencies = pdeps
+        @set! newsys.noise_eqs = fast_substitute(get_noise_eqs(sys), rules)
+        @set! newsys.costs = Vector{Union{Real, BasicSymbolic}}(fast_substitute(
+            get_costs(sys), rules))
+        @set! newsys.observed = fast_substitute(get_observed(sys), rules)
+        @set! newsys.initialization_eqs = fast_substitute(
+            get_initialization_eqs(sys), rules)
+        @set! newsys.constraints = fast_substitute(get_constraints(sys), rules)
+        @set! newsys.systems = map(s -> substitute(s, rules), get_systems(sys))
     else
         error("substituting symbols is not supported for $(typeof(sys))")
     end

src/systems/nonlinear/initializesystem.jl

@@ -41,7 +41,7 @@ end
 """
 $(TYPEDSIGNATURES)
 
-Generate `System` of nonlinear equations which initializes a problem from specified initial conditions of an `AbstractTimeDependentSystem`.
+Generate `System` of nonlinear equations which initializes a problem from specified initial conditions of a time-dependent `AbstractSystem`.
 """
 function generate_initializesystem_timevarying(sys::AbstractSystem;
         op = Dict(),
@@ -71,7 +71,6 @@ function generate_initializesystem_timevarying(sys::AbstractSystem;
     eqs_ics = Equation[]
     defs = copy(defaults(sys)) # copy so we don't modify sys.defaults
     additional_guesses = anydict(guesses)
-    additional_initialization_eqs = Vector{Equation}(initialization_eqs)
     guesses = merge(get_guesses(sys), additional_guesses)
     idxs_diff = isdiffeq.(eqs)
 
@@ -190,11 +189,13 @@ function generate_initializesystem_timevarying(sys::AbstractSystem;
     push!(pars, get_iv(sys))
 
     # 8) use observed equations for guesses of observed variables if not provided
+    guessed = Set(keys(defs)) # x(t), D(x(t)), ...
+    guessed = union(guessed, Set(default_toterm.(guessed))) # x(t), D(x(t)), xˍt(t), ...
     for eq in trueobs
-        haskey(defs, eq.lhs) && continue
-        any(x -> isequal(default_toterm(x), eq.lhs), keys(defs)) && continue
-
-        defs[eq.lhs] = eq.rhs
+        if !(eq.lhs in guessed)
+            defs[eq.lhs] = eq.rhs
+            #push!(guessed, eq.lhs) # should not encounter eq.lhs twice, so don't need to track it
+        end
     end
     append!(eqs_ics, trueobs)
 
@@ -216,7 +217,7 @@ end
 """
 $(TYPEDSIGNATURES)
 
-Generate `System` of nonlinear equations which initializes a problem from specified initial conditions of an `AbstractTimeDependentSystem`.
+Generate `System` of nonlinear equations which initializes a problem from specified initial conditions of a time-independent `AbstractSystem`.
 """
 function generate_initializesystem_timeindependent(sys::AbstractSystem;
         op = Dict(),
@@ -228,12 +229,10 @@ function generate_initializesystem_timeindependent(sys::AbstractSystem;
     eqs = equations(sys)
     trueobs, eqs = unhack_observed(observed(sys), eqs)
     vars = unique([unknowns(sys); getfield.(trueobs, :lhs)])
-    vars_set = Set(vars) # for efficient in-lookup
 
     eqs_ics = Equation[]
     defs = copy(defaults(sys)) # copy so we don't modify sys.defaults
     additional_guesses = anydict(guesses)
-    additional_initialization_eqs = Vector{Equation}(initialization_eqs)
     guesses = merge(get_guesses(sys), additional_guesses)
 
     # PREPROCESSING
@@ -738,6 +737,7 @@ function SciMLBase.late_binding_update_u0_p(
         end
         newp = setp_oop(sys, syms)(newp, vals)
     else
+        allsyms = nothing
         # if `p` is not provided or is symbolic
         p === missing || eltype(p) <: Pair || return newu0, newp
         (newu0 === nothing || isempty(newu0)) && return newu0, newp
@@ -755,6 +755,9 @@ function SciMLBase.late_binding_update_u0_p(
     if eltype(p) <: Pair
         syms = []
         vals = []
+        if allsyms === nothing
+            allsyms = all_symbols(sys)
+        end
         for (k, v) in p
             v === nothing && continue
             (symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v)) || continue

src/systems/system.jl

@@ -411,7 +411,7 @@ function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = [];
         end
         metadata = meta
     end
-    metadata = Base.ImmutableDict(metadata, MutableCacheKey => MutableCacheT())
+    metadata = refreshed_metadata(metadata)
     System(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)), eqs, noise_eqs, jumps, constraints,
         costs, consolidate, dvs, ps, brownians, iv, observed, Equation[],
         var_to_name, name, description, defaults, guesses, systems, initialization_eqs,
@@ -1097,3 +1097,52 @@ function supports_initialization(sys::System)
     return isempty(jumps(sys)) && _iszero(cost(sys)) &&
            isempty(constraints(sys))
 end
+
+safe_eachrow(::Nothing) = nothing
+safe_eachrow(x::AbstractArray) = eachrow(x)
+
+safe_issetequal(::Nothing, ::Nothing) = true
+safe_issetequal(::Nothing, x) = false
+safe_issetequal(x, ::Nothing) = false
+safe_issetequal(x, y) = issetequal(x, y)
+
+"""
+    $(TYPEDSIGNATURES)
+
+Check if two systems are about equal, to the extent that ModelingToolkit.jl supports. Note
+that if this returns `true`, the systems are not guaranteed to be exactly equivalent
+(unless `sysa === sysb`) but are highly likely to represent a similar mathematical problem.
+If this returns `false`, the systems are very likely to be different.
+"""
+function Base.isapprox(sysa::System, sysb::System)
+    sysa === sysb && return true
+    return nameof(sysa) == nameof(sysb) &&
+           isequal(get_iv(sysa), get_iv(sysb)) &&
+           issetequal(get_eqs(sysa), get_eqs(sysb)) &&
+           safe_issetequal(
+               safe_eachrow(get_noise_eqs(sysa)), safe_eachrow(get_noise_eqs(sysb))) &&
+           issetequal(get_jumps(sysa), get_jumps(sysb)) &&
+           issetequal(get_constraints(sysa), get_constraints(sysb)) &&
+           issetequal(get_costs(sysa), get_costs(sysb)) &&
+           isequal(get_consolidate(sysa), get_consolidate(sysb)) &&
+           issetequal(get_unknowns(sysa), get_unknowns(sysb)) &&
+           issetequal(get_ps(sysa), get_ps(sysb)) &&
+           issetequal(get_brownians(sysa), get_brownians(sysb)) &&
+           issetequal(get_observed(sysa), get_observed(sysb)) &&
+           issetequal(get_parameter_dependencies(sysa), get_parameter_dependencies(sysb)) &&
+           isequal(get_description(sysa), get_description(sysb)) &&
+           isequal(get_defaults(sysa), get_defaults(sysb)) &&
+           isequal(get_guesses(sysa), get_guesses(sysb)) &&
+           issetequal(get_initialization_eqs(sysa), get_initialization_eqs(sysb)) &&
+           issetequal(get_continuous_events(sysa), get_continuous_events(sysb)) &&
+           issetequal(get_discrete_events(sysa), get_discrete_events(sysb)) &&
+           isequal(get_connector_type(sysa), get_connector_type(sysb)) &&
+           isequal(get_assertions(sysa), get_assertions(sysb)) &&
+           isequal(get_metadata(sysa), get_metadata(sysb)) &&
+           isequal(get_is_dde(sysa), get_is_dde(sysb)) &&
+           issetequal(get_tstops(sysa), get_tstops(sysb)) &&
+           safe_issetequal(get_ignored_connections(sysa), get_ignored_connections(sysb)) &&
+           isequal(get_is_initializesystem(sysa), get_is_initializesystem(sysb)) &&
+           isequal(get_is_discrete(sysa), get_is_discrete(sysb)) &&
+           isequal(get_isscheduled(sysa), get_isscheduled(sysb))
+end

src/systems/systemstructure.jl

@@ -431,7 +431,7 @@ function TearingState(sys; quick_cancel = false, check = true, sort_eqs = true)
     if sort_eqs
         # sort equations lexicographically to reduce simplification issues
         # depending on order due to NP-completeness of tearing.
-        sortidxs = Base.sortperm(eqs, by = string)
+        sortidxs = Base.sortperm(string.(eqs)) # "by = string" creates more strings
         eqs = eqs[sortidxs]
         original_eqs = original_eqs[sortidxs]
         symbolic_incidence = symbolic_incidence[sortidxs]

src/utils.jl

@@ -606,7 +606,7 @@ end
 """
     $(TYPEDSIGNATURES)
 
-Indicate whether the given equation type (Equation, Pair, etc) supports `collect_vars!`. 
+Indicate whether the given equation type (Equation, Pair, etc) supports `collect_vars!`.
 Can be dispatched by higher-level libraries to indicate support.
 """
 eqtype_supports_collect_vars(eq) = false
@@ -791,7 +791,7 @@ Check if `T` is an appropriate symtype for a symbolic variable representing a fl
 point number or array of such numbers.
 """
 function is_floatingpoint_symtype(T::Type)
-    return T == Real || T == Number || T <: AbstractFloat ||
+    return T == Real || T == Number || T == Complex || T <: AbstractFloat ||
            T <: AbstractArray && is_floatingpoint_symtype(eltype(T))
 end
 

test/odesystem.jl

@@ -1574,3 +1574,28 @@ end
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0))
     @test prob.problem_type == "A"
 end
+
+@testset "`substitute` retains events and metadata" begin
+    @parameters p(t) = 1.0
+    @variables x(t) = 0.0
+    event = [0.5] => [p ~ Pre(t)]
+    event2 = [x ~ 0.75] => [p ~ 2 * Pre(t)]
+
+    struct TestMeta end
+
+    eq = [
+        D(x) ~ p
+    ]
+    @named sys = System(eq, t, [x], [p], discrete_events = [event],
+        continuous_events = [event2], metadata = Dict(TestMeta => "test"))
+
+    @variables x2(t) = 0.0
+    sys2 = substitute(sys, [x => x2])
+
+    @test length(ModelingToolkit.get_discrete_events(sys)) == 1
+    @test length(ModelingToolkit.get_discrete_events(sys2)) == 1
+    @test length(ModelingToolkit.get_continuous_events(sys)) == 1
+    @test length(ModelingToolkit.get_continuous_events(sys2)) == 1
+    @test getmetadata(sys, TestMeta, nothing) == "test"
+    @test getmetadata(sys2, TestMeta, nothing) == "test"
+end

test/serialization.jl

@@ -28,6 +28,7 @@ str = String(take!(io))
 
 sys = include_string(@__MODULE__, str)
 rc2 = expand_connections(rc_model)
+@test isapprox(sys, rc2)
 @test issetequal(equations(sys), equations(rc2))
 @test issetequal(unknowns(sys), unknowns(rc2))
 @test issetequal(parameters(sys), parameters(rc2))