new inference of D and tests and History file

Author: TorkelE

HISTORY.md

@@ -7,6 +7,29 @@
   (at the time the release is made). If you need a dependency version increased,
   please open an issue and we can update it and make a new Catalyst release once
   testing against the newer dependency version is complete. 
+- New formula for inferring variables from equations (declared using the `@equations` options) in the DSL. The order of inference of species/variables/parameters is now:
+    (1) Every symbol explicitly declared using `@species`, `@variables`, and `@parameters` are assigned to the correct category.
+    (2) Every symbol used as a reaction reactant is inferred as a species.
+    (3) Every symbol not declared in (1) or (2) that occurs in an expression provided after `@equations` is inferred as a variable.
+    (4) Every symbol not declared in (1), (2), or (3) that occurs either as a reaction rate or stoichiometric coefficient is inferred to be a parameter.
+E.g. in
+```julia
+@reaction_network begin
+    @equations V1 + S ~ V2^2
+    (p + S + V1), S --> 0
+end
+```
+`S` is inferred as a species, `V1` and `V2` as variables, and `p` as a parameter. The previous special cases for the `@observables`, `@compounds`, and `@differentials` options still hold. Finally, the `@require_declaration` options (described in more detail below) can now be used to require everything to be explicitly declared.
+- New formula for determining whether the default differentials have been used within an `@equations` option. Right now, if any expression `D(...)` is encountered (where `...`) can be anything, this is inferred as usage of the default differential D. E.g. in the following equations `D` is inferred as a differential with respect to the default independent variable:
+```julia
+@reaction_network begin
+  @equations D(V) + V ~ 1
+end
+@reaction_network begin
+  @equations D(D(V)) ~ 1
+end
+```
+Please note that this cannot be used at the same time as `D` is used to represent a species, variable, or parameter. 
 - Array symbolics support is more consistent with ModelingToolkit v9. Parameter
   arrays are no longer scalarized by Catalyst, while species and variables
   arrays still are (as in ModelingToolkit). As such, parameter arrays should now

src/reactionsystem.jl

@@ -524,8 +524,6 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
     # Extracts any species, variables, and parameters that occur in (non-reaction) equations.
     # Creates the new reactions + equations vector, `fulleqs` (sorted reactions first, equations next).
     if !isempty(eqs)
-        println(eqs)
-        println(iv)
         osys = ODESystem(eqs, iv; name = gensym())
         fulleqs = CatalystEqType[rxs; equations(osys)]
         union!(us, unknowns(osys))

test/dsl/dsl_options.jl

@@ -424,6 +424,74 @@ let
     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
+
+    # Species + variable.
+    @test_broken false #@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
+
+    # Species + differential.
+    @test_throws false Exception @eval @reaction_network begin
+        @species X(t)
+        @differentials X = Differential(t)
+    end
+
+    # Parameter + differential.
+    @test_throws false Exception @eval @reaction_network begin
+        @parameters X
+        @differentials X = Differential(t)
+    end
+
+    # Variable + differential.
+    @test_throws false Exception @eval @reaction_network begin
+        @variables X(t)
+        @differentials X = Differential(t)
+    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
+
+    # Species + compound.
+    @test_broken false #@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
+
+    # 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
+end
+
 ### Test Independent Variable Designations ###
 
 # Test ivs in DSL.
@@ -450,6 +518,118 @@ let
     @test issetequal(Catalyst.get_sivs(rn), [x])
 end
 
