Best practices for defining stochastic affect! functions for Jump Problems (RNG)

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 JumpProcesses.jl/JumpProcesses.jl at 23aa6ac0f78b05cff1fa573a2b8d3f06ff34d585 · SciML/JumpProcesses.jl · GitHub, and set the default value of the rng keyword to it, see JumpProcesses.jl/problem.jl at 23aa6ac0f78b05cff1fa573a2b8d3f06ff34d585 · SciML/JumpProcesses.jl · GitHub, 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.