@@ -9,25 +9,27 @@ function DiffEqBase.DiscreteProblem(lrs::LatticeReactionSystem, u0_in, tspan,
9
9
10
10
# Converts potential symmaps to varmaps.
11
11
u0_in = symmap_to_varmap (lrs, u0_in)
12
- p_in = symmap_to_varmap (lrs, p_in)
12
+ p_in = symmap_to_varmap (lrs, p_in)
13
13
14
14
# Converts u0 and p to their internal forms.
15
15
# u0 is simply a vector with all the species' initial condition values across all vertices.
16
16
# u0 is [spec 1 at vert 1, spec 2 at vert 1, ..., spec 1 at vert 2, ...].
17
- u0 = lattice_process_u0 (u0_in, species (lrs), lrs)
17
+ u0 = lattice_process_u0 (u0_in, species (lrs), lrs)
18
18
# vert_ps and `edge_ps` are vector maps, taking each parameter's Symbolics representation to its value(s).
19
19
# vert_ps values are vectors. Here, index (i) is a parameter's value in vertex i.
20
20
# edge_ps values are sparse matrices. Here, index (i,j) is a parameter's value in the edge from vertex i to vertex j.
21
21
# Uniform vertex/edge parameters store only a single value (a length 1 vector, or size 1x1 sparse matrix).
22
- vert_ps, edge_ps = lattice_process_p (p_in, vertex_parameters (lrs), edge_parameters (lrs), lrs)
22
+ vert_ps, edge_ps = lattice_process_p (p_in, vertex_parameters (lrs), edge_parameters (lrs),
23
+ lrs)
23
24
24
25
# Returns a DiscreteProblem (which basically just stores the processed input).
25
26
return DiscreteProblem (u0, tspan, [vert_ps; edge_ps], args... ; kwargs... )
26
27
end
27
28
28
29
# Builds a spatial JumpProblem from a DiscreteProblem containing a `LatticeReactionSystem`.
29
- function JumpProcesses. JumpProblem (lrs:: LatticeReactionSystem , dprob, aggregator, args... ; name = nameof (reactionsystem (lrs)),
30
- combinatoric_ratelaws = get_combinatoric_ratelaws (reactionsystem (lrs)), kwargs... )
30
+ function JumpProcesses. JumpProblem (lrs:: LatticeReactionSystem , dprob, aggregator,args... ;
31
+ combinatoric_ratelaws = get_combinatoric_ratelaws (reactionsystem (lrs)),
32
+ name = nameof (reactionsystem (lrs)), kwargs... )
31
33
# Error checks.
32
34
if ! isnothing (dprob. f. sys)
33
35
throw (ArgumentError (" Unexpected `DiscreteProblem` passed into `JumpProblem`. Was a `LatticeReactionSystem` used as input to the initial `DiscreteProblem`?" ))
@@ -40,12 +42,13 @@ function JumpProcesses.JumpProblem(lrs::LatticeReactionSystem, dprob, aggregator
40
42
# The non-spatial DiscreteProblem have a u0 matrix with entries for all combinations of species and vertexes.
41
43
hopping_constants = make_hopping_constants (dprob, lrs)
42
44
sma_jumps = make_spatial_majumps (dprob, lrs)
43
- non_spat_dprob = DiscreteProblem (reshape (dprob. u0, num_species (lrs), num_verts (lrs)), dprob. tspan, first .(dprob. p[1 ]))
45
+ non_spat_dprob = DiscreteProblem (reshape (dprob. u0, num_species (lrs), num_verts (lrs)),
46
+ dprob. tspan, first .(dprob. p[1 ]))
44
47
45
48
# Creates and returns a spatial JumpProblem (masked lattices are not supported by these).
46
49
spatial_system = has_masked_lattice (lrs) ? get_lattice_graph (lrs) : lattice (lrs)
47
- return JumpProblem (non_spat_dprob, aggregator, sma_jumps;
48
- hopping_constants, spatial_system , name, kwargs... )
50
+ return JumpProblem (non_spat_dprob, aggregator, sma_jumps;
51
+ hopping_constants, spatial_system, name, kwargs... )
49
52
end
50
53
51
54
# Creates the hopping constants from a discrete problem and a lattice reaction system.
@@ -59,18 +62,17 @@ function make_hopping_constants(dprob::DiscreteProblem, lrs::LatticeReactionSyst
59
62
# Creates an array (of the same size as the hopping constant array) containing all edges.
60
63
# First the array is a NxM matrix (number of species x number of vertices). Each element is a
61
64
# vector containing all edges leading out from that vertex (sorted by destination index).
62
- edge_array = [Pair{Int64,Int64}[] for _1 in 1 : num_species (lrs), _2 in 1 : num_verts (lrs)]
65
+ edge_array = [Pair{Int64, Int64}[] for _1 in 1 : num_species (lrs), _2 in 1 : num_verts (lrs)]
63
66
for e in edge_iterator (lrs), s_idx in 1 : num_species (lrs)
64
67
push! (edge_array[s_idx, e[1 ]], e)
65
68
end
66
69
foreach (e_vec -> sort! (e_vec; by = e -> e[2 ]), edge_array)
67
70
68
71
# Creates the hopping constants array. It has the same shape as the edge array, but each
69
72
# element is that species transportation rate along that edge
70
- hopping_constants = [
71
- [Catalyst. get_edge_value (all_diff_rates[s_idx], e) for e in edge_array[s_idx, src_idx]]
72
- for s_idx in 1 : num_species (lrs), src_idx in 1 : num_verts (lrs)
73
- ]
73
+ hopping_constants = [[Catalyst. get_edge_value (all_diff_rates[s_idx], e)
74
+ for e in edge_array[s_idx, src_idx]]
75
+ for s_idx in 1 : num_species (lrs), src_idx in 1 : num_verts (lrs)]
74
76
return hopping_constants
75
77
end
76
78
79
81
# Not sure if there is any form of performance improvement from that though. Likely not the case.
80
82
function make_spatial_majumps (dprob, lrs:: LatticeReactionSystem )
81
83
# Creates a vector, storing which reactions have spatial components.
82
- is_spatials = [has_spatial_vertex_component (rx. rate, dprob. p)
83
- for rx in reactions (reactionsystem (lrs))]
84
+ is_spatials = [has_spatial_vertex_component (rx. rate, dprob. p)
85
+ for rx in reactions (reactionsystem (lrs))]
84
86
85
87
# Creates templates for the rates (uniform and spatial) and the stoichiometries.
86
88
# We cannot fetch reactant_stoich and net_stoich from a (non-spatial) MassActionJump.
87
89
# The reason is that we need to re-order the reactions so that uniform appears first, and spatial next.
88
- u_rates = Vector {Float64} (undef, length (reactions (reactionsystem (lrs))) - count (is_spatials))
90
+ num_rxs = length (reactions (reactionsystem (lrs)))
91
+ u_rates = Vector {Float64} (undef, num_rxs - count (is_spatials))
89
92
s_rates = Matrix {Float64} (undef, count (is_spatials), num_verts (lrs))
90
- reactant_stoich = Vector {Vector{Pair{Int64, Int64}}} (undef, length ( reactions ( reactionsystem (lrs))) )
91
- net_stoich = Vector {Vector{Pair{Int64, Int64}}} (undef, length ( reactions ( reactionsystem (lrs))) )
93
+ reactant_stoich = Vector {Vector{Pair{Int64, Int64}}} (undef, num_rxs )
94
+ net_stoich = Vector {Vector{Pair{Int64, Int64}}} (undef, num_rxs )
92
95
93
96
# Loops through reactions with non-spatial rates, computes their rates and stoichiometries.
94
97
cur_rx = 1
@@ -101,9 +104,10 @@ function make_spatial_majumps(dprob, lrs::LatticeReactionSystem)
101
104
cur_rx += 1
102
105
end
103
106
# Loops through reactions with spatial rates, computes their rates and stoichiometries.
104
- for (is_spat, rx) in zip (is_spatials, reactions (reactionsystem (lrs)))
107
+ for (is_spat, rx) in zip (is_spatials, reactions (reactionsystem (lrs)))
105
108
is_spat || continue
106
- s_rates[cur_rx - length (u_rates), :] .= compute_vertex_value (rx. rate, lrs; ps = dprob. p)
109
+ s_rates[cur_rx - length (u_rates), :] .= compute_vertex_value (rx. rate, lrs;
110
+ ps = dprob. p)
107
111
substoich_map = Pair .(rx. substrates, rx. substoich)
108
112
reactant_stoich[cur_rx] = int_map (substoich_map, reactionsystem (lrs))
109
113
net_stoich[cur_rx] = int_map (rx. netstoich, reactionsystem (lrs))
0 commit comments