Merge branch 'master' into myb/obs

Author: YingboMa

.github/workflows/Downstream.yml

@@ -18,6 +18,7 @@ jobs:
         os: [ubuntu-latest]
         package:
           - {user: SciML, repo: Catalyst.jl, group: All}
+          - {user: SciML, repo: CellMLToolkit.jl, group: All}
           - {user: SciML, repo: NeuralPDE.jl, group: NNPDE}
           - {user: SciML, repo: DataDrivenDiffEq.jl, group: Standard}
 

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "5.13.6"
+version = "5.14.2"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -51,7 +51,7 @@ Distributions = "0.23, 0.24"
 DocStringExtensions = "0.7, 0.8"
 IfElse = "0.1"
 LabelledArrays = "1.3"
-Latexify = "0.11, 0.12, 0.13, 0.14"
+Latexify = "0.11, 0.12, 0.13, 0.14, 0.15"
 LightGraphs = "1.3"
 MacroTools = "0.5"
 NaNMath = "0.3"
@@ -65,8 +65,8 @@ SciMLBase = "1.3"
 Setfield = "0.7"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
-SymbolicUtils = "0.8.3, 0.9"
-Symbolics = "0.1.8"
+SymbolicUtils = "0.11.0"
+Symbolics = "0.1.14"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
 julia = "1.2"
@@ -81,7 +81,8 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "ForwardDiff", "GalacticOptim", "NonlinearSolve", "OrdinaryDiffEq", "Optim", "Random", "SteadyStateDiffEq", "Test", "StochasticDiffEq"]
+test = ["BenchmarkTools", "ForwardDiff", "GalacticOptim", "NonlinearSolve", "OrdinaryDiffEq", "Optim", "Random", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]

docs/src/basics/ContextualVariables.md

@@ -7,7 +7,7 @@ the `@variable` which is defined by
 [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl). For example:
 
 ```julia
-@variabes x y(x)
+@variables x y(x)
 ```
 
 This is used for the "normal" variable of a given system, like the states of a

docs/src/tutorials/ode_modeling.md

@@ -13,7 +13,7 @@ But if you want to just see some code and run, here's an example:
 ```julia
 using ModelingToolkit
 
-@variables t, x(t) RHS(t)  # independent and dependent variables
+@variables t x(t) RHS(t)  # independent and dependent variables
 @parameters τ       # parameters
 D = Differential(t) # define an operator for the differentiation w.r.t. time
 
@@ -51,7 +51,7 @@ first set the forcing function to a constant value.
 ```julia
 using ModelingToolkit
 
-@variables t, x(t)  # independent and dependent variables
+@variables t x(t)  # independent and dependent variables
 @parameters τ       # parameters
 D = Differential(t) # define an operator for the differentiation w.r.t. time
 

docs/src/tutorials/tearing_parallelism.md

@@ -37,7 +37,7 @@ function connect_heat(ps...)
 end
 
 # Basic electric components
-const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
+@parameters t
 const D = Differential(t)
 function Pin(;name)
     @variables v(t) i(t)

test/rc_model.jl renamed to examples/electrical_components.jl

@@ -2,19 +2,19 @@ using Test
 using ModelingToolkit, OrdinaryDiffEq
 
 # Basic electric components
-const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
+@parameters t
 function Pin(;name)
     @variables v(t) i(t)
     ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])
 end
 
-function Ground(name)
+function Ground(;name)
     @named g = Pin()
     eqs = [g.v ~ 0]
     ODESystem(eqs, t, [], [], systems=[g], name=name)
 end
 
-function ConstantVoltage(name; V = 1.0)
+function ConstantVoltage(;name, V = 1.0)
     val = V
     @named p = Pin()
     @named n = Pin()
@@ -26,7 +26,7 @@ function ConstantVoltage(name; V = 1.0)
     ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
 end
 
-function Resistor(name; R = 1.0)
+function Resistor(;name, R = 1.0)
     val = R
     @named p = Pin()
     @named n = Pin()
@@ -40,7 +40,7 @@ function Resistor(name; R = 1.0)
     ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
 end
 
-function Capacitor(name; C = 1.0)
+function Capacitor(;name, C = 1.0)
     val = C
     @named p = Pin()
     @named n = Pin()
@@ -55,15 +55,23 @@ function Capacitor(name; C = 1.0)
     ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
 end
 
