Merge remote-tracking branch 'origin/master' into graph-docs

Author: vyudu

Project.toml

@@ -46,7 +46,7 @@ BifurcationKit = "0.4.4"
 CairoMakie = "0.12"
 Combinatorics = "1.0.2"
 DataStructures = "0.18"
-DiffEqBase = "< 6.159.0"
+DiffEqBase = "6.159.0"
 DocStringExtensions = "0.8, 0.9"
 DynamicPolynomials = "0.5, 0.6"
 DynamicQuantities = "0.13.2, 1"
@@ -62,11 +62,11 @@ Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
-SciMLBase = "< 2.57.2"
+SciMLBase = "2.57.2"
 Setfield = "1"
 # StructuralIdentifiability = "0.5.8"
 SymbolicUtils = "2.1.2, 3.3.0"
-Symbolics = "5.30.1, 6"
+Symbolics = "6.22"
 Unitful = "1.12.4"
 julia = "1.10"
 

docs/Project.toml

@@ -50,7 +50,7 @@ BifurcationKit = "0.4.4"
 CairoMakie = "0.12"
 Catalyst = "14.4"
 DataFrames = "1.6"
-DiffEqBase = "< 6.159.0"
+DiffEqBase = "6.159.0"
 DiffEqParamEstim = "2.2"
 Distributions = "0.25"
 Documenter = "1.4.1"
@@ -87,4 +87,4 @@ SpecialFunctions = "2.4"
 StaticArrays = "1.9"
 SteadyStateDiffEq = "2.2"
 StochasticDiffEq = "6.65"
-Symbolics = "5.30.1, 6"
+Symbolics = "6.22"

docs/src/model_creation/dsl_advanced.md

@@ -272,23 +272,24 @@ sol = solve(oprob)
 plot(sol)
 ```
 
-### [Turning off species, parameter, and variable inferring](@id dsl_advanced_options_require_dec)
-In some cases it may be desirable for Catalyst to not infer species and parameters from the DSL, as in the case of reaction networks with very many variables, or as a sanity check that variable names are written correctly. To turn off inferring, simply add the `@require_declaration` macro to one of the lines of the `@reaction_network` declaration. Having this macro means that every single variable, species, or parameter will have to be explicitly declared using the `@variable`, `@species`, or `@parameter` macro. In the case that the DSL parser encounters an undeclared symbolic, it will error with an `UndeclaredSymbolicError` and print the reaction or equation that the undeclared symbolic was found in. 
+### [Turning off species, parameter, and variable inference](@id dsl_advanced_options_require_dec)
+In some cases it may be desirable for Catalyst to not infer species and parameters from the DSL, as in the case of reaction networks with very many variables, or as a sanity check that variable names are written correctly. To turn off inference, simply use the `@require_declaration` option when using the `@reaction_network` DSL. This will require every single variable, species, or parameter used within the DSL to be explicitly declared using the `@variable`, `@species`, or `@parameter` options. In the case that the DSL parser encounters an undeclared symbolic, it will error with an `UndeclaredSymbolicError` and print the reaction or equation that the undeclared symbolic was found in. 
 
-```@example dsl_advanced_no_infer
+```julia
 using Catalyst
-# The following case will throw an UndeclaredSymbolicError.
-try @macroexpand @reaction_network begin
+rn = @reaction_network begin
     @require_declaration
     (k1, k2), A <--> B
 end
-catch e
-    println(e.msg)
-end
 ```
-In order to avoid an error in this case all the relevant species and parameters will have to be declared.
+Running the code above will yield the following error: 
 ```
+LoadError: UndeclaredSymbolicError: Unrecognized variables A detected in reaction expression: "((k1, k2), A <--> B)". Since the flag @require_declaration is declared, all species must be explicitly declared with the @species macro.
+```
+In order to avoid the error in this case all the relevant species and parameters will have to be declared.
+```@example dsl_advanced_require_dec
 # The following case will not error. 
+using Catalyst
 t = default_t()
 rn = @reaction_network begin
     @require_declaration
@@ -299,11 +300,11 @@ end
 ```
 
 The following cases in which the DSL would normally infer variables will all throw errors if `@require_declaration` is set and the variables are not explicitly declared.
