From b1cfde7a5489c17419e20ac12729a49653fc89f5 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Fri, 10 Jul 2026 11:03:45 -0400
Subject: [PATCH 01/12] symbolic logging code

---
 lib/ModelingToolkitBase/Project.toml          |   2 +
 .../src/ModelingToolkitBase.jl                |   4 +
 lib/ModelingToolkitBase/src/debugging.jl      | 243 ++++++++++++++++++
 3 files changed, 249 insertions(+)

diff --git a/lib/ModelingToolkitBase/Project.toml b/lib/ModelingToolkitBase/Project.toml
index 711d78bbd2..a182e824e1 100644
--- a/lib/ModelingToolkitBase/Project.toml
+++ b/lib/ModelingToolkitBase/Project.toml
@@ -40,6 +40,7 @@ OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReadOnlyDicts = "795d4caa-f5a7-4580-b5d8-c01d53451803"
@@ -160,6 +161,7 @@ OrdinaryDiffEqTsit5 = "1, 2"
 Pkg = "1"
 PreallocationTools = "0.4.27, 1"
 PrecompileTools = "1.2.1"
+Printf = "1"
 Pyomo = "0.1.0"
 REPL = "1"
 Random = "1"
diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
index 34ff9ebf23..428a64214a 100644
--- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
+++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
@@ -21,6 +21,10 @@ using PrecompileTools, Reexport
     import BandedMatrices: BandedMatrices, BandedMatrix, bandwidths
 end
 
+import SciMLBase
+import SciMLBase: diagnose_symbolic_instability
+using Printf: @sprintf
+
 import SymbolicUtils
 import SymbolicUtils as SU
 import SymbolicUtils: iscall, arguments, operation, maketerm, promote_symtype,
diff --git a/lib/ModelingToolkitBase/src/debugging.jl b/lib/ModelingToolkitBase/src/debugging.jl
index 82e4c80627..6ad45f08e7 100644
--- a/lib/ModelingToolkitBase/src/debugging.jl
+++ b/lib/ModelingToolkitBase/src/debugging.jl
@@ -99,3 +99,246 @@ function get_assertions_expr(sys::AbstractSystem)
     end
     return term
 end
