I have a JumpProcesses.JumpProblem that I am stepping through piecewise (using integrator = init(...)/step!(integrator, dt, true)) and I would like to save states only during some of these intervals according to a fixed but complex schedule (depending only on t).
It seems like the following approach works correctly:
setting integrator.save_everystep = true right before step!ing intervals where I wish to save everything, and resetting it to false right afterwards
savevalues!(integrator, true) when I want to save only at specific time points after stepping there exactly using step!(integrator, ..., true)
Is there any reason why I shouldn’t be doing that? (This seems to indicate so.)
If so, how would I implement this correctly? Would I need to define multiple callbacks and switch them out manually by setting integrator.cb? (I cannot change save_positions since DiscreteCallback is immutable.)
I wouldn’t recommend changing save_everystep by hand. Create your problem by passing save_positions = (false,false) i.e. in the call to JumpProblem(...; save_positions = (false,false)).
By setting save_positions on the Problem I am controlling the initial saving behavior; my question is about how I can change that after simulating for some dt. For example:
integrator = init(jump_problem, SSAStepper(), save_start = false)
integrator.save_everystep = false # cannot set this on init
step!(integrator, 1.0, true) # don't save from 0 to 1
integrator.save_everystep = true
step!(integrator, 1.0, true) # save all changes from 1 to 2
integrator.save_everystep = false
step!(integrator, 1.0, true) # don't save from 2 to 3
savevalues!(integrator, true) # save state at 3
solve!(integrator) # solve to end, save final state
This seems to work, but it feels like I am using it incorrectly.
My understanding is that I would need to define one callback for each time interval that I want different saving behavior for, and then switch it out on the integrator – is that correct? And how would I do that, just by setting integrator.cb?
Just to clarify, do you want to save at discrete time points that you know / specify, or do you want to save the exact state (i.e. at all jump times) over some sub-intervals?
Both: I have a list of (non-overlapping, ordered, potentially zero-width) time intervals and I want to record all state changes (new states at all jumps) within the intervals, as well as the current state at the end of each of these intervals. (If an interval is zero-width, that means only saving the state at that time point.) I don’t want to save anything outside of these intervals.
(The reason is that I have another Gillespie implementation that I am configuring in this way, and now I am trying to build a SciML-based replacement which I would prefer to behave the same initially.)
I guess another option is to make a DiscreteCallback where the condition is true whenever you are at the desired discrete times or within the desired time intervals. You can set tstops in solve for the discrete time points, and for the saving interval endpoints, to make sure you step to them exactly. Then your affect! function just saves the state using savevalues!. That seems like it would do what you want (if I’m understanding correctly).
That said, what you are with the integrator is a reasonable approach too (though I think of save_everystep is internal, but as your problem shows perhaps this shouldn’t be the case since mutating it might be needed for some models).
Using a DiscreteCallback like this would imply checking the intersections at every step (though I guess that could be implemented in O(1) amortized by keeping a watermark).
I’ll keep using the current solution described above until I start hitting problems. Otherwise maybe another solution could be to schedule PresetTimeCallbacks or similar, and then manipulate integrator.save_everystep from there – that should be allowed according to the affect! documentation for DiscreteCallback (or at least that’s how I would interpret it).
PresetTimeCallbacks just generate DiscreteCallbacks, so I think they are still getting called every time step and checking if the current time is in their time list. But yes, if checking if the current time is a saving time is expensive relative to one-step of executing the simulation, hand-using the integrator interface will likely be better.