From fe2f49230ceb08724884dd4040fe8000a5363fd5 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 20 May 2024 16:26:12 -0400 Subject: [PATCH] handle weird reactions --- src/Catalyst.jl | 2 +- src/chemistry_functionality.jl | 35 ++++++++++- .../miscellaneous_tests/reaction_balancing.jl | 61 ++++++++++++++----- 3 files changed, 79 insertions(+), 19 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index b347f0ab31..fab49bd5ce 100644 --- a/src/Catalyst.jl +++ b/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 ### diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index f23b6692fe..63f2fdc182 100644 --- a/src/chemistry_functionality.jl +++ b/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 \ No newline at end of file +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 \ No newline at end of file diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index 8436e6154d..bb18512abd 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/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 \ No newline at end of file