Merge branch 'master' into add-matrices

README.md

@@ -92,7 +92,7 @@ an ordinary differential equation.
 
 ```julia
 # Fetch required packages.
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqDefault, Plots
 
 # Create model.
 model = @reaction_network begin

ext/CatalystGraphMakieExtension/rn_graph_plot.jl

@@ -60,7 +60,7 @@ function gen_distances(g::SRGraphWrap; inc = 0.2)
 end
 
 """
-    plot_network(rn::ReactionSystem; interactive=false)
+    plot_network(rn::ReactionSystem)
 
 Converts a [`ReactionSystem`](@ref) into a GraphMakie plot of the species reaction graph
 (or Petri net representation). Reactions correspond to small green circles, and 
@@ -76,7 +76,6 @@ Notes:
   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))
@@ -119,7 +118,7 @@ function Catalyst.plot_network(rn::ReactionSystem)
 end
 
 """
-    plot_complexes(rn::ReactionSystem; interactive=false)
+    plot_complexes(rn::ReactionSystem)
 
     Creates a GraphMakie plot of the [`ReactionComplex`](@ref)s in `rn`. Reactions
     correspond to arrows and reaction complexes to blue circles.

src/spatial_reaction_systems/lattice_reaction_systems.jl

@@ -8,7 +8,7 @@ const GridLattice{N, T} = Union{Array{Bool, N}, CartesianGridRej{N, T}}
 """
 $(TYPEDEF)
 
-A representation of a spatial system of chemical reactions on a discrete (lattice) space. 
+A representation of a spatial system of chemical reactions on a discrete (lattice) space.
 
 # Fields
 $(FIELDS)
@@ -20,7 +20,7 @@ Arguments:
 to which the non-spatial model is expanded.
 
 Keyword Arguments:
-- `diagonal_connections = false`: Only relevant for Cartesian and masked lattices. If `true`, 
+- `diagonal_connections = false`: Only relevant for Cartesian and masked lattices. If `true`,
 diagonally adjacent compartments are considered adjacent, and spatial reactions in between these
 are possible.
 
@@ -82,12 +82,12 @@ struct LatticeReactionSystem{Q, R, S, T} <: MT.AbstractTimeDependentSystem
     """
     parameters::Vector{Any}
     """
-    Parameters which values are tied to vertices, 
+    Parameters which values are tied to vertices,
     e.g. that possibly could have unique values at each vertex of the system.
     """
     vertex_parameters::Vector{Any}
     """
-    Parameters whose values are tied to edges (adjacencies), 
+    Parameters whose values are tied to edges (adjacencies),
     e.g. that possibly could have unique values at each edge of the system.
     """
     edge_parameters::Vector{Any}
@@ -108,23 +108,23 @@ struct LatticeReactionSystem{Q, R, S, T} <: MT.AbstractTimeDependentSystem
         if !iscomplete(rs)
             throw(ArgumentError("A non-complete `ReactionSystem` was used as input, this is not permitted."))
         end
-        if !isempty(MT.get_systems(rs))
-            throw(ArgumentError("A non-flattened (hierarchical) `ReactionSystem` was used as input. `LatticeReactionSystem`s can only be based on non-hierarchical `ReactionSystem`s."))
-        end
         if length(reactions(rs)) != length(equations(rs))
             throw(ArgumentError("The `ReactionSystem` used as input contain equations (in addition to reactions). This is not permitted."))
         end
         if length(species(rs)) != length(unknowns(rs))
             throw(ArgumentError("The `ReactionSystem` used as input contain variable unknowns (in addition to species unknowns). This is not permitted (the input `ReactionSystem` must contain species unknowns only)."))
         end
         if !isempty(MT.continuous_events(rs)) || !isempty(MT.discrete_events(rs))
-            throw(ArgumentError("The `ReactionSystem` used as input to `LatticeReactionSystem contain events. These will be ignored in any simulations based on the created `LatticeReactionSystem`."))
+            throw(ArgumentError("The `ReactionSystem` used as input to `LatticeReactionSystem` contain events. These will be ignored in any simulations based on the created `LatticeReactionSystem`."))
+        end
+        if !isnothing(MT.get_parent(rs)) && !isempty(MT.get_systems(MT.get_parent(rs)))
+            @warn "The `ReactionSystem` used as input to `LatticeReactionSystem` was originally created as a hierarchical model. While this won't necessarily result in errors, it has not been well-tested, and is not recommended."
         end
         if !isempty(observed(rs))
