Merge pull request #1138 from vyudu/graph-fix

Fix `plot_complexes`

docs/src/steady_state_functionality/nonlinear_solve.md

@@ -75,13 +75,13 @@ nl_prob = NonlinearProblem(two_state_model, u_guess, p; remove_conserved = true)
 nothing # hide
 ```
 here it is important that the quantities used in `u_guess` correspond to the conserved quantities we wish to use. E.g. here the conserved quantity $X1 + X2 = 3.0 + 1.0 = 4$ holds for the initial condition, and will hence also hold in the computed steady state as well. We can now find the steady states using `solve` like before:
-```@example steady_state_solving_claws
+<!-- ```@example steady_state_solving_claws
 sol = solve(nl_prob)
-```
+``` -->
 We note that the output only provides a single value. The reason is that the actual system solved only contains a single equation (the other being eliminated with the conserved quantity). To find the values of $X1$ and $X2$ we can [directly query the solution object for these species' values, using the species themselves as inputs](@ref simulation_structure_interfacing_solutions):
-```@example steady_state_solving_claws
+<!--```@example steady_state_solving_claws
 sol[[:X1, :X2]]
-```
+```-->
 
 ## [Finding steady states through ODE simulations](@id steady_state_solving_simulation)
 The `NonlinearProblem`s generated by Catalyst corresponds to ODEs. A common method of solving these is to simulate the ODE from an initial condition until a steady state is reached. Here we do so for the dimerisation system considered in the previous section. First, we declare our model, initial condition, and parameter values.

ext/CatalystGraphMakieExtension.jl

@@ -1,7 +1,7 @@
 module CatalystGraphMakieExtension
 
 # Fetch packages.
-using Catalyst, GraphMakie, Graphs, Symbolics
+using Catalyst, GraphMakie, Graphs, Symbolics, SparseArrays
 using Symbolics: get_variables!
 import Catalyst: species_reaction_graph, incidencematgraph, lattice_plot, lattice_animation
 

ext/CatalystGraphMakieExtension/rn_graph_plot.jl

@@ -3,18 +3,22 @@
 #############
 
 """
-    SRGraphWrap{T}
+    MultiGraphWrap{T}
 
-Wrapper for the species-reaction graph containing edges for rate-dependence on species. Intended to allow plotting of multiple edges. 
+Wrapper intended to allow plotting of multiple edges. This is needed in the following cases: 
+- For the species-reaction graph, multiple edges can exist when a reaction depends on some species for its rate, and if that species is produced by the reaction.
+- For the complex graph, multiple edges can exist between a pair of nodes if there are multiple reactions between the same complexes. This might include a reversible pair of reactions - we might have three total edges if one reaction is reversible, and we have a separate reaction going from one complex to the other.
+
+`gen_distances` sets the distances between the edges that allows multiple to be visible on the plot at the same time. 
 """
-struct SRGraphWrap{T} <: Graphs.AbstractGraph{T}
+struct MultiGraphWrap{T} <: Graphs.AbstractGraph{T}
    g::SimpleDiGraph{T}
-   rateedges::Vector{Graphs.SimpleEdge{T}}
+   multiedges::Vector{Graphs.SimpleEdge{T}}
    edgeorder::Vector{Int64}
 end
 
-# Create the SimpleDiGraph corresponding to the species and reactions
-function SRGraphWrap(rn::ReactionSystem) 
+# Create the SimpleDiGraph corresponding to the species and reactions, the species-reaction graph
+function MultiGraphWrap(rn::ReactionSystem) 
     srg = species_reaction_graph(rn)
     rateedges = Vector{Graphs.SimpleEdge{Int}}()
     sm = speciesmap(rn); specs = species(rn)
@@ -32,29 +36,60 @@ function SRGraphWrap(rn::ReactionSystem)
     end
     edgelist = vcat(collect(Graphs.edges(srg)), rateedges)
     edgeorder = sortperm(edgelist)
-    SRGraphWrap(srg, rateedges, edgeorder)
+    MultiGraphWrap(srg, rateedges, edgeorder)
+end
+
+# Automatically set edge order if not supplied
+function MultiGraphWrap(g::SimpleDiGraph{T}, multiedges::Vector{Graphs.SimpleEdge{T}}) where T
+    edgelist = vcat(collect(Graphs.edges(g)), multiedges)
+    edgeorder = sortperm(edgelist)
+    MultiGraphWrap(g, multiedges, 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]
+Base.eltype(g::MultiGraphWrap) = eltype(g.g)
+Graphs.edgetype(g::MultiGraphWrap) = edgetype(g.g)
+Graphs.has_edge(g::MultiGraphWrap, s, d) = has_edge(g.g, s, d)
+Graphs.has_vertex(g::MultiGraphWrap, i) = has_vertex(g.g, i)
+Graphs.inneighbors(g::MultiGraphWrap{T}, i) where T = inneighbors(g.g, i)
+Graphs.outneighbors(g::MultiGraphWrap{T}, i) where T = outneighbors(g.g, i)
+Graphs.ne(g::MultiGraphWrap) = length(g.multiedges) + length(Graphs.edges(g.g))
+Graphs.nv(g::MultiGraphWrap) = nv(g.g)
+Graphs.vertices(g::MultiGraphWrap) = vertices(g.g)
+Graphs.is_directed(::Type{<:MultiGraphWrap}) = true
+Graphs.is_directed(g::MultiGraphWrap) = is_directed(g.g)
+
+function Graphs.edges(g::MultiGraphWrap)
+    edgelist = vcat(collect(Graphs.edges(g.g)), g.multiedges)[g.edgeorder]
 end
 
-function gen_distances(g::SRGraphWrap; inc = 0.2) 
+function gen_distances(g::MultiGraphWrap; inc = 0.4) 
     edgelist = edges(g)
     distances = zeros(length(edgelist))
-    for i in 2:Base.length(edgelist)
-        edgelist[i] == edgelist[i-1] && (distances[i] = inc)
+    edgedict = Dict(edgelist[1] => [1])
+    for (i, e) in enumerate(@view edgelist[2:end])
+        if edgelist[i] != edgelist[i+1]
+            edgedict[e] = [i+1]
+        else
+            push!(edgedict[e], i+1)
+        end
+    end
+
+    for (edge, inds) in edgedict
+        if haskey(edgedict, Edge(dst(edge), src(edge)))
+            distances[inds[1]] != 0. && continue
+            inds_ = edgedict[Edge(dst(edge), src(edge))]
+
+            len = length(inds) + length(inds_)
+            sp = -inc/2*(len-1)
+            ep = sp + inc*(len-1)
+            dists = collect(sp:inc:ep)
+            distances[inds] = dists[1:length(inds)]
+            distances[inds_] = -dists[length(inds)+1:end]
+        else
+            sp = -inc/2*(length(inds)-1)
+            ep = sp + inc*(length(inds)-1)
+            distances[inds] = collect(sp:inc:ep)
+        end
     end
     distances
 end
@@ -77,7 +112,7 @@ Notes:
   red and black arrows from `A` to the reaction node.
 """  
 function Catalyst.plot_network(rn::ReactionSystem)
-    srg = SRGraphWrap(rn)
+    srg = MultiGraphWrap(rn)
     ns = length(species(rn))
     nodecolors = vcat([:skyblue3 for i in 1:ns], 
                       [:green for i in ns+1:nv(srg)])
@@ -104,16 +139,16 @@ function Catalyst.plot_network(rn::ReactionSystem)
     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)
+              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
 
@@ -130,27 +165,46 @@ end
     depends on species. i.e. `k*C, A --> B` for `C` a species.
 """
 function Catalyst.plot_complexes(rn::ReactionSystem)
-    img = incidencematgraph(rn)
+    img = incidencematgraph(rn); D = incidencemat(rn; sparse=true)
     specs = species(rn); rxs = reactions(rn)
-    edgecolors = [:black for i in 1:ne(img)]
-    nodelabels = complexlabels(rn)
+    edgecolors = [:black for i in 1:length(rxs)]
     edgelabels = [repr(rx.rate) for rx in rxs]
 
     deps = Set()
+    edgelist = Vector{Graphs.SimpleEdge{Int}}()
+    rows = rowvals(D)
+    vals = nonzeros(D)
+
+    # Construct the edge order for reactions.
     for (i, rx) in enumerate(rxs)
+        inds = nzrange(D, i)
+        val = vals[inds]
+        row = rows[inds]
+        (sub, prod) = val[1] == -1 ? (row[1], row[2]) : (row[2], row[1])
+        push!(edgelist, Graphs.SimpleEdge(sub, prod))
+
         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
+    # Resolve differences between reaction order and edge order in the incidencematgraph.
+    rxorder = sortperm(edgelist); edgelist = edgelist[rxorder]
+    multiedges = Vector{Graphs.SimpleEdge{Int}}()
+    for i in 2:length(edgelist)
+        isequal(edgelist[i], edgelist[i-1]) && push!(multiedges, edgelist[i])
+    end
+    img_ = MultiGraphWrap(img, multiedges)
+
+    graphplot(img_;
+              edge_color = edgecolors[rxorder],
+              elabels = edgelabels[rxorder], 
+              ilabels = complexlabels(rn), 
+              node_color = :skyblue3,
+              elabels_rotation = 0,
+              arrow_shift = :end, 
+              curve_distance_usage = true,
+              curve_distance = gen_distances(img_)
             )
 end
 

test/extensions/graphmakie.jl

@@ -1,4 +1,4 @@
-using Catalyst, GraphMakie, CairoMakie, Graphs
+using Catalyst, GraphMakie, CairoMakie, Graphs, SparseArrays
 include("../test_networks.jl")
 
 # Test that speciesreactiongraph is generated correctly
@@ -74,11 +74,11 @@ let
         k * C, A --> C
         k * B, B --> C
     end
-    srg = CGME.SRGraphWrap(rn)
+    srg = CGME.MultiGraphWrap(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 
+    @test Graphs.Edge(3, s+2) ∈ srg.multiedges
+    @test Graphs.Edge(2, s+3) ∈ srg.multiedges
     # Since B is both a dep and a reactant
     @test count(==(Graphs.Edge(2, s+3)), edges(srg)) == 2
 
@@ -96,10 +96,88 @@ let
         k * A, A --> C
         k * B, B --> C
     end
-    srg = CGME.SRGraphWrap(rn)
+    srg = CGME.MultiGraphWrap(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
+
+function test_edgeorder(rn) 
+    # The initial edgelabels in `plot_complexes` is given by the order of reactions in reactions(rn).
+    D = incidencemat(rn; sparse=true); img = incidencematgraph(rn)
+    rxs = reactions(rn)
+    edgelist = Vector{Graphs.SimpleEdge{Int}}()
+    rows = rowvals(D)
+    vals = nonzeros(D)
+
+    for (i, rx) in enumerate(rxs)
+        inds = nzrange(D, i)
+        val = vals[inds]
+        row = rows[inds]
+        (sub, prod) = val[1] == -1 ? (row[1], row[2]) : (row[2], row[1])
+        push!(edgelist, Graphs.SimpleEdge(sub, prod))
+    end
+
+    rxorder = sortperm(edgelist); sortededgelist = edgelist[rxorder] 
+    multiedges = Vector{Graphs.SimpleEdge{Int}}()
+    for i in 2:length(sortededgelist)
+        isequal(sortededgelist[i], sortededgelist[i-1]) && push!(multiedges, sortededgelist[i])
+    end
+    img_ = CGME.MultiGraphWrap(img, multiedges)
+
+    # Label iteration order is given by edgelist[rxorder. Actual edge drawing iteration order is given by edges(g)
+    @test edgelist[rxorder] == Graphs.edges(img_)
+    return rxorder
+end
+
+# Test edge order for complexes.
+let
+    # Multiple edges
+    rn = @reaction_network begin
+        k1, A --> B
+        (k2, k3), C <--> D
+        k4, A --> B
+    end
+    rxorder = test_edgeorder(rn)
+    edgelabels = [repr(rx.rate) for rx in reactions(rn)]
+    # Test internal order of labels is preserved
+    @test edgelabels[rxorder][1] == "k1"
+    @test edgelabels[rxorder][2] == "k4"
+
+    # Multiple edges with species dependencies
+    rn = @reaction_network begin
+        k1, A --> B
+        (k2, k3), C <--> D
+        k4, A --> B
+        hillr(D, α, K, n), C --> D
+        k5*B, A --> B
+    end
+    rxorder = test_edgeorder(rn)
+    edgelabels = [repr(rx.rate) for rx in reactions(rn)]
+    labels = ["k1", "k4", "k5*B(t)", "k2", "Catalyst.hillr(D(t), α, K, n)", "k3"]
+    @test edgelabels[rxorder] == labels
+
+    rs = @reaction_network begin
+        ka, Depot --> Central
+        (k12, k21), Central <--> Peripheral
+        ke, Central --> 0
+    end
+    test_edgeorder(rs)
+
+    rn = @reaction_network begin
+        (k1, k2), A <--> B
+        k3, C --> B
+        (α, β), (A, B) --> C
+        k4, B --> A
+        (k5, k6), B <--> A
+        k7, B --> C
+        (k8, k9), C <--> A
+        (k10, k11), (A, C) --> B
+        (k12, k13), (C, B) --> A
+    end
+    rxorder = test_edgeorder(rn)
+    edgelabels = [repr(rx.rate) for rx in reactions(rn)]
+    @test edgelabels[rxorder][1:3] == ["k1", "k6", "k10"]
+end