diff --git a/HISTORY.md b/HISTORY.md index 2bb952cd00..069e74ae5d 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -1,9 +1,32 @@ # Breaking updates and feature summaries across releases ## Catalyst unreleased (master branch) +- Simulation of spatial ODEs now supported. For full details, please see /~https://github.com/SciML/Catalyst.jl/pull/644 and upcoming documentation. Note that these methods are currently considered alpha, with the interface and approach changing even in non-breaking Catalyst releases. +- LatticeReactionSystem structure represents a spatial reaction network: + ```julia + rn = @reaction_network begin + (p,d), 0 <--> X + end + tr = @transport_reaction D X + lattice = Graphs.grid([5, 5]) + lrs = LatticeReactionSystem(rn, [tr], lattice) +``` +- Here, if a `u0` or `p` vector is given with scalar values: + ```julia + u0 = [:X => 1.0] + p = [:p => 1.0, :d => 0.5, :D => 0.1] + ``` + this value will be used across the entire system. If their values are instead vectors, different values are used across the spatial system. Here + ```julia + X0 = zeros(25) + X0[1] = 1.0 + u0 = [:X => X0] + ``` + X's value will be `1.0` in the first vertex, but `0.0` in the remaining one (the system have 25 vertexes in total). SInce th parameters `p` and `d` are part of the non-spatial reaction network, their values are tied to vertexes. However, if the `D` parameter (which governs diffusion between vertexes) is given several values, these will instead correspond to the specific edges (and transportation along those edges.) + ## Catalyst 13.5 -- Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reactin system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example: +- Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example: ```julia wilhelm_2009_model = @reaction_network begin k1, Y --> 2X diff --git a/Project.toml b/Project.toml index 6300c1427e..c43b053c36 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" @@ -40,6 +41,7 @@ ModelingToolkit = "8.66" Parameters = "0.12" Reexport = "0.2, 1.0" Requires = "1.0" +RuntimeGeneratedFunctions = "0.5.12" SymbolicUtils = "1.0.3" Symbolics = "5.0.3" Unitful = "1.12.4" diff --git a/src/Catalyst.jl b/src/Catalyst.jl index b33b6782d0..f10795d6e9 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -16,6 +16,10 @@ const MT = ModelingToolkit using Unitful @reexport using ModelingToolkit using Symbolics + +using RuntimeGeneratedFunctions +RuntimeGeneratedFunctions.init(@__MODULE__) + import Symbolics: BasicSymbolic import SymbolicUtils using ModelingToolkit: Symbolic, value, istree, get_states, get_ps, get_iv, get_systems, @@ -33,6 +37,7 @@ import ModelingToolkit: check_variables, import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show import MacroTools, Graphs +import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne import DataStructures: OrderedDict, OrderedSet import Parameters: @with_kw_noshow @@ -107,4 +112,22 @@ export balance_reaction function hc_steady_states end export hc_steady_states +### Spatial Reaction Networks ### + +# spatial reactions +include("spatial_reaction_systems/spatial_reactions.jl") +export TransportReaction, TransportReactions, @transport_reaction +export isedgeparameter + +# lattice reaction systems +include("spatial_reaction_systems/lattice_reaction_systems.jl") +export LatticeReactionSystem +export spatial_species, vertex_parameters, edge_parameters + +# variosu utility functions +include("spatial_reaction_systems/utility.jl") + +# spatial lattice ode systems. +include("spatial_reaction_systems/spatial_ODE_systems.jl") + end # module diff --git a/src/spatial_reaction_systems/lattice_reaction_systems.jl b/src/spatial_reaction_systems/lattice_reaction_systems.jl new file mode 100644 index 0000000000..79cb5d4fd0 --- /dev/null +++ b/src/spatial_reaction_systems/lattice_reaction_systems.jl @@ -0,0 +1,91 @@ +### Lattice Reaction Network Structure ### +# Describes a spatial reaction network over a graph. +# Adding the "<: MT.AbstractTimeDependentSystem" part messes up show, disabling me from creating LRSs. +struct LatticeReactionSystem{S,T} # <: MT.AbstractTimeDependentSystem + # Input values. + """The reaction system within each compartment.""" + rs::ReactionSystem{S} + """The spatial reactions defined between individual nodes.""" + spatial_reactions::Vector{T} + """The graph on which the lattice is defined.""" + lattice::SimpleDiGraph{Int64} + + # Derived values. + """The number of compartments.""" + num_verts::Int64 + """The number of edges.""" + num_edges::Int64 + """The number of species.""" + num_species::Int64 + """Whenever the initial input was a digraph.""" + init_digraph::Bool + """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). + """ + 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. + """ + 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. + """ + edge_parameters::Vector{BasicSymbolic{Real}} + + function LatticeReactionSystem(rs::ReactionSystem{S}, spatial_reactions::Vector{T}, + lattice::DiGraph; init_digraph = true) where {S, T} + # There probably some better way to ascertain that T has that type. Not sure how. + if !(T <: AbstractSpatialReaction) + error("The second argument must be a vector of AbstractSpatialReaction subtypes.") + end + + 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])) + rs_edge_parameters = filter(isedgeparameter, parameters(rs)) + if isempty(spatial_reactions) + srs_edge_parameters = Vector{BasicSymbolic{Real}}[] + else + srs_edge_parameters = setdiff(reduce(vcat, [parameters(sr) for sr in spatial_reactions]), parameters(rs)) + end + edge_parameters = unique([rs_edge_parameters; srs_edge_parameters]) + vertex_parameters = filter(!isedgeparameter, parameters(rs)) + # Ensures the parameter order begins similarly to in the non-spatial ReactionSystem. + ps = [parameters(rs); setdiff([edge_parameters; vertex_parameters], parameters(rs))] + + foreach(sr -> check_spatial_reaction_validity(rs, sr; edge_parameters=edge_parameters), spatial_reactions) + return new{S,T}(rs, spatial_reactions, lattice, nv(lattice), ne(lattice), num_species, + init_digraph, spat_species, ps, vertex_parameters, edge_parameters) + end +end +function LatticeReactionSystem(rs, srs, lat::SimpleGraph) + return LatticeReactionSystem(rs, srs, DiGraph(lat); init_digraph = false) +end + +### Lattice ReactionSystem Getters ### + +# Get all species. +species(lrs::LatticeReactionSystem) = unique([species(lrs.rs); lrs.spat_species]) +# Get all species that may be transported. +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). +vertex_parameters(lrs::LatticeReactionSystem) = lrs.vertex_parameters +# Get all parameters which values are tied to edges (adjacencies). +edge_parameters(lrs::LatticeReactionSystem) = lrs.edge_parameters + +# Gets the lrs name (same as rs name). +ModelingToolkit.nameof(lrs::LatticeReactionSystem) = nameof(lrs.rs) + +# Checks if a lattice reaction system is a pure (linear) transport reaction system. +is_transport_system(lrs::LatticeReactionSystem) = all(sr -> sr isa TransportReaction, lrs.spatial_reactions) diff --git a/src/spatial_reaction_systems/spatial_ODE_systems.jl b/src/spatial_reaction_systems/spatial_ODE_systems.jl new file mode 100644 index 0000000000..f0d71383ba --- /dev/null +++ b/src/spatial_reaction_systems/spatial_ODE_systems.jl @@ -0,0 +1,304 @@ +### Spatial ODE Functor Structures ### + +# Functor with information for the forcing function of a spatial ODE with spatial movement on a lattice. +struct LatticeTransportODEf{Q,R,S,T} + """The ODEFunction of the (non-spatial) reaction system which generated this function.""" + ofunc::Q + """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.""" + vert_ps::Vector{Vector{R}} + """ + 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. + """ + work_vert_ps::Vector{R} + """ + 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. + This field (specifically) is an enumerate over that array. + """ + v_ps_idx_types::Vector{Bool} + """ + A vector of pairs, 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. + """ + transport_rates::Vector{Pair{Int64, Vector{R}}} + """ + A matrix, NxM, where N is the number of species with transportation and M 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). + """ + leaving_rates::Matrix{R} + """An (enumerate'ed) iterator over all the edges of the lattice.""" + edges::S + """ + The edge parameters used to create the spatial ODEProblem. Currently unused, + but will be needed to support changing these (e.g. due to events). + Contain one vector for each edge parameter (length one if uniform, else one value for each edge). + """ + edge_ps::Vector{Vector{T}} + + function LatticeTransportODEf(ofunc::Q, vert_ps::Vector{Vector{R}}, transport_rates::Vector{Pair{Int64, Vector{R}}}, + edge_ps::Vector{Vector{T}}, lrs::LatticeReactionSystem) where {Q,R,T} + leaving_rates = zeros(length(transport_rates), lrs.num_verts) + for (s_idx, trpair) in enumerate(transport_rates) + rates = last(trpair) + for (e_idx, e) in enumerate(edges(lrs.lattice)) + # Updates the exit rate for species s_idx from vertex e.src + leaving_rates[s_idx, e.src] += get_component_value(rates, e_idx) + end + end + work_vert_ps = zeros(lrs.num_verts) + # 1 if ps are constant across the graph, 0 else. + v_ps_idx_types = map(vp -> length(vp) == 1, vert_ps) + eds = edges(lrs.lattice) + new{Q,R,typeof(eds),T}(ofunc, lrs.num_verts, lrs.num_species, vert_ps, work_vert_ps, + v_ps_idx_types, transport_rates, leaving_rates, eds, edge_ps) + end +end + +# Functor with information for the Jacobian function of a spatial ODE with spatial movement on a lattice. +struct LatticeTransportODEjac{Q,R,S,T} + """The ODEFunction of the (non-spatial) reaction system which generated this function.""" + ofunc::Q + """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.""" + vert_ps::Vector{Vector{R}} + """ + 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. + """ + work_vert_ps::Vector{R} + """ + 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. + """ + v_ps_idx_types::Vector{Bool} + """Whether the Jacobian is sparse or not.""" + sparse::Bool + """The transport rates. Can be a dense matrix (for non-sparse) or as the "nzval" field if sparse.""" + jac_transport::S + """ + The edge parameters used to create the spatial ODEProblem. Currently unused, + but will be needed to support changing these (e.g. due to events). + Contain one vector for each edge parameter (length one if uniform, else one value for each edge). + """ + edge_ps::Vector{Vector{T}} + + function LatticeTransportODEjac(ofunc::R, vert_ps::Vector{Vector{S}}, lrs::LatticeReactionSystem, + jac_transport::Union{Nothing, SparseMatrixCSC{Float64, Int64}}, + edge_ps::Vector{Vector{T}}, sparse::Bool) where {R,S,T} + work_vert_ps = zeros(lrs.num_verts) + v_ps_idx_types = map(vp -> length(vp) == 1, vert_ps) + new{R,S,typeof(jac_transport),T}(ofunc, lrs.num_verts, lrs.num_species, vert_ps, + work_vert_ps, v_ps_idx_types, sparse, jac_transport, edge_ps) + end +end + +### ODEProblem ### + +# Creates an ODEProblem from a LatticeReactionSystem. +function DiffEqBase.ODEProblem(lrs::LatticeReactionSystem, u0_in, tspan, + p_in = DiffEqBase.NullParameters(), args...; + jac = false, sparse = false, + name = nameof(lrs), include_zero_odes = true, + combinatoric_ratelaws = get_combinatoric_ratelaws(lrs.rs), + remove_conserved = false, checks = false, kwargs...) + if !is_transport_system(lrs) + 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. + u0_in = symmap_to_varmap(lrs, u0_in) + p_in = (p_in isa Tuple{<:Any,<:Any}) ? + (symmap_to_varmap(lrs, p_in[1]),symmap_to_varmap(lrs, p_in[2])) : + symmap_to_varmap(lrs, p_in) + + # Converts u0 and p to their internal forms. + # 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.num_verts) + # Both vert_ps and edge_ps becomes vectors of vectors. Each have 1 element for each parameter. + # These elements are length 1 vectors (if the parameter is uniform), + # or length num_verts/nE, with unique values for each vertex/edge (for vert_ps/edge_ps, respectively). + vert_ps, edge_ps = lattice_process_p(p_in, vertex_parameters(lrs), edge_parameters(lrs), lrs) + + # Creates ODEProblem. + ofun = build_odefunction(lrs, vert_ps, edge_ps, jac, sparse, name, include_zero_odes, + combinatoric_ratelaws, remove_conserved, checks) + return ODEProblem(ofun, u0, tspan, vert_ps, args...; kwargs...) +end + +# Builds an ODEFunction for a spatial ODEProblem. +function build_odefunction(lrs::LatticeReactionSystem, vert_ps::Vector{Vector{T}}, + edge_ps::Vector{Vector{T}}, 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). + transport_rates = make_sidxs_to_transrate_map(vert_ps, edge_ps, lrs) + + # 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. + 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) + if sparse + f = LatticeTransportODEf(ofunc_sparse, vert_ps, transport_rates, edge_ps, 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, edge_ps, true) + jac_prototype = jac_vals + else + f = LatticeTransportODEf(ofunc_dense, vert_ps, transport_rates, edge_ps, lrs) + J = LatticeTransportODEjac(ofunc_dense, vert_ps, lrs, jac_vals, edge_ps, false) + jac_prototype = nothing + end + else + if sparse + ofunc_sparse = ODEFunction(osys; jac = false, sparse = true) + f = LatticeTransportODEf(ofunc_sparse, vert_ps, transport_rates, edge_ps, lrs) + jac_prototype = build_jac_prototype(ofunc_sparse.jac_prototype, transport_rates, lrs; set_nonzero = false) + else + ofunc_dense = ODEFunction(osys; jac = false, sparse = false) + f = LatticeTransportODEf(ofunc_dense, vert_ps, transport_rates, edge_ps, lrs) + jac_prototype = nothing + end + J = nothing + end + + 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}, trans_rates, + lrs::LatticeReactionSystem; set_nonzero = false) + # Finds the indexes of the transport species, and the species with transport only (and no non-spatial dynamics). + trans_species = first.(trans_rates) + trans_only_species = filter(s_idx -> !Base.isstored(ns_jac_prototype, s_idx, s_idx), trans_species) + + # Finds the indexes of all terms in the non-spatial jacobian. + ns_jac_prototype_idxs = findnz(ns_jac_prototype) + ns_i_idxs = ns_jac_prototype_idxs[1] + ns_j_idxs = ns_jac_prototype_idxs[2] + + # Prepares vectors to store i and j indexes of Jacobian entries. + idx = 1 + num_entries = lrs.num_verts * length(ns_i_idxs) + + lrs.num_edges * (length(trans_only_species) + length(trans_species)) + i_idxs = Vector{Int}(undef, num_entries) + j_idxs = Vector{Int}(undef, num_entries) + + # Indexes of elements due to 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) + j_idxs[idx] = get_index(vert, ns_j_idxs[n], lrs.num_species) + idx += 1 + end + end + + # Indexes of elements due to spatial dynamics. + for e in edges(lrs.lattice) + # Indexes due to terms for a species leaves its current 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 + i_idxs[idx] = get_index(e.src, s_idx, lrs.num_species) + j_idxs[idx] = i_idxs[idx] + idx += 1 + end + # Indexes due to terms for species arriving into a new vertex. + for s_idx in trans_species + i_idxs[idx] = get_index(e.src, s_idx, lrs.num_species) + j_idxs[idx] = get_index(e.dst, s_idx, lrs.num_species) + idx += 1 + end + end + + # Create sparse jacobian prototype with 0-valued entries. + jac_prototype = sparse(i_idxs, j_idxs, zeros(num_entries)) + + # Set element values. + if set_nonzero + for (s, rates) in trans_rates, (e_idx, e) in enumerate(edges(lrs.lattice)) + idx_src = get_index(e.src, s, lrs.num_species) + idx_dst = get_index(e.dst, s, lrs.num_species) + val = get_component_value(rates, e_idx) + + # Term due to species leaving source vertex. + jac_prototype[idx_src, idx_src] -= val + + # Term due to species arriving to destination vertex. + jac_prototype[idx_src, idx_dst] += val + end + end + + return jac_prototype +end + +# Defines the forcing functor's effect on the (spatial) ODE system. +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 + idxs = get_indexes(vert_i, f_func.num_species) + + # vector of vertex ps at vert_i + vert_i_ps = view_vert_ps_vector!(f_func.work_vert_ps, p, vert_i, enumerate(f_func.v_ps_idx_types)) + + # evaluate reaction contributions to du at vert_i + f_func.ofunc((@view du[idxs]), (@view u[idxs]), vert_i_ps, t) + end + + # s_idx is species index among transport species, s is 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 + for vert_i in 1:(f_func.num_verts) + idx = get_index(vert_i, s, f_func.num_species) + du[idx] -= f_func.leaving_rates[s_idx, vert_i] * u[idx] + end + # Add rates for entering a given vertex via an incoming edge + for (e_idx, e) in enumerate(f_func.edges) + idx_dst = get_index(e.dst, s, f_func.num_species) + idx_src = get_index(e.src, s, f_func.num_species) + du[idx_dst] += get_component_value(rates, e_idx) * u[idx_src] + end + end +end + +# 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 + for vert_i in 1:(jac_func.num_verts) + idxs = get_indexes(vert_i, jac_func.num_species) + vert_ps = view_vert_ps_vector!(jac_func.work_vert_ps, p, vert_i, enumerate(jac_func.v_ps_idx_types)) + jac_func.ofunc.jac((@view J[idxs, idxs]), (@view u[idxs]), vert_ps, t) + end + + # Updates for the spatial reactions (adds the Jacobian values from the diffusion reactions). + J .+= jac_func.jac_transport +end \ No newline at end of file diff --git a/src/spatial_reaction_systems/spatial_reactions.jl b/src/spatial_reaction_systems/spatial_reactions.jl new file mode 100644 index 0000000000..c1a01a3888 --- /dev/null +++ b/src/spatial_reaction_systems/spatial_reactions.jl @@ -0,0 +1,152 @@ +### Spatial Reaction Structures ### + +# Abstract spatial reaction structures. +abstract type AbstractSpatialReaction end + +### EdgeParameter Metadata ### + +# Implements the edgeparameter metadata field. +struct EdgeParameter end +Symbolics.option_to_metadata_type(::Val{:edgeparameter}) = EdgeParameter + +# Implements the isedgeparameter check function. +isedgeparameter(x::Num, args...) = isedgeparameter(Symbolics.unwrap(x), args...) +function isedgeparameter(x, default = false) + p = Symbolics.getparent(x, nothing) + p === nothing || (x = p) + 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""" + rate::Any + """The species that is subject to diffusion.""" + species::BasicSymbolic{Real} + + # Creates a diffusion reaction. + function TransportReaction(rate, species) + if any(!ModelingToolkit.isparameter(var) for var in ModelingToolkit.get_variables(rate)) + error("TransportReaction rate contains variables: $(filter(var -> !ModelingToolkit.isparameter(var), ModelingToolkit.get_variables(rate))). The rate must consist of parameters only.") + end + new(rate, species.val) + end +end +# Creates a vector of TransportReactions. +function TransportReactions(transport_reactions) + [TransportReaction(tr[1], tr[2]) for tr in transport_reactions] +end + +# Macro for creating a transport reaction. +macro transport_reaction(rateex::ExprValues, species::ExprValues) + make_transport_reaction(MacroTools.striplines(rateex), species) +end +function make_transport_reaction(rateex, species) + # Handle interpolation of variables + rateex = esc_dollars!(rateex) + species = esc_dollars!(species) + + # Parses input expression. + parameters = ExprValues[] + find_parameters_in_rate!(parameters, rateex) + + # Checks for input errors. + forbidden_symbol_check(union([species], parameters)) + + # Creates expressions corresponding to actual code from the internal DSL representation. + sexprs = get_sexpr([species], Dict{Symbol, Expr}()) + pexprs = get_pexpr(parameters, Dict{Symbol, Expr}()) + iv = :(@variables $(DEFAULT_IV_SYM)) + trxexpr = :(TransportReaction($rateex, $species)) + + quote + $pexprs + $iv + $sexprs + $trxexpr + end +end + +# Gets the parameters in a transport reaction. +ModelingToolkit.parameters(tr::TransportReaction) = Symbolics.get_variables(tr.rate) + +# Gets the species in a transport reaction. +spatial_species(tr::TransportReaction) = [tr.species] + +# Checks that a transport reaction 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). + 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 + + # Checks that the rate does not depend on species. + rate_vars = ModelingToolkit.getname.(Symbolics.get_variables(tr.rate)) + if !isempty(intersect(ModelingToolkit.getname.(species(rs)), rate_vars)) + error("The following species were used in rates of a transport reactions: $(setdiff(ModelingToolkit.getname.(species(rs)), rate_vars)).") + end + + # 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.") + 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. + 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.") + end +end +equivalent_metadata(p1, p2) = isempty(setdiff(p1.metadata, p2.metadata, [Catalyst.EdgeParameter => true])) + +# Since MTK's "isequal" ignores metadata, we have to use a special function that accounts for this. +# This is important because whether something is an edge parameter is defined in metadata. +function isequivalent(sym1, sym2) + !isequal(sym1, sym2) && (return false) + (sym1.metadata != sym2.metadata) && (return false) + return true +end + +# Implements custom `==`. +""" + ==(rx1::TransportReaction, rx2::TransportReaction) + +Tests whether two [`TransportReaction`](@ref)s are identical. +""" +function (==)(tr1::TransportReaction, tr2::TransportReaction) + isequal(tr1.rate, tr2.rate) || return false + isequal(tr1.species, tr2.species) || return false + return true +end + +# Implements custom `hash`. +function hash(tr::TransportReaction, h::UInt) + h = Base.hash(tr.rate, h) + Base.hash(tr.species, h) +end + +### Utility ### +# Loops through a rate and extract all parameters. +function find_parameters_in_rate!(parameters, rateex::ExprValues) + if rateex isa Symbol + if rateex in [:t, :∅, :im, :nothing, CONSERVED_CONSTANT_SYMBOL] + error("Forbidden term $(rateex) used in transport reaction rate.") + elseif !(rateex in [:ℯ, :pi, :π]) + push!(parameters, rateex) + end + elseif rateex isa Expr + # Note, this (correctly) skips $(...) expressions + for i in 2:length(rateex.args) + find_parameters_in_rate!(parameters, rateex.args[i]) + end + end + nothing +end \ No newline at end of file diff --git a/src/spatial_reaction_systems/utility.jl b/src/spatial_reaction_systems/utility.jl new file mode 100644 index 0000000000..8e6543a2b3 --- /dev/null +++ b/src/spatial_reaction_systems/utility.jl @@ -0,0 +1,297 @@ +### Processes Input u0 & p ### + +# Defines _symbol_to_var, but where the system is a LRS. Required to make symmapt_to_varmap to work. +function _symbol_to_var(lrs::LatticeReactionSystem, sym) + # Checks if sym is a parameter. + p_idx = findfirst(sym==p_sym for p_sym in ModelingToolkit.getname.(parameters(lrs))) + isnothing(p_idx) || return parameters(lrs)[p_idx] + + # Checks if sym is a species. + s_idx = findfirst(sym==s_sym for s_sym in ModelingToolkit.getname.(species(lrs))) + isnothing(s_idx) || return species(lrs)[s_idx] + + 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, num_verts) + # u0 values can be given in various forms. This converts it to a Vector{Vector{}} form. + # Top-level vector: Contains one vector for each species. + # Second-level vector: contain one value if species uniform across lattice, else one value for each vertex). + u0 = lattice_process_input(u0_in, u0_syms, num_verts) + + # Perform various error checks on the (by the user provided) initial conditions. + check_vector_lengths(u0, length(u0_syms), num_verts) + + # Converts the Vector{Vector{}} format to a single Vector (with one values for each species and vertex). + expand_component_values(u0, num_verts) +end + +# From p input, splits it into diffusion parameters and compartment parameters. +# Store these in the desired internal format. +function lattice_process_p(p_in, p_vertex_syms, p_edge_syms, lrs::LatticeReactionSystem) + # If the user provided parameters as a single map (mixing vertex and edge parameters): + # Split into two separate vectors. + vert_ps_in, edge_ps_in = split_parameters(p_in, p_vertex_syms, p_edge_syms) + + # Parameter values can be given in various forms. This converts it to the Vector{Vector{}} form. + vert_ps = lattice_process_input(vert_ps_in, p_vertex_syms, lrs.num_verts) + + # Parameter values can be given in various forms. This converts it to the Vector{Vector{}} form. + edge_ps = lattice_process_input(edge_ps_in, p_edge_syms, lrs.num_edges) + + # If the lattice defined as (N edge) undirected graph, and we provides N/2 values for some edge parameter: + # Presume they want to expand that parameters value so it has the same value in both directions. + lrs.init_digraph || duplicate_trans_params!(edge_ps, lrs) + + # Perform various error checks on the (by the user provided) vertex and edge parameters. + check_vector_lengths(vert_ps, length(p_vertex_syms), lrs.num_verts) + check_vector_lengths(edge_ps, length(p_edge_syms), lrs.num_edges) + + return vert_ps, edge_ps +end + +# Splits parameters into those for the vertexes and those for the edges. + +# If they are already split, return that. +split_parameters(ps::Tuple{<:Any, <:Any}, args...) = ps +# Providing parameters to a spatial reaction system as a single vector of values (e.g. [1.0, 4.0, 0.1]) is not allowed. +# Either use tuple (e.g. ([1.0, 4.0], [0.1])) or map format (e.g. [A => 1.0, B => 4.0, D => 0.1]). +function split_parameters(ps::Vector{<:Number}, args...) + error("When providing parameters for a spatial system as a single vector, the paired form (e.g :D =>1.0) must be used.") +end +# Splitting is only done for Vectors of Pairs (where the first value is a Symbols, and the second a value). +function split_parameters(ps::Vector{<: Pair}, p_vertex_syms::Vector, p_edge_syms::Vector) + vert_ps_in = [p for p in ps if any(isequal(p[1]), p_vertex_syms)] + edge_ps_in = [p for p in ps if any(isequal(p[1]), p_edge_syms)] + + # Error check, in case some input parameters where neither recognised as vertex or edge parameters. + if (sum(length.([vert_ps_in, edge_ps_in])) != length(ps)) + error("These input parameters are not recognised: $(setdiff(first.(ps), vcat(first.([vert_ps_in, edge_ps_in]))))") + end + + return vert_ps_in, edge_ps_in +end + +# Input may have the following forms (after potential Symbol maps to Symbolic maps conversions): + # - A vector of values, where the i'th value corresponds to the value of the i'th + # initial condition value (for u0_in), vertex parameter value (for vert_ps_in), or edge parameter value (for edge_ps_in). + # - A vector of vectors of values. The same as previously, + # but here the species/parameter can have different values across the spatial structure. + # - A map of Symbols to values. These can either be a single value (if uniform across the spatial structure) + # or a vector (with different values for each vertex/edge). + # These can be mixed (e.g. [X => 1.0, Y => [1.0, 2.0, 3.0, 4.0]] is allowed). + # - A matrix. E.g. for initial conditions you can have a num_species * num_vertex matrix, + # indicating the value of each species at each vertex. + +# The lattice_process_input function takes input initial conditions/vertex parameters/edge parameters +# of whichever form the user have used, and converts them to the Vector{Vector{}} form used internally. +# E.g. for parameters the top-level vector contain one vector for each parameter (same order as in parameters(::ReactionSystem)). +# If a parameter is uniformly-values across the spatial structure, its vector has a single value. +# Else, it has a number of values corresponding to the number of vertexes/edges (for edge/vertex parameters). +# Initial conditions works similarly. + +# If the input is given in a map form, the vector needs sorting and the first value removed. +# The creates a Vector{Vector{Value}} or Vector{value} form, which is then again sent to lattice_process_input for reprocessing. +function lattice_process_input(input::Vector{<:Pair}, syms::Vector{BasicSymbolic{Real}}, args...) + if !isempty(setdiff(first.(input), syms)) + error("Some input symbols are not recognised: $(setdiff(first.(input), syms)).") + end + sorted_input = sort(input; by = p -> findfirst(isequal(p[1]), syms)) + return lattice_process_input(last.(sorted_input), syms, args...) +end +# If the input is a matrix: Processes the input and gives it in a form where it is a vector of vectors +# (some of which may have a single value). Sends it back to lattice_process_input for reprocessing. +function lattice_process_input(input::Matrix{<:Number}, args...) + lattice_process_input([vec(input[i, :]) for i in 1:size(input, 1)], args...) +end +# Possibly we want to support this type of input at some point. +function lattice_process_input(input::Array{<:Number, 3}, args...) + error("3 dimensional array parameter input currently not supported.") +end +# If the input is a Vector containing both vectors and single values, converts it to the Vector{<:Vector} form. +# Technically this last lattice_process_input is probably not needed. +function lattice_process_input(input::Vector{<:Any}, args...) + isempty(input) ? Vector{Vector{Float64}}() : + lattice_process_input([(val isa Vector{<:Number}) ? val : [val] for val in input], + args...) +end +# If the input is of the correct form already, return it. +lattice_process_input(input::Vector{<:Vector}, syms::Vector{BasicSymbolic{Real}}, n::Int64) = input + +# Checks that a value vector have the right length, as well as that of all its sub vectors. +# Error check if e.g. the user does not provide values for all species/parameters, +# or for one: provides a vector of values, but that has the wrong length +# (e.g providing 7 values for one species, but there are 8 vertexes). +function check_vector_lengths(input::Vector{<:Vector}, n_syms, n_locations) + if (length(input)!=n_syms) + error("Missing values for some initial conditions/parameters. Expected $n_syms values, got $(length(input)).") + end + if !isempty(setdiff(unique(length.(input)), [1, n_locations])) + error("Some inputs where given values of inappropriate length.") + end +end + +# For transport parameters, if the lattice was given as an undirected graph of size n: +# this is converted to a directed graph of size 2n. +# If transport parameters are given with n values, we want to use the same value for both directions. +# Since the order of edges in the new graph is non-trivial, this function +# distributes the n input values to a 2n length vector, putting the correct value in each position. +function duplicate_trans_params!(edge_ps::Vector{Vector{Float64}}, lrs::LatticeReactionSystem) + cum_adjacency_counts = [0;cumsum(length.(lrs.lattice.fadjlist[1:end-1]))] + for idx in 1:length(edge_ps) + # If the edge parameter already values for each directed edge, we can continue. + (2length(edge_ps[idx]) == lrs.num_edges) || continue # + + # This entire thing depends on the fact that, in the edges(lattice) iterator, the edges are sorted by: + # (1) Their source node + # (2) Their destination node. + + # A vector where we will put the edge parameters new values. + # Has the correct length (the number of directed edges in the lattice). + new_vals = Vector{Float64}(undef, lrs.num_edges) + # As we loop through the edges of the di-graph, this keeps track of each edge's index in the original graph. + original_edge_count = 0 + for edge in edges(lrs.lattice) # For each edge. + # The digraph conversion only adds edges so that src > dst. + (edge.src < edge.dst) ? (original_edge_count += 1) : continue + # For original edge i -> j, finds the index of i -> j in DiGraph. + idx_fwd = cum_adjacency_counts[edge.src] + findfirst(isequal(edge.dst),lrs.lattice.fadjlist[edge.src]) + # For original edge i -> j, finds the index of j -> i in DiGraph. + idx_bwd = cum_adjacency_counts[edge.dst] + findfirst(isequal(edge.src),lrs.lattice.fadjlist[edge.dst]) + new_vals[idx_fwd] = edge_ps[idx][original_edge_count] + new_vals[idx_bwd] = edge_ps[idx][original_edge_count] + end + # Replaces the edge parameters values with the updated value vector. + edge_ps[idx] = new_vals + end +end + +# For a set of input values on the given forms, and their symbolics, convert into a dictionary. +vals_to_dict(syms::Vector, vals::Vector{<:Vector}) = Dict(zip(syms, vals)) +# Produces a dictionary with all parameter values. +function param_dict(vert_ps, edge_ps, lrs) + merge(vals_to_dict(vertex_parameters(lrs), vert_ps), + vals_to_dict(edge_parameters(lrs), edge_ps)) +end + +# Computes the transport rates and stores them in a desired format +# (a Dictionary from species index to rates across all edges). +function compute_all_transport_rates(vert_ps::Vector{Vector{Float64}}, edge_ps::Vector{Vector{Float64}}, lrs::LatticeReactionSystem) + # Creates a dict, allowing us to access the values of wll parameters. + p_val_dict = param_dict(vert_ps, edge_ps, lrs) + + # 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.num_edges) + for s in spatial_species(lrs)] + + # Sorts all the species => rate pairs according to their species index in species(::ReactionSystem). + 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). +# We could allows several transportation reactions for one species and simply sum them though, easy change. +function get_transport_rate_law(s::BasicSymbolic{Real}, lrs::LatticeReactionSystem) + rates = filter(sr -> isequal(s, sr.species), lrs.spatial_reactions) + (length(rates) > 1) && error("Species $s have more than one diffusion reaction.") + return rates[1].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::Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Float64}}, num_edges::Int64) + # Finds parameters involved in rate and create a function evaluating teh rate law. + 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. `rates` becomes a vector with 1 value. + if all(length(p_val_dict[P]) == 1 for P in relevant_ps) + return [rate_law_func([p_val_dict[p][1] for p in relevant_ps]...)] + # If at least on parameter the rate depends on have a value varying across all edges, + # we have to compute one rate value for each edge. + else + return [rate_law_func([get_component_value(p_val_dict[p], idxE) for p in relevant_ps]...) + for idxE in 1:num_edges] + end +end + +# Creates a map, taking each species (with transportation) to its transportation rate. +# The species is represented by its index (in species(lrs). +# If the rate is uniform across all edges, the vector will be length 1 (with this value), +# else there will be a separate value for each edge. +# Pair{Int64, Vector{T}}[] is required in case vector is empty (otherwise it becomes Any[], causing type error later). +function make_sidxs_to_transrate_map(vert_ps::Vector{Vector{Float64}}, edge_ps::Vector{Vector{T}}, + lrs::LatticeReactionSystem) where T + transport_rates_speciesmap = compute_all_transport_rates(vert_ps, edge_ps, lrs) + return Pair{Int64, Vector{T}}[ + speciesmap(lrs.rs)[spat_rates[1]] => spat_rates[2] for spat_rates in transport_rates_speciesmap + ] +end + +### Accessing State & Parameter Array Values ### + +# Gets the index in the u array of species s in vertex vert (when their 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). +get_indexes(vert::Int64, num_species::Int64) = ((vert - 1) * num_species + 1):(vert * num_species) + +# For vectors of length 1 or n, we want to get value idx (or the one value, if length is 1). +# This function gets that. Here: +# - values is the vector with the values of the component across all locations +# (where the internal vectors may or may not be of size 1). +# - component_idx is the initial condition species/vertex parameter/edge parameters's index. +# This is predominantly used for parameters, for initial conditions, +# it is only used once (at initialisation) to re-process the input vector. +# - location_idx is the index of the vertex or edge for which we wish to access a initial condition or parameter values. +# The first two function takes the full value vector, and call the function of at the components specific index. +function get_component_value(values::Vector{<:Vector}, component_idx::Int64, + location_idx::Int64) + get_component_value(values[component_idx], location_idx) +end +# Sometimes we have pre-computed, for each component, whether it's vector is length 1 or not. +# This is stored in location_types. +function get_component_value(values::Vector{<:Vector}, component_idx::Int64, + location_idx::Int64, location_types::Vector{Bool}) + get_component_value(values[component_idx], location_idx, location_types[component_idx]) +end +# For a components value (which is a vector of either length 1 or some other length), retrieves its value. +function get_component_value(values::Vector{<:Number}, location_idx::Int64) + get_component_value(values, location_idx, length(values) == 1) +end +# Again, the location type (length of the value vector) may be pre-computed. +function get_component_value(values::Vector{<:Number}, location_idx::Int64, + location_type::Bool) + location_type ? values[1] : values[location_idx] +end + +# Converts a vector of vectors to a 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::Vector{<:Vector}, n) + vcat([get_component_value.(values, comp) for comp in 1:n]...) +end +function expand_component_values(values::Vector{<:Vector}, n, location_types::Vector{Bool}) + vcat([get_component_value.(values, comp, location_types) for comp in 1:n]...) +end + +# Creates a view of the vert_ps vector at a given location. +# Provides a work vector to which the converted vector is written. +function view_vert_ps_vector!(work_vert_ps, vert_ps, comp, enumerated_vert_ps_idx_types) + # Loops through all parameters. + for (idx,loc_type) in enumerated_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]) + end + return work_vert_ps +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) + reshape(expand_component_values(values, n), length(values), n) +end diff --git a/test/miscellaneous_tests/nonlinear_solve.jl b/test/miscellaneous_tests/nonlinear_solve.jl index f0e3a0b806..c2f74c462d 100644 --- a/test/miscellaneous_tests/nonlinear_solve.jl +++ b/test/miscellaneous_tests/nonlinear_solve.jl @@ -26,7 +26,7 @@ let # Solves it using standard algorithm and simulation based algorithm. sol1 = solve(nl_prob; abstol=1e-12, reltol=1e-12).u - sol2 = solve(nl_prob, DynamicSS(Rosenbrock23(); abstol=1e-12, reltol=1e-12); abstol=1e-12, reltol=1e-12).u + sol2 = solve(nl_prob, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12).u # Tests solutions are correct. @test isapprox(sol1[1], p[1] / p[2]; atol=1e-10) @@ -55,7 +55,7 @@ let # Solves it using standard algorithm and simulation based algorithm. sol1 = solve(nl_prob; abstol=1e-12, reltol=1e-12).u - sol2 = solve(nl_prob, DynamicSS(Rosenbrock23(); abstol=1e-12, reltol=1e-12); abstol=1e-12, reltol=1e-12).u + sol2 = solve(nl_prob, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12).u # Computes NonlinearFunction (manually and automatically). nfunc = NonlinearFunction(convert(NonlinearSystem, steady_state_network_2)) @@ -91,7 +91,7 @@ let # Solves it using standard algorithm and simulation based algorithm. sol1 = solve(nl_prob_1; abstol=1e-12, reltol=1e-12) - sol2 = solve(nl_prob_2, DynamicSS(Rosenbrock23(); abstol=1e-12, reltol=1e-12); abstol=1e-12, reltol=1e-12) + sol2 = solve(nl_prob_2, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12) # Checks output using NonlinearFunction. nfunc = NonlinearFunction(convert(NonlinearSystem, steady_state_network_3)) diff --git a/test/model_simulation/solve_steady_state_problems.jl b/test/model_simulation/solve_steady_state_problems.jl deleted file mode 100644 index 55c44b6115..0000000000 --- a/test/model_simulation/solve_steady_state_problems.jl +++ /dev/null @@ -1,78 +0,0 @@ -### Fetch Packages and Reaction Networks ### - -# Fetch packages. -using Catalyst, OrdinaryDiffEq, Random, SteadyStateDiffEq, Test - -# Sets rnd number. -using StableRNGs -rng = StableRNG(12345) - -# Fetch test networks. -include("../test_networks.jl") - -### Compares Solution to Known Steady States ### - -# Simple network. -let - steady_state_network_1 = @reaction_network begin - (k1, k2), ∅ ↔ X1 - (k3, k4), ∅ ↔ 3X2 - (k5, k6), ∅ ↔ X3 + X4 - end - - for factor in [1e-1, 1e0, 1e1] - u0 = factor * rand(rng, length(states(steady_state_network_1))) - u0[4] = u0[3] - p = 0.01 .+ factor * rand(rng, length(parameters(steady_state_network_1))) - prob = SteadyStateProblem(steady_state_network_1, u0, p) - sol = solve(prob, SSRootfind()).u - (minimum(sol[1:1]) > 1e-2) && (@test abs.(sol[1] - p[1] / p[2]) < 1e-8) - (minimum(sol[2:2]) > 1e-2) && - (@test abs.(sol[2]^3 / factorial(3) - p[3] / p[4]) < 1e-8) - (minimum(sol[3:4]) > 1e-2) && (@test abs.(sol[3] * sol[4] - p[5] / p[6]) < 1e-8) - end -end - -# These are disabled due to problem in SteadyStateProblem solution exactness. We do not recommend this method for finding steady states. - -# steady_state_network_2 = @reaction_network begin -# v / 10 + hill(X, v, K, n), ∅ → X -# d, X → ∅ -# end - -# for factor in [1e-1, 1e1, 1e1], repeat in 1:3 -# u0_small = factor * rand(rng, length(get_states(steady_state_network_2))) / 100 -# u0_large = factor * rand(rng, length(get_states(steady_state_network_2))) * 100 -# p = factor * rand(rng, length(get_ps(steady_state_network_2))) -# p[3] = round(p[3]) + 1 -# sol1 = solve(SteadyStateProblem(steady_state_network_2, u0_small, p), SSRootfind()).u[1] -# sol2 = solve(SteadyStateProblem(steady_state_network_2, u0_large, p), SSRootfind()).u[1] -# diff1 = abs(p[1] / 10 + p[1] * (sol1^p[3]) / (sol1^p[3] + p[2]^p[3]) - p[4] * sol1) -# diff2 = abs(p[1] / 10 + p[1] * (sol2^p[3]) / (sol2^p[3] + p[2]^p[3]) - p[4] * sol2) -# @test (diff1 < 1e-8) || (diff2 < 1e-8) -# end - -### For a couple of networks, test that the steady state solution is identical to the long term ODE solution. ### - -# steady_state_test_networks = [ -# reaction_networks_standard[8], -# reaction_networks_standard[10], -# reaction_networks_weird[1], -# ] -# for network in steady_state_test_networks, factor in [1e-1, 1e0, 1e1] -# u0 = factor * rand(rng, length(get_states(network))) -# p = factor * rand(rng, length(get_ps(network))) -# sol_ode = solve(ODEProblem(network, u0, (0.0, 1000000), p), Rosenbrock23()) -# sol_ss = solve(SteadyStateProblem(network, u0, p), SSRootfind()) -# -# @test all(abs.(sol_ode.u[end] .- sol_ss.u) .< 1e-4) -# end - -### No parameter test ### - -# no_param_network = @reaction_network begin (0.6, 3.2), ∅ ↔ X end -# for factor in [1e0, 1e1, 1e2] -# u0 = factor * rand(rng, length(get_states(no_param_network))) -# sol_ss = solve(SteadyStateProblem(no_param_network, u0), SSRootfind(), abstol = 1e-11) -# @test abs.(sol_ss.u[1] - 0.6 / 3.2) < 1e-8 -# end diff --git a/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl b/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl new file mode 100644 index 0000000000..c3c801c603 --- /dev/null +++ b/test/performance_benchmarks/lattice_reaction_systems_ODE_performance.jl @@ -0,0 +1,192 @@ +# Not actually run in CI, but useful for reference of ODE simulation performance across updates. + +### Preparations ### + +# Fetch packages. +using OrdinaryDiffEq +using Random, Statistics, SparseArrays, Test + +# Fetch test networks. +include("../spatial_test_networks.jl") + +### Runtime Checks ### +# Current timings are taken from the SciML CI server. +# Current not used, simply here for reference. +# Useful when attempting to optimise workflow. + +# using BenchmarkTools, Sundials +# runtime_reduction_margin = 10.0 + +# Small grid, small, non-stiff, system. +let + lrs = LatticeReactionSystem(SIR_system, SIR_srs_2, small_2d_grid) + 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) + @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) + + runtime_target = 0.00027 + runtime = minimum((@benchmark solve($oprob, Tsit5())).times) / 1000000000 + println("Small grid, small, non-stiff, system. Runtime: $(runtime), previous standard: $(runtime_target)") + @test runtime < runtime_reduction_margin * runtime_target +end + +# Large grid, small, non-stiff, system. +let + lrs = LatticeReactionSystem(SIR_system, SIR_srs_2, large_2d_grid) + 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) + @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) + + runtime_target = 0.12 + runtime = minimum((@benchmark solve($oprob, Tsit5())).times) / 1000000000 + println("Large grid, small, non-stiff, system. Runtime: $(runtime), previous standard: $(runtime_target)") + @test runtime < runtime_reduction_margin * runtime_target +end + +# Small grid, small, stiff, system. +let + lrs = LatticeReactionSystem(brusselator_system, brusselator_srs_1, small_2d_grid) + 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)) + @test SciMLBase.successful_retcode(solve(oprob, CVODE_BDF(linear_solver=:GMRES))) + + runtime_target = 0.013 + runtime = minimum((@benchmark solve($oprob, CVODE_BDF(linear_solver=:GMRES))).times) / 1000000000 + println("Small grid, small, stiff, system. Runtime: $(runtime), previous standard: $(runtime_target)") + @test runtime < runtime_reduction_margin * runtime_target +end + +# Large grid, small, stiff, system. +let + lrs = LatticeReactionSystem(brusselator_system, brusselator_srs_1, large_2d_grid) + 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)) + @test SciMLBase.successful_retcode(solve(oprob, CVODE_BDF(linear_solver=:GMRES))) + + runtime_target = 11. + runtime = minimum((@benchmark solve($oprob, CVODE_BDF(linear_solver=:GMRES))).times) / 1000000000 + println("Large grid, small, stiff, system. Runtime: $(runtime), previous standard: $(runtime_target)") + @test runtime < runtime_reduction_margin * runtime_target +end + +# Small grid, mid-sized, non-stiff, system. +let + lrs = LatticeReactionSystem(CuH_Amination_system, CuH_Amination_srs_2, + small_2d_grid) + u0 = [ + :CuoAc => 0.005 .+ rand_v_vals(lrs.lattice, 0.005), + :Ligand => 0.005 .+ rand_v_vals(lrs.lattice, 0.005), + :CuoAcLigand => 0.0, + :Silane => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :CuHLigand => 0.0, + :SilaneOAc => 0.0, + :Styrene => 0.16, + :AlkylCuLigand => 0.0, + :Amine_E => 0.39, + :AlkylAmine => 0.0, + :Cu_ELigand => 0.0, + :E_Silane => 0.0, + :Amine => 0.0, + :Decomposition => 0.0, + ] + 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) + @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) + + runtime_target = 0.0012 + runtime = minimum((@benchmark solve($oprob, Tsit5())).times) / 1000000000 + println("Small grid, mid-sized, non-stiff, system. Runtime: $(runtime), previous standard: $(runtime_target)") + @test runtime < runtime_reduction_margin * runtime_target +end + +# Large grid, mid-sized, non-stiff, system. +let + lrs = LatticeReactionSystem(CuH_Amination_system, CuH_Amination_srs_2, + large_2d_grid) + u0 = [ + :CuoAc => 0.005 .+ rand_v_vals(lrs.lattice, 0.005), + :Ligand => 0.005 .+ rand_v_vals(lrs.lattice, 0.005), + :CuoAcLigand => 0.0, + :Silane => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :CuHLigand => 0.0, + :SilaneOAc => 0.0, + :Styrene => 0.16, + :AlkylCuLigand => 0.0, + :Amine_E => 0.39, + :AlkylAmine => 0.0, + :Cu_ELigand => 0.0, + :E_Silane => 0.0, + :Amine => 0.0, + :Decomposition => 0.0, + ] + 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) + @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) + + runtime_target = 0.56 + runtime = minimum((@benchmark solve($oprob, Tsit5())).times) / 1000000000 + println("Large grid, mid-sized, non-stiff, system. Runtime: $(runtime), previous standard: $(runtime_target)") + @test runtime < runtime_reduction_margin * runtime_target +end + +# Small grid, mid-sized, stiff, system. +let + lrs = LatticeReactionSystem(sigmaB_system, sigmaB_srs_2, small_2d_grid) + u0 = [ + :w => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :w2 => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :w2v => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :v => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :w2v2 => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :vP => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :σB => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :w2σB => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :vPp => 0.0, + :phos => 0.4, + ] + pV = sigmaB_p + pE = [:DσB => 0.1, :Dw => 0.1, :Dv => 0.1] + 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 + runtime = minimum((@benchmark solve($oprob, CVODE_BDF(linear_solver=:GMRES))).times) / 1000000000 + println("Small grid, mid-sized, stiff, system. Runtime: $(runtime), previous standard: $(runtime_target)") + @test runtime < runtime_reduction_margin * runtime_target +end + +# Large grid, mid-sized, stiff, system. +let + lrs = LatticeReactionSystem(sigmaB_system, sigmaB_srs_2, large_2d_grid) + u0 = [ + :w => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :w2 => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :w2v => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :v => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :w2v2 => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :vP => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :σB => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :w2σB => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :vPp => 0.0, + :phos => 0.4, + ] + 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). + @test SciMLBase.successful_retcode(solve(oprob, CVODE_BDF(linear_solver=:GMRES))) + + runtime_target = 59. + runtime = minimum((@benchmark solve($oprob, CVODE_BDF(linear_solver=:GMRES))).times) / 1000000000 + println("Large grid, mid-sized, stiff, system. Runtime: $(runtime), previous standard: $(runtime_target)") + @test runtime < runtime_reduction_margin * runtime_target +end diff --git a/test/runtests.jl b/test/runtests.jl index 221489f7a2..0685c2fc81 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -39,8 +39,11 @@ using SafeTestsets @time @safetestset "U0 and Parameters Input Variants" begin include("model_simulation/u0_n_parameter_inputs.jl") end @time @safetestset "SDE System Simulations" begin include("model_simulation/simulate_SDEs.jl") end @time @safetestset "Jump System Simulations" begin include("model_simulation/simulate_jumps.jl") end - @time @safetestset "DiffEq Steady State Solving" begin include("model_simulation/solve_steady_state_problems.jl") end - @time @safetestset "PDE Systems Simulations" begin include("model_simulation/simulate_PDEs.jl") end + + ### Tests Spatial Network Simulations. ### + @time @safetestset "PDE Systems Simulations" begin include("spatial_reaction_systems/simulate_PDEs.jl") end + @time @safetestset "Lattice Reaction Systems" begin include("spatial_reaction_systems/lattice_reaction_systems.jl") end + @time @safetestset "ODE Lattice Systems Simulations" begin include("spatial_reaction_systems/lattice_reaction_systems_ODEs.jl") end ### Tests network visualization. ### @time @safetestset "Latexify" begin include("visualization/latexify.jl") end diff --git a/test/spatial_reaction_systems/lattice_reaction_systems.jl b/test/spatial_reaction_systems/lattice_reaction_systems.jl new file mode 100644 index 0000000000..00ded50494 --- /dev/null +++ b/test/spatial_reaction_systems/lattice_reaction_systems.jl @@ -0,0 +1,276 @@ +### Preparations ### + +# Fetch packages. +using Catalyst, Graphs, Test + +# Pre declares a grid. +grid = Graphs.grid([2, 2]) + +### Tests LatticeReactionSystem Getters Correctness ### +# Test case 1. +let + rs = @reaction_network begin + (p, 1), 0 <--> X + end + tr = @transport_reaction d X + lrs = LatticeReactionSystem(rs, [tr], grid) + + @unpack X, p = rs + d = edge_parameters(lrs)[1] + @test issetequal(species(lrs), [X]) + @test issetequal(spatial_species(lrs), [X]) + @test issetequal(parameters(lrs), [p, d]) + @test issetequal(vertex_parameters(lrs), [p]) + @test issetequal(edge_parameters(lrs), [d]) +end + +# Test case 2. +let + rs = @reaction_network begin + @parameters p1 p2 [edgeparameter=true] + end + @unpack p1, p2 = rs + + @test !isedgeparameter(p1) + @test isedgeparameter(p2) +end + +# Test case 3. +let + rs = @reaction_network begin + @parameters pX pY dX [edgeparameter=true] dY + (pX, 1), 0 <--> X + (pY, 1), 0 <--> Y + end + tr_1 = @transport_reaction dX X + tr_2 = @transport_reaction dY Y + lrs = LatticeReactionSystem(rs, [tr_1, tr_2], grid) + + @unpack X, Y, pX, pY, dX, dY = rs + @test issetequal(species(lrs), [X, Y]) + @test issetequal(spatial_species(lrs), [X, Y]) + @test issetequal(parameters(lrs), [pX, pY, dX, dY]) + @test issetequal(vertex_parameters(lrs), [pX, pY, dY]) + @test issetequal(edge_parameters(lrs), [dX]) +end + +# Test case 4. +let + rs = @reaction_network begin + @parameters dX p + (pX, 1), 0 <--> X + (pY, 1), 0 <--> Y + end + tr_1 = @transport_reaction dX X + lrs = LatticeReactionSystem(rs, [tr_1], grid) + + @unpack dX, p, X, Y, pX, pY = rs + @test issetequal(species(lrs), [X, Y]) + @test issetequal(spatial_species(lrs), [X]) + @test issetequal(parameters(lrs), [dX, p, pX, pY]) + @test issetequal(vertex_parameters(lrs), [dX, p, pX, pY]) + @test issetequal(edge_parameters(lrs), []) +end + +# Test case 5. +let + rs = @reaction_network begin + @species W(t) + @parameters pX pY dX [edgeparameter=true] dY + (pX, 1), 0 <--> X + (pY, 1), 0 <--> Y + (pZ, 1), 0 <--> Z + (pV, 1), 0 <--> V + end + @unpack dX, X, V = rs + @parameters dV dW + @variables t + @species W(t) + tr_1 = TransportReaction(dX, X) + tr_2 = @transport_reaction dY Y + tr_3 = @transport_reaction dZ Z + tr_4 = TransportReaction(dV, V) + tr_5 = TransportReaction(dW, W) + lrs = LatticeReactionSystem(rs, [tr_1, tr_2, tr_3, tr_4, tr_5], grid) + + @unpack pX, pY, pZ, pV, dX, dY, X, Y, Z, V = rs + dZ, dV, dW = edge_parameters(lrs)[2:end] + @test issetequal(species(lrs), [W, X, Y, Z, V]) + @test issetequal(spatial_species(lrs), [X, Y, Z, V, W]) + @test issetequal(parameters(lrs), [pX, pY, dX, dY, pZ, pV, dZ, dV, dW]) + @test issetequal(vertex_parameters(lrs), [pX, pY, dY, pZ, pV]) + @test issetequal(edge_parameters(lrs), [dX, dZ, dV, dW]) +end + +# Test case 6. +let + rs = @reaction_network customname begin + (p, 1), 0 <--> X + end + tr = @transport_reaction d X + lrs = LatticeReactionSystem(rs, [tr], grid) + + @test nameof(lrs) == :customname +end + +### Tests Spatial Reactions Getters Correctness ### + +# Test case 1. +let + tr_1 = @transport_reaction dX X + tr_2 = @transport_reaction dY1*dY2 Y + + # @test ModelingToolkit.getname.(species(tr_1)) == ModelingToolkit.getname.(spatial_species(tr_1)) == [:X] # species(::TransportReaction) currently not supported. + # @test ModelingToolkit.getname.(species(tr_2)) == ModelingToolkit.getname.(spatial_species(tr_2)) == [:Y] + @test ModelingToolkit.getname.(spatial_species(tr_1)) == [:X] + @test ModelingToolkit.getname.(spatial_species(tr_2)) == [:Y] + @test ModelingToolkit.getname.(parameters(tr_1)) == [:dX] + @test ModelingToolkit.getname.(parameters(tr_2)) == [:dY1, :dY2] + + # @test issetequal(species(tr_1), [tr_1.species]) + # @test issetequal(species(tr_2), [tr_2.species]) + @test issetequal(spatial_species(tr_1), [tr_1.species]) + @test issetequal(spatial_species(tr_2), [tr_2.species]) +end + +# Test case 2. +let + rs = @reaction_network begin + @species X(t) Y(t) + @parameters dX dY1 dY2 + end + @unpack X, Y, dX, dY1, dY2 = rs + tr_1 = TransportReaction(dX, X) + tr_2 = TransportReaction(dY1*dY2, Y) + # @test isequal(species(tr_1), [X]) + # @test isequal(species(tr_1), [X]) + @test issetequal(spatial_species(tr_2), [Y]) + @test issetequal(spatial_species(tr_2), [Y]) + @test issetequal(parameters(tr_1), [dX]) + @test issetequal(parameters(tr_2), [dY1, dY2]) +end + +### Tests Spatial Reactions Generation ### + +# Tests TransportReaction with non-trivial rate. +let + rs = @reaction_network begin + @parameters dV dE [edgeparameter=true] + (p,1), 0 <--> X + end + @unpack dV, dE, X = rs + + tr = TransportReaction(dV*dE, X) + @test isequal(tr.rate, dV*dE) +end + +# Tests transport_reactions function for creating TransportReactions. +let + rs = @reaction_network begin + @parameters d + (p,1), 0 <--> X + end + @unpack d, X = rs + trs = TransportReactions([(d, X), (d, X)]) + @test isequal(trs[1], trs[2]) +end + +# Test reactions with constants in rate. +let + @variables t + @species X(t) Y(t) + + tr_1 = TransportReaction(1.5, X) + tr_1_macro = @transport_reaction 1.5 X + @test isequal(tr_1.rate, tr_1_macro.rate) + @test isequal(tr_1.species, tr_1_macro.species) + + tr_2 = TransportReaction(π, Y) + tr_2_macro = @transport_reaction π Y + @test isequal(tr_2.rate, tr_2_macro.rate) + @test isequal(tr_2.species, tr_2_macro.species) +end + +### Test Interpolation ### + +# Does not currently work. The 3 tr_macro_ lines generate errors. +# Test case 1. +let + rs = @reaction_network begin + @species X(t) Y(t) Z(t) + @parameters dX dY1 dY2 dZ + end + @unpack X, Y, Z, dX, dY1, dY2, dZ = rs + rate1 = dX + rate2 = dY1*dY2 + species3 = Z + tr_1 = TransportReaction(dX, X) + tr_2 = TransportReaction(dY1*dY2, Y) + tr_3 = TransportReaction(dZ, Z) + tr_macro_1 = @transport_reaction $dX X + tr_macro_2 = @transport_reaction $(rate2) Y + # tr_macro_3 = @transport_reaction dZ $species3 # Currently does not work, something with meta programming. + + @test isequal(tr_1, tr_macro_1) + @test isequal(tr_2, tr_macro_2) # Unsure why these fails, since for components equality hold: `isequal(tr_1.species, tr_macro_1.species)` and `isequal(tr_1.rate, tr_macro_1.rate)` are both true. + # @test isequal(tr_3, tr_macro_3) +end + +### Tests Error generation ### + +# Test creation of TransportReaction with non-parameters in rate. +# Tests that it works even when rate is highly nested. +let + @variables t + @species X(t) Y(t) + @parameters D1 D2 D3 + @test_throws ErrorException TransportReaction(D1 + D2*(D3 + Y), X) + @test_throws ErrorException TransportReaction(Y, X) +end + +# Network where diffusion species is not declared in non-spatial network. +let + rs = @reaction_network begin + (p, d), 0 <--> X + end + tr = @transport_reaction D Y + @test_throws ErrorException LatticeReactionSystem(rs, [tr], grid) +end + +# Network where the rate depend on a species +let + rs = @reaction_network begin + @species Y(t) + (p, d), 0 <--> X + end + tr = @transport_reaction D*Y X + @test_throws ErrorException LatticeReactionSystem(rs, [tr], grid) +end + +# Network with edge parameter in non-spatial reaction rate. +let + rs = @reaction_network begin + @parameters p [edgeparameter=true] + (p, d), 0 <--> X + end + tr = @transport_reaction D X + @test_throws ErrorException LatticeReactionSystem(rs, [tr], grid) +end + +# Network where metadata has been added in rs (which is not seen in transport reaction). +let + rs = @reaction_network begin + @species X(t) [description="Species with added metadata"] + (p, d), 0 <--> X + end + tr = @transport_reaction D X + @test_throws ErrorException LatticeReactionSystem(rs, [tr], grid) + + rs = @reaction_network begin + @parameters D [description="Parameter with added metadata"] + (p, d), 0 <--> X + end + tr = @transport_reaction D X + @test_throws ErrorException LatticeReactionSystem(rs, [tr], grid) +end + diff --git a/test/spatial_reaction_systems/lattice_reaction_systems_ODEs.jl b/test/spatial_reaction_systems/lattice_reaction_systems_ODEs.jl new file mode 100644 index 0000000000..7355f34ef4 --- /dev/null +++ b/test/spatial_reaction_systems/lattice_reaction_systems_ODEs.jl @@ -0,0 +1,616 @@ +### Preparations ### + +# Fetch packages. +using OrdinaryDiffEq +using Random, Statistics, SparseArrays, Test + +# Fetch test networks. +include("../spatial_test_networks.jl") + +# Sets rnd number. +using StableRNGs +rng = StableRNG(12345) + +### Tests Simulations Don't Error ### +for grid in [small_2d_grid, short_path, small_directed_cycle] + # Non-stiff case + for srs in [Vector{TransportReaction}(), SIR_srs_1, SIR_srs_2] + lrs = LatticeReactionSystem(SIR_system, srs, grid) + u0_1 = [:S => 999.0, :I => 1.0, :R => 0.0] + u0_2 = [:S => 500.0 .+ 500.0 * rand_v_vals(lrs.lattice), :I => 1.0, :R => 0.0] + u0_3 = [ + :S => 950.0, + :I => 50 * rand_v_vals(lrs.lattice), + :R => 50 * rand_v_vals(lrs.lattice), + ] + u0_4 = [ + :S => 500.0 .+ 500.0 * rand_v_vals(lrs.lattice), + :I => 50 * rand_v_vals(lrs.lattice), + :R => 50 * rand_v_vals(lrs.lattice), + ] + u0_5 = make_u0_matrix(u0_3, vertices(lrs.lattice), + map(s -> Symbol(s.f), species(lrs.rs))) + for u0 in [u0_1, u0_2, u0_3, u0_4, u0_5] + p1 = [:α => 0.1 / 1000, :β => 0.01] + p2 = [:α => 0.1 / 1000, :β => 0.02 * rand_v_vals(lrs.lattice)] + p3 = [ + :α => 0.1 / 2000 * rand_v_vals(lrs.lattice), + :β => 0.02 * rand_v_vals(lrs.lattice), + ] + p4 = make_u0_matrix(p1, vertices(lrs.lattice), Symbol.(parameters(lrs.rs))) + for pV in [p1, p2, p3, p4] + pE_1 = map(sp -> sp => 0.01, spatial_param_syms(lrs)) + pE_2 = map(sp -> sp => 0.01, spatial_param_syms(lrs)) + pE_3 = map(sp -> sp => rand_e_vals(lrs.lattice, 0.01), + spatial_param_syms(lrs)) + pE_4 = make_u0_matrix(pE_3, edges(lrs.lattice), spatial_param_syms(lrs)) + for pE in [pE_1, pE_2, pE_3, pE_4] + oprob = ODEProblem(lrs, u0, (0.0, 500.0), (pV, pE)) + @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) + + oprob = ODEProblem(lrs, u0, (0.0, 10.0), (pV, pE); jac = false) + @test SciMLBase.successful_retcode(solve(oprob, Tsit5())) + end + end + end + end + + # Stiff case + for srs in [Vector{TransportReaction}(), brusselator_srs_1, brusselator_srs_2] + lrs = LatticeReactionSystem(brusselator_system, srs, grid) + u0_1 = [:X => 1.0, :Y => 20.0] + u0_2 = [:X => rand_v_vals(lrs.lattice, 10.0), :Y => 2.0] + u0_3 = [:X => rand_v_vals(lrs.lattice, 20), :Y => rand_v_vals(lrs.lattice, 10)] + u0_4 = make_u0_matrix(u0_3, vertices(lrs.lattice), + map(s -> Symbol(s.f), species(lrs.rs))) + for u0 in [u0_1, u0_2, u0_3, u0_4] + p1 = [:A => 1.0, :B => 4.0] + p2 = [:A => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), :B => 4.0] + p3 = [ + :A => 0.5 .+ rand_v_vals(lrs.lattice, 0.5), + :B => 4.0 .+ rand_v_vals(lrs.lattice, 1.0), + ] + p4 = make_u0_matrix(p2, vertices(lrs.lattice), Symbol.(parameters(lrs.rs))) + for pV in [p1, p2, p3, p4] + pE_1 = map(sp -> sp => 0.2, spatial_param_syms(lrs)) + pE_2 = map(sp -> sp => rand(rng), spatial_param_syms(lrs)) + pE_3 = map(sp -> sp => rand_e_vals(lrs.lattice, 0.2), + spatial_param_syms(lrs)) + pE_4 = make_u0_matrix(pE_3, edges(lrs.lattice), spatial_param_syms(lrs)) + for pE in [pE_1, pE_2, pE_3, pE_4] + oprob = ODEProblem(lrs, u0, (0.0, 10.0), (pV, pE)) + @test SciMLBase.successful_retcode(solve(oprob, QNDF())) + + oprob = ODEProblem(lrs, u0, (0.0, 10.0), (pV, pE); sparse = false) + @test SciMLBase.successful_retcode(solve(oprob, QNDF())) + end + end + end + end +end + +### Tests Simulation Correctness ### + +# Checks that non-spatial brusselator simulation is identical to all on an unconnected lattice. +let + lrs = LatticeReactionSystem(brusselator_system, brusselator_srs_1, unconnected_graph) + u0 = [:X => 2.0 + 2.0 * rand(rng), :Y => 10.0 + 10.0 * rand(rng)] + pV = brusselator_p + pE = [:dX => 0.2] + oprob_nonspatial = ODEProblem(brusselator_system, u0, (0.0, 100.0), pV) + oprob_spatial = ODEProblem(lrs, u0, (0.0, 100.0), (pV, pE)) + sol_nonspatial = solve(oprob_nonspatial, QNDF(); abstol = 1e-12, reltol = 1e-12) + sol_spatial = solve(oprob_spatial, QNDF(); abstol = 1e-12, reltol = 1e-12) + + for i in 1:nv(unconnected_graph) + @test all(isapprox.(sol_nonspatial.u[end], + sol_spatial.u[end][((i - 1) * 2 + 1):((i - 1) * 2 + 2)])) + end +end + +# Compares Jacobian and forcing functions of spatial system to analytically computed on. +let + # Creates LatticeReactionNetwork ODEProblem. + rs = @reaction_network begin + pX, 0 --> X + d, X --> 0 + pY*X, 0 --> Y + d, Y --> 0 + end + tr = @transport_reaction D X + lattice = path_graph(3) + lrs = LatticeReactionSystem(rs, [tr], lattice); + + D_vals = [0.2, 0.2, 0.3, 0.3] + u0 = [:X => [1.0, 2.0, 3.0], :Y => 1.0] + ps = [:pX => [2.0, 2.5, 3.0], :pY => 0.5, :d => 0.1, :D => D_vals] + oprob = ODEProblem(lrs, u0, (0.0, 0.0), ps; jac=true, sparse=true) + + # Creates manual f and jac functions. + function f_manual!(du, u, p, t) + X1, Y1, X2, Y2, X3, Y3 = u + pX, d, pY = p + pX1, pX2, pX3 = pX + pY, = pY + d, = d + D1, D2, D3, D4 = D_vals + du[1] = pX1 - d*X1 - D1*X1 + D2*X2 + du[2] = pY*X1 - d*Y1 + du[3] = pX2 - d*X2 + D1*X1 - (D2+D3)*X2 + D4*X3 + du[4] = pY*X2 - d*Y2 + du[5] = pX3 - d*X3 + D3*X2 - D4*X3 + du[6] = pY*X3 - d*Y3 + end + function jac_manual!(J, u, p, t) + X1, Y1, X2, Y2, X3, Y3 = u + pX, d, pY = p + pX1, pX2, pX3 = pX + pY, = pY + d, = d + D1, D2, D3, D4 = D_vals + + J .= 0.0 + + J[1,1] = - d - D1 + J[1,2] = 0 + J[2,1] = pY + J[2,2] = - d + + J[3,3] = - d - D2 - D3 + J[3,4] = 0 + J[4,3] = pY + J[4,4] = - d + + J[5,5] = - d - D4 + J[5,6] = 0 + J[6,5] = pY + J[6,6] = - d + + J[1,3] = D1 + J[3,1] = D2 + J[3,5] = D3 + J[5,3] = D4 + end + + # Sets test input values. + u = rand(rng, 6) + p = [rand(rng, 3), rand(rng, 1), rand(rng, 1)] + + # Tests forcing function. + du1 = fill(0.0, 6) + du2 = fill(0.0, 6) + oprob.f(du1, u, p, 0.0) + f_manual!(du2, u, p, 0.0) + @test du1 ≈ du2 + + # Tests Jacobian. + J1 = deepcopy(oprob.f.jac_prototype) + J2 = deepcopy(oprob.f.jac_prototype) + oprob.f.jac(J1, u, p, 0.0) + jac_manual!(J2, u, p, 0.0) + @test J1 ≈ J2 +end + +# Checks that result becomes homogeneous on a connected lattice. +let + lrs = LatticeReactionSystem(binding_system, binding_srs, undirected_cycle) + u0 = [ + :X => 1.0 .+ rand_v_vals(lrs.lattice), + :Y => 2.0 * rand_v_vals(lrs.lattice), + :XY => 0.5, + ] + oprob = ODEProblem(lrs, u0, (0.0, 1000.0), binding_p; tstops = 0.1:0.1:1000.0) + ss = solve(oprob, Tsit5()).u[end] + + @test all(isapprox.(ss[1:3:end], ss[1])) + @test all(isapprox.(ss[2:3:end], ss[2])) + @test all(isapprox.(ss[3:3:end], ss[3])) +end + +# Checks that various combinations of jac and sparse gives the same result. +let + lrs = LatticeReactionSystem(brusselator_system, brusselator_srs_1, small_2d_grid) + 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, 50.0), (pV, pE); jac = false, sparse = false) + oprob_sparse = ODEProblem(lrs, u0, (0.0, 50.0), (pV, pE); jac = false, sparse = true) + oprob_jac = ODEProblem(lrs, u0, (0.0, 50.0), (pV, pE); jac = true, sparse = false) + oprob_sparse_jac = ODEProblem(lrs, u0, (0.0, 50.0), (pV, pE); jac = true, sparse = true) + + ss = solve(oprob, Rosenbrock23(); abstol = 1e-10, reltol = 1e-10).u[end] + @test all(isapprox.(ss, + solve(oprob_sparse, Rosenbrock23(); abstol = 1e-10, reltol = 1e-10).u[end]; + rtol = 0.0001)) + @test all(isapprox.(ss, + solve(oprob_jac, Rosenbrock23(); abstol = 1e-10, reltol = 1e-10).u[end]; + rtol = 0.0001)) + @test all(isapprox.(ss, + solve(oprob_sparse_jac, Rosenbrock23(); abstol = 1e-10, reltol = 1e-10).u[end]; + rtol = 0.0001)) +end + +# Checks that, when non directed graphs are provided, the parameters are re-ordered correctly. +let + # Create the same lattice (one as digraph, one not). Algorithm depends on Graphs.jl reordering edges, hence the jumbled order. + lattice_1 = SimpleGraph(5) + lattice_2 = SimpleDiGraph(5) + + add_edge!(lattice_1, 5, 2) + add_edge!(lattice_1, 1, 4) + add_edge!(lattice_1, 1, 3) + add_edge!(lattice_1, 4, 3) + add_edge!(lattice_1, 4, 5) + + add_edge!(lattice_2, 4, 1) + add_edge!(lattice_2, 3, 4) + add_edge!(lattice_2, 5, 4) + add_edge!(lattice_2, 5, 2) + add_edge!(lattice_2, 4, 3) + add_edge!(lattice_2, 4, 5) + add_edge!(lattice_2, 3, 1) + add_edge!(lattice_2, 2, 5) + add_edge!(lattice_2, 1, 4) + add_edge!(lattice_2, 1, 3) + + lrs_1 = LatticeReactionSystem(SIR_system, SIR_srs_2, lattice_1) + lrs_2 = LatticeReactionSystem(SIR_system, SIR_srs_2, lattice_2) + + u0 = [:S => 990.0, :I => 20.0 * rand_v_vals(lrs_1.lattice), :R => 0.0] + pV = [:α => 0.1 / 1000, :β => 0.01] + + pE_1 = [:dS => [1.3, 1.4, 2.5, 3.4, 4.5], :dI => 0.01, :dR => 0.02] + pE_2 = [:dS => [1.3, 1.4, 2.5, 1.3, 3.4, 1.4, 3.4, 4.5, 2.5, 4.5], :dI => 0.01, :dR => 0.02] + ss_1 = solve(ODEProblem(lrs_1, u0, (0.0, 500.0), (pV, pE_1)), Tsit5()).u[end] + ss_2 = solve(ODEProblem(lrs_2, u0, (0.0, 500.0), (pV, pE_2)), Tsit5()).u[end] + @test all(isapprox.(ss_1, ss_2)) +end + +### Test Transport Reaction Types ### + +# Compares where spatial reactions are created with/without the macro. +let + @parameters dS dI + @unpack S, I = SIR_system + tr_1 = TransportReaction(dS, S) + tr_2 = TransportReaction(dI, I) + tr_macros_1 = @transport_reaction dS S + tr_macros_2 = @transport_reaction dI I + + lrs_1 = LatticeReactionSystem(SIR_system, [tr_1, tr_2], small_2d_grid) + lrs_2 = LatticeReactionSystem(SIR_system, [tr_macros_1, tr_macros_2], small_2d_grid) + u0 = [:S => 990.0, :I => 20.0 * rand_v_vals(lrs_1.lattice), :R => 0.0] + pV = [:α => 0.1 / 1000, :β => 0.01] + pE = [:dS => 0.01, :dI => 0.01] + ss_1 = solve(ODEProblem(lrs_1, u0, (0.0, 500.0), (pV, pE)), Tsit5()).u[end] + ss_2 = solve(ODEProblem(lrs_2, u0, (0.0, 500.0), (pV, pE)), Tsit5()).u[end] + @test all(isapprox.(ss_1, ss_2)) +end + +# Tries non-trivial diffusion rates. +let + SIR_tr_S_alt = @transport_reaction dS1+dS2 S + SIR_tr_I_alt = @transport_reaction dI1*dI2 I + SIR_tr_R_alt = @transport_reaction log(dR1)+dR2 R + SIR_srs_2_alt = [SIR_tr_S_alt, SIR_tr_I_alt, SIR_tr_R_alt] + lrs_1 = LatticeReactionSystem(SIR_system, SIR_srs_2, small_2d_grid) + lrs_2 = LatticeReactionSystem(SIR_system, SIR_srs_2_alt, small_2d_grid) + + u0 = [:S => 990.0, :I => 20.0 * rand_v_vals(lrs_1.lattice), :R => 0.0] + pV = [:α => 0.1 / 1000, :β => 0.01] + pE_1 = [:dS => 0.01, :dI => 0.01, :dR => 0.01] + pE_2 = [:dS1 => 0.005, :dS1 => 0.005, :dI1 => 2, :dI2 => 0.005, :dR1 => 1.010050167084168, :dR2 => 1.0755285551056204e-16] + + ss_1 = solve(ODEProblem(lrs_1, u0, (0.0, 500.0), (pV, pE_1)), Tsit5()).u[end] + ss_2 = solve(ODEProblem(lrs_2, u0, (0.0, 500.0), (pV, pE_2)), Tsit5()).u[end] + @test all(isapprox.(ss_1, ss_2)) +end + +# Tries various ways of creating TransportReactions. +let + CuH_Amination_system_alt_1 = @reaction_network begin + @species Newspecies1(t) Newspecies2(t) + @parameters dCuoAc [edgeparameter=true] dLigand dSilane dStyrene dCu_ELigand + 10.0^kp1, CuoAc + Ligand --> CuoAcLigand + 10.0^kp2, CuoAcLigand + Silane --> CuHLigand + SilaneOAc + 10.0^k1, CuHLigand + Styrene --> AlkylCuLigand + 10.0^k_1, AlkylCuLigand --> CuHLigand + Styrene + 10.0^k2, AlkylCuLigand + Amine_E --> AlkylAmine + Cu_ELigand + 10.0^k_2, AlkylAmine + Cu_ELigand --> AlkylCuLigand + Amine_E + 10.0^k3, Cu_ELigand + Silane --> CuHLigand + E_Silane + 10.0^kam, CuHLigand + Amine_E --> Amine + Cu_ELigand + 10.0^kdc, CuHLigand + CuHLigand --> Decomposition + end + @unpack dLigand, dSilane, Silane = CuH_Amination_system_alt_1 + @parameters dAmine_E dNewspecies1 + @variables t + @species Ligand(t) Amine_E(t) Newspecies1(t) + tr_alt_1_1 = TransportReaction(dLigand, Ligand) + tr_alt_1_2 = TransportReaction(dSilane, Silane) + tr_alt_1_3 = TransportReaction(dAmine_E, Amine_E) + tr_alt_1_4 = TransportReaction(dNewspecies1, Newspecies1) + tr_alt_1_5 = @transport_reaction dDecomposition Decomposition + tr_alt_1_6 = @transport_reaction dCu_ELigand Cu_ELigand + tr_alt_1_7 = @transport_reaction dNewspecies2 Newspecies2 + CuH_Amination_srs_alt_1 = [tr_alt_1_1, tr_alt_1_2, tr_alt_1_3, tr_alt_1_4, tr_alt_1_5, tr_alt_1_6, tr_alt_1_7] + lrs_1 = LatticeReactionSystem(CuH_Amination_system_alt_1, CuH_Amination_srs_alt_1, small_2d_grid) + + CuH_Amination_system_alt_2 = @reaction_network begin + @species Newspecies1(t) Newspecies2(t) + @parameters dCuoAc [edgeparameter=true] dLigand dSilane dStyrene dCu_ELigand + 10.0^kp1, CuoAc + Ligand --> CuoAcLigand + 10.0^kp2, CuoAcLigand + Silane --> CuHLigand + SilaneOAc + 10.0^k1, CuHLigand + Styrene --> AlkylCuLigand + 10.0^k_1, AlkylCuLigand --> CuHLigand + Styrene + 10.0^k2, AlkylCuLigand + Amine_E --> AlkylAmine + Cu_ELigand + 10.0^k_2, AlkylAmine + Cu_ELigand --> AlkylCuLigand + Amine_E + 10.0^k3, Cu_ELigand + Silane --> CuHLigand + E_Silane + 10.0^kam, CuHLigand + Amine_E --> Amine + Cu_ELigand + 10.0^kdc, CuHLigand + CuHLigand --> Decomposition + end + @unpack Decomposition, dCu_ELigand, Cu_ELigand = CuH_Amination_system_alt_2 + @parameters dNewspecies2 dDecomposition + @variables t + @species Newspecies2(t) + tr_alt_2_1 = @transport_reaction dLigand Ligand + tr_alt_2_2 = @transport_reaction dSilane Silane + tr_alt_2_3 = @transport_reaction dAmine_E Amine_E + tr_alt_2_4 = @transport_reaction dNewspecies1 Newspecies1 + tr_alt_2_5 = TransportReaction(dDecomposition, Decomposition) + tr_alt_2_6 = TransportReaction(dCu_ELigand, Cu_ELigand) + tr_alt_2_7 = TransportReaction(dNewspecies2, Newspecies2) + CuH_Amination_srs_alt_2 = [tr_alt_2_1, tr_alt_2_2, tr_alt_2_3, tr_alt_2_4, tr_alt_2_5, tr_alt_2_6, tr_alt_2_7] + lrs_2 = LatticeReactionSystem(CuH_Amination_system_alt_2, CuH_Amination_srs_alt_2, small_2d_grid) + + u0 = [CuH_Amination_u0; :Newspecies1 => 0.1; :Newspecies2 => 0.1] + pV = [CuH_Amination_p; :dLigand => 0.01; :dSilane => 0.01; :dCu_ELigand => 0.009; :dStyrene => -10000.0] + pE = [:dAmine_E => 0.011, :dNewspecies1 => 0.013, :dDecomposition => 0.015, :dNewspecies2 => 0.016, :dCuoAc => -10000.0] + + ss_1 = solve(ODEProblem(lrs_1, u0, (0.0, 500.0), (pV, pE)), Tsit5()).u[end] + ss_2 = solve(ODEProblem(lrs_2, u0, (0.0, 500.0), (pV, pE)), Tsit5()).u[end] + @test all(isequal.(ss_1, ss_2)) +end + +### Tests Special Cases ### + +# Create network with various combinations of graph/di-graph and parameters. +let + lrs_digraph = LatticeReactionSystem(SIR_system, SIR_srs_2, complete_digraph(3)) + lrs_graph = LatticeReactionSystem(SIR_system, SIR_srs_2, complete_graph(3)) + u0 = [:S => 990.0, :I => 20.0 * rand_v_vals(lrs_digraph.lattice), :R => 0.0] + pV = SIR_p + pE_digraph_1 = [:dS => [0.10, 0.12, 0.10, 0.14, 0.12, 0.14], :dI => 0.01, :dR => 0.01] + pE_digraph_2 = [[0.10, 0.12, 0.10, 0.14, 0.12, 0.14], 0.01, 0.01] + pE_digraph_3 = [0.10 0.12 0.10 0.14 0.12 0.14; 0.01 0.01 0.01 0.01 0.01 0.01; 0.01 0.01 0.01 0.01 0.01 0.01] + pE_graph_1 = [:dS => [0.10, 0.12, 0.14], :dI => 0.01, :dR => 0.01] + pE_graph_2 = [[0.10, 0.12, 0.14], 0.01, 0.01] + pE_graph_3 = [0.10 0.12 0.14; 0.01 0.01 0.01; 0.01 0.01 0.01] + oprob_digraph_1 = ODEProblem(lrs_digraph, u0, (0.0, 500.0), (pV, pE_digraph_1)) + oprob_digraph_2 = ODEProblem(lrs_digraph, u0, (0.0, 500.0), (pV, pE_digraph_2)) + oprob_digraph_3 = ODEProblem(lrs_digraph, u0, (0.0, 500.0), (pV, pE_digraph_3)) + oprob_graph_11 = ODEProblem(lrs_graph, u0, (0.0, 500.0), (pV, pE_digraph_1)) + oprob_graph_12 = ODEProblem(lrs_graph, u0, (0.0, 500.0), (pV, pE_graph_1)) + oprob_graph_21 = ODEProblem(lrs_graph, u0, (0.0, 500.0), (pV, pE_digraph_2)) + oprob_graph_22 = ODEProblem(lrs_graph, u0, (0.0, 500.0), (pV, pE_graph_2)) + oprob_graph_31 = ODEProblem(lrs_graph, u0, (0.0, 500.0), (pV, pE_digraph_3)) + oprob_graph_32 = ODEProblem(lrs_graph, u0, (0.0, 500.0), (pV, pE_graph_3)) + sim_end_digraph_1 = solve(oprob_digraph_1, Tsit5()).u[end] + sim_end_digraph_2 = solve(oprob_digraph_2, Tsit5()).u[end] + sim_end_digraph_3 = solve(oprob_digraph_3, Tsit5()).u[end] + sim_end_graph_11 = solve(oprob_graph_11, Tsit5()).u[end] + sim_end_graph_12 = solve(oprob_graph_12, Tsit5()).u[end] + sim_end_graph_21 = solve(oprob_graph_21, Tsit5()).u[end] + sim_end_graph_22 = solve(oprob_graph_22, Tsit5()).u[end] + sim_end_graph_31 = solve(oprob_graph_31, Tsit5()).u[end] + sim_end_graph_32 = solve(oprob_graph_32, Tsit5()).u[end] + + @test all(sim_end_digraph_1 .== sim_end_digraph_2 .== sim_end_digraph_3 .== + sim_end_graph_11 .== sim_end_graph_12 .== sim_end_graph_21 .== + sim_end_graph_22 .== sim_end_graph_31 .== sim_end_graph_32) +end + +# Creates networks where some species or parameters have no effect on the system. +let + binding_system_alt = @reaction_network begin + @species X(t) Y(t) XY(t) Z(t) V(t) W(t) + @parameters k1 k2 dX [edgeparameter = true] dXY [edgeparameter = true] dZ [edgeparameter = true] dV [edgeparameter = true] p1 p2 + (k1, k2), X + Y <--> XY + end + @unpack dX, dXY, dZ, dV, X, XY, Z, V = binding_system_alt + binding_srs_alt = [ + TransportReaction(dX, X), + TransportReaction(dXY, XY), + TransportReaction(dZ, Z), + TransportReaction(dV, V), + ] + lrs_alt = LatticeReactionSystem(binding_system_alt, binding_srs_alt, small_2d_grid) + u0_alt = [ + :X => 1.0, + :Y => 2.0 * rand_v_vals(lrs_alt.lattice), + :XY => 0.5, + :Z => 2.0 * rand_v_vals(lrs_alt.lattice), + :V => 0.5, + :W => 1.0, + ] + p_alt = [ + :k1 => 2.0, + :k2 => 0.1 .+ rand_v_vals(lrs_alt.lattice), + :dX => 1.0 .+ rand_e_vals(lrs_alt.lattice), + :dXY => 3.0, + :dZ => rand_e_vals(lrs_alt.lattice), + :dV => 0.2, + :p1 => 1.0, + :p2 => rand_v_vals(lrs_alt.lattice), + ] + oprob_alt = ODEProblem(lrs_alt, u0_alt, (0.0, 10.0), p_alt) + ss_alt = solve(oprob_alt, Tsit5()).u[end] + + binding_srs_main = [TransportReaction(dX, X), TransportReaction(dXY, XY)] + lrs = LatticeReactionSystem(binding_system, binding_srs_main, small_2d_grid) + u0 = u0_alt[1:3] + p = p_alt[1:4] + oprob = ODEProblem(lrs, u0, (0.0, 10.0), p) + ss = solve(oprob, Tsit5()).u[end] + + for i in 1:25 + @test isapprox(ss_alt[((i - 1) * 6 + 1):((i - 1) * 6 + 3)], + ss[((i - 1) * 3 + 1):((i - 1) * 3 + 3)]) < 1e-3 + end +end + +# Provides initial conditions and parameters in various different ways. +let + lrs = LatticeReactionSystem(SIR_system, SIR_srs_2, very_small_2d_grid) + u0_1 = [:S => 990.0, :I => [1.0, 3.0, 2.0, 5.0], :R => 0.0] + u0_2 = [990.0, [1.0, 3.0, 2.0, 5.0], 0.0] + u0_3 = [990.0 990.0 990.0 990.0; 1.0 3.0 2.0 5.0; 0.0 0.0 0.0 0.0] + pV_1 = [:α => 0.1 / 1000, :β => [0.01, 0.02, 0.01, 0.03]] + pV_2 = [0.1 / 1000, [0.01, 0.02, 0.01, 0.03]] + pV_3 = [0.1/1000 0.1/1000 0.1/1000 0.1/1000; 0.01 0.02 0.01 0.03] + pE_1 = [:dS => [0.01, 0.02, 0.03, 0.04], :dI => 0.01, :dR => 0.01] + pE_2 = [[0.01, 0.02, 0.03, 0.04], :0.01, 0.01] + pE_3 = [0.01 0.02 0.03 0.04; 0.01 0.01 0.01 0.01; 0.01 0.01 0.01 0.01] + + p1 = [ + :α => 0.1 / 1000, + :β => [0.01, 0.02, 0.01, 0.03], + :dS => [0.01, 0.02, 0.03, 0.04], + :dI => 0.01, + :dR => 0.01, + ] + ss_1_1 = solve(ODEProblem(lrs, u0_1, (0.0, 1.0), p1), Tsit5()).u[end] + for u0 in [u0_1, u0_2, u0_3], pV in [pV_1, pV_2, pV_3], pE in [pE_1, pE_2, pE_3] + ss = solve(ODEProblem(lrs, u0, (0.0, 1.0), (pV, pE)), Tsit5()).u[end] + @test all(isequal.(ss, ss_1_1)) + end +end + +# Confirms parameters can be provided in [pV; pE] and (pV, pE) form. +let + lrs = LatticeReactionSystem(SIR_system, SIR_srs_2, small_2d_grid) + u0 = [:S => 990.0, :I => 20.0 * rand_v_vals(lrs.lattice), :R => 0.0] + p1 = ([:α => 0.1 / 1000, :β => 0.01], [:dS => 0.01, :dI => 0.01, :dR => 0.01]) + p2 = [:α => 0.1 / 1000, :β => 0.01, :dS => 0.01, :dI => 0.01, :dR => 0.01] + oprob1 = ODEProblem(lrs, u0, (0.0, 500.0), p1; jac = false) + oprob2 = ODEProblem(lrs, u0, (0.0, 500.0), p2; jac = false) + + @test all(isapprox.(solve(oprob1, Tsit5()).u[end], solve(oprob2, Tsit5()).u[end])) +end + +### Compare to Hand-written Functions ### + +# Compares the brusselator for a line of cells. +let + function spatial_brusselator_f(du, u, p, t) + # Non-spatial + for i in 1:2:(length(u) - 1) + du[i] = p[1] + 0.5 * (u[i]^2) * u[i + 1] - u[i] - p[2] * u[i] + du[i + 1] = p[2] * u[i] - 0.5 * (u[i]^2) * u[i + 1] + end + + # Spatial + du[1] += p[3] * (u[3] - u[1]) + du[end - 1] += p[3] * (u[end - 3] - u[end - 1]) + for i in 3:2:(length(u) - 3) + du[i] += p[3] * (u[i - 2] + u[i + 2] - 2u[i]) + end + end + function spatial_brusselator_jac(J, u, p, t) + J .= 0 + # Non-spatial + for i in 1:2:(length(u) - 1) + J[i, i] = u[i] * u[i + 1] - 1 - p[2] + J[i, i + 1] = 0.5 * (u[i]^2) + J[i + 1, i] = p[2] - u[i] * u[i + 1] + J[i + 1, i + 1] = -0.5 * (u[i]^2) + end + + # Spatial + J[1, 1] -= p[3] + J[1, 3] += p[3] + J[end - 1, end - 1] -= p[3] + J[end - 1, end - 3] += p[3] + for i in 3:2:(length(u) - 3) + J[i, i] -= 2 * p[3] + J[i, i - 2] += p[3] + J[i, i + 2] += p[3] + end + end + function spatial_brusselator_jac_sparse(J, u, p, t) + # Spatial + J.nzval .= 0.0 + J.nzval[7:6:(end - 9)] .= -2p[3] + J.nzval[1] = -p[3] + J.nzval[end - 3] = -p[3] + J.nzval[3:3:(end - 4)] .= p[3] + + # Non-spatial + for i in 1:1:Int64(lenth(u) / 2 - 1) + j = 6(i - 1) + 1 + J.nzval[j] = u[i] * u[i + 1] - 1 - p[2] + J.nzval[j + 1] = 0.5 * (u[i]^2) + J.nzval[j + 3] = p[2] - u[i] * u[i + 1] + J.nzval[j + 4] = -0.5 * (u[i]^2) + end + J.nzval[end - 3] = u[end - 1] * u[end] - 1 - p[end - 1] + J.nzval[end - 2] = 0.5 * (u[end - 1]^2) + J.nzval[end - 1] = p[2] - u[end - 1] * u[end] + J.nzval[end] = -0.5 * (u[end - 1]^2) + end + function make_jac_prototype(u0) + jac_prototype_pre = zeros(length(u0), length(u0)) + for i in 1:2:(length(u0) - 1) + jac_prototype_pre[i, i] = 1 + jac_prototype_pre[i + 1, i] = 1 + jac_prototype_pre[i, i + 1] = 1 + jac_prototype_pre[i + 1, i + 1] = 1 + end + for i in 3:2:(length(u0) - 1) + jac_prototype_pre[i - 2, i] = 1 + jac_prototype_pre[i, i - 2] = 1 + end + return sparse(jac_prototype_pre) + end + + u0 = 2 * rand(rng, 10000) + p = [1.0, 4.0, 0.1] + tspan = (0.0, 100.0) + + ofun_hw_dense = ODEFunction(spatial_brusselator_f; jac = spatial_brusselator_jac) + ofun_hw_sparse = ODEFunction(spatial_brusselator_f; jac = spatial_brusselator_jac, + jac_prototype = make_jac_prototype(u0)) + + lrs = LatticeReactionSystem(brusselator_system, brusselator_srs_1, + path_graph(Int64(length(u0) / 2))) + u0V = [:X => u0[1:2:(end - 1)], :Y => u0[2:2:end]] + pV = [:A => p[1], :B => p[2]] + pE = [:dX => p[3]] + ofun_aut_dense = ODEProblem(lrs, u0V, tspan, (pV, pE); jac = true, sparse = false).f + ofun_aut_sparse = ODEProblem(lrs, u0V, tspan, (pV, pE); jac = true, sparse = true).f + + du_hw_dense = deepcopy(u0) + du_hw_sparse = deepcopy(u0) + 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) + + @test isapprox(du_hw_dense, du_aut_dense) + @test isapprox(du_hw_sparse, du_aut_sparse) + + J_hw_dense = deepcopy(zeros(length(u0), length(u0))) + J_hw_sparse = deepcopy(make_jac_prototype(u0)) + J_aut_dense = deepcopy(zeros(length(u0), length(u0))) + J_aut_sparse = deepcopy(make_jac_prototype(u0)) + + ofun_hw_dense.jac(J_hw_dense, u0, p, 0.0) + ofun_hw_sparse.jac(J_hw_sparse, u0, p, 0.0) + ofun_aut_dense.jac(J_aut_dense, u0, p, 0.0) + ofun_aut_sparse.jac(J_aut_sparse, u0, p, 0.0) + + @test isapprox(J_hw_dense, J_aut_dense) + @test isapprox(J_hw_sparse, J_aut_sparse) +end diff --git a/test/model_simulation/simulate_PDEs.jl b/test/spatial_reaction_systems/simulate_PDEs.jl similarity index 100% rename from test/model_simulation/simulate_PDEs.jl rename to test/spatial_reaction_systems/simulate_PDEs.jl diff --git a/test/spatial_test_networks.jl b/test/spatial_test_networks.jl new file mode 100644 index 0000000000..f2f9472e87 --- /dev/null +++ b/test/spatial_test_networks.jl @@ -0,0 +1,193 @@ +### Fetch packages ### +using Catalyst, Graphs + +# Sets rnd number. +using StableRNGs +rng = StableRNG(12345) + +### Helper Functions ### + +# Generates randomised initial condition or parameter values. +rand_v_vals(grid) = rand(rng, nv(grid)) +rand_v_vals(grid, x::Number) = rand_v_vals(grid) * x +rand_e_vals(grid) = rand(rng, ne(grid)) +rand_e_vals(grid, x::Number) = rand_e_vals(grid) * x +function make_u0_matrix(value_map, vals, symbols) + (length(symbols) == 0) && (return zeros(0, length(vals))) + d = Dict(value_map) + return [(d[s] isa Vector) ? d[s][v] : d[s] for s in symbols, v in 1:length(vals)] +end + +# Gets a symbol list of spatial parameters. +function spatial_param_syms(lrs::LatticeReactionSystem) + ModelingToolkit.getname.(edge_parameters(lrs)) +end + +# Converts to integer value (for JumpProcess simulations). +function make_values_int(values::Vector{<:Pair}) + [val[1] => round.(Int64, val[2]) for val in values] +end +make_values_int(values::Matrix{<:Number}) = round.(Int64, values) +make_values_int(values::Vector{<:Number}) = round.(Int64, values) +make_values_int(values::Vector{Vector}) = [round.(Int64, vals) for vals in values] + +### Declares Models ### + +# Small non-stiff system. +SIR_system = @reaction_network begin + α, S + I --> 2I + β, I --> R +end +SIR_p = [:α => 0.1 / 1000, :β => 0.01] +SIR_u0 = [:S => 999.0, :I => 1.0, :R => 0.0] + +SIR_tr_S = @transport_reaction dS S +SIR_tr_I = @transport_reaction dI I +SIR_tr_R = @transport_reaction dR R +SIR_srs_1 = [SIR_tr_S] +SIR_srs_2 = [SIR_tr_S, SIR_tr_I, SIR_tr_R] + +# Small non-stiff system. +binding_system = @reaction_network begin (k1, k2), X + Y <--> XY end +binding_tr_X = @transport_reaction dX X +binding_tr_Y = @transport_reaction dY Y +binding_tr_XY = @transport_reaction dXY XY +binding_srs = [binding_tr_X, binding_tr_Y, binding_tr_XY] +binding_u0 = [:X => 1.0, :Y => 2.0, :XY => 0.5] +binding_p = [:k1 => 2.0, :k2 => 0.1, :dX => 3.0, :dY => 5.0, :dXY => 2.0] + +# Mid-sized non-stiff system. +CuH_Amination_system = @reaction_network begin + 10.0^kp1, CuoAc + Ligand --> CuoAcLigand + 10.0^kp2, CuoAcLigand + Silane --> CuHLigand + SilaneOAc + 10.0^k1, CuHLigand + Styrene --> AlkylCuLigand + 10.0^k_1, AlkylCuLigand --> CuHLigand + Styrene + 10.0^k2, AlkylCuLigand + Amine_E --> AlkylAmine + Cu_ELigand + 10.0^k_2, AlkylAmine + Cu_ELigand --> AlkylCuLigand + Amine_E + 10.0^k3, Cu_ELigand + Silane --> CuHLigand + E_Silane + 10.0^kam, CuHLigand + Amine_E --> Amine + Cu_ELigand + 10.0^kdc, CuHLigand + CuHLigand --> Decomposition +end +CuH_Amination_p = [ + :kp1 => 1.2, + :kp2 => -0.72, + :k1 => 0.57, + :k_1 => -3.5, + :k2 => -0.35, + :k_2 => -0.77, + :k3 => -0.025, + :kam => -2.6, + :kdc => -3.0, +] +CuH_Amination_u0 = [ + :CuoAc => 0.0065, + :Ligand => 0.0072, + :CuoAcLigand => 0.0, + :Silane => 0.65, + :CuHLigand => 0.0, + :SilaneOAc => 0.0, + :Styrene => 0.16, + :AlkylCuLigand => 0.0, + :Amine_E => 0.39, + :AlkylAmine => 0.0, + :Cu_ELigand => 0.0, + :E_Silane => 0.0, + :Amine => 0.0, + :Decomposition => 0.0, +] + +CuH_Amination_tr_1 = @transport_reaction D1 CuoAc +CuH_Amination_tr_2 = @transport_reaction D2 Silane +CuH_Amination_tr_3 = @transport_reaction D3 Cu_ELigand +CuH_Amination_tr_4 = @transport_reaction D4 Amine +CuH_Amination_tr_5 = @transport_reaction D5 CuHLigand +CuH_Amination_srs_1 = [CuH_Amination_tr_1] +CuH_Amination_srs_2 = [ + CuH_Amination_tr_1, + CuH_Amination_tr_2, + CuH_Amination_tr_3, + CuH_Amination_tr_4, + CuH_Amination_tr_5, +] + +# Small stiff system. +brusselator_system = @reaction_network begin + A, ∅ → X + 1, 2X + Y → 3X + B, X → Y + 1, X → ∅ +end +brusselator_p = [:A => 1.0, :B => 4.0] + +brusselator_tr_x = @transport_reaction dX X +brusselator_tr_y = @transport_reaction dY Y +brusselator_srs_1 = [brusselator_tr_x] +brusselator_srs_2 = [brusselator_tr_x, brusselator_tr_y] + +# Mid-sized stiff system. +# Unsure about stiffness, but non-spatial version oscillates for this parameter set. +sigmaB_system = @reaction_network begin + kDeg, (w, w2, w2v, v, w2v2, vP, σB, w2σB) ⟶ ∅ + kDeg, vPp ⟶ phos + (kBw, kDw), 2w ⟷ w2 + (kB1, kD1), w2 + v ⟷ w2v + (kB2, kD2), w2v + v ⟷ w2v2 + kK1, w2v ⟶ w2 + vP + kK2, w2v2 ⟶ w2v + vP + (kB3, kD3), w2 + σB ⟷ w2σB + (kB4, kD4), w2σB + v ⟷ w2v + σB + (kB5, kD5), vP + phos ⟷ vPp + kP, vPp ⟶ v + phos + v0 * ((1 + F * σB) / (K + σB)), ∅ ⟶ σB + λW * v0 * ((1 + F * σB) / (K + σB)), ∅ ⟶ w + λV * v0 * ((1 + F * σB) / (K + σB)), ∅ ⟶ v +end +sigmaB_p = [:kBw => 3600, :kDw => 18, :kB1 => 3600, :kB2 => 3600, :kB3 => 3600, + :kB4 => 1800, :kB5 => 3600, + :kD1 => 18, :kD2 => 18, :kD3 => 18, :kD4 => 1800, :kD5 => 18, :kK1 => 36, :kK2 => 6, + :kP => 180, :kDeg => 0.7, + :v0 => 0.4, :F => 30, :K => 0.2, :λW => 4, :λV => 4.5] +sigmaB_u0 = [ + :w => 1.0, + :w2 => 1.0, + :w2v => 1.0, + :v => 1.0, + :w2v2 => 1.0, + :vP => 1.0, + :σB => 1.0, + :w2σB => 1.0, + :vPp => 0.0, + :phos => 0.4, +] + +sigmaB_tr_σB = @transport_reaction DσB σB +sigmaB_tr_w = @transport_reaction Dw w +sigmaB_tr_v = @transport_reaction Dv v +sigmaB_srs_1 = [sigmaB_tr_σB] +sigmaB_srs_2 = [sigmaB_tr_σB, sigmaB_tr_w, sigmaB_tr_v] + +### Declares Lattices ### + +# Grids. +very_small_2d_grid = Graphs.grid([2, 2]) +small_2d_grid = Graphs.grid([5, 5]) +medium_2d_grid = Graphs.grid([20, 20]) +large_2d_grid = Graphs.grid([100, 100]) + +small_3d_grid = Graphs.grid([5, 5, 5]) +medium_3d_grid = Graphs.grid([20, 20, 20]) +large_3d_grid = Graphs.grid([100, 100, 100]) + +# Paths. +short_path = path_graph(100) +long_path = path_graph(1000) + +# Unconnected graphs. +unconnected_graph = SimpleGraph(10) + +# Undirected cycle. +undirected_cycle = cycle_graph(49) + +# Directed cycle. +small_directed_cycle = cycle_graph(100) +large_directed_cycle = cycle_graph(1000) \ No newline at end of file