From f0fab91e22432760996fd28fe0fd6edad8741cd7 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 6 Feb 2024 19:49:28 -0500 Subject: [PATCH] "finish" remake. --- .../lattice_reaction_systems.jl | 91 ++++---- .../spatial_ODE_systems.jl | 170 ++++++++------- .../spatial_reactions.jl | 27 +-- src/spatial_reaction_systems/utility.jl | 203 ++++++++++-------- ...attice_reaction_systems_ODE_performance.jl | 16 +- .../lattice_reaction_systems_ODEs.jl | 11 +- 6 files changed, 276 insertions(+), 242 deletions(-) diff --git a/src/spatial_reaction_systems/lattice_reaction_systems.jl b/src/spatial_reaction_systems/lattice_reaction_systems.jl index 907fbee47d..a525757753 100644 --- a/src/spatial_reaction_systems/lattice_reaction_systems.jl +++ b/src/spatial_reaction_systems/lattice_reaction_systems.jl @@ -4,57 +4,60 @@ # Adding the "<: MT.AbstractTimeDependentSystem" part messes up show, disabling me from creating LRSs. Should be fixed some time. struct LatticeReactionSystem{Q,R,S,T} # <: MT.AbstractTimeDependentSystem # Input values. - """The reaction system within each compartment.""" + """The (non-spatial) reaction system within each vertexes.""" rs::ReactionSystem{Q} - """The spatial reactions defined between individual nodes.""" + """The spatial reactions defined between individual vertexes.""" spatial_reactions::Vector{R} - """The graph on which the lattice is defined.""" + """The lattice on which the (discrete) spatial system is defined.""" lattice::S # Derived values. - """The number of compartments.""" + """The number of vertexes (compartments).""" num_verts::Int64 """The number of edges.""" num_edges::Int64 """The number of species.""" num_species::Int64 - """Species that may move spatially.""" + """List of species that may move spatially.""" spat_species::Vector{BasicSymbolic{Real}} """ All parameters related to the lattice reaction system - (both with spatial and non-spatial effects). + (both those whose values are tied to vertexes and edges). """ parameters::Vector{BasicSymbolic{Real}} """ - Parameters which values are tied to vertexes (adjacencies), - e.g. (possibly) have a unique value at each vertex of the system. + Parameters which values are tied to vertexes, + e.g. that possibly could have unique values at each vertex of the system. """ vertex_parameters::Vector{BasicSymbolic{Real}} """ - Parameters which values are tied to edges (adjacencies), - e.g. (possibly) have a unique value at each edge of the system. + 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{BasicSymbolic{Real}} """ - An iterator over all the edges on the lattice. - The format depends on the type of lattice (Cartesian grid, grid, or graph). + An iterator over all the lattice's edges. Currently, the format is always a Vector{Pair{Int64,Int64}}. + However, in the future, different types could potentially be used for different types of lattice + (E.g. for a Cartesian grid, we do not technically need to enumerate each edge) """ edge_iterator::T function LatticeReactionSystem(rs::ReactionSystem{Q}, spatial_reactions::Vector{R}, lattice::S, - num_verts, num_edges, edge_iterator::T) where {Q,R, S, T} + num_verts::Int64, num_edges::Int64, edge_iterator::T) where {Q,R,S,T} # Error checks. if !(R <: AbstractSpatialReaction) error("The second argument must be a vector of AbstractSpatialReaction subtypes.") end - # Computes derived values spatial species. Counts the total number of species. + # Computes the species which are parts of spatial reactions. Also counts total number of + # species types. if isempty(spatial_reactions) spat_species = Vector{BasicSymbolic{Real}}[] else spat_species = unique(reduce(vcat, [spatial_species(sr) for sr in spatial_reactions])) end + num_species = length(unique([species(rs); spat_species])) # Computes the sets of vertex, edge, and all, parameters. rs_edge_parameters = filter(isedgeparameter, parameters(rs)) @@ -66,10 +69,6 @@ struct LatticeReactionSystem{Q,R,S,T} # <: MT.AbstractTimeDependentSystem edge_parameters = unique([rs_edge_parameters; srs_edge_parameters]) vertex_parameters = filter(!isedgeparameter, parameters(rs)) - # Counts the number of species types. The number of vertexes and compartments is given in input. - # `num_edges` cannot be computed here, because how it is computed depends on the lattice and `typeof(edge_iterator)`. - num_species = length(unique([species(rs); spat_species])) - # Ensures the parameter order begins similarly to in the non-spatial ReactionSystem. ps = [parameters(rs); setdiff([edge_parameters; vertex_parameters], parameters(rs))] @@ -86,25 +85,26 @@ function LatticeReactionSystem(rs, srs, lattice_in::CartesianGridRej{S,T}; diago # Error checks. (length(lattice_in.dims) > 3) && error("Grids of higher dimension than 3 is currently not supported.") - # Ensures that the matrix has a 3d form (for intermediary computations, original is passed to constructor). + # Ensures that the matrix has a 3d form (used for intermediary computations only, + # the original is passed to the constructor). lattice = CartesianGrid((lattice_in.dims..., fill(1, 3-length(lattice_in.dims))...)) # Counts vertexes and edges. The `num_edges` count formula counts the number of internal, side, - # edge, and corner vertexes (on the grid). The number of edges from each depend on whether diagonal + # edge, and corner vertexes (on the grid). The number of edges from each depends on whether diagonal # connections are allowed. The formula holds even if l, m, and/or n are 1. l,m,n = lattice.dims num_verts = l * m * n (ni, ns, ne, nc) = diagonal_connections ? (26,17,11,7) : (6,5,4,3) - num_edges = ni*(l-2)*(m-2)*(n-2) + # Internal vertexes. - ns*(2(l-2)*(m-2) + 2(l-2)*(n-2) + 2(m-2)*(n-2)) + # Side vertexes. - ne*(4(l-2) + 4(m-2) + 4(n-2)) + # Edge vertexes. - nc*8 # Corner vertexes. + num_edges = ni*(l-2)*(m-2)*(n-2) + # Edges from internal vertexes. + ns*(2(l-2)*(m-2) + 2(l-2)*(n-2) + 2(m-2)*(n-2)) + # Edges from side vertexes. + ne*(4(l-2) + 4(m-2) + 4(n-2)) + # Edges from edge vertexes. + nc*8 # Edges from corner vertexes. # Creates an iterator over all edges. # Currently creates a full vector. Future version might be for efficient. edge_iterator = Vector{Pair{Int64,Int64}}(undef, num_edges) # Loops through, simultaneously, the coordinates of each position in the grid, as well as that - # coordinate's (scalar) flat index. For each grid point, loops through all potential neighbours. + # coordinate's (scalar) flat index. For each grid point, loops through all potential neighbours. indices = [(L, M, N) for L in 1:l, M in 1:m, N in 1:n] flat_indices = 1:num_verts next_vert = 0 @@ -117,7 +117,7 @@ function LatticeReactionSystem(rs, srs, lattice_in::CartesianGridRej{S,T}; diago !diagonal_connections && (count([L==LL, M==MM, N==NN]) == 2) || continue diagonal_connections && (L==LL) && (M==MM) && (N==NN) && continue - # Computes the neighbours scalar index. Add that connection to `edge_iterator`. + # Computes the neighbour's flat (scalar) index. Add the edge to edge_iterator. neighbour_idx = LL + (MM - 1) * l + (NN - 1) * m * l edge_iterator[next_vert += 1] = (idx => neighbour_idx) end @@ -132,17 +132,16 @@ function LatticeReactionSystem(rs, srs, lattice_in::Array{Bool, T}; diagonal_con dims = size(lattice_in) (length(dims) > 3) && error("Grids of higher dimension than 3 is currently not supported.") - # Ensures that the matrix has a 3d form (for intermediary computations, original is passed to constructor). - # Also gets some basic lattice information. + # Ensures that the matrix has a 3d form (used for intermediary computations only, + # the original is passed to the constructor). lattice = reshape(lattice_in, [dims...; fill(1, 3-length(dims))]...) # Counts vertexes (edges have to be counted after the iterator have been created). num_verts = count(lattice) - # Creates an iterator over all edges.Currently a full vector of all edges (as pairs). - edge_iterator = Vector{Pair{Int64,Int64}}() - # Makes a template matrix to store each vertex's index. The matrix is 0 where there are no vertex. + # Makes a template matrix to store each vertex's index. The matrix is 0 where there is no vertex. + # The template is used in the next step. idx_matrix = fill(0, size(lattice_in)) cur_vertex_idx = 0 for flat_idx in 1:length(lattice) @@ -151,8 +150,10 @@ function LatticeReactionSystem(rs, srs, lattice_in::Array{Bool, T}; diagonal_con end end - # Loops through, simultaneously, the coordinates of each position in the grid, as well as that - # coordinate's (scalar) flat index. For each grid point, loops through all potential neighbours. + # Creates an iterator over all edges. A vector with pairs of each edge's source to its destination. + edge_iterator = Vector{Pair{Int64,Int64}}() + # Loops through, the coordinates of each position in the grid. + # For each grid point, loops through all potential neighbours and adds edges to edge_iterator. l, m, n = size(lattice) indices = [(L, M, N) for L in 1:l, M in 1:m, N in 1:n] for (L, M, N) in indices @@ -169,7 +170,7 @@ function LatticeReactionSystem(rs, srs, lattice_in::Array{Bool, T}; diagonal_con !diagonal_connections && (count([L==LL, M==MM, N==NN]) == 2) || continue diagonal_connections && (L==LL) && (M==MM) && (N==NN) && continue - # Computes the neighbours scalar index. Add that connection to `edge_iterator`. + # Computes the neighbour's scalar index. Add that connection to `edge_iterator`. push!(edge_iterator, idx_matrix[L,M,N] => idx_matrix[LL,MM,NN]) end end @@ -188,6 +189,7 @@ end # Creates a LatticeReactionSystem from a (undirected) Graph lattice (graph grid). LatticeReactionSystem(rs, srs, lattice::SimpleGraph) = LatticeReactionSystem(rs, srs, DiGraph(lattice)) + ### Lattice ReactionSystem Getters ### # Get all species. @@ -197,9 +199,9 @@ spatial_species(lrs::LatticeReactionSystem) = lrs.spat_species # Get all parameters. ModelingToolkit.parameters(lrs::LatticeReactionSystem) = lrs.parameters -# Get all parameters which values are tied to vertexes (compartments). +# Get all parameters whose values are tied to vertexes (compartments). vertex_parameters(lrs::LatticeReactionSystem) = lrs.vertex_parameters -# Get all parameters which values are tied to edges (adjacencies). +# Get all parameters whose values are tied to edges (adjacencies). edge_parameters(lrs::LatticeReactionSystem) = lrs.edge_parameters # Gets the lrs name (same as rs name). @@ -211,35 +213,39 @@ is_transport_system(lrs::LatticeReactionSystem) = all(sr -> sr isa TransportReac """ has_cartesian_lattice(lrs::LatticeReactionSystem) -Returns `true` if `lrs` was created using a cartesian grid lattice (e.g. created via `CartesianGrid(5,5)`). Else, returns `false`. +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) = lrs.lattice isa CartesianGridRej{S,T} where {S,T} """ has_masked_lattice(lrs::LatticeReactionSystem) -Returns `true` if `lrs` was created using a masked grid lattice (e.g. created via `[true true; true false]`). Else, returns `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) = lrs.lattice isa Array{Bool, T} where T """ has_grid_lattice(lrs::LatticeReactionSystem) -Returns `true` if `lrs` was created using a cartesian or masked grid lattice. Else, returns `false`. +Returns `true` if `lrs` was created using a cartesian or masked grid lattice. Otherwise, returns `false`. """ has_grid_lattice(lrs::LatticeReactionSystem) = (has_cartesian_lattice(lrs) || has_masked_lattice(lrs)) """ has_graph_lattice(lrs::LatticeReactionSystem) -Returns `true` if `lrs` was created using a graph grid lattice (e.g. created via `path_graph(5)`). Else, returns `false`. +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) = lrs.lattice isa SimpleDiGraph """ grid_size(lrs::LatticeReactionSystem) -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. +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. """ function grid_size(lrs::LatticeReactionSystem) has_cartesian_lattice(lrs) && (return lrs.lattice.dims) @@ -250,6 +256,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). The output is either `1`, `2`, or `3`. +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) = length(grid_size(lrs)) \ No newline at end of file diff --git a/src/spatial_reaction_systems/spatial_ODE_systems.jl b/src/spatial_reaction_systems/spatial_ODE_systems.jl index d073061ca4..28fed3aee4 100644 --- a/src/spatial_reaction_systems/spatial_ODE_systems.jl +++ b/src/spatial_reaction_systems/spatial_ODE_systems.jl @@ -2,43 +2,46 @@ # Functor with information for the forcing function of a spatial ODE with spatial movement on a lattice. struct LatticeTransportODEf{S,T} - """The ODEFunction of the (non-spatial) reaction system which generated this function.""" + """The ODEFunction of the (non-spatial) ReactionSystem that generated this LatticeTransportODEf instance.""" ofunc::S """The number of vertices.""" num_verts::Int64 """The number of species.""" num_species::Int64 - """The values of the parameters which values are tied to vertexes.""" + """The values of the parameters that are tied to vertexes.""" vert_ps::Vector{Vector{T}} """ - Temporary vector. For parameters which values are identical across the lattice, - at some point these have to be converted of a length num_verts vector. - To avoid re-allocation they are written to this vector. + Vector for storing temporary values. Repeatedly during simulations, we need to retrieve the + parameter values in a certain vertex. However, since most parameters (likely) are uniform + (and hence only have 1 value stored), we need to create a new vector each time we need to retrieve + the parameter values in a new vertex. To avoid relocating these values repeatedly, we write them + to this vector. """ work_vert_ps::Vector{T} """ - For each parameter in vert_ps, its value is a vector with length either num_verts or 1. - To know whenever a parameter's value need expanding to the work_vert_ps array, its length needs checking. - This check is done once, and the value stored to this array. True means a uniform value. + For each parameter in vert_ps, its value is a vector with a length of either num_verts or 1. + To know whenever a parameter's value needs expanding to the work_vert_ps array, its length needs checking. + This check is done once, and the value is stored in this array. True means a uniform value. """ v_ps_idx_types::Vector{Bool} """ - A vector of sparse, with a value for each species with transportation. - The first value is the species index (in the species(::ReactionSystem) vector), - and the second is a vector with its transport rate values. - If the transport rate is uniform (across all edges), that value is the only value in the vector. - Else, there is one value for each edge in the lattice. + A vector that stores, for each species with transportation, its transportation rate(s). + Each entry is a pair from (the index of) the transported species (in the `species(lrs)` vector) + to its transportation rate (each species only has a single transportation rate, the sum of all + its transportation reactions' rates). If the transportation rate is uniform across all edges, + stores a single value (in a size (1,1) sparse matrix). Otherwise, stores these in a sparse matrix + where value (i,j) is the species transportation rate from vertex i to vertex j. """ transport_rates::Vector{Pair{Int64,SparseMatrixCSC{T, Int64}}} """ - For each transport rate in transport_rates, its value is a (sparse) matrix with size either - (num_verts,num_verts) or (1,1). In the second case, that transportation rate is uniform across all edges. - To know how to access transport rate's value (without checking sizes), we can use this vector directly. - True means a uniform value. + For each transport rate in transport_rates, its value is a (sparse) matrix with a size of either + (num_verts,num_verts) or (1,1). In the second case, the transportation rate is uniform across + all edges. To avoid having to check which case holds for each transportation rate, we store the + corresponding case in this value. `true` means that a species has a uniform transportation rate. """ t_rate_idx_types::Vector{Bool} """ - A matrix, NxM, where N is the number of species with transportation and M the number of vertexes. + A matrix, NxM, where N is the number of species with transportation and M is the number of vertexes. Each value is the total rate at which that species leaves that vertex (e.g. for a species with constant diffusion rate D, in a vertex with n neighbours, this value is n*D). """ @@ -46,27 +49,28 @@ struct LatticeTransportODEf{S,T} """An iterator over all the edges of the lattice.""" edge_iterator::Vector{Pair{Int64, Int64}} - function LatticeTransportODEf(ofunc::S, vert_ps::Vector{Pair{BasicSymbolic{Real},Vector{T}}}, transport_rates::Vector{Pair{Int64, SparseMatrixCSC{T, Int64}}}, + function LatticeTransportODEf(ofunc::S, vert_ps::Vector{Pair{BasicSymbolic{Real},Vector{T}}}, + transport_rates::Vector{Pair{Int64, SparseMatrixCSC{T, Int64}}}, lrs::LatticeReactionSystem) where {S,T} - # Records. which parameters and rates are uniform and not. + # Records which parameters and rates are uniform and which are not. v_ps_idx_types = map(vp -> length(vp[2]) == 1, vert_ps) t_rate_idx_types = map(tr -> size(tr[2]) == (1,1), transport_rates) - # Input `vert_ps` is a vector map taking each parameter symbolic to its value (potentially a vector). - # This vector is sorted according to the parameters order. Here, we simply extract its values only. + # Input `vert_ps` is a vector map taking each parameter symbolic to its value (potentially a + # vector). This vector is already sorted according to the order of the parameters. Here, we extract + # its values only and put them into `vert_ps`. vert_ps = [vp[2] for vp in vert_ps] - # Computes the leaving rates. + # Computes the leaving rate matrix. leaving_rates = zeros(length(transport_rates), lrs.num_verts) - for (s_idx, trpair) in enumerate(transport_rates) - t_rate = trpair[2] + for (s_idx, tr_pair) in enumerate(transport_rates) for e in lrs.edge_iterator - # Updates the exit rate for species s_idx from vertex e.src - leaving_rates[s_idx, e[1]] += get_transport_rate(t_rate, e, t_rate_idx_types[s_idx]) + # Updates the exit rate for species s_idx from vertex e.src. + leaving_rates[s_idx, e[1]] += get_transport_rate(tr_pair[2], e, t_rate_idx_types[s_idx]) end end - # Declares `work_vert_ps` (used as storage during computation) and an iterator over the edges. + # Declares `work_vert_ps` (used as storage during computation) and the edge iterator. work_vert_ps = zeros(length(vert_ps)) edge_iterator = lrs.edge_iterator new{S,T}(ofunc, lrs.num_verts, lrs.num_species, vert_ps, work_vert_ps, @@ -76,25 +80,26 @@ end # Functor with information for the Jacobian function of a spatial ODE with spatial movement on a lattice. struct LatticeTransportODEjac{R,S,T} - """The ODEFunction of the (non-spatial) reaction system which generated this function.""" + """The ODEFunction of the (non-spatial) ReactionSystem that generated this LatticeTransportODEf instance.""" ofunc::R """The number of vertices.""" num_verts::Int64 """The number of species.""" num_species::Int64 - """The values of the parameters which values are tied to vertexes.""" + """The values of the parameters that are tied to vertexes.""" vert_ps::Vector{Vector{S}} """ - Temporary vector. For parameters which values are identical across the lattice, - at some point these have to be converted of a length(num_verts) vector. - To avoid re-allocation they are written to this vector. + Vector for storing temporary values. Repeatedly during simulations, we need to retrieve the + parameter values in a certain vertex. However, since most parameters (likely) are uniform + (and hence only have 1 value stored), we need to create a new vector each time we need to retrieve + the parameter values in a new vertex. To avoid relocating these values repeatedly, we write them + to this vector. """ - work_vert_ps::Vector{S} + work_vert_ps::Vector{S} """ - For each parameter in vert_ps, it either have length num_verts or 1. - To know whenever a parameter's value need expanding to the work_vert_ps array, - its length needs checking. This check is done once, and the value stored to this array. - This field (specifically) is an enumerate over that array. + For each parameter in vert_ps, its value is a vector with a length of either num_verts or 1. + To know whenever a parameter's value needs expanding to the work_vert_ps array, its length needs checking. + This check is done once, and the value is stored in this array. True means a uniform value. """ v_ps_idx_types::Vector{Bool} """Whether the Jacobian is sparse or not.""" @@ -102,11 +107,12 @@ struct LatticeTransportODEjac{R,S,T} """The transport rates. This is a dense or sparse matrix (depending on what type of Jacobian is used).""" jac_transport::T - function LatticeTransportODEjac(ofunc::R, vert_ps::Vector{Pair{BasicSymbolic{Real},Vector{S}}}, lrs::LatticeReactionSystem, + function LatticeTransportODEjac(ofunc::R, vert_ps::Vector{Pair{BasicSymbolic{Real},Vector{S}}}, jac_transport::Union{Nothing, SparseMatrixCSC{Float64, Int64}}, - sparse::Bool) where {R,S} - # Input `vert_ps` is a vector map taking each parameter symbolic to its value (potentially a vector). - # This vector is sorted according to the parameters order. Here, we simply extract its values only. + lrs::LatticeReactionSystem, sparse::Bool) where {R,S} + # Input `vert_ps` is a vector map taking each parameter symbolic to its value (potentially a + # vector). This vector is already sorted according to the order of the parameters. Here, we extract + # its values only and put them into `vert_ps`. vert_ps = [vp[2] for vp in vert_ps] work_vert_ps = zeros(lrs.num_verts) @@ -129,26 +135,24 @@ function DiffEqBase.ODEProblem(lrs::LatticeReactionSystem, u0_in, tspan, error("Currently lattice ODE simulations are only supported when all spatial reactions are TransportReactions.") end - # Converts potential symmaps to varmaps - # Vertex and edge parameters may be given in a tuple, or in a common vector, making parameter case complicated. + # Converts potential symmaps to varmaps. u0_in = symmap_to_varmap(lrs, u0_in) p_in = symmap_to_varmap(lrs, p_in) # Converts u0 and p to their internal forms. - # u0 is simply a vector with all the species initial condition values across all vertexes. + # u0 is simply a vector with all the species' initial condition values across all vertexes. # u0 is [spec 1 at vert 1, spec 2 at vert 1, ..., spec 1 at vert 2, ...]. u0 = lattice_process_u0(u0_in, species(lrs), lrs) - # vert_ps and `edge_ps` are vector maps, taking each parameter's Symbolics to its value(s). - # vert_ps values are vectors. Here, index (i) is a parameters value in vertex i. - # edge_ps becomes sparse matrix. Here, index (i,j) is a parameters value in the edge from vertex i to vertex j. - # Uniform vertex/edge parameters stores only a single value (in a length 1 vector, or size 1x1 sparse matrix). - # This is the parameters single value. + # vert_ps and `edge_ps` are vector maps, taking each parameter's Symbolics representation to its value(s). + # vert_ps values are vectors. Here, index (i) is a parameter's value in vertex i. + # edge_ps becomes a sparse matrix. Here, index (i,j) is a parameter's value in the edge from vertex i to vertex j. + # Uniform vertex/edge parameters store only a single value (a length 1 vector, or size 1x1 sparse matrix). # In the `ODEProblem` vert_ps and edge_ps are merged (but for building the ODEFunction, they are separate). vert_ps, edge_ps = lattice_process_p(p_in, vertex_parameters(lrs), edge_parameters(lrs), lrs) - # Creates ODEProblem. + # Creates the ODEFunction. ofun = build_odefunction(lrs, vert_ps, edge_ps, jac, sparse, name, include_zero_odes, - combinatoric_ratelaws, remove_conserved, checks) + combinatoric_ratelaws, remove_conserved, checks) # Combines `vert_ps` and `edge_ps` to a single vector with values only (not a map). Creates ODEProblem. ps = [p[2] for p in [vert_ps; edge_ps]] @@ -156,34 +160,35 @@ function DiffEqBase.ODEProblem(lrs::LatticeReactionSystem, u0_in, tspan, end # Builds an ODEFunction for a spatial ODEProblem. -function build_odefunction(lrs::LatticeReactionSystem, vert_ps::Vector{Pair{A,Vector{T}}}, - edge_ps::Vector{Pair{B,SparseMatrixCSC{T, Int64}}}, jac::Bool, sparse::Bool, - name, include_zero_odes, combinatoric_ratelaws, remove_conserved, checks) where {A,B,T} +function build_odefunction(lrs::LatticeReactionSystem, vert_ps::Vector{Pair{BasicSymbolic{Real},Vector{T}}}, + edge_ps::Vector{Pair{BasicSymbolic{Real},SparseMatrixCSC{T, Int64}}}, + jac::Bool, sparse::Bool, name, include_zero_odes, combinatoric_ratelaws, + remove_conserved, checks) where {T} if remove_conserved error("Removal of conserved quantities is currently not supported for `LatticeReactionSystem`s") end - # Creates a map, taking (the index in species(lrs) each species (with transportation) - # to its transportation rate (uniform or one value for each edge). + # Creates a map, taking (the index in species(lrs) of) each species (with transportation) + # to its transportation rate (uniform or one value for each edge). The rates are sparse matrices. transport_rates = make_sidxs_to_transrate_map(vert_ps, edge_ps, lrs) - # Prepares the Jacobian and forcing functions (depending on jacobian and sparsity selection). + # Prepares the Jacobian and forcing functions (depending on Jacobian and sparsity selection). osys = convert(ODESystem, lrs.rs; name, combinatoric_ratelaws, include_zero_odes, checks) if jac # `build_jac_prototype` currently assumes a sparse (non-spatial) Jacobian. Hence compute this. # `LatticeTransportODEjac` currently assumes a dense (non-spatial) Jacobian. Hence compute this. - # Long term we could write separate version of these functions for generic input. + # Long term we could write separate versions of these functions for generic input. ofunc_dense = ODEFunction(osys; jac = true, sparse = false) ofunc_sparse = ODEFunction(osys; jac = true, sparse = true) - jac_vals = build_jac_prototype(ofunc_sparse.jac_prototype, transport_rates, lrs; set_nonzero = true) + jac_transport = build_jac_prototype(ofunc_sparse.jac_prototype, transport_rates, lrs; set_nonzero = true) if sparse f = LatticeTransportODEf(ofunc_sparse, vert_ps, transport_rates, lrs) - jac_vals = build_jac_prototype(ofunc_sparse.jac_prototype, transport_rates, lrs; set_nonzero = true) - J = LatticeTransportODEjac(ofunc_dense, vert_ps, lrs, jac_vals, true) - jac_prototype = jac_vals + jac_transport = build_jac_prototype(ofunc_sparse.jac_prototype, transport_rates, lrs; set_nonzero = true) + J = LatticeTransportODEjac(ofunc_dense, vert_ps, jac_transport, lrs, true) + jac_prototype = jac_transport else f = LatticeTransportODEf(ofunc_dense, vert_ps, transport_rates, lrs) - J = LatticeTransportODEjac(ofunc_dense, vert_ps, lrs, jac_vals, false) + J = LatticeTransportODEjac(ofunc_dense, vert_ps, jac_transport, lrs, false) jac_prototype = nothing end else @@ -202,10 +207,13 @@ function build_odefunction(lrs::LatticeReactionSystem, vert_ps::Vector{Pair{A,Ve return ODEFunction(f; jac = J, jac_prototype = jac_prototype) end -# Builds a jacobian prototype. If requested, populate it with the Jacobian's (constant) values as well. -function build_jac_prototype(ns_jac_prototype::SparseMatrixCSC{Float64, Int64}, transport_rates::Vector{Pair{Int64,SparseMatrixCSC{T, Int64}}}, - lrs::LatticeReactionSystem; set_nonzero = false) where T - # Finds the indexes of the transport species, and the species with transport only (and no non-spatial dynamics). +# Builds a jacobian prototype. +# If requested, populate it with the constant values of the Jacobian's transportation part. +function build_jac_prototype(ns_jac_prototype::SparseMatrixCSC{Float64, Int64}, + transport_rates::Vector{Pair{Int64,SparseMatrixCSC{T, Int64}}}, + lrs::LatticeReactionSystem; set_nonzero = false) where {T} + # Finds the indexes of both the transport species, + # and the species with transport only (that is, with no non-spatial dynamics but with spatial dynamics). trans_species = [tr[1] for tr in transport_rates] trans_only_species = filter(s_idx -> !Base.isstored(ns_jac_prototype, s_idx, s_idx), trans_species) @@ -221,7 +229,7 @@ function build_jac_prototype(ns_jac_prototype::SparseMatrixCSC{Float64, Int64}, i_idxs = Vector{Int}(undef, num_entries) j_idxs = Vector{Int}(undef, num_entries) - # Indexes of elements due to non-spatial dynamics. + # Indexes of elements caused by non-spatial dynamics. for vert in 1:lrs.num_verts for n in 1:length(ns_i_idxs) i_idxs[idx] = get_index(vert, ns_i_idxs[n], lrs.num_species) @@ -230,9 +238,9 @@ function build_jac_prototype(ns_jac_prototype::SparseMatrixCSC{Float64, Int64}, end end - # Indexes of elements due to spatial dynamics. + # Indexes of elements caused by spatial dynamics. for e in lrs.edge_iterator - # Indexes due to terms for a species leaves its current vertex (but does not have + # Indexes due to terms for a species leaving its source vertex (but does not have # non-spatial dynamics). If the non-spatial Jacobian is fully dense, these would already # be accounted for. for s_idx in trans_only_species @@ -240,7 +248,7 @@ function build_jac_prototype(ns_jac_prototype::SparseMatrixCSC{Float64, Int64}, j_idxs[idx] = i_idxs[idx] idx += 1 end - # Indexes due to terms for species arriving into a new vertex. + # Indexes due to terms for species arriving into a destination vertex. for s_idx in trans_species i_idxs[idx] = get_index(e[1], s_idx, lrs.num_species) j_idxs[idx] = get_index(e[2], s_idx, lrs.num_species) @@ -248,7 +256,7 @@ function build_jac_prototype(ns_jac_prototype::SparseMatrixCSC{Float64, Int64}, end end - # Create sparse jacobian prototype with 0-valued entries. + # Create a sparse Jacobian prototype with 0-valued entries. jac_prototype = sparse(i_idxs, j_idxs, zeros(num_entries)) # Set element values. @@ -273,25 +281,25 @@ end function (f_func::LatticeTransportODEf)(du, u, p, t) # Updates for non-spatial reactions. for vert_i in 1:(f_func.num_verts) - # Gets the indices of species at vertex i. + # Gets the indexes of all the species at vertex i. idxs = get_indexes(vert_i, f_func.num_species) - # Updates the vector which contains the vertex parameter values for vertex vert_i. + # Updates the work vector to contain the vertex parameter values for vertex vert_i. update_work_vert_ps!(f_func, p, vert_i) # Evaluate reaction contributions to du at vert_i. f_func.ofunc((@view du[idxs]), (@view u[idxs]), f_func.work_vert_ps, t) end - # s_idx is species index among transport species, s is index among all species. + # s_idx is the species index among transport species, s is the index among all species. # rates are the species' transport rates. for (s_idx, (s, rates)) in enumerate(f_func.transport_rates) - # Rate for leaving vert_i + # Rate for leaving source vertex vert_i. for vert_i in 1:(f_func.num_verts) idx_src = get_index(vert_i, s, f_func.num_species) du[idx_src] -= f_func.leaving_rates[s_idx, vert_i] * u[idx_src] end - # Add rates for entering a given vertex via an incoming edge. + # Add rates for entering a destination vertex via an incoming edge. for e in f_func.edge_iterator idx_src = get_index(e[1], s, f_func.num_species) idx_dst = get_index(e[2], s, f_func.num_species) @@ -300,17 +308,17 @@ function (f_func::LatticeTransportODEf)(du, u, p, t) end end -# Defines the jacobian functor's effect on the (spatial) ODE system. +# Defines the Jacobian functor's effect on the (spatial) ODE system. function (jac_func::LatticeTransportODEjac)(J, u, p, t) J .= 0.0 - # Update the Jacobian from reaction terms. + # Update the Jacobian from non-spatial reaction terms. for vert_i in 1:(jac_func.num_verts) idxs = get_indexes(vert_i, jac_func.num_species) update_work_vert_ps!(jac_func, p, vert_i) jac_func.ofunc.jac((@view J[idxs, idxs]), (@view u[idxs]), jac_func.work_vert_ps, t) end - # Updates for the spatial reactions (adds the Jacobian values from the diffusion reactions). + # Updates for the spatial reactions (adds the Jacobian values from the transportation reactions). J .+= jac_func.jac_transport end diff --git a/src/spatial_reaction_systems/spatial_reactions.jl b/src/spatial_reaction_systems/spatial_reactions.jl index c1a01a3888..a70a788784 100644 --- a/src/spatial_reaction_systems/spatial_reactions.jl +++ b/src/spatial_reaction_systems/spatial_reactions.jl @@ -3,7 +3,7 @@ # Abstract spatial reaction structures. abstract type AbstractSpatialReaction end -### EdgeParameter Metadata ### +### Edge Parameter Metadata ### # Implements the edgeparameter metadata field. struct EdgeParameter end @@ -17,12 +17,13 @@ function isedgeparameter(x, default = false) Symbolics.getmetadata(x, EdgeParameter, default) end + ### Transport Reaction Structures ### # A transport reaction. These are simple to handle, and should cover most types of spatial reactions. # Only permit constant rates (possibly consisting of several parameters). struct TransportReaction <: AbstractSpatialReaction - """The rate function (excluding mass action terms). Currently only constants supported""" + """The rate function (excluding mass action terms). Currently, only constants supported""" rate::Any """The species that is subject to diffusion.""" species::BasicSymbolic{Real} @@ -40,7 +41,7 @@ function TransportReactions(transport_reactions) [TransportReaction(tr[1], tr[2]) for tr in transport_reactions] end -# Macro for creating a transport reaction. +# Macro for creating a TransportReactions. macro transport_reaction(rateex::ExprValues, species::ExprValues) make_transport_reaction(MacroTools.striplines(rateex), species) end @@ -70,17 +71,17 @@ function make_transport_reaction(rateex, species) end end -# Gets the parameters in a transport reaction. +# Gets the parameters in a TransportReactions. ModelingToolkit.parameters(tr::TransportReaction) = Symbolics.get_variables(tr.rate) -# Gets the species in a transport reaction. +# Gets the species in a TransportReactions. spatial_species(tr::TransportReaction) = [tr.species] -# Checks that a transport reaction is valid for a given reaction system. +# Checks that a TransportReactions is valid for a given reaction system. function check_spatial_reaction_validity(rs::ReactionSystem, tr::TransportReaction; edge_parameters=[]) # Checks that the species exist in the reaction system. # (ODE simulation code becomes difficult if this is not required, - # as non-spatial jacobian and f function generated from rs is of wrong size). + # as non-spatial jacobian and f function generated from rs are of the wrong size). if !any(isequal(tr.species), species(rs)) error("Currently, species used in TransportReactions must have previously been declared within the non-spatial ReactionSystem. This is not the case for $(tr.species).") end @@ -93,16 +94,16 @@ function check_spatial_reaction_validity(rs::ReactionSystem, tr::TransportReacti # Checks that the species does not exist in the system with different metadata. if any([isequal(tr.species, s) && !isequivalent(tr.species, s) for s in species(rs)]) - error("A transport reaction used a species, $(tr.species), with metadata not matching its lattice reaction system. Please fetch this species from the reaction system and used in transport reaction creation.") + error("A transport reaction used a species, $(tr.species), with metadata not matching its lattice reaction system. Please fetch this species from the reaction system and use it in the TransportReaction's creation.") end if any([isequal(rs_p, tr_p) && !equivalent_metadata(rs_p, tr_p) for rs_p in parameters(rs), tr_p in Symbolics.get_variables(tr.rate)]) error("A transport reaction used a parameter with metadata not matching its lattice reaction system. Please fetch this parameter from the reaction system and used in transport reaction creation.") end - # Checks that no edge parameter occur among rates of non-spatial reactions. + # Checks that no edge parameter occurs among rates of non-spatial reactions. if any([!isempty(intersect(Symbolics.get_variables(r.rate), edge_parameters)) for r in reactions(rs)]) - error("Edge paramter(s) were found as a rate of a non-spatial reaction.") + error("Edge parameter(s) were found as a rate of a non-spatial reaction.") end end equivalent_metadata(p1, p2) = isempty(setdiff(p1.metadata, p2.metadata, [Catalyst.EdgeParameter => true])) @@ -133,8 +134,10 @@ function hash(tr::TransportReaction, h::UInt) Base.hash(tr.species, h) end + ### Utility ### -# Loops through a rate and extract all parameters. + +# Loops through a rate and extracts all parameters. function find_parameters_in_rate!(parameters, rateex::ExprValues) if rateex isa Symbol if rateex in [:t, :∅, :im, :nothing, CONSERVED_CONSTANT_SYMBOL] @@ -143,7 +146,7 @@ function find_parameters_in_rate!(parameters, rateex::ExprValues) push!(parameters, rateex) end elseif rateex isa Expr - # Note, this (correctly) skips $(...) expressions + # Note, this (correctly) skips $(...) expressions. for i in 2:length(rateex.args) find_parameters_in_rate!(parameters, rateex.args[i]) end diff --git a/src/spatial_reaction_systems/utility.jl b/src/spatial_reaction_systems/utility.jl index eeab03a78c..0b0fa55d63 100644 --- a/src/spatial_reaction_systems/utility.jl +++ b/src/spatial_reaction_systems/utility.jl @@ -13,45 +13,46 @@ function _symbol_to_var(lrs::LatticeReactionSystem, sym) error("Could not find property parameter/species $sym in lattice reaction system.") end -# From u0 input, extracts their values and store them in the internal format. -# Internal format: a vector on the form [spec 1 at vert 1, spec 2 at vert 1, ..., spec 1 at vert 2, ...]). -function lattice_process_u0(u0_in, u0_syms, lrs::LatticeReactionSystem) +# From u0 input, extract their values and store them in the internal format. +# Internal format: a vector on the form [spec 1 at vert 1, spec 2 at vert 1, ..., spec 1 at vert 2, ...]). +function lattice_process_u0(u0_in, u0_syms::Vector{BasicSymbolic{Real}}, lrs::LatticeReactionSystem) # u0 values can be given in various forms. This converts it to a Vector{Pair{Symbolics,...}} form. # Top-level vector: Maps each species to its value(s). u0 = lattice_process_input(u0_in, u0_syms) # Species' initial condition values can be given in different forms (also depending on the lattice). - # This converts each species's values to a Vector. For species with uniform initial conditions, - # The value vector holds that value only. For spatially heterogeneous initial conditions, - # the vector have teh same length as the number of vertexes (with one value for each). + # This converts each species's values to a Vector. In it, for species with uniform initial conditions, + # it holds that value only. For spatially heterogeneous initial conditions, + # the vector has the same length as the number of vertexes (storing one value for each). u0 = vertex_value_map(u0, lrs) - # Converts the initial condition to a single Vector (with one values for each species and vertex). + # Converts the initial condition to a single Vector (with one value for each species and vertex). return expand_component_values([entry[2] for entry in u0], lrs.num_verts) end -# From p input, splits it into diffusion parameters and compartment parameters. +# From a parameter input, split it into vertex parameters and edge parameters. # Store these in the desired internal format. -function lattice_process_p(ps_in, ps_vertex_syms, ps_edge_syms, lrs::LatticeReactionSystem) +function lattice_process_p(ps_in, ps_vertex_syms::Vector{BasicSymbolic{Real}}, + ps_edge_syms::Vector{BasicSymbolic{Real}}, lrs::LatticeReactionSystem) # p values can be given in various forms. This converts it to a Vector{Pair{Symbolics,...}} form. # Top-level vector: Maps each parameter to its value(s). # Second-level: Contains either a vector (vertex parameters) or a sparse matrix (edge parameters). - # For uniform parameters these have size 1/1x1. Else, they have size num_verts/num_vertsxnum_verts. + # For uniform parameters these have size 1/(1,1). Else, they have size num_verts/(num_verts,num_verts). ps = lattice_process_input(ps_in, [ps_vertex_syms; ps_edge_syms]) # Split the parameter vector into one for vertex parameters and one for edge parameters. - # Next, converts the values to the correct form (vectors for vert_ps and sparse matrices for edge_ps). + # Next, convert their values to the correct form (vectors for vert_ps and sparse matrices for edge_ps). vert_ps, edge_ps = split_parameters(ps, ps_vertex_syms, ps_edge_syms) vert_ps = vertex_value_map(vert_ps, lrs) edge_ps = edge_value_map(edge_ps, lrs) - + return vert_ps, edge_ps end # The input (parameters or initial conditions) may either be a dictionary (symbolics to value(s).) # or a map (in vector or tuple form) from symbolics to value(s). This converts the input to a # (Vector) map from symbolics to value(s), where the entries have the same order as `syms`. -function lattice_process_input(input::Dict{BasicSymbolic{Real}, <:Any}, syms::Vector{BasicSymbolic{Real}}) +function lattice_process_input(input::Dict{BasicSymbolic{Real}, T}, syms::Vector{BasicSymbolic{Real}}) where {T} # Error checks if !isempty(setdiff(keys(input), syms)) error("You have provided values for the following unrecognised parameters/initial conditions: $(setdiff(keys(input), syms)).") @@ -62,40 +63,56 @@ function lattice_process_input(input::Dict{BasicSymbolic{Real}, <:Any}, syms::Ve return [sym => input[sym] for sym in syms] end -function lattice_process_input(input, syms::Vector{BasicSymbolic{Real}}) - lattice_process_input(Dict(input), syms) +function lattice_process_input(input, syms::Vector{BasicSymbolic{Real}}) + if ((input isa Vector) || (input isa Vector)) && all(entry isa Pair for entry in input) + return lattice_process_input(Dict(input), syms) + end + error("Input parameters/initial conditions have the wrong format ($(typeof(input))). These should either be a Dictionary, or a Tuple or a Vector (where each entry is a Pair taking a parameter/species to its value).") end # Splits parameters into vertex and edge parameters. -#function split_parameters(ps::Vector{<: Pair}, p_vertex_syms::Vector, p_edge_syms::Vector) -function split_parameters(ps, p_vertex_syms, p_edge_syms) +# function split_parameters(ps::Vector{<: Pair}, p_vertex_syms::Vector, p_edge_syms::Vector) +function split_parameters(ps, p_vertex_syms::Vector{BasicSymbolic{Real}}, p_edge_syms::Vector{BasicSymbolic{Real}}) vert_ps = [p for p in ps if any(isequal(p[1]), p_vertex_syms)] edge_ps = [p for p in ps if any(isequal(p[1]), p_edge_syms)] return vert_ps, edge_ps end -# Converts the values for initial condition/vertex parameters to the correct form: -# Map from symbolics to to vectors of length either 1 (for uniform values) or num_verts. -function vertex_value_map(values, lrs) +# Converts the values for the initial conditions/vertex parameters to the correct form: +# A map vector from symbolics to vectors of either length 1 (for uniform values) or num_verts. +function vertex_value_map(values, lrs::LatticeReactionSystem) isempty(values) && (return Pair{BasicSymbolic{Real}, Vector{Float64}}[]) return [entry[1] => vertex_value_form(entry[2], lrs, entry[1]) for entry in values] end -# Converts the values for a specific component (species/parameter) to the correct vector form. -function vertex_value_form(values, lrs::LatticeReactionSystem, sym) + +# Converts the values for an individual species/vertex parameter to its correct vector form. +function vertex_value_form(values, lrs::LatticeReactionSystem, sym::BasicSymbolic{Real}) + # If the value is a scalar (i.e. uniform across the lattice), return it in vector form. (values isa AbstractArray) || (return [values]) + + # If the value is a vector (something all three lattice types accept). if values isa Vector + # For the case where we have a 1d (Cartesian or masked) grid, and the vector's values + # correspond to individual grid points. if has_grid_lattice(lrs) && (size(values) == grid_size(lrs)) vertex_value_form(values, lrs.num_verts, lrs.lattice, sym) end + + # For the case where the i'th value of the vector corresponds to the value in the i'th vertex. + # This is the only (non-uniform) case possible for graph grids. if (length(values) != lrs.num_verts) error("You have provided ($(length(values))) values for $sym. This is not equal to the number of vertexes ($(lrs.num_verts)).") end return values end + + # (2d and 3d) Cartesian and masked grids can take non-vector, non-scalar, values input. return vertex_value_form(values, lrs.num_verts, lrs.lattice, sym) end -# Converts values to correct vector form for a Cartesian grid lattice. -function vertex_value_form(values::AbstractArray, num_verts::Int64, lattice::CartesianGridRej{S,T}, sym) where {S,T} + +# Converts values to the correct vector form for a Cartesian grid lattice. +function vertex_value_form(values::AbstractArray, num_verts::Int64, lattice::CartesianGridRej{S,T}, + sym::BasicSymbolic{Real}) where {S,T} if size(values) != lattice.dims error("The values for $sym did not have the same format as the lattice. Expected a $(lattice.dims) array, got one of size $(size(values))") end @@ -104,32 +121,40 @@ function vertex_value_form(values::AbstractArray, num_verts::Int64, lattice::Car end return [values[flat_idx] for flat_idx in 1:num_verts] end -# Converts values to correct vector form for a masked grid lattice. -function vertex_value_form(values::AbstractArray, num_verts::Int64, lattice::Array{Bool, T}, sym) where {T} + +# Converts values to the correct vector form for a masked grid lattice. +function vertex_value_form(values::AbstractArray, num_verts::Int64, lattice::Array{Bool,T}, + sym::BasicSymbolic{Real}) where {T} if size(values) != size(lattice) error("The values for $sym did not have the same format as the lattice. Expected a $(size(lattice)) array, got one of size $(size(values))") end + + # Pre-declares a vector with the values in each vertex (return_values). + # Loops through the lattice and the values, adding these to the return_values. return_values = Vector{typeof(values[1])}(undef, num_verts) cur_idx = 0 - for i = 1:length(lattice) - lattice[i] || continue - return_values[cur_idx += 1] = values[i] + for (idx,val) = enumerate(values) + lattice[idx] || continue + return_values[cur_idx += 1] = val end + + # Checks that the correct number of values was provided, and returns the values. if (length(return_values) != num_verts) error("You have provided ($(length(return_values))) values for $sym. This is not equal to the number of vertexes ($(num_verts)).") end return return_values end -# Converts the values for initial condition/vertex parameters to the correct form: -# Map from symbolics to to vectors of length either 1 (for uniform values) or num_verts. -function edge_value_map(values, lrs) +# Converts the values for the edge parameters to the correct form: +# A map vector from symbolics to sparse matrices of size either (1,1) or (num_verts,num_verts). +function edge_value_map(values, lrs::LatticeReactionSystem) isempty(values) && (return Pair{BasicSymbolic{Real}, SparseMatrixCSC{Float64, Int64}}[]) return [entry[1] => edge_value_form(entry[2], lrs, entry[1]) for entry in values] end -# Converts the values for a specific component (species/parameter) to the correct vector form. + +# Converts the values for an individual edge parameter to its correct sparse matrix form. function edge_value_form(values, lrs::LatticeReactionSystem, sym) - # If a scalar have been given, converts it to a size (1,1) sparse matrix. + # If the value is a scalar (i.e. uniform across the lattice), return it in sparse matrix form. (values isa SparseMatrixCSC) || (return sparse([1], [1], [values])) # Error checks. @@ -140,6 +165,8 @@ function edge_value_form(values, lrs::LatticeReactionSystem, sym) error("Values was not provided for some edges for edge parameter $sym.") end + # Unlike initial conditions/vertex parameters, (unless uniform) edge parameters' values are + # always provided in the same (sparse matrix) form. return values end @@ -147,12 +174,14 @@ end # The species is represented by its index (in species(lrs). # If the rate is uniform across all edges, the transportation rate will be a size (1,1) sparse matrix. # Else, the rate will be a size (num_verts,num_verts) sparse matrix. -# In the first step, computes a map from species symbolics form to value(s). -# Second step converts to map from species index to value(s). function make_sidxs_to_transrate_map(vert_ps::Vector{Pair{BasicSymbolic{Real},Vector{T}}}, edge_ps::Vector{Pair{BasicSymbolic{Real},SparseMatrixCSC{T, Int64}}}, - lrs::LatticeReactionSystem) where T + lrs::LatticeReactionSystem) where {T} + # Creates a dictionary with each parameter's value(s). p_val_dict = Dict(vcat(vert_ps, edge_ps)) + + # First, compute a map from species in their symbolics form to their values. + # Next, convert to map from species index to values. transport_rates_speciesmap = compute_all_transport_rates(p_val_dict, lrs) return Pair{Int64,SparseMatrixCSC{T, Int64}}[ speciesmap(lrs.rs)[spat_rates[1]] => spat_rates[2] for spat_rates in transport_rates_speciesmap @@ -160,38 +189,30 @@ function make_sidxs_to_transrate_map(vert_ps::Vector{Pair{BasicSymbolic{Real},Ve end # Computes the transport rates for all species with transportation rates. Output is a map -# taking each species; symbolics form to its transportation rates across all edges. +# taking each species' symbolics form to its transportation rates across all edges. function compute_all_transport_rates(p_val_dict, lrs::LatticeReactionSystem) # For all species with transportation, compute their transportation rate (across all edges). # This is a vector, pairing each species to these rates. - unsorted_rates = [s => compute_transport_rates(get_transport_rate_law(s, lrs), p_val_dict, lrs) - for s in spatial_species(lrs)] + unsorted_rates = [s => compute_transport_rates(s, p_val_dict, lrs) for s in spatial_species(lrs)] - # Sorts all the species => rate pairs according to their species index in species(::ReactionSystem). + # Sorts all the species => rate pairs according to their species index in species(lrs). return sort(unsorted_rates; by = rate -> findfirst(isequal(rate[1]), species(lrs))) end -# For a species, retrieves the symbolic expression for its transportation rate -# (likely only a single parameter, such as `D`, but could be e.g. L*D, where L and D are parameters). -# If there are several transportation reactions for the species, their sum is used. -#function get_transport_rate_law(s::BasicSymbolic{Real}, lrs::LatticeReactionSystem) -function get_transport_rate_law(s, lrs) - rates = filter(sr -> isequal(s, sr.species), lrs.spatial_reactions) - return sum(getfield.(rates, :rate)) -end -# For the numeric expression describing the rate of transport (likely only a single parameter, e.g. `D`), -# and the values of all our parameters, computes the transport rate(s). -# If all parameters the rate depend on are uniform all edges, this becomes a length 1 vector. -# Else a vector with each value corresponding to the rate at one specific edge. -#function compute_transport_rates(rate_law::Num, p_val_dict, lrs::LatticeReactionSystem) -function compute_transport_rates(rate_law, p_val_dict, lrs) - # Finds parameters involved in rate and create a function evaluating the rate law. + +# For the expression describing the rate of transport (likely only a single parameter, e.g. `D`), +# and the values of all our parameters, compute the transport rate(s). +# If all parameters that the rate depends on are uniform across all edges, this becomes a length-1 vector. +# Else it becomes a vector where each value corresponds to the rate at one specific edge. +function compute_transport_rates(s::BasicSymbolic{Real}, p_val_dict, lrs::LatticeReactionSystem) + # Find parameters involved in the rate and create a function evaluating the rate law. + rate_law = get_transport_rate_law(s, lrs) relevant_ps = Symbolics.get_variables(rate_law) rate_law_func = drop_expr(@RuntimeGeneratedFunction(build_function(rate_law, relevant_ps...))) - # If all these parameters are spatially uniform, the rates becomes a size (1,1) sparse matrix. - # Else, the rates becomes a size (num_verts,num_verts) sparse matrix. + # If all these parameters are spatially uniform, the rates become a size (1,1) sparse matrix. + # Else, the rates become a size (num_verts,num_verts) sparse matrix. if all(size(p_val_dict[p]) == (1,1) for p in relevant_ps) - relevant_p_vals = [get_edge_value(p_val_dict[p], [1 => 1]) for p in relevant_ps] + relevant_p_vals = [get_edge_value(p_val_dict[p], 1 => 1) for p in relevant_ps] return sparse([1],[1],rate_law_func(relevant_p_vals...)) else transport_rates = spzeros(lrs.num_verts, lrs.num_verts) @@ -203,76 +224,70 @@ function compute_transport_rates(rate_law, p_val_dict, lrs) end end -# Produces a dictionary with all parameters' values. Vertex parameters have their values converted to -# a sparse matrix (with one value for each edge, always using the source vertex's value) -function param_dict(vert_ps, edge_ps, lrs) - return merge(Dict(zip(vertex_parameters(lrs), vert_ps)), Dict(zip(edge_parameters(lrs), edge_ps))) +# For a species, retrieve the symbolic expression for its transportation rate +# (likely only a single parameter, such as `D`, but could be e.g. L*D, where L and D are parameters). +# If there are several transportation reactions for the species, their sum is used. +function get_transport_rate_law(s::BasicSymbolic{Real}, lrs::LatticeReactionSystem) + rates = filter(sr -> isequal(s, sr.species), lrs.spatial_reactions) + return sum(getfield.(rates, :rate)) end -### Accessing State & Parameter Array Values ### +### Accessing State & Parameter Array Values ### -# Converts a vector of vectors to a long vector. +# Converts a vector of vectors to a single, long, vector. # These are used when the initial condition is converted to a single vector (from vector of vector form). -function expand_component_values(values, num_verts) +function expand_component_values(values::Vector{Vector{T}}, num_verts::Int64) where {T} vcat([get_vertex_value.(values, vert) for vert in 1:num_verts]...) end -# Gets the index in the u array of species s in vertex vert (when their are num_species species). +# Gets the index in the u array of species s in vertex vert (when there are num_species species). get_index(vert::Int64, s::Int64, num_species::Int64) = (vert - 1) * num_species + s -# Gets the indexes in the u array of all species in vertex vert (when their are num_species species). +# Gets the indexes in the u array of all species in vertex vert (when there are num_species species). get_indexes(vert::Int64, num_species::Int64) = ((vert - 1) * num_species + 1):(vert * num_species) -# Returns the value of a parameter in an edge. For vertex parameters, uses their values in the source. -function get_edge_value(values::Vector{T}, edge) where {T} +# Returns the value of a parameter in an edge. For vertex parameters, use their values in the source. +function get_edge_value(values::Vector{T}, edge::Pair{Int64,Int64}) where {T} return (length(values) == 1) ? values[1] : values[edge[1]] end -function get_edge_value(values::SparseMatrixCSC{T, Int64}, edge) where {T} +function get_edge_value(values::SparseMatrixCSC{T, Int64}, edge::Pair{Int64,Int64}) where {T} return (size(values) == (1,1)) ? values[1,1] : values[edge[1],edge[2]] end -# Returns the value of an initial condition of parameter in a vertex. -function get_vertex_value(values::Vector{T}, vert_idx) where {T} +# Returns the value of an initial condition of vertex parameter in a vertex. +function get_vertex_value(values::Vector{T}, vert_idx::Int64) where {T} return (length(values) == 1) ? values[1] : values[vert_idx] end - - - - - - - - -# Finds the transport rate of a parameter going from a source vertex to a destination vertex. -function get_transport_rate(transport_rate, edge::Pair{Int64,Int64}, t_rate_idx_types::Bool) +# Finds the transport rate of a parameter along a specific edge. +function get_transport_rate(transport_rate::SparseMatrixCSC{T, Int64}, edge::Pair{Int64,Int64}, + t_rate_idx_types::Bool) where {T} return t_rate_idx_types ? transport_rate[1,1] : transport_rate[edge[1],edge[2]] end -# Finds the transportation rate for a specific species and a `LatticeTransportODEf` struct. +# Finds the transportation rate for a specific species, LatticeTransportODEf struct, and edge. function get_transport_rate(trans_s_idx::Int64, f_func::LatticeTransportODEf, edge::Pair{Int64,Int64}) get_transport_rate(f_func.transport_rates[trans_s_idx][2], edge, f_func.t_rate_idx_types[trans_s_idx]) end - - -# Updates the internal work_vert_ps vector for a given location. -# To this vector, we write the systems parameter values at a specific vertex. -function update_work_vert_ps!(work_vert_ps, vert_ps, comp, vert_ps_idx_types) +# Updates the internal work_vert_ps vector for a given vertex. +# To this vector, we write the system's parameter values at the specific vertex. +function update_work_vert_ps!(work_vert_ps::Vector{S}, all_ps::Vector{T}, vert::Int64, + vert_ps_idx_types::Vector{Bool}) where {S,T} # Loops through all parameters. for (idx,loc_type) in enumerate(vert_ps_idx_types) # If the parameter is uniform across the spatial structure, it will have a length-1 value vector # (which value we write to the work vector). # Else, we extract it value at the specific location. - work_vert_ps[idx] = (loc_type ? vert_ps[idx][1] : vert_ps[idx][comp]) + work_vert_ps[idx] = (loc_type ? all_ps[idx][1] : all_ps[idx][vert]) end end -# Input is always either a LatticeTransportODEf or LatticeTransportODEjac function (which fields we then pass on). -function update_work_vert_ps!(lt_ode_func, vert_ps, comp) - return update_work_vert_ps!(lt_ode_func.work_vert_ps, vert_ps, comp, lt_ode_func.v_ps_idx_types) +# Input is either a LatticeTransportODEf or LatticeTransportODEjac function (which fields we pass on). +function update_work_vert_ps!(lt_ode_func, all_ps::Vector{T}, vert::Int64) where {T} + return update_work_vert_ps!(lt_ode_func.work_vert_ps, all_ps, vert, lt_ode_func.v_ps_idx_types) end # Expands a u0/p information stored in Vector{Vector{}} for to Matrix form # (currently only used in Spatial Jump systems). -function matrix_expand_component_values(values::Vector{<:Vector}, n) +function matrix_expand_component_values(values::Vector{<:Vector}, n::Int64) reshape(expand_component_values(values, n), length(values), n) end diff --git a/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl b/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl index 992d864b20..74328d6979 100644 --- a/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl +++ b/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl @@ -23,7 +23,7 @@ let u0 = [:S => 990.0, :I => 20.0 * rand_v_vals(lrs.lattice), :R => 0.0] pV = SIR_p pE = [:dS => 0.01, :dI => 0.01, :dR => 0.01] - oprob = ODEProblem(lrs, u0, (0.0, 500.0), (pV, pE); jac = false) + oprob = ODEProblem(lrs, u0, (0.0, 500.0), [pV; pE]; jac = false) @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) runtime_target = 0.00027 @@ -38,7 +38,7 @@ let u0 = [:S => 990.0, :I => 20.0 * rand_v_vals(lrs.lattice), :R => 0.0] pV = SIR_p pE = [:dS => 0.01, :dI => 0.01, :dR => 0.01] - oprob = ODEProblem(lrs, u0, (0.0, 500.0), (pV, pE); jac = false) + oprob = ODEProblem(lrs, u0, (0.0, 500.0), [pV; pE]; jac = false) @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) runtime_target = 0.12 @@ -53,7 +53,7 @@ let u0 = [:X => rand_v_vals(lrs.lattice, 10), :Y => rand_v_vals(lrs.lattice, 10)] pV = brusselator_p pE = [:dX => 0.2] - oprob = ODEProblem(lrs, u0, (0.0, 100.0), (pV, pE)) + oprob = ODEProblem(lrs, u0, (0.0, 100.0), [pV; pE]) @test SciMLBase.successful_retcode(solve(oprob, CVODE_BDF(linear_solver=:GMRES))) runtime_target = 0.013 @@ -68,7 +68,7 @@ let u0 = [:X => rand_v_vals(lrs.lattice, 10), :Y => rand_v_vals(lrs.lattice, 10)] pV = brusselator_p pE = [:dX => 0.2] - oprob = ODEProblem(lrs, u0, (0.0, 100.0), (pV, pE)) + oprob = ODEProblem(lrs, u0, (0.0, 100.0), [pV; pE]) @test SciMLBase.successful_retcode(solve(oprob, CVODE_BDF(linear_solver=:GMRES))) runtime_target = 11. @@ -99,7 +99,7 @@ let ] pV = CuH_Amination_p pE = [:D1 => 0.1, :D2 => 0.1, :D3 => 0.1, :D4 => 0.1, :D5 => 0.1] - oprob = ODEProblem(lrs, u0, (0.0, 10.0), (pV, pE); jac = false) + oprob = ODEProblem(lrs, u0, (0.0, 10.0), [pV; pE]; jac = false) @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) runtime_target = 0.0012 @@ -130,7 +130,7 @@ let ] pV = CuH_Amination_p pE = [:D1 => 0.1, :D2 => 0.1, :D3 => 0.1, :D4 => 0.1, :D5 => 0.1] - oprob = ODEProblem(lrs, u0, (0.0, 10.0), (pV, pE); jac = false) + oprob = ODEProblem(lrs, u0, (0.0, 10.0), [pV; pE]; jac = false) @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) runtime_target = 0.56 @@ -156,7 +156,7 @@ let ] pV = sigmaB_p pE = [:DσB => 0.1, :Dw => 0.1, :Dv => 0.1] - oprob = ODEProblem(lrs, u0, (0.0, 50.0), (pV, pE)) + oprob = ODEProblem(lrs, u0, (0.0, 50.0), [pV; pE]) @test SciMLBase.successful_retcode(solve(oprob, CVODE_BDF(linear_solver=:GMRES))) runtime_target = 0.61 @@ -182,7 +182,7 @@ let ] pV = sigmaB_p pE = [:DσB => 0.1, :Dw => 0.1, :Dv => 0.1] - oprob = ODEProblem(lrs, u0, (0.0, 10.0), (pV, pE)) # Time reduced from 50.0 (which casues Julai to crash). + oprob = ODEProblem(lrs, u0, (0.0, 10.0), [pV; pE]) # Time reduced from 50.0 (which casues Julai to crash). @test SciMLBase.successful_retcode(solve(oprob, CVODE_BDF(linear_solver=:GMRES))) runtime_target = 59. diff --git a/test/spatial_reaction_systems/lattice_reaction_systems_ODEs.jl b/test/spatial_reaction_systems/lattice_reaction_systems_ODEs.jl index 4397658cee..050044109f 100644 --- a/test/spatial_reaction_systems/lattice_reaction_systems_ODEs.jl +++ b/test/spatial_reaction_systems/lattice_reaction_systems_ODEs.jl @@ -300,7 +300,8 @@ let end u0 = 2 * rand(rng, 10000) - p = [1.0, 4.0, 0.1] + p_hw = [1.0, 4.0, 0.1] + p_aut = [[1.0], [4.0], sparse([1], [1], [0.1])] tspan = (0.0, 100.0) ofun_hw_dense = ODEFunction(spatial_brusselator_f; jac = spatial_brusselator_jac) @@ -320,10 +321,10 @@ let du_aut_dense = deepcopy(u0) du_aut_sparse = deepcopy(u0) - ofun_hw_dense(du_hw_dense, u0, p, 0.0) - ofun_hw_sparse(du_hw_sparse, u0, p, 0.0) - ofun_aut_dense(du_aut_dense, u0, p, 0.0) - ofun_aut_sparse(du_aut_sparse, u0, p, 0.0) + ofun_hw_dense(du_hw_dense, u0, p_hw, 0.0) + ofun_hw_sparse(du_hw_sparse, u0, p_hw, 0.0) + ofun_aut_dense(du_aut_dense, u0, p_aut, 0.0) + ofun_aut_sparse(du_aut_sparse, u0, p_aut, 0.0) @test isapprox(du_hw_dense, du_aut_dense) @test isapprox(du_hw_sparse, du_aut_sparse)