Merge pull request #988 from isaacsas/t_as_parameter_fix

make `t` a parameter now

Author: isaacsas

docs/src/model_creation/model_file_loading_and_export.md

@@ -30,7 +30,7 @@ Here, `include` is used to execute the Julia code from any file. This means that
 let
 
 # Independent variable:
-@variables t
+@parameters t
 
 # Parameters:
 ps = @parameters kB kD kP
@@ -52,7 +52,7 @@ complete(rs)
 end
 ```
 !!! note
-    The code that `save_reactionsystem` prints uses [programmatic modelling](@ref programmatic_CRN_construction) to generate the written model. 
+    The code that `save_reactionsystem` prints uses [programmatic modelling](@ref programmatic_CRN_construction) to generate the written model.
 
 In addition to transferring models between Julia sessions, the `save_reactionsystem` function can also be used or print a model to a text file where you can easily inspect its components.
 
@@ -76,7 +76,7 @@ A general-purpose format for storing CRN models is so-called .net files. These c
 using ReactionNetworkImporters
 prn = loadrxnetwork(BNGNetwork(), "repressilator.net")
 ```
-Here, .net files not only contain information regarding the reaction network itself, but also the numeric values (initial conditions and parameter values) required for simulating it. Hence, `loadrxnetwork` generates a `ParsedReactionNetwork` structure, containing all this information. You can access the model as `prn.rn`, the initial conditions as `prn.u0`, and the parameter values as `prn.p`. Furthermore, these initial conditions and parameter values are also made [*default* values](@ref dsl_advanced_options_default_vals) of the model. 
+Here, .net files not only contain information regarding the reaction network itself, but also the numeric values (initial conditions and parameter values) required for simulating it. Hence, `loadrxnetwork` generates a `ParsedReactionNetwork` structure, containing all this information. You can access the model as `prn.rn`, the initial conditions as `prn.u0`, and the parameter values as `prn.p`. Furthermore, these initial conditions and parameter values are also made [*default* values](@ref dsl_advanced_options_default_vals) of the model.
 
 A parsed reaction network's content can then be provided to various problem types for simulation. E.g. here we perform an ODE simulation of our repressilator model:
 ```julia
@@ -128,7 +128,7 @@ Above, we described how to use SBMLImporter to import SBML files. Alternatively,
 While CRN models can be represented through various file formats, they can also be represented in various matrix forms. E.g. a CRN with $m$ species and $n$ reactions (and with constant rates) can be represented with either
 - An $mxn$ substrate matrix (with each species's substrate stoichiometry in each reaction) and an $nxm$ product matrix (with each species's product stoichiometry in each reaction).
 
-Or 
+Or
 - An $mxn$ complex stoichiometric matrix (...) and a $2mxn$ incidence matrix (...).
 
 The advantage of these forms is that they offer a compact and very general way to represent a large class of CRNs. ReactionNetworkImporters have the functionality for converting matrices of these forms directly into Catalyst `ReactionSystem` models. Instructions on how to do this are available in [ReactionNetworkImporter's documentation](https://docs.sciml.ai/ReactionNetworkImporters/stable/#Loading-a-matrix-representation).

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

@@ -1,17 +1,17 @@
 ### Structural Identifiability ODE Creation ###
 
-# For a reaction system, list of measured quantities and known parameters, generate a StructuralIdentifiability compatible ODE. 
+# For a reaction system, list of measured quantities and known parameters, generate a StructuralIdentifiability compatible ODE.
 """
 make_si_ode(rs::ReactionSystem; measured_quantities=observed(rs), known_p = [], ignore_no_measured_warn=false)
 
-Creates a reaction rate equation ODE system of the form used within the StructuralIdentifiability.jl package. The output system is compatible with all StructuralIdentifiability functions. 
+Creates a reaction rate equation ODE system of the form used within the StructuralIdentifiability.jl package. The output system is compatible with all StructuralIdentifiability functions.
 
 Arguments:
 - `rs::ReactionSystem`; The reaction system we wish to convert to an ODE.
-- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.x`, or the symbol `:x`). 
+- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.x`, or the symbol `:x`).
 - `known_p = []`: List of parameters for which their values are assumed to be known.
-- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty. 
-- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly). 
+- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty.
+- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly).
 
 Example:
 ```julia
@@ -45,10 +45,10 @@ Applies StructuralIdentifiability.jl's `assess_local_identifiability` function t
 
 Arguments:
 - `rs::ReactionSystem`; The reaction system for which we wish to compute structural identifiability of the associated reaction rate equation ODE model.
-- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.s`, or the symbol `:x`). 
-- `known_p = []`: List of parameters which values are known. 
-- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty. 
-- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly). 
+- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.s`, or the symbol `:x`).
+- `known_p = []`: List of parameters which values are known.
+- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty.
+- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly).
 
 Example:
 ```julia
@@ -85,10 +85,10 @@ Applies StructuralIdentifiability.jl's `assess_identifiability` function to a Ca
 
 Arguments:
 - `rs::ReactionSystem`; The reaction system we wish to compute structural identifiability for.
-- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.s`, or the symbol `:x`). 
-- `known_p = []`: List of parameters which values are known. 
-- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty. 
-- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly). 
+- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.s`, or the symbol `:x`).
+- `known_p = []`: List of parameters which values are known.
+- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty.
+- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly).
 
 Example:
 ```julia
@@ -125,10 +125,10 @@ Applies StructuralIdentifiability.jl's `find_identifiable_functions` function to
 
 Arguments:
 - `rs::ReactionSystem`; The reaction system we wish to compute structural identifiability for.
-- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.s`, or the symbol `:x`). 
-- `known_p = []`: List of parameters which values are known. 
-- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty. 
-- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly). 
+- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.s`, or the symbol `:x`).
+- `known_p = []`: List of parameters which values are known.
+- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty.
+- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly).
 
 Example:
 ```julia
@@ -158,7 +158,7 @@ end
 
 ### Helper Functions ###
 
-# From a reaction system, creates the corresponding MTK-style ODESystem for SI application 
+# From a reaction system, creates the corresponding MTK-style ODESystem for SI application
 # Also compute the, later needed, conservation law equations and list of system symbols (unknowns and parameters).
 function make_osys(rs::ReactionSystem; remove_conserved = true)
     # Creates the ODESystem corresponding to the ReactionSystem (expanding functions and flattening it).
@@ -200,7 +200,7 @@ function make_measured_quantities(
 
     # Creates one internal observation variable for each measured quantity (`___internal_observables`).
     # Creates a vector of equations, setting each measured quantity equal to one observation variable.
-    @variables t (___internal_observables(Catalyst.get_iv(rs)))[1:length(mqs)]
+    @variables (___internal_observables(Catalyst.get_iv(rs)))[1:length(mqs)]
     return Equation[(q isa Equation) ? q : (___internal_observables[i] ~ q)
                     for (i, q) in enumerate(mqs)]
 end

src/dsl.jl

@@ -324,10 +324,10 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     if haskey(options, :ivs)
         ivs = Tuple(extract_syms(options, :ivs))
         ivexpr = copy(options[:ivs])
-        ivexpr.args[1] = Symbol("@", "variables")
+        ivexpr.args[1] = Symbol("@", "parameters")
     else
         ivs = (DEFAULT_IV_SYM,)
-        ivexpr = :(@variables $(DEFAULT_IV_SYM))
+        ivexpr = :($(DEFAULT_IV_SYM) = default_t())
     end
     tiv = ivs[1]
     sivs = (length(ivs) > 1) ? Expr(:vect, ivs[2:end]...) : nothing
@@ -877,7 +877,7 @@ function make_reaction(ex::Expr)
     sexprs = get_sexpr(species, Dict{Symbol, Expr}())
     pexprs = get_pexpr(parameters, Dict{Symbol, Expr}())
     rxexpr = get_rxexprs(reaction)
-    iv = :(@variables $(DEFAULT_IV_SYM))
+    iv = :($(DEFAULT_IV_SYM) = default_t())
 
     # Returns the rephrased expression.
     quote

src/reactionsystem_serialisation/serialise_fields.jl

@@ -8,7 +8,7 @@ end
 # Extract a string which declares the system's independent variable.
 function get_iv_string(rn::ReactionSystem)
     iv_dec = MT.get_iv(rn)
-    return "@variables $(iv_dec)"
+    return "@parameters $(iv_dec)"
 end
 
 # Creates an annotation for the system's independent variable.
@@ -92,12 +92,12 @@ function handle_us_n_ps(file_text::String, rn::ReactionSystem, annotate::Bool,
 
         # Pre-declares the sets with written/remaining parameters/species/variables.
         # Whenever all/none are written depends on whether there were any initial dependencies.
-        # `deepcopy` is required as these get mutated by `dependency_split!`. 
+        # `deepcopy` is required as these get mutated by `dependency_split!`.
         remaining_ps = (p_deps ? deepcopy(ps_all) : [])
         remaining_sps = (sp_deps ? deepcopy(sps_all) : [])
         remaining_vars = (var_deps ? deepcopy(vars_all) : [])
 
-        # Iteratively loops through all parameters, species, and/or variables. In each iteration, 
+        # Iteratively loops through all parameters, species, and/or variables. In each iteration,
         # adds the declaration of those that can still be declared.
         while !(isempty(remaining_ps) && isempty(remaining_sps) && isempty(remaining_vars))
             # Checks which parameters/species/variables can be written. The `dependency_split`
@@ -145,7 +145,7 @@ function seri_has_parameters(rn::ReactionSystem)
     return !isempty(get_ps(rn))
 end
 
-# Extract a string which declares the system's parameters. Uses multiline declaration (a 
+# Extract a string which declares the system's parameters. Uses multiline declaration (a
 # `begin ... end` block) if more than 3 parameters have a "complicated" declaration (if they
 # have metadata, default value, or type designation).
 function get_parameters_string(ps)
@@ -167,7 +167,7 @@ function seri_has_species(rn::ReactionSystem)
     return !isempty(get_species(rn))
 end
 
-# Extract a string which declares the system's species. Uses multiline declaration (a 
+# Extract a string which declares the system's species. Uses multiline declaration (a
 # `begin ... end` block) if more than 3 species have a "complicated" declaration (if they
 # have metadata, default value, or type designation).
 function get_species_string(sps)
@@ -189,7 +189,7 @@ function seri_has_variables(rn::ReactionSystem)
     return length(get_unknowns(rn)) > length(get_species(rn))
 end
 
-# Extract a string which declares the system's variables. Uses multiline declaration (a 
+# Extract a string which declares the system's variables. Uses multiline declaration (a
 # `begin ... end` block) if more than 3 variables have a "complicated" declaration (if they
 # have metadata, default value, or type designation).
 function get_variables_string(vars)
@@ -222,7 +222,7 @@ function get_reactions_string(rn::ReactionSystem)
         return "rxs = [$(reaction_string(get_rxs(rn)[1], strip_call_dict))]"
     end
 
-    # Creates the string corresponding to the code which generates the system's reactions. 
+    # Creates the string corresponding to the code which generates the system's reactions.
     rxs_string = "rxs = ["
     for rx in get_rxs(rn)
         @string_append! rxs_string "\n\t"*reaction_string(rx, strip_call_dict) ","
@@ -286,7 +286,7 @@ function get_equations_string(rn::ReactionSystem)
         return "eqs = [$(expression_2_string(get_eqs(rn)[end]; strip_call_dict))]"
     end
 
-    # Creates the string corresponding to the code which generates the system's reactions. 
+    # Creates the string corresponding to the code which generates the system's reactions.
     eqs_string = "eqs = ["
     for eq in get_eqs(rn)[(length(get_rxs(rn)) + 1):end]
         @string_append! eqs_string "\n\t" expression_2_string(eq; strip_call_dict) ","
@@ -392,7 +392,7 @@ function get_continuous_events_string(rn::ReactionSystem)
         return "continuous_events = [$(continuous_event_string(MT.get_continuous_events(rn)[1], strip_call_dict))]"
     end
 
-    # Creates the string corresponding to the code which generates the system's reactions. 
+    # Creates the string corresponding to the code which generates the system's reactions.
     continuous_events_string = "continuous_events = ["
     for continuous_event in MT.get_continuous_events(rn)
         @string_append! continuous_events_string "\n\t" continuous_event_string(
@@ -450,7 +450,7 @@ function get_discrete_events_string(rn::ReactionSystem)
         return "discrete_events = [$(discrete_event_string(MT.get_discrete_events(rn)[1], strip_call_dict))]"
     end
 
-    # Creates the string corresponding to the code which generates the system's reactions. 
+    # Creates the string corresponding to the code which generates the system's reactions.
     discrete_events_string = "discrete_events = ["
     for discrete_event in MT.get_discrete_events(rn)
         @string_append! discrete_events_string "\n\t" discrete_event_string(
@@ -493,7 +493,7 @@ DISCRETE_EVENTS_FS = (seri_has_discrete_events, get_discrete_events_string,
 ### Handles Systems ###
 
 # Specific `push_field` function, which is used for the system field (where the annotation option
-# must be passed to the `get_component_string` function). Since non-ReactionSystem systems cannot be 
+# must be passed to the `get_component_string` function). Since non-ReactionSystem systems cannot be
 # written to file, this function throws an error if any such systems are encountered.
 function push_systems_field(file_text::String, rn::ReactionSystem, annotate::Bool,
         top_level::Bool)

src/spatial_reaction_systems/spatial_reactions.jl

@@ -59,7 +59,7 @@ function make_transport_reaction(rateex, species)
     # Creates expressions corresponding to actual code from the internal DSL representation.
     sexprs = get_sexpr([species], Dict{Symbol, Expr}())
     pexprs = get_pexpr(parameters, Dict{Symbol, Expr}())
-    iv = :(@variables $(DEFAULT_IV_SYM))
+    iv = :($(DEFAULT_IV_SYM) = default_t())
     trxexpr = :(TransportReaction($rateex, $species))
 
     # Appends `edgeparameter` metadata to all declared parameters.
@@ -86,12 +86,12 @@ function check_spatial_reaction_validity(rs::ReactionSystem, tr::TransportReacti
         edge_parameters = [])
     # Checks that the species exist in the reaction system.
     # (ODE simulation code becomes difficult if this is not required,
-    # as non-spatial jacobian and f function generated from rs are of the wrong size).  
+    # as non-spatial jacobian and f function generated from rs are of the wrong size).
     if !any(isequal(tr.species), species(rs))
         error("Currently, species used in TransportReactions must have previously been declared within the non-spatial ReactionSystem. This is not the case for $(tr.species).")
     end
 
-    # Checks that the rate does not depend on species.    
+    # Checks that the rate does not depend on species.
     rate_vars = ModelingToolkit.getname.(Symbolics.get_variables(tr.rate))
     if !isempty(intersect(ModelingToolkit.getname.(species(rs)), rate_vars))
         error("The following species were used in rates of a transport reactions: $(setdiff(ModelingToolkit.getname.(species(rs)), rate_vars)).")

test/dsl/dsl_advanced_model_construction.jl

@@ -154,7 +154,7 @@ end
 # Tests reaction metadata accessor functions.
 let
     # Creates reactions directly.
-    @variables t
+    t = default_t()
     @parameters k η
     @species X(t) X2(t)
 
@@ -195,7 +195,7 @@ let
     @test Catalyst.hasmetadata(rxs[2], :md_1)
     @test !Catalyst.hasmetadata(rxs[3], :noise_scaling)
     @test !Catalyst.hasmetadata(rxs[3], :md_1)
-    
+
     @test isequal(Catalyst.getmetadata(rxs[1], :noise_scaling), η)
     @test isequal(Catalyst.getmetadata(rxs[2], :md_1), 1.0)
 
@@ -211,11 +211,11 @@ let
 end
 
 # Checks that repeated metadata throws errors.
-let 
-    @test_throws LoadError @eval @reaction k, 0 --> X, [md1=1.0, md1=2.0] 
+let
+    @test_throws LoadError @eval @reaction k, 0 --> X, [md1=1.0, md1=2.0]
     @test_throws LoadError @eval @reaction_network begin
-        k, 0 --> X, [md1=1.0, md1=1.0] 
-    end 
+        k, 0 --> X, [md1=1.0, md1=1.0]
+    end
 end
 
 # Tests for nested metadata.
@@ -233,14 +233,14 @@ let
         k8, Y7 --> X7, [md7="Hello"]
         k8, Y8 --> X7, [md7="Hi",md8="Yo"]
     end
-    
+
     rn2 = @reaction_network reactions begin
         (k1,k2,k3), (X1,X2,X3) --> (Y1,Y2,Y3), ([md1=1.0,md2=2.0],[md1=0.0],[md3="Hello world"])
         k, (X4,X5) --> (Y4,Y5), [md4=:sym]
         (k6, k6), X6 <--> Y6, ([md6=0.0],[md6=2.0])
         (k7,k8), X7 <--> (Y7,Y8), ([md7="Hi"],([md7="Hello"],[md7="Hi",md8="Yo"]))
     end
-    
+
     @test isequal(rn1, rn2)
 end
 
@@ -254,7 +254,7 @@ let
         k5, 3X5 --> Z5, [unnecessary_metadata=[1,2,3]]
         k6, X6 => Z6, [unnecessary_metadata=true]
     end
-    
+
     rn2 = @reaction_network reactions begin
         k1*X1, X1 + 2Y1 --> Z1, [only_use_rate=false]
         k2, 4X2 --> Z2 + W3, [only_use_rate=true]
@@ -263,7 +263,7 @@ let
         k5, 3X5 --> Z5, [only_use_rate=false, unnecessary_metadata=[1,2,3]]
         k6, X6 --> Z6, [only_use_rate=true, unnecessary_metadata=true]
     end
-    
+
     @test isequal(rn1,rn2)
 end