Merge remote-tracking branch 'origin' into validate-inputs-odesys

merge

Author: vyudu

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.60.0"
+version = "9.61.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -123,7 +123,7 @@ NonlinearSolve = "4.3"
 OffsetArrays = "1"
 OrderedCollections = "1"
 OrdinaryDiffEq = "6.82.0"
-OrdinaryDiffEqCore = "1.13.0"
+OrdinaryDiffEqCore = "1.15.0"
 OrdinaryDiffEqDefault = "1.2"
 OrdinaryDiffEqNonlinearSolve = "1.3.0"
 PrecompileTools = "1"
@@ -132,7 +132,7 @@ RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
 SCCNonlinearSolve = "1.0.0"
-SciMLBase = "2.68.1"
+SciMLBase = "2.71.1"
 SciMLStructures = "1.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"

docs/pages.jl

@@ -32,6 +32,7 @@ pages = [
         "basics/InputOutput.md",
         "basics/MTKLanguage.md",
         "basics/Validation.md",
+        "basics/Debugging.md",
         "basics/DependencyGraphs.md",
         "basics/Precompilation.md",
         "basics/FAQ.md"],

docs/src/basics/Debugging.md

@@ -0,0 +1,40 @@
+# Debugging
+
+Every (mortal) modeler writes models that contain mistakes or are susceptible to numerical errors in their hunt for the perfect model.
+Debugging such errors is part of the modeling process, and ModelingToolkit includes some functionality that helps with this.
+
+For example, consider an ODE model with "dangerous" functions (here `√`):
+
+```@example debug
+using ModelingToolkit, OrdinaryDiffEq
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+@variables u1(t) u2(t) u3(t)
+eqs = [D(u1) ~ -√(u1), D(u2) ~ -√(u2), D(u3) ~ -√(u3)]
+defaults = [u1 => 1.0, u2 => 2.0, u3 => 3.0]
+@named sys = ODESystem(eqs, t; defaults)
+sys = structural_simplify(sys)
+```
+
+This problem causes the ODE solver to crash:
+
+```@repl debug
+prob = ODEProblem(sys, [], (0.0, 10.0), []);
+sol = solve(prob, Tsit5());
+```
+
+This suggests *that* something went wrong, but not exactly *what* went wrong and *where* it did.
+In such situations, the `debug_system` function is helpful:
+
+```@repl debug
+dsys = debug_system(sys; functions = [sqrt]);
+dprob = ODEProblem(dsys, [], (0.0, 10.0), []);
+dsol = solve(dprob, Tsit5());
+```
+
+Now we see that it crashed because `u1` decreased so much that it became negative and outside the domain of the `√` function.
+We could have figured that out ourselves, but it is not always so obvious for more complex models.
+
+```@docs
+debug_system
+```

src/debugging.jl

@@ -1,36 +1,44 @@
-const LOGGED_FUN = Set([log, sqrt, (^), /, inv])
-is_legal(::typeof(/), a, b) = is_legal(inv, b)
-is_legal(::typeof(inv), a) = !iszero(a)
-is_legal(::Union{typeof(log), typeof(sqrt)}, a) = a isa Complex || a >= zero(a)
-is_legal(::typeof(^), a, b) = a isa Complex || b isa Complex || isinteger(b) || a >= zero(a)
-
+struct LoggedFunctionException <: Exception
+    msg::String
+end
 struct LoggedFun{F}
     f::F
     args::Any
+    error_nonfinite::Bool
+end
+function LoggedFunctionException(lf::LoggedFun, args, msg)
+    LoggedFunctionException(
+        "Function $(lf.f)($(join(lf.args, ", "))) " * msg * " with input" *
+        join("\n  " .* string.(lf.args .=> args)) # one line for each "var => val" for readability
+    )
 end
+Base.showerror(io::IO, err::LoggedFunctionException) = print(io, err.msg)
 Base.nameof(lf::LoggedFun) = nameof(lf.f)
 SymbolicUtils.promote_symtype(::LoggedFun, Ts...) = Real
 function (lf::LoggedFun)(args...)
-    f = lf.f
-    symbolic_args = lf.args
-    if is_legal(f, args...)
-        f(args...)
-    else
-        args_str = join(string.(symbolic_args .=> args), ", ", ", and ")
-        throw(DomainError(args, "$(lf.f) errors with input(s): $args_str"))
+    val = try
+        lf.f(args...) # try to call with numerical input, as usual
+    catch err
+        throw(LoggedFunctionException(lf, args, "errors")) # Julia automatically attaches original error message
     end
+    if lf.error_nonfinite && !isfinite(val)
+        throw(LoggedFunctionException(lf, args, "output non-finite value $val"))
+    end
+    return val
 end
 
-function logged_fun(f, args...)
+function logged_fun(f, args...; error_nonfinite = true) # remember to update error_nonfinite in debug_system() docstring
     # Currently we don't really support complex numbers