+### Test Symbolic Variable Inference ###
+
+# Basic checks that that symbolic variables not explicitly declared are correctly inferred.
+let
+    # Case 1 (a reaction only).
+    rn1 = @reaction_network begin
+        (p1/(S1+p2) + S2), S1 --> S2
+    end
+    @test issetequal(species(rn1), [rn1.S1, rn1.S2])
+    @test issetequal(parameters(rn1), [rn1.p1, rn1.p2])
+
+    # Case 2 (reactions and equations).
+    rn2 = @reaction_network begin
+        @equations V1 + log(V2 + S1) ~ V2^2
+        (p1/V1 + S1 + log(S2 + V2 + p2)), S1 --> S2
+    end
+    @test issetequal(species(rn2), [rn2.S1, rn2.S2])
+    @test issetequal(nonspecies(rn2), [rn2.V1, rn2.V2])
+    @test issetequal(parameters(rn2), [rn2.p1, rn2.p2])
+
+    # Case 3 (reaction and equations with a differential).
+    rn3 = @reaction_network begin
+        @equations begin
+            D(V1) ~ S1 + V1
+            V2 + S2 ~ V1^2 + V2^2
+        end
+        (p1/V1 + S1 + log(S2 + V2 + p2)), S1 --> S2
+    end
+    @test issetequal(species(rn3), [rn3.S1, rn3.S2])
+    @test issetequal(nonspecies(rn3), [rn3.V1, rn3.V2])
+    @test issetequal(parameters(rn3), [rn3.p1, rn3.p2])
+
+    # Case 4 (reactions and equations with a pre-declared parameter).
+    rn4 = @reaction_network begin
+        @parameters p1
+        @equations V1 + sin(p1 + S1) ~ S2*V2
+        (p1+p2+V1+V2+S1+S2), S1 --> S2
+    end
+    @test issetequal(species(rn4), [rn2.S1, rn2.S2])
+    @test issetequal(nonspecies(rn4), [rn2.V1, rn2.V2])
+    @test issetequal(parameters(rn4), [rn2.p1, rn2.p2])
+
+    # Case 5 (algebraic equation containing D, which is pre-declared as a species).
+    rn5 = @reaction_network begin
+        @species D(t)
+        @equations D * (S1 + V1 + V2) ~ S2
+        (p1 + p2*(D + V1 + V2 + S2 + S2)), S1 --> S2 + D
+    end
+    @test issetequal(species(rn5), [rn5.S1, rn5.S2, rn5.D])
+    @test issetequal(nonspecies(rn5), [rn5.V1, rn5.V2])
+    @test issetequal(parameters(rn5), [rn5.p1, rn5.p2])
+
+    # Case 6 (algebraic equation containing D, which is pre-declared as a parameter).
+    rn6 = @reaction_network begin
+        @parameters D
+        @equations D * (S1 + V1 + V2) ~ S2
+        (p1 + p2*(D + V1 + V2 + S2 + S2)), S1 --> S2
+    end
+    @test issetequal(species(rn6), [rn6.S1, rn6.S2])
+    @test issetequal(nonspecies(rn6), [rn6.V1, rn6.V2])
+    @test issetequal(parameters(rn6), [rn6.p1, rn6.p2, rn6.D])
+
+    # Case 7 (algebraic equation containing D, which is pre-declared as a variable).
+    rn7 = @reaction_network begin
+        @variables D(t)
+        @equations D * (S1 + V1 + V2) ~ S2
+        (p1 + p2*(D + V1 + V2 + S2 + S2)), S1 --> S2
+    end
+    @test issetequal(species(rn7), [rn7.S1, rn7.S2])
+    @test issetequal(nonspecies(rn7), [rn7.V1, rn7.V2, rn7.D])
+    @test issetequal(parameters(rn7), [rn7.p1, rn7.p2])
+
+    # Case 8 (reactions, equations, and a custom differential).
+    rn8 = @reaction_network begin
+        @differentials Δ = Differential(t)
+        @equations Δ(V1) + Δ(V2) + log(V2 + S1) ~  S2
+        (p1/V1 + S1 + log(S2 + V2 + p2)), S1 --> S2
+    end
+    @test issetequal(species(rn8), [rn8.S1, rn8.S2])
+    @test issetequal(nonspecies(rn8), [rn8.V1, rn8.V2])
+    @test issetequal(parameters(rn8), [rn8.p1, rn8.p2])
+end
+
+# Checks that various cases where symbolic variables cannot (or shouldn't) be inferred generate errors.
+let
+    # Species/variables/parameter named after default differential used as function call.
+    # In the future, species/variables should be usable this way (designating a time delay).
+    @test_throws Exception @eval @reaction_network begin
+        @equations D(V) ~ 1 - V
+        d, D --> 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
+        d, X --> 0
+    end
+
+    # Symbol only occurring in events.
+    @test_throws Exception @eval @reaction_network begin
+        @discrete_event (X > 1.0) => [V => V/2]
+        d, X --> 0
+    end
+    @test_throws Exception @eval @reaction_network begin
+        @continuous_event [X > 1.0] => [V => V/2]
+        d, X --> 0
+    end
+end
 
 ### Observables ###
 
