Merge branch 'master' into fix_doc_project_versions

docs/pages.jl

@@ -18,7 +18,7 @@ pages = Any[
         #"model_creation/network_analysis.md",
         "model_creation/chemistry_related_functionality.md",
         "Model creation examples" => Any[
-            "model_creation/examples/basic_CRN_library.md",
+            #"model_creation/examples/basic_CRN_library.md",
             "model_creation/examples/programmatic_generative_linear_pathway.md",
             #"model_creation/examples/hodgkin_huxley_equation.md",
             #"model_creation/examples/smoluchowski_coagulation_equation.md"
@@ -28,7 +28,7 @@ pages = Any[
         "model_simulation/simulation_introduction.md",
         "model_simulation/simulation_plotting.md",
         "model_simulation/simulation_structure_interfacing.md",
-        "model_simulation/ensemble_simulations.md",
+        #"model_simulation/ensemble_simulations.md",
         # Stochastic simulation statistical analysis.
         "model_simulation/ode_simulation_performance.md",
         # SDE Performance considerations/advice.
@@ -37,7 +37,7 @@ pages = Any[
     ],
     "Steady state analysis" => Any[
         "steady_state_functionality/homotopy_continuation.md",
-        "steady_state_functionality/nonlinear_solve.md",
+        #"steady_state_functionality/nonlinear_solve.md",
         "steady_state_functionality/steady_state_stability_computation.md",        
         "steady_state_functionality/bifurcation_diagrams.md",
         "steady_state_functionality/dynamical_systems.md"

src/Catalyst.jl

@@ -6,9 +6,9 @@ module Catalyst
 using DocStringExtensions
 using SparseArrays, DiffEqBase, Reexport, Setfield
 using LaTeXStrings, Latexify, Requires
-using JumpProcesses: JumpProcesses,
-                     JumpProblem, MassActionJump, ConstantRateJump,
-                     VariableRateJump
+using JumpProcesses: JumpProcesses, JumpProblem, 
+                     MassActionJump, ConstantRateJump, VariableRateJump,
+                     SpatialMassActionJump
 
 # ModelingToolkit imports and convenience functions we use
 using ModelingToolkit
@@ -165,20 +165,21 @@ export make_si_ode
 
 ### Spatial Reaction Networks ###
 
-# spatial reactions
+# Spatial reactions.
 include("spatial_reaction_systems/spatial_reactions.jl")
 export TransportReaction, TransportReactions, @transport_reaction
 export isedgeparameter
 
-# lattice reaction systems
+# Lattice reaction systems
 include("spatial_reaction_systems/lattice_reaction_systems.jl")
 export LatticeReactionSystem
 export spatial_species, vertex_parameters, edge_parameters
 
-# variosu utility functions
+# Various utility functions
 include("spatial_reaction_systems/utility.jl")
 
-# spatial lattice ode systems.
+# Specific spatial problem types.
 include("spatial_reaction_systems/spatial_ODE_systems.jl")
+include("spatial_reaction_systems/lattice_jump_systems.jl")
 
 end # module

src/spatial_reaction_systems/lattice_jump_systems.jl

@@ -0,0 +1,136 @@
+### JumpProblem ###
+
+# Builds a spatial DiscreteProblem from a Lattice Reaction System.
+function DiffEqBase.DiscreteProblem(lrs::LatticeReactionSystem, u0_in, tspan, p_in = DiffEqBase.NullParameters(), args...; kwargs...)
+    if !is_transport_system(lrs) 
+        error("Currently lattice Jump simulations only supported when all spatial reactions are transport reactions.")
+    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)   
+
+    # Returns a DiscreteProblem.
+    # Previously, a Tuple was used for (vert_ps, edge_ps), but this was converted to a Vector internally.
+    return DiscreteProblem(u0, tspan, [vert_ps, edge_ps], args...; kwargs...)
+end
+
+# Builds a spatial JumpProblem from a DiscreteProblem containg a Lattice Reaction System.
+function JumpProcesses.JumpProblem(lrs::LatticeReactionSystem, dprob, aggregator, args...; name = nameof(lrs.rs), 
+                                   combinatoric_ratelaws = get_combinatoric_ratelaws(lrs.rs), kwargs...)
+    # Error checks.
+    if !isnothing(dprob.f.sys)
+        error("Unexpected `DiscreteProblem` passed into `JumpProblem`. Was a `LatticeReactionSystem` used as input to the initial `DiscreteProblem`?")
+    end
+
+    # Computes hopping constants and mass action jumps (requires some internal juggling).
+    # Currently, JumpProcesses requires uniform vertex parameters (hence `p=first.(dprob.p[1])`).
+    # Currently, the resulting JumpProblem does not depend on parameters (no way to incorporate these).
+    # Hence the parameters of this one does nto actually matter. If at some point JumpProcess can 
+    # handle parameters this can be updated and improved.
+    # The non-spatial DiscreteProblem have a u0 matrix with entries for all combinations of species and vertexes.
+    hopping_constants = make_hopping_constants(dprob, lrs)
+    sma_jumps = make_spatial_majumps(dprob, lrs)
+    non_spat_dprob = DiscreteProblem(reshape(dprob.u0, lrs.num_species, lrs.num_verts), dprob.tspan, first.(dprob.p[1]))
+
+    return JumpProblem(non_spat_dprob, aggregator, sma_jumps; 
+                       hopping_constants, spatial_system = lrs.lattice, name, kwargs...)
+end
+
+# Creates the hopping constants from a discrete problem and a lattice reaction system.
+function make_hopping_constants(dprob::DiscreteProblem, lrs::LatticeReactionSystem)
+    # Creates the all_diff_rates vector, containing for each species, its transport rate across all edges.
+    # If transport rate is uniform for one species, the vector have a single element, else one for each edge.
+    spatial_rates_dict = Dict(compute_all_transport_rates(dprob.p[1], dprob.p[2], lrs))
+    all_diff_rates = [haskey(spatial_rates_dict, s) ? spatial_rates_dict[s] : [0.0] for s in species(lrs)]
+
+    # Creates the hopping constant Matrix. It contains one element for each combination of species and vertex.
+    # Each element is a Vector, containing the outgoing hopping rates for that species, from that vertex, on that edge.
+    hopping_constants = [Vector{Float64}(undef, length(lrs.lattice.fadjlist[j])) 
+                         for i in 1:(lrs.num_species), j in 1:(lrs.num_verts)]
+
+    # For each edge, finds each position in `hopping_constants`.
+    for (e_idx, e) in enumerate(edges(lrs.lattice))
+        dst_idx = findfirst(isequal(e.dst), lrs.lattice.fadjlist[e.src])      
+        # For each species, sets that hopping rate.  
+        for s_idx in 1:(lrs.num_species)
+            hopping_constants[s_idx, e.src][dst_idx] = get_component_value(all_diff_rates[s_idx], e_idx)
+        end
+    end
+
+    return hopping_constants
+end
+
+# Creates a SpatialMassActionJump struct from a (spatial) DiscreteProblem and a LatticeReactionSystem.
+# Could implementation a version which, if all reaction's rates are uniform, returns a MassActionJump.
+# Not sure if there is any form of performance improvement from that though. Possibly is not the case.
+function make_spatial_majumps(dprob, lrs::LatticeReactionSystem)
+    # Creates a vector, storing which reactions have spatial components.
+    is_spatials = [Catalyst.has_spatial_vertex_component(rx.rate, lrs; vert_ps = dprob.p[1]) for rx in reactions(lrs.rs)]
+
+    # Creates templates for the rates (uniform and spatial) and the stoichiometries.
+    # We cannot fetch reactant_stoich and net_stoich from a (non-spatial) MassActionJump.
+    # The reason is that we need to re-order the reactions so that uniform appears first, and spatial next.
+    u_rates = Vector{Float64}(undef, length(reactions(lrs.rs)) - count(is_spatials))
+    s_rates = Matrix{Float64}(undef, count(is_spatials), lrs.num_verts)
+    reactant_stoich = Vector{Vector{Pair{Int64, Int64}}}(undef, length(reactions(lrs.rs)))
+    net_stoich = Vector{Vector{Pair{Int64, Int64}}}(undef, length(reactions(lrs.rs)))
+
+    # Loops through reactions with non-spatial rates, computes their rates and stoichiometries.
+    cur_rx = 1;
+    for (is_spat, rx) in zip(is_spatials, reactions(lrs.rs))
+        is_spat && continue
+        u_rates[cur_rx] = compute_vertex_value(rx.rate, lrs; vert_ps = dprob.p[1])[1]
+        substoich_map = Pair.(rx.substrates, rx.substoich)
+        reactant_stoich[cur_rx] = int_map(substoich_map, lrs.rs)
+        net_stoich[cur_rx] = int_map(rx.netstoich, lrs.rs)
+        cur_rx += 1
+    end
+    # Loops through reactions with spatial rates, computes their rates and stoichiometries.
+    for (is_spat, rx) in zip(is_spatials, reactions(lrs.rs))    
+        is_spat || continue
+        s_rates[cur_rx-length(u_rates),:] = compute_vertex_value(rx.rate, lrs; vert_ps = dprob.p[1])
+        substoich_map = Pair.(rx.substrates, rx.substoich)
+        reactant_stoich[cur_rx] = int_map(substoich_map, lrs.rs)
+        net_stoich[cur_rx] = int_map(rx.netstoich, lrs.rs)
+        cur_rx += 1
+    end
+    # SpatialMassActionJump expects empty rate containers to be nothing.
+    isempty(u_rates) && (u_rates = nothing)
+    (count(is_spatials)==0) && (s_rates = nothing)
+
+    return SpatialMassActionJump(u_rates, s_rates, reactant_stoich, net_stoich)
+end
+
+### Extra ###
+
+# Temporary. Awaiting implementation in SII, or proper implementation withinCatalyst (with more general functionality).
+function int_map(map_in, sys) where {T,S}
+    return [ModelingToolkit.variable_index(sys, pair[1]) => pair[2] for pair in map_in]
+end
+
+# Currently unused. If we want to create certain types of MassActionJumps (instead of SpatialMassActionJumps) we can take this one back.
+# Creates the (non-spatial) mass action jumps from a (non-spatial) DiscreteProblem (and its Reaction System of origin).
+# function make_majumps(non_spat_dprob, rs::ReactionSystem)
+#     # Computes various required inputs for assembling the mass action jumps.
+#     js = convert(JumpSystem, rs)
+#     statetoid = Dict(ModelingToolkit.value(state) => i for (i, state) in enumerate(states(rs)))
+#     eqs = equations(js)
+#     invttype = non_spat_dprob.tspan[1] === nothing ? Float64 : typeof(1 / non_spat_dprob.tspan[2])
+# 
+#     # Assembles the non-spatial mass action jumps.
+#     p = (non_spat_dprob.p isa DiffEqBase.NullParameters || non_spat_dprob.p === nothing) ? Num[] : non_spat_dprob.p    
+#     majpmapper = ModelingToolkit.JumpSysMajParamMapper(js, p; jseqs = eqs, rateconsttype = invttype)
+#     return ModelingToolkit.assemble_maj(eqs.x[1], statetoid, majpmapper)
+# end

src/spatial_reaction_systems/utility.jl

@@ -295,3 +295,99 @@ end
 function matrix_expand_component_values(values::Vector{<:Vector}, n)
     reshape(expand_component_values(values, n), length(values), n)
 end
+
+# For an expression, computes its values using the provided state and parameter vectors.
+# The expression is assumed to be valid in edges (and can have edges parameter components).
+# If some component is non-uniform, output is a vector of length equal to the number of vertexes.
+# If all components are uniform, the output is a length one vector.
+function compute_edge_value(exp, lrs::LatticeReactionSystem, edge_ps)
+    # Finds the symbols in the expression. Checks that all correspond to edge parameters.
+    relevant_syms = Symbolics.get_variables(exp)
+    if !all(any(isequal(sym, p) for p in edge_parameters(lrs)) for sym in relevant_syms) 
+        error("An non-edge parameter was encountered in expressions: $exp. Here, only edge parameters are expected.")
+    end
+
+    # Creates a Function tha computes the expressions value for a parameter set.
+    exp_func = drop_expr(@RuntimeGeneratedFunction(build_function(exp, relevant_syms...)))
+    # Creates a dictionary with the value(s) for all edge parameters.
+    sym_val_dict = vals_to_dict(edge_parameters(lrs), edge_ps)
+
+    # If all values are uniform, compute value once. Else, do it at all edges.
+    if !has_spatial_edge_component(exp, lrs, edge_ps)
+        return [exp_func([sym_val_dict[sym][1] for sym in relevant_syms]...)]
+    end
+    return [exp_func([get_component_value(sym_val_dict[sym], idxE) for sym in relevant_syms]...) 
+                                                                            for idxE in 1:lrs.num_edges]
+end
+
+# For an expression, computes its values using the provided state and parameter vectors.
+# The expression is assumed to be valid in vertexes (and can have vertex parameter and state components).
+# If at least one component is non-uniform, output is a vector of length equal to the number of vertexes.
+# If all components are uniform, the output is a length one vector.
+function compute_vertex_value(exp, lrs::LatticeReactionSystem; u=nothing, vert_ps=nothing)
+    # Finds the symbols in the expression. Checks that all correspond to states or vertex parameters.
+    relevant_syms = Symbolics.get_variables(exp)
+    if any(any(isequal(sym) in edge_parameters(lrs)) for sym in relevant_syms) 
+        error("An edge parameter was encountered in expressions: $exp. Here, on vertex-based components are expected.")
+    end
+    # Creates a Function tha computes the expressions value for a parameter set.
+    exp_func = drop_expr(@RuntimeGeneratedFunction(build_function(exp, relevant_syms...)))
+    # Creates a dictionary with the value(s) for all edge parameters.
+    if !isnothing(u) && !isnothing(vert_ps)
+        all_syms = [species(lrs); vertex_parameters(lrs)]
+        all_vals = [u; vert_ps]
+    elseif !isnothing(u) && isnothing(vert_ps)
+        all_syms = species(lrs)
+        all_vals = u
+
+    elseif isnothing(u) && !isnothing(vert_ps)
+        all_syms = vertex_parameters(lrs)
+        all_vals = vert_ps
+    else
+        error("Either u or vertex_ps have to be provided to has_spatial_vertex_component.")
+    end
+    sym_val_dict = vals_to_dict(all_syms, all_vals)
+
+    # If all values are uniform, compute value once. Else, do it at all edges.
+    if !has_spatial_vertex_component(exp, lrs; u, vert_ps)
+        return [exp_func([sym_val_dict[sym][1] for sym in relevant_syms]...)]
+    end
+    return [exp_func([get_component_value(sym_val_dict[sym], idxV) for sym in relevant_syms]...) 
+                                                                            for idxV in 1:lrs.num_verts]
+end
+
+### System Property Checks ###
+
+# For a Symbolic expression, a LatticeReactionSystem, and a parameter list of the internal format:
+# Checks if any edge parameter in the expression have a spatial component (that is, is not uniform).
+function has_spatial_edge_component(exp, lrs::LatticeReactionSystem, edge_ps)
+    # Finds the edge parameters in the expression. Computes their indexes.
+    exp_syms = Symbolics.get_variables(exp)
+    exp_edge_ps = filter(sym -> any(isequal(sym), edge_parameters(lrs)), exp_syms)
+    p_idxs = [findfirst(isequal(sym, edge_p) for edge_p in edge_parameters(lrs)) for sym in exp_syms]
+    # Checks if any of the corresponding value vectors have length != 1 (that is, is not uniform).
+    return any(length(edge_ps[p_idx]) != 1 for p_idx in p_idxs)
+end
+
+# For a Symbolic expression, a LatticeReactionSystem, and a parameter list of the internal format (vector of vectors):
+# Checks if any vertex parameter in the expression have a spatial component (that is, is not uniform).
+function has_spatial_vertex_component(exp, lrs::LatticeReactionSystem; u=nothing, vert_ps=nothing)
+    # Finds all the symbols in the expression.
+    exp_syms = Symbolics.get_variables(exp)
+
+    # If vertex parameter values where given, checks if any of these have non-uniform values.
+    if !isnothing(vert_ps)
+        exp_vert_ps = filter(sym -> any(isequal(sym), vertex_parameters(lrs)), exp_syms)
+        p_idxs = [ModelingToolkit.parameter_index(lrs.rs, sym) for sym in exp_vert_ps]
+        any(length(vert_ps[p_idx]) != 1 for p_idx in p_idxs) && return true
+    end
+
+    # If states values where given, checks if any of these have non-uniform values.
+    if !isnothing(u)
+        exp_u = filter(sym -> any(isequal(sym), species(lrs)), exp_syms)
+        u_idxs = [ModelingToolkit.variable_index(lrs.rs, sym) for sym in exp_u]
+        any(length(u[u_idx]) != 1 for u_idx in u_idxs) && return true
+    end
+
+    return false
+end

test/runtests.jl

@@ -61,6 +61,7 @@ using SafeTestsets, Test
         @time @safetestset "PDE Systems Simulations" begin include("spatial_modelling/simulate_PDEs.jl") end
         @time @safetestset "Lattice Reaction Systems" begin include("spatial_modelling/lattice_reaction_systems.jl") end
         @time @safetestset "ODE Lattice Systems Simulations" begin include("spatial_modelling/lattice_reaction_systems_ODEs.jl") end
+        @time @safetestset "Jump Lattice Systems Simulations" begin include("spatial_reaction_systems/lattice_reaction_systems_jumps.jl") end
     #end
 
     #if GROUP == "All" || GROUP == "Visualisation-Extensions"