DifferentialEquations.jl: marked point processes (or storing user data)

Specific Domains Numerics
1

I’m trying to simulate a simple marked point process with DifferentialEquations. As far as I can tell, this is not directly supported (the documentation does mention a mark in the context of RegularJump, but this seems to be a feature specific to Gillespie-type models).
I want to be able to extract jump time and marks from the solution/integrator.
One possible, albeit somewhat backward, solution is to store the last time and mark directly in the state of the underlying ODE/SDE problem, and extract them after solving with unique. So a simple example would look like this:

time_and_mark(t) = (t, randn())
jump!(integrator) = integrator.u .= time_and_mark(integrator.t)
ode = ODEProblem((du,u,p,t)->(), [0.0,0.0], (0.0, 1.0))
prob = JumpProblem( ode, Direct(), ConstantRateJump( (u,p,t)->10.0, jump! ) )
sol = solve(prob,Tsit5())
times, marks = [[x[i] for x in u] for i=1:2 for u=(unique(x->x[1],sol.u)[2:end],)]

It get even more awkward when the marks are not scalar. I was hoping for a more elegant solution.

The docs mention that the integrator interface includes a p field for “user-provided data”, and that it can be initialized in init, but I couldn’t get it to work, for example:

jump!(integrator) = (@show integrator.p)
ode = ODEProblem((du,u,p,t)->(), [0.0], (0.0, 1.0))
prob = JumpProblem( ode, Direct(), ConstantRateJump( (u,p,t)->5.0, jump! ) )
integrator = init(prob, Tsit5(), p = Float64[])
sol = solve!(integrator)

shows that integrator.p is nothing rather than the expected empty vector. Trying to store anything into it inside jump! results in error since it tries to convert it nothing.

Any advice on how to store user data with the integrator or another solution for modeling marked point processes would be appreciated.

2

Your process seems quite simple. So why don’t you simulate it directly?

3

The code examples are not the actual process I’m trying to simulate. They were simplified to illustrate my problem without unnecessary detail.
I encountered this problem a few times with different mixed models (diffusion process + marked point process), so I would prefer a general solution. So far I worked around the problem by putting the times and jumps into the diffusion process as illustrated in the first example.

4

Hey,
Yes, I’ve somewhat ignored marks since I’m not really sure what they are used for and wanted to learn more before trying to make them formal or document an example with them. Is it just counting the number of times a given jump has occurred? And vector when there are multiple types of jumps? How are you using this in a model? I’m curious because I haven’t had to make use of it and want to know what you need with them in order to get it implemented :slight_smile:.

As for doing this in non-regular jumps, yes using the parameters is probably a good way to do it. Or use closures. Either way is fine. For parameters, they go on the ODEProblem: it’s an optional last argument. ODEProblem(f,u0,tspan,p). So your example is

time_and_mark(t) = (t, randn())
jump!(integrator) = integrator.u .= time_and_mark(integrator.t)
ode = ODEProblem((du,u,p,t)->(), [0.0,0.0], (0.0, 1.0),Float64[])
prob = JumpProblem( ode, Direct(), ConstantRateJump( (u,p,t)->10.0, jump! ) )
sol = solve(prob,Tsit5())

and then that empty array of Float64’s shows up as integrator.p. Then of course you can put whatever you want in there and modify it as you please.

5

Somehow I didn’t make the connection between this p argument and the problem parameters, and thought this is some additional feature. The documentation at http://docs.juliadiffeq.org/latest/basics/integrator.html#Handing-Integrators-1 is a bit confusing.
I guess I’ll use parameters for now, though it feels like a bit of a hack to use them to store extra data generated while solving. I suppose it’s also a bit error-prone since they won’t be automatically initialized in solve, so some care will be needed with the Monte Carlo interface, but I might just skip it and manage the collection of statistics myself.
Incidentally, the docs also mention “user data” here http://docs.juliadiffeq.org/stable/basics/common_solver_opts.html#User-Data-1, but it’s not very clear how to use it. Is it meant for something of this sort (storing extra data generated while solving)?

The context is modeling a large population of sensory neurons. Neurons emit spikes at random times influenced by external stimulus and possibly influencing other brain processes (which have a more abstracted model so they live in the SDE part). There may be additional information for each spike, which is where marks come in. In my case the mark attached to each spike has parameters of the spiking neuron (you can also formulate this as an unmarked process for each neuron rather than a single marked process, however the marked process approach allows some simplification by approximating a large population using a continuous distribution over neuron parameters). I imagine other models of spiking neural population might also have marks specific to spikes (not just neurons).
I can’t say much unfortunately about how others use MPP models, but I expect support for equations like this should cover most cases, without requiring any new algorithms:
du_t = f(u_t,p,t)dt + Σgᵢ(u_t,p,t)dWⁱ + \int_y h(u_t,p,t,y)N(dt,dy)
where the last term is an integral w.r.p the MPP N viewed as a random discrete measure over [t_0,t_1]\times\mathcal{Y}, with \mathcal{Y} being the space of marks. There could also be a sum over such integrals with different processes N^i, though strictly speaking it doesn’t increase the generality.
Mathematically, a marked point process is characterized by an “intensity kernel” \lambda_t(dy) (relative to some relevant filtration \mathcal{F}_t, which for this purpose could just be “all the history”), such that the rate of points with marks in Y\subset\mathcal{Y} at time t (conditioned on \mathcal{F}_t) is \int_Y\lambda_t(dy). Practically, that just means you generate the times as a point process with rate \int_\mathcal{Y}\lambda_t(dy) and sample each mark from \lambda_t(dy)/\int_\mathcal{Y}\lambda_t(dy). So I imagine an interface like this

# total rate (∫λ(u,p,t,dy), rate of the unmarked process)
total_rate(u,p,t) = ...
# generate a random mark for a point at time t with state u
# (sample from (λ(u,p,t,dy) / ∫λ(u,p,t,dy))
sample_mark(u,p,t) = ...
# affect the state u following a point at time t with mark m
# (the effect of ∫h(...)N(dt,dy))
point_affect!(du,u,p,t,m) = ...

mpp = MarkedPointProcess( total_rate, sample_mark, point_affect! )
prob = JumpProblem( ode_or_sde, mpp )
sol = solve(prob)
times, marks = get_point_events(sol)

which would just mean something like

function jump!(du,u,p,t)
    m = sample_mark(u,p,t)
    store_time_and_mark_somewhere(t,m,somewhere)
    point_affect!(du,u,p,t,m)
end
jump = VariableRateJump( total_rate, jump! )
prob = JumpProblem( ode_or_sde, jump )

(Not sure about the names. There’s a bit of a discrepancy between the “marked points” and “jumps” terminology.)

6

Thanks for the explanation! Yes, I came across marks in Platen’s book on jump diffusions but haven’t gotten to doing anything with them. Could you file an issue in DiffEqJump.jl? We should definitely get this supported in the main jump definitions, likely using Distributions.jl to define a mark distribution.

Edit: No, Distributions.jl wouldn’t do well because it should be a function of parameters. You’re right: allowing users to pass a mark sampler function is the key.

7

Finally got round to it: