save progress

src/Catalyst.jl

@@ -41,6 +41,7 @@ import ModelingToolkit: check_variables,
 
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
 import MacroTools, Graphs
+using MacroTools: striplines
 import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne
 import DataStructures: OrderedDict, OrderedSet
 import Parameters: @with_kw_noshow

src/chemistry_functionality.jl

@@ -63,11 +63,11 @@ t = default_t()
 @compound CO2(t) ~ C + 2O
 ```
 
-Notes: 
+Notes:
 - The component species must be defined before using the `@compound` macro.
 """
 macro compound(expr)
-    make_compound(MacroTools.striplines(expr))
+    make_compound(striplines(expr))
 end
 
 # Declares compound error messages:
@@ -83,12 +83,12 @@ function make_compound(expr)
         error(COMPOUND_CREATION_ERROR_BAD_SEPARATOR)
 
     # Loops through all components, add the component and the coefficients to the corresponding vectors
-    # Cannot extract directly using e.g. "getfield.(composition, :reactant)" because then 
+    # Cannot extract directly using e.g. "getfield.(composition, :reactant)" because then
     # we get something like :([:C, :O]), rather than :([C, O]).
     composition = Catalyst.recursive_find_reactants!(expr.args[3], 1,
         Vector{ReactantStruct}(undef, 0))
-    components = :([])                                      # Becomes something like :([C, O]).                                         
-    coefficients = :([])                                    # Becomes something like :([1, 2]). 
+    components = :([])                                      # Becomes something like :([C, O]).
+    coefficients = :([])                                    # Becomes something like :([1, 2]).
     for comp in composition
         push!(components.args, comp.reactant)
         push!(coefficients.args, comp.stoichiometry)
@@ -110,13 +110,13 @@ function make_compound(expr)
                for comp in $components])))
 
     # Creates the found expressions that will create the compound species.
-    # The `Expr(:escape, :(...))` is required so that the expressions are evaluated in 
-    # the scope the users use the macro in (to e.g. detect already exiting species).     
+    # The `Expr(:escape, :(...))` is required so that the expressions are evaluated in
+    # the scope the users use the macro in (to e.g. detect already exiting species).
     # Creates something like (where `compound_ivs` and `component_ivs` evaluates to all the compound's and components' ivs):
     #   `@species CO2(..)`
     #   `isempty([])` && length(component_ivs) && error("When ...)
-    #   `CO2 = CO2(component_ivs..)` 
-    #   `issetequal(compound_ivs, component_ivs) || error("The ...)` 
+    #   `CO2 = CO2(component_ivs..)`
+    #   `issetequal(compound_ivs, component_ivs) || error("The ...)`
     #   `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)`
     #   `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [C, O])`
     #   `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [1, 2])`
@@ -159,7 +159,7 @@ Macro that creates several compound species, which each is composed of smaller c
 Example:
 ```julia
 t = default_t()
-@species C(t) H(t) O(t) 
+@species C(t) H(t) O(t)
 @compounds
     CH4(t) = C + 4H
     O2(t) = 2O
@@ -168,11 +168,11 @@ t = default_t()
 end
 ```
 
-Notes: 
+Notes:
 - The component species must be defined before using the `@compound` macro.
 """
 macro compounds(expr)
-    make_compounds(MacroTools.striplines(expr))
+    make_compounds(striplines(expr))
 end
 
 # Function managing the @compound macro.
@@ -183,7 +183,7 @@ function make_compounds(expr)
     # For each compound in `expr`, creates the set of 7 compound creation lines (using `make_compound`).
     # Next, loops through all 7*[Number of compounds] lines and add them to compound_declarations.
     compound_calls = [Catalyst.make_compound(line) for line in expr.args]
-    for compound_call in compound_calls, line in MacroTools.striplines(compound_call).args
+    for compound_call in compound_calls, line in striplines(compound_call).args
         push!(compound_declarations.args, line)
     end
 
@@ -249,7 +249,7 @@ brxs = balance_reaction(rx) # No solution.
 
 Notes:
 - Balancing reactions that contain compounds of compounds is currently not supported.
-- A reaction may not always yield a single solution; it could have an infinite number of solutions or none at all. When there are multiple solutions, a vector of all possible `Reaction` objects is returned. However, substrates and products may be interchanged as we currently do not solve for a linear combination that maintains the set of substrates and products. 
+- A reaction may not always yield a single solution; it could have an infinite number of solutions or none at all. When there are multiple solutions, a vector of all possible `Reaction` objects is returned. However, substrates and products may be interchanged as we currently do not solve for a linear combination that maintains the set of substrates and products.
 - If the reaction cannot be balanced, an empty `Reaction` vector is returned.
 """
 function balance_reaction(reaction::Reaction)
@@ -369,16 +369,16 @@ From a system, creates a new system where each reaction is a balanced version of
 reaction of the original system. For more information, consider the `balance_reaction` function
 (which is internally applied to each system reaction).
 
-Arguments 
+Arguments
 - `rs`: The reaction system that should be balanced.
 
 Notes:
 - If any reaction in the system cannot be balanced, throws an error.
-- If any reaction in the system have an infinite number of potential reactions, throws an error. 
+- If any reaction in the system have an infinite number of potential reactions, throws an error.
 Here, it would be possible to generate a valid reaction, however, no such routine is currently
 implemented in `balance_system`.
 - `balance_system` will not modify reactions of subsystems to the input system. It is recommended
-not to apply `balance_system` to non-flattened systems. 
+not to apply `balance_system` to non-flattened systems.
 """
 function balance_system(rs::ReactionSystem)
     @set! rs.eqs = CatalystEqType[get_balanced_reaction(eq) for eq in get_eqs(rs)]
@@ -391,7 +391,7 @@ end
 function get_balanced_reaction(rx::Reaction)
     brxs = balance_reaction(rx)
 
-    # In case there are no, or multiple, solutions to the balancing problem. 
+    # In case there are no, or multiple, solutions to the balancing problem.
     if isempty(brxs)
         error("Could not balance reaction `$rx`, unable to create a balanced `ReactionSystem`.")
     end

src/dsl.jl

@@ -141,17 +141,17 @@ ReactionSystems generated through `@reaction_network` are complete.
 """
 macro reaction_network(name::Symbol, ex::Expr)
     :(complete($(make_reaction_system(
-        MacroTools.striplines(ex); name = :($(QuoteNode(name)))))))
+        striplines(ex); name = :($(QuoteNode(name)))))))
 end
 
 # Allows @reaction_network $name begin ... to interpolate variables storing a name.
 macro reaction_network(name::Expr, ex::Expr)
     :(complete($(make_reaction_system(
-        MacroTools.striplines(ex); name = :($(esc(name.args[1])))))))
+        striplines(ex); name = :($(esc(name.args[1])))))))
 end
 
 macro reaction_network(ex::Expr)
-    ex = MacroTools.striplines(ex)
+    ex = striplines(ex)
 
     # no name but equations: @reaction_network begin ... end ...
     if ex.head == :block
@@ -179,16 +179,16 @@ Equivalent to `@reaction_network` except the generated `ReactionSystem` is not m
 complete.
 """
 macro network_component(name::Symbol, ex::Expr)
-    make_reaction_system(MacroTools.striplines(ex); name = :($(QuoteNode(name))))
+    make_reaction_system(striplines(ex); name = :($(QuoteNode(name))))
 end
 
 # Allows @network_component $name begin ... to interpolate variables storing a name.
 macro network_component(name::Expr, ex::Expr)
-    make_reaction_system(MacroTools.striplines(ex); name = :($(esc(name.args[1]))))
+    make_reaction_system(striplines(ex); name = :($(esc(name.args[1]))))
 end
 
 macro network_component(ex::Expr)
-    ex = MacroTools.striplines(ex)
+    ex = striplines(ex)
 
     # no name but equations: @network_component begin ... end ...
     if ex.head == :block
@@ -315,7 +315,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters_declared = extract_syms(options, :parameters)
     variables_declared = extract_syms(options, :variables)
 
-    # Reads equations. 
+    # Reads equations.
     vars_extracted, add_default_diff, equations = read_equations_options(
         options, variables_declared)
     variables = vcat(variables_declared, vars_extracted)
@@ -343,7 +343,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     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, including ones inferred from the reactions.
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
         variables)))
     species_extracted, parameters_extracted = extract_species_and_parameters!(
@@ -384,7 +384,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
@@ -435,7 +435,7 @@ function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(u
             push_reactions!(reactions, reaction.args[3], reaction.args[2],
                 rate, metadata, arrow)
         else
-            throw("Malformed reaction, invalid arrow type used in: $(MacroTools.striplines(line))")
+            throw("Malformed reaction, invalid arrow type used in: $(striplines(line))")
         end
     end
     return reactions
@@ -651,7 +651,7 @@ function read_events_option(options, event_type::Symbol)
         error("Trying to read an unsupported event type.")
     end
     events_input = haskey(options, event_type) ? options[event_type].args[3] :
-                   MacroTools.striplines(:(begin end))
+                   striplines(:(begin end))
     events_input = option_block_form(events_input)
 
     # Goes through the events, checks for errors, and adds them to the output vector.
@@ -727,7 +727,7 @@ function create_differential_expr(options, add_default_diff, used_syms, tiv)
     # 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.
     diffexpr = (haskey(options, :differentials) ? options[:differentials].args[3] :
-                MacroTools.striplines(:(begin end)))
+                striplines(:(begin end)))
     diffexpr = option_block_form(diffexpr)
 
     # Goes through all differentials, checking that they are correctly formatted and their symbol is not used elsewhere.
@@ -783,11 +783,6 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted)
             # For Observables that have already been declared using @species/@variables,
             # or are interpolated, this parts should not be carried out.
             if !((obs_name in species_n_vars_declared) || is_escaped_expr(obs_eq.args[2]))
-                # Appends (..) to the observable (which is later replaced with the extracted ivs).
-                # Adds the observable to the first line of the output expression (starting with `@variables`).
-                obs_expr = insert_independent_variable(obs_eq.args[2], :(..))
-                push!(observed_vars.args[1].args, obs_expr)
-
                 # Adds a line to the `observed_vars` expression, setting the ivs for this observable.
                 # Cannot extract directly using e.g. "getfield.(dependants_structs, :reactant)" because
                 # then we get something like :([:X1, :X2]), rather than :([X1, X2]).
@@ -798,8 +793,11 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted)
                            for dep in $dep_var_expr])))
                 ivs_get_expr_sorted = :(sort($(ivs_get_expr);
                     by = iv -> findfirst(MT.getname(iv) == ivs for ivs in $ivs_sorted)))
-                push!(observed_vars.args,
-                    :($obs_name = $(obs_name)($(ivs_get_expr_sorted)...)))
+
+                ivs_get_expr_sorted = :([τ, x])
+                println(ivs_get_expr_sorted)
+                obs_expr = insert_independent_variable(obs_eq.args[2], :(([τ, x])...))
+                push!(observed_vars.args[1].args, obs_expr)
             end
 
             # In case metadata was given, this must be cleared from `observed_eqs`.
@@ -813,7 +811,8 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted)
         end
 
         # If nothing was added to `observed_vars`, it has to be modified not to throw an error.
-        (length(observed_vars.args) == 1) && (observed_vars = :())
+        (striplines(observed_vars) == striplines(Expr(:block, :(@variables)))) &&
+            (observed_vars = :())
     else
         # If option is not used, return empty expression and vector.
         observed_vars = :()
@@ -917,7 +916,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/spatial_reaction_systems/spatial_reactions.jl

@@ -43,7 +43,7 @@ end
 
 # Macro for creating a `TransportReaction`.
 macro transport_reaction(rateex::ExprValues, species::ExprValues)
-    make_transport_reaction(MacroTools.striplines(rateex), species)
+    make_transport_reaction(striplines(rateex), species)
 end
 function make_transport_reaction(rateex, species)
     # Handle interpolation of variables

test/dsl/dsl_options.jl

@@ -565,6 +565,35 @@ let
     @test sol[:X][1] == u0[:X1]^2 + ps[:op_1]*(u0[:X2] + 2*u0[:X3]) + u0[:X1]*u0[:X4]/ps[:op_2] + ps[:p]
 end
 
+# Checks that models created w/o specifying `@variables` for observables are identical.
+let
+    # With default ivs.
+    rn1 = @reaction_network rn begin
+        @variables X(t) X1(t) X2(t)
+        @observables X ~ X1 + X2
+    end
+    rn2 = @reaction_network rn begin
+        @variables X1(t) X2(t)
+        @observables X ~ X1 + X2
+    end
+    @test isequal(rn1, rn2)
+    @test isequal(rn1.X, rn2.X)
+
+    # With non-default ivs.
+    rn3 = @reaction_network rn begin
+        @ivs τ x
+        @variables X(τ,x) X1(τ,x) X2(τ,x)
+        @observables X ~ X1 + X2
+    end
+    rn4 = @reaction_network rn begin
+        @ivs τ x
+        @variables X1(τ,x) X2(τ,x)
+        @observables X ~ X1 + X2
+    end
+    @test isequal(rn3, rn4)
+    @test isequal(rn3.X, rn4.X)
+end
+
 # Checks that ivs are correctly found.
 let
     rn = @reaction_network begin
@@ -951,7 +980,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 +994,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 +1007,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 +1038,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