Merge remote-tracking branch 'origin/master' into event_discpars_autoinfer

Author: TorkelE

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.26.1"
+version = "10.28.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -26,6 +26,7 @@ DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 ExprTools = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
@@ -44,6 +45,7 @@ NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -114,6 +116,7 @@ DynamicQuantities = "^0.11.2, 0.12, 0.13, 1"
 EnumX = "1.0.4"
 ExprTools = "0.1.10"
 FMI = "0.14"
+FillArrays = "1.13.0"
 FindFirstFunctions = "1"
 ForwardDiff = "0.10.3, 1"
 FunctionWrappers = "1.1"
@@ -141,6 +144,7 @@ OrdinaryDiffEq = "6.82.0"
 OrdinaryDiffEqCore = "1.34.0"
 OrdinaryDiffEqDefault = "1.2"
 OrdinaryDiffEqNonlinearSolve = "1.5.0"
+PreallocationTools = "0.4.27"
 PrecompileTools = "1"
 Pyomo = "0.1.0"
 REPL = "1"

README.md

@@ -11,7 +11,7 @@
 [![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)
 
 ModelingToolkit.jl is a modeling framework for high-performance symbolic-numeric computation
-in scientific computing  and scientific machine learning.
+in scientific computing and scientific machine learning.
 It allows for users to give a high-level description of a model for
 symbolic preprocessing to analyze and enhance the model. ModelingToolkit can
 automatically generate fast functions for model components like Jacobians
@@ -36,88 +36,95 @@ lower it to a first order system, symbolically generate the Jacobian function
 for the numerical integrator, and solve it.
 
 ```julia
-using OrdinaryDiffEqDefault, ModelingToolkit
+using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
+# Defines a ModelingToolkit `System` model.
 @parameters σ ρ β
 @variables x(t) y(t) z(t)
-
-eqs = [D(D(x)) ~ σ * (y - x),
+eqs = [
+    D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
-    D(z) ~ x * y - β * z]
-
+    D(z) ~ x * y - β * z
+]
 @mtkcompile sys = System(eqs, t)
 
-u0 = [D(x) => 2.0,
+# Simulate the model for a specific condition (initial condition and parameter values).
+using OrdinaryDiffEqDefault
+sim_cond = [
+    D(x) => 2.0,
     x => 1.0,
     y => 0.0,
-    z => 0.0]
-
-p = [σ => 28.0,
+    z => 0.0,
+    σ => 28.0,
     ρ => 10.0,
-    β => 8 / 3]
-
-tspan = (0.0, 100.0)
-prob = ODEProblem(sys, u0, tspan, p, jac = true)
+    β => 8 / 3
+]
+tend = 100.0
+prob = ODEProblem(sys, sim_cond, tend; jac = true)
 sol = solve(prob)
+
+# Plot the solution in phase-space.
 using Plots
 plot(sol, idxs = (x, y))
 ```
 
-![Lorenz2](https://user-images.githubusercontent.com/1814174/79118645-744eb580-7d5c-11ea-9c37-13c4efd585ca.png)
+![Lorenz2](https://github.com/user-attachments/assets/e82fb2ce-97b7-4f56-b272-85653c88bdb3)
 
-This automatically will have generated fast Jacobian functions, making
+This will have automatically generated fast Jacobian functions, making
 it more optimized than directly building a function. In addition, we can then
 use ModelingToolkit to compose multiple ODE subsystems. Now, let's define two
 interacting Lorenz equations and simulate the resulting Differential-Algebraic
 Equation (DAE):
 
 ```julia
-using DifferentialEquations, ModelingToolkit
+using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
-@parameters σ ρ β
-@variables x(t) y(t) z(t)
-
-eqs = [D(x) ~ σ * (y - x),
+# Defines two lorenz system models.
+eqs = [
+    D(x) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
-    D(z) ~ x * y - β * z]
-
+    D(z) ~ x * y - β * z
+]
 @named lorenz1 = System(eqs, t)
 @named lorenz2 = System(eqs, t)
 
+# Connect the two models, creating a single model.
 @variables a(t)
 @parameters γ
 connections = [0 ~ lorenz1.x + lorenz2.y + a * γ]
-@mtkcompile connected = System(connections, t, systems = [lorenz1, lorenz2])
+@mtkcompile connected_lorenz = System(connections, t; systems = [lorenz1, lorenz2])
 
-u0 = [lorenz1.x => 1.0,
+# Simulate the model for a specific condition (initial condition and parameter values).
+using OrdinaryDiffEqDefault
+sim_cond = [
+    lorenz1.x => 1.0,
     lorenz1.y => 0.0,
     lorenz1.z => 0.0,
     lorenz2.x => 0.0,
-    lorenz2.y => 1.0,
     lorenz2.z => 0.0,
-    a => 2.0]
-
-p = [lorenz1.σ => 10.0,
+    a => 2.0,
+    lorenz1.σ => 10.0,
     lorenz1.ρ => 28.0,
     lorenz1.β => 8 / 3,
     lorenz2.σ => 10.0,
     lorenz2.ρ => 28.0,
     lorenz2.β => 8 / 3,
-    γ => 2.0]
-
-tspan = (0.0, 100.0)
-prob = ODEProblem(connected, u0, tspan, p)
+    γ => 2.0
+]
+tend = 100.0
+prob = ODEProblem(connected_lorenz, sim_cond, tend)
 sol = solve(prob)
 
+# Plot the solution in phase-space.
 using Plots
 plot(sol, idxs = (a, lorenz1.x, lorenz2.z))
 ```
 
-![](https://user-images.githubusercontent.com/17304743/187790221-528046c3-dbdb-4853-b977-799596c147f3.png)
+![LorenzConnected](https://github.com/user-attachments/assets/ef65d812-c10e-42c6-945d-e61515e3b6a1)
 
-# Citation
+## Citation
 
 If you use ModelingToolkit.jl in your research, please cite [this paper](https://arxiv.org/abs/2103.05244):
 

docs/src/API/System.md

@@ -214,7 +214,7 @@ ModelingToolkit.eqtype_supports_collect_vars
 ModelingToolkit.modified_unknowns!
 ```
 
-## Namepsace manipulation
+## Namespace manipulation
 
 ModelingToolkit namespaces variables from subsystems when using them in a parent system to
 disambiguate from identically named variables in other subsystems or the parent system. The

docs/src/API/codegen.md

@@ -23,6 +23,9 @@ ModelingToolkit.build_explicit_observed_function
 ModelingToolkit.generate_control_function
 ModelingToolkit.generate_update_A
 ModelingToolkit.generate_update_b
+ModelingToolkit.generate_semiquadratic_functions
+ModelingToolkit.generate_semiquadratic_jacobian
+ModelingToolkit.get_semiquadratic_W_sparsity
 ```
 
 For functions such as jacobian calculation which require symbolic computation, there

docs/src/API/problems.md

@@ -17,6 +17,8 @@ SciMLBase.ODEFunction
 SciMLBase.ODEProblem
 SciMLBase.DAEFunction
 SciMLBase.DAEProblem
+ModelingToolkit.SemilinearODEFunction
+ModelingToolkit.SemilinearODEProblem
 SciMLBase.SDEFunction
 SciMLBase.SDEProblem
 SciMLBase.DDEFunction

docs/src/basics/Composition.md

@@ -57,7 +57,7 @@ p = [decay1.a => 0.1
      decay2.a => 0.2]
 
 using OrdinaryDiffEq
-prob = ODEProblem(simplified_sys, x0, (0.0, 100.0), p)
+prob = ODEProblem(simplified_sys, [x0; p], (0.0, 100.0))
 sol = solve(prob, Tsit5())
 sol[decay2.f]
 ```

docs/src/basics/Events.md

@@ -418,7 +418,7 @@ ev = ModelingToolkit.SymbolicDiscreteCallback(
 @mtkcompile sys = System(
     D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [ev])
 
-prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0])
+prob = ODEProblem(sys, [x => 0.0, c => 1.0], (0.0, 2pi))
 sol = solve(prob, Tsit5())
 sol[c]
 ```
@@ -439,7 +439,7 @@ will be saved. If we repeat the above example with `c` not a `discrete_parameter
 @mtkcompile sys = System(
     D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ Pre(c) + 1]])
 
-prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0])
+prob = ODEProblem(sys, [x => 0.0, c => 1.0], (0.0, 2pi))
 sol = solve(prob, Tsit5())
 sol.ps[c] # sol[c] will error, since `c` is not a timeseries value
 ```

docs/src/tutorials/ode_modeling.md

@@ -36,7 +36,7 @@ end
 
 using OrdinaryDiffEq
 @mtkcompile fol = FOL()
-prob = ODEProblem(fol, [], (0.0, 10.0), [])
+prob = ODEProblem(fol, [], (0.0, 10.0))
 sol = solve(prob)
 
 using Plots

docs/src/tutorials/programmatically_generating.md

@@ -49,7 +49,7 @@ eqs = [D(x) ~ (h - x) / τ] # create an array of equations
 # Perform the standard transformations and mark the model complete
 # Note: Complete models cannot be subsystems of other models!
 fol = mtkcompile(model)
-prob = ODEProblem(fol, [], (0.0, 10.0), [])
+prob = ODEProblem(fol, [], (0.0, 10.0))
 using OrdinaryDiffEq
 sol = solve(prob)
 

src/ModelingToolkit.jl

@@ -104,6 +104,9 @@ const DQ = DynamicQuantities
 import DifferentiationInterface as DI
 using ADTypes: AutoForwardDiff
 import SciMLPublic: @public
+import PreallocationTools
+import PreallocationTools: DiffCache
+import FillArrays
 
 export @derivatives
 
@@ -262,6 +265,7 @@ export IntervalNonlinearProblem
 export OptimizationProblem, constraints
 export SteadyStateProblem
 export JumpProblem
+export SemilinearODEFunction, SemilinearODEProblem
 export alias_elimination, flatten
 export connect, domain_connect, @connector, Connection, AnalysisPoint, Flow, Stream,
        instream