handle weird reactions

Author: TorkelE

src/Catalyst.jl

@@ -147,7 +147,7 @@ export Graph, savegraph, complexgraph
 include("chemistry_functionality.jl")
 export @compound, @compounds
 export iscompound, components, coefficients, component_coefficients
-export balance_reaction
+export balance_reaction, balance_system
 
 ### Extensions ###
 

src/chemistry_functionality.jl

@@ -353,9 +353,38 @@ end
 From a system, creates a new system where each reaction is a balanced version of the corresponding
 reaction of the original system. For more information, consider the `balance_reaction` function
 (which is internally applied to each system reaction).
+
+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. 
+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. 
 """
 function balance_system(rs::ReactionSystem)
-    @set! rs.eqs = [(eq isa Reaction) ? balance_reaction(eq) : eq for eq in get_eqs(rs)]
-    @set! rs.rxs = [balance_reaction(rx) for rx in get_rxs(rs)]
+    @set! rs.eqs = CatalystEqType[get_balanced_reaction(eq) for eq in get_eqs(rs)]
+    @set! rs.rxs = [get_balanced_reaction(rx) for rx in get_rxs(rs)]
     return rs
-end
+end
+
+# Selects a balanced version of an input reaction. Handles potential problems when there are no,
+# or several, balanced alternatives.
+function get_balanced_reaction(rx::Reaction)
+    brxs = balance_reaction(rx)
+
+    # 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
+    if length(brxs) > 1
+        error("Infinite number of balanced reactions possible for reaction ($rx) are possible. No method for automatically generating a valid reaction is currently implemented in `balance_system`.")
+    end
+
+    return only(brxs)
+end
+# For non-`Reaction` equations, returns the original equation.
+get_balanced_reaction(eq::Equation) = eq

test/miscellaneous_tests/reaction_balancing.jl

@@ -433,15 +433,15 @@ end
 
 # Simple example with multiple reactions.
 let
-    # Creates a reaction system and tis balanced version.
+    # Creates an unbalanced `ReactionSystem`.
     rs = @reaction_network begin
         @species C(t) O(t) H(t)
         @compounds begin
-            H2(t) = 2H
-            CH4(t) = C + 4H
-            O2(t) = 2O
-            CO2(t) = C + 2O
-            H2O(t) = 2H + O
+            H2(t) ~ 2H
+            CH4(t) ~ C + 4H
+            O2(t) ~ 2O
+            CO2(t) ~ C + 2O
+            H2O(t) ~ 2H + O
         end
         1.0, C + H2 --> CH4
         2.0, CH4 + O2 --> CO2 + H2O
@@ -452,23 +452,23 @@ let
     # substrates/products for granted).
     rx1 = filter(rx -> isequal(rx.rate, 1.0), reactions(rs_balanced))[1]
     rx2 = filter(rx -> isequal(rx.rate, 2.0), reactions(rs_balanced))[1]
-    @test issetequal(rx1.substoich, [1,2])
+    @test issetequal(rx1.substoich, [1, 2])
     @test issetequal(rx1.prodstoich, [1])
-    @test issetequal(rx2.substoich, [1,2])
-    @test issetequal(rx2.prodstoich, [1,2])
+    @test issetequal(rx2.substoich, [1, 2])
+    @test issetequal(rx2.prodstoich, [1, 2])
 end
 
 # Checks for system with non-reaction equations.
 let
-    # Creates a reaction system and tis balanced version.
+    # Creates an unbalanced `ReactionSystem`.
     rs = @reaction_network begin
         @species O(t) H(t)
         @compounds begin
-            H2(t) = 2H
-            O2(t) = 2O
-            H2O(t) = 2H + O
+            H2(t) ~ 2H
+            O2(t) ~ 2O
+            H2O(t) ~ 2H + O
         end
-        @equation begin
+        @equations begin
             D(V) ~ 1 - V
         end
         1.0, H2 + O2 --> H2O
@@ -477,6 +477,37 @@ let
 
     # Checks that the system is correctly balanced (while not taking order of reactions or of
     # substrates/products for granted).
-    @test issetequal(reactions(rs_balanced)[1].substoich, [2,1])
+    @test issetequal(reactions(rs_balanced)[1].substoich, [2, 1])
     @test issetequal(reactions(rs_balanced)[1].prodstoich, [2])
+end
+
+# Checks that systems problematic reactions yield errors.
+let
+    # Check system with reaction with an infinite number of balanced versions.
+    rs1 = @reaction_network begin
+        @species C(t) H(t) O(t)
+        @compounds begin
+            CO ~ C + O
+            CO2 ~ C + 2O
+            H2 ~ 2H
+            CH4 ~ C + 4H
+            H2O ~ 2H + O
+        end
+        k, CO + CO2 + H2 --> CH4 + H20
+    end
+    @test_throws Exception balance_system(rs1)
+    
+    # Check system with reaction without any balanced versions.
+    rs2 = @reaction_network begin
+        @species Fe(t) S(t) O(t) H(t) N(t)
+        @compounds begin
+            FeS2 ~ Fe + 2S
+            HNO3 ~ H + N + 3O
+            Fe2S3O12 ~ 2Fe + 3S + 12O
+            NO ~ N + O
+            H2SO4 ~ 2H + S + 4O
+        end
+        k, FeS2 + HNO3 --> Fe2S3O12 + NO + H2SO4
+    end
+    @test_throws Exception balance_system(rs2)
 end