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 <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "10.5.0"
+version = "10.6.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

benchmark/benchmarks.jl

@@ -4,6 +4,7 @@ using ModelingToolkitStandardLibrary.Electrical
 using ModelingToolkitStandardLibrary.Mechanical.Rotational
 using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEqDefault
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 const SUITE = BenchmarkGroup()
 
@@ -45,12 +46,33 @@ end
 
 @named model = DCMotor()
 
+# first call
+mtkcompile(model)
 SUITE["mtkcompile"] = @benchmarkable mtkcompile($model)
 
 model = mtkcompile(model)
 u0 = unknowns(model) .=> 0.0
 tspan = (0.0, 6.0)
-SUITE["ODEProblem"] = @benchmarkable ODEProblem($model, $u0, $tspan)
 
 prob = ODEProblem(model, u0, tspan)
+SUITE["ODEProblem"] = @benchmarkable ODEProblem($model, $u0, $tspan)
+
+# first call
+init(prob)
 SUITE["init"] = @benchmarkable init($prob)
+
+large_param_init = SUITE["large_parameter_init"] = BenchmarkGroup()
+
+N = 25
+@variables x(t)[1:N]
+@parameters A[1:N, 1:N]
+
+defval = collect(x) * collect(x)'
+@mtkcompile model = System(
+    [D(x) ~ x], t, [x], [A]; defaults = [A => defval], guesses = [A => fill(NaN, N, N)])
+
+u0 = [x => rand(N)]
+prob = ODEProblem(model, u0, tspan)
+large_param_init["ODEProblem"] = @benchmarkable ODEProblem($model, $u0, $tspan)
+
+large_param_init["init"] = @benchmarkable init($prob)

src/bipartite_graph.jl

@@ -535,13 +535,39 @@ function set_neighbors!(g::BipartiteGraph, i::Integer, new_neighbors)
     end
 end
 
-function delete_srcs!(g::BipartiteGraph, srcs)
+function delete_srcs!(g::BipartiteGraph{I}, srcs; rm_verts = false) where {I}
     for s in srcs
         set_neighbors!(g, s, ())
     end
+    if rm_verts
+        old_to_new_idxs = collect(one(I):I(nsrcs(g)))
+        for s in srcs
+            old_to_new_idxs[s] = zero(I)
+        end
+        offset = zero(I)
+        for i in eachindex(old_to_new_idxs)
+            if iszero(old_to_new_idxs[i])
+                offset += one(I)
+                continue
+            end
+            old_to_new_idxs[i] -= offset
+        end
+
+        if g.badjlist isa AbstractVector
+            for i in 1:ndsts(g)
+                for j in eachindex(g.badjlist[i])
+                    g.badjlist[i][j] = old_to_new_idxs[g.badjlist[i][j]]
+                end
+                filter!(!iszero, g.badjlist[i])
+            end
+        end
+        deleteat!(g.fadjlist, srcs)
+    end
     g
 end
-delete_dsts!(g::BipartiteGraph, srcs) = delete_srcs!(invview(g), srcs)
+function delete_dsts!(g::BipartiteGraph, srcs; rm_verts = false)
+    delete_srcs!(invview(g), srcs; rm_verts)
+end
 
 ###
 ### Edges iteration

src/problems/initializationproblem.jl

@@ -39,7 +39,7 @@ All other keyword arguments are forwarded to the wrapped nonlinear problem const
     for k in keys(op)
         has_u0_ics |= is_variable(sys, k) || isdifferential(k) ||
                       symbolic_type(k) == ArraySymbolic() &&
-                      is_sized_array_symbolic(k) && is_variable(sys, first(collect(k)))
+                      is_sized_array_symbolic(k) && is_variable(sys, unwrap(first(wrap(k))))
     end
     if !has_u0_ics && get_initializesystem(sys) !== nothing
         isys = get_initializesystem(sys; initialization_eqs, check_units)
@@ -79,7 +79,7 @@ All other keyword arguments are forwarded to the wrapped nonlinear problem const
         @warn errmsg
     end
 
-    uninit = setdiff(unknowns(sys), [unknowns(isys); observables(isys)])
+    uninit = setdiff(unknowns(sys), unknowns(isys), observables(isys))
 
     # TODO: throw on uninitialized arrays
     filter!(x -> !(x isa Symbolics.Arr), uninit)

