diff --git a/docs/src/api.md b/docs/src/api.md index f78a5f70e4..6efabcd08b 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -86,9 +86,9 @@ ReactionSystem ``` ## [Options for the `@reaction_network` DSL](@id api_dsl_options) -We have [previously described](@ref dsl_advanced_options) how options permits the user to supply non-reaction information to [`ReactionSystem`](@ref) created through the DSL. Here follows a list +We have [previously described](@ref dsl_advanced_options) how options permit the user to supply non-reaction information to [`ReactionSystem`](@ref) created through the DSL. Here follows a list of all options currently available. -- [`parameters`]:(@ref dsl_advanced_options_declaring_species_and_parameters) Allows the designation of a set of symbols as system parameter. +- [`parameters`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system parameters. - [`species`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system species. - [`variables`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system non-species variables. - [`ivs`](@ref dsl_advanced_options_ivs): Allows the designation of a set of symbols as system independent variables. @@ -100,6 +100,7 @@ of all options currently available. - [`continuous_events`](@ref constraint_equations_events): Allows the creation of continuous events. - [`discrete_events`](@ref constraint_equations_events): Allows the creation of discrete events. - [`combinatoric_ratelaws`](@ref faq_combinatoric_ratelaws): Takes a single option (`true` or `false`), which sets whether to use combinatorial rate laws. +- [`require_declaration`](@ref dsl_advanced_options_require_dec): Turns off all inference of parameters, species, variables, the default differential, and observables (requiring these to be explicitly declared using e.g. `@species`). ## [ModelingToolkit and Catalyst accessor functions](@id api_accessor_functions) A [`ReactionSystem`](@ref) is an instance of a diff --git a/docs/src/model_creation/dsl_advanced.md b/docs/src/model_creation/dsl_advanced.md index af806864b6..bac449cf94 100644 --- a/docs/src/model_creation/dsl_advanced.md +++ b/docs/src/model_creation/dsl_advanced.md @@ -204,7 +204,6 @@ ModelingToolkit.getdescription(two_state_system.kA) ``` ### [Designating constant-valued/fixed species parameters](@id dsl_advanced_options_constant_species) - Catalyst enables the designation of parameters as `constantspecies`. These parameters can be used as species in reactions, however, their values are not changed by the reaction and remain constant throughout the simulation (unless changed by e.g. the [occurrence of an event](@ref constraint_equations_events). Practically, this is done by setting the parameter's `isconstantspecies` metadata to `true`. Here, we create a simple reaction where the species `X` is converted to `Xᴾ` at rate `k`. By designating `X` as a constant species parameter, we ensure that its quantity is unchanged by the occurrence of the reaction. ```@example dsl_advanced_constant_species using Catalyst # hide @@ -575,6 +574,45 @@ nothing # hide !!! note When using interpolation, expressions like `2$spec` won't work; the multiplication symbol must be explicitly included like `2*$spec`. +## [Creating individual reaction using the `@reaction` macro](@id dsl_advanced_options_reaction_macro) +Catalyst exports a macro `@reaction`, which can be used to generate a singular [`Reaction`](@ref) object of the same type which is stored within the [`ReactionSystem`](@ref) structure (which in turn can be generated by `@reaction_network`). Generally, `@reaction` follows [identical rules to those of `@reaction_network`](@ref dsl_description_reactions) for writing and interpreting reactions (however, bi-directional reactions are not permitted). E.g. here we create a simple dimerisation reaction: +```@example dsl_advanced_reaction_macro +using Catalyst # hide +rx_dimerisation = @reaction kD, 2X --> X2 +``` +Here, `@reaction` is followed by a single line consisting of three parts: +- A rate (at which the reaction occurs). +- Any number of substrates (which are consumed by the reaction). +- Any number of products (which are produced by the reaction). + +In the next example, we first create a simple [SIR model](@ref basic_CRN_library_sir). Next, we instead create its individual reaction components using the `@reaction` macro. Finally, we confirm that these are identical to those stored in the initial model (using the [`reactions`](@ref) function). +```@example dsl_advanced_reaction_macro +sir_model = @reaction_network begin + α, S + I --> 2I + β, I --> R +end +infection_rx = @reaction α, S + I --> 2I +recovery_rx = @reaction β, I --> R +sir_rxs = [infection_rx, recovery_rx] +issetequal(reactions(sir_model), sir_rxs) +``` +One of the primary uses of the `@reaction` macro is to provide some of the convenience of the DSL to [*programmatic modelling](@ref programmatic_CRN_construction). E.g. here we can combine our reactions to create a `ReactionSystem` directly, and also confirm that this is identical to the model created through the DSL: +```@example dsl_advanced_reaction_macro +sir_programmatic = complete(ReactionSystem(sir_rxs, default_t(); name = :sir)) +sir_programmatic == sir_model +``` + +During programmatic modelling, it can be good to keep in mind that already declared symbolic variables can be [*interpolated*](@ref dsl_advanced_options_symbolics_and_DSL_interpolation). E.g. here we create two production reactions both depending on the same Michaelis-Menten function: +```@example dsl_advanced_reaction_macro +t = default_t() +@species X(t) +@parameters v K +mm_term = Catalyst.mm(X, v, K) +rx1 = @reaction $mm_term, 0 --> P1 +rx2 = @reaction $mm_term, 0 --> P2 +nothing # hide +``` + ## [Disabling mass action for reactions](@id dsl_advanced_options_disable_ma) As [described previously](@ref math_models_in_catalyst_rre_odes), Catalyst uses *mass action kinetics* to generate ODEs from reactions. Here, each reaction generates a term for each of its reactants, which consists of the reaction's rate, substrates, and the reactant's stoichiometry. E.g. the following reaction: diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 52daa0cec5..7c0cbbc4ee 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -72,11 +72,6 @@ const CONSERVED_CONSTANT_SYMBOL = :Γ const forbidden_symbols_skip = Set([:ℯ, :pi, :π, :t, :∅]) const forbidden_symbols_error = union(Set([:im, :nothing, CONSERVED_CONSTANT_SYMBOL]), forbidden_symbols_skip) -const forbidden_variables_error = let - fvars = copy(forbidden_symbols_error) - delete!(fvars, :t) - fvars -end ### Package Main ### diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 548ce06b4a..0576004691 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -86,7 +86,7 @@ function make_compound(expr) # 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)) + Vector{DSLReactant}(undef, 0)) components = :([]) # Becomes something like :([C, O]). coefficients = :([]) # Becomes something like :([1, 2]). for comp in composition diff --git a/src/dsl.jl b/src/dsl.jl index ac1bb88de9..902984668a 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -1,64 +1,4 @@ -""" -Macro that inputs an expression corresponding to a reaction network and outputs -a `ReactionNetwork` that can be used as input to generation of ODE, SDE, and -Jump problems. - -Most arrows accepted (both right, left, and bi-drectional arrows). Note that -while --> is a correct arrow, neither <-- nor <--> works. Using non-filled -arrows (⇐, ⟽, ⇒, ⟾, ⇔, ⟺) will disable mass kinetics and let you cutomize -reaction rates yourself. Use 0 or ∅ for degradation/creation to/from nothing. - -Example systems: - - ### Basic Usage ### - rn = @reaction_network begin # Creates a ReactionSystem. - 2.0, X + Y --> XY # This will have reaction rate corresponding to 2.0*[X][Y] - 2.0, XY ← X + Y # Identical to 2.0, X + Y --> XY - end - - ### Manipulating Reaction Rates ### - rn = @reaction_network begin - 2.0, X + Y ⟾ XY # Ignores mass kinetics. This will have reaction rate corresponding to 2.0. - 2.0X, X + Y --> XY # Reaction rate needs not be constant. This will have reaction rate corresponding to 2.0*[X]*[X]*[Y]. - XY+log(X)^2, X + Y --> XY # Reaction rate accepts quite complicated expressions. - hill(XY,2,2,2), X + Y --> XY # Reaction inis activated by XY according to a hill function. hill(x,v,K,N). - mm(XY,2,2), X + Y --> XY # Reaction inis activated by XY according to a michaelis menten function. mm(x,v,K). - end - - ### Multiple Reactions on a Single Line ### - rn = @reaction_network begin - (2.0,1.0), X + Y ↔ XY # Identical to reactions (2.0, X + Y --> XY) and (1.0, XY --> X + Y). - 2.0, (X,Y) --> 0 # This corresponds to both X and Y degrading at rate 2.0. - (2.0, 1.0), (X,Y) --> 0 # This corresponds to X and Y degrading at rates 2.0 and 1.0, respectively. - 2.0, (X1,Y1) --> (X2,Y2) # X1 and Y1 becomes X2 and Y2, respectively, at rate 2.0. - end - - ### Adding Parameters ### - kB = 2.0; kD = 1.0 - p = [kB, kD] - p = [] - rn = @reaction_network begin - (kB, kD), X + Y ↔ XY # Lets you define parameters outside on network. Parameters can be changed without recalling the network. - end - - ### Defining New Functions ### - my_hill_repression(x, v, k, n) = v*k^n/(k^n+x^n) - - # may be necessary to - # @register_symbolic my_hill_repression(x, v, k, n) - # see https://docs.sciml.ai/ModelingToolkit/stable/basics/Validation/#User-Defined-Registered-Functions-and-Types - - r = @reaction_network MyReactionType begin - my_hill_repression(x, v_x, k_x, n_x), 0 --> x - end - - ### Simulating Reaction Networks ### - probODE = ODEProblem(rn, args...; kwargs...) # Using multiple dispatch the reaction network can be used as input to create ODE, SDE and Jump problems. - probSDE = SDEProblem(rn, args...; kwargs...) - probJump = JumpProblem(prob,aggregator::Direct,rn) -""" - -### Constant Declarations ### +### Constants Declarations ### # Declare various arrow types symbols used for the empty set (also 0). const empty_set = Set{Symbol}([:∅]) @@ -110,159 +50,187 @@ end """ @reaction_network -Generates a [`ReactionSystem`](@ref dsl_description) that encodes a chemical reaction -network. +Macro for generating chemical reaction network models (Catalyst `ReactionSystem`s). See the +following two section ([DSL introduction](https://docs.sciml.ai/Catalyst/stable/model_creation/dsl_basics/) +and [advantage usage](https://docs.sciml.ai/Catalyst/stable/model_creation/dsl_advanced/)) of +the Catalyst documentation for more details on the domain-specific language (DSL) that the +macro implements. The macro's output (a `ReactionSystem` structure) is central to Catalyst +and its functionality. How to e.g. simulate these is described in the [Catalyst documentation](https://docs.sciml.ai/Catalyst/stable/). -See [The Reaction DSL](@ref dsl_description) documentation for details on -parameters to the macro. +Returns: +- A Catalyst `ReactionSystem`, i.e. a symbolic model for the reaction network. The returned +system is marked `complete`. To obtain a `ReactionSystem` that is not marked complete, for +example to then use in compositional modelling, see the otherwise equivalent `@network_component` +macro. Examples: +Here we create a basic SIR model. It contains two reactions (infection and recovery): ```julia -# a basic SIR model, with name SIR -sir_model = @reaction_network SIR begin - c1, s + i --> 2i - c2, i --> r -end - -# a basic SIR model, with random generated name sir_model = @reaction_network begin - c1, s + i --> 2i - c2, i --> r + c1, S + I --> 2I + c2, I --> R end +``` -# an empty network with name empty -emptyrn = @reaction_network empty - -# an empty network with random generated name -emptyrn = @reaction_network +Next, we create a self-activation loop. Here, a single component (`X`) activates its own production +with a Michaelis-Menten function: +```julia +sa_loop = @reaction_network begin + mm(X,v,K), 0 --> X + d, X --> 0 +end ``` +This model also contains production and degradation reactions, where `0` denotes that there are +either no substrates or no products in a reaction. -ReactionSystems generated through `@reaction_network` are complete. +Options: +In addition to reactions, the macro also supports "option" inputs (permitting e.g. the addition +of observables). Each option is designated by a tag starting with a `@` followed by its input. +A list of options can be found [here](https://docs.sciml.ai/Catalyst/stable/api/#api_dsl_options). """ -macro reaction_network(name::Symbol, ex::Expr) - :(complete($(make_reaction_system( - striplines(ex); name = :($(QuoteNode(name))))))) +macro reaction_network(name::Symbol, network_expr::Expr) + make_rs_expr(QuoteNode(name), network_expr) end -# Allows @reaction_network $name begin ... to interpolate variables storing a name. -macro reaction_network(name::Expr, ex::Expr) - :(complete($(make_reaction_system( - striplines(ex); name = :($(esc(name.args[1]))))))) +# The case where the name contains an interpolation. +macro reaction_network(name::Expr, network_expr::Expr) + make_rs_expr(esc(name.args[1]), network_expr) end -macro reaction_network(ex::Expr) - ex = striplines(ex) - - # no name but equations: @reaction_network begin ... end ... - if ex.head == :block - :(complete($(make_reaction_system(ex)))) - else # empty but has interpolated name: @reaction_network $name - networkname = :($(esc(ex.args[1]))) - return Expr(:block, :(@parameters t), - :(complete(ReactionSystem(Reaction[], t, [], []; name = $networkname)))) - end +# The case where nothing, or only a name, is provided. +macro reaction_network(name::Symbol = gensym(:ReactionSystem)) + make_rs_expr(QuoteNode(name)) end -# Returns a empty network (with, or without, a declared name). -macro reaction_network(name::Symbol = gensym(:ReactionSystem)) - return Expr(:block, :(@parameters t), - :(complete(ReactionSystem(Reaction[], t, [], []; name = $(QuoteNode(name)))))) +# Handles two disjoint cases. +macro reaction_network(expr::Expr) + # Case 1: The input is a name with interpolation. + (expr.head != :block) && return make_rs_expr(esc(expr.args[1])) + # Case 2: The input is a reaction network (and no name is provided). + return make_rs_expr(:(gensym(:ReactionSystem)), expr) end -# Ideally, it would have been possible to combine the @reaction_network and @network_component macros. -# However, this issue: https://github.com/JuliaLang/julia/issues/37691 causes problem with interpolations -# if we make the @reaction_network macro call the @network_component macro. +# Ideally, it would have been possible to combine the @reaction_network and @network_component +# macros. However, this issue: https://github.com/JuliaLang/julia/issues/37691 causes problem +# with interpolations if we make the @reaction_network macro call the @network_component macro. +# Instead, these use the same input, but pass `complete = false` to `make_rs_expr`. """ @network_component Equivalent to `@reaction_network` except the generated `ReactionSystem` is not marked as complete. """ -macro network_component(name::Symbol, ex::Expr) - make_reaction_system(striplines(ex); name = :($(QuoteNode(name)))) +macro network_component(name::Symbol, network_expr::Expr) + make_rs_expr(QuoteNode(name), network_expr; complete = false) end -# Allows @network_component $name begin ... to interpolate variables storing a name. -macro network_component(name::Expr, ex::Expr) - make_reaction_system(striplines(ex); name = :($(esc(name.args[1])))) +# The case where the name contains an interpolation. +macro network_component(name::Expr, network_expr::Expr) + make_rs_expr(esc(name.args[1]), network_expr; complete = false) end -macro network_component(ex::Expr) - ex = striplines(ex) +# The case where nothing, or only a name, is provided. +macro network_component(name::Symbol = gensym(:ReactionSystem)) + make_rs_expr(QuoteNode(name); complete = false) +end - # no name but equations: @network_component begin ... end ... - if ex.head == :block - make_reaction_system(ex) - else # empty but has interpolated name: @network_component $name - networkname = :($(esc(ex.args[1]))) - return Expr(:block, :(@parameters t), - :(ReactionSystem(Reaction[], t, [], []; name = $networkname))) - end +# Handles two disjoint cases. +macro network_component(expr::Expr) + # Case 1: The input is a name with interpolation. + (expr.head != :block) && return make_rs_expr(esc(expr.args[1]); complete = false) + # Case 2: The input is a reaction network (and no name is provided). + return make_rs_expr(:(gensym(:ReactionSystem)), expr; complete = false) end -# Returns a empty network (with, or without, a declared name). -macro network_component(name::Symbol = gensym(:ReactionSystem)) - return Expr(:block, :(@parameters t), - :(ReactionSystem(Reaction[], t, [], []; name = $(QuoteNode(name))))) +### DSL Macros Helper Functions ### + +# For when no reaction network has been designated. Generates an empty network. +function make_rs_expr(name; complete = true) + rs_expr = :(ReactionSystem(Reaction[], t, [], []; name = $name)) + complete && (rs_expr = :(complete($rs_expr))) + return Expr(:block, :(@parameters t), rs_expr) +end + +# When both a name and a network expression are generated, dispatch these to the internal +# `make_reaction_system` function. +function make_rs_expr(name, network_expr; complete = true) + rs_expr = make_reaction_system(striplines(network_expr), name) + complete && (rs_expr = :(complete($rs_expr))) + return rs_expr end ### Internal DSL Structures ### -# Structure containing information about one reactant in one reaction. -struct ReactantStruct +# Internal structure containing information about one reactant in one reaction. +struct DSLReactant reactant::Union{Symbol, Expr} stoichiometry::ExprValues end -# Structure containing information about one Reaction. Contain all its substrates and products as well as its rate. Contains a specialized constructor. -struct ReactionStruct - substrates::Vector{ReactantStruct} - products::Vector{ReactantStruct} +# Internal structure containing information about one Reaction. Contains all its substrates and +# products as well as its rate and potential metadata. Uses a specialized constructor. +struct DSLReaction + substrates::Vector{DSLReactant} + products::Vector{DSLReactant} rate::ExprValues metadata::Expr rxexpr::Expr - function ReactionStruct(sub_line::ExprValues, prod_line::ExprValues, rate::ExprValues, - metadata_line::ExprValues, rx_line::Expr) - sub = recursive_find_reactants!(sub_line, 1, Vector{ReactantStruct}(undef, 0)) - prod = recursive_find_reactants!(prod_line, 1, Vector{ReactantStruct}(undef, 0)) + function DSLReaction(sub_line::ExprValues, prod_line::ExprValues, + rate::ExprValues, metadata_line::ExprValues, rx_line::Expr) + subs = recursive_find_reactants!(sub_line, 1, Vector{DSLReactant}(undef, 0)) + prods = recursive_find_reactants!(prod_line, 1, Vector{DSLReactant}(undef, 0)) metadata = extract_metadata(metadata_line) - new(sub, prod, rate, metadata, rx_line) + new(subs, prods, rate, metadata, rx_line) end end # Recursive function that loops through the reaction line and finds the reactants and their -# stoichiometry. Recursion makes it able to handle weird cases like 2(X+Y+3(Z+XY)). +# stoichiometry. Recursion makes it able to handle weird cases like 2(X + Y + 3(Z + XY)). The +# reactants are stored in the `reactants` vector. As the expression tree is parsed, the +# stoichiometry is updated and new reactants are added. function recursive_find_reactants!(ex::ExprValues, mult::ExprValues, - reactants::Vector{ReactantStruct}) - if typeof(ex) != Expr || (ex.head == :escape) || (ex.head == :ref) + reactants::Vector{DSLReactant}) + # We have reached the end of the expression tree and can finalise and return the reactants. + if (typeof(ex) != Expr) || (ex.head == :escape) || (ex.head == :ref) + # The final bit of the expression is not a relevant reactant, no additions are required. (ex == 0 || in(ex, empty_set)) && (return reactants) - if any(ex == reactant.reactant for reactant in reactants) - idx = findall(x -> x == ex, getfield.(reactants, :reactant))[1] - reactants[idx] = ReactantStruct(ex, - processmult(+, mult, - reactants[idx].stoichiometry)) + + # If the expression corresponds to a reactant on our list, increase its multiplicity. + idx = findfirst(r.reactant == ex for r in reactants) + if !isnothing(idx) + newmult = processmult(+, mult, reactants[idx].stoichiometry) + reactants[idx] = DSLReactant(ex, newmult) + + # If the expression corresponds to a new reactant, add it to the list. else - push!(reactants, ReactantStruct(ex, mult)) + push!(reactants, DSLReactant(ex, mult)) end + + # If we have encountered a multiplication (i.e. a stoichiometry and a set of reactants). elseif ex.args[1] == :* + # The normal case (e.g. 3*X or 3*(X+Y)). Update the current multiplicity and continue. if length(ex.args) == 3 - recursive_find_reactants!(ex.args[3], processmult(*, mult, ex.args[2]), - reactants) + newmult = processmult(*, mult, ex.args[2]) + recursive_find_reactants!(ex.args[3], newmult, reactants) + # More complicated cases (e.g. 2*3*X). Yes, `ex.args[1:(end - 1)]` should start at 1 (not 2). else newmult = processmult(*, mult, Expr(:call, ex.args[1:(end - 1)]...)) recursive_find_reactants!(ex.args[end], newmult, reactants) end + # If we have encountered a sum of different reactants, apply recursion on each. elseif ex.args[1] == :+ for i in 2:length(ex.args) recursive_find_reactants!(ex.args[i], mult, reactants) end else - throw("Malformed reaction, bad operator: $(ex.args[1]) found in stochiometry expression $ex.") + throw("Malformed reaction, bad operator: $(ex.args[1]) found in stoichiometry expression $ex.") end reactants end +# Helper function for updating the multiplicity throughout recursion (handles e.g. parametric +# stoichiometries). The `op` argument is an operation (e.g. `*`, but could also e.g. be `+`). function processmult(op, mult, stoich) if (mult isa Number) && (stoich isa Number) op(mult, stoich) @@ -271,7 +239,8 @@ function processmult(op, mult, stoich) end end -# Finds the metadata from a metadata expression (`[key=val, ...]`) and returns as a Vector{Pair{Symbol,ExprValues}}. +# Finds the metadata from a metadata expression (`[key=val, ...]`) and returns as a +# `Vector{Pair{Symbol,ExprValues}}``. function extract_metadata(metadata_line::Expr) metadata = :([]) for arg in metadata_line.args @@ -284,8 +253,7 @@ function extract_metadata(metadata_line::Expr) return metadata end - - +### Specialised @require_declaration Option Error ### struct UndeclaredSymbolicError <: Exception msg::String end @@ -298,149 +266,114 @@ end ### DSL Internal Master Function ### # Function for creating a ReactionSystem structure (used by the @reaction_network macro). -function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) +function make_reaction_system(ex::Expr, name) - # Handle interpolation of variables + # Handle interpolation of variables in the input. ex = esc_dollars!(ex) - # Read lines with reactions and options. + # Extract the lines with reactions, the lines with options, and the options. Check for input errors. reaction_lines = Expr[x for x in ex.args if x.head == :tuple] option_lines = Expr[x for x in ex.args if x.head == :macrocall] + options = Dict(Symbol(String(arg.args[1])[2:end]) => arg for arg in option_lines) + allunique(arg.args[1] for arg in option_lines) || + error("Some options where given multiple times.") + (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.") + any(!in(option_keys), keys(options)) && + error("The following unsupported options were used: $(filter(opt_in->!in(opt_in,option_keys), keys(options)))") - # Get macro options. - if length(unique(arg.args[1] for arg in option_lines)) < length(option_lines) - error("Some options were given multiple times.") - end - options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg, - option_lines)) - - # Reads options. - default_reaction_metadata = :([]) - check_default_noise_scaling!(default_reaction_metadata, options) - compound_expr, compound_species = read_compound_options(options) - continuous_events_expr = read_events_option(options, :continuous_events) - discrete_events_expr = read_events_option(options, :discrete_events) + # Read options that explicitly declare some symbol (e.g. `@species`). Compiles a list of + # all declared symbols and checks that there have been no double-declarations. + sps_declared = extract_syms(options, :species) + ps_declared = extract_syms(options, :parameters) + vs_declared = extract_syms(options, :variables) + tiv, sivs, ivs, ivsexpr = read_ivs_option(options) + cmpexpr_init, cmps_declared = read_compound_option(options) + diffsexpr, diffs_declared = read_differentials_option(options) + syms_declared = collect(Iterators.flatten((cmps_declared, sps_declared, ps_declared, + vs_declared, ivs, diffs_declared))) + if !allunique(syms_declared) + nonunique_syms = [s for s in syms_declared if count(x -> x == s, syms_declared) > 1] + error("The following symbols $(unique(nonunique_syms)) have explicitly been declared as multiple types of components (e.g. occur in at least two of the `@species`, `@parameters`, `@variables`, `@ivs`, `@compounds`, `@differentials`). This is not allowed.") + end + + # Reads the reactions and equation. From these, infer species, variables, and parameters. requiredec = haskey(options, :require_declaration) - - # Parses reactions, species, and parameters. reactions = get_reactions(reaction_lines) - species_declared = [extract_syms(options, :species); compound_species] - parameters_declared = extract_syms(options, :parameters) - variables_declared = extract_syms(options, :variables) - - # Handle independent variables - if haskey(options, :ivs) - ivs = Tuple(extract_syms(options, :ivs)) - ivexpr = copy(options[:ivs]) - ivexpr.args[1] = Symbol("@", "parameters") - else - ivs = (DEFAULT_IV_SYM,) - ivexpr = :($(DEFAULT_IV_SYM) = default_t()) - end - tiv = ivs[1] - sivs = (length(ivs) > 1) ? Expr(:vect, ivs[2:end]...) : nothing - all_ivs = (isnothing(sivs) ? [tiv] : [tiv; sivs.args]) - - if haskey(options, :combinatoric_ratelaws) - combinatoric_ratelaws = options[:combinatoric_ratelaws].args[end] - else - combinatoric_ratelaws = true - end - - # Collect species and parameters, including ones inferred from the reactions. - declared_syms = Set(Iterators.flatten((parameters_declared, species_declared, - variables_declared))) - species_extracted, parameters_extracted = extract_species_and_parameters!( - reactions, declared_syms; requiredec) - - # Reads equations (and infers potential variables). - # Excludes any parameters already extracted (if they also was a variable). - declared_syms = union(declared_syms, species_extracted) - vars_extracted, add_default_diff, equations = read_equations_options( - options, declared_syms, parameters_extracted; requiredec) - variables = vcat(variables_declared, vars_extracted) - parameters_extracted = setdiff(parameters_extracted, vars_extracted) - - # Creates the finalised parameter and species lists. - species = vcat(species_declared, species_extracted) - parameters = vcat(parameters_declared, parameters_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; requiredec) - - # 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.") - any(!in(opt_in, option_keys) for opt_in in keys(options)) && - error("The following unsupported options were used: $(filter(opt_in->!in(opt_in,option_keys), keys(options)))") - forbidden_symbol_check(union(species, parameters)) - forbidden_variable_check(variables) - unique_symbol_check(union(species, parameters, variables, ivs)) + sps_inferred, ps_pre_inferred = extract_sps_and_ps(reactions, syms_declared; requiredec) + vs_inferred, diffs_inferred, equations = read_equations_option!(diffsexpr, options, + union(syms_declared, sps_inferred), tiv; requiredec) + ps_inferred = setdiff(ps_pre_inferred, vs_inferred, diffs_inferred) + syms_inferred = union(sps_inferred, ps_inferred, vs_inferred, diffs_inferred) + all_syms = union(syms_declared, syms_inferred) + obsexpr, obs_eqs, obs_syms = read_observed_option(options, ivs, + union(sps_declared, vs_declared), all_syms; requiredec) + + # Read options not related to the declaration or inference of symbols. + continuous_events_expr = read_events_option(options, :continuous_events) + discrete_events_expr = read_events_option(options, :discrete_events) + default_reaction_metadata = read_default_noise_scaling_option(options) + combinatoric_ratelaws = read_combinatoric_ratelaws_option(options) # 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; iv_symbols = ivs) - pexprs = get_pexpr(parameters_extracted, options) - ps, pssym = assign_expr_to_var(!isempty(parameters), pexprs, "ps") - vars, varssym = assign_expr_to_var(!isempty(variables), vexprs, "vars"; - scalarize = true) - sps, spssym = assign_expr_to_var(!isempty(species), sexprs, "specs"; scalarize = true) - comps, compssym = assign_expr_to_var(!isempty(compound_species), compound_expr, - "comps"; scalarize = true) - rxexprs = :(CatalystEqType[]) - for reaction in reactions - push!(rxexprs.args, get_rxexprs(reaction)) - end - for equation in equations - equation = escape_equation_RHS!(equation) - push!(rxexprs.args, equation) - end - - # Output code corresponding to the reaction system. - quote - $ivexpr - $ps - $vars - $sps - $observed_vars - $comps - $diffexpr - - sivs_vec = $sivs - rx_eq_vec = $rxexprs - vars = setdiff(union($spssym, $varssym, $compssym), $obs_syms) - obseqs = $observed_eqs - cevents = $continuous_events_expr - devents = $discrete_events_expr + psexpr_init = get_psexpr(ps_inferred, options) + spsexpr_init = get_usexpr(sps_inferred, options; ivs) + vsexpr_init = get_usexpr(vs_inferred, options, :variables; ivs) + psexpr, psvar = assign_var_to_symvar_declaration(psexpr_init, "ps") + spsexpr, spsvar = assign_var_to_symvar_declaration(spsexpr_init, "specs") + vsexpr, vsvar = assign_var_to_symvar_declaration(vsexpr_init, "vars") + cmpsexpr, cmpsvar = assign_var_to_symvar_declaration(cmpexpr_init, "comps") + rxsexprs = get_rxexprs(reactions, equations, all_syms) + + # Assemblies the full expression that declares all required symbolic variables, and + # then the output `ReactionSystem`. + MacroTools.flatten(striplines(quote + # Inserts the expressions which generate the `ReactionSystem` input. + $ivsexpr + $psexpr + $spsexpr + $vsexpr + $obsexpr + $cmpsexpr + $diffsexpr + + # Stores each kwarg in a variable. Not necessary, but useful when debugging generated code. + name = $name + spatial_ivs = $sivs + rx_eq_vec = $rxsexprs + us = setdiff(union($spsvar, $vsvar, $cmpsvar), $obs_syms) + observed = $obs_eqs + continuous_events = $continuous_events_expr + discrete_events = $discrete_events_expr + combinatoric_ratelaws = $combinatoric_ratelaws + default_reaction_metadata = $default_reaction_metadata remake_ReactionSystem_internal( - make_ReactionSystem_internal( - rx_eq_vec, $tiv, vars, $pssym; - name = $name, spatial_ivs = sivs_vec, observed = obseqs, - continuous_events = cevents, discrete_events = devents, - combinatoric_ratelaws = $combinatoric_ratelaws); - default_reaction_metadata = $default_reaction_metadata) - end + make_ReactionSystem_internal(rx_eq_vec, $tiv, us, $psvar; name, spatial_ivs, + observed, continuous_events, discrete_events, combinatoric_ratelaws); + default_reaction_metadata) + end)) end ### DSL Reaction Reading Functions ### -# Generates a vector containing a number of reaction structures, each containing the information about one reaction. -function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(undef, 0)) +# Generates a vector of reaction structures, each containing information about one reaction. +function get_reactions(exprs::Vector{Expr}) + # Declares an array to which we add all found reactions. + reactions = Vector{DSLReaction}(undef, 0) + + # Loops through each line of reactions. Extracts and adds each line's reactions to `reactions`. for line in exprs # Reads core reaction information. arrow, rate, reaction, metadata = read_reaction_line(line) - # Checks the type of arrow used, and creates the corresponding reaction(s). Returns them in an array. + # Checks which type of line is used, and calls `push_reactions!` on the processed line. if in(arrow, double_arrows) - if typeof(rate) != Expr || rate.head != :tuple + (typeof(rate) != Expr || rate.head != :tuple) && error("Error: Must provide a tuple of reaction rates when declaring a bi-directional reaction.") - end + (typeof(metadata) != Expr || metadata.head != :tuple) && + error("Error: Must provide a tuple of reaction metadata when declaring a bi-directional reaction.") + push_reactions!(reactions, reaction.args[2], reaction.args[3], rate.args[1], metadata.args[1], arrow, line) push_reactions!(reactions, reaction.args[3], reaction.args[2], @@ -458,10 +391,10 @@ function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(u return reactions end -# Extracts the rate, reaction, and metadata fields (the last one optional) from a reaction line. +# Extract the rate, reaction, and metadata fields (the last one optional) from a reaction line. function read_reaction_line(line::Expr) - # Handles rate, reaction, and arrow. - # Special routine required for the`-->` case, which creates different expression from all other cases. + # Handles rate, reaction, and arrow. A special routine is required for the`-->` case, + # which creates an expression different from what the other arrows create. rate = line.args[1] reaction = line.args[2] if reaction.head == :--> @@ -481,22 +414,23 @@ function read_reaction_line(line::Expr) return arrow, rate, reaction, metadata end -# Takes a reaction line and creates reaction(s) from it and pushes those to the reaction array. -# Used to create multiple reactions from, for instance, `k, (X,Y) --> 0`. -function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, - prod_line::ExprValues, rate::ExprValues, metadata::ExprValues, arrow::Symbol, line::Expr) - # The rates, substrates, products, and metadata may be in a tupple form (e.g. `k, (X,Y) --> 0`). - # This finds the length of these tuples (or 1 if not in tuple forms). Errors if lengs inconsistent. - lengs = (tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate), tup_leng(metadata)) - if any(!(leng == 1 || leng == maximum(lengs)) for leng in lengs) - throw("Malformed reaction, rate=$rate, subs=$sub_line, prods=$prod_line, metadata=$metadata.") - end - - # Loops through each reaction encoded by the reaction composites. Adds the reaction to the reactions vector. - for i in 1:maximum(lengs) - # If the `only_use_rate` metadata was not provided, this has to be infered from the arrow used. +# Takes a reaction line and creates reaction(s) from it and pushes those to the reaction vector. +# Used to create multiple reactions from bundled reactions (like `k, (X,Y) --> 0`). +function push_reactions!(reactions::Vector{DSLReaction}, subs::ExprValues, + prods::ExprValues, rate::ExprValues, metadata::ExprValues, arrow::Symbol, line::Expr) + # The rates, substrates, products, and metadata may be in a tuple form (e.g. `k, (X,Y) --> 0`). + # This finds these tuples' lengths (or 1 for non-tuple forms). Inconsistent lengths yield error. + lengs = (tup_leng(subs), tup_leng(prods), tup_leng(rate), tup_leng(metadata)) + maxl = maximum(lengs) + any(!(leng == 1 || leng == maxl) for leng in lengs) && + error("Malformed reaction, rate: $rate, subs: $subs, prods: $prods, metadata: $metadata.") + + # Loops through each reaction encoded by the reaction's different components. + # Creates a `DSLReaction` representation and adds it to `reactions`. + for i in 1:maxl + # If the `only_use_rate` metadata was not provided, this must be inferred from the arrow. metadata_i = get_tup_arg(metadata, i) - if all(arg -> arg.args[1] != :only_use_rate, metadata_i.args) + if all(arg.args[1] != :only_use_rate for arg in metadata_i.args) push!(metadata_i.args, :(only_use_rate = $(in(arrow, pure_rate_arrows)))) end @@ -505,9 +439,9 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues error("Some reaction metadata fields where repeated: $(metadata_entries)") end - push!(reactions, - ReactionStruct(get_tup_arg(sub_line, i), - get_tup_arg(prod_line, i), get_tup_arg(rate, i), metadata_i, line)) + # Extracts substrates, products, and rates for the i'th reaction. + subs_i, prods_i, rate_i = get_tup_arg.((subs, prods, rate), i) + push!(reactions, DSLReaction(subs_i, prods_i, rate_i, metadata_i, line)) end end @@ -516,155 +450,197 @@ end # When the user have used the @species (or @parameters) options, extract species (or # parameters) from its input. function extract_syms(opts, vartype::Symbol) - if haskey(opts, vartype) + # If the corresponding option has been used, use `Symbolics._parse_vars` to find all + # variable within it (returning them in a vector). + return if haskey(opts, vartype) ex = opts[vartype] vars = Symbolics._parse_vars(vartype, Real, ex.args[3:end]) - syms = Vector{Union{Symbol, Expr}}(vars.args[end].args) + Vector{Union{Symbol, Expr}}(vars.args[end].args) else - syms = Union{Symbol, Expr}[] + Union{Symbol, Expr}[] end - return syms end # Function looping through all reactions, to find undeclared symbols (species or -# parameters), and assign them to the right category. -function extract_species_and_parameters!(reactions, excluded_syms; requiredec = false) +# parameters) and assign them to the right category. +function extract_sps_and_ps(reactions, excluded_syms; requiredec = false) + # Loops through all reactants and extract undeclared ones as species. species = OrderedSet{Union{Symbol, Expr}}() for reaction in reactions for reactant in Iterators.flatten((reaction.substrates, reaction.products)) add_syms_from_expr!(species, reactant.reactant, excluded_syms) - (!isempty(species) && requiredec) && throw(UndeclaredSymbolicError( - "Unrecognized variables $(join(species, ", ")) detected in reaction expression: \"$(string(reaction.rxexpr))\". Since the flag @require_declaration is declared, all species must be explicitly declared with the @species macro.")) end + (!isempty(species) && requiredec) && + throw(UndeclaredSymbolicError("Unrecognized reactant $(join(species, ", ")) detected in reaction expression: \"$(string(reaction.rxexpr))\". Since the flag @require_declaration is declared, all species must be explicitly declared with the @species option.")) end + union!(excluded_syms, species) - foreach(s -> push!(excluded_syms, s), species) + # Loops through all rates and stoichiometries, extracting used symbols as parameters. parameters = OrderedSet{Union{Symbol, Expr}}() for reaction in reactions add_syms_from_expr!(parameters, reaction.rate, excluded_syms) - (!isempty(parameters) && requiredec) && throw(UndeclaredSymbolicError( - "Unrecognized parameter $(join(parameters, ", ")) detected in rate expression: $(reaction.rate) for the following reaction expression: \"$(string(reaction.rxexpr))\". Since the flag @require_declaration is declared, all parameters must be explicitly declared with the @parameters macro.")) + (!isempty(parameters) && requiredec) && + throw(UndeclaredSymbolicError("Unrecognized symbol $(join(parameters, ", ")) detected in rate expression: $(reaction.rate) for the following reaction expression: \"$(string(reaction.rxexpr))\". Since the flag @require_declaration is declared, all parameters must be explicitly declared with the @parameters option.")) for reactant in Iterators.flatten((reaction.substrates, reaction.products)) add_syms_from_expr!(parameters, reactant.stoichiometry, excluded_syms) - (!isempty(parameters) && requiredec) && throw(UndeclaredSymbolicError( - "Unrecognized parameters $(join(parameters, ", ")) detected in the stoichiometry for reactant $(reactant.reactant) in the following reaction expression: \"$(string(reaction.rxexpr))\". Since the flag @require_declaration is declared, all parameters must be explicitly declared with the @parameters macro.")) + (!isempty(parameters) && requiredec) && + throw(UndeclaredSymbolicError("Unrecognized symbol $(join(parameters, ", ")) detected in the stoichiometry for reactant $(reactant.reactant) in the following reaction expression: \"$(string(reaction.rxexpr))\". Since the flag @require_declaration is declared, all parameters must be explicitly declared with the @parameters option.")) end end collect(species), collect(parameters) end -# Function called by extract_species_and_parameters!, recursively loops through an +# Function called by `extract_sps_and_ps`, recursively loops through an # expression and find symbols (adding them to the push_symbols vector). -function add_syms_from_expr!(push_symbols::AbstractSet, rateex::ExprValues, excluded_syms) - if rateex isa Symbol - if !(rateex in forbidden_symbols_skip) && !(rateex in excluded_syms) - push!(push_symbols, rateex) +function add_syms_from_expr!(push_symbols::AbstractSet, expr::ExprValues, excluded_syms) + # If we have encountered a Symbol in the recursion, we can try extracting it. + if expr isa Symbol + if !(expr in forbidden_symbols_skip) && !(expr in excluded_syms) + push!(push_symbols, expr) end - elseif rateex isa Expr + elseif expr isa Expr # note, this (correctly) skips $(...) expressions - for i in 2:length(rateex.args) - add_syms_from_expr!(push_symbols, rateex.args[i], excluded_syms) + for i in 2:length(expr.args) + add_syms_from_expr!(push_symbols, expr.args[i], excluded_syms) end end - nothing end ### DSL Output Expression Builders ### -# Given the species that were extracted from the reactions, and the options dictionary, creates the @species ... expression for the macro output. -function get_sexpr(species_extracted, options, key = :species; - iv_symbols = (DEFAULT_IV_SYM,)) - if haskey(options, key) - sexprs = options[key] - elseif isempty(species_extracted) - sexprs = :() +# Given the extracted species (or variables) and the option dictionary, create the +# `@species ...` (or `@variables ..`) expression which would declare these. +# If `key = :variables`, does this for variables (and not species). +function get_usexpr(us_extracted, options, key = :species; ivs = (DEFAULT_IV_SYM,)) + usexpr = if haskey(options, key) + options[key] + elseif isempty(us_extracted) + :() else - sexprs = Expr(:macrocall, Symbol("@", key), LineNumberNode(0)) + Expr(:macrocall, Symbol("@", key), LineNumberNode(0)) end - foreach(s -> (s isa Symbol) && push!(sexprs.args, Expr(:call, s, iv_symbols...)), - species_extracted) - sexprs + for u in us_extracted + u isa Symbol && push!(usexpr.args, Expr(:call, u, ivs...)) + end + usexpr end -# Given the parameters that were extracted from the reactions, and the options dictionary, creates the @parameters ... expression for the macro output. -function get_pexpr(parameters_extracted, options) - pexprs = (haskey(options, :parameters) ? options[:parameters] : - (isempty(parameters_extracted) ? :() : :(@parameters))) +# Given the parameters that were extracted from the reactions, and the options dictionary, +# creates the `@parameters ...` expression for the macro output. +function get_psexpr(parameters_extracted, options) + pexprs = if haskey(options, :parameters) + options[:parameters] + elseif isempty(parameters_extracted) + :() + else + :(@parameters) + end foreach(p -> push!(pexprs.args, p), parameters_extracted) pexprs end -# Creates the reaction vector declaration statement. -function get_rxexprs(rxstruct) - subs_init = isempty(rxstruct.substrates) ? nothing : :([]) +# From the system reactions (as `DSLReaction`s) and equations (as expressions), +# creates the expression that evaluates to the reaction (+ equations) vector. +function get_rxexprs(reactions, equations, all_syms) + rxsexprs = :(Catalyst.CatalystEqType[]) + foreach(rx -> push!(rxsexprs.args, get_rxexpr(rx)), reactions) + foreach(eq -> push!(rxsexprs.args, escape_equation!(eq, all_syms)), equations) + return rxsexprs +end + +# From a `DSLReaction` struct, creates the expression which evaluates to the creation +# of the corresponding reaction. +function get_rxexpr(rx::DSLReaction) + # Initiates the `Reaction` expression. + rate = recursive_escape_functions!(rx.rate) + subs_init = isempty(rx.substrates) ? nothing : :([]) subs_stoich_init = deepcopy(subs_init) - prod_init = isempty(rxstruct.products) ? nothing : :([]) + prod_init = isempty(rx.products) ? nothing : :([]) prod_stoich_init = deepcopy(prod_init) - reaction_func = :(Reaction($(recursive_escape_functions!(rxstruct.rate)), $subs_init, - $prod_init, $subs_stoich_init, $prod_stoich_init, - metadata = $(rxstruct.metadata))) - for sub in rxstruct.substrates - push!(reaction_func.args[3].args, sub.reactant) - push!(reaction_func.args[5].args, sub.stoichiometry) + rx_constructor = :(Reaction($rate, $subs_init, $prod_init, $subs_stoich_init, + $prod_stoich_init; metadata = $(rx.metadata))) + + # Loops through all products and substrates, and adds them (and their stoichiometries) + # to the `Reaction` expression. + for sub in rx.substrates + push!(rx_constructor.args[4].args, sub.reactant) + push!(rx_constructor.args[6].args, sub.stoichiometry) end - for prod in rxstruct.products - push!(reaction_func.args[4].args, prod.reactant) - push!(reaction_func.args[6].args, prod.stoichiometry) + for prod in rx.products + push!(rx_constructor.args[5].args, prod.reactant) + push!(rx_constructor.args[7].args, prod.stoichiometry) end - reaction_func + return rx_constructor end -# takes a ModelingToolkit declaration macro like @parameters and returns an expression -# that calls the macro and saves it in a variable given by namesym based on name. -# scalarizes if desired -function assign_expr_to_var(nonempty, ex, name; scalarize = false) +# Takes a ModelingToolkit declaration macro (like @parameters ...) and return and expression: +# That calls the macro and then scalarizes all the symbols created into a vector of Nums. +# stores the created symbolic variables in a variable (whose name is generated from `name`). +# It will also return the name used for the variable that stores the symbolic variables. +function assign_var_to_symvar_declaration(expr_init, name) + # Generates a random variable name which (in generated code) will store the produced + # symbolic variables (e.g. `var"##ps#384"`). namesym = gensym(name) - if nonempty - if scalarize - symvec = gensym() - ex = quote - $symvec = $ex - $namesym = reduce(vcat, Symbolics.scalarize($symvec)) - end - else - ex = quote - $namesym = $ex - end + + # If the input expression is non-empty, wrap it with additional information. + if expr_init != :(()) + symvec = gensym() + expr = quote + $symvec = $expr_init + $namesym = reduce(vcat, Symbolics.scalarize($symvec)) end else - ex = :($namesym = Num[]) + expr = :($namesym = Num[]) end - ex, namesym + + return expr, namesym +end + +# Recursively escape functions within equations of an equation written using user-defined functions. +# Does not escape special function calls like "hill(...)" and differential operators. Does +# also not escape stuff corresponding to e.g. species or parameters (required for good error +# for when e.g. a species is used as a differential, or for time delays in the future). +function escape_equation!(eqexpr::Expr, all_syms) + eqexpr.args[2] = recursive_escape_functions!(eqexpr.args[2], all_syms) + eqexpr.args[3] = recursive_escape_functions!(eqexpr.args[3], all_syms) + eqexpr end ### DSL Option Handling ### -# Checks if the `@default_noise_scaling` option is used. If so, read its input and adds it as a -# default metadata value to the `default_reaction_metadata` vector. -function check_default_noise_scaling!(default_reaction_metadata, options) +# Returns the `default_reaction_metadata` output. Technically Catalyst's code could have been made +# more generic to account for other default reaction metadata. Practically, this will likely +# be the only relevant reaction metadata to have a default value via the DSL. If another becomes +# relevant, the code can be rewritten to take this into account. +# Checks if the `@default_noise_scaling` option is used. If so, use it as the default value of +# the `default_noise_scaling` reaction metadata, otherwise, returns an empty vector. +function read_default_noise_scaling_option(options) if haskey(options, :default_noise_scaling) - if (length(options[:default_noise_scaling].args) != 3) # Because of how expressions are, 3 is used. - error("@default_noise_scaling should only have a single input, this appears not to be the case: \"$(options[:default_noise_scaling])\"") - end - push!(default_reaction_metadata.args, - :(:noise_scaling => $(options[:default_noise_scaling].args[3]))) - end -end - -# When compound species are declared using the "@compound begin ... end" option, get a list of the compound species, and also the expression that crates them. -function read_compound_options(opts) - # If the compound option is used retrieve a list of compound species (need to be added to the reaction system's species), and the option that creates them (used to declare them as compounds at the end). - if haskey(opts, :compounds) - compound_expr = opts[:compounds] - # Find compound species names, and append the independent variable. - compound_species = [find_varinfo_in_declaration(arg.args[2])[1] - for arg in compound_expr.args[3].args] + (length(options[:default_noise_scaling].args) != 3) && + error("@default_noise_scaling should only have a single expression as its input, this appears not to be the case: \"$(options[:default_noise_scaling])\"") + return :([:noise_scaling => $(options[:default_noise_scaling].args[3])]) + end + return :([]) +end + +# When compound species are declared using the "@compound begin ... end" option, get a list +# of the compound species, and also the expression that creates them. +function read_compound_option(options) + # If the compound option is used, retrieve a list of compound species and the option line + # that creates them (used to declare them as compounds at the end). Due to some expression + # handling, in the case of a single compound we must change to the `@compound` macro. + if haskey(options, :compounds) + cmpexpr_init = options[:compounds] + cmpexpr_init.args[3] = option_block_form(get_block_option(cmpexpr_init)) + cmps_declared = [find_varinfo_in_declaration(arg.args[2])[1] + for arg in cmpexpr_init.args[3].args] + (length(cmps_declared) == 1) && (cmpexpr_init.args[1] = Symbol("@compound")) else # If option is not used, return empty vectors and expressions. - compound_expr = :() - compound_species = Union{Symbol, Expr}[] + cmpexpr_init = :() + cmps_declared = Union{Symbol, Expr}[] end - return compound_expr, compound_species + return cmpexpr_init, cmps_declared end # Read the events (continuous or discrete) provided as options to the DSL. Returns an expression which evaluates to these. @@ -673,7 +649,7 @@ function read_events_option(options, event_type::Symbol) if event_type ∉ [:continuous_events, :discrete_events] error("Trying to read an unsupported event type.") end - events_input = haskey(options, event_type) ? options[event_type].args[3] : + events_input = haskey(options, event_type) ? get_block_option(options[event_type]) : striplines(:(begin end)) events_input = option_block_form(events_input) @@ -691,7 +667,7 @@ function read_events_option(options, event_type::Symbol) error("The condition part of continuous events (the left-hand side) must be a vector. This is not the case for: $(arg).") end if (arg isa Expr) && (arg.args[3] isa Expr) && (arg.args[3].head != :vect) - error("The affect part of all events (the righ-hand side) must be a vector. This is not the case for: $(arg).") + error("The affect part of all events (the right-hand side) must be a vector. This is not the case for: $(arg).") end # Adds the correctly formatted event to the event creation expression. @@ -701,14 +677,13 @@ function read_events_option(options, event_type::Symbol) return events_expr end -# Reads the variables options. Outputs: -# `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, syms_declared, parameters_extracted; requiredec = false) - # Prepares the equations. First, extracts equations from provided option (converting to block form if required). +# Reads the variables options. Outputs a list of the variables inferred from the equations, +# as well as the equation vector. If the default differential was used, update the `diffsexpr` +# expression so that this declares this as well. +function read_equations_option!(diffsexpr, options, syms_unavailable, tiv; requiredec = false) + # Prepares the equations. First, extract equations from the 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) + eqs_input = haskey(options, :equations) ? get_block_option(options[:equations]) : :(begin end) eqs_input = option_block_form(eqs_input) equations = Expr[] ModelingToolkit.parse_equations!(Expr(:block), equations, @@ -717,31 +692,38 @@ function read_equations_options(options, syms_declared, parameters_extracted; re # 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 = OrderedSet{Union{Symbol, Expr}}() + vs_inferred = OrderedSet{Union{Symbol, Expr}}() 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 \"~\".") + error("Malformed equation: \"$eq\". Equation's left hand and right-hand sides should be separated by a \"~\".") end - # If the default differential (`D`) is used, record that it should be decalred later on. - if (:D ∉ union(syms_declared, parameters_extracted)) && find_D_call(eq) + # If the default differential (`D`) is used, record that it should be declared later on. + if (:D ∉ syms_unavailable) && find_D_call(eq) requiredec && throw(UndeclaredSymbolicError( "Unrecognized symbol D was used as a differential in an equation: \"$eq\". Since the @require_declaration flag is set, all differentials in equations must be explicitly declared using the @differentials option.")) add_default_diff = true - push!(syms_declared, :D) + push!(syms_unavailable, :D) end - # Any undecalred symbolic variables encountered should be extracted as variables. - add_syms_from_expr!(vars_extracted, eq, syms_declared) - (!isempty(vars_extracted) && requiredec) && throw(UndeclaredSymbolicError( - "Unrecognized symbolic variables $(join(vars_extracted, ", ")) detected in equation expression: \"$(string(eq))\". Since the flag @require_declaration is declared, all symbolic variables must be explicitly declared with the @species, @variables, and @parameters options.")) + # Any undeclared symbolic variables encountered should be extracted as variables. + add_syms_from_expr!(vs_inferred, eq, syms_unavailable) + (!isempty(vs_inferred) && requiredec) && throw(UndeclaredSymbolicError( + "Unrecognized symbol $(join(vs_inferred, ", ")) detected in equation expression: \"$(string(eq))\". Since the flag @require_declaration is declared, all symbolic variables must be explicitly declared with the @species, @variables, and @parameters options.")) + end + + # If `D` differential is used, add it to the differential expression and inferred differentials list. + diffs_inferred = Union{Symbol, Expr}[] + if add_default_diff && !any(diff_dec.args[1] == :D for diff_dec in diffsexpr.args) + diffs_inferred = [:D] + push!(diffsexpr.args, :(D = Differential($(tiv)))) end - return collect(vars_extracted), add_default_diff, equations + return collect(vs_inferred), diffs_inferred, equations end -# Searches an expresion `expr` and returns true if it have any subexpression `D(...)` (where `...` can be anything). +# Searches an expression `expr` and returns true if it has any subexpression `D(...)` (where `...` can be anything). # Used to determine whether the default differential D has been used in any equation provided to `@equations`. function find_D_call(expr) return if Base.isexpr(expr, :call) && expr.args[1] == :D @@ -755,71 +737,66 @@ end # 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) +function read_differentials_option(options) # Creates the differential expression. - # If differentials was provided as options, this is used as the initial expression. + # If differentials were 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] : - 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. - for dexpr in diffexpr.args + diffsexpr = (haskey(options, :differentials) ? + get_block_option(options[:differentials]) : striplines(:(begin end))) + diffsexpr = option_block_form(diffsexpr) + + # Goes through all differentials, checking that they are correctly formatted. Adds their + # symbol to the list of declared differential symbols. + diffs_declared = Union{Symbol, Expr}[] + for dexpr in diffsexpr.args (dexpr.head != :(=)) && error("Differential declaration must have form like D = Differential(t), instead \"$(dexpr)\" was given.") (dexpr.args[1] isa Symbol) || error("Differential left-hand side must be a single symbol, instead \"$(dexpr.args[1])\" was given.") - in(dexpr.args[1], used_syms) && - error("Differential name ($(dexpr.args[1])) is also a species, variable, or parameter. This is ambiguous and not allowed.") in(dexpr.args[1], forbidden_symbols_error) && error("A forbidden symbol ($(dexpr.args[1])) was used as a differential name.") + push!(diffs_declared, dexpr.args[1]) end - # If the default differential D has been used, but not pre-declared using the @differentials - # 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($(tiv)))) - end - - return diffexpr + return diffsexpr, diffs_declared end -# Reads the observables options. Outputs an expression ofr creating the observable variables, and a vector of observable equations. -function read_observed_options(options, species_n_vars_declared, ivs_sorted; requiredec = false) +# Reads the observables options. Outputs an expression for creating the observable variables, +# a vector containing the observable equations, and a list of all observable symbols (this +# list contains both those declared separately or inferred from the `@observables` option` input`). +function read_observed_option(options, all_ivs, us_declared, all_syms; requiredec = false) + syms_unavailable = setdiff(all_syms, us_declared) if haskey(options, :observables) # Gets list of observable equations and prepares variable declaration expression. # (`options[:observables]` includes `@observables`, `.args[3]` removes this part) - observed_eqs = make_observed_eqs(options[:observables].args[3]) - observed_vars = Expr(:block, :(@variables)) + obs_eqs = make_obs_eqs(get_block_option(options[:observables])) + obsexpr = Expr(:block, :(@variables)) obs_syms = :([]) - for (idx, obs_eq) in enumerate(observed_eqs.args) - # Extract the observable, checks errors, and continues the loop if the observable has been declared. + for (idx, obs_eq) in enumerate(obs_eqs.args) + # Extract the observable, checks for errors. obs_name, ivs, defaults, metadata = find_varinfo_in_declaration(obs_eq.args[2]) - if (requiredec && !in(obs_name, species_n_vars_declared)) - throw(UndeclaredSymbolicError( - "An undeclared variable ($obs_name) was declared as an observable in the following observable equation: \"$obs_eq\". Since the flag @require_declaration is set, all variables must be declared with the @species, @parameters, or @variables macros.")) - end + + # Error checks. + (requiredec && !in(obs_name, us_declared)) && + throw(UndeclaredSymbolicError("An undeclared symbol ($obs_name) was used as an observable in the following observable equation: \"$obs_eq\". Since the flag @require_declaration is set, all observables must be declared with either the @species or @variables options.")) isempty(ivs) || error("An observable ($obs_name) was given independent variable(s). These should not be given, as they are inferred automatically.") isnothing(defaults) || error("An observable ($obs_name) was given a default value. This is forbidden.") - (obs_name in forbidden_symbols_error) && + in(obs_name, forbidden_symbols_error) && error("A forbidden symbol ($(obs_eq.args[2])) was used as an observable name.") - - # Error checks. - if (obs_name in species_n_vars_declared) && is_escaped_expr(obs_eq.args[2]) + in(obs_name, syms_unavailable) && + error("An observable ($obs_name) uses a name that already have been already been declared or inferred as another model property.") + (obs_name in us_declared) && is_escaped_expr(obs_eq.args[2]) && error("An interpolated observable have been used, which has also been explicitly declared within the system using either @species or @variables. This is not permitted.") - end - if ((obs_name in species_n_vars_declared) || is_escaped_expr(obs_eq.args[2])) && - !isnothing(metadata) + ((obs_name in us_declared) || is_escaped_expr(obs_eq.args[2])) && !isnothing(metadata) && error("Metadata was provided to observable $obs_name in the `@observables` macro. However, the observable was also declared separately (using either @species or @variables). When this is done, metadata should instead be provided within the original @species or @variable declaration.") - end - # This bits adds the observables to the @variables vector which is given as output. + # This bit adds the observables to the @variables vector which is given as output. # 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])) + if !((obs_name in us_declared) || is_escaped_expr(obs_eq.args[2])) # Creates an expression which extracts the ivs of the species & variables the # observable depends on, and splats them out in the correct order. dep_var_expr = :(filter(!MT.isparameter, @@ -828,15 +805,15 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted; req vcat, [sorted_arguments(MT.unwrap(dep)) 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))) + by = iv -> findfirst(MT.getname(iv) == ivs for ivs in $all_ivs))) obs_expr = insert_independent_variable(obs_eq.args[2], :($ivs_get_expr_sorted...)) - push!(observed_vars.args[1].args, obs_expr) + push!(obsexpr.args[1].args, obs_expr) end - # In case metadata was given, this must be cleared from `observed_eqs`. + # In case metadata was given, this must be cleared from `obs_eqs`. # For interpolated observables (I.e. $X ~ ...) this should and cannot be done. - is_escaped_expr(obs_eq.args[2]) || (observed_eqs.args[idx].args[2] = obs_name) + is_escaped_expr(obs_eq.args[2]) || (obs_eqs.args[idx].args[2] = obs_name) # Adds the observable to the list of observable names. # This is required for filtering away so these are not added to the ReactionSystem's species list. @@ -844,67 +821,107 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted; req is_escaped_expr(obs_eq.args[2]) || push!(obs_syms.args, obs_name) end - # If nothing was added to `observed_vars`, it has to be modified not to throw an error. - (striplines(observed_vars) == striplines(Expr(:block, :(@variables)))) && - (observed_vars = :()) + # If nothing was added to `obsexpr`, it has to be modified not to throw an error. + (striplines(obsexpr) == striplines(Expr(:block, :(@variables)))) && + (obsexpr = :()) else - # If option is not used, return empty expression and vector. - observed_vars = :() - observed_eqs = :([]) + # If observables option is not used, return empty expression and vector. + obsexpr = :() + obs_eqs = :([]) obs_syms = :([]) end - return observed_vars, observed_eqs, obs_syms + + return obsexpr, obs_eqs, obs_syms end -# From the input to the @observables options, creates a vector containing one equation for each observable. -# Checks separate cases for "@observables O ~ ..." and "@observables begin ... end". Other cases errors. -function make_observed_eqs(observables_expr) - if observables_expr.head == :call - return :([$(observables_expr)]) - elseif observables_expr.head == :block - observed_eqs = :([]) - for arg in observables_expr.args - push!(observed_eqs.args, arg) - end - return observed_eqs +# From the input to the @observables options, create a vector containing one equation for +# each observable. `option_block_form` handles if single line declaration of `@observables`, +# i.e. without a `begin ... end` block, was used. +function make_obs_eqs(observables_expr) + observables_expr = option_block_form(observables_expr) + obs_eqs = :([]) + foreach(arg -> push!(obs_eqs.args, arg), observables_expr.args) + return obs_eqs +end + +# Reads the combinatorial ratelaw options, which determines if a combinatorial rate law should +# be used or not. If not provided, use the default (true). +function read_combinatoric_ratelaws_option(options) + return haskey(options, :combinatoric_ratelaws) ? + get_block_option(options[:combinatoric_ratelaws]) : true +end + +# Finds the time independent variable, and any potential spatial independent variables. +# Returns these (individually and combined), as well as an expression for declaring them. +function read_ivs_option(options) + # Creates the independent variables expressions (depends on whether the `ivs` option was used). + if haskey(options, :ivs) + ivs = Tuple(extract_syms(options, :ivs)) + ivsexpr = copy(options[:ivs]) + ivsexpr.args[1] = Symbol("@", "parameters") else - error("Malformed observables option usage: $(observables_expr).") + ivs = (DEFAULT_IV_SYM,) + ivsexpr = :($(DEFAULT_IV_SYM) = default_t()) end + + # Extracts the independent variables symbols (time and spatial), and returns the output. + tiv = ivs[1] + sivs = (length(ivs) > 1) ? Expr(:vect, ivs[2:end]...) : nothing + return tiv, sivs, ivs, ivsexpr end ### `@reaction` Macro & its Internals ### -@doc raw""" -@reaction +""" + @reaction + +Macro for generating a single [`Reaction`](@ref) object using a similar syntax as the `@reaction_network` +macro (but permitting only a single reaction). A more detailed introduction to the syntax can +be found in the description of `@reaction_network`. -Generates a single [`Reaction`](@ref) object. +The `@reaction` macro is followed by a single line consisting of three parts: +- A rate (at which the reaction occurs). +- Any number of substrates (which are consumed by the reaction). +- Any number of products (which are produced by the reaction). + +The output is a reaction (just like created using the `Reaction` constructor). Examples: +Here we create a simple binding reaction and store it in the variable rx: +```julia +rx = @reaction k, X + Y --> XY +``` +The macro will automatically deduce `X`, `Y`, and `XY` to be species (as these occur as reactants) +and `k` as a parameter (as it does not occur as a reactant). + +The `@reaction` macro provides a more concise notation to the `Reaction` constructor. I.e. here +we create the same reaction using both approaches, and also confirm that they are identical. ```julia +# Creates a reaction using the `@reaction` macro. rx = @reaction k*v, A + B --> C + D -# is equivalent to +# Creates a reaction using the `Reaction` constructor. t = default_t() @parameters k v @species A(t) B(t) C(t) D(t) -rx == Reaction(k*v, [A,B], [C,D]) +rx2 = Reaction(k*v, [A, B], [C, D]) + +# Confirms that the two approaches yield identical results: +rx1 == rx2 ``` -Here `k` and `v` will be parameters and `A`, `B`, `C` and `D` will be variables. -Interpolation of existing parameters/variables also works +Interpolation of already declared symbolic variables into `@reaction` is possible: ```julia t = default_t() @parameters k b @species A(t) ex = k*A^2 + t -rx = @reaction b*$ex*$A, $A --> C +rx = @reaction b*\$ex*\$A, \$A --> C ``` Notes: -- Any symbols arising in the rate expression that aren't interpolated are treated as -parameters. In the reaction part (`α*A + B --> C + D`), coefficients are treated as -parameters, e.g. `α`, and rightmost symbols as species, e.g. `A,B,C,D`. -- Works with any *single* arrow types supported by [`@reaction_network`](@ref). -- Interpolation of Julia variables into the macro works similar to the `@reaction_network` +- `@reaction` does not support bi-directional type reactions (using `<-->`) or reaction bundling +(e.g. `d, (X,Y) --> 0`). +- Interpolation of Julia variables into the macro works similarly to the `@reaction_network` macro. See [The Reaction DSL](@ref dsl_description) tutorial for more details. """ macro reaction(ex) @@ -914,68 +931,61 @@ end # Function for creating a Reaction structure (used by the @reaction macro). function make_reaction(ex::Expr) - # Handle interpolation of variables + # Handle interpolation of variables in the input. ex = esc_dollars!(ex) - # Parses reactions, species, and parameters. + # Parses reactions. Extracts species and parameters within it. reaction = get_reaction(ex) - species, parameters = extract_species_and_parameters!([reaction], []) + species, parameters = extract_sps_and_ps([reaction], []) - # Checks for input errors. + # Checks for input errors. Needed here but not in `@reaction_network` as `ReactionSystem` performs this check but `Reaction` doesn't. forbidden_symbol_check(union(species, parameters)) - # 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}()) - rxexpr = get_rxexprs(reaction) + # Creates expressions corresponding to code for declaring the parameters, species, and reaction. + spexprs = get_usexpr(species, Dict{Symbol, Expr}()) + pexprs = get_psexpr(parameters, Dict{Symbol, Expr}()) + rxexpr = get_rxexpr(reaction) iv = :($(DEFAULT_IV_SYM) = default_t()) - # Returns the rephrased expression. + # Returns a rephrased expression which generates the `Reaction`. quote $pexprs $iv - $sexprs + $spexprs $rxexpr end end -# Reads a single line and creates the corresponding ReactionStruct. +# Reads a single line and creates the corresponding DSLReaction. function get_reaction(line) reaction = get_reactions([line]) - if (length(reaction) != 1) + (length(reaction) != 1) && error("Malformed reaction. @reaction macro only creates a single reaction. E.g. double arrows, such as `<-->` are not supported.") - end - return reaction[1] + return only(reaction) end ### Generic Expression Manipulation ### -# 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) +# 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, syms_skip = []) (typeof(expr) != Expr) && (return expr) - foreach(i -> expr.args[i] = recursive_escape_functions!(expr.args[i]), + foreach(i -> expr.args[i] = recursive_escape_functions!(expr.args[i], syms_skip), 1:length(expr.args)) - if expr.head == :call - !isdefined(Catalyst, expr.args[1]) && (expr.args[1] = esc(expr.args[1])) + if (expr.head == :call) && !isdefined(Catalyst, expr.args[1]) && expr.args[1] ∉ syms_skip + expr.args[1] = esc(expr.args[1]) end expr end -# Recursively escape functions in the right-hand-side of an equation written using user-defined functions. Special function calls like "hill(...)" are not expanded. -function escape_equation_RHS!(eqexpr::Expr) - rhs = recursive_escape_functions!(eqexpr.args[3]) - eqexpr.args[3] = rhs - eqexpr -end - -# Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical). +# Returns the length of an expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical). function tup_leng(ex::ExprValues) (typeof(ex) == Expr && ex.head == :tuple) && (return length(ex.args)) return 1 end -# Gets the ith element in a expression tuple, or returns the input itself if it is not an expression tuple -# (probably a Symbol/Numerical). +# Gets the ith element in an expression tuple, or returns the input itself if it is not an expression tuple +# (probably a Symbol/Numerical). This is used to handle bundled reactions (like `d, (X,Y) --> 0`). function get_tup_arg(ex::ExprValues, i::Int) (tup_leng(ex) == 1) && (return ex) return ex.args[i] diff --git a/src/expression_utils.jl b/src/expression_utils.jl index c1faa8ca8c..b839844fcf 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -2,13 +2,15 @@ # Function that handles variable interpolation. function esc_dollars!(ex) - if ex isa Expr - if ex.head == :$ - return esc(:($(ex.args[1]))) - else - for i in 1:length(ex.args) - ex.args[i] = esc_dollars!(ex.args[i]) - end + # If we do not have an expression: recursion has finished and we return the input. + (ex isa Expr) || (return ex) + + # If we have encountered an interpolation, perform the appropriate modification, else recur. + if ex.head == :$ + return esc(:($(ex.args[1]))) + else + for i in eachindex(ex.args) + ex.args[i] = esc_dollars!(ex.args[i]) end end ex @@ -22,33 +24,28 @@ end ### Parameters/Species/Variables Symbols Correctness Checking ### # Throws an error when a forbidden symbol is used. -function forbidden_symbol_check(v) - !isempty(intersect(forbidden_symbols_error, v)) && - error("The following symbol(s) are used as species or parameters: " * - ((map(s -> "'" * string(s) * "', ", - intersect(forbidden_symbols_error, v))...)) * - "this is not permitted.") - nothing -end - -# Throws an error when a forbidden variable is used (a forbidden symbol that is not `:t`). -function forbidden_variable_check(v) - !isempty(intersect(forbidden_variables_error, v)) && - error("The following symbol(s) are used as variables: " * - ((map(s -> "'" * string(s) * "', ", - intersect(forbidden_variables_error, v))...)) * - "this is not permitted.") +function forbidden_symbol_check(syms) + used_forbidden_syms = intersect(forbidden_symbols_error, syms) + isempty(used_forbidden_syms) || + error("The following symbol(s) are used as species or parameters: $used_forbidden_syms, this is not permitted.") end -function unique_symbol_check(syms) - allunique(syms) || - error("Reaction network independent variables, parameters, species, and variables must all have distinct names, but a duplicate has been detected. ") - nothing -end ### Catalyst-specific Expressions Manipulation ### -# Some options takes input on form that is either `@option ...` or `@option begin ... end`. +# Many option inputs can be on a form `@option input` or `@option begin ... end`. In both these +# cases we want to retrieve the third argument in the option expression. Further more, we wish +# to throw an error if there is more inputs (suggesting e.g. multiple inputs on a single line). +# Note that there are only some options for which we wish to make this check. +function get_block_option(expr) + (length(expr.args) < 3) && + error("The $(expr.args[1]) option's input was misformatted (full declaration: `$expr`). It seems that it has no inputs, whereas some input is expected.") + (length(expr.args) > 3) && + error("The $(expr.args[1]) option's input was misformatted (full declaration: `$expr`). Potentially, it has multiple inputs on a single line, in which case these should be split across multiple lines using a `begin ... end` block.") + return expr.args[3] +end + +# Some options takes input on form that is either `@option ...` or `@option begin ... end`. # This transforms input of the latter form to the former (with only one line in the `begin ... end` block) function option_block_form(expr) (expr.head == :block) && return expr @@ -74,12 +71,12 @@ function find_varinfo_in_declaration(expr) # Case: X (expr isa Symbol) && (return expr, [], nothing, nothing) - # Case: X(t) + # Case: X(t) (expr.head == :call) && (return expr.args[1], expr.args[2:end], nothing, nothing) if expr.head == :(=) # Case: X = 1.0 (expr.args[1] isa Symbol) && (return expr.args[1], [], expr.args[2], nothing) - # Case: X(t) = 1.0 + # Case: X(t) = 1.0 (expr.args[1].head == :call) && (return expr.args[1].args[1], expr.args[1].args[2:end], expr.args[2].args[1], nothing) @@ -116,7 +113,7 @@ end # (In this example the independent variable :t was inserted). # Here, the iv is a iv_expr, which can be anything, which is inserted function insert_independent_variable(expr_in, iv_expr) - # If expr is a symbol, just attach the iv. If not we have to create a new expr and mutate it. + # If expr is a symbol, just attach the iv. If not we have to create a new expr and mutate it. # Because Symbols (a possible input) cannot be mutated, this function cannot mutate the input # (would have been easier if Expr input was guaranteed). (expr_in isa Symbol) && (return Expr(:call, expr_in, iv_expr)) diff --git a/src/spatial_reaction_systems/spatial_reactions.jl b/src/spatial_reaction_systems/spatial_reactions.jl index 3b846d3049..421435cf98 100644 --- a/src/spatial_reaction_systems/spatial_reactions.jl +++ b/src/spatial_reaction_systems/spatial_reactions.jl @@ -57,9 +57,10 @@ function make_transport_reaction(rateex, species) # Checks for input errors. forbidden_symbol_check(union([species], parameters)) + # 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}()) + sexprs = get_usexpr([species], Dict{Symbol, Expr}()) + pexprs = get_psexpr(parameters, Dict{Symbol, Expr}()) iv = :($(DEFAULT_IV_SYM) = default_t()) trxexpr = :(TransportReaction($rateex, $species)) diff --git a/test/dsl/dsl_advanced_model_construction.jl b/test/dsl/dsl_advanced_model_construction.jl index 092f7bd0be..6dd389f637 100644 --- a/test/dsl/dsl_advanced_model_construction.jl +++ b/test/dsl/dsl_advanced_model_construction.jl @@ -75,7 +75,6 @@ end ### Test Interpolation Within the DSL ### - # Declares parameters and species used across the test. @parameters α k k1 k2 @species A(t) B(t) C(t) D(t) @@ -123,19 +122,6 @@ let @test rn == rn2 end -# Creates a reaction network using `eval` and internal function. -let - ex = quote - (Ka, Depot --> Central) - (CL / Vc, Central --> 0) - end - # Line number nodes aren't ignored so have to be manually removed - Base.remove_linenums!(ex) - exsys = Catalyst.make_reaction_system(ex) - sys = @eval Catalyst $exsys - @test sys isa ReactionSystem -end - # Miscellaneous interpolation tests. Unsure what they do here (not related to DSL). let rx = @reaction k*h, A + 2*B --> 3*C + D @@ -214,8 +200,8 @@ end # Checks that repeated metadata throws errors. let - @test_throws LoadError @eval @reaction k, 0 --> X, [md1=1.0, md1=2.0] - @test_throws LoadError @eval @reaction_network begin + @test_throws Exception @eval @reaction k, 0 --> X, [md1=1.0, md1=2.0] + @test_throws Exception @eval @reaction_network begin k, 0 --> X, [md1=1.0, md1=1.0] end end @@ -269,6 +255,24 @@ let @test isequal(rn1,rn2) end +# Tests that erroneous metadata declarations yields errors. +let + # Malformed metadata/value separator. + @test_throws Exception @eval @reaction_network begin + d, X --> 0, [misc=>"Metadata should use `=`, not `=>`."] + end + + # Malformed lhs + @test_throws Exception @eval @reaction_network begin + d, X --> 0, [misc,description=>"Metadata lhs should be a single symbol."] + end + + # Malformed metadata separator. + @test_throws Exception @eval @reaction_network begin + d, X --> 0, [misc=>:misc; description="description"] + end +end + ### Other Tests ### # Test floating point stoichiometry work. @@ -352,7 +356,7 @@ let @species (X(t))[1:2] (k[1],k[2]), X[1] <--> X[2] end - + @parameters k[1:2] @species (X(t))[1:2] rx1 = Reaction(k[1], [X[1]], [X[2]]) diff --git a/test/dsl/dsl_basic_model_construction.jl b/test/dsl/dsl_basic_model_construction.jl index b79c229c6e..ae2a638ae2 100644 --- a/test/dsl/dsl_basic_model_construction.jl +++ b/test/dsl/dsl_basic_model_construction.jl @@ -55,6 +55,65 @@ end ## Run Tests ### +# Tests the various network constructors. Test for `@network_component` and `@network_component`. +# Tests for combinations of reactions/no reactions, no name/name/interpolated name. +let + # Declare comparison networks programmatically. + @parameters d + @species X(t) + rx = Reaction(d, [X], []) + + rs_empty = ReactionSystem([], t; name = :name) + rs = ReactionSystem([rx], t; name = :name) + rs_empty_comp = complete(rs_empty) + rs_comp = complete(rs) + + # Declare empty networks. + name_sym = :name + rs_empty_1 = @network_component + rs_empty_2 = @network_component name + rs_empty_3 = @network_component $name_sym + rs_empty_comp_1 = @reaction_network + rs_empty_comp_2 = @reaction_network name + rs_empty_comp_3 = @reaction_network $name_sym + + # Check that empty networks are correct. + isequivalent(rs_empty_1, rs_empty) + rs_empty_2 == rs_empty + rs_empty_3 == rs_empty + isequivalent(rs_empty_comp_1, rs_empty_comp) + rs_empty_comp_2 == rs_empty_comp + rs_empty_comp_3 == rs_empty_comp + + # Declare non-empty networks. + rs_1 = @network_component begin + d, X --> 0 + end + rs_2 = @network_component name begin + d, X --> 0 + end + rs_3 = @network_component $name_sym begin + d, X --> 0 + end + rs_comp_1 = @reaction_network begin + d, X --> 0 + end + rs_comp_2 = @reaction_network name begin + d, X --> 0 + end + rs_comp_3 = @reaction_network $name_sym begin + d, X --> 0 + end + + # Check that non-empty networks are correct. + isequivalent(rs_1, rs) + rs_2 == rs + rs_3 == rs + isequivalent(rs_empty_1, rs_empty) + rs_empty_2 == rs_empty + rs_empty_3 == rs_empty +end + # Test basic properties of networks. let basic_test(reaction_networks_standard[1], 10, [:X1, :X2, :X3], @@ -130,7 +189,7 @@ let u0 = rnd_u0(networks[1], rng; factor) p = rnd_ps(networks[1], rng; factor) t = rand(rng) - + @test f_eval(networks[1], u0, p, t) ≈ f_eval(networks[2], u0, p, t) @test jac_eval(networks[1], u0, p, t) ≈ jac_eval(networks[2], u0, p, t) @test g_eval(networks[1], u0, p, t) ≈ g_eval(networks[2], u0, p, t) @@ -148,7 +207,7 @@ let (l3, l4), Y2 ⟷ Y3 (l5, l6), Y3 ⟷ Y4 c, Y4 → ∅ - end + end # Checks that the networks' functions evaluates equally for various randomised inputs. @unpack X1, X2, X3, X4, p, d, k1, k2, k3, k4, k5, k6 = network @@ -156,10 +215,10 @@ let u0_1 = Dict(rnd_u0(network, rng; factor)) p_1 = Dict(rnd_ps(network, rng; factor)) u0_2 = [:Y1 => u0_1[X1], :Y2 => u0_1[X2], :Y3 => u0_1[X3], :Y4 => u0_1[X4]] - p_2 = [:q => p_1[p], :c => p_1[d], :l1 => p_1[k1], :l2 => p_1[k2], :l3 => p_1[k3], + p_2 = [:q => p_1[p], :c => p_1[d], :l1 => p_1[k1], :l2 => p_1[k2], :l3 => p_1[k3], :l4 => p_1[k4], :l5 => p_1[k5], :l6 => p_1[k6]] t = rand(rng) - + @test f_eval(network, u0_1, p_1, t) ≈ f_eval(differently_written_5, u0_2, p_2, t) @test jac_eval(network, u0_1, p_1, t) ≈ jac_eval(differently_written_5, u0_2, p_2, t) @test g_eval(network, u0_1, p_1, t) ≈ g_eval(differently_written_5, u0_2, p_2, t) @@ -212,7 +271,7 @@ let u0 = rnd_u0(networks[1], rng; factor) p = rnd_ps(networks[1], rng; factor) t = rand(rng) - + @test f_eval(networks[1], u0, p, t) ≈ f_eval(networks[2], u0, p, t) @test jac_eval(networks[1], u0, p, t) ≈ jac_eval(networks[2], u0, p, t) @test g_eval(networks[1], u0, p, t) ≈ g_eval(networks[2], u0, p, t) @@ -234,7 +293,7 @@ let (sqrt(3.7), exp(1.9)), X4 ⟷ X1 + X2 end push!(identical_networks_3, reaction_networks_standard[9] => no_parameters_9) - push!(parameter_sets, [:p1 => 1.5, :p2 => 1, :p3 => 2, :d1 => 0.01, :d2 => 2.3, :d3 => 1001, + push!(parameter_sets, [:p1 => 1.5, :p2 => 1, :p3 => 2, :d1 => 0.01, :d2 => 2.3, :d3 => 1001, :k1 => π, :k2 => 42, :k3 => 19.9, :k4 => 999.99, :k5 => sqrt(3.7), :k6 => exp(1.9)]) no_parameters_10 = @reaction_network begin @@ -246,14 +305,14 @@ let 1.0, X5 ⟶ ∅ end push!(identical_networks_3, reaction_networks_standard[10] => no_parameters_10) - push!(parameter_sets, [:p => 0.01, :k1 => 3.1, :k2 => 3.2, :k3 => 0.0, :k4 => 2.1, :k5 => 901.0, + push!(parameter_sets, [:p => 0.01, :k1 => 3.1, :k2 => 3.2, :k3 => 0.0, :k4 => 2.1, :k5 => 901.0, :k6 => 63.5, :k7 => 7, :k8 => 8, :d => 1.0]) for (networks, p_1) in zip(identical_networks_3, parameter_sets) for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3] u0 = rnd_u0(networks[1], rng; factor) t = rand(rng) - + @test f_eval(networks[1], u0, p_1, t) ≈ f_eval(networks[2], u0, [], t) @test jac_eval(networks[1], u0, p_1, t) ≈ jac_eval(networks[2], u0, [], t) @test g_eval(networks[1], u0, p_1, t) ≈ g_eval(networks[2], u0, [], t) @@ -324,7 +383,7 @@ let τ = rand(rng) u = rnd_u0(reaction_networks_conserved[1], rng; factor) p_2 = rnd_ps(time_network, rng; factor) - p_1 = [p_2; reaction_networks_conserved[1].k1 => τ; + p_1 = [p_2; reaction_networks_conserved[1].k1 => τ; reaction_networks_conserved[1].k4 => τ; reaction_networks_conserved[1].k5 => τ] @test f_eval(reaction_networks_conserved[1], u, p_1, τ) ≈ f_eval(time_network, u, p_2, τ) @@ -404,7 +463,7 @@ let @test rn1 == rn2 end -# Tests arrow variants in `@reaction`` macro . +# Tests arrow variants in `@reaction` macro. let @test isequal((@reaction k, 0 --> X), (@reaction k, X <-- 0)) @test isequal((@reaction k, 0 --> X), (@reaction k, X ⟻ 0)) @@ -412,8 +471,8 @@ let @test isequal((@reaction k, 0 --> X), (@reaction k, 0 ⥟ X)) end -# Test that symbols with special meanings, or that are forbidden, are handled properly. -let +# Test that symbols with special meanings are handled properly. +let test_network = @reaction_network begin t * k, X --> ∅ end @test length(species(test_network)) == 1 @test length(parameters(test_network)) == 1 @@ -432,11 +491,74 @@ let @test length(species(test_network)) == 1 @test length(parameters(test_network)) == 0 @test reactions(test_network)[1].rate == ℯ +end + +### Error Test ### + +# Erroneous `@reaction` usage. +let + # Bi-directional reaction using the `@reaction` macro. + @test_throws Exception @eval @reaction (k1,k2), X1 <--> X2 - @test_throws LoadError @eval @reaction im, 0 --> B - @test_throws LoadError @eval @reaction nothing, 0 --> B - @test_throws LoadError @eval @reaction k, 0 --> im - @test_throws LoadError @eval @reaction k, 0 --> nothing + # Bundles reactions. + @test_throws Exception @eval @reaction k, (X1,X2) --> 0 end +# Tests that malformed reactions yields errors. +let + # Checks that malformed combinations of entries yields errors. + @test_throws Exception @eval @reaction_network begin + d, X --> 0, Y --> 0 + end + @test_throws Exception @eval @reaction_network begin + d, X --> 0, [misc="Ok metadata"], [description="Metadata in (erroneously) extra []."] + end + # Checks that incorrect bundling yields error. + @test_throws Exception @eval @reaction_network begin + (k1,k2,k3), (X1,X2) --> 0 + end + + # Checks that incorrect stoichiometric expression yields an error. + @test_throws Exception @eval @reaction_network begin + k, X^Y --> XY + end +end + +# Check that forbidden symbols correctly generates errors. +let + # @reaction macro, symbols that cannot be in the rate. + @test_throws Exception @eval @reaction im, 0 --> X + @test_throws Exception @eval @reaction nothing, 0 --> X + @test_throws Exception @eval @reaction Γ, 0 --> X + @test_throws Exception @eval @reaction ∅, 0 --> X + + # @reaction macro, symbols that cannot be a reactant. + @test_throws Exception @eval @reaction 1, 0 --> im + @test_throws Exception @eval @reaction 1, 0 --> nothing + @test_throws Exception @eval @reaction 1, 0 --> Γ + @test_throws Exception @eval @reaction 1, 0 --> ℯ + @test_throws Exception @eval @reaction 1, 0 --> pi + @test_throws Exception @eval @reaction 1, 0 --> π + @test_throws Exception @eval @reaction 1, 0 --> t + + # @reaction_network macro, symbols that cannot be in the rate. + @test_throws Exception @eval @reaction_network begin im, 0 --> X end + @test_throws Exception @eval @reaction_network begin nothing, 0 --> X end + @test_throws Exception @eval @reaction_network begin Γ, 0 --> X end + @test_throws Exception @eval @reaction_network begin ∅, 0 --> X end + + # @reaction_network macro, symbols that cannot be a reactant. + @test_throws Exception @eval @reaction_network begin 1, 0 --> im end + @test_throws Exception @eval @reaction_network begin 1, 0 --> nothing end + @test_throws Exception @eval @reaction_network begin 1, 0 --> Γ end + @test_throws Exception @eval @reaction_network begin 1, 0 --> ℯ end + @test_throws Exception @eval @reaction_network begin 1, 0 --> pi end + @test_throws Exception @eval @reaction_network begin 1, 0 --> π end + @test_throws Exception @eval @reaction_network begin 1, 0 --> t end + + # Checks that non-supported arrow type usage yields error. + @test_throws Exception @eval @reaction_network begin + d, X ⇻ 0 + end +end \ No newline at end of file diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl index 16ef93d1f7..71c8c4417a 100644 --- a/test/dsl/dsl_options.jl +++ b/test/dsl/dsl_options.jl @@ -415,36 +415,45 @@ let spcs = (A, B1, B2, C) @test issetequal(unknowns(rn), sts) @test issetequal(species(rn), spcs) +end - @test_throws ArgumentError begin - rn = @reaction_network begin - @variables K - k, K*A --> B - end +# Tests error when disallowed name is used for variable. +let + @test_throws Exception @eval @reaction_network begin + @variables π(t) + end + @test_throws Exception @eval @reaction_network begin + @variables Γ(t) + end + @test_throws Exception @eval @reaction_network begin + @variables ∅(t) + end + @test_throws Exception @eval @reaction_network begin + @ivs s + @variables t(s) end end # Tests that explicitly declaring a single symbol as several things does not work. # Several of these are broken, but note sure how to test broken-ness on `@test_throws false Exception @eval`. -# Relevant issue: https://github.com/SciML/Catalyst.jl/issues/1173 let # Species + parameter. - @test_broken false #@test_throws Exception @eval @reaction_network begin - #@species X(t) - #@parameters X - #end + @test_throws Exception @eval @reaction_network begin + @species X(t) + @parameters X + end # Species + variable. - @test_broken false #@test_throws Exception @eval @reaction_network begin - #@species X(t) - #@variables X(t) - #end + @test_throws Exception @eval @reaction_network begin + @species X(t) + @variables X(t) + end # Variable + parameter. - @test_broken false #@test_throws Exception @eval @reaction_network begin - #@variables X(t) - #@parameters X - #end + @test_throws Exception @eval @reaction_network begin + @variables X(t) + @parameters X + end # Species + differential. @test_throws Exception @eval @reaction_network begin @@ -465,31 +474,31 @@ let end # Parameter + observable (species/variable + observable is OK, as this e.g. provide additional observables information). - @test_broken false #@test_throws Exception @eval @reaction_network begin - #@species Y(t) - #@parameters X - #@observables X ~ Y - #end + @test_throws Exception @eval @reaction_network begin + @species Y(t) + @parameters X + @observables X ~ Y + end # Species + compound. - @test_broken false #@test_throws Exception @eval @reaction_network begin - #@species X(t) O(t) - #@compounds begin X(t) ~ 2O end - #end + @test_throws Exception @eval @reaction_network begin + @species X(t) O(t) + @compounds begin X(t) ~ 2O end + end # Parameter + compound. - @test_broken false #@test_throws Exception @eval @reaction_network begin - #@species O(t) - #@parameters X - #@compounds begin X(t) ~ 2O end - #end + @test_throws Exception @eval @reaction_network begin + @species O(t) + @parameters X + @compounds begin X(t) ~ 2O end + end # Variable + compound. - @test_broken false #@test_throws Exception @eval @reaction_network begin - #@species O(t) - #@variables X(t) - #@compounds begin X(t) ~ 2O end - #end + @test_throws Exception @eval @reaction_network begin + @species O(t) + @variables X(t) + @compounds begin X(t) ~ 2O end + end end ### Test Independent Variable Designations ### @@ -609,11 +618,11 @@ let @equations D(V) ~ 1 - V d, D --> 0 end - @test_broken false # @test_throws Exception @eval @reaction_network begin - #@variables D(t) - #@equations D(V) ~ 1 - V - #d, X --> 0 - #end + @test_throws Exception @eval @reaction_network begin + @variables D(t) + @equations D(V) ~ 1 - V + d, X --> 0 + end @test_throws Exception @eval @reaction_network begin @parameters D @equations D(V) ~ 1 - V @@ -675,7 +684,7 @@ let @test plot(sol; idxs=:X).series_list[1].plotattributes[:y][end] ≈ 10.0 @test plot(sol; idxs=[X, Y]).series_list[2].plotattributes[:y][end] ≈ 3.0 @test plot(sol; idxs=[rn.X, rn.Y]).series_list[2].plotattributes[:y][end] ≈ 3.0 - @test plot(sol; idxs=[:X, :Y]).series_list[2].plotattributes[:y][end] ≈ 3.0 # (https://github.com/SciML/ModelingToolkit.jl/issues/2778) + @test plot(sol; idxs=[:X, :Y]).series_list[2].plotattributes[:y][end] ≈ 3.0 end # Compares programmatic and DSL system with observables. @@ -960,6 +969,14 @@ let @observables $X ~ X1 + X2 (k1,k2), X1 <--> X2 end + + # Observable metadata provided twice. + @test_throws Exception @eval @reaction_network begin + @species X2 [description="Twice the amount of X"] + @observables (X2, [description="X times two."]) ~ 2X + d, X --> 0 + end + end @@ -1128,55 +1145,184 @@ let @equations X = 1 - S (p,d), 0 <--> S end + + # Differential equation using a forbidden variable (in the DSL). + @test_throws Exception @eval @reaction_network begin + @equations D(π) ~ -1 + end + + # Algebraic equation using a forbidden variable (in the DSL). + @test_throws Exception @eval @reaction_network begin + @equations Γ ~ 1 + 3(Γ^2 + Γ) + end +end + +### Other DSL Option Tests ### + +# Test that various options can be provided in block and single line form. +# Also checks that the single line form takes maximally one argument. +let + # The `@equations` option. + rn11 = @reaction_network rn1 begin + @equations D(V) ~ 1 - V + end + rn12 = @reaction_network rn1 begin + @equations begin + D(V) ~ 1 - V + end + end + @test isequal(rn11, rn12) + @test_throws Exception @eval @reaction_network begin + @equations D(V) ~ 1 - V D(W) ~ 1 - W + end + + # The `@observables` option. + rn21 = @reaction_network rn1 begin + @species X(t) + @observables X2 ~ 2X + end + rn22 = @reaction_network rn1 begin + @species X(t) + @observables begin + X2 ~ 2X + end + end + @test isequal(rn21, rn22) + @test_throws Exception @eval @reaction_network begin + @species X(t) + @observables X2 ~ 2X X3 ~ 3X + end + + # The `@compounds` option. + rn31 = @reaction_network rn1 begin + @species X(t) + @compounds X2 ~ 2X + end + rn32 = @reaction_network rn1 begin + @species X(t) + @compounds begin + X2 ~ 2X + end + end + @test isequal(rn31, rn32) + @test_throws Exception @eval @reaction_network begin + @species X(t) + @compounds X2 ~ 2X X3 ~ 3X + end + + # The `@differentials` option. + rn41 = @reaction_network rn1 begin + @differentials D = Differential(t) + end + rn42 = @reaction_network rn1 begin + @differentials begin + D = Differential(t) + end + end + @test isequal(rn41, rn42) + @test_throws Exception @eval @reaction_network begin + @differentials D = Differential(t) Δ = Differential(t) + end + + # The `@continuous_events` option. + rn51 = @reaction_network rn1 begin + @species X(t) + @continuous_events [X ~ 3.0] => [X ~ X - 1] + end + rn52 = @reaction_network rn1 begin + @species X(t) + @continuous_events begin + [X ~ 3.0] => [X ~ X - 1] + end + end + @test isequal(rn51, rn52) + @test_throws Exception @eval @reaction_network begin + @species X(t) + @continuous_events [X ~ 3.0] => [X ~ X - 1] [X ~ 1.0] => [X ~ X + 1] + end + + # The `@discrete_events` option. + rn61 = @reaction_network rn1 begin + @species X(t) + @discrete_events [X > 3.0] => [X ~ X - 1] + end + rn62 = @reaction_network rn1 begin + @species X(t) + @discrete_events begin + [X > 3.0] => [X ~ X - 1] + end + end + @test isequal(rn61, rn62) + @test_throws Exception @eval @reaction_network begin + @species X(t) + @discrete_events [X > 3.0] => [X ~ X - 1] [X < 1.0] => [X ~ X + 1] + end end # test combinatoric_ratelaws DSL option let + # Test for `@combinatoric_ratelaws false`. rn = @reaction_network begin @combinatoric_ratelaws false (k1,k2), 2A <--> B - end + end combinatoric_ratelaw = Catalyst.get_combinatoric_ratelaws(rn) @test combinatoric_ratelaw == false rl = oderatelaw(reactions(rn)[1]; combinatoric_ratelaw) @unpack k1, A = rn @test isequal(rl, k1*A^2) + # Test for `@combinatoric_ratelaws true`. rn2 = @reaction_network begin @combinatoric_ratelaws true (k1,k2), 2A <--> B - end + end combinatoric_ratelaw = Catalyst.get_combinatoric_ratelaws(rn2) @test combinatoric_ratelaw == true rl = oderatelaw(reactions(rn2)[1]; combinatoric_ratelaw) @unpack k1, A = rn2 @test isequal(rl, k1*A^2/2) + # Test for interpolation into `@combinatoric_ratelaws`. crl = false rn3 = @reaction_network begin @combinatoric_ratelaws $crl (k1,k2), 2A <--> B - end + end combinatoric_ratelaw = Catalyst.get_combinatoric_ratelaws(rn3) @test combinatoric_ratelaw == crl rl = oderatelaw(reactions(rn3)[1]; combinatoric_ratelaw) @unpack k1, A = rn3 @test isequal(rl, k1*A^2) + + # Test for erroneous inputs (to few, to many, wrong type). + @test_throws Exception @eval @reaction_network begin + @combinatoric_ratelaws + d, 3X --> 0 + end + @test_throws Exception @eval @reaction_network begin + @combinatoric_ratelaws false false + d, 3X --> 0 + end + @test_throws Exception @eval @reaction_network begin + @combinatoric_ratelaws "false" + d, 3X --> 0 + end end # Test whether user-defined functions are properly expanded in equations. let f(A, t) = 2*A*t + g(A) = 2*A + 2 - # Test user-defined function + # Test user-defined function (on both lhs and rhs). rn = @reaction_network begin - @equations D(A) ~ f(A, t) + @equations D(A) + g(A) ~ f(A, t) end @test length(equations(rn)) == 1 @test equations(rn)[1] isa Equation @species A(t) - @test isequal(equations(rn)[1], D(A) ~ 2*A*t) - + @test isequal(equations(rn)[1], D(A) + 2*A + 2 ~ 2*A*t) # Test whether expansion happens properly for unregistered/registered functions. hill_unregistered(A, v, K, n) = v*(A^n) / (A^n + K^n) @@ -1201,7 +1347,6 @@ let @parameters v K n @test isequal(equations(rn2r)[1], D(A) ~ hill2(A, v, K, n)) - rn3 = @reaction_network begin @species Iapp(t) @equations begin @@ -1223,7 +1368,6 @@ let rn3_sym = complete(rn3_sym) @test isequivalent(rn3, rn3_sym) - # 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 @@ -1236,10 +1380,18 @@ let @test isequal(Catalyst.expand_registered_functions(equations(rn4)[1]), D(A) ~ v*(A^n)) end -### test that @no_infer properly throws errors when undeclared variables are written ### +# Erroneous `@default_noise_scaling` declaration (other noise scaling tests are mostly in the SDE file). +let + # Default noise scaling with multiple entries. + @test_throws Exception @eval @reaction_network begin + @default_noise_scaling η1 η2 + end +end -import Catalyst: UndeclaredSymbolicError +# test that @require_declaration properly throws errors when undeclared variables are written. let + import Catalyst: UndeclaredSymbolicError + # Test error when species are inferred @test_throws UndeclaredSymbolicError @macroexpand @reaction_network begin @require_declaration