From a4e8b6516214d8521f34e899d9a2dee25ca65a17 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 19 Jun 2024 12:17:01 -0400 Subject: [PATCH 1/8] JuliaFormatter --- src/Catalyst.jl | 3 ++- src/network_analysis.jl | 8 ++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index cec9cca52f..84b1a31b71 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -127,7 +127,8 @@ export @reaction_network, @network_component, @reaction, @species include("network_analysis.jl") export reactioncomplexmap, reactioncomplexes, incidencemat export complexstoichmat -export complexoutgoingmat, incidencematgraph, linkageclasses, stronglinkageclasses, terminallinkageclasses, deficiency, subnetworks +export complexoutgoingmat, incidencematgraph, linkageclasses, stronglinkageclasses, + terminallinkageclasses, deficiency, subnetworks export linkagedeficiencies, isreversible, isweaklyreversible export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants diff --git a/src/network_analysis.jl b/src/network_analysis.jl index b6f2954ce0..4db0913bf7 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -365,7 +365,7 @@ linkageclasses(incidencegraph) = Graphs.connected_components(incidencegraph) Return the strongly connected components of a reaction network's incidence graph (i.e. sub-groups of reaction complexes such that every complex is reachable from every other one in the sub-group). """ -function stronglinkageclasses(rn::ReactionSystem) +function stronglinkageclasses(rn::ReactionSystem) nps = get_networkproperties(rn) if isempty(nps.stronglinkageclasses) nps.stronglinkageclasses = stronglinkageclasses(incidencematgraph(rn)) @@ -381,17 +381,17 @@ stronglinkageclasses(incidencegraph) = Graphs.strongly_connected_components(inci Return the terminal strongly connected components of a reaction network's incidence graph (i.e. sub-groups of reaction complexes that are 1) strongly connected and 2) every reaction in the component produces a complex in the component). """ -function terminallinkageclasses(rn::ReactionSystem) +function terminallinkageclasses(rn::ReactionSystem) nps = get_networkproperties(rn) if isempty(nps.terminallinkageclasses) slcs = stronglinkageclasses(rn) - tslcs = filter(lc->isterminal(lc, rn), slcs) + tslcs = filter(lc -> isterminal(lc, rn), slcs) nps.terminallinkageclasses = tslcs end nps.terminallinkageclasses end -function isterminal(lc::Vector, rn::ReactionSystem) +function isterminal(lc::Vector, rn::ReactionSystem) imat = incidencemat(rn) for r in 1:size(imat, 2) From 9b868d04054e3479297cf975e2d2860a21c3ce97 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 25 Jun 2024 13:20:44 -0400 Subject: [PATCH 2/8] commented code --- src/network_analysis.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 4db0913bf7..b8166c0360 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -378,7 +378,7 @@ stronglinkageclasses(incidencegraph) = Graphs.strongly_connected_components(inci """ terminallinkageclasses(rn::ReactionSystem) - Return the terminal strongly connected components of a reaction network's incidence graph (i.e. sub-groups of reaction complexes that are 1) strongly connected and 2) every reaction in the component produces a complex in the component). + Return the terminal strongly connected components of a reaction network's incidence graph (i.e. sub-groups of reaction complexes that are 1) strongly connected and 2) every outgoing reaction from a complex in the component produces a complex also in the component). """ function terminallinkageclasses(rn::ReactionSystem) @@ -391,11 +391,16 @@ function terminallinkageclasses(rn::ReactionSystem) nps.terminallinkageclasses end + +# Check whether a given linkage class in a reaction network is terminal, i.e. all outgoing reactions from complexes in the linkage class produce a complex also in hte linkage class function isterminal(lc::Vector, rn::ReactionSystem) imat = incidencemat(rn) for r in 1:size(imat, 2) + # Find the index of the reactant complex for a given reaction s = findfirst(==(-1), @view imat[:, r]) + + # If the reactant complex is in the linkage class, check whether the product complex is also in the linkage class. If any of them are not, return false. if s in Set(lc) p = findfirst(==(1), @view imat[:, r]) p in Set(lc) ? continue : return false From 2728402e415109b10d00c5247f8b622a9ed607a1 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 25 Jun 2024 13:32:10 -0400 Subject: [PATCH 3/8] commented tests --- test/network_analysis/network_properties.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index baee75e32a..811d6418cc 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -328,6 +328,8 @@ end ### STRONG LINKAGE CLASS TESTS + +# a) Checks that strong/terminal linkage classes are correctly found. Should identify the (A, B+C) linkage class as non-terminal, since B + C produces D let rn = @reaction_network begin (k1, k2), A <--> B + C @@ -347,6 +349,7 @@ let @test issubset([[3,4,5], [6,7]], tslcs) end +# b) Makes the D + E --> G reaction irreversible. Thus, (D+E) becomes a non-terminal linkage class. Checks whether correctly identifies both (A, B+C) and (D+E) as non-terminal let rn = @reaction_network begin (k1, k2), A <--> B + C @@ -366,6 +369,7 @@ let @test issubset([[3,4,5], [7]], tslcs) end +# From a), makes the B + C <--> D reaction reversible. Thus, the non-terminal (A, B+C) linkage class gets absorbed into the terminal (A, B+C, D, E, 2F) linkage class, and the terminal linkage classes and strong linkage classes coincide. let rn = @reaction_network begin (k1, k2), A <--> B + C @@ -385,6 +389,7 @@ let @test issubset([[1,2,3,4,5], [6,7]], tslcs) end +# Simple test for strong and terminal linkage classes let rn = @reaction_network begin (k1, k2), A <--> 2B From 9544511e230c2127e02f1d47828003d2cc98367c Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 11 Jul 2024 11:00:15 -0400 Subject: [PATCH 4/8] up --- src/network_analysis.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index b8166c0360..621bbe8977 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -392,7 +392,8 @@ function terminallinkageclasses(rn::ReactionSystem) end -# Check whether a given linkage class in a reaction network is terminal, i.e. all outgoing reactions from complexes in the linkage class produce a complex also in hte linkage class +# Helper function for terminallinkageclasses. Given a linkage class and a reaction network, say whether the linkage class is terminal, +# i.e. all outgoing reactions from complexes in the linkage class produce a complex also in the linkage class function isterminal(lc::Vector, rn::ReactionSystem) imat = incidencemat(rn) From ed8fda6ca8fddbe7ec0d992a04566af076cfd109 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 11 Jul 2024 11:07:00 -0400 Subject: [PATCH 5/8] formatter --- src/network_analysis.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 621bbe8977..03c648edd8 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -391,7 +391,6 @@ function terminallinkageclasses(rn::ReactionSystem) nps.terminallinkageclasses end - # Helper function for terminallinkageclasses. Given a linkage class and a reaction network, say whether the linkage class is terminal, # i.e. all outgoing reactions from complexes in the linkage class produce a complex also in the linkage class function isterminal(lc::Vector, rn::ReactionSystem) From 66277600a2f000eafe7365f3a5cb13fad8188308 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 11 Jul 2024 16:11:35 -0400 Subject: [PATCH 6/8] fixed reset --- src/reactionsystem.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index c45231f502..5944e895fc 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -119,6 +119,7 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V} nps.isempty && return nps.netstoichmat = Matrix{Int}(undef, 0, 0) nps.conservationmat = Matrix{I}(undef, 0, 0) + nps.cyclemat = Matrix{Int}(undef, 0, 0) empty!(nps.col_order) nps.rank = 0 nps.nullity = 0 @@ -134,6 +135,8 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V} nps.complexoutgoingmat = Matrix{Int}(undef, 0, 0) nps.incidencegraph = Graphs.DiGraph() empty!(nps.linkageclasses) + empty!(nps.stronglinkageclasses) + empty!(nps.terminallinkageclasses) nps.deficiency = 0 # this needs to be last due to setproperty! setting it to false From 1b3cdd7f6d3aad240a6c89ec914ad39c95615141 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 12 Jul 2024 12:03:26 -0400 Subject: [PATCH 7/8] NPS fix --- src/reactionsystem.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index ffee08fcb9..9a24863bb8 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -78,6 +78,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re isempty::Bool = true netstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) conservationmat::Matrix{I} = Matrix{I}(undef, 0, 0) + cyclemat::Matrix{I} = Matrix{I}(undef, 0, 0) col_order::Vector{Int} = Int[] rank::Int = 0 nullity::Int = 0 From 28be9ebc5a3a5b64a18c2339b98f85e3cfe1eab0 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 12 Jul 2024 12:57:43 -0400 Subject: [PATCH 8/8] NPS fix --- src/reactionsystem.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 9a24863bb8..7bd09f555b 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -94,6 +94,8 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re complexoutgoingmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) incidencegraph::Graphs.SimpleDiGraph{Int} = Graphs.DiGraph() linkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0) + stronglinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0) + terminallinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0) deficiency::Int = 0 end #! format: on