Merge pull request #851 from SciML/full_system_balancing

Full system balancing

Author: TorkelE

docs/src/model_creation/chemistry_related_functionality.md

@@ -136,4 +136,23 @@ balanced_reaction = balance_reaction(unbalanced_reaction)[1]
 Reactions declared as a part of a `ReactionSystem` (e.g. using the DSL) can be retrieved for balancing using the [`reactions`](@ref) function. Please note that balancing these will not mutate the `ReactionSystem`, but a new reaction system will need to be created using the balanced reactions.
 
 !!! note
-    Reaction balancing is currently not supported for reactions involving compounds of compounds.
+    Reaction balancing is currently not supported for reactions involving compounds of compounds.
+
+### Balancing full systems
+It is possible to balance all the reactions of a reaction system simultaneously using the `balance_system` function. Here, the output is a new system, where all reactions are balanced. E.g. We can use it to balance this system of methane formation/combustion:
+```@example chem2
+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
+    end
+    1.0, C + H2 --> CH4
+    2.0, CH4 + O2 --> CO2 + H2O
+end
+rs_balanced = balance_system(rs)
+```
+Except for the modified reaction stoichiometries, the new system is identical to the previous one.

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

@@ -259,7 +259,7 @@ function balance_reaction(reaction::Reaction)
     end
 
     isempty(balancedrxs) && (@warn "Unable to balance reaction.")
-    (length(balancedrxs) > 1) && (@warn "Infinite balanced reactions from ($reaction) are possible, returning a basis for them. Note that we do not check if they preserve the set of substrates and products from the original reaction.")
+    (length(balancedrxs) > 1) && (@warn "The space of possible balanced versions of the reaction ($reaction) is greater than one-dimension. This prevents the selection of a single appropriate balanced reaction. Instead, a basis for balanced reactions is returned. Note that we do not check if they preserve the set of substrates and products from the original reaction.")
     return balancedrxs
 end
 
@@ -346,3 +346,45 @@ function create_matrix(reaction::Catalyst.Reaction)
 
     return A
 end
+
+"""
+    balance_system(rs::ReactionSystem)
+
+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 = 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
+
+# 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("The space of possible balanced versions of the reaction ($reaction) is greater than one-dimension. This prevents the selection of a single appropriate balanced reaction. No method to, in this case, automatically generate 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

@@ -375,7 +375,7 @@ let
 
     rx = Reaction(1.0,[CO,CO2,H2],[CH4,H2O])
     brxs = balance_reaction(rx)
-    @test_logs (:warn, r"Infinite balanced reactions*") match_mode=:any balance_reaction(rx)
+    @test_logs (:warn, r"The space of possible balanced *") match_mode=:any balance_reaction(rx)
     @test length(brxs) == 2
 end
 
@@ -425,5 +425,89 @@ let
     @test isequal([rn.CO2,  rn.H2O], brxs.substrates)
     @test isequal([rn.C6H12O6,  rn.O2], brxs.products)
     @test isequal([6, 6], brxs.substoich)
-    @test isequal([1,  6], brxs.prodstoich)    
+    @test isequal([1, 6], brxs.prodstoich)    
+end
+
+
+### Reaction System Balancing ###
+
+# Simple example with multiple reactions.
+let
+    # 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
+        end
+        1.0, C + H2 --> CH4
+        2.0, CH4 + O2 --> CO2 + H2O
+    end
+    rs_balanced = balance_system(rs)
+
+    # Checks that the system is correctly balanced (while not taking order of reactions or of
+    # 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.prodstoich, [1])
+    @test issetequal(rx2.substoich, [1, 2])
+    @test issetequal(rx2.prodstoich, [1, 2])
+end
+
+# Checks for system with non-reaction equations.
+let
+    # 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
+        end
+        @equations begin
+            D(V) ~ 1 - V
+        end
+        1.0, H2 + O2 --> H2O
+    end
+    rs_balanced = balance_system(rs)
+
+    # 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].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