Merge branch 'master' into test-SDEFunctionExpr

Author: ChrisRackauckas

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"
@@ -143,7 +143,7 @@ StaticArrays = "0.10, 0.11, 0.12, 1.0"
 StochasticDiffEq = "6.72.1"
 StochasticDelayDiffEq = "1.8.1"
 SymbolicIndexingInterface = "0.3.36"
-SymbolicUtils = "3.7"
+SymbolicUtils = "3.10"
 Symbolics = "6.22.1"
 URIs = "1"
 UnPack = "0.1, 1.0"

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/ModelingToolkit.jl

@@ -54,7 +54,8 @@ import SCCNonlinearSolve
 using Reexport
 using RecursiveArrayTools
 import Graphs: SimpleDiGraph, add_edge!, incidence_matrix
-import BlockArrays: BlockedArray, Block, blocksize, blocksizes
+import BlockArrays: BlockArray, BlockedArray, Block, blocksize, blocksizes, blockpush!,
+                    undef_blocks, blocks
 import CommonSolve
 import EnumX
 
@@ -145,8 +146,8 @@ include("systems/index_cache.jl")
 include("systems/parameter_buffer.jl")
 include("systems/abstractsystem.jl")
 include("systems/model_parsing.jl")
-include("systems/analysis_points.jl")
 include("systems/connectors.jl")
+include("systems/analysis_points.jl")
 include("systems/imperative_affect.jl")
 include("systems/callbacks.jl")
 include("systems/problem_utils.jl")

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

@@ -1875,6 +1875,13 @@ Equivalent to `length(equations(expand_connections(sys))) - length(filter(eq ->
 function n_expanded_connection_equations(sys::AbstractSystem)
     # TODO: what about inputs?
     isconnector(sys) && return length(get_unknowns(sys))
+    sys = remove_analysis_points(sys)
+    n_variable_connect_eqs = 0
+    for eq in equations(sys)
+        is_causal_variable_connection(eq.rhs) || continue
+        n_variable_connect_eqs += length(get_systems(eq.rhs)) - 1
+    end
+
     sys, (csets, _) = generate_connection_set(sys)
     ceqs, instream_csets = generate_connection_equations_and_stream_connections(csets)
     n_outer_stream_variables = 0
@@ -1897,7 +1904,7 @@ function n_expanded_connection_equations(sys::AbstractSystem)
     #    n_toplevel_unused_flows += count(x->get_connection_type(x) === Flow && !(x in toplevel_flows), get_unknowns(m))
     #end
 
-    nextras = n_outer_stream_variables + length(ceqs)
+    nextras = n_outer_stream_variables + length(ceqs) + n_variable_connect_eqs
 end
 
 function Base.show(
@@ -1919,7 +1926,7 @@ function Base.show(
     nrows > 0 && hint && print(io, " see hierarchy($name)")
     for i in 1:nrows
         sub = subs[i]
-        name = String(nameof(sub))
+        local name = String(nameof(sub))
         print(io, "\n  ", name)
         desc = description(sub)
         if !isempty(desc)
@@ -2274,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 = debug_system(complete(sys))
 
-julia> sys = complete(sys);
+julia> prob = ODEProblem(sys, [0.0, 2.0], (0.0, 1.0))
 
-julia> prob = ODEProblem(sys, [], (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/analysis_points.jl

@@ -208,6 +208,14 @@ function Symbolics.connect(in::AbstractSystem, name::Symbol, out, outs...; verbo
     return AnalysisPoint() ~ AnalysisPoint(in, name, [out; collect(outs)]; verbose)
 end
 
+function Symbolics.connect(
+        in::ConnectableSymbolicT, name::Symbol, out::ConnectableSymbolicT,
+        outs::ConnectableSymbolicT...; verbose = true)
+    allvars = (in, out, outs...)
+    validate_causal_variables_connection(allvars)
+    return AnalysisPoint() ~ AnalysisPoint(in, name, [out; collect(outs)]; verbose)
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -240,7 +248,7 @@ connection. This is the variable named `u` if present, and otherwise the only
 variable in the system. If the system does not have a variable named `u` and
 contains multiple variables, throw an error.
 """
-function ap_var(sys)
+function ap_var(sys::AbstractSystem)
     if hasproperty(sys, :u)
         return sys.u
     end
@@ -249,6 +257,15 @@ function ap_var(sys)
     error("Could not determine the analysis-point variable in system $(nameof(sys)). To use an analysis point, apply it to a connection between causal blocks which have a variable named `u` or a single unknown of the same size.")
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+For an `AnalysisPoint` involving causal variables. Simply return the variable.
+"""
+function ap_var(var::ConnectableSymbolicT)
+    return var
+end
+
 """
     $(TYPEDEF)