up

Author: vyudu

ext/CatalystGraphMakieExtension/rn_graph_plot.jl

@@ -12,6 +12,26 @@ struct SRGraphWrap{T} <: AbstractGraph{T}
    rateedges::Vector{SimpleEdge}
 end
 
+# Create the SimpleDiGraph corresponding to the species and reactions
+function SRGraphWrap(rn::ReactionSystem) 
+    srg = speciesreactiongraph(rn)
+    rateedges = Vector{SimpleEdge}()
+    sm = speciesmap(rn); specs = species(rn)
+
+    deps = Set()
+    for (i, rx) in enumerate(reactions(rn)) 
+        empty!(deps)
+        get_variables!(deps, rx.rate, specs)
+        if !isempty(deps)
+            for spec in deps 
+                specidx = sm[spec]
+                push!(rateedges, SimpleEdge(specidx, i + length(specs)))
+            end
+        end
+    end
+    SRGraphWrap(srg, rateedges)
+end
+
 Base.eltype(g::SRGraphWrap) = eltype(g.g)
 Graphs.edgetype(g::SRGraphWrap) = edgetype(g.g)
 Graphs.has_edge(g::SRGraphWrap, s, d) = has_edge(g.g, s, d)
@@ -25,8 +45,7 @@ Graphs.is_directed(g::SRGraphWrap) = is_directed(g.g)
 
 function Graphs.edges(g::SRGraphWrap)
     edgelist = vcat(collect(Graphs.edges(g.g)), g.rateedges)
-    edgeorder = sortperm(edgelist)
-    edgelist = edgelist[edgeorder]
+    edgelist = sort!(edgelist)
 end
 
 function gen_distances(g::SRGraphWrap; inc = 0.2) 
@@ -114,26 +133,6 @@ function plot_speciesreaction_graph(rn::ReactionSystem; interactive = false)
     f
 end
 
-# Create the SimpleDiGraph corresponding to the species and reactions
-function SRGraphWrap(rn::ReactionSystem) 
-    srg = speciesreactiongraph(rn)
-    rateedges = Vector{SimpleEdge}()
-    sm = speciesmap(rn); specs = species(rn)
-
-    deps = Set()
-    for (i, rx) in enumerate(reactions(rn)) 
-        empty!(deps)
-        get_variables!(deps, rx.rate, specs)
-        if !isempty(deps)
-            for spec in deps 
-                specidx = sm[spec]
-                push!(rateedges, SimpleEdge(specidx, i + length(specs)))
-            end
-        end
-    end
-    SRGraphWrap(srg, rateedges)
-end
-
 """
     ComplexGraph(rn::ReactionSystem; interactive=false)
 

src/network_analysis.jl

@@ -316,37 +316,37 @@ end
 
 
 """
-    speciesreactiongraph(rn::ReactionSystem)
+    species_reaction_graph(rn::ReactionSystem)
 
 Construct a directed simple graph where there are two types of nodes: species and reactions. 
-An edge from a species to reaction indicates that it is a reactant, and an edge from a reaction
-to a species indicates that it is a product. By default, the species vertices are listed
-first, so the first n indices correspond to species nodes. 
+An edge from a species *s* to reaction *r* indicates that *s* is a reactant for *r*, and an edge from a reaction
+r to a species s indicates that *s* is a product of *r*. By default, the species vertices are listed
+first, so the first *n* indices correspond to species nodes. 
 
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-speciesreactiongraph(sir)
+species_reaction_graph(sir)
 """
 function species_reaction_graph(rn::ReactionSystem) 
     specs = species(rn)
     rxs = reactions(rn)
     sm = speciesmap(rn)
-    s = length(specs); r = length(rxs); nv = s + r
+    s = length(specs)
 
-    adjmat = zeros(Int64, nv, nv)
+    edgelist = Graphs.Edge[]
     for (i, rx) in enumerate(rxs) 
-        for (spec, stoich) in zip(rx.substrates, rx.substoich)
-            adjmat[sm[spec], s+i] = stoich
+        for spec in rx.substrates
+            push!(edgelist, Graphs.Edge(sm[spec], s+i))
         end
-        for (spec, stoich) in zip(rx.products, rx.prodstoich)
-            adjmat[s+i, sm[spec]] = stoich
+        for spec in rx.products
+            push!(edgelist, Graphs.Edge(s+i, sm[spec]))
         end
     end
-    srg = SimpleDiGraph(adjmat)
+    srg = Graphs.SimpleDiGraphFromIterator(edgelist)
 end
 
 ### Linkage, Deficiency, Reversibility ###

test/extensions/graphmakie.jl

@@ -0,0 +1,63 @@
+using Catalyst, GraphMakie, GLMakie
+include("../test_networks.jl")
+# Test that speciesreactiongraph is generated correctly
+let
+    brusselator = @reaction_network begin
+        A, ∅ --> X
+        1, 2X + Y --> 3X
+        B, X --> Y
+        1, X --> ∅
+    end
+
+    srg = Catalyst.species_reaction_graph(brusselator)
+    s = length(species(brusselator))
+    edgel = Edge.([(s+1, 1),
+                   (1, s+2),
+                   (2, s+2),
+                   (s+2, 1),
+                   (s+3, 2),
+                   (1, s+3),
+                   (1, s+4)])
+    @test all(∈(collect(Graphs.edges(srg))), edgel)
+
+    MAPK = @reaction_network MAPK begin
+        (k₁, k₂),KKK + E1 <--> KKKE1
+        k₃, KKKE1 --> KKK_ + E1
+        (k₄, k₅), KKK_ + E2 <--> KKKE2
+        k₆, KKKE2 --> KKK + E2
+        (k₇, k₈), KK + KKK_ <--> KK_KKK_
+        k₉, KK_KKK_ --> KKP + KKK_
+        (k₁₀, k₁₁), KKP + KKK_ <--> KKPKKK_
+        k₁₂, KKPKKK_ --> KKPP + KKK_
+        (k₁₃, k₁₄), KKP + KKPase <--> KKPKKPase
+        k₁₅, KKPPKKPase --> KKP + KKPase
+        k₁₆,KKPKKPase --> KK + KKPase
+        (k₁₇, k₁₈), KKPP + KKPase <--> KKPPKKPase
+        (k₁₉, k₂₀), KKPP + K <--> KKPPK
+        k₂₁, KKPPK --> KKPP + KP
+        (k₂₂, k₂₃), KKPP + KP <--> KPKKPP
+        k₂₄, KPKKPP --> KPP + KKPP
+        (k₂₅, k₂₆), KP + KPase <--> KPKPase
+        k₂₇, KKPPKPase --> KP + KPase
+        k₂₈, KPKPase --> K + KPase
+        (k₂₉, k₃₀), KPP + KPase <--> KKPPKPase
+    end
+    srg = Catalyst.species_reaction_graph(MAPK)
+    @test nv(srg) == length(species(MAPK)) + length(reactions(MAPK))
+    @test ne(srg) == 90
+end
+
+# Test that rate edges are inferred correctly
+let
+    rn = @reaction_network begin
+        k, A --> B
+        k * C, A --> C
+        k * B, B --> C
+    end
+    srg = SRGraphWrap(rn)
+    s = length(species(rn))
+    @test Edge(3, s+2) ∈ srg.rateedges 
+    @test Edge(2, s+3) ∈ srg.rateedges 
+    # Since B is both a dep and a reactant
+    @test count(==(Edge(2, s+3)), edges(srg)) == 2
+end