-- Inferring a species in a reaction, as in the example above
-- Inferring a parameter in a reaction rate expression, as in the reaction line `k*n, A --> B`
-- Inferring a parameter in the stoichiometry of a species, as in the reaction line `k, n*A --> B`
-- Inferring a differential variable on the LHS of a coupled differential equation, as in `A` in `@equations D(A) ~ A^2`
-- Inferring an [observable](@ref dsl_advanced_options_observables) that is declared using `@observables`
+- Occurrence of an undeclared species in a reaction, as in the example above.
+- Occurrence of an undeclared parameter in a reaction rate expression, as in the reaction line `k*n, A --> B`.
+- Occurrence of an undeclared parameter in the stoichiometry of a species, as in the reaction line `k, n*A --> B`.
+- Occurrence of an undeclared differential variable on the LHS of a coupled differential equation, as in `A` in `@equations D(A) ~ A^2`.
+- Occurrence of an undeclared [observable](@ref dsl_advanced_options_observables) in an `@observables` expression, such as `@observables X1 ~ A + B`.
 
 ## [Naming reaction networks](@id dsl_advanced_options_naming)
 Each reaction network model has a name. It can be accessed using the `nameof` function. By default, some generic name is used:

src/Catalyst.jl

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

src/chemistry_functionality.jl

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

src/dsl.jl

@@ -141,17 +141,17 @@ ReactionSystems generated through `@reaction_network` are complete.
 """
 macro reaction_network(name::Symbol, ex::Expr)
     :(complete($(make_reaction_system(
-        MacroTools.striplines(ex); name = :($(QuoteNode(name)))))))
+        striplines(ex); name = :($(QuoteNode(name)))))))
 end
 
 # Allows @reaction_network $name begin ... to interpolate variables storing a name.
 macro reaction_network(name::Expr, ex::Expr)
     :(complete($(make_reaction_system(
-        MacroTools.striplines(ex); name = :($(esc(name.args[1])))))))
+        striplines(ex); name = :($(esc(name.args[1])))))))
 end
 
 macro reaction_network(ex::Expr)
-    ex = MacroTools.striplines(ex)
+    ex = striplines(ex)
 
     # no name but equations: @reaction_network begin ... end ...
     if ex.head == :block
@@ -179,16 +179,16 @@ Equivalent to `@reaction_network` except the generated `ReactionSystem` is not m
 complete.
 """
 macro network_component(name::Symbol, ex::Expr)
-    make_reaction_system(MacroTools.striplines(ex); name = :($(QuoteNode(name))))
+    make_reaction_system(striplines(ex); name = :($(QuoteNode(name))))
 end
 
 # Allows @network_component $name begin ... to interpolate variables storing a name.
 macro network_component(name::Expr, ex::Expr)
-    make_reaction_system(MacroTools.striplines(ex); name = :($(esc(name.args[1]))))
+    make_reaction_system(striplines(ex); name = :($(esc(name.args[1]))))
 end
 
 macro network_component(ex::Expr)
-    ex = MacroTools.striplines(ex)
+    ex = striplines(ex)
 
     # no name but equations: @network_component begin ... end ...
     if ex.head == :block
@@ -328,7 +328,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters_declared = extract_syms(options, :parameters)
     variables_declared = extract_syms(options, :variables)
 
-    # Reads equations. 
+    # Reads equations.
     vars_extracted, add_default_diff, equations = read_equations_options(
         options, variables_declared; requiredec)
     variables = vcat(variables_declared, vars_extracted)
@@ -356,7 +356,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     observed_vars, observed_eqs, obs_syms = read_observed_options(
         options, [species_declared; variables], all_ivs; requiredec)
 
