Comments (14)
Also if want to check with closures method, here is the code:
function test_closures()
p = (μ = 0.01, σ = 0.1, N = 300) # if all constant
T = 10.0 # maximum time length
x_iv = rand(p.N) # just draws from the inital condition
prob = SDEProblem(μ_SDE, σ_SDE, x_iv ,(0.0, T), p)
rate(u,p,t) = 0.2
affect_index!(integrator, index) = (integrator.u[index] = max(integrator.u[index], integrator.u[rand(1:integrator.p.N)]) )
jumps = [ConstantRateJump(rate,AffectIndex(affect_index!, i)) for i in 1:p.N]
jump_prob = JumpProblem(prob,Direct(),jumps...)
@btime solve($jump_prob);
end
test_closures()
from jumpprocesses.jl.
Thanks, I’ll play around with this when I’m back in Boston next week.
from jumpprocesses.jl.
Thanks. We will investigate other approaches for when the jumps are all of identical frequency (as discussed before), so the main purpose for this would be to have a framework to modify for when the jumps have different rates
from jumpprocesses.jl.
But also, we just aren't sure what to consider swapping out for the Direct()
and SRIW1()
in general. Our understanding is that changjng Direct to DirectFW in the closure example is the best choice, but otherwise we don't know permutations of either.
from jumpprocesses.jl.
Direct()
and SRIW1
are completely different things and cannot be swapped. SSAStepper is a DiscreteProblem-only fast stepper which is an alternative to SRIW1 if there's no differential equation. Direct is an aggregator, which would need to be swapped with a different jump aggregator.
from jumpprocesses.jl.
OK. So it sounds like the following is true: Direct -> DirectFW if we are wrapping a whole bunch of things of different types, otherwise it shouldn't matter.
And SRIW1 is the only alternative with this setup (since we aren't a pure jump process).
from jumpprocesses.jl.
SOSRI is a bit better.
from jumpprocesses.jl.
@jlperla Just an update. I've been playing around with this, and I am seeing that even with save_positions=(false,false)
and saving at a small number of time points the performance is pretty slow (N=300
). Unfortunately the current brokenness of Cairo.jl means that I haven't been able to actually visualize any profiling data to see what is going on.
Did you or @PooyaFa ever profile the code and see where it is actually spending its time? (i.e. is the major amount of time actually being spent in Direct
, or is it in another part of the DiffEq library code?)
from jumpprocesses.jl.
Thanks! I think based on what we discussed with Chris, tau leaping seems like the only way to really jack up the number of points. However, even without the jump size issues, the current regular jump code only does pure jump processes. Alas, we need jump diffusions.
from jumpprocesses.jl.
I don't think tau-leaping is needed, and it can potentially introduce other issues. One should certainly be able to simulate thousands of jump processes coupled to diffusion processes without needing tau-leaping. The problem is determining where the bottleneck in the current code base is. It could be in Direct
, but it could also be in how mixed jump-diffusion systems are being handled by the timestepper. Profiling, if Cairo.jl ever gets fixed, would help in determining where the bottleneck is occurring.
from jumpprocesses.jl.
If I've understood your problem / code correctly, there are no "dependencies" between the jumps. That is, when a given jump occurs it does not change u
in such a way that we would need to recalculate the rate for other jumps (since none of the jump rates actually depend on u
). In that case NRM
should be faster and work well.
Could you try this code out, which uses NRM
, and see if it is actually giving results that look accurate to you?
function test_SRIW1()
p = (μ = 0.01, σ = 0.1, N = 500) # if all constant
T = 10.0 # maximum time length
x_iv = rand(p.N) # just draws from the inital condition
prob = SDEProblem(μ_SDE, σ_SDE, x_iv ,(0.0, T), p)
rate(u,p,t) = 0.2
affect_index!(integrator, index) = (integrator.u[index] = max(integrator.u[index], integrator.u[rand(1:integrator.p.N)]) )
jumps = [ConstantRateJump(rate,AffectIndex(affect_index!, i)) for i in 1:p.N]
jump_prob = JumpProblem(prob,NRM(),JumpSet((), jumps, nothing, nothing),save_positions=(false,false),dep_graph=[Int[] for i=1:p.N])
sol = solve(jump_prob, SRIW1(), saveat=T);
#@btime begin solve($jump_prob, $SRIW1(), saveat=$T); end;
return sol;
end
from jumpprocesses.jl.
Note that it currently only saves at the final time point, but you could add more save points using saveat
.
from jumpprocesses.jl.
Also, did you ever actually try DirectFW
with save_positions=(false,false)
and specifying some saveat
? I'm finding just this is much faster on my computer. Direct
is 12.241 s, DirectFW
is 91.373 ms and NRM is 74.197 ms
.
from jumpprocesses.jl.
I am going to close this out for now, as I think that the combination of not saving everything and using DirectFW (or NRM) closes this particular issue.
from jumpprocesses.jl.
Related Issues (20)
- 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
- figure out why CI is now taking a ridiculous amount of time HOT 2
- Uniformization Direct Methods
- network free solvers
- Better documentation of `remake` mutating original problems
- Discrete events with repeated dosage fails for JumpSystems (fine ofr e.g. ODEs) HOT 7
- Automatic selection of solver method (e.g. `SSAStepper()`) HOT 2
- (Hierarchical) Exact Accelerated SSA HOT 1
- Seeding ensemble problems HOT 8
- `FunctionMap` not defined HOT 7
- add option to convert MassActionJumps to ConstantRateJumps HOT 3
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.