updates

Author: TorkelE

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -107,7 +107,7 @@ function filter_negative_f(sols; neg_thres=-1e-20)
     return filter(sol -> all(>=(0), sol), sols)
 end
 
-# Sometimes (when polynomials are created from hybrid CRN/DAEs), the steady state polynomial have the wrong type.
+# Sometimes (when polynomials are created from coupled CRN/DAEs), the steady state polynomial have the wrong type.
 # This converts it to the correct type, which homotopy continuation can handle.
 const WRONG_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}}}
 const CORRECT_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}, Float64}}

src/Catalyst.jl

@@ -30,7 +30,8 @@ using ModelingToolkit: Symbolic, value, istree, get_unknowns, get_ps, get_iv, ge
 
 import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_variables!,
                         modified_unknowns!, validate, namespace_variables,
-                        namespace_parameters, rename, renamespace, getname, flatten
+                        namespace_parameters, rename, renamespace, getname, flatten,
+                        is_alg_equation, is_diff_equation
 
 # internal but needed ModelingToolkit functions
 import ModelingToolkit: check_variables,

src/reaction_network.jl

@@ -401,7 +401,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters = vcat(parameters_declared, parameters_extracted)
 
     # Create differential expression.
-    diffexpr = create_differential_expr(options, add_default_diff, [species; parameters; variables])
+    diffexpr = create_differential_expr(options, add_default_diff, [species; parameters; variables], tiv)
 
     # Checks for input errors.
     (sum(length.([reaction_lines, option_lines])) != length(ex.args)) &&
@@ -414,7 +414,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
 
     # Creates expressions corresponding to actual code from the internal DSL representation.
     sexprs = get_sexpr(species_extracted, options; iv_symbols = ivs)
-    vexprs = get_sexpr(vars_extracted, options, :variables)
+    vexprs = get_sexpr(vars_extracted, options, :variables; iv_symbols = ivs)
     pexprs = get_pexpr(parameters_extracted, options)
     ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps")
     vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars")
@@ -852,12 +852,13 @@ function read_equations_options(options, variables_declared)
             in(diff_var, variables_declared) || push!(vars_extracted, diff_var)
         end
     end
-
+    
     return vars_extracted, add_default_diff, equations
 end
 
-# Creates an expression declaring differentials.
-function create_differential_expr(options, add_default_diff, used_syms)
+# Creates an expression declaring differentials. Here, `tiv` is the time independent variables, 
+# which is used by the default differential (if it is used).
+function create_differential_expr(options, add_default_diff, used_syms, tiv)
     # Creates the differential expression.
     # If differentials was provided as options, this is used as the initial expression.
     # If the default differential (D(...)) was used in equations, this is added to the expression.
@@ -875,7 +876,7 @@ function create_differential_expr(options, add_default_diff, used_syms)
     # If the default differential D has been used, but not pre-declared using the @differenitals
     # options, add this declaration to the list of declared differentials.
     if add_default_diff && !any(diff_dec.args[1] == :D for diff_dec in diffexpr.args) 
-        push!(diffexpr.args, :(D = Differential($(DEFAULT_IV_SYM))))
+        push!(diffexpr.args, :(D = Differential($(tiv))))
     end
 
     return diffexpr

src/reactionsystem.jl

@@ -312,8 +312,8 @@ end
 ### Reaction Acessor Functions ###
 
 # Overwrites functions in ModelingToolkit to give the correct input.
-is_diff_equation(rx::Reaction) = false
-is_alg_equation(rx::Reaction) = false
+ModelingToolkit.is_diff_equation(rx::Reaction) = false
+ModelingToolkit.is_alg_equation(rx::Reaction) = false
 
 ################################## Reaction Complexes ####################################
 
@@ -725,6 +725,11 @@ function findvars!(ps, us, exprtosearch, ivs, vars)
     end
     empty!(vars)
 end
+# Special dispatch for equations, applied `findvars!` to left-hand and right-hand sides.
+function findvars!(ps, us, eq_to_search::Equation, ivs, vars)
+    findvars!(ps, us, eq_to_search.lhs, ivs, vars)
+    findvars!(ps, us, eq_to_search.lhs, ivs, vars)
+end
 
 # Called internally (whether DSL-based or programmtic model creation is used). 
 # Creates a sorted reactions + equations vector, also ensuring reaction is first in this vector.
@@ -813,22 +818,14 @@ function find_event_vars!(ps, us, event, ivs, vars)
     # If not, it is a vector of conditions and we must check each.
     if conds isa Vector
         for cond in conds
-            # For continious events the conditions are equations (with lhs and rhs).
-            # For discrete events, they are single expressions.
-            if cond isa Equation
-                findvars!(ps, us, cond.lhs, ivs, vars)
-                findvars!(ps, us, cond.rhs, ivs, vars)
-            else
-                findvars!(ps, us, cond, ivs, vars)
-            end
+            findvars!(ps, us, cond, ivs, vars)
         end
     else
         findvars!(ps, us, conds, ivs, vars)
     end
     # The affects is always a vector of equations. Here, we handle the lhs and rhs separately.
     for affect in affects
-        findvars!(ps, us, affect.lhs, ivs, vars)
-        findvars!(ps, us, affect.rhs, ivs, vars)
+        findvars!(ps, us, cond, ivs, vars)
     end
 end
 """
@@ -1616,6 +1613,9 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
                          as_odes = false, include_zero_odes)
     eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
 
+    # Throws a warning if there are differential equations in non-standard format.
+    # Next, sets all differential terms to `0`.
+    nonlinear_convert_differentials_check(eqs, nonspecies(rs), get_iv(rs))
     eqs = [remove_diffs(eq.lhs) ~ remove_diffs(eq.rhs) for eq in eqs]
 
 
@@ -1637,7 +1637,32 @@ function remove_diffs(expr)
 end
 diff_2_zero(expr) = (Symbolics.is_derivative(expr) ? 0.0 : expr)
 
-    
+# Ideally, when `ReactionSystem`s are converted to `NonlinearSystem`s, any coupled ODEs should be 
+# on the form D(X) ~ ..., where lhs is the time derivative w.r.t. a single variable, and the rhs
+# 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(eqs, vars, iv)
+    for eq in eqs
+        # For each equation (that is not a reaction), 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(iv)) || 
+                (length(arguments(eq.lhs)) != 1) ||
+                !any(isequal(arguments(eq.lhs)[1]), vars)
+            @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: 
+                    (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
+        end
+    end
+end
 
 """
 ```julia
@@ -1820,7 +1845,7 @@ function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0, tspan::Tuple,
     return DiscreteProblem(jsys, u0map, tspan, pmap, args...; kwargs...)
 end
 
-# JumpProblem from AbstractReactionNetworkg
+# JumpProblem from AbstractReactionNetwork
 function JumpProcesses.JumpProblem(rs::ReactionSystem, prob, aggregator, args...;
                                    name = nameof(rs),
                                    combinatoric_ratelaws = get_combinatoric_ratelaws(rs),

test/dsl/dsl_options.jl

@@ -687,7 +687,7 @@ let
 end
 
 
-### Hybrid CRN/Equations Models ###
+### Coupled CRN/Equations Models ###
 
 # Checks creation of basic network.
 # Check indexing of output solution.

test/extensions/bifurcation_kit.jl

@@ -180,7 +180,7 @@ end
 
 ### Other Tests ###
 
-# Checks that bifurcation diagrams can be computed for hybrid CRN/DAE systems.
+# Checks that bifurcation diagrams can be computed for coupled CRN/DAE systems.
 let
     # Prepares the model (production/degradation of X, with equations for volume and X concentration).
     rs = @reaction_network begin

test/extensions/homotopy_continuation.jl

@@ -109,7 +109,7 @@ end
 
 ### Other Tests ###
 
-# Tests that homotopy continuation can be applied to hybrid DAE/CRN systems.
+# Tests that homotopy continuation can be applied to coupled DAE/CRN systems.
 let
     # Prepares the model (production/degradation of X, with equations for volume and X concentration).
     rs = @reaction_network begin

test/extensions/structural_identifiability.jl

@@ -312,7 +312,7 @@ end
 
 ### Other Tests ###
 
-# Checks that identifiability can be assessed for hybrid CRN/DAE systems.
+# Checks that identifiability can be assessed for coupled CRN/DAE systems.
 let
     rs = @reaction_network begin
         @parameters k c1 c2

test/programmatic_model_creation/component_based_model_creation.jl

@@ -363,7 +363,6 @@ let
     rn2 = extend(osys, rn)
     rn2 = complete(rn2)
     rnodes = complete(convert(ODESystem, rn2))
-    # @test_throws ErrorException convert(NonlinearSystem, rn2)
 
     # Ensure right number of equations are generated.
     @variables G(t)