src/structural_transformation/symbolics_tearing.jl

@@ -960,7 +960,8 @@ function update_simplified_system!(
         obs_sub[eq.lhs] = eq.rhs
     end
     # TODO: compute the dependency correctly so that we don't have to do this
-    obs = [fast_substitute(observed(sys), obs_sub); solved_eqs]
+    obs = [fast_substitute(observed(sys), obs_sub); solved_eqs;
+           fast_substitute(state.additional_observed, obs_sub)]
 
     unknown_idxs = filter(
         i -> diff_to_var[i] === nothing && ispresent(i) && !(fullvars[i] in solved_vars), eachindex(state.fullvars))

src/systems/abstractsystem.jl

@@ -805,7 +805,23 @@ has_equations(::AbstractSystem) = true
 
 Invalidate cached jacobians, etc.
 """
-invalidate_cache!(sys::AbstractSystem) = sys
+function invalidate_cache!(sys::AbstractSystem)
+    has_metadata(sys) || return sys
+    empty!(getmetadata(sys, MutableCacheKey, nothing))
+    return sys
+end
+
+# `::MetadataT` but that is defined later
+function refreshed_metadata(meta::Base.ImmutableDict)
+    newmeta = MetadataT()
+    for (k, v) in meta
+        if k === MutableCacheKey
+            v = MutableCacheT()
+        end
+        newmeta = Base.ImmutableDict(newmeta, k => v)
+    end
+    return newmeta
+end
 
 function Setfield.get(obj::AbstractSystem, ::Setfield.PropertyLens{field}) where {field}
     getfield(obj, field)
@@ -815,6 +831,8 @@ end
         args = map(fieldnames(obj)) do fn
             if fn in fieldnames(patch)
                 :(patch.$fn)
+            elseif fn == :metadata
+                :($refreshed_metadata(getfield(obj, $(Meta.quot(fn)))))
             else
                 :(getfield(obj, $(Meta.quot(fn))))
             end
@@ -2507,7 +2525,15 @@ function extend(sys::AbstractSystem, basesys::AbstractSystem;
     cevs = union(get_continuous_events(basesys), get_continuous_events(sys))
     devs = union(get_discrete_events(basesys), get_discrete_events(sys))
     defs = merge(get_defaults(basesys), get_defaults(sys)) # prefer `sys`
-    meta = merge(get_metadata(basesys), get_metadata(sys))
+    meta = MetadataT()
+    for kvp in get_metadata(basesys)
+        kvp[1] == MutableCacheKey && continue
+        meta = Base.ImmutableDict(meta, kvp)
+    end
+    for kvp in get_metadata(sys)
+        kvp[1] == MutableCacheKey && continue
+        meta = Base.ImmutableDict(meta, kvp)
+    end
     syss = union(get_systems(basesys), get_systems(sys))
     args = length(ivs) == 0 ? (eqs, sts, ps) : (eqs, ivs[1], sts, ps)
     kwargs = (observed = obs, continuous_events = cevs,
@@ -2705,7 +2731,9 @@ function process_parameter_equations(sys::AbstractSystem)
                     is_sized_array_symbolic(sym) &&
                     all(Base.Fix1(is_parameter, sys), collect(sym))
         end
-            if !isparameter(eq.lhs)
+            # Everything in `varsbuf` is a parameter, so this is a cheap `is_parameter`
+            # check.
+            if !(eq.lhs in varsbuf)
                 throw(ArgumentError("""
                 LHS of parameter dependency equation must be a single parameter. Found \
                 $(eq.lhs).

src/systems/connectors.jl

@@ -872,7 +872,8 @@ function expand_connections(sys::AbstractSystem; tol = 1e-10)
     eqs = [equations(sys); ceqs; stream_eqs]
     # substitute `instream(..)` expressions with their new values
     for i in eachindex(eqs)
-        eqs[i] = fixpoint_sub(eqs[i], instream_subs; maxiters = length(instream_subs))
+        eqs[i] = fixpoint_sub(
+            eqs[i], instream_subs; maxiters = max(length(instream_subs), 10))
     end
     # get the defaults for domain networks
     d_defs = domain_defaults(sys, domain_csets)