diff --git a/HISTORY.md b/HISTORY.md index f3f72c5d10..902c6a597b 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -81,6 +81,7 @@ end plot_network(brusselator) ``` +- The letter Ø (used in Danish/Norwegian alphabet) is now conisdred the same as ∅ (empty set). It can no longer be used as a species/parameter. ## Catalyst 14.4.1 - Support for user-defined functions on the RHS when providing coupled equations diff --git a/Project.toml b/Project.toml index a3846c6ee9..a2bb6789da 100644 --- a/Project.toml +++ b/Project.toml @@ -32,6 +32,7 @@ BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a" StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" @@ -39,7 +40,7 @@ StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" [extensions] CatalystBifurcationKitExtension = "BifurcationKit" CatalystCairoMakieExtension = "CairoMakie" -CatalystGraphMakieExtension = ["GraphMakie", "NetworkLayout"] +CatalystGraphMakieExtension = ["GraphMakie", "NetworkLayout", "Makie"] CatalystHomotopyContinuationExtension = "HomotopyContinuation" CatalystStructuralIdentifiabilityExtension = "StructuralIdentifiability" @@ -57,8 +58,9 @@ Graphs = "1.4" HomotopyContinuation = "2.9" JumpProcesses = "9.13.2" LaTeXStrings = "1.3.0" -Latexify = "0.16.5" +Latexify = "0.16.6" MacroTools = "0.5.5" +Makie = "0.22.1" ModelingToolkit = "< 9.60" NetworkLayout = "0.4.7" Parameters = "0.12" diff --git a/docs/Project.toml b/docs/Project.toml index 9affe0ea77..c28df9ce52 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -62,7 +62,7 @@ HomotopyContinuation = "2.9" IncompleteLU = "0.2" JumpProcesses = "9.13.2" Latexify = "0.16.5" -LinearSolve = "2.30" +LinearSolve = "2.30, 3" ModelingToolkit = "< 9.60" NetworkLayout = "0.4" NonlinearSolve = "3.12, 4" diff --git a/docs/src/api.md b/docs/src/api.md index f78a5f70e4..f099341206 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 allow one to supply additional information to a [`ReactionSystem`](@ref) created with 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..b03b22d4c2 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). Then, we specify the same model by instead creating 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/docs/src/model_creation/model_visualisation.md b/docs/src/model_creation/model_visualisation.md index d378deb3a6..7555929944 100644 --- a/docs/src/model_creation/model_visualisation.md +++ b/docs/src/model_creation/model_visualisation.md @@ -88,20 +88,18 @@ In this section we demonstrate some of the ways that plot objects can be manipul f, ax, p = plot_complexes(brusselator, show_rate_labels = true) ``` -It seems like a bit of the top node is cut off. Let's hide the tick marks and grid and increase the top and bottom margins by increasing `yautolimitmargin`. +It seems like a bit of the top node is cut off. Let's increase the top and bottom margins by increasing `yautolimitmargin`. ```@example visualisation_graphs -hidedecorations!(ax) -ax.yautolimitmargin = (0.1, 0.1) # defaults to (0.05, 0.05) +ax.yautolimitmargin = (0.3, 0.3) # defaults to (0.15, 0.15) ax.aspect = DataAspect() +f ``` -There are many keyword arguments that can be passed to `plot_network` or `plot_complexes` to change the look of the graph (which get passed to the `graphplot` Makie recipe). Let's change the color of the nodes and make the inner labels a bit smaller. As before, we hide the tick marks and grid. Let's also give the plot a title. +There are many keyword arguments that can be passed to `plot_network` or `plot_complexes` to change the look of the graph (which get passed to the `graphplot` Makie recipe). Let's change the color of the nodes and make the inner labels a bit smaller. Let's also give the plot a title. ```@example visualisation_graphs f, ax, p = plot_complexes(brusselator, show_rate_labels = true, node_color = :yellow, ilabels_fontsize = 10) -hidedecorations!(ax) -ax.yautolimitmargin = (0.1, 0.1) # defaults to (0.05, 0.05) -ax.aspect = DataAspect() ax.title = "Brusselator" +f ``` Most of the kwargs that modify the nodes or edges will also accept a vector with the same length as the number of nodes or edges, respectively. See [here](https://graph.makie.org/stable/#The-graphplot-Recipe) for a full list of keyword arguments to `graph_plot`. Note that `plot_complexes` and `plot_network` default to `layout = Stress()` rather than `layout = Spring()`, since `Stress()` is better at generating plots with fewer edge crossings. More layout options and customizations (such as pinning nodes to certain positions) can be found in the [`NetworkLayout` documentation](https://juliagraphs.org/NetworkLayout.jl/stable/). @@ -109,6 +107,7 @@ Most of the kwargs that modify the nodes or edges will also accept a vector with Once a graph is already created we can also change the keyword arguments by modifying the fields of the `Plot` object `p`. ```@example visualisation_graphs p.node_color = :orange +f ``` Custom node positions can also be given, if the automatic layout is unsatisfactory. @@ -116,6 +115,7 @@ Custom node positions can also be given, if the automatic layout is unsatisfacto fixedlayout = [(0,0), (1,0), (0,1), (1,1), (2,0)] p.layout = fixedlayout autolimits!(ax) +f ``` Makie graph plots can be made to be interactive, allowing one to drag nodes and edges. To do this, we retrieve the axis from the GraphMakie plot, and then register the interactions. **Note that this can only be done if `GLMakie` is the installed Makie backend. See the [GraphMakie docs](https://graph.makie.org/stable/#Predefined-Interactions) for more information about the types of interactions one can register.** Below is a non-interactive code example that shows how to do this: diff --git a/docs/src/spatial_modelling/lattice_reaction_systems.md b/docs/src/spatial_modelling/lattice_reaction_systems.md index 6182c2b59c..ea107b8034 100644 --- a/docs/src/spatial_modelling/lattice_reaction_systems.md +++ b/docs/src/spatial_modelling/lattice_reaction_systems.md @@ -60,9 +60,9 @@ lat_getu(sol, :X1, lrs) and plot the simulation using ```@example spatial_intro_basics import CairoMakie -lattice_animation(sol, :X1, lrs, "lattice_simulation_2d.mp4") +lattice_animation(sol, :X1, lrs, "lattice_simulation.mp4") ``` -![](./lattice_simulation_2d.mp4) +![](./lattice_simulation.mp4) More information on how to retrieve values from spatial simulations can be found [here](@ref lattice_simulation_structure_interaction_simulation_species), and for plotting them, [here](@ref lattice_simulation_plotting). Finally, a list of functions for querying `LatticeReactionSystems` for various properties can be found [here](@ref api_lattice_simulations). ## [Spatial reactions](@id spatial_lattice_modelling_intro_spatial_reactions) diff --git a/ext/CatalystCairoMakieExtension/cairo_makie_extension_spatial_modelling.jl b/ext/CatalystCairoMakieExtension/cairo_makie_extension_spatial_modelling.jl index 7dc1dd3cad..5b719558f8 100644 --- a/ext/CatalystCairoMakieExtension/cairo_makie_extension_spatial_modelling.jl +++ b/ext/CatalystCairoMakieExtension/cairo_makie_extension_spatial_modelling.jl @@ -1,6 +1,6 @@ ### 1d Lattice Simulation Plots/Animations ### -# Internal dispatch for the plotting of a lattice simulation on a 1d lattice (Cartesian or masked). +# Internal dispatch for the plotting of a lattice simulation on a 1d lattice (Cartesian or masked). function lattice_plot( sol, sp, lrs::LatticeReactionSystem{Q, R, <:Catalyst.GridLattice{1, S}, T}; t = sol.t[end], markersize = 20, kwargs...) where {Q, R, S, T} @@ -12,7 +12,7 @@ function lattice_plot( markersize = markersize, kwargs...) end -# Internal dispatch for the animation of a lattice simulation on a 1d lattice (Cartesian or masked). +# Internal dispatch for the animation of a lattice simulation on a 1d lattice (Cartesian or masked). function lattice_animation( sol, sp, lrs::LatticeReactionSystem{Q, R, <:Catalyst.GridLattice{1, S}, T}, filename::String; @@ -40,7 +40,7 @@ function lattice_animation( return nothing end -# Internal dispatch for the kymographs of a lattice simulation on a 1d lattice (Cartesian or masked). +# Internal dispatch for the kymographs of a lattice simulation on a 1d lattice (Cartesian or masked). function lattice_kymograph( sol, sp, lrs::LatticeReactionSystem{Q, R, <:Catalyst.GridLattice{1, S}, T}; colormap = :BuGn_7, @@ -64,7 +64,7 @@ end ### 2d Lattice Simulation Plots/Animations ### -# Internal dispatch for the plotting of a lattice simulation on a 2d lattice (Cartesian or masked). +# Internal dispatch for the plotting of a lattice simulation on a 2d lattice (Cartesian or masked). function lattice_plot(sol, sp, lrs::LatticeReactionSystem{Q, R, <:Catalyst.GridLattice{2, S}, T}; t = sol.t[end], colormap = :BuGn_7, plot_min = nothing, plot_max = nothing, @@ -80,14 +80,14 @@ function lattice_plot(sol, sp, return heatmap(x_vals, y_vals, vals; - axis = (xlabel = "Time", ylabel = "Compartment", + axis = (xlabel = "Compartment", ylabel = "Compartment", xgridvisible = false, ygridvisible = false), colormap, colorrange = (plot_min, plot_max), kwargs...) end -# Internal dispatch for the animation of a lattice simulation on a 2d lattice (Cartesian or masked). +# Internal dispatch for the animation of a lattice simulation on a 2d lattice (Cartesian or masked). function lattice_animation( sol, sp, lrs::LatticeReactionSystem{Q, R, <:Catalyst.GridLattice{2, S}, T}, filename::String; @@ -101,7 +101,7 @@ function lattice_animation( # Creates the base figure (which is modified in the animation). fig, ax, hm = heatmap(x_vals, y_vals, vals[1]; - axis = (xgridvisible = false, ygridvisible = false), + axis = (xgridvisible = false, ygridvisible = false, xlabel = "Compartment", ylabel = "Compartment"), colormap, colorrange = (plot_min, plot_max), kwargs...) ttitle && (ax.title = "Time: $(round(t[1]; sigdigits = 3))") @@ -116,14 +116,14 @@ end ### 3d Lattice Simulation Plots/Animations (Errors Only) ### -# Internal dispatch for the plotting of a lattice simulation on a 3d lattice (Cartesian or masked). +# Internal dispatch for the plotting of a lattice simulation on a 3d lattice (Cartesian or masked). function lattice_plot( sol, sp, lrs::LatticeReactionSystem{Q, R, <:Catalyst.GridLattice{3, S}, T}; kwargs...) where {Q, R, S, T} throw(ArgumentError("The `lattice_plot` function does not support 3d Cartesian/masked lattices.")) end -# Internal dispatch for the animation of a lattice simulation on a 3d lattice (Cartesian or masked). +# Internal dispatch for the animation of a lattice simulation on a 3d lattice (Cartesian or masked). function lattice_animation( sol, sp, lrs::LatticeReactionSystem{Q, R, <:Catalyst.GridLattice{3, S}, T}, filename::String; kwargs...) where {Q, R, S, T} diff --git a/ext/CatalystGraphMakieExtension.jl b/ext/CatalystGraphMakieExtension.jl index a87dcf9af3..0a63fa1415 100644 --- a/ext/CatalystGraphMakieExtension.jl +++ b/ext/CatalystGraphMakieExtension.jl @@ -1,7 +1,7 @@ module CatalystGraphMakieExtension # Fetch packages. -using Catalyst, GraphMakie, Graphs, Symbolics, SparseArrays, NetworkLayout +using Catalyst, GraphMakie, Graphs, Symbolics, SparseArrays, NetworkLayout, Makie using Symbolics: get_variables! import Catalyst: species_reaction_graph, incidencematgraph, lattice_plot, lattice_animation diff --git a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl index 26271e2b38..be83f6f8e6 100644 --- a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl +++ b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl @@ -209,6 +209,9 @@ function Catalyst.plot_network(rn::ReactionSystem; kwargs...) f.axis.xautolimitmargin = (0.15, 0.15) f.axis.yautolimitmargin = (0.15, 0.15) + hidedecorations!(f.axis) + hidespines!(f.axis) + f.axis.aspect = DataAspect() f end @@ -260,9 +263,13 @@ function Catalyst.plot_complexes(rn::ReactionSystem; show_rate_labels = false, k curve_distance = gen_distances(cg), kwargs... ) + f.axis.xautolimitmargin = (0.15, 0.15) f.axis.yautolimitmargin = (0.15, 0.15) - + hidedecorations!(f.axis) + hidespines!(f.axis) + f.axis.aspect = DataAspect() + f end diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 52daa0cec5..c8c2af25c6 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -69,14 +69,9 @@ const ExprValues = Union{Expr, Symbol, Float64, Int, Bool} const CONSERVED_CONSTANT_SYMBOL = :Γ # Declares symbols which may neither be used as parameters nor unknowns. -const forbidden_symbols_skip = Set([:ℯ, :pi, :π, :t, :∅]) +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..7207edcf7b 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -1,67 +1,7 @@ -""" -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}([:∅]) +const empty_set = Set{Symbol}([:∅, :Ø]) const fwd_arrows = Set{Symbol}([:>, :(=>), :→, :↣, :↦, :⇾, :⟶, :⟼, :⥟, :⥟, :⇀, :⇁, :⇒, :⟾]) const bwd_arrows = Set{Symbol}([:<, :(<=), :←, :↢, :↤, :⇽, :⟵, :⟻, :⥚, :⥞, :↼, :↽, :⇐, :⟽, Symbol("<--")]) @@ -110,159 +50,183 @@ 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 +([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/)) sections 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. """ @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, :(t = default_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 +235,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 +249,7 @@ function extract_metadata(metadata_line::Expr) return metadata end - - +### Specialised @require_declaration Option Error ### struct UndeclaredSymbolicError <: Exception msg::String end @@ -298,149 +262,116 @@ 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] + allunique(arg.args[1] for arg in option_lines) || + error("Some options where given multiple times.") + numlines = length(reaction_lines) + length(option_lines) + (numlines != length(ex.args)) && + error("@reaction_network input contain $(length(ex.args) - numlines) malformed lines.") + options = Dict(Symbol(String(arg.args[1])[2:end]) => arg for arg in option_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", scalarize = false) + 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 + $vsexpr + $spsexpr + $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 = _observed, continuous_events = _continuous_events, + discrete_events = _discrete_events, combinatoric_ratelaws = _combinatoric_ratelaws); + default_reaction_metadata = _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 +389,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 +412,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)) + maxlen = maximum(lengs) + any(!(leng == 1 || leng == maxlen) 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:maxlen + # 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 +437,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,55 +448,59 @@ 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 + excluded_syms = 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 -# 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 called by `extract_sps_and_ps`, recursively loops through an expression and find +# symbols (adding them to the push_symbols vector). Returns `nothing` to ensure type stability. +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 @@ -572,99 +508,144 @@ 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 + for u in us_extracted + u isa Symbol && push!(usexpr.args, Expr(:call, u, ivs...)) end - foreach(s -> (s isa Symbol) && push!(sexprs.args, Expr(:call, s, iv_symbols...)), - species_extracted) - sexprs + 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 : :([]) - subs_stoich_init = deepcopy(subs_init) - prod_init = isempty(rxstruct.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) - end - for prod in rxstruct.products - push!(reaction_func.args[4].args, prod.reactant) - push!(reaction_func.args[6].args, prod.stoichiometry) - end - reaction_func +# 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 -# 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) +# 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(rx.products) ? nothing : :([]) + prod_stoich_init = deepcopy(prod_init) + 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 rx.products + push!(rx_constructor.args[5].args, prod.reactant) + push!(rx_constructor.args[7].args, prod.stoichiometry) + end + return rx_constructor +end + +# 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. +# If requested, performs scalarization. +function assign_var_to_symvar_declaration(expr_init, name; scalarize = true) + # 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 the input expression is non-empty, wrap it with additional information. + if expr_init != :(()) if scalarize symvec = gensym() - ex = quote - $symvec = $ex + expr = quote + $symvec = $expr_init $namesym = reduce(vcat, Symbolics.scalarize($symvec)) end else - ex = quote - $namesym = $ex + expr = quote + $namesym = $expr_init end 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 +654,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 +672,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 +682,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 +697,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 +742,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 +810,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 +826,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("@", "independent_variables") 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`. + +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). -Generates a single [`Reaction`](@ref) object. +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 +936,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/latexify_recipes.jl b/src/latexify_recipes.jl index ecbd42f6d1..4759f76937 100644 --- a/src/latexify_recipes.jl +++ b/src/latexify_recipes.jl @@ -14,14 +14,14 @@ const LATEX_DEFS = CatalystLatexParams() @latexrecipe function f(rs::ReactionSystem; form = :reactions, expand_functions = true) expand_functions && (rs = expand_registered_functions(rs)) if form == :reactions # Returns chemical reaction network code. - cdot --> false + mult_symbol --> "" env --> :chem return rs elseif form == :ode # Returns ODE system code. - cdot --> false + mult_symbol --> "" return convert(ODESystem, rs) elseif form == :sde # Returns SDE system code. - cdot --> false + mult_symbol --> "" return convert(SDESystem, rs) end error("Unrecognised form argument given: $form. This should be either reactions (default), :ode, or :sde.") diff --git a/src/spatial_reaction_systems/spatial_reactions.jl b/src/spatial_reaction_systems/spatial_reactions.jl index 3b846d3049..fb861efcc6 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)) @@ -154,7 +155,7 @@ end # Loops through a rate and extracts all parameters. function find_parameters_in_rate!(parameters, rateex::ExprValues) if rateex isa Symbol - if rateex in [:t, :∅, :im, :nothing, CONSERVED_CONSTANT_SYMBOL] + if rateex in [:t, :∅, :Ø, :im, :nothing, CONSERVED_CONSTANT_SYMBOL] error("Forbidden term $(rateex) used in transport reaction rate.") elseif !(rateex in [:ℯ, :pi, :π]) push!(parameters, rateex) diff --git a/test/dsl/dsl_advanced_model_construction.jl b/test/dsl/dsl_advanced_model_construction.jl index 092f7bd0be..29a205a88c 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 @@ -150,6 +136,20 @@ let @test rx == Reaction(b+ex, [A,C], nothing, [2,1], nothing) 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) + name = QuoteNode(:rs) + exsys = Catalyst.make_reaction_system(ex, name) + sys = @eval Catalyst $exsys + @test sys isa ReactionSystem +end + ### Tests Reaction Metadata ### # Tests construction for various types of reaction metadata. @@ -214,8 +214,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 +269,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 +370,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..ccff3cd81f 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, τ) @@ -385,6 +444,35 @@ let @test any(isequal(I), unknowns(rn)) end +# Test that Ø (Danish/Norwegian letter), ∅ (empty set), and 0 (zero) are equivalent. +let + rn1 = @reaction_network rn begin + p, Ø --> X + d, X --> Ø + end + rn2 = @reaction_network rn begin + p, ∅ --> X + d, X --> ∅ + end + rn3 = @reaction_network rn begin + p, 0 --> X + d, X --> 0 + end + rn4 = @reaction_network rn begin + p, Ø --> X + d, X --> ∅ + end + rn5 = @reaction_network rn begin + p, Ø --> X + d, X --> 0 + end + rn6 = @reaction_network rn begin + p, ∅ --> X + d, X --> 0 + end + @test rn1 == rn2 == rn3 == rn4 == rn5 == rn6 +end + # Tests backwards and bi-directional arrows. let rn1 = @reaction_network arrowtest begin @@ -404,7 +492,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,7 +500,7 @@ 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. +# 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 @@ -432,11 +520,74 @@ let @test length(species(test_network)) == 1 @test length(parameters(test_network)) == 0 @test reactions(test_network)[1].rate == ℯ +end - @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 +### Error Test ### + +# Erroneous `@reaction` usage. +let + # Bi-directional reaction using the `@reaction` macro. + @test_throws Exception @eval @reaction (k1,k2), X1 <--> X2 + + # 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 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 diff --git a/test/reactionsystem_core/custom_crn_functions.jl b/test/reactionsystem_core/custom_crn_functions.jl index 2cfe60ae65..e5f2387761 100644 --- a/test/reactionsystem_core/custom_crn_functions.jl +++ b/test/reactionsystem_core/custom_crn_functions.jl @@ -49,7 +49,7 @@ let ps = rnd_ps(custom_function_network_1, rng; factor) t = rand(rng) @test f_eval(custom_function_network_1, u0, ps, t) ≈ f_eval(custom_function_network_2, u0, ps, t) - @test jac_eval(custom_function_network_1, u0, ps, t) ≈ jac_eval(custom_function_network_2, u0, ps, t) + @test_broken jac_eval(custom_function_network_1, u0, ps, t) ≈ jac_eval(custom_function_network_2, u0, ps, t) # Weird error dur to some SciML update. Reported in https://github.com/SciML/ModelingToolkit.jl/issues/3447. @test g_eval(custom_function_network_1, u0, ps, t) ≈ g_eval(custom_function_network_2, u0, ps, t) end end @@ -101,20 +101,20 @@ end # Tests `ReactionSystem`s. let @species x(t) y(t) - @parameters k v n + @parameters k v n rs1 = @reaction_network rs begin mm(x, v, k), 0 --> x mmr(x, v, k), 0 --> x hill(x, v, k, n), 0 --> x hillr(x, v, k, n), 0 --> x - hillar(x, y, v, k, n), 0 --> x + hillar(x, y, v, k, n), 0 --> x end rs2 = @reaction_network rs begin v * x / (x + k), 0 --> x v * k / (x + k), 0 --> x v * (x^n) / (x^n + k^n), 0 --> x v * (k^n) / (x^n + k^n), 0 --> x - v * (x^n) / (x^n + y^n + k^n), 0 --> x + v * (x^n) / (x^n + y^n + k^n), 0 --> x end @test Catalyst.expand_registered_functions(rs1) == rs2 @@ -123,14 +123,14 @@ end # Tests `Reaction`s. let @species x(t) y(t) - @parameters k v n - + @parameters k v n + r1 = @reaction mm(x, v, k), 0 --> x r2 = @reaction mmr(x, v, k), 0 --> x r3 = @reaction hill(x, v, k, n), 0 --> x r4 = @reaction hillr(x, v, k, n), 0 --> x r5 = @reaction hillar(x, y, v, k, n), 0 --> x + y - + @test isequal(Catalyst.expand_registered_functions(r1).rate, v * x / (x + k)) @test isequal(Catalyst.expand_registered_functions(r2).rate, v * k / (x + k)) @test isequal(Catalyst.expand_registered_functions(r3).rate, v * (x^n) / (x^n + k^n)) @@ -143,14 +143,14 @@ end let @parameters T @variables X(T) Y(T) - @parameters K V N - + @parameters K V N + eq1 = 0 ~ mm(X, V, K) eq2 = 0 ~ mmr(X, V, K) eq3 = 0 ~ hill(X, V, K, N) eq4 = 0 ~ hillr(X, V, K, N) eq5 = 0 ~ hillar(X, Y, V, K, N) - + @test isequal(Catalyst.expand_registered_functions(eq1), 0 ~ V * X / (X + K)) @test isequal(Catalyst.expand_registered_functions(eq2), 0 ~ V * K / (X + K)) @test isequal(Catalyst.expand_registered_functions(eq3), 0 ~ V * (X^N) / (X^N + K^N)) @@ -166,7 +166,7 @@ let @parameters v K eqs = [ Reaction(mm(X,v,K), [], [X]), - mm(V,v,K) ~ V + 1 + mm(V,v,K) ~ V + 1 ] @named rs = ReactionSystem(eqs, t) @@ -211,7 +211,7 @@ let (v * (X^n) / (X^n + K^n) > 1000.0) => [X ~ v * (K^n) / (X^n + K^n) + 2] ] continuous_events = ModelingToolkit.SymbolicContinuousCallback.(continuous_events) - discrete_events = ModelingToolkit.SymbolicDiscreteCallback.(discrete_events) + discrete_events = ModelingToolkit.SymbolicDiscreteCallback.(discrete_events) @test isequal(only(Catalyst.get_rxs(rs_expanded)).rate, v0 + v * (X^n) / (X^n + Y^n + K^n)) @test isequal(get_continuous_events(rs_expanded), continuous_events) @test isequal(get_discrete_events(rs_expanded), discrete_events)