Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] LatticeReactionSystem internal update #756

Merged
merged 50 commits into from
Jul 9, 2024
Merged

Conversation

TorkelE
Copy link
Member

@TorkelE TorkelE commented Dec 31, 2023

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.

@isaacsas
Copy link
Member

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?

@TorkelE
Copy link
Member Author

TorkelE commented Jan 1, 2024

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
Copy link
Member Author

TorkelE commented Jan 18, 2024

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

@TorkelE TorkelE changed the title LatticeReactionSystem internal update [WIP] LatticeReactionSystem internal update Jan 26, 2024
@TorkelE TorkelE changed the base branch from master to Catalyst_version_14 February 7, 2024 01:07
@TorkelE TorkelE force-pushed the spatial_internal_update branch from f0fab91 to afef75e Compare February 7, 2024 13:36
@TorkelE TorkelE changed the title [WIP] LatticeReactionSystem internal update [v14 ready] LatticeReactionSystem internal update Feb 7, 2024
@TorkelE
Copy link
Member Author

TorkelE commented Feb 8, 2024

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 #737 should be merged before this one anyway, I will look into it further after that.

Base automatically changed from Catalyst_version_14 to master May 15, 2024 17:37
@TorkelE TorkelE changed the title [v14 ready] LatticeReactionSystem internal update [WIP] LatticeReactionSystem internal update Jun 11, 2024
@TorkelE TorkelE merged commit 2228734 into master Jul 9, 2024
5 of 6 checks passed
@TorkelE TorkelE deleted the spatial_internal_update branch July 9, 2024 22:24
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants