-
-
Notifications
You must be signed in to change notification settings - Fork 34
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #219 from Vilin97/spatially-varying-rx-rates-WIP
working on spatially varying reaction rates, WIP
- Loading branch information
Showing
10 changed files
with
248 additions
and
10 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,71 @@ | ||
| struct SpatialMassActionJump{A<:Union{AbstractVector,Nothing},B<:Union{AbstractMatrix,Nothing},S,U,V} <: AbstractMassActionJump | ||
| uniform_rates::A # reactions that are uniform in space | ||
| spatial_rates::B # reactions whose rate depends on the site | ||
| reactant_stoch::S | ||
| net_stoch::U | ||
| param_mapper::V | ||
|
|
||
| """ | ||
| uniform rates go first in ordering | ||
| """ | ||
| function SpatialMassActionJump{A,B,S,U,V}(uniform_rates::A, spatial_rates::B, reactant_stoch::S, net_stoch::U, param_mapper::V, scale_rates::Bool, useiszero::Bool, nocopy::Bool) where {A<:Union{AbstractVector,Nothing},B<:Union{AbstractMatrix,Nothing}, S, U, V} | ||
| uniform_rates = (nocopy || isnothing(uniform_rates)) ? uniform_rates : copy(uniform_rates) | ||
| spatial_rates = (nocopy || isnothing(spatial_rates)) ? spatial_rates : copy(spatial_rates) | ||
| reactant_stoch = nocopy ? reactant_stoch : copy(reactant_stoch) | ||
| for i in eachindex(reactant_stoch) | ||
| if useiszero && (length(reactant_stoch[i]) == 1) && iszero(reactant_stoch[i][1][1]) | ||
| reactant_stoch[i] = typeof(reactant_stoch[i])() | ||
| end | ||
| end | ||
| num_unif_rates = isnothing(uniform_rates) ? 0 : length(uniform_rates) | ||
| if scale_rates && num_unif_rates > 0 | ||
| scalerates!(uniform_rates, reactant_stoch) | ||
| end | ||
| if scale_rates && !isnothing(spatial_rates) && !isempty(spatial_rates) | ||
| scalerates!(spatial_rates, reactant_stoch[num_unif_rates+1:end]) | ||
| end | ||
| new(uniform_rates, spatial_rates, reactant_stoch, net_stoch, param_mapper) | ||
| end | ||
|
|
||
| end | ||
|
|
||
| ################ Constructors ################## | ||
|
|
||
| SpatialMassActionJump(urates::A, srates::B, rs::S, ns::U, pmapper::V; scale_rates = true, useiszero = true, nocopy=false) where {A<:Union{AbstractVector,Nothing},B<:Union{AbstractMatrix,Nothing},S,U,V} = SpatialMassActionJump{A,B,S,U,V}(urates, srates, rs, ns, pmapper, scale_rates, useiszero, nocopy) | ||
| SpatialMassActionJump(urates::A, srates::B, rs, ns; scale_rates = true, useiszero = true, nocopy=false) where {A<:Union{AbstractVector,Nothing},B<:Union{AbstractMatrix,Nothing}} = SpatialMassActionJump(urates, srates, rs, ns, nothing; scale_rates=scale_rates, useiszero=useiszero, nocopy=nocopy) | ||
|
|
||
| SpatialMassActionJump(srates::B, rs, ns, pmapper; scale_rates = true, useiszero = true, nocopy=false) where {B<:Union{AbstractMatrix,Nothing}} = SpatialMassActionJump(nothing, srates, rs, ns, pmapper; scale_rates = scale_rates, useiszero = useiszero, nocopy=nocopy) | ||
| SpatialMassActionJump(srates::B, rs, ns; scale_rates = true, useiszero = true, nocopy=false) where {B<:Union{AbstractMatrix,Nothing}} = SpatialMassActionJump(nothing, srates, rs, ns, nothing; scale_rates = scale_rates, useiszero = useiszero, nocopy=nocopy) | ||
|
|
||
| SpatialMassActionJump(urates::A, rs, ns, pmapper; scale_rates = true, useiszero = true, nocopy=false) where {A<:Union{AbstractVector,Nothing}} = SpatialMassActionJump(urates, nothing, rs, ns, pmapper; scale_rates = scale_rates, useiszero = useiszero, nocopy=nocopy) | ||
| SpatialMassActionJump(urates::A, rs, ns; scale_rates = true, useiszero = true, nocopy=false) where {A<:Union{AbstractVector,Nothing}} = SpatialMassActionJump(urates, nothing, rs, ns, nothing; scale_rates = scale_rates, useiszero = useiszero, nocopy=nocopy) | ||
|
|
||
| SpatialMassActionJump(ma_jumps::MassActionJump{T,S,U,V}; scale_rates = true, useiszero = true, nocopy=false) where {T,S,U,V} = SpatialMassActionJump(ma_jumps.scaled_rates, ma_jumps.reactant_stoch, ma_jumps.net_stoch, ma_jumps.param_mapper; scale_rates = scale_rates, useiszero = useiszero, nocopy=nocopy) | ||
|
|
||
| ############################################## | ||
|
|
||
| get_num_majumps(spatial_majump::SpatialMassActionJump{Nothing,Nothing,S,U,V}) where {S,U,V} = 0 | ||
| get_num_majumps(spatial_majump::SpatialMassActionJump{Nothing,B,S,U,V}) where {B,S,U,V} = size(spatial_majump.spatial_rates, 1) | ||
| get_num_majumps(spatial_majump::SpatialMassActionJump{A,Nothing,S,U,V}) where {A,S,U,V} = length(spatial_majump.uniform_rates) | ||
| get_num_majumps(spatial_majump::SpatialMassActionJump{A,B,S,U,V}) where {A<:AbstractVector,B<:AbstractMatrix,S,U,V} = length(spatial_majump.uniform_rates) + size(spatial_majump.spatial_rates, 1) | ||
| using_params(spatial_majump::SpatialMassActionJump) = false | ||
|
|
||
| rate_at_site(rx, site, spatial_majump::SpatialMassActionJump{Nothing,B,S,U,V}) where {B,S,U,V} = spatial_majump.spatial_rates[rx, site] | ||
| rate_at_site(rx, site, spatial_majump::SpatialMassActionJump{A,Nothing,S,U,V}) where {A,S,U,V} = spatial_majump.uniform_rates[rx] | ||
| function rate_at_site(rx, site, spatial_majump::SpatialMassActionJump{A,B,S,U,V}) where {A<:AbstractVector,B<:AbstractMatrix,S,U,V} | ||
| num_unif_rxs = length(spatial_majump.uniform_rates) | ||
| rx <= num_unif_rxs ? spatial_majump.uniform_rates[rx] : spatial_majump.spatial_rates[rx-num_unif_rxs, site] | ||
| end | ||
|
|
||
| function evalrxrate(speciesmat::AbstractMatrix{T}, rxidx::S, majump::SpatialMassActionJump, site::Int) where {T,S} | ||
| val = one(T) | ||
| @inbounds for specstoch in majump.reactant_stoch[rxidx] | ||
| specpop = speciesmat[specstoch[1], site] | ||
| val *= specpop | ||
| @inbounds for k = 2:specstoch[2] | ||
| specpop -= one(specpop) | ||
| val *= specpop | ||
| end | ||
| end | ||
| @inbounds return val * rate_at_site(rxidx, site, majump) | ||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,130 @@ | ||
| using DiffEqJump, DiffEqBase, OrdinaryDiffEq | ||
| using Test, Graphs, LinearAlgebra | ||
|
|
||
| reltol = 0.05 | ||
| Nsims = 10^4 | ||
|
|
||
| dim = 1 | ||
| linear_size = 5 | ||
| dims = Tuple(repeat([linear_size], dim)) | ||
| num_nodes = prod(dims) | ||
| center_site = trunc(Int,(linear_size^dim + 1)/2) | ||
| u0 = zeros(Int, 1, num_nodes) | ||
| end_time = 100.0 | ||
| diffusivity = 1.0 | ||
| death_rate = 0.1 | ||
|
|
||
| num_species = 1 | ||
| reactstoch = [Pair{Int64, Int64}[], [1 => 1]] | ||
| netstoch = [[1 => 1],[1 => -1]] | ||
| uniform_rates = ones(2, num_nodes) | ||
| uniform_rates[2,:] *= death_rate | ||
| non_uniform_rates = zeros(2, num_nodes) | ||
| non_uniform_rates[:,center_site] = uniform_rates[:,center_site] | ||
|
|
||
| # DiscreteProblem setup | ||
| tspan = (0.0, end_time) | ||
| prob = DiscreteProblem(u0, tspan, uniform_rates) | ||
|
|
||
| # spatial system | ||
| grid = Graphs.grid(dims) | ||
| hopping_constants = [diffusivity for i in u0] | ||
|
|
||
| # majumps | ||
| uniform_majumps_1 = SpatialMassActionJump(uniform_rates[:,1], reactstoch, netstoch) | ||
| uniform_majumps_2 = SpatialMassActionJump(uniform_rates, reactstoch, netstoch) | ||
| uniform_majumps_3 = SpatialMassActionJump([1.], reshape(uniform_rates[2,:], 1, num_nodes), reactstoch, netstoch) # hybrid | ||
| uniform_majumps_4 = SpatialMassActionJump(MassActionJump(uniform_rates[:,1], reactstoch, netstoch)) | ||
| uniform_majumps = [uniform_majumps_1, uniform_majumps_2, uniform_majumps_3, uniform_majumps_4] | ||
|
|
||
| non_uniform_majumps_1 = SpatialMassActionJump(non_uniform_rates, reactstoch, netstoch) # reactions are zero outside of center site | ||
| non_uniform_majumps_2 = SpatialMassActionJump([1.], reshape(non_uniform_rates[2,:], 1, num_nodes), reactstoch,netstoch) # birth everywhere, death only at center site | ||
| non_uniform_majumps_3 = SpatialMassActionJump([1. 0. 0. 0. 0.; 0. 0. 0. 0. death_rate], reactstoch,netstoch) # birth on the left, death on the right | ||
| non_uniform_majumps = [non_uniform_majumps_1, non_uniform_majumps_2, non_uniform_majumps_3] | ||
|
|
||
| # put together the JumpProblem's | ||
| uniform_jump_problems = JumpProblem[JumpProblem(prob, NSM(), majump, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(false,false)) for majump in uniform_majumps] | ||
| # flattenned | ||
| append!(uniform_jump_problems, JumpProblem[JumpProblem(prob, NRM(), majump, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(false,false)) for majump in uniform_majumps]) | ||
|
|
||
| # non-uniform | ||
| non_uniform_jump_problems = JumpProblem[JumpProblem(prob, NSM(), majump, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(false,false)) for majump in non_uniform_majumps] | ||
|
|
||
|
|
||
| # testing | ||
| function get_mean_end_state(jump_prob, Nsims) | ||
| end_state = zeros(size(jump_prob.prob.u0)) | ||
| for i in 1:Nsims | ||
| sol = solve(jump_prob, SSAStepper()) | ||
| end_state .+= sol.u[end] | ||
| end | ||
| end_state/Nsims | ||
| end | ||
|
|
||
| function discrete_laplacian_from_spatial_system(spatial_system, hopping_rate) | ||
| sites = 1:num_sites(spatial_system) | ||
| laplacian = zeros(length(sites), length(sites)) | ||
| for site in sites | ||
| laplacian[site,site] = -outdegree(spatial_system, site) | ||
| for nb in neighbors(spatial_system, site) | ||
| laplacian[site, nb] = 1 | ||
| end | ||
| end | ||
| laplacian .*= hopping_rate | ||
| laplacian | ||
| end | ||
| L = discrete_laplacian_from_spatial_system(grid, diffusivity) | ||
|
|
||
| # birth and death everywhere | ||
| f(u,p,t) = L*u - death_rate*u + uniform_rates[1,:] | ||
| ode_prob = ODEProblem(f, zeros(num_nodes), tspan) | ||
| sol = solve(ode_prob, Tsit5()) | ||
|
|
||
| for spatial_jump_prob in uniform_jump_problems | ||
| solution = solve(spatial_jump_prob, SSAStepper()) | ||
| mean_end_state = get_mean_end_state(spatial_jump_prob, Nsims) | ||
| mean_end_state = reshape(mean_end_state, num_nodes) | ||
| diff = mean_end_state - sol.u[end] | ||
| for (i,d) in enumerate(diff) | ||
| @test abs(d) < reltol*sol.u[end][i] | ||
| end | ||
| end | ||
|
|
||
| # birth and death zero outside of center site | ||
| f(u,p,t) = L*u - diagm([0., 0., death_rate, 0., 0.])*u + [0., 0., 1., 0., 0.] | ||
| ode_prob = ODEProblem(f, zeros(num_nodes), tspan) | ||
| sol = solve(ode_prob, Tsit5()) | ||
|
|
||
| solution = solve(non_uniform_jump_problems[1], SSAStepper()) | ||
| mean_end_state = get_mean_end_state(non_uniform_jump_problems[1], Nsims) | ||
| mean_end_state = reshape(mean_end_state, num_nodes) | ||
| diff = mean_end_state - sol.u[end] | ||
| for (i,d) in enumerate(diff) | ||
| @test abs(d) < reltol*sol.u[end][i] | ||
| end | ||
|
|
||
| # birth everywhere, death only at center site | ||
| f(u,p,t) = L*u - diagm([0., 0., death_rate, 0., 0.])*u + ones(num_nodes) | ||
| ode_prob = ODEProblem(f, zeros(num_nodes), tspan) | ||
| sol = solve(ode_prob, Tsit5()) | ||
|
|
||
| solution = solve(non_uniform_jump_problems[2], SSAStepper()) | ||
| mean_end_state = get_mean_end_state(non_uniform_jump_problems[2], Nsims) | ||
| mean_end_state = reshape(mean_end_state, num_nodes) | ||
| diff = mean_end_state - sol.u[end] | ||
| for (i,d) in enumerate(diff) | ||
| @test abs(d) < reltol*sol.u[end][i] | ||
| end | ||
|
|
||
| # birth on left end, death on right end | ||
| f(u,p,t) = L*u - diagm([0., 0., 0., 0., death_rate])*u + [1., 0., 0., 0., 0.] | ||
| ode_prob = ODEProblem(f, zeros(num_nodes), tspan) | ||
| sol = solve(ode_prob, Tsit5()) | ||
|
|
||
| solution = solve(non_uniform_jump_problems[3], SSAStepper()) | ||
| mean_end_state = get_mean_end_state(non_uniform_jump_problems[3], Nsims) | ||
| mean_end_state = reshape(mean_end_state, num_nodes) | ||
| diff = mean_end_state - sol.u[end] | ||
| for (i,d) in enumerate(diff) | ||
| @test abs(d) < reltol*sol.u[end][i] | ||
| end |