Merge pull request #1115 from vyudu/GraphMakie

GraphMakie update

ext/CatalystGraphMakieExtension.jl

@@ -1,11 +1,11 @@
 module CatalystGraphMakieExtension
 
 # Fetch packages.
-using Catalyst, GraphMakie
-import Catalyst: lattice_plot, lattice_animation, extract_vals
-import Graphs: AbstractGraph, SimpleGraph
+using Catalyst, GraphMakie, Graphs, Symbolics
+using Symbolics: get_variables!
+import Catalyst: species_reaction_graph, incidencematgraph, lattice_plot, lattice_animation
 
-# Creates and exports hc_steady_states function.
+# Creates and exports graph plotting functions.
 include("CatalystGraphMakieExtension/graph_makie_extension_spatial_modelling.jl")
-
+include("CatalystGraphMakieExtension/rn_graph_plot.jl")
 end

ext/CatalystGraphMakieExtension/rn_graph_plot.jl

@@ -0,0 +1,185 @@
+#############
+# Adapted from https://github.com/MakieOrg/GraphMakie.jl/issues/52#issuecomment-1018527479
+#############
+
+"""
+    SRGraphWrap{T}
+
+Wrapper for the species-reaction graph containing edges for rate-dependence on species. Intended to allow plotting of multiple edges. 
+"""
+struct SRGraphWrap{T} <: Graphs.AbstractGraph{T}
+   g::SimpleDiGraph{T}
+   rateedges::Vector{Graphs.SimpleEdge{T}}
+   edgeorder::Vector{Int64}
+end
+
+# Create the SimpleDiGraph corresponding to the species and reactions
+function SRGraphWrap(rn::ReactionSystem) 
+    srg = species_reaction_graph(rn)
+    rateedges = Vector{Graphs.SimpleEdge{Int}}()
+    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, Graphs.SimpleEdge(specidx, i + length(specs)))
+            end
+        end
+    end
+    edgelist = vcat(collect(Graphs.edges(srg)), rateedges)
+    edgeorder = sortperm(edgelist)
+    SRGraphWrap(srg, rateedges, edgeorder)
+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)
+Graphs.has_vertex(g::SRGraphWrap, i) = has_vertex(g.g, i)
+Graphs.inneighbors(g::SRGraphWrap{T}, i) where T = inneighbors(g.g, i)
+Graphs.outneighbors(g::SRGraphWrap{T}, i) where T = outneighbors(g.g, i)
+Graphs.ne(g::SRGraphWrap) = length(g.rateedges) + length(Graphs.edges(g.g))
+Graphs.nv(g::SRGraphWrap) = nv(g.g)
+Graphs.vertices(g::SRGraphWrap) = vertices(g.g)
+Graphs.is_directed(g::SRGraphWrap) = is_directed(g.g)
+
+function Graphs.edges(g::SRGraphWrap)
+    edgelist = vcat(collect(Graphs.edges(g.g)), g.rateedges)[g.edgeorder]
+end
+
+function gen_distances(g::SRGraphWrap; inc = 0.2) 
+    edgelist = edges(g)
+    distances = zeros(length(edgelist))
+    for i in 2:Base.length(edgelist)
+        edgelist[i] == edgelist[i-1] && (distances[i] = inc)
+    end
+    distances
+end
+
+"""
+    plot_network(rn::ReactionSystem; interactive=false)
+
+Converts a [`ReactionSystem`](@ref) into a GraphMakie plot of the species reaction graph
+(or Petri net representation). Reactions correspond to small green circles, and 
+species to blue circles.
+
+Notes:
+- Black arrows from species to reactions indicate reactants, and are labelled
+  with their input stoichiometry.
+- Black arrows from reactions to species indicate products, and are labelled
+  with their output stoichiometry.
+- Red arrows from species to reactions indicate that species is used within the
+  rate expression. For example, in the reaction `k*A, B --> C`, there would be a
+  red arrow from `A` to the reaction node. In `k*A, A+B --> C`, there would be
+  red and black arrows from `A` to the reaction node.
+"""  
+# TODO: update docs for interacting with plots. The `interactive` flag sets the ability to interactively drag nodes and edges in the generated plot. Only allowed if `GLMakie` is the loaded Makie backend.
+function Catalyst.plot_network(rn::ReactionSystem)
+    srg = SRGraphWrap(rn)
+    ns = length(species(rn))
+    nodecolors = vcat([:skyblue3 for i in 1:ns], 
+                      [:green for i in ns+1:nv(srg)])
+    ilabels = vcat(map(s -> String(tosymbol(s, escape=false)), species(rn)),
+                   fill("", nv(srg.g) - ns))
+    nodesizes = vcat([30 for i in 1:ns],
+                     [10 for i in ns+1:nv(srg)])
+
+    ssm = substoichmat(rn); psm = prodstoichmat(rn)
+    # Get stoichiometry of reaction
+    edgelabels = map(Graphs.edges(srg.g)) do e
+        string(src(e) > ns ? 
+            psm[dst(e), src(e)-ns] :
+            ssm[src(e), dst(e)-ns]) 
+    end 
+    edgecolors = [:black for i in 1:ne(srg)]
+
+    num_e = ne(srg.g)
+    for i in 1:length(srg.edgeorder)
+        if srg.edgeorder[i] > num_e
+            edgecolors[i] = :red
+            insert!(edgelabels, i, "")
+        end
+    end
+
+    graphplot(srg; 
+             edge_color = edgecolors,
+             elabels = edgelabels,
+             elabels_rotation = 0,
+             ilabels = ilabels,
+             node_color = nodecolors,
+             node_size = nodesizes,
+             arrow_shift = :end,
+             arrow_size = 20,
+             curve_distance_usage = true,
+             curve_distance = gen_distances(srg)
+            )
+end
+
+"""
+    plot_complexes(rn::ReactionSystem; interactive=false)
+
+    Creates a GraphMakie plot of the [`ReactionComplex`](@ref)s in `rn`. Reactions
+    correspond to arrows and reaction complexes to blue circles.
+
+    Notes:
+    - Black arrows from complexes to complexes indicate reactions whose rate is a
+      parameter or a `Number`. i.e. `k, A --> B`.
+    - Red arrows from complexes to complexes indicate reactions whose rate
+    depends on species. i.e. `k*C, A --> B` for `C` a species.
+"""
+function Catalyst.plot_complexes(rn::ReactionSystem)
+    img = incidencematgraph(rn)
+    specs = species(rn); rxs = reactions(rn)
+    edgecolors = [:black for i in 1:ne(img)]
+    nodelabels = complexlabels(rn)
+    edgelabels = [repr(rx.rate) for rx in rxs]
+
+    deps = Set()
+    for (i, rx) in enumerate(rxs)
+        empty!(deps)
+        get_variables!(deps, rx.rate, specs)
+        (!isempty(deps)) && (edgecolors[i] = :red)
+    end
+
+    graphplot(img; 
+             edge_color = edgecolors,
+             elabels = edgelabels, 
+             ilabels = complexlabels(rn), 
+             node_color = :skyblue3,
+             elabels_rotation = 0,
+             arrow_shift = :end, 
+             curve_distance = 0.2
+            )
+end
+
+function complexelem_tostr(e::Catalyst.ReactionComplexElement, specstrs) 
+    if e.speciesstoich == 1
+        return "$(specstrs[e.speciesid])"  
+    else
+        return "$(e.speciesstoich)$(specstrs[e.speciesid])"
+    end
+end
+
+# Get the strings corresponding to the reaction complexes
+function complexlabels(rn::ReactionSystem)
+    labels = String[]
+
+    specstrs = map(s -> String(tosymbol(s, escape=false)), species(rn))
+    complexes, B = reactioncomplexes(rn)
+
+    for complex in complexes
+        if isempty(complex) 
+            push!(labels, "∅")
+        elseif length(complex) == 1
+            push!(labels, complexelem_tostr(complex[1], specstrs))
+        else
+            elems = map(c -> complexelem_tostr(c, specstrs), complex)
+            str = reduce((e1, e2) -> *(e1, " + ", e2), @view elems[2:end]; init = elems[1])
+            push!(labels, str)
+        end
+    end
+    labels
+end

