Comments (14)
In addition, the different model spec given in the docs throws an error:
using DifferentialEquations
prob = DiscreteProblem([999,1,0],(0.0,250.0))
rate = (t,u) -> (0.1/1000.0)*u[1]*u[2]
affect! = function (integrator)
integrator.u[1] -= 1
integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!)
rate = (t,u) -> 0.01*u[2]
affect! = function (integrator)
integrator.u[2] -= 1
integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,affect!)
jump_prob = JumpProblem(prob,jump,jump2)
sol = solve(jump_prob,Discrete(apply_map=false))
MethodError: no method matching DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::DiffEqBase.DiscreteProblem{Array{Int64,1},Float64,true,DiffEqBase.##130#131}, ::DiffEqJump.ConstantRateJump{##17#18,##19#20}, ::DiffEqJump.ConstantRateJump{##21#22,##23#24})
Closest candidates are:
DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::Any, ::DiffEqJump.AbstractAggregatorAlgorithm, ::DiffEqJump.ConstantRateJump{F1,F2}; kwargs...) at /home/simon/juliapro/JuliaPro-0.5.1.1/JuliaPro/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:9
DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::Any, ::DiffEqJump.AbstractAggregatorAlgorithm, ::Any...; kwargs...) at /home/simon/juliapro/JuliaPro-0.5.1.1/JuliaPro/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:11
DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::P, ::A, ::J<:Union{DiffEqJump.AbstractJumpAggregator,Void}, ::C, ::J2) at /home/simon/juliapro/JuliaPro-0.5.1.1/JuliaPro/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:2
from jumpprocesses.jl.
Second problem, that needs a jump aggregation method:
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
For the first problem, the first set of random numbers should be good. But, hopefully this version of the loop clarifies what's up:
nums = Int[]
@time for i in 1:1000
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,Discrete(apply_map=false))
push!(nums,sol[end][3])
end
mean(nums)
(Note: I checked the means against your Gillespie.jl, so unless something has changed, it should end up with the same distribution)
Since the first jump is only calculated when the jump_prob
is constructed, a funny thing is happening. First of all, if it It's using the last states of the previous run, causing it to want to die out more. But secondly, it's not going to "redo" the 1st-jump calculation unless you call jump_prob = JumpProblem(prob,Direct(),jump,jump2)
. So in order to test different runs, you need to do:
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,Discrete(apply_map=false))
You need to run both of those lines! If you do that, it should work.
I don't know when this got introduced, but this is indeed awful behavior. This is because the jump's code isn't ever touched until it reaches the first timepoint of the jump, and so there isn't a way to input "setup the first jump again". But, if there was an initialization phase for callbacks:
then this would be fixed.
from jumpprocesses.jl.
And thank you for testing this out! I would like to hear some feedback from you. This setup using the DiscreteCallbacks
in the ODE solver is more flexible because it can be attached to any DE, but it is still not as efficient as doing the straight loop. I think I got like 4x from Gillespie.jl last time I checked? The only-jumps version might need to be pulled out and not use the ODE solver logic in order to get more efficient. But I am also wondering if the fact that I am using tuples of tuples is making a difference here. For small equations, tuples of tuples to store the reaction rates probably slows things down considerably over using a dense matrix like Gillespie.jl, but should have an advantage in large sparse models. I don't have a good way of testing how much of a difference is directly due to that choice though.
from jumpprocesses.jl.
For reference, 590cf3f pulled out the initialization from the solve
command. But that should go directly into the DE solvers. I'll see if I can get this done right now.
from jumpprocesses.jl.
Adding in the initialization step worked, but the latest build throws an error.
using DiffEqJump, DiffEqBase, OrdinaryDiffEq
using Base.Test
rate = (t,u) -> (0.1/1000.0)*u[1]*u[2]
affect! = function (integrator)
integrator.u[1] -= 1
integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!;save_positions=(false,true))
rate = (t,u) -> 0.01u[2]
affect! = function (integrator)
integrator.u[2] -= 1
integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,affect!;save_positions=(false,true))
prob = DiscreteProblem([999.0,1.0,0.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,Discrete(apply_map=false))
nums = Int[]
@time for i in 1:1000
sol = solve(jump_prob,Discrete(apply_map=false))
push!(nums,sol[end][1])
end
mean(nums)
MethodError: no method matching DiffEqBase.DiscreteCallback{F1,F2,F3}(::DiffEqJump.DirectJumpAggregation{Float64,Tuple{##1#2,##5#6},Tuple{##3#4,##7#8}}, ::DiffEqJump.DirectJumpAggregation{Float64,Tuple{##1#2,##5#6},Tuple{##3#4,##7#8}}, ::Tuple{Bool,Bool})
Closest candidates are:
DiffEqBase.DiscreteCallback{F1,F2,F3}{F1,F2,F3}(::F1, ::F2, ::F3, ::Tuple{Bool,Bool}) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqBase/src/callbacks.jl:45
DiffEqBase.DiscreteCallback{F1,F2,F3}(::Any, ::Any; initialize, save_positions) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqBase/src/callbacks.jl:50
DiffEqBase.DiscreteCallback{F1,F2,F3}(::Any, ::Any, ::DiffEqJump.ConstantRateJump{F1,F2}, ::Any) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqJump/src/callbacks.jl:31
...
in #JumpProblem#10(::Tuple{Bool,Bool}, ::Type{T}, ::DiffEqBase.DiscreteProblem{Array{Float64,1},Float64,true,DiffEqBase.##130#131}, ::DiffEqJump.Direct, ::DiffEqJump.JumpSet{Tuple{},Tuple{DiffEqJump.ConstantRateJump{##1#2,##3#4},DiffEqJump.ConstantRateJump{##5#6,##7#8}}}) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:23
in #JumpProblem#9(::Array{Any,1}, ::Type{T}, ::DiffEqBase.DiscreteProblem{Array{Float64,1},Float64,true,DiffEqBase.##130#131}, ::DiffEqJump.Direct, ::DiffEqJump.ConstantRateJump{##1#2,##3#4}, ::Vararg{Any,N}) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:11
in DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2}(::DiffEqBase.DiscreteProblem{Array{Float64,1},Float64,true,DiffEqBase.##130#131}, ::DiffEqJump.Direct, ::DiffEqJump.ConstantRateJump{##1#2,##3#4}, ::DiffEqJump.ConstantRateJump{##5#6,##7#8}) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/v0.5/DiffEqJump/src/problem.jl:11
from jumpprocesses.jl.
In addition, with the old commit (I get the same error as above with the latest build) I get different results using the Reaction syntax.
using DifferentialEquations
infection = Reaction(0.1/1000,[1,2],[(1,-1),(2,1)])
recovery = Reaction(0.01,[2],[(2,-1),(3,1)])
srand(1234)
nums = Int[]
@time for i in 1:1000
sir_prob = DiscreteProblem([999,1,0],(0.0,250.0))
sir_jump_prob = GillespieProblem(sir_prob,Direct(),infection,recovery)
sir_sol = solve(sir_jump_prob,Discrete())
push!(nums,sir_sol[end][3])
end
mean(nums)
from jumpprocesses.jl.
Right now that master branch is:
Pkg.checkout("OrdinaryDiffEq")
Pkg.checkout("DiffEqJump")
(along with Stochastic and Delay DiffEq, all for the initialize change, along with retcodes). I'll be pushing tags when tests pass.
My test case is:
using DifferentialEquations
infection = Reaction(0.1/1000,[1,2],[(1,-1),(2,1)])
recovery = Reaction(0.01,[2],[(2,-1),(3,1)])
sir_prob = DiscreteProblem([999,1,0],(0.0,250.0))
sir_jump_prob = GillespieProblem(sir_prob,Direct(),infection,recovery)
sir_sol = solve(sir_jump_prob,Discrete())
using Plots; plot(sir_sol)
srand(1234)
nums = Int[]
@time for i in 1:100000
sir_sol = solve(sir_jump_prob,Discrete())
push!(nums,sir_sol[end][3])
end
println("Reaction DSL: $(mean(nums))")
using DiffEqJump, DiffEqBase, OrdinaryDiffEq
using Base.Test
rate = (t,u) -> (0.1/1000.0)*u[1]*u[2]
affect! = function (integrator)
integrator.u[1] -= 1
integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!;save_positions=(false,true))
rate = (t,u) -> 0.01u[2]
affect! = function (integrator)
integrator.u[2] -= 1
integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,affect!;save_positions=(false,true))
prob = DiscreteProblem([999.0,1.0,0.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,Discrete(apply_map=false))
using Plots; plot(sol)
nums = Int[]
@time for i in 1:100000
sol = solve(jump_prob,Discrete(apply_map=false))
push!(nums,sol[end][3])
end
println("Direct Jumps: $(mean(nums))")
using Gillespie
function F(x,parms)
(S,I,R) = x
(beta,gamma) = parms
infection = beta*S*I
recovery = gamma*I
[infection,recovery]
end
x0 = [999,1,0]
nu = [[-1 1 0];[0 -1 1]]
parms = [0.1/1000.0,0.01]
tf = 250.0
srand(1234)
nums = Int[]
@time for i in 1:100000
result = ssa(x0,F,nu,parms,tf)
data = ssa_data(result)
push!(nums,data[:x3][end])
end
println("Gillespie: $(mean(nums))")
and I get for the means:
Reaction DSL: 725.36465
Direct Jumps: 725.03103
Gillepsie: 725.53742
For the record, the times are:
Reaction DSL: 137.620577 seconds (1.97 G allocations: 87.976 GB, 23.05% gc time)
Direct Jumps: 116.954541 seconds (1.81 G allocations: 63.778 GB, 16.63% gc time)
Gillespie: 25.244438 seconds (1.48 G allocations: 57.974 GB, 34.18% gc time)
so there seems to be a performance regression somewhere changing 2x-4x to 5x-6x. Performance will be another issue though. A lot of it might have to do with the splatting penalty, which seems to keep growing:
The easiest way around all of this would be generated functions, but in theory there's a way without them
from jumpprocesses.jl.
Dear @ChrisRackauckas
Interesting. Looking forward to the tagged releases; still, it's good to know that my implementation gets something for being less friendly. May I suggest using a fixed random number seed and testing on the basis of the output (which is what I do in Gillespie.jl to avoid messing things up)?
from jumpprocesses.jl.
May I suggest using a fixed random number seed and testing on the basis of the output (which is what I do in Gillespie.jl to avoid messing things up)?
Yeah, that is something a little finicky too. Juno actually uses some random numbers when setting up its progress bar, so many times when I am benchmarking I get subtle differences when I go to a large problem and turn the progress bar on.
Interesting. Looking forward to the tagged releases; still, it's good to know that my implementation gets something for being less friendly.
Your implementation is much easier to check the correctness of, and it's designed as a straight loop and so it's quite easy to see that it's "nearly optimal". I'm really glad you made it, because it's the perfect test case!
I think that in the end we can have something more extensive here, but that also means that the barrier to entry for "how this works" is much higher, and much more testing has to be done to make it robust.
from jumpprocesses.jl.
I think all of the releases are tagged now. Those test case should run.
from jumpprocesses.jl.
Dear @ChrisRackauckas
Looks OK now...thanks!
from jumpprocesses.jl.
Cool. Double checking: what performance numbers do you get? Similar to what I posted above?
from jumpprocesses.jl.
Yes, pretty similar - 3-4x performance hit on a Mac Pro too.
from jumpprocesses.jl.
Good to know. Sounds like a good universal test case to work on. Thanks for your input. This should make the whole jump much easier to use.
from jumpprocesses.jl.
Related Issues (20)
- formatting needs to be rerun
- `save_positions` not working with callbacks
- `Coevolve` can call `urate` with a time beyond the time span HOT 1
- Check if wrapping affects to always return nothing has perf impact
- Efficency of ExtendedJumpArray broadcasting in ode_interpolant HOT 10
- Improve docstrings
- Simulating jump problems require non-integer tspan values HOT 3
- Interface for Jump simulations seems non-ideally messy? HOT 7
- possible interpolation bug in hybrid systems HOT 1
- Including global variables that participate in MassActionJumps in spatial systems HOT 10
- Remarks about docs [JuliaCon Proceedings review] HOT 4
- Remarks about paper [JuliaCon Proceedings review] HOT 1
- Remarks about package [JuliaCon Proceedings review] HOT 7
- Dynamical ODE and SDE solvers fail when they are combined with VariableRateJump. HOT 1
- Method ambiguities reported by Aqua HOT 2
- mutating parameter vector needs to call reset_aggregated_jumps! HOT 7
- Reenable formatter once SciMLStyle is fixed.
- Certain Continuous Callbacks does not work for JumpProcessess HOT 14
- Remarks about paper - episode 2 [JuliaCon Proceedings review] HOT 5
- Add / use cache-optimized indexed priority queue in NRM-based methods. HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from jumpprocesses.jl.