From d854bfd01c9ccb47efbca7dea7074f68ea07abab Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 25 Nov 2024 16:16:15 -0500 Subject: [PATCH 1/5] init --- ext/CatalystGraphMakieExtension.jl | 3 +- .../rn_graph_plot.jl | 124 +++++++++++------- 2 files changed, 81 insertions(+), 46 deletions(-) diff --git a/ext/CatalystGraphMakieExtension.jl b/ext/CatalystGraphMakieExtension.jl index a47ef7f5e5..d924c2fc0b 100644 --- a/ext/CatalystGraphMakieExtension.jl +++ b/ext/CatalystGraphMakieExtension.jl @@ -1,9 +1,10 @@ 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 +import SparseArrays: rowvals, nonzeros, nzrange # Creates and exports graph plotting functions. include("CatalystGraphMakieExtension/graph_makie_extension_spatial_modelling.jl") diff --git a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl index 64e3afe50b..3682611bee 100644 --- a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl +++ b/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,33 +36,46 @@ 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.25) edgelist = edges(g) distances = zeros(length(edgelist)) - for i in 2:Base.length(edgelist) - edgelist[i] == edgelist[i-1] && (distances[i] = inc) + for (i,e) in enumerate(edgelist) + fwd = findall(==(e), edgelist); rev = findall(==(reverse(e)), edgelist) + distances[i] != 0 && continue + distances[fwd] = collect(0:inc:inc*(length(fwd)-1)) + distances[rev] = collect(inc:inc:inc*(length(rev))) end distances end +Base.copy(g::MultiGraphWrap) = MultiGraphWrap(g.g, g.multiedges, g.edgeorder) + """ plot_network(rn::ReactionSystem) @@ -77,7 +94,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 +121,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 +147,44 @@ 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 = SparseArrays.rowvals(D); vals = SparseArrays.nonzeros(D) + + # Construct the edge order for reactions. for (i, rx) in enumerate(rxs) + inds = SparseArrays.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 From a5105a7b3e279121da2dfb54b2c49eee8bb6dc4c Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 25 Nov 2024 18:03:33 -0500 Subject: [PATCH 2/5] fix --- test/extensions/graphmakie.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/extensions/graphmakie.jl b/test/extensions/graphmakie.jl index ecba9d3f7e..de5d5e34fe 100644 --- a/test/extensions/graphmakie.jl +++ b/test/extensions/graphmakie.jl @@ -74,7 +74,7 @@ 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 @@ -96,7 +96,7 @@ 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 From 4d0b9afdbf4fa5a1ca57ef00bf33b71858b0b6a8 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 26 Nov 2024 00:01:23 -0500 Subject: [PATCH 3/5] add tests --- ext/CatalystGraphMakieExtension.jl | 1 - .../rn_graph_plot.jl | 42 +++++++--- test/extensions/graphmakie.jl | 84 ++++++++++++++++++- 3 files changed, 112 insertions(+), 15 deletions(-) diff --git a/ext/CatalystGraphMakieExtension.jl b/ext/CatalystGraphMakieExtension.jl index d924c2fc0b..f6f6665813 100644 --- a/ext/CatalystGraphMakieExtension.jl +++ b/ext/CatalystGraphMakieExtension.jl @@ -4,7 +4,6 @@ module CatalystGraphMakieExtension using Catalyst, GraphMakie, Graphs, Symbolics, SparseArrays using Symbolics: get_variables! import Catalyst: species_reaction_graph, incidencematgraph, lattice_plot, lattice_animation -import SparseArrays: rowvals, nonzeros, nzrange # Creates and exports graph plotting functions. include("CatalystGraphMakieExtension/graph_makie_extension_spatial_modelling.jl") diff --git a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl index 3682611bee..af83f8b1cb 100644 --- a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl +++ b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl @@ -62,20 +62,38 @@ function Graphs.edges(g::MultiGraphWrap) edgelist = vcat(collect(Graphs.edges(g.g)), g.multiedges)[g.edgeorder] end -function gen_distances(g::MultiGraphWrap; inc = 0.25) +function gen_distances(g::MultiGraphWrap; inc = 0.4) edgelist = edges(g) distances = zeros(length(edgelist)) - for (i,e) in enumerate(edgelist) - fwd = findall(==(e), edgelist); rev = findall(==(reverse(e)), edgelist) - distances[i] != 0 && continue - distances[fwd] = collect(0:inc:inc*(length(fwd)-1)) - distances[rev] = collect(inc:inc:inc*(length(rev))) + edgedict = Dict(edgelist[1] => [1]) + for (i, e) in enumerate(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 -Base.copy(g::MultiGraphWrap) = MultiGraphWrap(g.g, g.multiedges, g.edgeorder) - """ plot_network(rn::ReactionSystem) @@ -154,12 +172,14 @@ function Catalyst.plot_complexes(rn::ReactionSystem) deps = Set() edgelist = Vector{Graphs.SimpleEdge{Int}}() - rows = SparseArrays.rowvals(D); vals = SparseArrays.nonzeros(D) + rows = rowvals(D) + vals = nonzeros(D) # Construct the edge order for reactions. for (i, rx) in enumerate(rxs) - inds = SparseArrays.nzrange(D, i) - val = vals[inds]; row = rows[inds] + 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)) diff --git a/test/extensions/graphmakie.jl b/test/extensions/graphmakie.jl index de5d5e34fe..a7fe1d1487 100644 --- a/test/extensions/graphmakie.jl +++ b/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 @@ -77,8 +77,8 @@ let 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 @@ -103,3 +103,81 @@ let @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 From f607e44129eb38f6b1904b1a06ac6beb682010ca Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 26 Nov 2024 14:08:25 -0500 Subject: [PATCH 4/5] comment out doc --- docs/src/steady_state_functionality/nonlinear_solve.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/steady_state_functionality/nonlinear_solve.md b/docs/src/steady_state_functionality/nonlinear_solve.md index 4a84198ee9..fb56fc0d3b 100644 --- a/docs/src/steady_state_functionality/nonlinear_solve.md +++ b/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 + 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 + ## [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. From 5ec4511dbcf9c9aedb9f4209159ccc68031169b3 Mon Sep 17 00:00:00 2001 From: Vincent Du <54586336+vyudu@users.noreply.github.com> Date: Wed, 27 Nov 2024 10:09:08 -0500 Subject: [PATCH 5/5] Update ext/CatalystGraphMakieExtension/rn_graph_plot.jl Co-authored-by: Sam Isaacson --- ext/CatalystGraphMakieExtension/rn_graph_plot.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl index af83f8b1cb..4420124907 100644 --- a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl +++ b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl @@ -66,7 +66,7 @@ function gen_distances(g::MultiGraphWrap; inc = 0.4) edgelist = edges(g) distances = zeros(length(edgelist)) edgedict = Dict(edgelist[1] => [1]) - for (i, e) in enumerate(edgelist[2:end]) + for (i, e) in enumerate(@view edgelist[2:end]) if edgelist[i] != edgelist[i+1] edgedict[e] = [i+1] else