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

**URL:** <https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [January 9, 2019, 9:39pm UTC](https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453 "2019-01-09T21:39:22Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![yha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yha/32/3502_2.png) [@yha](https://discourse.julialang.org/u/yha)\
**Post date:** [January 9, 2019, 9:39pm UTC](https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453/1 "2019-01-09T21:39:22Z")

</div>

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:

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [January 9, 2019, 11:04pm UTC](https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453/2 "2019-01-09T23:04:05Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![yha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yha/32/3502_2.png) [@yha](https://discourse.julialang.org/u/yha)\
**Post date:** [January 9, 2019, 11:18pm UTC](https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453/3 "2019-01-09T23:18:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 10, 2019, 7:28am UTC](https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453/4 "2019-01-10T07:28:38Z")

</div>

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 🙂.

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![yha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yha/32/3502_2.png) [@yha](https://discourse.julialang.org/u/yha)\
**Post date:** [January 14, 2019, 1:24am UTC](https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453/5 "2019-01-14T01:24:39Z")

</div>

> [@ChrisRackauckas](#):
>
> yes using the parameters is probably a good way to do it

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](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](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)?

> [@ChrisRackauckas](#):
>
> 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 🙂.

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

```julia
# 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

```julia
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.)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 15, 2019, 6:06pm UTC](https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453/6 "2019-01-15T18:06:55Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![yha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yha/32/3502_2.png) [@yha](https://discourse.julialang.org/u/yha)\
**Post date:** [January 27, 2019, 4:27pm UTC](https://discourse.julialang.org/t/differentialequations-jl-marked-point-processes-or-storing-user-data/19453/7 "2019-01-27T16:27:19Z")

</div>

> [@ChrisRackauckas](#):
>
> Could you file an issue in DiffEqJump.jl?

Finally got round to it:

> <https://github.com/SciML/JumpProcesses.jl/issues/63>
>
> Following the discussion at https://discourse.julialang.org/t/differentialequati…ons-jl-marked-point-processes-or-storing-user-data/19453, I suggest adding support for marked point processes. Specifically, for simulating the process alongside an SDE affected by it, and retrieving times and marks from the solution.
> 
> As a suggested interface (slightly different than the one I suggested in the discourse thread above), the constructors for \`VariableRateJump\` and \`ConstantRateJump\` can have an additional optional (keyword?) argument \`mark\` specifying a mark sampling function. The \`affect!\` function would then take an additional argument for the mark.
> \`\`\`
> total\_rate(u,p,t) = ...
> sample\_mark(u,p,t) = ...
> affect!(du,u,p,t,y) = ...
> mpp = VariableRateJump( total\_rate, affect!, mark=sample\_mark )
> prob = JumpProblem( ode\_or\_sde, mpp )
> sol = solve(prob)
> times, marks = get\_point\_events(sol)
> \`\`\`
> which would mean something like
> \`\`\`
> function \_affect!(du,u,p,t)
> y = sample\_mark(u,p,t)
> store\_time\_and\_mark\_somewhere!(t,y,somewhere)
> affect!(du,u,p,t,y)
> end
> jump = VariableRateJump( total\_rate, \_affect! )
> prob = JumpProblem( ode\_or\_sde, jump )
> \`\`\`
> 
> Mathematically, this is enough to support equations of the form
> du\_t=\\\<drift and diffusion terms\\\> +∫\_y h(u\_t,p,t,y)N(dt,dy)
> where \`N\` is a MPP with an intensity kernel of the form λₜ(dy)=λ(u\_t,p,t,dy) (\`total\_rate\` is ∫λₜ(dy), \`sample\_mark\` samples \`y\` from λₜ(dy)/∫λₜ(dy))