Merge branch 'master' into document-fully-determined

.github/workflows/Downstream.yml

@@ -32,7 +32,7 @@ jobs:
           - {user: SciML, repo: NeuralPDE.jl, group: NNPDE}
           - {user: SciML, repo: DataDrivenDiffEq.jl, group: Downstream}
           - {user: SciML, repo: StructuralIdentifiability.jl, group: All}
-          - {user: SciML, repo: ModelingToolkitStandardLibrary.jl}
+          - {user: SciML, repo: ModelingToolkitStandardLibrary.jl, group: Core}
           - {user: SciML, repo: ModelOrderReduction.jl, group: All}
           - {user: SciML, repo: MethodOfLines.jl, group: Interface}
           - {user: SciML, repo: MethodOfLines.jl, group: 2D_Diffusion}

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 = "9.8.0"
+version = "9.25.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -32,7 +32,8 @@ Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -63,10 +64,12 @@ MTKDeepDiffsExt = "DeepDiffs"
 [compat]
 AbstractTrees = "0.3, 0.4"
 ArrayInterface = "6, 7"
+BifurcationKit = "0.3"
 Combinatorics = "1"
 Compat = "3.42, 4"
 ConstructionBase = "1"
 DataStructures = "0.17, 0.18"
+DeepDiffs = "1"
 DiffEqBase = "6.103.0"
 DiffEqCallbacks = "2.16, 3"
 DiffRules = "0.1, 1.0"
@@ -89,7 +92,9 @@ Libdl = "1"
 LinearAlgebra = "1"
 MLStyle = "0.4.17"
 NaNMath = "0.3, 1"
-OrdinaryDiffEq = "6.73.0"
+NonlinearSolve = "3.12"
+OrderedCollections = "1"
+OrdinaryDiffEq = "6.82.0"
 PrecompileTools = "1"
 RecursiveArrayTools = "2.3, 3"
 Reexport = "0.2, 1"
@@ -102,9 +107,9 @@ SimpleNonlinearSolve = "0.1.0, 1"
 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.11"
-SymbolicUtils = "1.0"
-Symbolics = "5.26"
+SymbolicIndexingInterface = "0.3.12"
+SymbolicUtils = "2.1"
+Symbolics = "5.32"
 URIs = "1"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
@@ -115,14 +120,17 @@ AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
+DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
@@ -136,4 +144,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg"]
+test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"]

docs/Project.toml

@@ -13,6 +13,7 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
@@ -27,14 +28,15 @@ Distributions = "0.25"
 Documenter = "1"
 DynamicQuantities = "^0.11.2, 0.12"
 ModelingToolkit = "8.33, 9"
-NonlinearSolve = "0.3, 1, 2, 3"
+NonlinearSolve = "3"
 Optim = "1.7"
 Optimization = "3.9"
 OptimizationOptimJL = "0.1"
 OrdinaryDiffEq = "6.31"
 Plots = "1.36"
+SciMLStructures = "1.1"
 StochasticDiffEq = "6"
 SymbolicIndexingInterface = "0.3.1"
-SymbolicUtils = "1"
+SymbolicUtils = "2"
 Symbolics = "5"
 Unitful = "1.12"

docs/pages.jl

@@ -17,7 +17,8 @@ pages = [
         "Basic Examples" => Any["examples/higher_order.md",
             "examples/spring_mass.md",
             "examples/modelingtoolkitize_index_reduction.md",
-            "examples/parsing.md"],
+            "examples/parsing.md",
+            "examples/remake.md"],
         "Advanced Examples" => Any["examples/tearing_parallelism.md",
             "examples/sparse_jacobians.md",
             "examples/perturbation.md"]],