-    term(LoggedFun(f, args), args..., type = Real)
+    term(LoggedFun(f, args, error_nonfinite), args..., type = Real)
 end
 
-debug_sub(eq::Equation) = debug_sub(eq.lhs) ~ debug_sub(eq.rhs)
-function debug_sub(ex)
+function debug_sub(eq::Equation, funcs; kw...)
+    debug_sub(eq.lhs, funcs; kw...) ~ debug_sub(eq.rhs, funcs; kw...)
+end
+function debug_sub(ex, funcs; kw...)
     iscall(ex) || return ex
     f = operation(ex)
-    args = map(debug_sub, arguments(ex))
-    f in LOGGED_FUN ? logged_fun(f, args...) :
+    args = map(ex -> debug_sub(ex, funcs; kw...), arguments(ex))
+    f in funcs ? logged_fun(f, args...; kw...) :
     maketerm(typeof(ex), f, args, metadata(ex))
 end

src/systems/abstractsystem.jl

@@ -2281,37 +2281,37 @@ macro mtkbuild(exprs...)
 end
 
 """
-$(SIGNATURES)
+    debug_system(sys::AbstractSystem; functions = [log, sqrt, (^), /, inv, asin, acos], error_nonfinite = true)
 
-Replace functions with singularities with a function that errors with symbolic
-information. E.g.
+Wrap `functions` in `sys` so any error thrown in them shows helpful symbolic-numeric
+information about its input. If `error_nonfinite`, functions that output nonfinite
+values (like `Inf` or `NaN`) also display errors, even though the raw function itself
+does not throw an exception (like `1/0`). For example:
 
 ```julia-repl
-julia> sys = debug_system(sys);
-
-julia> sys = complete(sys);
+julia> sys = debug_system(complete(sys))
 
-julia> prob = ODEProblem(sys, [], (0, 1.0));
+julia> prob = ODEProblem(sys, [0.0, 2.0], (0.0, 1.0))
 
-julia> du = zero(prob.u0);
-
-julia> prob.f(du, prob.u0, prob.p, 0.0)
-ERROR: DomainError with (-1.0,):
-log errors with input(s): -cos(Q(t)) => -1.0
-Stacktrace:
-  [1] (::ModelingToolkit.LoggedFun{typeof(log)})(args::Float64)
-  ...
+julia> prob.f(prob.u0, prob.p, 0.0)
+ERROR: Function /(1, sin(P(t))) output non-finite value Inf with input
+  1 => 1
+  sin(P(t)) => 0.0
 ```
 """
-function debug_system(sys::AbstractSystem)
+function debug_system(
+        sys::AbstractSystem; functions = [log, sqrt, (^), /, inv, asin, acos], kw...)
+    if !(functions isa Set)
+        functions = Set(functions) # more efficient "in" lookup
+    end
     if has_systems(sys) && !isempty(get_systems(sys))
-        error("debug_system only works on systems with no sub-systems!")
+        error("debug_system(sys) only works on systems with no sub-systems! Consider flattening it with flatten(sys) or structural_simplify(sys) first.")
     end
     if has_eqs(sys)
-        @set! sys.eqs = debug_sub.(equations(sys))
+        @set! sys.eqs = debug_sub.(equations(sys), Ref(functions); kw...)
     end
     if has_observed(sys)
-        @set! sys.observed = debug_sub.(observed(sys))
+        @set! sys.observed = debug_sub.(observed(sys), Ref(functions); kw...)
     end
     return sys
 end

src/systems/diffeqs/abstractodesystem.jl

@@ -1335,7 +1335,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem,
 
     # TODO: throw on uninitialized arrays
     filter!(x -> !(x isa Symbolics.Arr), uninit)
-    if !isempty(uninit)
+    if is_time_dependent(sys) && !isempty(uninit)
         allow_incomplete || throw(IncompleteInitializationError(uninit))
         # for incomplete initialization, we will add the missing variables as parameters.
         # they will be updated by `update_initializeprob!` and `initializeprobmap` will
@@ -1373,6 +1373,21 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem,
     end
 
     u0map = merge(ModelingToolkit.guesses(sys), todict(guesses), todict(u0map))
+
+    # Replace dummy derivatives in u0map: D(x) -> x_t etc.
+    if has_schedule(sys)
+        schedule = get_schedule(sys)
+        if !isnothing(schedule)
+            for (var, val) in u0map
+                dvar = get(schedule.dummy_sub, var, var) # with dummy derivatives
+                if dvar !== var # then replace it
+                    delete!(u0map, var)
+                    push!(u0map, dvar => val)
+                end
+            end
+        end
+    end
+
     fullmap = merge(u0map, parammap)
     u0T = Union{}
     for sym in unknowns(isys)

src/systems/diffeqs/odesystem.jl

@@ -374,6 +374,8 @@ function Base.:(==)(sys1::ODESystem, sys2::ODESystem)
         _eq_unordered(get_eqs(sys1), get_eqs(sys2)) &&
         _eq_unordered(get_unknowns(sys1), get_unknowns(sys2)) &&
         _eq_unordered(get_ps(sys1), get_ps(sys2)) &&