-    # Collect species and parameters, including ones inferred from the reactions. 
+    # Collect species and parameters, including ones inferred from the reactions.
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
         variables)))
     species_extracted, parameters_extracted = extract_species_and_parameters!(
@@ -397,7 +397,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         push!(rxexprs.args, equation)
     end
 
-    # Output code corresponding to the reaction system. 
+    # Output code corresponding to the reaction system.
     quote
         $ivexpr
         $ps
@@ -448,7 +448,7 @@ function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(u
             push_reactions!(reactions, reaction.args[3], reaction.args[2],
                 rate, metadata, arrow, line)
         else
-            throw("Malformed reaction, invalid arrow type used in: $(MacroTools.striplines(line))")
+            throw("Malformed reaction, invalid arrow type used in: $(striplines(line))")
         end
     end
     return reactions
@@ -670,7 +670,7 @@ function read_events_option(options, event_type::Symbol)
         error("Trying to read an unsupported event type.")
     end
     events_input = haskey(options, event_type) ? options[event_type].args[3] :
-                   MacroTools.striplines(:(begin end))
+                   striplines(:(begin end))
     events_input = option_block_form(events_input)
 
     # Goes through the events, checks for errors, and adds them to the output vector.
@@ -750,7 +750,7 @@ function create_differential_expr(options, add_default_diff, used_syms, tiv)
     # If differentials was provided as options, this is used as the initial expression.
     # If the default differential (D(...)) was used in equations, this is added to the expression.
     diffexpr = (haskey(options, :differentials) ? options[:differentials].args[3] :
-                MacroTools.striplines(:(begin end)))
+                striplines(:(begin end)))
     diffexpr = option_block_form(diffexpr)
 
     # Goes through all differentials, checking that they are correctly formatted and their symbol is not used elsewhere.
@@ -810,23 +810,18 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted; req
             # For Observables that have already been declared using @species/@variables,
             # or are interpolated, this parts should not be carried out.
             if !((obs_name in species_n_vars_declared) || is_escaped_expr(obs_eq.args[2]))
-                # Appends (..) to the observable (which is later replaced with the extracted ivs).
-                # Adds the observable to the first line of the output expression (starting with `@variables`).
-                obs_expr = insert_independent_variable(obs_eq.args[2], :(..))
-                push!(observed_vars.args[1].args, obs_expr)
-
-                # Adds a line to the `observed_vars` expression, setting the ivs for this observable.
-                # Cannot extract directly using e.g. "getfield.(dependants_structs, :reactant)" because
-                # then we get something like :([:X1, :X2]), rather than :([X1, X2]).
+                # 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,
                     Symbolics.get_variables($(obs_eq.args[3]))))
                 ivs_get_expr = :(unique(reduce(
                     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)))
-                push!(observed_vars.args,
-                    :($obs_name = $(obs_name)($(ivs_get_expr_sorted)...)))
+
+                obs_expr = insert_independent_variable(obs_eq.args[2], :($ivs_get_expr_sorted...))
+                push!(observed_vars.args[1].args, obs_expr)
             end
 
             # In case metadata was given, this must be cleared from `observed_eqs`.
@@ -840,7 +835,8 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted; req
         end
 
         # If nothing was added to `observed_vars`, it has to be modified not to throw an error.
-        (length(observed_vars.args) == 1) && (observed_vars = :())
+        (striplines(observed_vars) == striplines(Expr(:block, :(@variables)))) &&
+            (observed_vars = :())
     else
         # If option is not used, return empty expression and vector.
         observed_vars = :()
@@ -944,7 +940,7 @@ end
 
 ### Generic Expression Manipulation ###
 
-# Recursively traverses an expression and escapes all the user-defined functions. Special function calls like "hill(...)" are not expanded. 
+# Recursively traverses an expression and escapes all the user-defined functions. Special function calls like "hill(...)" are not expanded.
 function recursive_escape_functions!(expr::ExprValues)
     (typeof(expr) != Expr) && (return expr)
     foreach(i -> expr.args[i] = recursive_escape_functions!(expr.args[i]),

src/spatial_reaction_systems/spatial_reactions.jl

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