Merge branch 'master' into sde_throws

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 = "9.19.0"
+version = "9.23.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -108,8 +108,8 @@ 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.12"
-SymbolicUtils = "2"
-Symbolics = "5.30.1"
+SymbolicUtils = "2.1"
+Symbolics = "5.32"
 URIs = "1"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"

docs/pages.jl

@@ -31,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/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.

src/discretedomain.jl

@@ -4,6 +4,8 @@ struct SampleTime <: Operator
     SampleTime() = SymbolicUtils.term(SampleTime, type = Real)
 end
 SymbolicUtils.promote_symtype(::Type{<:SampleTime}, t...) = Real
+Base.nameof(::SampleTime) = :SampleTime
+SymbolicUtils.isbinop(::SampleTime) = false
 
 # Shift
 
@@ -32,6 +34,9 @@ struct Shift <: Operator
 end
 Shift(steps::Int) = new(nothing, steps)
 normalize_to_differential(s::Shift) = Differential(s.t)^s.steps
+Base.nameof(::Shift) = :Shift
+SymbolicUtils.isbinop(::Shift) = false
+
 function (D::Shift)(x, allow_zero = false)
     !allow_zero && D.steps == 0 && return x
     Term{symtype(x)}(D, Any[x])
@@ -108,6 +113,8 @@ Sample(x) = Sample()(x)
 (D::Sample)(x) = Term{symtype(x)}(D, Any[x])
 (D::Sample)(x::Num) = Num(D(value(x)))
 SymbolicUtils.promote_symtype(::Sample, x) = x
+Base.nameof(::Sample) = :Sample
+SymbolicUtils.isbinop(::Sample) = false
 
 Base.show(io::IO, D::Sample) = print(io, "Sample(", D.clock, ")")
 
@@ -137,6 +144,8 @@ end
 (D::Hold)(x) = Term{symtype(x)}(D, Any[x])
 (D::Hold)(x::Num) = Num(D(value(x)))
 SymbolicUtils.promote_symtype(::Hold, x) = x
+Base.nameof(::Hold) = :Hold
+SymbolicUtils.isbinop(::Hold) = false
 
 Hold(x) = Hold()(x)
 

src/inputoutput.jl

@@ -195,6 +195,8 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
         disturbance_inputs = disturbances(sys);
         implicit_dae = false,
         simplify = false,
+        eval_expression = false,
+        eval_module = @__MODULE__,
         kwargs...)
     isempty(inputs) && @warn("No unbound inputs were found in system.")
 
@@ -240,7 +242,8 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
     end
     process = get_postprocess_fbody(sys)
     f = build_function(rhss, args...; postprocess_fbody = process,
-        expression = Val{false}, kwargs...)
+        expression = Val{true}, kwargs...)
+    f = eval_or_rgf.(f; eval_expression, eval_module)
     (; f, dvs, ps, io_sys = sys)
 end
 
@@ -395,7 +398,7 @@ model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.i
 
 `f_oop` will have an extra state corresponding to the integrator in the disturbance model. This state will not be affected by any input, but will affect the dynamics from where it enters, in this case it will affect additively from `model.torque.tau.u`.
 """
-function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing)
+function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing; kwargs...)
     t = get_iv(sys)
     @variables d(t)=0 [disturbance = true]
     @variables u(t)=0 [input = true] # New system input
@@ -418,6 +421,6 @@ function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing)
     augmented_sys = extend(augmented_sys, sys)
 
     (f_oop, f_ip), dvs, p = generate_control_function(augmented_sys, all_inputs,
-        [d])
+        [d]; kwargs...)
     (f_oop, f_ip), augmented_sys, dvs, p
 end