@@ -790,7 +970,7 @@ end
 # Check that DAE is solved correctly.
 let
     rn = @reaction_network rn begin
-        @parameters k
+        @parameters k d
         @variables X(t) Y(t)
         @equations begin
             X + 5 ~ k*S
@@ -940,10 +1120,10 @@ let
     @test 5*sol[:Y][end] ≈ sol[:S][end] + sol[:X][end]
 end
 
-# Tests that the correct symbolic variables are infered as species, variables, and paraemters.
+# Tests that the correct symbolic variables are inferred as species, variables, and parameters.
 let
     rn = @reaction_network begin
-        @parameters p1 p2
+        @parameters p1 p2 k1 k2 k3 k4
         @species X1(t) X2(t)
         @variables W(t)
         @equations begin
@@ -968,12 +1148,6 @@ let
         @equations X = 1 - S
         (p,d), 0 <--> S
     end
-
-    # Equation with component undeclared elsewhere.
-    @test_throws Exception @eval @reaction_network begin
-        @equations X ~ p - S
-        (P,D), 0 <--> S
-    end
 end
 
 # test combinatoric_ratelaws DSL option
@@ -1082,7 +1256,7 @@ 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
+### test that @no_infer properly throws errors when undeclared variables are written ###
 
 import Catalyst: UndeclaredSymbolicError
 let

test/reactionsystem_core/coupled_equation_crn_systems.jl

@@ -899,38 +899,6 @@ end
 
 # Checks that various misformatted declarations yield errors.
 let
-    # Symbol in equation not appearing elsewhere (1).
-    @test_throws Exception @eval @reaction_network begin
-        @equations D(V) ~ -X
-    end
-
-    # Symbol in equation not appearing elsewhere (2).
-    @test_throws Exception @eval @reaction_network begin
-        @equations 1 + log(x) ~ 2X
-    end
-
-    # Attempting to infer differential variable not isolated on lhs (1).
-    @test_throws Exception @eval @reaction_network begin
-        @equations D(V) + 1 ~ 0
-    end
-
-    # Attempting to infer differential variable not isolated on lhs (2).
-    @test_throws Exception @eval @reaction_network begin
-        @equations -1.0 ~ D(V)
-    end
-
-    # Attempting to infer differential operator not isolated on lhs (1).
-    @test_throws Exception @eval @reaction_network begin
-        @variables V(t)
-        @equations D(V) + 1 ~ 0
-    end
-
-    # Attempting to infer a variable when using a non-default differential.
-    @test_throws Exception @eval @reaction_network begin
-        @differentials Δ = Differential(t)
-        @equations Δ(V) ~ -1,0
-    end
-
     # Attempting to create a new differential from an unknown iv.
     @test_throws Exception @eval @reaction_network begin
         @differentials D = Differential(τ)
@@ -944,14 +912,8 @@ let
 
     # Several equations without `begin ... end` block.
     @test_throws Exception @eval @reaction_network begin
-        @variables V(t)
         @equations D(V) + 1 ~ - 1.0
-    end
-
-    # Undeclared differential.
-    @test_throws Exception @eval @reaction_network begin
-        @species V
-        @equations Δ(V) ~ -1.0
+        @equations D(W) + 1 ~ - 1.0
     end
 
     # System using multiple ivs.