save prorgess

Author: TorkelE

src/dsl.jl

@@ -315,11 +315,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters_declared = extract_syms(options, :parameters)
     variables_declared = extract_syms(options, :variables)
 
-    # Reads equations. 
-    vars_extracted, add_default_diff, equations = read_equations_options(
-        options, variables_declared)
-    variables = vcat(variables_declared, vars_extracted)
-
     # Handle independent variables
     if haskey(options, :ivs)
         ivs = Tuple(extract_syms(options, :ivs))
@@ -339,23 +334,29 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         combinatoric_ratelaws = true
     end
 
-    # Reads observables.
-    observed_vars, observed_eqs, obs_syms = read_observed_options(
-        options, [species_declared; variables], all_ivs)
-
-    # Collect species and parameters, including ones inferred from the reactions. 
+    # Collect species and parameters.
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
-        variables)))
+        variables_declared)))
     species_extracted, parameters_extracted = extract_species_and_parameters!(
         reactions, declared_syms)
 
     species = vcat(species_declared, species_extracted)
     parameters = vcat(parameters_declared, parameters_extracted)
 
+    # Reads equations.
+    designated_syms = [species; parameters; variables_declared]
+    vars_extracted, add_default_diff, equations = read_equations_options(
+        options, designated_syms)
+    variables = vcat(variables_declared, vars_extracted)
+
     # Create differential expression.
     diffexpr = create_differential_expr(
         options, add_default_diff, [species; parameters; variables], tiv)
 
+    # Reads observables.
+    observed_vars, observed_eqs, obs_syms = read_observed_options(
+        options, [species_declared; variables], all_ivs)
+
     # Checks for input errors.
     (sum(length.([reaction_lines, option_lines])) != length(ex.args)) &&
         error("@reaction_network input contain $(length(ex.args) - sum(length.([reaction_lines,option_lines]))) malformed lines.")
@@ -384,7 +385,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         push!(rxexprs.args, equation)
     end
 
-    # Output code corresponding to the reaction system. 
+    # Output code corresponding to the reaction system.
     quote
         $ivexpr
         $ps
@@ -682,7 +683,7 @@ end
 # `vars_extracted`: A vector with extracted variables (lhs in pure differential equations only).
 # `dtexpr`: If a differential equation is defined, the default derivative (D ~ Differential(t)) must be defined.
 # `equations`: a vector with the equations provided.
-function read_equations_options(options, variables_declared)
+function read_equations_options(options, syms_declared)
     # Prepares the equations. First, extracts equations from provided option (converting to block form if required).
     # Next, uses MTK's `parse_equations!` function to split input into a vector with the equations.
     eqs_input = haskey(options, :equations) ? options[:equations].args[3] : :(begin end)
@@ -694,30 +695,35 @@ 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 = Vector{Symbol}()
+    vars_extracted = OrderedSet{Union{Symbol, Expr}}()
+    excluded_syms = syms_declared
     add_default_diff = false
     for eq in equations
         if (eq.head != :call) || (eq.args[1] != :~)
             error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")
         end
 
         # Checks if the equation have the format D(X) ~ ... (where X is a symbol). This means that the
-        # default differential has been used. X is added as a declared variable to the system, and
-        # we make a note that a differential D = Differential(iv) should be made as well.
+        # default differential has been used and we make a note that it should be decalred in the DSL output.
         lhs = eq.args[2]
-        # if lhs: is an expression. Is a function call. The function's name is D. Calls a single symbol.
+        # If lhs: is an expression. Is a function call. The function's name is D. It has a single argument.
         if (lhs isa Expr) && (lhs.head == :call) && (lhs.args[1] == :D) &&
            (lhs.args[2] isa Symbol)
             diff_var = lhs.args[2]
             if in(diff_var, forbidden_symbols_error)
                 error("A forbidden symbol ($(diff_var)) was used as an variable in this differential equation: $eq")
             end
-            add_default_diff = true
-            in(diff_var, variables_declared) || push!(vars_extracted, diff_var)
+            if !add_default_diff
+                add_default_diff = true
+                excluded_syms = [excluded_syms; :D]
+            end
         end
+
+        # Any undecalred symbolic variables encountered should be extracted as variables.
+        add_syms_from_expr!(vars_extracted, eq, excluded_syms)
     end
 
-    return vars_extracted, add_default_diff, equations
+    return collect(vars_extracted), add_default_diff, equations
 end
 
 # Creates an expression declaring differentials. Here, `tiv` is the time independent variables,
