Merge pull request #2726 from SciML/dw/backport_odefunctionexpr

Backport #2719 to MTK v8

.github/workflows/ci.yml

@@ -1,8 +1,6 @@
 name: CI
 on:
   pull_request:
-    branches:
-      - master
     paths-ignore:
       - 'docs/**'
   push:

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 = "8.75.0"
+version = "8.76.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -99,7 +99,7 @@ SparseArrays = "1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 SymbolicIndexingInterface = "0.3.1"
-SymbolicUtils = "1.0"
+SymbolicUtils = "1.0 - 1.5.1"
 Symbolics = "5.7"
 URIs = "1"
 UnPack = "0.1, 1.0"

docs/make.jl

@@ -16,7 +16,7 @@ mathengine = MathJax3(Dict(:loader => Dict("load" => ["[tex]/require", "[tex]/ma
             "ams",
             "autoload",
             "mathtools",
-            "require",
+            "require"
         ])))
 
 makedocs(sitename = "ModelingToolkit.jl",

docs/pages.jl

@@ -10,7 +10,8 @@ pages = [
         "tutorials/parameter_identifiability.md",
         "tutorials/bifurcation_diagram_computation.md",
         "tutorials/domain_connections.md"],
-    "Examples" => Any["Basic Examples" => Any["examples/higher_order.md",
+    "Examples" => Any[
+        "Basic Examples" => Any["examples/higher_order.md",
             "examples/spring_mass.md",
             "examples/modelingtoolkitize_index_reduction.md",
             "examples/parsing.md"],
@@ -35,5 +36,5 @@ pages = [
         "systems/DiscreteSystem.md",
         "systems/PDESystem.md"],
     "comparison.md",
-    "internals.md",
+    "internals.md"
 ]

docs/src/basics/Composition.md

@@ -22,7 +22,7 @@ function decay(; name)
     @variables x(t) f(t)
     D = Differential(t)
     ODESystem([
-            D(x) ~ -a * x + f,
+            D(x) ~ -a * x + f
         ];
         name = name)
 end
@@ -32,8 +32,9 @@ end
 
 @parameters t
 D = Differential(t)
-connected = compose(ODESystem([decay2.f ~ decay1.x
-            D(decay1.f) ~ 0], t; name = :connected), decay1, decay2)
+connected = compose(
+    ODESystem([decay2.f ~ decay1.x
+               D(decay1.f) ~ 0], t; name = :connected), decay1, decay2)
 
 equations(connected)
 
@@ -52,10 +53,10 @@ Now we can solve the system:
 
 ```@example composition
 x0 = [decay1.x => 1.0
-    decay1.f => 0.0
-    decay2.x => 1.0]
+      decay1.f => 0.0
+      decay2.x => 1.0]
 p = [decay1.a => 0.1
-    decay2.a => 0.2]
+     decay2.a => 0.2]
 
 using DifferentialEquations
 prob = ODEProblem(simplified_sys, x0, (0.0, 100.0), p)
@@ -96,7 +97,7 @@ in their namespaced form. For example:
 
 ```julia
 u0 = [x => 2.0
-    subsys.x => 2.0]
+      subsys.x => 2.0]
 ```
 
 Note that any default values within the given subcomponent will be
@@ -210,19 +211,27 @@ N = S + I + R
 @named ieqn = ODESystem([D(I) ~ β * S * I / N - γ * I])
 @named reqn = ODESystem([D(R) ~ γ * I])
 
-sir = compose(ODESystem([
+sir = compose(
+    ODESystem(
+        [
             S ~ ieqn.S,
             I ~ seqn.I,
             R ~ ieqn.R,
             ieqn.S ~ seqn.S,
             seqn.I ~ ieqn.I,
             seqn.R ~ reqn.R,
             ieqn.R ~ reqn.R,
-            reqn.I ~ ieqn.I], t, [S, I, R], [β, γ],
+            reqn.I ~ ieqn.I],
+        t,
+        [S, I, R],
+        [β, γ],
         defaults = [seqn.β => β
-            ieqn.β => β
-            ieqn.γ => γ
-            reqn.γ => γ], name = :sir), seqn, ieqn, reqn)
+                    ieqn.β => β
+                    ieqn.γ => γ
+                    reqn.γ => γ], name = :sir),
+    seqn,
+    ieqn,
+    reqn)
 ```
 
 Note that the states are forwarded by an equality relationship, while
@@ -244,7 +253,7 @@ u0 = [seqn.S => 990.0,
     reqn.R => 0.0]
 
 p = [β => 0.5
-    γ => 0.25]
+     γ => 0.25]
 
 tspan = (0.0, 40.0)
 prob = ODEProblem(sireqn_simple, u0, tspan, p, jac = true)

docs/src/basics/Events.md

@@ -69,7 +69,7 @@ function UnitMassWithFriction(k; name)
     @variables t x(t)=0 v(t)=0
     D = Differential(t)
     eqs = [D(x) ~ v
-        D(v) ~ sin(t) - k * sign(v)]
+           D(v) ~ sin(t) - k * sign(v)]
     ODESystem(eqs, t; continuous_events = [v ~ 0], name) # when v = 0 there is a discontinuity
 end
 @named m = UnitMassWithFriction(0.7)
