Merge remote-tracking branch 'origin/master' into isaacsas-uncap-bases

Author: isaacsas

HISTORY.md

@@ -25,6 +25,22 @@
   StructuralIdentifiability has with Julia 1.10.5 and 1.11.
 - A tutorial on making interactive plot displays using Makie has been added.
 - The BifurcationKit extension has been updated to v.4.
+- There is a new DSL option `@require_declaration` that will turn off automatic inferring for species, parameters, and variables in the DSL. For example, the following will now error:
+  ```julia
+  rn = @reaction_network begin
+    @require_declaration
+    (k1, k2), A <--> B
+  end
+  ```
+  When this flag is set, all symbolics must be explicitly declared. 
+  ```julia
+  rn = @reaction_network begin
+    @species A(t) B(t)
+    @parameters k1 k2
+    @require_declaration
+    (k1, k2), A <--> B
+  end
+  ```
 
 ## Catalyst 14.4.1
 - Support for user-defined functions on the RHS when providing coupled equations 

Project.toml

@@ -66,7 +66,7 @@ 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"
 
@@ -77,7 +77,11 @@ Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
+OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -93,7 +97,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["DiffEqCallbacks", "DomainSets", "Graphviz_jll", "Logging", "NonlinearSolve", 
-        "OrdinaryDiffEq", "Pkg", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", 
-        "StableRNGs", "StaticArrays", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", 
-        "Test", "Unitful"]
+test = ["DiffEqCallbacks", "DomainSets", "Graphviz_jll", "Logging", "NonlinearSolve", "OrdinaryDiffEqBDF", "OrdinaryDiffEqDefault", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner", "Pkg", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "StaticArrays", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]

docs/Project.toml

@@ -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/faqs.md

@@ -309,3 +309,16 @@ end
 In some cases, it may be necessary or desirable to register functions with
 Symbolics.jl before their use in Catalyst, see the discussion
 [here](https://symbolics.juliasymbolics.org/stable/manual/functions/).
+
+## How can I turn off automatic inferring of species and parameters when using the DSL?
+This option can be set using the `@require_declaration` option inside `@reaction_network`. In this case all the species, parameters, and variables in the system must be pre-declared using one of the `@species`, `@parameters`, or `@variables` macros. For more information about what is inferred automatically and not, please see the section on [`@require_declaration`](@ref dsl_advanced_options_require_dec).
+
+```@example faq9
+using Catalyst
+rn = @reaction_network begin
+    @require_declaration
+    @species A(t) B(t)
+    @parameters k1 k2
+    (k1, k2), A <--> B
+end
+```

docs/src/model_creation/dsl_advanced.md

@@ -205,7 +205,7 @@ 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.
+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
 rn = @reaction_network begin
@@ -272,6 +272,40 @@ sol = solve(oprob)
 plot(sol)
 ```
 
+### [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. 
+
+```julia
+using Catalyst
+rn = @reaction_network begin
+    @require_declaration
+    (k1, k2), A <--> B
+end
+```
+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
+    @species A(t) B(t)
+    @parameters k1 k2
+    (k1, k2), A <--> B
+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.
+- 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:
 ```@example dsl_advanced_names

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