SciML / JumpProcesses.jl

Build and simulate jump equations like Gillespie simulations and jump diffusions with constant and state-dependent rates and mix with differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/JumpProcesses/stable/
Other
136 stars 35 forks source link

Setfields for jump problems #310

Closed TorkelE closed 1 year ago

TorkelE commented 1 year ago

Two things:

enables setdinex! and getindex for `JumpProblems. E.g.

using Catalyst, JumpProcesses
rn = @reaction_network begin
    (p,d), 0 <--> X
end

u0 = [:X => 0]
p = [:p => 1.0, :d => 0.1]
dp = DiscreteProblem(rn,u0,(0.0,10.0),p)
jp = JumpProblem(rn,dp,Direct())

jp[:X]
jp[:p]
jp[:p] = 1.5
jp[:X] = 15

It also adds a check so that if setindex! is called on a JumpProblem then update_parameters!(prob.massaction_jump, prob.prob.p) is also called

TorkelE commented 1 year ago

I would add some additional tests, but that would require adding MTK as a test dependency. Would that be preferred?

coveralls commented 1 year ago

Pull Request Test Coverage Report for Build 4612486657

Warning: This coverage report may be inaccurate.

This pull request's base commit is no longer the HEAD commit of its target branch. This means it includes changes from outside the original pull request, including, potentially, unrelated coverage changes.

Details


Changes Missing Coverage Covered Lines Changed/Added Lines %
src/problem.jl 2 6 33.33%
<!-- Total: 2 6 33.33% -->
Totals Coverage Status
Change from base Build 4549721708: 1.1%
Covered Lines: 2114
Relevant Lines: 2397

💛 - Coveralls
isaacsas commented 1 year ago

I don't think we want to have a test dependency on ModelingToolkit here, the dependencies / tests for this should be comparable to what is done in SciMLBase I'd think.

isaacsas commented 1 year ago

Maybe just add a simple test that you use to check this is working by directly building a JumpProblem instead of going through MTK?

TorkelE commented 1 year ago

I'd avoid the explicit import of the SciMLBase function, but otherwise looks good.

This is a good point

TorkelE commented 1 year ago

Maybe just add a simple test that you use to check this is working by directly building a JumpProblem instead of going through MTK?

Will try to make some simple test not requiring MTK

TorkelE commented 1 year ago

Actually, is that possible?

This is only for the cases when a prob has a f which has a sys. However, is there some non-MTK based way of creating systems? Even in SciMLBase we have to make downstream tests that depend on these.

This is only relevant when you have symbols like :X and :p associated with variables/parameters, but I don't think this is ever the case in JumpProcesses before e.g. MTK gets involved? I will look at it a bit more, but not sure if I will be able to figure it out.

isaacsas commented 1 year ago
using JumpProcesses
rate1(u,p,t) = 1.0
rate2(u,p,t) = 1.0
affect1!(integ) = (integ.u[1] +=1)
affect2!(integ) = (integ.u[2] += 1)
crj1 = ConstantRateJump(rate1, affect1!)
crj2 = ConstantRateJump(rate2, affect2!)
g = DiscreteFunction((du,u,p,t) -> nothing; syms = [:a, :b])
dprob = DiscreteProblem(g, [0, 0], (0.0, 10.0))
jprob = JumpProblem(dprob, Direct(), crj1, crj2)
sol = solve(jprob, SSAStepper())

Can be indexed with :a and :b:

julia> sol[:a]
22-element Vector{Int64}:
 0
 1
 1
 1
 2
 2
 2
 3
 3
 3
 3
 3
 4
 5
 6
 6
 7
 8
 8
 8
 8
 8
TorkelE commented 1 year ago

Thanks a lot! Will use this