Extend RegularJump interface to allow flexible matrix-free jump sizes
The current interface for regular jumps is https://docs.juliadiffeq.org/latest/types/jump_types.html#Defining-a-Regular-Jump-1
where c(dc,u,p,t,mark) is the Stoichiometry matrix. If it simulates N jumps for a particular jump process, then it finds the jump sizes in this matrix and the process jumps that N times in the dt interface. For some applications, this has two issues: (1) it might be more efficient to have a matrix-free calculation instead of the c(...) interface; and (2) for some applications you may want the jump sizes to change if there are multiple arrivals.
For the implementation in the algorithm, I don't think this would end up very complicated. Basically, instead of https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/simple_regular_solve.jl#L42-L44
!rj.constant_c && c(dc,uprev,p,tprev,mark)
mul!(update,dc,counts)
u[i] = uprev .+ update
it would be something like
c!(u, counts, p, tprev, mark)
Or something like that, and it would modify u in place given the counts. The user could copy the u themselves if they need to. The constant_c flag is not really necessary if it matrix-free, though it could be passed for the user to deal wtih.