@@ -30,6 +31,7 @@ pages = [
         "basics/MTKModel_Connector.md",
         "basics/Validation.md",
         "basics/DependencyGraphs.md",
+        "basics/Precompilation.md",
         "basics/FAQ.md"],
     "System Types" => Any["systems/ODESystem.md",
         "systems/SDESystem.md",

docs/src/basics/FAQ.md

@@ -186,7 +186,7 @@ p, replace, alias = SciMLStructures.canonicalize(Tunable(), prob.p)
 
 This error can come up after running `structural_simplify` on a system that generates dummy derivatives (i.e. variables with `ˍt`).  For example, here even though all the variables are defined with initial values, the `ODEProblem` generation will throw an error that defaults are missing from the variable map.
 
-```
+```julia
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
@@ -197,13 +197,13 @@ eqs = [x1 + x2 + 1 ~ 0
        2 * D(D(x1)) + D(D(x2)) + D(D(x3)) + D(x4) + 4 ~ 0]
 @named sys = ODESystem(eqs, t)
 sys = structural_simplify(sys)
-prob = ODEProblem(sys, [], (0,1))
+prob = ODEProblem(sys, [], (0, 1))
 ```
 
 We can solve this problem by using the `missing_variable_defaults()` function
 
-```
-prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0,1))
+```julia
+prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 1))
 ```
 
 This function provides 0 for the default values, which is a safe assumption for dummy derivatives of most models.  However, the 2nd argument allows for a different default value or values to be used if needed.
@@ -221,12 +221,26 @@ julia> ModelingToolkit.missing_variable_defaults(sys, [1,2,3])
 Use the `u0_constructor` keyword argument to map an array to the desired
 container type. For example:
 
-```
+```julia
 using ModelingToolkit, StaticArrays
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
-sts = @variables x1(t)=0.0
+sts = @variables x1(t) = 0.0
 eqs = [D(x1) ~ 1.1 * x1]
 @mtkbuild sys = ODESystem(eqs, t)
-prob = ODEProblem{false}(sys, [], (0,1); u0_constructor = x->SVector(x...))
+prob = ODEProblem{false}(sys, [], (0, 1); u0_constructor = x -> SVector(x...))
+```
+
+## Using a custom independent variable
+
+When possible, we recommend `using ModelingToolkit: t_nounits as t, D_nounits as D` as the independent variable and its derivative.
+However, if you want to use your own, you can do so:
+
+```julia
+using ModelingToolkit
+
+@independent_variables x
+D = Differential(x)
+@variables y(x)
+@named sys = ODESystem([D(y) ~ x], x)
 ```

docs/src/basics/MTKModel_Connector.md

@@ -39,6 +39,7 @@ Let's explore these in more detail with the following example:
 
 ```@example mtkmodel-example
 using ModelingToolkit
+using ModelingToolkit: t
 
 @mtkmodel ModelA begin
     @parameters begin
@@ -191,6 +192,7 @@ getdefault(model_c3.model_a.k_array[2])
 
 ```@example mtkmodel-example
 using ModelingToolkit
+using ModelingToolkit: t
 
 @mtkmodel M begin
     @parameters begin
@@ -262,6 +264,7 @@ A simple connector can be defined with syntax similar to following example:
 
 ```@example connector
 using ModelingToolkit
+using ModelingToolkit: t
 
 @connector Pin begin
     v(t) = 0.0, [description = "Voltage"]
