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

Spatial constant rate jump #343

Draft
wants to merge 15 commits into
base: master
Choose a base branch
from
Draft

Conversation

Vilin97
Copy link
Contributor

@Vilin97 Vilin97 commented Sep 3, 2023

Add spatial constant rate jump support.

NB: I have not implemented flattening for constant rate jumps, because it's unclear how to flatten constant rate jumps.

Implementation details:
In order to use constant rate jumps, the user is expected to pass in a JumpSet and two dependency graphs (reaction-reaction and reaction-species). The syntax is as follows:

rate1(u,p,t,site) = ...
rate2(u,p,t,site) = ...
affect1!(integrator,site) = ...
affect2!(integrator,site) =...
crjumps = JumpSet([ConstantRateJump(rate1, affect1!), ConstantRateJump(rate2, affect2!)])
# define prob, hopping_constants, spatial_system, and dependency_graphs
jump_prob = JumpProblem(prob, NSM(), crjumps, hopping_constants = hopping_constants,
            spatial_system = spatial_system , save_positions = (false, false), dep_graph = dep_graph, jumptovars_map = jumptovars_map)

@Vilin97 Vilin97 changed the title [WIP] Spatial constant rate jump 2 Spatial constant rate jump 2 Sep 3, 2023
@Vilin97 Vilin97 marked this pull request as ready for review September 3, 2023 06:03
Conflicts:
	src/spatial/reaction_rates.jl
	test/spatial/reaction_rates.jl
@coveralls
Copy link

coveralls commented Sep 9, 2023

Pull Request Test Coverage Report for Build 6180688297

  • 31 of 38 (81.58%) changed or added relevant lines in 7 files are covered.
  • No unchanged relevant lines lost coverage.
  • Overall coverage increased (+0.1%) to 89.981%

Changes Missing Coverage Covered Lines Changed/Added Lines %
src/problem.jl 2 3 66.67%
src/spatial/directcrdirect.jl 3 6 50.0%
src/spatial/nsm.jl 3 6 50.0%
Totals Coverage Status
Change from base Build 6167199515: 0.1%
Covered Lines: 2326
Relevant Lines: 2585

💛 - Coveralls

@Vilin97 Vilin97 changed the title Spatial constant rate jump 2 Spatial constant rate jump Sep 12, 2023
Copy link
Member

@isaacsas isaacsas left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of passing the site into the affect and rate functions, we could store it in the aggregator and make an interface for users to access it via the integrator (which stores the aggregator). Something like get_current_site(integrator). Then this would work with old-style constant rate jump signatures. I don't know if we want that though.

The other big issue is I think there are going to be type stability issues without using FunctionWrappers on the rate and affect functions (since each function is its own type, every crj will be a different type).

@@ -39,6 +39,9 @@ function DirectCRDirectJumpAggregation(nj::SpatialJump{J}, njt::T, et::T, rx_rat

# a dependency graph is needed
if dep_graph === nothing
if length(rx_rates.cr_jumps) != 0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if length(rx_rates.cr_jumps) != 0
if !isempty(rx_rates.cr_jumps)

@@ -54,6 +57,9 @@ function DirectCRDirectJumpAggregation(nj::SpatialJump{J}, njt::T, et::T, rx_rat
end

if jumptovars_map === nothing
if length(rx_rates.cr_jumps) != 0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if length(rx_rates.cr_jumps) != 0
if !isempty(rx_rates.cr_jumps)

@@ -32,6 +32,9 @@ function NSMJumpAggregation(nj::SpatialJump{J}, njt::T, et::T, rx_rates::RX, hop

# a dependency graph is needed
if dep_graph === nothing
if length(rx_rates.cr_jumps) != 0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same as above

@@ -47,6 +50,9 @@ function NSMJumpAggregation(nj::SpatialJump{J}, njt::T, et::T, rx_rates::RX, hop
end

if jumptovars_map === nothing
if length(rx_rates.cr_jumps) != 0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same as above

"""

### spatial rx rates ###
struct RxRates{F, M}
struct RxRates{F, M, C}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For future consideration, do you know if storing the jumps in RxRates instead of within the aggregator directly results in an extra layer of indirection when accessing them (i.e. an extra pointer call)? If so, we might want to consider moving them up a level if that shows any performance benefit.

end
RxRates(num_sites::Int, ma_jumps::M) where {M<:AbstractMassActionJump} = RxRates(num_sites, ma_jumps, ConstantRateJump[])
RxRates(num_sites::Int, cr_jumps::C) where {C} = RxRates(num_sites, SpatialMassActionJump(), cr_jumps)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
RxRates(num_sites::Int, cr_jumps::C) where {C} = RxRates(num_sites, SpatialMassActionJump(), cr_jumps)
RxRates(num_sites::Int, cr_jumps::C) where {C} = RxRates(num_sites, SpatialMassActionJump(), cr_jumps)

I'd specify a type for C here or not include C at all. Generally a generic type like this would be suggested to just be dropped since it isn't being used (it may even give a warning in tests).

Comment on lines +69 to +70
cr_jump = rx_rates.cr_jumps[rx - get_num_majumps(ma_jumps)]
set_rx_rate_at_site!(rx_rates, site, rx, cr_jump.rate(u, integrator.p, integrator.t, site))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This isn't type stable is it? That is why we split the ConstantRateJump rates and affects in the non-spatial solvers and wrap them inside FunctionWrappers.

Comment on lines +96 to +97
cr_jump = rx_rates.cr_jumps[rx - get_num_majumps(rx_rates.ma_jumps)]
cr_jump.affect!(integrator, site)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This isn't type stable is it? That is why we split the ConstantRateJump rates and affects in the non-spatial solvers and wrap them inside FunctionWrappers.

@Vilin97 Vilin97 marked this pull request as draft January 25, 2024 22:02
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.

3 participants