JumpProblem hopping_constants callback?

There’s a new feature in JumpProcesses.jl (formerly DiffEqJump.jl) where you can model spatial jump processes like in this example post. One of the inputs is a hopping_constants matrix which controls how long a species will occupy a point in space. My question is: Can the hopping_constants matrix be modified in a callback?

Having different values in the hopping_constants matrix can mimic a boundary where the rate of diffusion differs on either side. By modifying the hopping_constants matrix you could simulate that boundary changing over time.

Here’s an example of some code I would like to add such a callback to:

using DifferentialEquations, Graphs, Plots, Printf

dims = (50,50)
numberOfSpecies = 1
numberOfNodes = prod(dims) # number of sites
grid = Graphs.grid(dims) # planar grid graph

#Put 25 molecules at each site less than 10 units from grid center
center = LinearIndices(dims)[dims .÷2...]
startingState = zeros(Int,numberOfSpecies, numberOfNodes)
startingState[gdistances(grid, center) .≤ 10] .= 25

#Create a Discrete problem 
tspan = (0.0, 10.0)
rates = [1.1, 1.0] 
prob = DiscreteProblem(startingState, tspan, rates)

#Encode the 2 reactions (could use Catalyst.jl to do this automatically)
reactstoch = [[1 => 1],[1 => 1]]
netstoch = [[1 => 1],[1 => -1]]
majumps = MassActionJump(rates, reactstoch, netstoch)

#Faster diffusion outside boundary
hopConstants = ones(numberOfSpecies, numberOfNodes)
hopConstants[gdistances(grid, center) .> 10] .= 2.0

alg = DirectCRDirect()
jump_prob = JumpProblem(prob,
                        spatial_system = grid,
                        save_positions=(false, false))

sol = solve(jump_prob, SSAStepper(), saveat = tspan[2]/300)

Code to animate the solution

#maxium value obtained
maxu = mapreduce(maximum,max,sol.u)

anim = @animate for t in range(tspan..., 300)
    currTime = @sprintf "Time: %.2f" t
        reshape(sol(t), dims),
        clims = (0,maxu),
        framestyle = :box,

gif(anim, "output.gif", fps = 60)


You could modify it in a callback, but we don’t have a version of reset_aggregated_jumps! that will communicate it has been changed and reset the internal data structures.

If you wouldn’t mind opening an issue on JumpProcesses we could look into adding that.

Oh for sure, I’ll open up an issue here. Thanks!