@@ -917,7 +923,7 @@ end
 
 ### Generic Expression Manipulation ###
 
-# Recursively traverses an expression and escapes all the user-defined functions. Special function calls like "hill(...)" are not expanded. 
+# Recursively traverses an expression and escapes all the user-defined functions. Special function calls like "hill(...)" are not expanded.
 function recursive_escape_functions!(expr::ExprValues)
     (typeof(expr) != Expr) && (return expr)
     foreach(i -> expr.args[i] = recursive_escape_functions!(expr.args[i]),

src/reactionsystem.jl

@@ -97,9 +97,9 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     stronglinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
     terminallinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
 
-    checkedrobust::Bool = false 
+    checkedrobust::Bool = false
     robustspecies::Vector{Int} = Vector{Int}(undef, 0)
-    deficiency::Int = -1 
+    deficiency::Int = -1
 end
 #! format: on
 
@@ -215,11 +215,11 @@ end
 
 ### ReactionSystem Structure ###
 
-""" 
+"""
 WARNING!!!
 
-The following variable is used to check that code that should be updated when the `ReactionSystem` 
-fields are updated has in fact been updated. Do not just blindly update this without first checking 
+The following variable is used to check that code that should be updated when the `ReactionSystem`
+fields are updated has in fact been updated. Do not just blindly update this without first checking
 all such code and updating it appropriately (e.g. serialization). Please use a search for
 `reactionsystem_fields` throughout the package to ensure all places which should be updated, are updated.
 """
@@ -318,7 +318,7 @@ struct ReactionSystem{V <: NetworkProperties} <:
     """
     discrete_events::Vector{MT.SymbolicDiscreteCallback}
     """
-    Metadata for the system, to be used by downstream packages. 
+    Metadata for the system, to be used by downstream packages.
     """
     metadata::Any
     """
@@ -480,10 +480,10 @@ function ReactionSystem(iv; kwargs...)
     ReactionSystem(Reaction[], iv, [], []; kwargs...)
 end
 
-# Called internally (whether DSL-based or programmatic model creation is used). 
+# Called internally (whether DSL-based or programmatic model creation is used).
 # Creates a sorted reactions + equations vector, also ensuring reaction is first in this vector.
-# Extracts potential species, variables, and parameters from the input (if not provided as part of 
-# the model creation) and creates the corresponding vectors. 
+# Extracts potential species, variables, and parameters from the input (if not provided as part of
+# the model creation) and creates the corresponding vectors.
 # While species are ordered before variables in the unknowns vector, this ordering is not imposed here,
 # but carried out at a later stage.
 function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
@@ -495,7 +495,7 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
     any(in(obs_vars), us_in) &&
         error("Found an observable in the list of unknowns. This is not allowed.")
 
-    # Creates a combined iv vector (iv and sivs). This is used later in the function (so that 
+    # Creates a combined iv vector (iv and sivs). This is used later in the function (so that
     # independent variables can be excluded when encountered quantities are added to `us` and `ps`).
     t = value(iv)
     ivs = Set([t])
@@ -524,6 +524,8 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
     # Extracts any species, variables, and parameters that occur in (non-reaction) equations.
     # Creates the new reactions + equations vector, `fulleqs` (sorted reactions first, equations next).
     if !isempty(eqs)
+        println(eqs)
+        println(iv)
         osys = ODESystem(eqs, iv; name = gensym())
         fulleqs = CatalystEqType[rxs; equations(osys)]
         union!(us, unknowns(osys))
@@ -560,7 +562,7 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
     end
     psv = collect(new_ps)
 
-    # Passes the processed input into the next `ReactionSystem` call.    
+    # Passes the processed input into the next `ReactionSystem` call.
     ReactionSystem(fulleqs, t, usv, psv; spatial_ivs, continuous_events,
         discrete_events, observed, kwargs...)
 end
@@ -1062,8 +1064,8 @@ end
 
 ### General `ReactionSystem`-specific Functions ###
 
-# Checks if the `ReactionSystem` structure have been updated without also updating the 
-# `reactionsystem_fields` constant. If this is the case, returns `false`. This is used in 
+# Checks if the `ReactionSystem` structure have been updated without also updating the
+# `reactionsystem_fields` constant. If this is the case, returns `false`. This is used in
 # certain functionalities which would break if the `ReactionSystem` structure is updated without
 # also updating these functionalities.
 function reactionsystem_uptodate_check()
@@ -1241,7 +1243,7 @@ end
 ### `ReactionSystem` Remaking ###
 
 """
-    remake_ReactionSystem_internal(rs::ReactionSystem; 
+    remake_ReactionSystem_internal(rs::ReactionSystem;
         default_reaction_metadata::Vector{Pair{Symbol, T}} = Vector{Pair{Symbol, Any}}()) where {T}
 
 Takes a `ReactionSystem` and remakes it, returning a modified `ReactionSystem`. Modifications depend
@@ -1274,7 +1276,7 @@ function set_default_metadata(rs::ReactionSystem; default_reaction_metadata = []
     # Currently, `noise_scaling` is the only relevant metadata supported this way.
     drm_dict = Dict(default_reaction_metadata)
     if haskey(drm_dict, :noise_scaling)
-        # Finds parameters, species, and variables in the noise scaling term.       
+        # Finds parameters, species, and variables in the noise scaling term.
         ns_expr = drm_dict[:noise_scaling]
         ns_syms = [Symbolics.unwrap(sym) for sym in get_variables(ns_expr)]
         ns_ps = Iterators.filter(ModelingToolkit.isparameter, ns_syms)
@@ -1414,7 +1416,7 @@ function ModelingToolkit.compose(sys::ReactionSystem, systems::AbstractArray; na
         MT.collect_scoped_vars!(newunknowns, newparams, ssys, iv)
     end
 
-    if !isempty(newunknowns) 
+    if !isempty(newunknowns)
         @set! sys.unknowns = union(get_unknowns(sys), newunknowns)
         sort!(get_unknowns(sys), by = !isspecies)
         @set! sys.species = filter(isspecies, get_unknowns(sys))

test/dsl/dsl_options.jl

@@ -901,6 +901,26 @@ let
     @test 5*sol[:Y][end] ≈ sol[:S][end] + sol[:X][end]
 end
 
+# Tests that the correct symbolic variables are infered as species, variables, and paraemters.
+let
+    rn = @reaction_network begin
+        @parameters p1 p2
+        @species X1(t) X2(t)
+        @variables W(t)
+        @equations begin
+            D(V1) ~ p1 * X1 - k1 * V1 + W
+            k2 * V2 ~ D(V2) + p2 * X2
+            V3 + X3 ~ V1^2 + X2^2
+        end
+        (k1, k2), X1 <--> X2
+        (k3, k4), X2 <--> X3
+    end
+
+    @test issetequal(species(rn), @species X1(t) X2(t) X3(t))
+    @test issetequal(parameters(rn), @parameters p1 p2 k1 k2 k3 k4)
+    @test issetequal(nonspecies(rn), @variables V1(t) V2(t) V3(t) W(t))
+end
+
 # Tests that various erroneous declarations throw errors.
 let
     # Using = instead of ~ (for equation).
@@ -951,7 +971,7 @@ let
     @test isequal(rl, k1*A^2)
 end
 
-# Test whether user-defined functions are properly expanded in equations. 
+# Test whether user-defined functions are properly expanded in equations.
 let
     f(A, t) = 2*A*t
 
@@ -965,7 +985,7 @@ let
     @test isequal(equations(rn)[1], D(A) ~ 2*A*t)
 
 
-    # Test whether expansion happens properly for unregistered/registered functions. 
+    # Test whether expansion happens properly for unregistered/registered functions.
     hill_unregistered(A, v, K, n) = v*(A^n) / (A^n + K^n)
     rn2 = @reaction_network begin
         @parameters v K n
@@ -978,7 +998,7 @@ let
 
     hill2(A, v, K, n) = v*(A^n) / (A^n + K^n)
     @register_symbolic hill2(A, v, K, n)
-    # Registered symbolic function should not expand. 
+    # Registered symbolic function should not expand.
     rn2r = @reaction_network begin
         @parameters v K n
         @equations D(A) ~ hill2(A, v, K, n)
@@ -1009,9 +1029,9 @@ let
     @named rn3_sym = ReactionSystem(eq, t)
     rn3_sym = complete(rn3_sym)
     @test isequivalent(rn3, rn3_sym)
-    
-    
-    # Test more complicated expression involving both registered function and a user-defined function. 
+
+
+    # Test more complicated expression involving both registered function and a user-defined function.
     g(A, K, n) = A^n + K^n
     rn4 = @reaction_network begin
         @parameters v K n