Disallow intermediate equations in DiffEqSystem

Author: HarrisonGrodin

src/systems/diffeqs/diffeqsystem.jl

@@ -4,6 +4,12 @@ mutable struct DiffEqSystem <: AbstractSystem
     dvs::Vector{Variable}
     ps::Vector{Variable}
     jac::Matrix{Expression}
+    function DiffEqSystem(eqs, ivs, dvs, ps, jac)
+        all(!isintermediate, eqs) ||
+            throw(ArgumentError("no intermediate equations permitted in DiffEqSystem"))
+
+        new(eqs, ivs, dvs, ps, jac)
+    end
 end
 
 DiffEqSystem(eqs, ivs, dvs, ps) = DiffEqSystem(eqs, ivs, dvs, ps, Matrix{Expression}(undef,0,0))
@@ -20,6 +26,9 @@ function DiffEqSystem(eqs, ivs)
     DiffEqSystem(eqs, ivs, dvs, ps, Matrix{Expression}(undef,0,0))
 end
 
+isintermediate(eq::Equation) = !(isa(eq.lhs, Operation) && isa(eq.lhs.op, Differential))
+
+
 function generate_ode_function(sys::DiffEqSystem;version = ArrayFunction)
     var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in eachindex(sys.dvs)]
     param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in eachindex(sys.ps)]
@@ -42,8 +51,6 @@ function generate_ode_function(sys::DiffEqSystem;version = ArrayFunction)
     end
 end
 
-isintermediate(eq::Equation) = !(isa(eq.lhs, Operation) && isa(eq.lhs.op, Differential))
-
 function build_equals_expr(eq::Equation)
     @assert !isintermediate(eq)
 
@@ -52,13 +59,7 @@ function build_equals_expr(eq::Equation)
 end
 
 function calculate_jacobian(sys::DiffEqSystem, simplify=true)
-    calcs, diff_exprs = partition(isintermediate, sys.eqs)
-    rhs = [eq.rhs for eq in diff_exprs]
-
-    # Handle intermediate calculations by substitution
-    for calc ∈ calcs
-        find_replace!.(rhs, calc.lhs, calc.rhs)
-    end
+    rhs = [eq.rhs for eq in sys.eqs]
 
     sys_exprs = calculate_jacobian(rhs, sys.dvs)
     sys_exprs = Expression[expand_derivatives(expr) for expr in sys_exprs]
@@ -68,7 +69,6 @@ end
 function generate_ode_jacobian(sys::DiffEqSystem, simplify=true)
     var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in eachindex(sys.dvs)]
     param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in eachindex(sys.ps)]
-    diff_exprs = filter(!isintermediate, sys.eqs)
     jac = calculate_jacobian(sys, simplify)
     sys.jac = jac
     jac_exprs = [:(J[$i,$j] = $(convert(Expr, jac[i,j]))) for i in 1:size(jac,1), j in 1:size(jac,2)]
@@ -80,7 +80,6 @@ end
 function generate_ode_iW(sys::DiffEqSystem, simplify=true)
     var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in eachindex(sys.dvs)]
     param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in eachindex(sys.ps)]
-    diff_exprs = filter(!isintermediate, sys.eqs)
     jac = sys.jac
 
     gam = Parameter(:gam)

src/systems/diffeqs/first_order_transform.jl

@@ -19,13 +19,11 @@ function ode_order_lowering(sys::DiffEqSystem; kwargs...)
 end
 ode_order_lowering(eqs; naming_scheme = "_") = ode_order_lowering!(deepcopy(eqs), naming_scheme)
 function ode_order_lowering!(eqs, naming_scheme)
-    ind = findfirst(x->!(isintermediate(x)), eqs)
-    idv = extract_idv(eqs[ind])
+    idv = extract_idv(eqs[1])
     D   = Differential(idv, 1)
     sym_order = Dict{Symbol, Int}()
     dv_name = eqs[1].lhs.subtype
     for eq in eqs
-        isintermediate(eq) && continue
         sym, maxorder = extract_symbol_order(eq)
         maxorder == 1 && continue # fast pass
         if maxorder > get(sym_order, sym, 0)