-R = 1.0
-C = 1.0
-V = 1.0
-resistor = Resistor(:resistor, R=R)
-capacitor = Capacitor(:capacitor, C=C)
-source = ConstantVoltage(:source, V=V)
-ground = Ground(:ground)
+function Inductor(; name, L = 1.0)
+    val = L
+    @named p = Pin()
+    @named n = Pin()
+    @variables v(t) i(t)
+    @parameters L
+    D = Differential(t)
+    eqs = [
+           v ~ p.v - n.v
+           0 ~ p.i + n.i
+           i ~ p.i
+           D(i) ~ v / L
+          ]
+    ODESystem(eqs, t, [v, i], [L], systems=[p, n], defaults=Dict(L => val), name=name)
+end
 
-function connect(ps...)
+function connect_pins(ps...)
     eqs = [
            0 ~ sum(p->p.i, ps) # KCL
           ]
@@ -74,10 +82,3 @@ function connect(ps...)
 
     return eqs
 end
-rc_eqs = [
-          connect(source.p, resistor.p)
-          connect(resistor.n, capacitor.p)
-          connect(capacitor.n, source.n, ground.g)
-         ]
-
-rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground], name=:rc)

examples/rc_model.jl

@@ -0,0 +1,17 @@
+include("electrical_components.jl")
+
+R = 1.0
+C = 1.0
+V = 1.0
+@named resistor = Resistor(R=R)
+@named capacitor = Capacitor(C=C)
+@named source = ConstantVoltage(V=V)
+@named ground = Ground()
+
+rc_eqs = [
+          connect_pins(source.p, resistor.p)
+          connect_pins(resistor.n, capacitor.p)
+          connect_pins(capacitor.n, source.n, ground.g)
+         ]
+
+@named rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground])

examples/serial_inductor.jl

@@ -0,0 +1,16 @@
+include("electrical_components.jl")
+
+@named source = ConstantVoltage(V=10.0)
+@named resistor = Resistor(R=1.0)
+@named inductor1 = Inductor(L=1.0e-2)
+@named inductor2 = Inductor(L=2.0e-2)
+@named ground = Ground()
+
+eqs = [
+       connect_pins(source.p, resistor.p)
+       connect_pins(resistor.n, inductor1.p)
+       connect_pins(inductor1.n, inductor2.p)
+       connect_pins(source.n, inductor2.n, ground.g)
+      ]
+
+@named ll_model = ODESystem(eqs, t, systems=[source, resistor, inductor1, inductor2, ground])

src/ModelingToolkit.jl

@@ -137,6 +137,7 @@ include("structural_transformation/StructuralTransformations.jl")
 @reexport using .StructuralTransformations
 
 export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr
+export DAEFunctionExpr, DAEProblemExpr
 export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
 export SystemStructure
 export JumpSystem
@@ -157,6 +158,7 @@ export Differential, expand_derivatives, @derivatives
 export IntervalDomain, ProductDomain, ⊗, CircleDomain
 export Equation, ConstrainedEquation
 export Term, Sym
+export SymScope, LocalScope, ParentScope, GlobalScope
 export independent_variable, states, parameters, equations, controls, observed, structure
 export structural_simplify
 

src/bipartite_graph.jl

@@ -66,13 +66,13 @@ badjlist = [[1,2,5,6],[3,4,6]]
 bg = BipartiteGraph(7, fadjlist, badjlist)
 ```
 """
-mutable struct BipartiteGraph{I<:Integer,M} <: LightGraphs.AbstractGraph{I}
+mutable struct BipartiteGraph{I<:Integer,F<:Vector{Vector{I}},B<:Union{Vector{Vector{I}},I},M} <: LightGraphs.AbstractGraph{I}
     ne::Int
-    fadjlist::Vector{Vector{I}} # `fadjlist[src] => dsts`
-    badjlist::Vector{Vector{I}} # `badjlist[dst] => srcs`
+    fadjlist::F # `fadjlist[src] => dsts`
+    badjlist::B # `badjlist[dst] => srcs` or `ndsts`
     metadata::M
 end
-BipartiteGraph(ne::Integer, fadj::AbstractVector, badj::AbstractVector) = BipartiteGraph(ne, fadj, badj, nothing)
+BipartiteGraph(ne::Integer, fadj::AbstractVector, badj::Union{AbstractVector,Integer}=maximum(maximum, fadj); metadata=nothing) = BipartiteGraph(ne, fadj, badj, metadata)
 
 """
 ```julia
