Catalyst.jl icon indicating copy to clipboard operation
Catalyst.jl copied to clipboard

[WIP] LatticeReactionSystem internal update

Open TorkelE opened this issue 1 year ago • 4 comments

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 of Bool, (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:
    1. 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.
    2. 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.
    3. 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) and make_directed_edge_values.

Minor stuff

  • The internal view_vert_ps_vector! function is now called directly on LatticeTransportODEf/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 the ODEProblem 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.

TorkelE avatar Dec 31 '23 13:12 TorkelE

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?

isaacsas avatar Dec 31 '23 18:12 isaacsas

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.

TorkelE avatar Jan 01 '24 01:01 TorkelE

Ok, ordering os variables/compartments doesn't really seem to be thing in MoL, so I will skip that step here.

TorkelE avatar Jan 18 '24 00:01 TorkelE

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.

TorkelE avatar Feb 08 '24 15:02 TorkelE