I have created a small module based on DifferentialEquations to help me define a certain kind of stochastic differential equations where the stochastic processes are Lévy processes. As the jumps are stochastic, the kind of “affect functions” that act on the integrator object at jump times end up being similar to:
function affect!(integrator)
integrator.u[1] += some_function(rand())
return nothing
end
I currently make use of the rand function but I don’t even include Random in my module’s dependencies, nor do I declare the RNG anywhere. Should I do so? If so, what’s the cleanest way of doing this? I am assuming that it’s not a good idea to declare the RNG as a global variable akin to
module MyModule
import Random: Xoshiro
rng = Xoshiro()
function affect!(integrator)
integrator.u[1] += some_function(rand(rng))
return nothing
end
end
You could make it an optional parameter with a default value in whatever function serves as the starting point for using your package / setting up problems with it.
Note also that JumpProblems can take rng as a keyword argument, which sets the internal random number generator used for simulating the jump process. So you could configure your package to ensure they use the same rng. (Not sure if SDEProblems allow it to be set too.)
Don’t I need to a define a particular instance of the generator?
I thought I needed to do something akin to
rng = Xoshiro()
JumpProblem(..., rng=rng)
and that the excerpt below wouldn’t work:
JumpProblem(..., rng = Xoshiro())
Or do you mean instantiating it on startup? I use instantiating functions for a few minor configurations, but those are simple deterministic things whose purpose is to make finding files and data easier. I feel uneasy about having an init function export it as a global variable, and if I just define a function that gets read by every simulation function, it will always create a new instance of the rng everytime it’s called.
So you could configure your package to ensure they use the same rng.
By this do you mean passing the same rng on both calls? In other words:
rand(rng) and JumpProblem(..., rng=rng).
We just define a default generator instance in JumpProcesses, see https://github.com/SciML/JumpProcesses.jl/blob/23aa6ac0f78b05cff1fa573a2b8d3f06ff34d585/src/JumpProcesses.jl#L33, and set the default value of the rng keyword to it, see https://github.com/SciML/JumpProcesses.jl/blob/23aa6ac0f78b05cff1fa573a2b8d3f06ff34d585/src/problem.jl#L164, when calling JumpProblem.
So there is a module-wide default that is created and set to the default Julia generator, but it can be overridden if the user sets the keyword argument when calling JumpProblem.
By this do you mean passing the same rng on both calls? In other words:
rand(rng) and JumpProblem(..., rng=rng).
Yes.
If I understand correctly, there is a global variable defined within the module, then. Does this guarantee that the rng used for the jump part ends up being the same as the one used in the diffusion component (for example)?
To follow your suggestion of using the same rng for the affect function, would I need to import the cons from the package?
Something like
import DifferentialEquations.JumpProcesses: DEFAULT_RNG
function affect!(integrator)
integrator.u[1] += some_function(rand(DEFAULT_RNG))
return nothing
end
I wouldn’t suggest using JumpProcesses.DEFAULT_RNG as that is internal to the package. I was suggesting you could define an analogous default within your module, and an analogous way for users of your module to pass in an alternative rng. You could then pass whatever rng has been selected in your module down into JumpProcesses and other sub-libraries.
Unfortunately, I don’t think StochasticDiffEq currently let’s you change the rng in the default usage via SDEProblem, but @ChrisRackauckas would know better about that.
It does, you just have to define a noise process and then pass it as noise.
2 Likes
Alright, thank you both. I’ll try to see if it’s possible to refactor what I have using a custom noise process, and go from there.