diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 1757e0241d..72cb17c229 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -8,6 +8,7 @@ using SparseArrays, DiffEqBase, Reexport using LaTeXStrings, Latexify, Requires using JumpProcesses: JumpProcesses, JumpProblem, MassActionJump, ConstantRateJump, VariableRateJump +using SplitApplyCombine # ModelingToolkit imports and convenience functions we use using ModelingToolkit @@ -30,7 +31,9 @@ import ModelingToolkit: check_variables, check_parameters, _iszero, _merge, chec get_unit, check_equations import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show -import MacroTools, Graphs +import MacroTools +import Graphs, Graphs.DiGraph, Graphs.SimpleGraph, Graphs.vertices, Graphs.edges, + Graphs.add_vertices!, Graphs.nv, Graphs.ne import DataStructures: OrderedDict, OrderedSet import Parameters: @with_kw_noshow @@ -68,6 +71,12 @@ export @reaction_network, @add_reactions, @reaction, @species include("registered_functions.jl") export mm, mmr, hill, hillr, hillar +# spatial reaction networks +include("lattice_reaction_system_diffusion.jl") +export DiffusionReaction, diffusion_reactions, isdiffusionparameter +export LatticeReactionSystem +export compartment_parameters, diffusion_parameters, diffusion_species + # functions to query network properties include("networkapi.jl") export species, nonspecies, reactionparams, reactions, speciesmap, paramsmap diff --git a/src/compound.jl b/src/compound.jl index 9f2ac5bca2..01c9081232 100644 --- a/src/compound.jl +++ b/src/compound.jl @@ -18,12 +18,12 @@ macro compound(species_expr, arr_expr...) # Construct the expressions that define the species species_expr = Expr(:macrocall, Symbol("@species"), LineNumberNode(0), - Expr(:call, species_name, species_arg)) + Expr(:call, species_name, species_arg)) # Construct the expression to set the iscompound metadata setmetadata_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), - Catalyst.CompoundSpecies, - true)) + Catalyst.CompoundSpecies, + true)) # Ensure the expressions are evaluated in the correct scope by escaping them escaped_species_expr = esc(species_expr) @@ -49,22 +49,22 @@ macro compound(species_expr, arr_expr...) # Construct the expression to set the components metadata setcomponents_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), - Catalyst.CompoundComponents, - $species_expr)) + Catalyst.CompoundComponents, + $species_expr)) # Ensure the expression is evaluated in the correct scope by escaping it escaped_setcomponents_expr = esc(setcomponents_expr) # Construct the expression to set the coefficients metadata setcoefficients_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), - Catalyst.CompoundCoefficients, - $coeffs_expr)) + Catalyst.CompoundCoefficients, + $coeffs_expr)) escaped_setcoefficients_expr = esc(setcoefficients_expr) # Return a block that contains the escaped expressions return Expr(:block, escaped_species_expr, escaped_setmetadata_expr, - escaped_setcomponents_expr, escaped_setcoefficients_expr) + escaped_setcomponents_expr, escaped_setcoefficients_expr) end # Check if a species is a compound diff --git a/src/lattice_reaction_system_diffusion.jl b/src/lattice_reaction_system_diffusion.jl new file mode 100644 index 0000000000..014e871cd2 --- /dev/null +++ b/src/lattice_reaction_system_diffusion.jl @@ -0,0 +1,266 @@ +### Diffusion Reaction Structure. ### + +# Implements the diffusionparameter metadata field. +struct DiffusionParameter end +Symbolics.option_to_metadata_type(::Val{:diffusionparameter}) = DiffusionParameter + +isdiffusionparameter(x::Num, args...) = isdiffusionparameter(Symbolics.unwrap(x), args...) +function isdiffusionparameter(x, default = false) + p = Symbolics.getparent(x, nothing) + p === nothing || (x = p) + Symbolics.getmetadata(x, DiffusionParameter, default) +end + +# Abstract spatial reaction structures. +abstract type AbstractSpatialReaction end + +# A diffusion reaction. These are simple to hanlde, and should cover most types of spatial reactions. +# Currently only permit constant rates. +struct DiffusionReaction <: AbstractSpatialReaction + """The rate function (excluding mass action terms). Currently only constants supported""" + rate::Num + """The species that is subject to difusion.""" + species::Num + """A symbol representation of the species that is subject to difusion.""" + species_sym::Symbol # Required for identification in certain cases. + + # Creates a diffusion reaction. + function DiffusionReaction(rate::Num, species::Num) + new(rate, species, ModelingToolkit.getname(species)) + end + function DiffusionReaction(rate::Number, species::Num) + new(Num(rate), species, ModelingToolkit.getname(species)) + end + function DiffusionReaction(rate::Symbol, species::Num) + new(Symbolics.variable(rate), species, ModelingToolkit.getname(species)) + end + function DiffusionReaction(rate::Num, species::Symbol) + new(rate, Symbolics.variable(species), species) + end + function DiffusionReaction(rate::Number, species::Symbol) + new(Num(rate), Symbolics.variable(species), species) + end + function DiffusionReaction(rate::Symbol, species::Symbol) + new(Symbolics.variable(rate), Symbolics.variable(species), species) + end +end +# Creates a vector of DiffusionReactions. +function diffusion_reactions(diffusion_reactions) + [DiffusionReaction(dr[1], dr[2]) for dr in diffusion_reactions] +end +# Gets the parameters in a diffusion reaction. +ModelingToolkit.parameters(dr::DiffusionReaction) = Symbolics.get_variables(dr.rate) + +### Lattice Reaction Network Structure ### +# Desribes a spatial reaction network over a graph. +struct LatticeReactionSystem # <: MT.AbstractTimeDependentSystem # Adding this part messes up show, disabling me from creating LRSs + """The reaction system within each comaprtment.""" + rs::ReactionSystem + """The spatial reactions defined between individual nodes.""" + spatial_reactions::Vector{<:AbstractSpatialReaction} + """The graph on which the lattice is defined.""" + lattice::DiGraph + + # Derrived values. + """The number of compartments.""" + nC::Int64 + """The number of edges.""" + nE::Int64 + """The number of species.""" + nS::Int64 + """Whenever the initial input was a di graph.""" + init_digraph::Bool + + function LatticeReactionSystem(rs::ReactionSystem, + spatial_reactions::Vector{<:AbstractSpatialReaction}, + lattice::DiGraph; init_digraph = true) + return new(rs, spatial_reactions, lattice, nv(lattice), ne(lattice), + length(species(rs)), init_digraph) + end + function LatticeReactionSystem(rs::ReactionSystem, + spatial_reactions::Vector{<:AbstractSpatialReaction}, + lattice::SimpleGraph) + return LatticeReactionSystem(rs, spatial_reactions, graph_to_digraph(lattice); + init_digraph = false) + end + function LatticeReactionSystem(rs::ReactionSystem, + spatial_reaction::AbstractSpatialReaction, + lattice::Graphs.AbstractGraph) + return LatticeReactionSystem(rs, [spatial_reaction], lattice) + end +end +# Covnerts a graph to a digraph (in a way where we know where the new edges are in teh edge vector). +function graph_to_digraph(g1) + g2 = Graphs.SimpleDiGraphFromIterator(reshape(permutedims(hcat(collect(edges(g1)), + reverse.(edges(g1)))), :, + 1)[:]) + add_vertices!(g2, nv(g1) - nv(g2)) + return g2 +end +# Gets the species of a lattice reaction system. +species(lrs::LatticeReactionSystem) = species(lrs.rs) +function diffusion_species(lrs::LatticeReactionSystem) + filter(s -> ModelingToolkit.getname(s) in getfield.(lrs.spatial_reactions, :species_sym), + species(lrs.rs)) +end + +# Gets the parameters in a lattice reaction system. +function ModelingToolkit.parameters(lrs::LatticeReactionSystem) + unique(vcat(parameters(lrs.rs), + Symbolics.get_variables.(getfield.(lrs.spatial_reactions, :rate))...)) +end +function compartment_parameters(lrs::LatticeReactionSystem) + filter(p -> !is_spatial_param(p, lrs), parameters(lrs)) +end +function diffusion_parameters(lrs::LatticeReactionSystem) + filter(p -> is_spatial_param(p, lrs), parameters(lrs)) +end + +# Checks whenever a parameter is a spatial parameter or not. +function is_spatial_param(p, lrs) + hasmetadata(p, DiffusionParameter) && getmetadata(p, DiffusionParameter) && + (return true) # Wanted to just depend on metadata, but seems like we cannot implement that trivially. + return (any(isequal.(p, parameters(lrs.rs))) ? false : true) +end + +### Processes Input u0 & p ### + +# From u0 input, extracts their values and store them in the internal format. +function lattice_process_u0(u0_in, u0_symbols, nC) + u0 = lattice_process_input(u0_in, u0_symbols, nC) + check_vector_lengths(u0, nC) + expand_component_values(u0, nC) +end + +# From p input, splits it into diffusion parameters and compartment parameters, and store these in the desired internal format. +function lattice_process_p(p_in, p_comp_symbols, p_diff_symbols, lrs::LatticeReactionSystem) + pC_in, pD_in = split_parameters(p_in, p_comp_symbols, p_diff_symbols) + pC = lattice_process_input(pC_in, p_comp_symbols, lrs.nC) + pD = lattice_process_input(pD_in, p_diff_symbols, lrs.nE) + lrs.init_digraph || foreach(idx -> duplicate_diff_params!(pD, idx, lrs), 1:length(pD)) + check_vector_lengths(pC, lrs.nC) + check_vector_lengths(pD, lrs.nE) + return pC, pD +end + +# Splits parameters into those for the compartments and those for the connections. +split_parameters(ps::Tuple{<:Any, <:Any}, args...) = ps +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 +function split_parameters(ps::Vector{<:Pair}, p_comp_symbols::Vector, + p_diff_symbols::Vector) + pC_in = [p for p in ps if Symbol(p[1]) in p_comp_symbols] + pD_in = [p for p in ps if Symbol(p[1]) in p_diff_symbols] + (sum(length.([pC_in, pD_in])) != length(ps)) && + error("These input parameters are not recongised: $(setdiff(first.(ps), vcat(first.([pC_in, pE_in]))))") + return pC_in, pD_in +end + +# If the input is given in a map form, teh vector needs sorting and the first value removed. +function lattice_process_input(input::Vector{<:Pair}, symbols::Vector{Symbol}, args...) + (length(setdiff(Symbol.(first.(input)), symbols)) != 0) && + error("Some input symbols are not recognised: $(setdiff(Symbol.(first.(input)), symbols)).") + sorted_input = sort(input; + by = p -> findfirst(ModelingToolkit.getname(p[1]) .== symbols)) + return lattice_process_input(last.(sorted_input), symbols, args...) +end +# Processes the input and gvies it in a form where it is a vector of vectors (some of which may have a single value). +function lattice_process_input(input::Matrix{<:Number}, args...) + lattice_process_input([vec(input[i, :]) for i in 1:size(input, 1)], args...) +end +function lattice_process_input(input::Array{<:Number, 3}, args...) + error("3 dimensional array parameter inpur currently not supported.") +end +function lattice_process_input(input::Vector{<:Any}, args...) + isempty(input) && return Vector{Vector{Float64}}() + lattice_process_input([(val isa Vector{<:Number}) ? val : [val] for val in input], + args...) +end +lattice_process_input(input::Vector{<:Vector}, symbols::Vector{Symbol}, n::Int64) = input +function check_vector_lengths(input::Vector{<:Vector}, n) + isempty(setdiff(unique(length.(input)), [1, n])) || + error("Some inputs where given values of inappropriate length.") +end + +# For diffusion parameters, if the graph was given as an undirected graph of length n, and the paraemter have n values, exapnd so that the same value are given for both values on the edge. +function duplicate_diff_params!(pD::Vector{Vector{Float64}}, idx::Int64, + lrs::LatticeReactionSystem) + (2length(pD[idx]) == lrs.nE) && (pD[idx] = [p_val for p_val in pD[idx] for _ in 1:2]) +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(pC, pD, lrs) + merge(vals_to_dict(compartment_parameters(lrs), pC), + vals_to_dict(diffusion_parameters(lrs), pD)) +end + +# Computes the diffusion rates and stores them in a format (Dictionary of species index to rates across all edges). +function compute_all_diffusion_rates(pC::Vector{Vector{Float64}}, + pD::Vector{Vector{Float64}}, + lrs::LatticeReactionSystem) + param_value_dict = param_dict(pC, pD, lrs) + return [s => Symbolics.value.(compute_diffusion_rates(get_diffusion_rate_law(s, lrs), + param_value_dict, lrs.nE)) + for s in diffusion_species(lrs)] +end +function get_diffusion_rate_law(s::Symbolics.BasicSymbolic, lrs::LatticeReactionSystem) + rates = filter(sr -> isequal(ModelingToolkit.getname(s), sr.species_sym), + lrs.spatial_reactions) + (length(rates) > 1) && error("Species $s have more than one diffusion reaction.") # We could allows several and simply sum them though, easy change. + return rates[1].rate +end +function compute_diffusion_rates(rate_law::Num, + param_value_dict::Dict{Any, Vector{Float64}}, nE::Int64) + relevant_parameters = Symbolics.get_variables(rate_law) + if all(length(param_value_dict[P]) == 1 for P in relevant_parameters) + return [ + substitute(rate_law, + Dict(p => param_value_dict[p][1] for p in relevant_parameters)), + ] + end + return [substitute(rate_law, + Dict(p => get_component_value(param_value_dict[p], idxE) + for p in relevant_parameters)) for idxE in 1:nE] +end + +### Accessing State & Parameter Array Values ### + +# Gets the index in the u array of species s in compartment comp (when their are nS species). +get_index(comp::Int64, s::Int64, nS::Int64) = (comp - 1) * nS + s +# Gets the indexes in the u array of all species in comaprtment comp (when their are nS species). +get_indexes(comp::Int64, nS::Int64) = ((comp - 1) * nS + 1):(comp * nS) + +# We have many vectors of length 1 or n, for which we want to get value idx (or the one value, if length is 1), this function gets that. +function get_component_value(values::Vector{<:Vector}, component_idx::Int64, + location_idx::Int64) + get_component_value(values[component_idx], location_idx) +end +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 +function get_component_value(values::Vector{<:Number}, location_idx::Int64) + get_component_value(values, location_idx, length(values) == 1) +end +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. +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 pC vector at a given comaprtment. +function view_pC_vector(pC, comp, pC_location_types, pC_idxs) + mapview(p_idx -> pC_location_types[p_idx] ? pC[p_idx][1] : pC[p_idx][comp], pC_idxs) +end +# Expands a u0/p information stored in Vector{Vector{}} for to Matrix form (currently 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/runtests.jl b/test/runtests.jl index 2bbfa84495..f0f500fa26 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,7 +36,11 @@ using SafeTestsets @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_system_base.jl") end + #@time @safetestset "Lattice Systems Simulations" begin include("spatial_reaction_systems/lattice_reaction_systems.jl") end ### Tests network visualization. ### @time @safetestset "Latexify" begin include("visualization/latexify.jl") end diff --git a/test/spatial_reaction_systems/lattice_reaction_system_base.jl b/test/spatial_reaction_systems/lattice_reaction_system_base.jl new file mode 100644 index 0000000000..998552756b --- /dev/null +++ b/test/spatial_reaction_systems/lattice_reaction_system_base.jl @@ -0,0 +1,47 @@ +# Fetch packages. +using Catalyst, JumpProcesses, OrdinaryDiffEq +using Random, Statistics, SparseArrays, Test +using Graphs + +### Creates a Model ### +SIR_system = @reaction_network begin + α, S + I --> 2I + β, I --> R +end + +SIR_srs_1 = [DiffusionReaction(:dS, :S)] +SIR_srs_2 = diffusion_reactions([(:dS, :S), (:dI, :I), (:dR, :R)]) + +grid = Graphs.grid([2, 2]) + +lrs1 = LatticeReactionSystem(SIR_system, SIR_srs_1, grid) +lrs2 = LatticeReactionSystem(SIR_system, SIR_srs_2, grid) + +@test length(compartment_parameters(lrs1)) == length(compartment_parameters(lrs2)) == 2 +@test length(diffusion_parameters(lrs1)) == length(diffusion_species(lrs1)) == 1 +@test length(diffusion_parameters(lrs2)) == length(diffusion_species(lrs2)) == 3 + +### Check Parameters and u0 Handling ### +l = Graphs.nv(grid) +SIR_pC = [:α => 0.1 / 1000, :β => 0.01] +SIR_pD = [:dS => 0.01] +SIR_u0 = [:S => 999.0, :I => 1.0, :R => 0.0] +pC_syms = [:α, :β] +pD_syms = [:dS] +u0_syms = [:S, :I, :R] + +# Split parameters. +@test Catalyst.split_parameters((SIR_pC, SIR_pD)) == (SIR_pC, SIR_pD) +@test Catalyst.split_parameters([SIR_pC; SIR_pD], [:α, :β], [:dS]) == (SIR_pC, SIR_pD) + +# Convert input to vector of vector (of values). +@test Catalyst.lattice_process_input([:S => 1, :I => 2, :R => 3], u0_syms, l) == [[1], [2], [3]] +@test Catalyst.lattice_process_input([:S => 1, :I => [2, 2, 2, 2], :R => 3], u0_syms, l) == [[1], [2, 2, 2, 2], [3]] +@test Catalyst.lattice_process_input([1 1 1 1; 2 2 2 2; 3 3 3 3], u0_syms, l) == [[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3]] +@test Catalyst.lattice_process_input([1, 2, 3], u0_syms, l) == [[1], [2], [3]] +@test Catalyst.lattice_process_input([1, [2, 2, 2, 2], 3], u0_syms, l) == [[1], [2, 2, 2, 2], [3]] +@test Catalyst.lattice_process_input([1, [1,2,3,4],1], u0_syms, l) + +# Fetching teh diffusion rates. +@test isequal(Catalyst.compute_all_diffusion_rates([[1.0], [2.0, 2.0, 2.0, 2.0], [3.0]], [[0.1]], lrs1), Pair.(diffusion_species(lrs1),[[0.1]])) +@test isequal(Catalyst.compute_all_diffusion_rates([[1.0], [2.0, 2.0, 2.0, 2.0], [3.0]], [[0.1], [0.2], [0.3]], lrs2), Pair.(diffusion_species(lrs2),[[0.1], [0.2], [0.3]])) 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