+
+function SciMLBase.diagnose_symbolic_instability(integrator::SciMLBase.DEIntegrator)
+    sys = integrator.f.sys
+    u = integrator.u
+    uprev = integrator.uprev
+    diagnosis = String[]
+
+    #check for assertion failures
+    unks = unknowns(sys)
+    curr_substitution_map = Dict(zip(unks, u))
+    prev_substitution_map = Dict(zip(unknowns(sys), uprev))
+
+    for (cond, msg) in assertions(sys)
+        subclauses = String[]
+        find_failing_subterms(cond, prev_substitution_map, curr_substitution_map, subclauses)
+        if !isempty(subclauses)
+            push!(diagnosis, "\n\nAssertion violated: $cond - \"$msg\"")
+            append!(diagnosis, subclauses)
+        end
+    end
+
+    #find singularity causes in equations
+    singularities = String[]
+    for eq in equations(sys)
+        find_singular_subterms(eq, eq.rhs, prev_substitution_map, singularities)
+    end
+    if !isempty(singularities)
+        push!(diagnosis, "\nSymbolic Analysis of MTK System:")
+        append!(diagnosis, singularities)
+    end
+
+    return isempty(diagnosis) ? "" : join(diagnosis, "\n")
+end
+
+function find_singular_subterms(eq, expr, sub_map, diagnosis)
+    expr = Symbolics.unwrap(expr)
+    !SymbolicUtils.iscall(expr) && return diagnosis
+    op = SymbolicUtils.operation(expr)
+    args = SymbolicUtils.arguments(expr)
+
+    if op === (/) #division, singular if we divide by small thing
+        d = Symbolics.value(Symbolics.substitute(args[2], sub_map))
+        if d isa Number && abs(d) < 1e-10
+            push!(diagnosis, "in equation $eq: division by very small value $(args[2]) ≈ $(@sprintf("%.4g", d)) leads to singularity.")
+        end
+    elseif op === log #singular if we log small thing
+        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        if x isa Number && x <= 1e-10
+            push!(diagnosis, "in equation $eq: log of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
+        end
+    elseif op === sqrt 
+        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        if x isa Number && x < 1e-10
+            push!(diagnosis, "in equation $eq: sqrt of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
+        end
+    elseif op === (^)
+        e = Symbolics.value(Symbolics.substitute(args[2], sub_map))
+        b = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        if e isa Number && b isa Number #two cases
+            if e < 0 && abs(b) < 1e-10
+                push!(diagnosis, "in equation $eq: ($(args[1])) raised to power $e with base ≈ $(@sprintf("%.4g", b)) going to 0; result diverges.")
+            elseif e > 0 && abs(b) > 1
+                push!(diagnosis, "in equation $eq: ($(args[1]) ≈ $(@sprintf("%.4g", b))) raised to power $e - base magnitude is large and being amplified.")
+            end
+        end
+    end
+
+    for arg in args
+        find_singular_subterms(eq, arg, sub_map, diagnosis)
+    end
+    return diagnosis
+end
+
+function find_failing_subterms(cond, prev_map, curr_map, diagnosis)
+    c = Symbolics.unwrap(cond)
+    !SymbolicUtils.iscall(c) && return diagnosis
+    op = SymbolicUtils.operation(c)
+    args = SymbolicUtils.arguments(c)
+
+    if (op === (<) || op === (>) || op === (<=) || op === (>=)) && length(args) == 2
+        #compare using previous non-nan values to find violating subclauses, then output current values
+        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
+        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
+        if lhs isa Number && rhs isa Number
+            # small margin -> violated
+            margin = (op === (<) || op === (<=)) ? rhs - lhs : lhs - rhs
+            if margin <= 1e-6
+                push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
+            end
+        end
+    elseif op === (!=) && length(args) == 2
+        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
+        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
+        if lhs isa Number && rhs isa Number && abs(lhs - rhs) <= 1e-6
+            push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
+        end
+    elseif op === (==) && length(args) == 2
+        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
+        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
+        if lhs isa Number && rhs isa Number && abs(lhs - rhs) > 1e-6
+            push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
+        end
+    else #recurse
+        for arg in args
+            find_failing_subterms(arg, prev_map, curr_map, diagnosis)
+        end
+    end
+    return diagnosis
+end
+
+function clause_values(c, curr_map)
+    parts = String[]
+    for v in Symbolics.get_variables(c)
+        val = Symbolics.value(Symbolics.substitute(v, curr_map))
+        push!(parts, val isa Number ? "$v = $(@sprintf("%.4g", val))" : "$v = $val")
+    end
+    return join(parts, ", ")
+end
+
+#=
+
+module OrdinaryDiffEqModelingToolkitExt
+
+using OrdinaryDiffEqCore, ModelingToolkit
+using Printf: @sprintf
+
+function OrdinaryDiffEqCore.system_singularity_rootcause(sys, u, uprev)
+    diagnosis = String[]
+
+    #check for assertion failures
+    unks = unknowns(sys)
+    curr_substitution_map = Dict(zip(unks, u))
+    prev_substitution_map = Dict(zip(unknowns(sys), uprev))
+
+    for (cond, msg) in ModelingToolkit.assertions(sys)
+        subclauses = String[]
+        find_failing_subterms(cond, prev_substitution_map, curr_substitution_map, subclauses)
+        if !isempty(subclauses)
+            push!(diagnosis, "\n\nAssertion violated: $cond - \"$msg\"")
+            append!(diagnosis, subclauses)
+        end
+    end
+
+    #find singularity causes in equations
+    singularities = String[]
+    for eq in equations(sys)
+        find_singular_subterms(eq, eq.rhs, prev_substitution_map, singularities)
+    end
+    if !isempty(singularities)
+        push!(diagnosis, "\nSymbolic Analysis of MTK System:")
+        append!(diagnosis, singularities)
+    end
+
+    return diagnosis
+end
+
+function find_singular_subterms(eq, expr, sub_map, diagnosis)
+    expr = Symbolics.unwrap(expr)
+    !SymbolicUtils.iscall(expr) && return diagnosis
+    op = SymbolicUtils.operation(expr)
+    args = SymbolicUtils.arguments(expr)
+
+    if op === (/) #division, singular if we divide by small thing
+        d = Symbolics.value(Symbolics.substitute(args[2], sub_map))
+        if d isa Number && abs(d) < 1e-10
+            push!(diagnosis, "in equation $eq: division by very small value $(args[2]) ≈ $(@sprintf("%.4g", d)) leads to singularity.")
+        end
+    elseif op === log #singular if we log small thing
+        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        if x isa Number && x <= 1e-10
+            push!(diagnosis, "in equation $eq: log of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
+        end
+    elseif op === sqrt 
+        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        if x isa Number && x < 1e-10
+            push!(diagnosis, "in equation $eq: sqrt of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
+        end
+    elseif op === (^)
+        e = Symbolics.value(Symbolics.substitute(args[2], sub_map))
+        b = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        if e isa Number && b isa Number #two cases
+            if e < 0 && abs(b) < 1e-10
+                push!(diagnosis, "in equation $eq: ($(args[1])) raised to power $e with base ≈ $(@sprintf("%.4g", b)) going to 0; result diverges.")
+            elseif e > 0 && abs(b) > 1
+                push!(diagnosis, "in equation $eq: ($(args[1]) ≈ $(@sprintf("%.4g", b))) raised to power $e - base magnitude is large and being amplified.")
+            end
+        end
+    end
+
+    for arg in args
+        find_singular_subterms(eq, arg, sub_map, diagnosis)
+    end
+    return diagnosis
+end
+
+function find_failing_subterms(cond, prev_map, curr_map, diagnosis)
+    c = Symbolics.unwrap(cond)
+    !SymbolicUtils.iscall(c) && return diagnosis
+    op = SymbolicUtils.operation(c)
+    args = SymbolicUtils.arguments(c)
+
+    if (op === (<) || op === (>) || op === (<=) || op === (>=)) && length(args) == 2
+        #compare using previous non-nan values to find violating subclauses, then output current values
+        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
+        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
+        if lhs isa Number && rhs isa Number
+            # small margin -> violated
+            margin = (op === (<) || op === (<=)) ? rhs - lhs : lhs - rhs
+            if margin <= 1e-6
+                push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
+            end
+        end
+    elseif op === (!=) && length(args) == 2
+        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
+        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
+        if lhs isa Number && rhs isa Number && abs(lhs - rhs) <= 1e-6
+            push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
+        end
+    elseif op === (==) && length(args) == 2
+        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
+        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
+        if lhs isa Number && rhs isa Number && abs(lhs - rhs) > 1e-6
+            push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
+        end
+    else #recurse
+        for arg in args
+            find_failing_subterms(arg, prev_map, curr_map, diagnosis)
+        end
+    end
+    return diagnosis
+end
+
+function clause_values(c, curr_map)
+    parts = String[]
+    for v in Symbolics.get_variables(c)
+        val = Symbolics.value(Symbolics.substitute(v, curr_map))
+        push!(parts, val isa Number ? "$v = $(@sprintf("%.4g", val))" : "$v = $val")
+    end
+    return join(parts, ", ")
+end
+
+end
+ =#
\ No newline at end of file

From 6c46654b17596368fc082b8178c3598da8d37874 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Fri, 10 Jul 2026 11:08:27 -0400
Subject: [PATCH 02/12] remove dead code

---
 lib/ModelingToolkitBase/src/debugging.jl | 125 -----------------------
 1 file changed, 125 deletions(-)

diff --git a/lib/ModelingToolkitBase/src/debugging.jl b/lib/ModelingToolkitBase/src/debugging.jl
index 6ad45f08e7..ab7b4b103d 100644
--- a/lib/ModelingToolkitBase/src/debugging.jl
+++ b/lib/ModelingToolkitBase/src/debugging.jl
@@ -217,128 +217,3 @@ function clause_values(c, curr_map)
     end
     return join(parts, ", ")
 end
-
-#=
-
-module OrdinaryDiffEqModelingToolkitExt
-
-using OrdinaryDiffEqCore, ModelingToolkit
-using Printf: @sprintf
-
-function OrdinaryDiffEqCore.system_singularity_rootcause(sys, u, uprev)
-    diagnosis = String[]
-
-    #check for assertion failures
-    unks = unknowns(sys)
-    curr_substitution_map = Dict(zip(unks, u))
-    prev_substitution_map = Dict(zip(unknowns(sys), uprev))
-
-    for (cond, msg) in ModelingToolkit.assertions(sys)
-        subclauses = String[]
-        find_failing_subterms(cond, prev_substitution_map, curr_substitution_map, subclauses)
-        if !isempty(subclauses)
-            push!(diagnosis, "\n\nAssertion violated: $cond - \"$msg\"")
-            append!(diagnosis, subclauses)
-        end
-    end
-
-    #find singularity causes in equations
-    singularities = String[]
-    for eq in equations(sys)
-        find_singular_subterms(eq, eq.rhs, prev_substitution_map, singularities)
-    end
-    if !isempty(singularities)
-        push!(diagnosis, "\nSymbolic Analysis of MTK System:")
-        append!(diagnosis, singularities)
-    end
-
-    return diagnosis
-end
-
-function find_singular_subterms(eq, expr, sub_map, diagnosis)
-    expr = Symbolics.unwrap(expr)
-    !SymbolicUtils.iscall(expr) && return diagnosis
-    op = SymbolicUtils.operation(expr)
-    args = SymbolicUtils.arguments(expr)
-
-    if op === (/) #division, singular if we divide by small thing
-        d = Symbolics.value(Symbolics.substitute(args[2], sub_map))
-        if d isa Number && abs(d) < 1e-10
-            push!(diagnosis, "in equation $eq: division by very small value $(args[2]) ≈ $(@sprintf("%.4g", d)) leads to singularity.")
-        end
-    elseif op === log #singular if we log small thing
-        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
-        if x isa Number && x <= 1e-10
-            push!(diagnosis, "in equation $eq: log of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
-        end
-    elseif op === sqrt 
-        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
-        if x isa Number && x < 1e-10
-            push!(diagnosis, "in equation $eq: sqrt of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
-        end
-    elseif op === (^)
-        e = Symbolics.value(Symbolics.substitute(args[2], sub_map))
-        b = Symbolics.value(Symbolics.substitute(args[1], sub_map))
-        if e isa Number && b isa Number #two cases
-            if e < 0 && abs(b) < 1e-10
-                push!(diagnosis, "in equation $eq: ($(args[1])) raised to power $e with base ≈ $(@sprintf("%.4g", b)) going to 0; result diverges.")
-            elseif e > 0 && abs(b) > 1
-                push!(diagnosis, "in equation $eq: ($(args[1]) ≈ $(@sprintf("%.4g", b))) raised to power $e - base magnitude is large and being amplified.")
-            end
-        end
-    end
-
-    for arg in args
-        find_singular_subterms(eq, arg, sub_map, diagnosis)
-    end
-    return diagnosis
-end
-
-function find_failing_subterms(cond, prev_map, curr_map, diagnosis)
-    c = Symbolics.unwrap(cond)
-    !SymbolicUtils.iscall(c) && return diagnosis
-    op = SymbolicUtils.operation(c)
-    args = SymbolicUtils.arguments(c)
-
-    if (op === (<) || op === (>) || op === (<=) || op === (>=)) && length(args) == 2
-        #compare using previous non-nan values to find violating subclauses, then output current values
-        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
-        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
-        if lhs isa Number && rhs isa Number
-            # small margin -> violated
-            margin = (op === (<) || op === (<=)) ? rhs - lhs : lhs - rhs
-            if margin <= 1e-6
-                push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
-            end
-        end
-    elseif op === (!=) && length(args) == 2
-        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
-        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
-        if lhs isa Number && rhs isa Number && abs(lhs - rhs) <= 1e-6
-            push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
-        end
-    elseif op === (==) && length(args) == 2
-        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
-        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
-        if lhs isa Number && rhs isa Number && abs(lhs - rhs) > 1e-6
-            push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
-        end
-    else #recurse
-        for arg in args
-            find_failing_subterms(arg, prev_map, curr_map, diagnosis)
-        end
-    end
-    return diagnosis
-end
-
-function clause_values(c, curr_map)
-    parts = String[]
-    for v in Symbolics.get_variables(c)
-        val = Symbolics.value(Symbolics.substitute(v, curr_map))
-        push!(parts, val isa Number ? "$v = $(@sprintf("%.4g", val))" : "$v = $val")
-    end
-    return join(parts, ", ")
-end
-
-end
- =#
\ No newline at end of file

From f8344c22d1ab7083e5012437abf5075f93f1300f Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <142109039+Shreyas-Ekanathan@users.noreply.github.com>
Date: Mon, 13 Jul 2026 11:17:22 -0400
Subject: [PATCH 03/12] Update lib/ModelingToolkitBase/src/debugging.jl

Co-authored-by: Aayush Sabharwal <aayush.sabharwal@gmail.com>
---
 lib/ModelingToolkitBase/src/debugging.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/ModelingToolkitBase/src/debugging.jl b/lib/ModelingToolkitBase/src/debugging.jl
index ab7b4b103d..d371841c4f 100644
--- a/lib/ModelingToolkitBase/src/debugging.jl
+++ b/lib/ModelingToolkitBase/src/debugging.jl
@@ -140,7 +140,7 @@ function find_singular_subterms(eq, expr, sub_map, diagnosis)
     args = SymbolicUtils.arguments(expr)
 
     if op === (/) #division, singular if we divide by small thing
-        d = Symbolics.value(Symbolics.substitute(args[2], sub_map))
+        d = Symbolics.value(SymbolicUtils.substitute(args[2], sub_map))
         if d isa Number && abs(d) < 1e-10
             push!(diagnosis, "in equation $eq: division by very small value $(args[2]) ≈ $(@sprintf("%.4g", d)) leads to singularity.")
         end

From 89b86b51a2c6d53de3d520085b1a5444b2aeb0b2 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <142109039+Shreyas-Ekanathan@users.noreply.github.com>
Date: Mon, 13 Jul 2026 11:17:29 -0400
Subject: [PATCH 04/12] Update lib/ModelingToolkitBase/src/debugging.jl

Co-authored-by: Aayush Sabharwal <aayush.sabharwal@gmail.com>
---
 lib/ModelingToolkitBase/src/debugging.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/ModelingToolkitBase/src/debugging.jl b/lib/ModelingToolkitBase/src/debugging.jl
index d371841c4f..8ddc30e3ff 100644
--- a/lib/ModelingToolkitBase/src/debugging.jl
+++ b/lib/ModelingToolkitBase/src/debugging.jl
@@ -134,7 +134,7 @@ function SciMLBase.diagnose_symbolic_instability(integrator::SciMLBase.DEIntegra
 end
 
 function find_singular_subterms(eq, expr, sub_map, diagnosis)
-    expr = Symbolics.unwrap(expr)
+    expr = unwrap(expr)
     !SymbolicUtils.iscall(expr) && return diagnosis
     op = SymbolicUtils.operation(expr)
     args = SymbolicUtils.arguments(expr)

From bb9f9f0a82504b8de70e0c90fec1ad7d0a6f972a Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Tue, 14 Jul 2026 10:18:39 -0400
Subject: [PATCH 05/12] some optimizations

---
 lib/ModelingToolkitBase/src/debugging.jl | 31 ++++++++++++------------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/lib/ModelingToolkitBase/src/debugging.jl b/lib/ModelingToolkitBase/src/debugging.jl
index 8ddc30e3ff..a5087e803e 100644
--- a/lib/ModelingToolkitBase/src/debugging.jl
+++ b/lib/ModelingToolkitBase/src/debugging.jl
@@ -100,16 +100,13 @@ function get_assertions_expr(sys::AbstractSystem)
     return term
 end
 
-function SciMLBase.diagnose_symbolic_instability(integrator::SciMLBase.DEIntegrator)
-    sys = integrator.f.sys
-    u = integrator.u
-    uprev = integrator.uprev
+function SciMLBase.diagnose_symbolic_instability(sys::AbstractSystem, u, uprev)
     diagnosis = String[]
 
     #check for assertion failures
     unks = unknowns(sys)
-    curr_substitution_map = Dict(zip(unks, u))
-    prev_substitution_map = Dict(zip(unknowns(sys), uprev))
+    curr_substitution_map = Dict(zip(unks, Symbolics.unwrap.(Symbolics.Num.(u))))
+    prev_substitution_map = Dict(zip(unknowns(sys), Symbolics.unwrap.(Symbolics.Num.(uprev))))
 
     for (cond, msg) in assertions(sys)
         subclauses = String[]
@@ -122,8 +119,10 @@ function SciMLBase.diagnose_symbolic_instability(integrator::SciMLBase.DEIntegra
 
     #find singularity causes in equations
     singularities = String[]
-    for eq in equations(sys)
-        find_singular_subterms(eq, eq.rhs, prev_substitution_map, singularities)
+    visited = IdDict()
+    subber = SymbolicUtils.IRSubstituter{true}(get_irstructure(sys), prev_substitution_map)
+    for eq in full_equations(sys)
+        find_singular_subterms(eq, eq.rhs, subber, singularities, visited)
     end
     if !isempty(singularities)
         push!(diagnosis, "\nSymbolic Analysis of MTK System:")
@@ -133,30 +132,32 @@ function SciMLBase.diagnose_symbolic_instability(integrator::SciMLBase.DEIntegra
     return isempty(diagnosis) ? "" : join(diagnosis, "\n")
 end
 
-function find_singular_subterms(eq, expr, sub_map, diagnosis)
+function find_singular_subterms(eq, expr, sub_map, diagnosis, visited)
     expr = unwrap(expr)
     !SymbolicUtils.iscall(expr) && return diagnosis
     op = SymbolicUtils.operation(expr)
     args = SymbolicUtils.arguments(expr)
+    haskey(visited, expr) && return diagnosis  
+    visited[expr] = nothing
 
     if op === (/) #division, singular if we divide by small thing
-        d = Symbolics.value(SymbolicUtils.substitute(args[2], sub_map))
+        d = Symbolics.value(sub_map(args[2]))
         if d isa Number && abs(d) < 1e-10
             push!(diagnosis, "in equation $eq: division by very small value $(args[2]) ≈ $(@sprintf("%.4g", d)) leads to singularity.")
         end
     elseif op === log #singular if we log small thing
-        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        x = Symbolics.value(sub_map(args[1]))
         if x isa Number && x <= 1e-10
             push!(diagnosis, "in equation $eq: log of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
         end
     elseif op === sqrt 
-        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        x = Symbolics.value(sub_map(args[1]))
         if x isa Number && x < 1e-10
             push!(diagnosis, "in equation $eq: sqrt of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
         end
     elseif op === (^)
-        e = Symbolics.value(Symbolics.substitute(args[2], sub_map))
-        b = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        e = Symbolics.value(sub_map(args[2]))
+        b = Symbolics.value(sub_map(args[1]))
         if e isa Number && b isa Number #two cases
             if e < 0 && abs(b) < 1e-10
                 push!(diagnosis, "in equation $eq: ($(args[1])) raised to power $e with base ≈ $(@sprintf("%.4g", b)) going to 0; result diverges.")
@@ -167,7 +168,7 @@ function find_singular_subterms(eq, expr, sub_map, diagnosis)
     end
 
     for arg in args
-        find_singular_subterms(eq, arg, sub_map, diagnosis)
+        find_singular_subterms(eq, arg, sub_map, diagnosis, visited)
     end
     return diagnosis
 end

From bef8677a71442a54ace286a14508f6b102294d7e Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 15 Jul 2026 10:39:19 +0530
Subject: [PATCH 06/12] Apply suggestions from code review

Co-authored-by: Aayush Sabharwal <aayush.sabharwal@gmail.com>
---
 lib/ModelingToolkitBase/src/debugging.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/lib/ModelingToolkitBase/src/debugging.jl b/lib/ModelingToolkitBase/src/debugging.jl
index a5087e803e..2dce5bc3b3 100644
--- a/lib/ModelingToolkitBase/src/debugging.jl
+++ b/lib/ModelingToolkitBase/src/debugging.jl
@@ -105,8 +105,8 @@ function SciMLBase.diagnose_symbolic_instability(sys::AbstractSystem, u, uprev)
 
     #check for assertion failures
     unks = unknowns(sys)
-    curr_substitution_map = Dict(zip(unks, Symbolics.unwrap.(Symbolics.Num.(u))))
-    prev_substitution_map = Dict(zip(unknowns(sys), Symbolics.unwrap.(Symbolics.Num.(uprev))))
+    curr_substitution_map = Dict(zip(unks, u))
+    prev_substitution_map = Dict(zip(unknowns(sys), uprev))
 
     for (cond, msg) in assertions(sys)
         subclauses = String[]
@@ -119,7 +119,7 @@ function SciMLBase.diagnose_symbolic_instability(sys::AbstractSystem, u, uprev)
 
     #find singularity causes in equations
     singularities = String[]
-    visited = IdDict()
+    visited = IdDict{SymbolicT, Nothing}()
     subber = SymbolicUtils.IRSubstituter{true}(get_irstructure(sys), prev_substitution_map)
     for eq in full_equations(sys)
         find_singular_subterms(eq, eq.rhs, subber, singularities, visited)

From 147487be8a37bf800be31a5e1eed8617c626b45a Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 15 Jul 2026 09:41:35 -0400
Subject: [PATCH 07/12] bug fix

---
 lib/ModelingToolkitBase/src/debugging.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lib/ModelingToolkitBase/src/debugging.jl b/lib/ModelingToolkitBase/src/debugging.jl
index 2dce5bc3b3..f7dd9a7cd4 100644
--- a/lib/ModelingToolkitBase/src/debugging.jl
+++ b/lib/ModelingToolkitBase/src/debugging.jl
@@ -105,8 +105,8 @@ function SciMLBase.diagnose_symbolic_instability(sys::AbstractSystem, u, uprev)
 
     #check for assertion failures
     unks = unknowns(sys)
-    curr_substitution_map = Dict(zip(unks, u))
-    prev_substitution_map = Dict(zip(unknowns(sys), uprev))
+    curr_substitution_map = Dict(zip(unks, Symbolics.unwrap.(Symbolics.Num.(u))))
+    prev_substitution_map = Dict(zip(unknowns(sys), Symbolics.unwrap.(Symbolics.Num.(uprev))))
 
     for (cond, msg) in assertions(sys)
         subclauses = String[]

From 6c994ef923f888f8cd12610a4d6b2d4165f82e3a Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 15 Jul 2026 10:53:13 -0400
Subject: [PATCH 08/12] type dictionary concretely

---
 lib/ModelingToolkitBase/src/debugging.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lib/ModelingToolkitBase/src/debugging.jl b/lib/ModelingToolkitBase/src/debugging.jl
index f7dd9a7cd4..bbc410f7a4 100644
--- a/lib/ModelingToolkitBase/src/debugging.jl
+++ b/lib/ModelingToolkitBase/src/debugging.jl
@@ -105,8 +105,8 @@ function SciMLBase.diagnose_symbolic_instability(sys::AbstractSystem, u, uprev)
 
     #check for assertion failures
     unks = unknowns(sys)
-    curr_substitution_map = Dict(zip(unks, Symbolics.unwrap.(Symbolics.Num.(u))))
-    prev_substitution_map = Dict(zip(unknowns(sys), Symbolics.unwrap.(Symbolics.Num.(uprev))))
+    curr_substitution_map = Dict{SymbolicT, SymbolicT}(zip(unks, u))
+    prev_substitution_map = Dict{SymbolicT, SymbolicT}(zip(unknowns(sys), uprev))
 
     for (cond, msg) in assertions(sys)
         subclauses = String[]

From f12f936ee04525ecab2fd63d4b6f6034e3ca2967 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 29 Jul 2026 09:37:08 -0400
Subject: [PATCH 09/12] compat bump

---
 .DS_Store                            | Bin 0 -> 8196 bytes
 lib/ModelingToolkitBase/Project.toml |   2 +-
 2 files changed, 1 insertion(+), 1 deletion(-)
 create mode 100644 .DS_Store

diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..8188a722d69989703696071490a97679e1188bc9
GIT binary patch
literal 8196
zcmeHMYitx%6u#fIz#Y24v=k_8v0GNC!6kj9AO*yw+Y))UbX&T8K)bs$!pL-{?96V_
zf~h7N4dUaM#8<*E<BLBk(WsF`CB(-cgBbXsi9ZY`niw_s$3)M}o#oNe#2>~e&Q0z;
z_ndprnK@^^bMM}}j4`z4)kemu7-O8OOO-NeZcyBDw@xXNFDWGo(r5N)#;`KPy^|T8
zGCL3i5eOm>L?DPj5P={9|Aq+Ap6wPs&Au;2gEojj5P|<?1l;u@Le-^lpH2wsEgjT&
zM*yPk5x^_dSKJQKh)?4_oe-2ps6eO^6{_%$7!c|tkNETAKAjL$s5AJ75BN`pe?o!x
zbh01urZdC^4cZ_AK?J5oz)hbLW-y&)GGo`C-woY%8XC?+Dk+^Ya~7Y)%lMw;erqi0
zq@BFon$8WlT(@NznPNYk*{y41LrP`BGM%(;YB^6qS2db4*sGhWHP&k7P1W{vlVX6&
zTpm&ukB>JsuZ=~T8rL3)MaCPOTUN&+Yg$^49Fh5gX#M)Gp##Q<X&n_Fz+VN}oJk%%
ze$^erE+oWtL`f<$D%`nPbVr_D7Mc@I_4N-7hK7`i0=;*~sH0n^y<5xKM2D1lHu<+D
zYnh!n%^mN}T1I}vw1<@PY}&}?jkKd}GK{2tzeciob~a}jMyF-#u9?1+qaAYGu_;5T
ze%s193C~zcC;MWSlDD-kjg0CAt&L3Em6Ke9d~)uhnp<yM)pW;(O`Esw>^e73DX$1s
zh8Fe~bX(6DTF;oSI{W+5c2+Z0-Q3qt0X7}|9$m}H<t0buZZkciovT{Bq#|4$zGXTs
zu0(k%pV1HIb>}|ujM{KiR^H-+D!tbZOLZ-(-blDMT*vu<D^KJ!YK1toG8~hYH1D=4
z0K!Cls8LpIx%&X2K&Uo{)^J6YdyI5OGlZxmw4N*5<lch#212+o6zAoL$*HV0GD7~?
zLYxS<bLA}G>ppvX+IBidHM96uiSQQA`|bVdoYpZ$Ff3~9;N1JJBzINXRo$_2A|SQl
zU7Yv1#x|)c*{N8(q$=4R4i71cZWeR`J$+SyedP<zd9?U-qVbiY8$GQKdp$4A(RhOp
z%e1jAWlgM&ZD+$Q$Hv()cATAHZ?ccs8TKXnntjK9V3z?%n1PunLnW%Q3^A-iJsQx0
zRwU4l%}8Ppd$1Sx!iIwa#&8%9;}JZH$M6K6!Ero`=WqhA;8nbilQ@MBa2g-tBYckY
z_zD+r5#Qry{DNQchr}gWnk`jJ%cT`kozy6;lXgfQQm3>_>XHT|RT`0U(ztLzDrVbW
z!53TUF!lV!gR;%P{FC0eakkamy=CjR?f;!@yD-<yulbb=7A}e`Us>P0_RiKT5EP&-
z;`Cbldyu~_Q9Q_dn3#-J%EWN>(ojtg<%D)<a`qPUQ<>$WpSv&(*DedyC^IN)1iPY+
z*qM|ug1tRfQ=>?fsDf>1tO+S4l!mgWyE;ZJe5PexD5R89LW;4CcM;3nHA3)-HnQZV
z*WU*83+!8ViT%WWC!o*8d{m*1z`X_=2-rK(N#O2A4}p6a8Dybi1Qre=M^HbE2QYyL
z@eq#U7@ovac$$FzJYK+ycnL2P*k2>KpTgUC7w_SHoWWUqf=}_82m5bu2|wcUBp~lB
zLRe7zb`{8lWX?1!b00~Hu)lut&BK)lSE}6efBW?J|2L1rfhR!(f(T5B07~1EZLK7p
zZozZU+7YV5RNdmvn-J7fq2{`a0OGIwVMzT5O?7qRKAjMhBvk(KivaI8`Cs_89qj+X
L{_i35=4}21%IIt#

literal 0
HcmV?d00001

diff --git a/lib/ModelingToolkitBase/Project.toml b/lib/ModelingToolkitBase/Project.toml
index a182e824e1..b04c7ff5fa 100644
--- a/lib/ModelingToolkitBase/Project.toml
+++ b/lib/ModelingToolkitBase/Project.toml
@@ -172,7 +172,7 @@ ReferenceTests = "0.10"
 RuntimeGeneratedFunctions = "0.5.12"
 SCCNonlinearSolve = "1.13"
 SafeTestsets = "0.1"
-SciMLBase = "3.18"
+SciMLBase = "3.40.1"
 SciMLPublic = "1.0.0"
 SciMLStructures = "1.7"
 Serialization = "1"

From bc4b9d7d30af4290c83179bf55dbc10337cc8fa3 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <142109039+Shreyas-Ekanathan@users.noreply.github.com>
Date: Wed, 29 Jul 2026 09:40:03 -0400
Subject: [PATCH 10/12] Delete .DS_Store

---
 .DS_Store | Bin 8196 -> 0 bytes
 1 file changed, 0 insertions(+), 0 deletions(-)
 delete mode 100644 .DS_Store

diff --git a/.DS_Store b/.DS_Store
deleted file mode 100644
index 8188a722d69989703696071490a97679e1188bc9..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 8196
zcmeHMYitx%6u#fIz#Y24v=k_8v0GNC!6kj9AO*yw+Y))UbX&T8K)bs$!pL-{?96V_
zf~h7N4dUaM#8<*E<BLBk(WsF`CB(-cgBbXsi9ZY`niw_s$3)M}o#oNe#2>~e&Q0z;
z_ndprnK@^^bMM}}j4`z4)kemu7-O8OOO-NeZcyBDw@xXNFDWGo(r5N)#;`KPy^|T8
zGCL3i5eOm>L?DPj5P={9|Aq+Ap6wPs&Au;2gEojj5P|<?1l;u@Le-^lpH2wsEgjT&
zM*yPk5x^_dSKJQKh)?4_oe-2ps6eO^6{_%$7!c|tkNETAKAjL$s5AJ75BN`pe?o!x
zbh01urZdC^4cZ_AK?J5oz)hbLW-y&)GGo`C-woY%8XC?+Dk+^Ya~7Y)%lMw;erqi0
zq@BFon$8WlT(@NznPNYk*{y41LrP`BGM%(;YB^6qS2db4*sGhWHP&k7P1W{vlVX6&
zTpm&ukB>JsuZ=~T8rL3)MaCPOTUN&+Yg$^49Fh5gX#M)Gp##Q<X&n_Fz+VN}oJk%%
ze$^erE+oWtL`f<$D%`nPbVr_D7Mc@I_4N-7hK7`i0=;*~sH0n^y<5xKM2D1lHu<+D
zYnh!n%^mN}T1I}vw1<@PY}&}?jkKd}GK{2tzeciob~a}jMyF-#u9?1+qaAYGu_;5T
ze%s193C~zcC;MWSlDD-kjg0CAt&L3Em6Ke9d~)uhnp<yM)pW;(O`Esw>^e73DX$1s
zh8Fe~bX(6DTF;oSI{W+5c2+Z0-Q3qt0X7}|9$m}H<t0buZZkciovT{Bq#|4$zGXTs
zu0(k%pV1HIb>}|ujM{KiR^H-+D!tbZOLZ-(-blDMT*vu<D^KJ!YK1toG8~hYH1D=4
z0K!Cls8LpIx%&X2K&Uo{)^J6YdyI5OGlZxmw4N*5<lch#212+o6zAoL$*HV0GD7~?
zLYxS<bLA}G>ppvX+IBidHM96uiSQQA`|bVdoYpZ$Ff3~9;N1JJBzINXRo$_2A|SQl
zU7Yv1#x|)c*{N8(q$=4R4i71cZWeR`J$+SyedP<zd9?U-qVbiY8$GQKdp$4A(RhOp
z%e1jAWlgM&ZD+$Q$Hv()cATAHZ?ccs8TKXnntjK9V3z?%n1PunLnW%Q3^A-iJsQx0
zRwU4l%}8Ppd$1Sx!iIwa#&8%9;}JZH$M6K6!Ero`=WqhA;8nbilQ@MBa2g-tBYckY
z_zD+r5#Qry{DNQchr}gWnk`jJ%cT`kozy6;lXgfQQm3>_>XHT|RT`0U(ztLzDrVbW
z!53TUF!lV!gR;%P{FC0eakkamy=CjR?f;!@yD-<yulbb=7A}e`Us>P0_RiKT5EP&-
z;`Cbldyu~_Q9Q_dn3#-J%EWN>(ojtg<%D)<a`qPUQ<>$WpSv&(*DedyC^IN)1iPY+
z*qM|ug1tRfQ=>?fsDf>1tO+S4l!mgWyE;ZJe5PexD5R89LW;4CcM;3nHA3)-HnQZV
z*WU*83+!8ViT%WWC!o*8d{m*1z`X_=2-rK(N#O2A4}p6a8Dybi1Qre=M^HbE2QYyL
z@eq#U7@ovac$$FzJYK+ycnL2P*k2>KpTgUC7w_SHoWWUqf=}_82m5bu2|wcUBp~lB
zLRe7zb`{8lWX?1!b00~Hu)lut&BK)lSE}6efBW?J|2L1rfhR!(f(T5B07~1EZLK7p
zZozZU+7YV5RNdmvn-J7fq2{`a0OGIwVMzT5O?7qRKAjMhBvk(KivaI8`Cs_89qj+X
L{_i35=4}21%IIt#


From f10b5c49915d307c2a6ff66fb210561829786a60 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 29 Jul 2026 10:26:12 -0400
Subject: [PATCH 11/12] fix compat edit

---
 lib/ModelingToolkitBase/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/ModelingToolkitBase/Project.toml b/lib/ModelingToolkitBase/Project.toml
index 5ed0f70a9b..915c9e207b 100644
--- a/lib/ModelingToolkitBase/Project.toml
+++ b/lib/ModelingToolkitBase/Project.toml
@@ -174,7 +174,7 @@ ReferenceTests = "0.10"
 RuntimeGeneratedFunctions = "0.5.12"
 SCCNonlinearSolve = "1.13"
 SafeTestsets = "0.1"
-SciMLBase = "3.40.1"
+SciMLBase = "3.38"
 SciMLPublic = "1.0.0"
 SciMLStructures = "1.7"
 Serialization = "1"

From 4a4ca8677258ab2027fb0a14252a874f8c8c9611 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 29 Jul 2026 15:59:33 -0400
Subject: [PATCH 12/12] fix few failing tests

---
 docs/Project.toml                                      | 2 +-
 lib/ModelingToolkitBase/test/optimization/Project.toml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 39d5ddb615..df14347ec6 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -45,7 +45,7 @@ CairoMakie = "0.15"
 CommonSolve = "0.2"
 ControlSystemsBase = "1.20"
 ControlSystemsMTK = "2.6"
-DataInterpolations = "8"
+DataInterpolations = "8, 9"
 Distributions = "0.25"
 Documenter = "1"
 DynamicQuantities = "1"
diff --git a/lib/ModelingToolkitBase/test/optimization/Project.toml b/lib/ModelingToolkitBase/test/optimization/Project.toml
index 7c89fb85f8..fa13aadfdc 100644
--- a/lib/ModelingToolkitBase/test/optimization/Project.toml
+++ b/lib/ModelingToolkitBase/test/optimization/Project.toml
@@ -32,7 +32,7 @@ ModelingToolkitBase = {path = "../.."}
 
 [compat]
 CasADi = "1.0.7"
-DataInterpolations = "8.8"
+DataInterpolations = "8.8, 9"
 OrdinaryDiffEqExplicitTableaus = "2"
 OrdinaryDiffEqImplicitTableaus = "2"
 SafeTestsets = "0.1, 1"