+        _eq_unordered(continuous_events(sys1), continuous_events(sys2)) &&
+        _eq_unordered(discrete_events(sys1), discrete_events(sys2)) && 
         all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2)))
 end
 
@@ -630,7 +632,15 @@ function build_explicit_observed_function(sys, ts;
     oop_fn = Func(args, [],
                  pre(Let(obsexprs,
                      return_value,
-                     false))) |> array_wrapper[1] |> oop_mtkp_wrapper |> toexpr
+                     false)), [Expr(:meta, :propagate_inbounds)]) |> array_wrapper[1] |> oop_mtkp_wrapper |> toexpr
+
+    if !checkbounds
+        oop_fn.args[end] = quote
+            @inbounds begin
+                $(oop_fn.args[end])
+            end
+        end
+    end
     oop_fn = expression ? oop_fn : eval_or_rgf(oop_fn; eval_expression, eval_module)
 
     if !isscalar
@@ -640,6 +650,18 @@ function build_explicit_observed_function(sys, ts;
             wrap_code = mtkparams_wrapper .∘ array_wrapper .∘
                         wrap_assignments(isscalar, obsexprs),
             expression = Val{true})[2]
+        if !checkbounds
+            iip_fn.args[end] = quote
+                @inbounds begin
+                    $(iip_fn.args[end])
+                end
+            end
+        end
+        iip_fn.args[end] = quote
+            $(Expr(:meta, :propagate_inbounds))
+            $(iip_fn.args[end])
+        end
+
         if !expression
             iip_fn = eval_or_rgf(iip_fn; eval_expression, eval_module)
         end

src/systems/diffeqs/sdesystem.jl

@@ -280,11 +280,13 @@ function Base.:(==)(sys1::SDESystem, sys2::SDESystem)
     iv2 = get_iv(sys2)
     isequal(iv1, iv2) &&
         isequal(nameof(sys1), nameof(sys2)) &&
-        isequal(get_eqs(sys1), get_eqs(sys2)) &&
-        isequal(get_noiseeqs(sys1), get_noiseeqs(sys2)) &&
+        _eq_unordered(get_eqs(sys1), get_eqs(sys2)) &&
+        _eq_unordered(get_noiseeqs(sys1), get_noiseeqs(sys2)) &&
         isequal(get_is_scalar_noise(sys1), get_is_scalar_noise(sys2)) &&
         _eq_unordered(get_unknowns(sys1), get_unknowns(sys2)) &&
         _eq_unordered(get_ps(sys1), get_ps(sys2)) &&
+        _eq_unordered(continuous_events(sys1), continuous_events(sys2)) &&
+        _eq_unordered(discrete_events(sys1), discrete_events(sys2)) && 
         all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2)))
 end
 
@@ -738,6 +740,7 @@ function DiffEqBase.SDEProblem{iip, specialize}(
     if !iscomplete(sys)
         error("A completed `SDESystem` is required. Call `complete` or `structural_simplify` on the system before creating an `SDEProblem`")
     end
+
     f, u0, p = process_SciMLProblem(
         SDEFunction{iip, specialize}, sys, u0map, parammap; check_length,
         t = tspan === nothing ? nothing : tspan[1], kwargs...)
@@ -768,6 +771,15 @@ function DiffEqBase.SDEProblem{iip, specialize}(
         noise_rate_prototype = noise_rate_prototype, kwargs...)
 end
 
+function DiffEqBase.SDEProblem(sys::ODESystem, args...; kwargs...)
+
+    if any(ModelingToolkit.isbrownian, unknowns(sys))
+        error("SDESystem constructed by defining Brownian variables with @brownian must be simplified by calling `structural_simplify` before a SDEProblem can be constructed.")
+    else
+        error("Cannot construct SDEProblem from a normal ODESystem.")
+    end
+end
+
 """
 ```julia
 DiffEqBase.SDEProblem{iip}(sys::SDESystem, u0map, tspan, p = parammap;

src/systems/jumps/jumpsystem.jl

@@ -426,14 +426,15 @@ function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple,
         error("The passed in JumpSystem contains `Equation`s or continuous events, please use a problem type that supports these features, such as ODEProblem.")
     end
 
-    _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap;
+    _f, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap;
         t = tspan === nothing ? nothing : tspan[1], use_union, tofloat = false, check_length = false)
     f = DiffEqBase.DISCRETE_INPLACE_DEFAULT
 
     observedfun = ObservedFunctionCache(
         sys; eval_expression, eval_module, checkbounds = get(kwargs, :checkbounds, false))
 
-    df = DiscreteFunction{true, true}(f; sys = sys, observed = observedfun)
+    df = DiscreteFunction{true, true}(f; sys = sys, observed = observedfun,
+        initialization_data = get(_f.kwargs, :initialization_data, nothing))
     DiscreteProblem(df, u0, tspan, p; kwargs...)
 end