@@ -93,16 +93,16 @@ $(SIGNATURES)
 
 Build an empty `BipartiteGraph` with `nsrcs` sources and `ndsts` destinations.
 """
-function BipartiteGraph(nsrcs::T, ndsts::T; metadata=nothing) where T
+function BipartiteGraph(nsrcs::T, ndsts::T, backedge::Val{B}=Val(true); metadata=nothing) where {T,B}
     fadjlist = map(_->T[], 1:nsrcs)
-    badjlist = map(_->T[], 1:ndsts)
+    badjlist = B ? map(_->T[], 1:ndsts) : ndsts
     BipartiteGraph(0, fadjlist, badjlist, metadata)
 end
 
 Base.eltype(::Type{<:BipartiteGraph{I}}) where I = I
 function Base.empty!(g::BipartiteGraph)
     foreach(empty!, g.fadjlist)
-    foreach(empty!, g.badjlist)
+    g.badjlist isa AbstractVector && foreach(empty!, g.badjlist)
     g.ne = 0
     if g.metadata !== nothing
         foreach(empty!, g.metadata)
@@ -111,17 +111,22 @@ function Base.empty!(g::BipartiteGraph)
 end
 Base.length(::BipartiteGraph) = error("length is not well defined! Use `ne` or `nv`.")
 
+@noinline throw_no_back_edges() = throw(ArgumentError("The graph has no back edges."))
+
 if isdefined(LightGraphs, :has_contiguous_vertices)
     LightGraphs.has_contiguous_vertices(::Type{<:BipartiteGraph}) = false
 end
 LightGraphs.is_directed(::Type{<:BipartiteGraph}) = false
 LightGraphs.vertices(g::BipartiteGraph) = (𝑠vertices(g), 𝑑vertices(g))
 𝑠vertices(g::BipartiteGraph) = axes(g.fadjlist, 1)
-𝑑vertices(g::BipartiteGraph) = axes(g.badjlist, 1)
+𝑑vertices(g::BipartiteGraph) = g.badjlist isa AbstractVector ? axes(g.badjlist, 1) : Base.OneTo(g.badjlist)
 has_𝑠vertex(g::BipartiteGraph, v::Integer) = v in 𝑠vertices(g)
 has_𝑑vertex(g::BipartiteGraph, v::Integer) = v in 𝑑vertices(g)
 𝑠neighbors(g::BipartiteGraph, i::Integer, with_metadata::Val{M}=Val(false)) where M = M ? zip(g.fadjlist[i], g.metadata[i]) : g.fadjlist[i]
-𝑑neighbors(g::BipartiteGraph, j::Integer, with_metadata::Val{M}=Val(false)) where M = M ? zip(g.badjlist[j], (g.metadata[i][j] for i in g.badjlist[j])) : g.badjlist[j]
+function 𝑑neighbors(g::BipartiteGraph, j::Integer, with_metadata::Val{M}=Val(false)) where M
+    g.badjlist isa AbstractVector || throw_no_back_edges()
+    M ? zip(g.badjlist[j], (g.metadata[i][j] for i in g.badjlist[j])) : g.badjlist[j]
+end
 LightGraphs.ne(g::BipartiteGraph) = g.ne
 LightGraphs.nv(g::BipartiteGraph) = sum(length, vertices(g))
 LightGraphs.edgetype(g::BipartiteGraph{I}) where I = BipartiteEdge{I}
@@ -145,7 +150,6 @@ const NO_METADATA = NoMetadata()
 LightGraphs.add_edge!(g::BipartiteGraph, i::Integer, j::Integer, md=NO_METADATA) = add_edge!(g, BipartiteEdge(i, j), md)
 function LightGraphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge, md=NO_METADATA)
     @unpack fadjlist, badjlist = g
-    verts = vertices(g)
     s, d = src(edge), dst(edge)
     (has_𝑠vertex(g, s) && has_𝑑vertex(g, d)) || error("edge ($edge) out of range.")
     @inbounds list = fadjlist[s]
@@ -157,15 +161,21 @@ function LightGraphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge, md=NO_MET
     end
 
     g.ne += 1
-    @inbounds list = badjlist[d]
-    index = searchsortedfirst(list, s)
-    insert!(list, index, s)
+    if badjlist isa AbstractVector
+        @inbounds list = badjlist[d]
+        index = searchsortedfirst(list, s)
+        insert!(list, index, s)
+    end
     return true  # edge successfully added
 end
 
 function LightGraphs.add_vertex!(g::BipartiteGraph{T}, type::VertType) where T
     if type === DST
-        push!(g.badjlist, T[])
+        if g.badjlist isa AbstractVector
+            push!(g.badjlist, T[])
+        else
+            g.badjlist += 1
+        end
     elseif type === SRC
         push!(g.fadjlist, T[])
     else