-            @warn "The `ReactionSystem` used as input to `LatticeReactionSystem contain observables. It will not be possible to access these from the created `LatticeReactionSystem`."
+            @warn "The `ReactionSystem` used as input to `LatticeReactionSystem` contain observables. It will not be possible to access these from the created `LatticeReactionSystem`."
         end
 
-        # Computes the species which are parts of spatial reactions. Also counts the total number of 
+        # Computes the species which are parts of spatial reactions. Also counts the total number of
         # species types.
         if isempty(spatial_reactions)
             spat_species = Vector{BasicSymbolic{Real}}[]
@@ -171,7 +171,7 @@ function LatticeReactionSystem(rs, srs, lattice::SimpleGraph)
     LatticeReactionSystem(rs, srs, DiGraph(lattice))
 end
 
-# Creates a LatticeReactionSystem from a CartesianGrid lattice (cartesian grid) or a Boolean Array 
+# Creates a LatticeReactionSystem from a CartesianGrid lattice (cartesian grid) or a Boolean Array
 # lattice (masked grid). These two are quite similar, so much code can be reused in a single interface.
 function LatticeReactionSystem(rs, srs, lattice::GridLattice{N, T};
         diagonal_connections = false) where {N, T}
@@ -226,7 +226,7 @@ function count_edges(grid::CartesianGridRej{N, T};
     return num_edges
 end
 
-# Counts and edges on a masked grid. Does so by looping through all the vertices of the grid, 
+# Counts and edges on a masked grid. Does so by looping through all the vertices of the grid,
 # finding their neighbours, and updating the edge count accordingly.
 function count_edges(grid::Array{Bool, N}; diagonal_connections = false) where {N}
     g_size = grid_size(grid)
@@ -239,7 +239,7 @@ function count_edges(grid::Array{Bool, N}; diagonal_connections = false) where {
     return num_edges
 end
 
-# For a (1d, 2d, or 3d) (Cartesian or masked) grid, returns a vector and an array, permitting the 
+# For a (1d, 2d, or 3d) (Cartesian or masked) grid, returns a vector and an array, permitting the
 # conversion between a vertex's flat (scalar) and grid indices. E.g. for a 2d grid, if grid point (3,2)
 # corresponds to the fifth vertex, then `flat_to_grid_idx[5] = (3,2)` and `grid_to_flat_idx[3,2] = 5`.
 function get_index_converters(grid::GridLattice{N, T}, num_verts) where {N, T}
@@ -335,7 +335,7 @@ num_verts(lrs::LatticeReactionSystem) = getfield(lrs, :num_verts)
 """
     num_edges(lrs::LatticeReactionSystem)
 
-Returns the number of edges (i.e. connections between vertices) in the lattice stored in a 
+Returns the number of edges (i.e. connections between vertices) in the lattice stored in a
 `LatticeReactionSystem`.
 """
 num_edges(lrs::LatticeReactionSystem) = getfield(lrs, :num_edges)
@@ -392,7 +392,7 @@ end
 """
     has_cartesian_lattice(lrs::LatticeReactionSystem)
 
-Returns `true` if `lrs` was created using a cartesian grid lattice (e.g. created via `CartesianGrid(5,5)`). 
+Returns `true` if `lrs` was created using a cartesian grid lattice (e.g. created via `CartesianGrid(5,5)`).
 Otherwise, returns `false`.
 """
 has_cartesian_lattice(lrs::LatticeReactionSystem) = lattice(lrs) isa
@@ -401,7 +401,7 @@ has_cartesian_lattice(lrs::LatticeReactionSystem) = lattice(lrs) isa
 """
     has_masked_lattice(lrs::LatticeReactionSystem)
 
-Returns `true` if `lrs` was created using a masked grid lattice (e.g. created via `[true true; true false]`). 
+Returns `true` if `lrs` was created using a masked grid lattice (e.g. created via `[true true; true false]`).
 Otherwise, returns `false`.
 """
 has_masked_lattice(lrs::LatticeReactionSystem) = lattice(lrs) isa Array{Bool, N} where {N}
@@ -418,15 +418,15 @@ end
 """
     has_graph_lattice(lrs::LatticeReactionSystem)
 
-Returns `true` if `lrs` was created using a graph grid lattice (e.g. created via `path_graph(5)`). 
+Returns `true` if `lrs` was created using a graph grid lattice (e.g. created via `path_graph(5)`).
 Otherwise, returns `false`.
 """
 has_graph_lattice(lrs::LatticeReactionSystem) = lattice(lrs) isa SimpleDiGraph
 
 """
     grid_size(lrs::LatticeReactionSystem)
 
-Returns the size of `lrs`'s lattice (only if it is a cartesian or masked grid lattice). 
+Returns the size of `lrs`'s lattice (only if it is a cartesian or masked grid lattice).
 E.g. for a lattice `CartesianGrid(4,6)`, `(4,6)` is returned.
 """
 grid_size(lrs::LatticeReactionSystem) = grid_size(lattice(lrs))
@@ -439,7 +439,7 @@ end
 """
     grid_dims(lrs::LatticeReactionSystem)
 
-Returns the number of dimensions of `lrs`'s lattice (only if it is a cartesian or masked grid lattice). 
+Returns the number of dimensions of `lrs`'s lattice (only if it is a cartesian or masked grid lattice).
 The output is either `1`, `2`, or `3`.
 """
 grid_dims(lrs::LatticeReactionSystem) = grid_dims(lattice(lrs))
@@ -506,19 +506,19 @@ end
 """
     make_edge_p_values(lrs::LatticeReactionSystem, make_edge_p_value::Function)
 
-Generates edge parameter values for a lattice reaction system. Only work for (Cartesian or masked) 
-grid lattices (without diagonal adjacencies). 
+Generates edge parameter values for a lattice reaction system. Only work for (Cartesian or masked)
+grid lattices (without diagonal adjacencies).
 
 Input:
 - `lrs`: The lattice reaction system for which values should be generated.
     - `make_edge_p_value`: a function describing a rule for generating the edge parameter values.
 
 Output:
-    - `ep_vals`: A sparse matrix of size (num_verts,num_verts) (where num_verts is the number of 
-    vertices in `lrs`). Here, `eps[i,j]` is filled only if there is an edge going from vertex i to 
+    - `ep_vals`: A sparse matrix of size (num_verts,num_verts) (where num_verts is the number of
+    vertices in `lrs`). Here, `eps[i,j]` is filled only if there is an edge going from vertex i to
     vertex j. The value of `eps[i,j]` is determined by `make_edge_p_value`.
 
-Here, `make_edge_p_value` should take two arguments, `src_vert` and `dst_vert`, which correspond to 
+Here, `make_edge_p_value` should take two arguments, `src_vert` and `dst_vert`, which correspond to
 the grid indices of an edge's source and destination vertices, respectively. It outputs a single value,
 which is the value assigned to that edge.
 
@@ -568,34 +568,34 @@ function make_edge_p_values(lrs::LatticeReactionSystem, make_edge_p_value::Funct
 end
 
 """
-    make_directed_edge_values(lrs::LatticeReactionSystem, x_vals::Tuple{T,T}, y_vals::Tuple{T,T} = (undef,undef), 
+    make_directed_edge_values(lrs::LatticeReactionSystem, x_vals::Tuple{T,T}, y_vals::Tuple{T,T} = (undef,undef),
                      z_vals::Tuple{T,T} = (undef,undef)) where {T}
 
-Generates edge parameter values for a lattice reaction system. Only work for (Cartesian or masked) 
-grid lattices (without diagonal adjacencies). Each dimension (x, and possibly y and z), and 
+Generates edge parameter values for a lattice reaction system. Only work for (Cartesian or masked)
+grid lattices (without diagonal adjacencies). Each dimension (x, and possibly y and z), and
 direction has assigned its own constant edge parameter value.
 
 Input:
     - `lrs`: The lattice reaction system for which values should be generated.
-    - `x_vals::Tuple{T,T}`: The values in the increasing (from a lower x index to a higher x index) 
+    - `x_vals::Tuple{T,T}`: The values in the increasing (from a lower x index to a higher x index)
     and decreasing (from a higher x index to a lower x index) direction along the x dimension.
-    - `y_vals::Tuple{T,T}`: The values in the increasing and decreasing direction along the y dimension. 
+    - `y_vals::Tuple{T,T}`: The values in the increasing and decreasing direction along the y dimension.
     Should only be used for 2 and 3-dimensional grids.
-    - `z_vals::Tuple{T,T}`: The values in the increasing and decreasing direction along the z dimension. 
+    - `z_vals::Tuple{T,T}`: The values in the increasing and decreasing direction along the z dimension.
     Should only be used for 3-dimensional grids.
 
 Output:
-    - `ep_vals`: A sparse matrix of size (num_verts,num_verts) (where num_verts is the number of 
-    vertices in `lrs`). Here, `eps[i,j]` is filled only if there is an edge going from vertex i to 
-    vertex j. The value of `eps[i,j]` is determined by the `x_vals`, `y_vals`, and `z_vals` Tuples, 
+    - `ep_vals`: A sparse matrix of size (num_verts,num_verts) (where num_verts is the number of
+    vertices in `lrs`). Here, `eps[i,j]` is filled only if there is an edge going from vertex i to
+    vertex j. The value of `eps[i,j]` is determined by the `x_vals`, `y_vals`, and `z_vals` Tuples,
     and vertices i and j's relative position in the grid.
 
-It should be noted that two adjacent vertices will always be different in exactly a single dimension 
+It should be noted that two adjacent vertices will always be different in exactly a single dimension
 (x, y, or z). The corresponding tuple determines which value is assigned.
 
 Example:
     In the following example, we wish to have diffusion in the x dimension, but a constant flow from
-    low y values to high y values (so not transportation from high to low y). We achieve it in the 
+    low y values to high y values (so not transportation from high to low y). We achieve it in the
     following manner:
 ```julia
 using Catalyst

test/extensions/Project.toml

@@ -1,12 +1,13 @@
 [deps]
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
-Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
-# StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
+# StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"

test/extensions/graphmakie.jl

@@ -1,4 +1,4 @@
-using Catalyst, GraphMakie, GLMakie, Graphs
+using Catalyst, GraphMakie, CairoMakie, Graphs
 include("../test_networks.jl")
 
 # Test that speciesreactiongraph is generated correctly
@@ -103,4 +103,3 @@ let
     @test count(==(Graphs.Edge(1, s+2)), edges(srg)) == 2
     @test count(==(Graphs.Edge(2, s+3)), edges(srg)) == 2
 end
-

test/extensions/lattice_simulation_plotting.jl

@@ -1,7 +1,8 @@
 ### Preparations ###
 
 # Fetch packages.
-using Catalyst, CairoMakie, GraphMakie, Graphs, JumpProcesses, OrdinaryDiffEqDefault, OrdinaryDiffEq, Test
+using Catalyst, CairoMakie, GraphMakie, Graphs
+using JumpProcesses, OrdinaryDiffEqDefault, OrdinaryDiffEqTsit5, Test
 
 
 ### Checks Basic Plot Cases ###

test/runtests.jl

@@ -86,6 +86,7 @@ end
     if GROUP == "All" || GROUP == "Extensions"
         activate_extensions_env()
 
+        @time @safetestset "Graph visualization" begin include("extensions/graphmakie.jl") end
         @time @safetestset "BifurcationKit Extension" begin include("extensions/bifurcation_kit.jl") end
         @time @safetestset "HomotopyContinuation Extension" begin include("extensions/homotopy_continuation.jl") end
 

test/spatial_modelling/lattice_reaction_systems.jl

@@ -248,7 +248,7 @@ end
     @test_logs (:warn, r"The `ReactionSystem` used as input to `LatticeReactionSystem contain observables. It *") match_mode=:any LatticeReactionSystem(rs4, [tr], short_path)
 end
 
-# Tests for hierarchical input system.
+# Tests for hierarchical input system (should yield a warning).
 let
     t = default_t()
     @parameters d D
@@ -257,7 +257,7 @@ let
     @named rs1 = ReactionSystem(rxs, t)
     @named rs2 = ReactionSystem(rxs, t; systems = [rs1])
     rs2 = complete(rs2)
-    @test_throws ArgumentError LatticeReactionSystem(rs2, [TransportReaction(D, X)], CartesianGrid((2,2)))
+    @test_logs (:warn, r"The `ReactionSystem` used as input to `LatticeReactionSystem` was originally created as a hierarchical model. While *") match_mode=:any LatticeReactionSystem(rs2, [TransportReaction(D, X)], CartesianGrid((2,2)))
 end
 
 # Tests for non-complete input `ReactionSystem`.