@@ -344,6 +347,7 @@ The if-elseif-else statements can be used inside `@equations`, `@parameters`,
 
 ```@example branches-in-components
 using ModelingToolkit
+using ModelingToolkit: t
 
 @mtkmodel C begin end
 
@@ -404,12 +408,13 @@ The conditional parts are reflected in the `structure`. For `BranchOutsideTheBlo
 
 ```julia
 julia> BranchOutsideTheBlock.structure
-Dict{Symbol, Any} with 6 entries:
+Dict{Symbol, Any} with 7 entries:
   :components            => Any[(:if, :flag, Vector{Union{Expr, Symbol}}[[:sys1, :C]], Any[])]
   :kwargs                => Dict{Symbol, Dict}(:flag=>Dict{Symbol, Bool}(:value=>1))
   :structural_parameters => Dict{Symbol, Dict}(:flag=>Dict{Symbol, Bool}(:value=>1))
   :independent_variable  => t
-  :parameters            => Dict{Symbol, Dict{Symbol, Any}}(:a2 => Dict(:type => AbstractArray{Real}, :condition => (:if, :flag, Dict{Symbol, Any}(:kwargs => Dict{Any, Any}(:a1 => Dict{Symbol, Union{Nothing, DataType}}(:value => nothing, :type => Real)), :parameters => Any[Dict{Symbol, Dict{Symbol, Any}}(:a1 => Dict(:type => AbstractArray{Real}))]), Dict{Symbol, Any}(:variables => Any[Dict{Symbol, Dict{Symbol, Any}}()], :kwargs => Dict{Any, Any}(:a2 => Dict{Symbol, Union{Nothing, DataType}}(:value => nothing, :type => Real)), :parameters => Any[Dict{Symbol, Dict{Symbol, Any}}(:a2 => Dict(:type => AbstractArray{Real}))]))), :a1 => Dict(:type => AbstractArray{Real}, :condition => (:if, :flag, Dict{Symbol, Any}(:kwargs => Dict{Any, Any}(:a1 => Dict{Symbol, Union{Nothing, DataType}}(:value => nothing, :type => Real)), :parameters => Any[Dict{Symbol, Dict{Symbol, Any}}(:a1 => Dict(:type => AbstractArray{Real}))]), Dict{Symbol, Any}(:variables => Any[Dict{Symbol, Dict{Symbol, Any}}()], :kwargs => Dict{Any, Any}(:a2 => Dict{Symbol, Union{Nothing, DataType}}(:value => nothing, :type => Real)), :parameters => Any[Dict{Symbol, Dict{Symbol, Any}}(:a2 => Dict(:type => AbstractArray{Real}))]))))
+  :parameters            => Dict{Symbol, Dict{Symbol, Any}}(:a2 => Dict(:type=>AbstractArray{Real}, :condition=>(:if, :flag, Dict{Symbol, Any}(:kwargs=>Dict{Any, Any}(:a1=>Dict{Symbol, Union{Nothing, DataType}}(:value=>nothing, :type=>Real)), :parameters=>Any[Dict{Symbol, Dict{Symbol, Any}}(:a1=>Dict(:type=>AbstractArray{Real}))]), Dict{Symbol, Any}(:variables=>Any[Dict{Symbol, Dict{Symbol, Any}}()], :kwargs=>Dict{Any, Any}(:a2=>Dict{Symbol, Union{Nothing, DataType}}(:value=>nothing, :type=>Real)), :parameters=>Any[Dict{Symbol, Dict{Symbol, Any}}(:a2=>Dict(:type=>AbstractArray{Real}))]))), :a1 => Dict(:type=>AbstractArray{Real}, :condition=>(:if, :flag, Dict{Symbol, Any}(:kwargs=>Dict{Any, Any}(:a1=>Dict{Symbol, Union{Nothing, DataType}}(:value=>nothing, :type=>Real)), :parameters=>Any[Dict{Symbol, Dict{Symbol, Any}}(:a1=>Dict(:type=>AbstractArray{Real}))]), Dict{Symbol, Any}(:variables=>Any[Dict{Symbol, Dict{Symbol, Any}}()], :kwargs=>Dict{Any, Any}(:a2=>Dict{Symbol, Union{Nothing, DataType}}(:value=>nothing, :type=>Real)), :parameters=>Any[Dict{Symbol, Dict{Symbol, Any}}(:a2=>Dict(:type=>AbstractArray{Real}))]))))
+  :defaults              => Dict{Symbol, Any}(:a1=>10)
   :equations             => Any[(:if, :flag, ["a1 ~ 0"], ["a2 ~ 0"])]
 ```
 

docs/src/basics/Precompilation.md