src/Catalyst.jl

@@ -167,6 +167,11 @@ export hc_steady_states
 function make_si_ode end
 export make_si_ode
 
+# GraphMakie
+function plot_network end
+function plot_complexes end
+export plot_network, plot_complexes
+
 ### Spatial Reaction Networks ###
 
 # Spatial reactions.

src/network_analysis.jl

@@ -314,6 +314,44 @@ function incidencematgraph(incidencemat::SparseMatrixCSC{Int, Int})
     return graph
 end
 
+
+"""
+    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 *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. 
+
+Note: this is equivalent to the Petri net representation of a chemical reaction network.
+
+For example,
+```julia
+sir = @reaction_network SIR begin
+    β, S + I --> 2I
+    ν, I --> R
+end
+species_reaction_graph(sir)
+"""
+function species_reaction_graph(rn::ReactionSystem) 
+    specs = species(rn)
+    rxs = reactions(rn)
+    sm = speciesmap(rn)
+    s = length(specs)
+
+    edgelist = Graphs.Edge[]
+    for (i, rx) in enumerate(rxs) 
+        for spec in rx.substrates
+            push!(edgelist, Graphs.Edge(sm[spec], s+i))
+        end
+        for spec in rx.products
+            push!(edgelist, Graphs.Edge(s+i, sm[spec]))
+        end
+    end
+    srg = Graphs.SimpleDiGraphFromIterator(edgelist)
+end
+
 ### Linkage, Deficiency, Reversibility ###
 
 """

test/extensions/graphmakie.jl

@@ -0,0 +1,106 @@
+using Catalyst, GraphMakie, GLMakie, Graphs
+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 = Graphs.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
+
+    # Test that figures are generated properly. 
+    f = plot_network(MAPK)
+    save("fig.png", f)
+    @test isfile("fig.png")
+    rm("fig.png")
+    f = plot_network(brusselator)
+    save("fig.png", f)
+    @test isfile("fig.png")
+    rm("fig.png")
+
+    f = plot_complexes(MAPK); save("fig.png", f)
+    @test isfile("fig.png")
+    rm("fig.png")
+    f = plot_complexes(brusselator); save("fig.png", f)
+    @test isfile("fig.png")
+    rm("fig.png")
+end
+
+CGME = Base.get_extension(parentmodule(ReactionSystem), :CatalystGraphMakieExtension)
+# Test that rate edges are inferred correctly. We should see two for the following reaction network. 
+let
+    # Two rate edges, one to species and one to product
+    rn = @reaction_network begin
+        k, A --> B
+        k * C, A --> C
+        k * B, B --> C
+    end
+    srg = CGME.SRGraphWrap(rn)
+    s = length(species(rn))
+    @test ne(srg) == 8
+    @test Graphs.Edge(3, s+2) ∈ srg.rateedges 
+    @test Graphs.Edge(2, s+3) ∈ srg.rateedges 
+    # Since B is both a dep and a reactant
+    @test count(==(Graphs.Edge(2, s+3)), edges(srg)) == 2
+
+    f = plot_network(rn)
+    save("fig.png", f)
+    @test isfile("fig.png")
+    rm("fig.png")
+    f = plot_complexes(rn); save("fig.png", f)
+    @test isfile("fig.png")
+    rm("fig.png")
+
+    # Two rate edges, both to reactants 
+    rn = @reaction_network begin
+        k, A --> B
+        k * A, A --> C
+        k * B, B --> C
+    end
+    srg = CGME.SRGraphWrap(rn)
+    s = length(species(rn))
+    @test ne(srg) == 8
+    # Since A, B is both a dep and a reactant
+    @test count(==(Graphs.Edge(1, s+2)), edges(srg)) == 2
+    @test count(==(Graphs.Edge(2, s+3)), edges(srg)) == 2
+end
+