up

Author: TorkelE

src/Catalyst.jl

@@ -38,9 +38,6 @@ import ModelingToolkit: check_variables,
                         check_parameters, _iszero, _merge, check_units,
                         get_unit, check_equations, iscomplete
 
-# Event-related MTK stuff.
-import ModelingToolkit: PresetTimeCallback
-
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
 import MacroTools, Graphs
 import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne

src/reaction_network.jl

@@ -831,7 +831,7 @@ function read_equations_options(options, variables_declared)
     # Loops through all equations, checks for lhs of the form `D(X) ~ ...`.
     # When this is the case, the variable X and differential D are extracted (for automatic declaration).
     # Also performs simple error checks.
-    vars_extracted = []
+    vars_extracted = Vector{Symbol}()
     add_default_diff = false
     for eq in equations
         if (eq.head != :call) || (eq.args[1] != :~)

src/reactionsystem.jl

@@ -1596,7 +1596,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
                       combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
                       include_zero_odes = true, remove_conserved = false, checks = false,
                       default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
-                      kwargs...)
+                      all_differentials_permitted = false, kwargs...)
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, NonlinearSystem)
     fullrs = Catalyst.flatten(rs)
@@ -1608,7 +1608,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
 
     # Throws a warning if there are differential equations in non-standard format.
     # Next, sets all differential terms to `0`.
-    nonlinear_convert_differentials_check(rs)
+    all_differentials_permitted || nonlinear_convert_differentials_check(rs)
     eqs = [remove_diffs(eq.lhs) ~ remove_diffs(eq.rhs) for eq in eqs]
 
 
@@ -1635,24 +1635,26 @@ diff_2_zero(expr) = (Symbolics.is_derivative(expr) ? 0.0 : expr)
 # does not contain any differentials. If this is not teh case, we throw a warning to let the user
 # know that they should be careful.
 function nonlinear_convert_differentials_check(rs::ReactionSystem)
-    for eq in equations(rs)
-        # For each equation (that is not a reaction), checks in order: 
+    for eq in filter(is_diff_equation, equations(rs))
+        # For each differential equation, checks in order: 
         # If there is a differential on the right hand side.
         # If the lhs is not on the form D(...).
         # If the lhs upper level function is not a differential w.r.t. time.
         # If the contenct of the differential is not a variable (and nothing more).
         # If either of this is a case, throws the warning.
-        (eq isa Reaction) && continue
         if Symbolics._occursin(Symbolics.is_derivative, eq.rhs) ||
                 !Symbolics.istree(eq.lhs) ||
                 !isequal(Symbolics.operation(eq.lhs), Differential(get_iv(rs))) || 
                 (length(arguments(eq.lhs)) != 1) ||
                 !any(isequal(arguments(eq.lhs)[1]), nonspecies(rs))
-            @warn "You are attempting to convert a `ReactionSystem` coupled with differential equations. However, some of these differentials are not of the form `D(x) ~ ...` where: 
+            error("You are attempting to convert a `ReactionSystem` coupled with differential equations to a `NonlinearSystem`. However, some of these differentials are not of the form `D(x) ~ ...` where: 
                     (1) The left-hand side is a differential of a single variable with respect to the time independent variable, and
                     (2) The right-hand side does not contain any differentials.
-                  This is permitted (and all differential will be set to `0`), however, it is recommended to proceed with caution to ensure that the produced nonlinear equation is the desired one."
-            return
+                This is generally not permitted.
+                
+                If you still would like to perform this conversions, please use the `all_differentials_permitted = true` option. In this case, all differential will be set to `0`. 
+                However, it is recommended to proceed with caution to ensure that the produced nonlinear equation makes sense for you intended application."
+                )
         end
     end
 end
@@ -1790,13 +1792,14 @@ function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0,
                                      name = nameof(rs), include_zero_odes = true,
                                      combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
                                      remove_conserved = false, checks = false,
-                                     check_length = false, kwargs...)
+                                     check_length = false, all_differentials_permitted = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     nlsys = convert(NonlinearSystem, rs; name, combinatoric_ratelaws, include_zero_odes,
-                    checks, remove_conserved)
+                    checks, all_differentials_permitted, remove_conserved)
     nlsys = complete(nlsys)
-    return NonlinearProblem(nlsys, u0map, pmap, args...; check_length, kwargs...)
+    return NonlinearProblem(nlsys, u0map, pmap, args...; check_length, 
+                            kwargs...)
 end
 
 # SDEProblem from AbstractReactionNetwork

test/reactionsystem_structure/coupled_equation_reaction_systems.jl

@@ -437,48 +437,49 @@ end
 
 ### Coupled NonlinearSystems Tests ###
 
-# Checks that systems with weird differential equations yield appropriate warnings.
+# Checks that systems with weird differential equations yield errors.
 let 
-    # This oen is normal, and should not yield a warning.
+    # This one is normal, and should not yield an error.
     begin
         rs = @reaction_network begin
             @equations D(V) ~ 1.0 - V
         end
         @test_nowarn convert(NonlinearSystem, rs)
     end
 
-    # Higher-order differential on the lhs, should yield a warning.
+    # Higher-order differential on the lhs, should yield an error.
     begin
         rs = @reaction_network begin
             @differentials D = Differential(t)
             @variables V(t)
             @equations D(D(V)) ~ 1.0 - V
             (p,d), 0 <--> X
         end
-        @test_logs (:warn, r".*") convert(NonlinearSystem, rs)
+        @test_throws Exception convert(NonlinearSystem, rs)
     end
 
-    # Differential on the rhs, should yield a warning.
+    # Differential on the rhs, should yield an error.
     begin
         rs = @reaction_network begin
             @variables U(t)
             @equations D(V) ~ 1.0 - V + D(U)
             (p,d), 0 <--> X
         end
-        @test_logs (:warn, r".*") convert(NonlinearSystem, rs)
+        @test_throws Exception convert(NonlinearSystem, rs)
     end
 
-    # Non-differential term on the lhs, should yield a warning.
+    # Non-differential term on the lhs, should yield an error.
     begin
         rs = @reaction_network begin
             @differentials D = Differential(t)
             @variables V(t)
             @equations D(V) + V ~ 1.0 - V
             (p,d), 0 <--> X
         end
-        @test_logs (:warn, r".*") convert(NonlinearSystem, rs)
+        @test_throws Exception convert(NonlinearSystem, rs)
     end
 end
+
 ### Unusual Differentials Tests ###
 
 # Tests that coupled CRN/DAEs with higher order differentials can be created.
@@ -522,7 +523,7 @@ let
     # Here `B => 0.1` has to be provided as well (and it shouldn't for the 2nd order ODE), hence the 
     # separate `u0` declaration.
     u0 = [X => 1.0, A => 2.0, D(A) => 1.0, B => 0.1]
-    nlprob = NonlinearProblem(coupled_rs, u0, ps; structural_simplify = true)
+    nlprob = NonlinearProblem(coupled_rs, u0, ps; structural_simplify = true, all_differentials_permitted = true)
     nlsol = solve(nlprob)
     @test nlsol[X][end] ≈ 2.0
     @test nlsol[A][end] ≈ 0.0