@@ -0,0 +1,117 @@
+# Working with Precompilation and Binary Building
+
+## tl;dr, I just want precompilation to work
+
+The tl;dr is, if you want to make precompilation work then instead of
+
+```julia
+ODEProblem(sys, u0, tspan, p)
+```
+
+use:
+
+```julia
+ODEProblem(sys, u0, tspan, p, eval_module = @__MODULE__, eval_expression = true)
+```
+
+As a full example, here's an example of a module that would precompile effectively:
+
+```julia
+module PrecompilationMWE
+using ModelingToolkit
+
+@variables x(ModelingToolkit.t_nounits)
+@named sys = ODESystem([ModelingToolkit.D_nounits(x) ~ -x + 1], ModelingToolkit.t_nounits)
+prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [],
+    eval_expression = true, eval_module = @__MODULE__)
+
+end
+```
+
+If you use that in your package's code then 99% of the time that's the right answer to get
+precompilation working.
+
+## I'm doing something fancier and need a bit more of an explanation
+
+Oh you dapper soul, time for the bigger explanation. Julia's `eval` function evaluates a
+function into a module at a specified world-age. If you evaluate a function within a function
+and try to call it from within that same function, you will hit a world-age error. This looks like:
+
+```julia
+function worldageerror()
+    f = eval(:((x) -> 2x))
+    f(2)
+end
+```
+
+```
+julia> worldageerror()
+ERROR: MethodError: no method matching (::var"#5#6")(::Int64)
+
+Closest candidates are:
+  (::var"#5#6")(::Any) (method too new to be called from this world context.)
+   @ Main REPL[12]:2
+```
+
+This is done for many reasons, in particular if the code that is called within a function could change
+at any time, then Julia functions could not ever properly optimize because the meaning of any function
+or dispatch could always change and you would lose performance by guarding against that. For a full
+discussion of world-age, see [this paper](https://arxiv.org/abs/2010.07516).
+
+However, this would be greatly inhibiting to standard ModelingToolkit usage because then something as
+simple as building an ODEProblem in a function and then using it would get a world age error:
+
+```julia
+function wouldworldage()
+    prob = ODEProblem(sys, [], (0.0, 1.0))
+    sol = solve(prob)
+end
+```
+
+The reason is because `prob.f` would be constructed via `eval`, and thus `prob.f` could not be called
+in the function, which means that no solve could ever work in the same function that generated the
+problem. That does mean that:
+
+```julia
+function wouldworldage()
+    prob = ODEProblem(sys, [], (0.0, 1.0))
+end
+sol = solve(prob)
+```
+
+is fine, or putting
+
+```julia
+prob = ODEProblem(sys, [], (0.0, 1.0))
+sol = solve(prob)
+```
+
+at the top level of a module is perfectly fine too. They just cannot happen in the same function.
+
+This would be a major limitation to ModelingToolkit, and thus we developed
+[RuntimeGeneratedFunctions](https://github.com/SciML/RuntimeGeneratedFunctions.jl) to get around
+this limitation. It will not be described beyond that, it is dark art and should not be investigated.
+But it does the job. But that does mean that it plays... oddly with Julia's compilation.
+
+There are ways to force RuntimeGeneratedFunctions to perform their evaluation and caching within
+a given module, but that is not recommended because it does not play nicely with Julia v1.9's
+introduction of package images for binary caching.
+
+Thus when trying to make things work with precompilation, we recommend using `eval`. This is
+done by simply adding `eval_expression=true` to the problem constructor. However, this is not
+a silver bullet because the moment you start using eval, all potential world-age restrictions
+apply, and thus it is recommended this is simply used for evaluating at the top level of modules
+for the purpose of precompilation and ensuring binaries of your MTK functions are built correctly.
+
+However, there is one caveat that `eval` in Julia works depending on the module that it is given.
+If you have `MyPackage` that you are precompiling into, or say you are using `juliac` or PackageCompiler
+or some other static ahead-of-time (AOT) Julia compiler, then you don't want to accidentally `eval`
+that function to live in ModelingToolkit and instead want to make sure it is `eval`'d to live in `MyPackage`
+(since otherwise it will not cache into the binary). ModelingToolkit cannot know that in advance, and thus
+you have to pass in the module you wish for the functions to "live" in. This is done via the `eval_module`
+argument.
+
+Hence `ODEProblem(sys, u0, tspan, p, eval_module=@__MODULE__, eval_expression=true)` will work if you
+are running this expression in the scope of the module you wish to be precompiling. However, if you are
+attempting to AOT compile a different module, this means that `eval_module` needs to be appropriately
+chosen. And, because `eval_expression=true`, all caveats of world-age apply.