Merge pull request #1852 from SciML/myb/debugging

Add `debug_system` that gives a more helpful error message

Author: YingboMa

src/ModelingToolkit.jl

@@ -147,6 +147,7 @@ include("systems/dependency_graphs.jl")
 include("systems/systemstructure.jl")
 using .SystemStructures
 
+include("debugging.jl")
 include("systems/alias_elimination.jl")
 include("structural_transformation/StructuralTransformations.jl")
 
@@ -210,5 +211,6 @@ export build_function
 export modelingtoolkitize
 export @variables, @parameters
 export @named, @nonamespace, @namespace, extend, compose
+export debug_system
 
 end # module

src/debugging.jl

@@ -0,0 +1,36 @@
+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 LoggedFun{F}
+    f::F
+    args::Any
+end
+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"))
+    end
+end
+
+function logged_fun(f, args...)
+    # Currently we don't really support complex numbers
+    term(LoggedFun(f, args), args..., type = Real)
+end
+
+debug_sub(eq::Equation) = debug_sub(eq.lhs) ~ debug_sub(eq.rhs)
+function debug_sub(ex)
+    istree(ex) || return ex
+    f = operation(ex)
+    args = map(debug_sub, arguments(ex))
+    f in LOGGED_FUN ? logged_fun(f, args...) :
+    similarterm(ex, f, args, metadata = metadata(ex))
+end

src/systems/abstractsystem.jl

@@ -953,6 +953,40 @@ end
 """
 $(SIGNATURES)
 
+Replace functions with singularities with a function that errors with symbolic
+information. E.g.
+
+```julia-repl
+julia> sys = debug_system(sys);
+
+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)
+  ...
+```
+"""
+function debug_system(sys::AbstractSystem)
+    if has_systems(sys) && !isempty(get_systems(sys))
+        error("debug_system only works on systems with no sub-systems!")
+    end
+    if has_eqs(sys)
+        @set! sys.eqs = debug_sub.(equations(sys))
+    end
+    if has_observed(sys)
+        @set! sys.observed = debug_sub.(observed(sys))
+    end
+    return sys
+end
+
+"""
+$(SIGNATURES)
+
 Structurally simplify algebraic equations in a system and compute the
 topological sort of the observed equations. When `simplify=true`, the `simplify`
 function will be applied during the tearing process. It also takes kwargs

test/odesystem.jl

@@ -889,3 +889,20 @@ eqs = [D(q) ~ -p / L - F
 testdict = Dict([:name => "test"])
 @named sys = ODESystem(eqs, t, metadata = testdict)
 @test get_metadata(sys) == testdict
+
+@variables t P(t)=0 Q(t)=2
+∂t = Differential(t)
+
+eqs = [∂t(Q) ~ 1 / sin(P)
+       ∂t(P) ~ log(-cos(Q))]
+@named sys = ODESystem(eqs, t, [P, Q], [])
+sys = debug_system(sys);
+prob = ODEProblem(sys, [], (0, 1.0));
+du = zero(prob.u0);
+if VERSION < v"1.8"
+    @test_throws DomainError prob.f(du, [1, 0], prob.p, 0.0)
+    @test_throws DomainError prob.f(du, [0, 2], prob.p, 0.0)
+else
+    @test_throws "-cos(Q(t))" prob.f(du, [1, 0], prob.p, 0.0)
+    @test_throws "sin(P(t))" prob.f(du, [0, 2], prob.p, 0.0)
+end