Catalyst.jl
Catalyst.jl copied to clipboard
[WIP] LatticeReactionSystem internal update
This is not complete yet, but I will initiate the PR so that I can summarise what I have lined up for the internal update of the LatticeReactionStructure
. The idea is to include all the needed stuff (except for using PDEBase to update the indexing of the solution).
This is an example of the new workflow:
using Catalyst, OrdinaryDiffEq
rn = @reaction_network begin
i, S + I --> 2I
r, I --> R
end
tr = @transport_reaction D S
lattice = CartesianGrid((20,20))
lrs = LatticeReactionSystem(rn, [tr], lattice)
I_values = zeros(20,20)
I_values[1,1] = 1
u0 = [:S => 1000.0, :I => I_values, :R => 0.0]
tspan = (0.0, 100.0)
ps = [:i => 0.01, :r => 0.001, :D => 0.05]
oprob = ODEProblem(lrs, u0, tspan, ps)
sol = solve(oprob, Tsit5())
This is what I will include/have included:
More or less Major
- Generalise how the lattice is stored. Now the lattice can be given as (1) A
CartesianGrid
, (2) an Array ofBool
, (3) An undirected graph, and (4) a directed graph (for the cartesian grid and array, I currently only allowed dimensions up to 3. - Change how edges (connections) are iterated through. Currently, we just use the list of edges from the graphs. Now I will instead generate a
edge_iterator
. This is an iterator over `Pair{Int64,Int64} representing the edges on the lattice. This is required to handle the 3 possible types of lattices now permitted. - Change how the edge parameter values are stored internally (in e.g.
ODEProblem
). Previously, for each parameter, we stored a vector with its values. Now we store a sparse matrix (where value i,j is its value from edge i to j). If the parameter is uniform, the sparse matrix is size 1x1, storing the parameters only value. - Input parameter and initial condition values are now only permitted in map (
[p => 1.0, d => 2.0]
) form (or corresponding dictionary form). - How non-uniform initial condition/parameter values are given has now been changed/improved:
- vertex parameter and initial condition values can now be a vector (of the same length as the number of vertexes). This is the same as before.
- For Cartesian/masked grids. They can also be given in array form (where the value in index i,j,k of the input array corresponds to the same value in the original lattice.
- Edge parameters should now be a sparse matrix with its value across all connections. Here, index (i,j) is the edge from vertex i to vertex j. Note that this means that flat/scalar indexes are used for Cartesian/masked grids. One could have allowed multi-dimensional sparse arrays (denoting the source and destination vertices' indexes on the grid), but this could result in 6d sparse arrays as input, which felt too messy.
- Added two functions to help the user generate edge parameter values (
make_edge_p_values
) andmake_directed_edge_values
.
Minor stuff
- The internal
view_vert_ps_vector!
function is now called directly onLatticeTransportODEf
/LatticeTransportODEjac
(instead of directly on their field). - Redone some of the workflow in utility.jl, i.e. where input initial conditions/parameters are processed to the correct format (to account for how this was changed).
- The
p
vector stored in theODEProblem
not contain both vertex and edge parameters (vertex parameters first, then edge parameters). - Added some new getter functions for fields in
LatticeReactionSystem
(getting information on the lattice).
Other notes
- I have currently not changed the order of storing species to [s1v1, s1v2, ...,s1vn, s2v1, ..., smvn] (where s=species and v=vertex). Currently, it is still [s1v1, s2v1, snv1, s1v2, ..., smvn].
- Enabling of indexing of the solution object will be done in a follow-up PR.
Edge parameter value generation functions
Since setting edge parameter values can be a fit fiddly, I have created two functions to help the user do this for (Cartesian and masked) grid lattices.
make_edge_p_values(lrs::LatticeReactionSystem, make_edge_p_value::Function)
enables the user to provide a function (make_edge_p_value
) which takes two arguments (src_vert
and dst_vert
) and from these determines a value along that edge. E.g. in the following example, we assign the value 0.1
to all edges, except for the one leading from vertex (1,1) to vertex (1,2), to which we assign the value 1.0
:
using Catalyst
rn = @reaction_network begin
(p,d), 0 <--> X
end
tr = @transport_reaction D X
lattice = CartesianGrid((5,5))
lrs = LatticeReactionSystem(rn, [tr], lattice)
function make_edge_p_value(src_vert, dst_vert)
if src_vert == (1,1) && dst_vert == (1,2)
return 1.0
else
return 0.1
end
end
D_vals = make_edge_p_values(lrs, make_edge_p_value)
make_directed_edge_values(lrs, ...)
allows the user to provide a number of two-values Tuple's (equal to the grid's number of dimensions). These tuples determines the edges values according to that dimension and that direction. E.g. in the following example, we wish to have diffusion in the x dimension, but a constant flow from low y values to high y values (so not transportation from high to low y). We achieve it in the following manner:
using Catalyst
rn = @reaction_network begin
(p,d), 0 <--> X
end
tr = @transport_reaction D X
lattice = CartesianGrid((5,5))
lrs = LatticeReactionSystem(rn, [tr], lattice)
D_vals = make_directed_edge_values(lrs, (0.1, 0.1), (0.1, 0.0))
Here, since we have a 2d grid, we only provide the first two Tuples to make_directed_edge_values
.
Change which order with which variables/vertexes are stored (so new order is first all variables at site 1, then all variables at site 2, and so on).
Have you verified what order MOL uses?
Still waiting for Alex to confirm, think he is probably on vacation. Holding off one actually changing anything until I've gotten in contact with him. If we need to chang the order I will do so here.
Ok, ordering os variables/compartments doesn't really seem to be thing in MoL, so I will skip that step here.
This PR is ready now. There is one error that I only get when I run CI tests locally/here (but now when I run it in a normal environment). I wouldn't be surprised if it is related to the new MTK changes that are lying around. Since https://github.com/SciML/Catalyst.jl/pull/737 should be merged before this one anyway, I will look into it further after that.