@@ -94,7 +94,7 @@ root_eqs = [x ~ 0]  # the event happens at the ground x(t) = 0
 affect = [v ~ -v] # the effect is that the velocity changes sign
 
 @named ball = ODESystem([D(x) ~ v
-        D(v) ~ -9.8], t; continuous_events = root_eqs => affect) # equation => affect
+                         D(v) ~ -9.8], t; continuous_events = root_eqs => affect) # equation => affect
 
 ball = structural_simplify(ball)
 
@@ -114,13 +114,14 @@ Multiple events? No problem! This example models a bouncing ball in 2D that is e
 D = Differential(t)
 
 continuous_events = [[x ~ 0] => [vx ~ -vx]
-    [y ~ -1.5, y ~ 1.5] => [vy ~ -vy]]
+                     [y ~ -1.5, y ~ 1.5] => [vy ~ -vy]]
 
-@named ball = ODESystem([
+@named ball = ODESystem(
+    [
         D(x) ~ vx,
         D(y) ~ vy,
         D(vx) ~ -9.8 - 0.1vx, # gravity + some small air resistance
-        D(vy) ~ -0.1vy,
+        D(vy) ~ -0.1vy
     ], t; continuous_events)
 
 ball = structural_simplify(ball)
@@ -189,7 +190,7 @@ affect interface:
 sts = @variables x(t), v(t)
 par = @parameters g = 9.8
 bb_eqs = [D(x) ~ v
-    D(v) ~ -g]
+          D(v) ~ -g]
 
 function bb_affect!(integ, u, p, ctx)
     integ.u[u.v] = -integ.u[u.v]

docs/src/basics/Linearization.md

@@ -26,8 +26,8 @@ using ModelingToolkit
 D = Differential(t)
 
 eqs = [u ~ kp * (r - y) # P controller
-    D(x) ~ -x + u    # First-order plant
-    y ~ x]           # Output equation
+       D(x) ~ -x + u    # First-order plant
+       y ~ x]           # Output equation
 
 @named sys = ODESystem(eqs, t)
 matrices, simplified_sys = linearize(sys, [r], [y]) # Linearize from r to y

docs/src/basics/Variable_metadata.md

