Avoiding Plotting to dense plots, saveat not working

Specific Domains Modelling & Simulations
1

I am running simulations of SDEs and plotting the results. I keep having problem that Jupyter crashes, despite that the simulations shouldn’t be that intense. I think hat I have tracked the problems to the plots being to big. I have earlier had this issue and got help with that. However the use of saveat (as suggested there) only seems to get me that far.

I have been investigating and it seems like the plots still are much larger than saveat indicates. Here I have tried to make a minimal(ish) example:

using Plots
gr();
using DifferentialEquations

function positive_domain()
    condition(u,t,integrator) = (any(u .< 0))
    affect!(integrator) = integrator.u .= integrator.uprev
    return DiscreteCallback(condition,affect!) 
end;

rn = @reaction_network rnType begin
    0.01, (X) → ∅
end 
prob_sde = SDEProblem(rn,[1.],(0.,2000.))
sol = solve(prob_sde,ImplicitEM(),dt=0.001,callback=PositiveDomain(),saveat=1.);
length(sol)

using saveat=1. I would expect to get a solution about the size of 2000, in fact its length is 2002002, which is quite a lot. I can also change to plotly(), zoom in and confirm that a plot(sol) shows things at a very small scale. I have tried other saveat values, like 1 and 0.5, but success.

Am I using saveat wrong somehow? I am having stochastic simulations over a time of about 2000, but not really interested in anything of a timescale less than 1, is there a good way to only plot about 2000 timepoints of my solution?

2

You are using the default values for save_positions in the callback.

http://docs.juliadiffeq.org/latest/features/callback_functions.html#DiscreteCallback-1

Notice that it says that the defaults are (true,true), so it’ll trigger a save just before and after each callback call. Of course, this is because if you want the full solution then this is required for proper discontinuity handling. But if you don’t want those saves, set save_positions = (false,false).

3

Thank you, that works.

Reading further, both you and the documents give warnings:

required for proper discontinuity handling

For discontinuous changes like a modification to u to be handled correctly (without error), one should set save_positions=(true,true)

Does this mean that I might get faulty solutions? After all I do not actually want to do anything with the solver, only what it outputs. Intuitively what I want to do (receive/plot a solution taking up slightly less space) should not be a problem for accurately simulating the system?

4

It doesn’t matter for simulating the system. It matters for plotting though. If you don’t do this and you have an event that causes a discontinuity (say, adds 5 to every variable), then you might notice that the plot looks incorrect because instead of increasing vertically at that point it will do some weird connection between the previous time step and the value after the discontinuity. For SDEs it’s not bad because it just interpolates the line in there, but technically that line is incorrect because it should be linear to the point before the continuity, and then vertical at the continuity, then continue. So for plotting purposes, along with the post-solution interpolation, this needs discontinuity needs to be handled exactly. For intermediate interpolation (like saveat and event handling), the integrator handles this just fine. So you don’t need to worry about correctness, but you may want to be careful with the interpretation of results around “large enough” discontinuities.

Hope that helps.

P.S.

If you use the GR backend, you can use gr(fmt=:png) to make the plots into PNGs, then they can have millions of points without a problem in Jupyter notebooks. Of course, that doesn’t fix your original problem so I didn’t mention it, but that’s something to keep in your back pocket if you really do need to plot a million points.

5

Thank, I see. So if I do not expect any large sudden event than I should not be worried about the plotting (although I maybe should be worried about my conviction that that is the case, but that is another issue).

Thanks for the gr(fmt=:png) tips, that might come in handy at some point.

6

Hi Gaussia,

I am new to Julia and am struggling with the non-negative solution of a system of ODEs.

Could you tell me why you defined a function called positive_domain() when you are using a pre-defined callback PositiveDomain()?
And in your defined function, after you go back to the values of the previous step what would happen? I see you have fixed the step size.

Thanks

7

There are a few problems in the positivity enforcement done earlier in this thread. For one, it has the issues you discussed. But more importantly, it seems that the SDE itself is actually trying to go negative (in cases, of course it’s stochastic) which of course would mess with it. These issues were found later on and handled differently, so I wouldn’t use this as a guide for how to handle non-negativity handling. PositiveDomain() or setting isoutofdomain would be better.

8

Hey, Chris.

Thank you for the reply.
I have just posted a question that is troubling me with these implementations. I must be doing something very silly. Can you look over the post if you find time?

9

Just re-iterating what Chris said that there are several irregularities with the above solution and you should not not spend to much time thinking on what is going on there.