@@ -142,7 +142,7 @@ Dₜ = Differential(t)
 @parameters T [tunable = true, bounds = (0, Inf)]
 @parameters k [tunable = true, bounds = (0, Inf)]
 eqs = [Dₜ(x) ~ (-x + k * u) / T # A first-order system with time constant T and gain k
-    y ~ x]
+       y ~ x]
 sys = ODESystem(eqs, t, name = :tunable_first_order)
 ```
 

docs/src/examples/parsing.md

@@ -7,8 +7,8 @@ symbolic forms from that? For example, say we had the following system we wanted
 
 ```@example parsing
 ex = [:(y ~ x)
-    :(y ~ -2x + 3 / z)
-    :(z ~ 2)]
+      :(y ~ -2x + 3 / z)
+      :(z ~ 2)]
 ```
 
 We can use the function `parse_expr_to_symbolic` from Symbolics.jl to generate the symbolic

docs/src/examples/spring_mass.md

@@ -26,7 +26,7 @@ end
 
 function connect_spring(spring, a, b)
     [spring.x ~ norm(scalarize(a .- b))
-        scalarize(spring.dir .~ scalarize(a .- b))]
+     scalarize(spring.dir .~ scalarize(a .- b))]
 end
 
 function spring_force(spring)
@@ -43,7 +43,7 @@ g = [0.0, -9.81]
 @named spring = Spring(k = k, l = l)
 
 eqs = [connect_spring(spring, mass.pos, center)
-    scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)]
+       scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)]
 
 @named _model = ODESystem(eqs, t, [spring.x; spring.dir; mass.pos], [])
 @named model = compose(_model, mass, spring)
@@ -96,7 +96,7 @@ We now define functions that help construct the equations for a mass-spring syst
 ```@example component
 function connect_spring(spring, a, b)
     [spring.x ~ norm(scalarize(a .- b))
-        scalarize(spring.dir .~ scalarize(a .- b))]
+     scalarize(spring.dir .~ scalarize(a .- b))]
 end
 ```
 
@@ -125,7 +125,7 @@ We can now create the equations describing this system, by connecting `spring` t
 
 ```@example component
 eqs = [connect_spring(spring, mass.pos, center)
-    scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)]
+       scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)]
 ```
 
 Finally, we can build the model using these equations and components.

docs/src/examples/tearing_parallelism.md

@@ -32,7 +32,7 @@ function ConstantVoltage(; name, V = 1.0)
     @named n = Pin()
     @parameters V = V
     eqs = [V ~ p.v - n.v
-        0 ~ p.i + n.i]
+           0 ~ p.i + n.i]
     compose(ODESystem(eqs, t, [], [V], name = name), p, n)
 end
 
@@ -48,10 +48,10 @@ function HeatingResistor(; name, R = 1.0, TAmbient = 293.15, alpha = 1.0)
     @variables v(t) RTherm(t)
     @parameters R=R TAmbient=TAmbient alpha=alpha
     eqs = [RTherm ~ R * (1 + alpha * (h.T - TAmbient))
-        v ~ p.i * RTherm
-        h.Q_flow ~ -v * p.i # -LossPower
-        v ~ p.v - n.v
-        0 ~ p.i + n.i]
+           v ~ p.i * RTherm
+           h.Q_flow ~ -v * p.i # -LossPower
+           v ~ p.v - n.v
+           0 ~ p.i + n.i]
     compose(ODESystem(eqs, t, [v, RTherm], [R, TAmbient, alpha],
             name = name), p, n, h)
 end
@@ -61,7 +61,7 @@ function HeatCapacitor(; name, rho = 8050, V = 1, cp = 460, TAmbient = 293.15)
     C = rho * V * cp
     @named h = HeatPort()
     eqs = [
-        D(h.T) ~ h.Q_flow / C,
+        D(h.T) ~ h.Q_flow / C
     ]
     compose(ODESystem(eqs, t, [], [rho, V, cp],
             name = name), h)
@@ -73,8 +73,8 @@ function Capacitor(; name, C = 1.0)
     @variables v(t) = 0.0
     @parameters C = C
     eqs = [v ~ p.v - n.v
-        0 ~ p.i + n.i
-        D(v) ~ p.i / C]
+           0 ~ p.i + n.i
+           D(v) ~ p.i / C]
     compose(ODESystem(eqs, t, [v], [C],
             name = name), p, n)
 end
@@ -85,9 +85,9 @@ function parallel_rc_model(i; name, source, ground, R, C)
     heat_capacitor = HeatCapacitor(name = Symbol(:heat_capacitor, i))
 
     rc_eqs = [connect(source.p, resistor.p)
-        connect(resistor.n, capacitor.p)
-        connect(capacitor.n, source.n, ground.g)
-        connect(resistor.h, heat_capacitor.h)]
+              connect(resistor.n, capacitor.p)
+              connect(capacitor.n, source.n, ground.g)
+              connect(resistor.h, heat_capacitor.h)]
 
     compose(ODESystem(rc_eqs, t, name = Symbol(name, i)),
         [resistor, capacitor, source, ground, heat_capacitor])
@@ -113,7 +113,7 @@ end;
 @variables E(t) = 0.0
 eqs = [
     D(E) ~ sum(((i, sys),) -> getproperty(sys, Symbol(:resistor, i)).h.Q_flow,
-        enumerate(rc_systems)),
+    enumerate(rc_systems))
 ]
 @named _big_rc = ODESystem(eqs, t, [E], [])
 @named big_rc = compose(_big_rc, rc_systems)
@@ -155,7 +155,8 @@ ts = TearingState(expand_connections(big_rc))
 inc_org = BipartiteGraphs.incidence_matrix(ts.structure.graph)
 blt_org = StructuralTransformations.sorted_incidence_matrix(ts, only_algeqs = true,
     only_algvars = true)
-blt_reduced = StructuralTransformations.sorted_incidence_matrix(ModelingToolkit.get_tearing_state(sys),
+blt_reduced = StructuralTransformations.sorted_incidence_matrix(
+    ModelingToolkit.get_tearing_state(sys),
     only_algeqs = true,
     only_algvars = true)
 ```

docs/src/tutorials/acausal_components.md

@@ -102,7 +102,7 @@ end
 
 @mtkbuild rc_model = RCModel(resistor.R = 2.0)
 u0 = [
-    rc_model.capacitor.v => 0.0,
+    rc_model.capacitor.v => 0.0
 ]
 prob = ODEProblem(rc_model, u0, (0, 10.0))
 sol = solve(prob)
@@ -316,7 +316,7 @@ This is done as follows:
 
 ```@example acausal
 u0 = [rc_model.capacitor.v => 0.0
-    rc_model.capacitor.p.i => 0.0]
+      rc_model.capacitor.p.i => 0.0]
 
 prob = ODEProblem(rc_model, u0, (0, 10.0))
 sol = solve(prob)