Representation of multivariate processes

[JoseKling]

Following PR #51, multivariate processes are represented as marked processes with i.i.d. integer marks, each mark representing one marginal process. This is fine for Poisson processes, since the marks do not influence the intensity of the process.

For Hawkes processes, however, marks may influence the intensity, so a univariate Hawkes process with integer marks is not equivalent to a multivariate Hawkes process. Currently I am treating any integer marks as really meaning the marginal process an event occurs, but I think this is not great.

This was shortly discussed in #43. I am not currently considering processes with vector-valued marks, so the options I see are:

1. Define that this package deals only with univariate processes, but because of the structure of the Poisson process, we can interpret this univariate process as a multivariate one.

2. Fields times and marks of History hold vectors of vectors, instead of vectors.

3. Add another field that holds a vector containing the marginal process corresponding to each event. If this field contains a vector filled with nothing, than it is a univariate process.

4. Have multiple subtypes of History, which may or may not have the fields marks and marginal_process (or whatever name).

I would not consider (1), since I started working on mutually exciting Hawkes processes. Between (2) and (3), I guess (2) is a more straightforward way to represent multivariate processes, but (3) has the benefit that the event times are already ordered, which is needed for pretty much every method related with multivariate Hawkes processes. (4) is possibly the most efficient (?), but a similar functionality can be obtained using aliases.

Thoughts? Other alternatives?

[JoseKling]

What I came up with, which should not change the interface too much, so we can still use what we have with small modifications is this (without the complete implementation of the constructors, so it is not too long)

struct History{N,T<:Real,M,D} # N is the number of dimensions, D can be Int or Nothing
    times::Vector{T}
    tmin::T
    tmax::T
    marks::Vector{M}
    dims::Vector{D}     # In which of the 'N' marginal process the event occurs

    # Add optional arguments check
    # Inner constructors are the standard for the methods to create a `History`
    function History{N}(times::Vector{<:Real}, tmin::Real, tmax::Real, marks::Vector, dims::Vector{<:Int}) where {N}
        new{N,T,M,D}(times, tmin, tmax, marks, dims)
    end

    function History{1}(times::Vector{<:Real}, tmin::Real, tmax::Real, marks::Vector, ::Any)
        new{1,T,M,Nothing}(times, tmin, tmax, marks, dims)
    end
end

const UnivariateHistory{T<:Real,M} = History{1,T,M,Nothing}

# Outer constructors are useful for inputing data
function History(times::Vector{<:Vector{<:Real}}, tmin::Real, tmax::Real, marks::Vector{<:Vector}) 
    # From a vector of vectors to the representation used in the `History` struct
    return History{N}(times_vec, tmin, tmax, marks_vec, dims)
end

function History(times::Vector{<:Real}, tmin::Real, tmax::Real)
    History([times], tmin, tmax; check=check)
end

function History(times::Vector{<:Real}, tmin::Real, tmax::Real, marks::Vector)
    History([times], tmin, tmax, [marks]; check=check)
end

function History(times::Vector{<:Vector{<:Real}}, tmin::Real, tmax::Real)
    marks = [fill(nothing, length(times[i])) for i in 1:length(times)]
    History(times, tmin, tmax, marks; check=check)
end

For univariate processes, I guess nothing changes. The call to the constructor and the fields are the same
, and methods should all work (maybe with some minor adjustments for the added dims)

[rsenne]

This seems like a good approach to me. I like that it won'i interfere with the current implementation too much. If I understand correctly we could do something this for a bivariate process?

h = History{2}(times=[1.0, 1.5, 2.5, 3.7], tmin=0.0, tmax=10.0,
               marks=[m1, m2, m3, m4], dims=[1, 2, 1, 1])

[JoseKling]

Exactly. We have the two alternatives

History{2}(times=[1.0, 1.5, 2.5, 3.7], tmin=0.0, tmax=10.0, marks=[m1, m2, m3, m4], dims=[1, 2, 1, 1])
# or
History(times=[[1.0, 2.5, 3.7], [1.5]], tmin=0.0, tmax=10.0, marks=[[m1, m3, m4], [m2]])

And if times and marks are vectors instead of vectors of vectors in the second alternative, it constructs a univariate process.

[rsenne]

Yeah i think this is good i say move forward with it

[JoseKling]

So, this was done in PR #77. All tests pass without changing the implementation, since, as expected, this does not change anything for univariate processes.

The issue now is that this requires changing how processes are implemented. Instead of using the Categorical distribution in mark_dist to mock a multivariate process, we can make them actually multivariate, so we need to change the implementation.

I guess the following is the obvious answer. The logic of the methods won't change much, but I guess everything would need to be adjusted.

abstract type AbstractPointProcess{N} end

# Both univariate and multivariate processes hold a vector of parameters, even for N=1
struct PoissonProcess{N,R<:Real} <: AbstractPointProcess{N}
λ::Vector{R}
mark_dist::Vector{<:Distribution}
end

[rsenne]

hmm okay. What are the downstream ramifications you’re expecting if we switch to AbstractPointProcess{N}?

My worry is that to keep things type-stable we’d need to thread N through a lot of APIs (or enforce invariants everywhere at runtime), and that feels a bit cumbersome. Not saying that we shouldn't pursue--just thinking out loud/wondering how xpansive this would be

[jucheval]

Do we really need to pass N as a type parameter ?
If I remember correctly, one of @gdalle's arguments against this implementation is that JIT is applied each time you modify the dimension of your processes. I guess that it is related to @rsenne's comment too.

I think it would be possible to pass N as a field of the struct instead. Don't you think ?

[gdalle]

Definitely don't put N as a type parameter, you'll just strain the compiler needlessly

[JoseKling]

hmm okay. What are the downstream ramifications you’re expecting if we switch to AbstractPointProcess{N}?

My worry is that to keep things type-stable we’d need to thread N through a lot of APIs (or enforce invariants everywhere at runtime), and that feels a bit cumbersome. Not saying that we shouldn't pursue--just thinking out loud/wondering how xpansive this would be

I guess this will require changing a lot of stuff. As I mentioned, the logic will remain basically the same, but I would say all methods involving an AbstractPointProcess will need to be modified in some way.

Do we really need to pass N as a type parameter?
If I remember correctly, one of @gdalle's arguments against this implementation is that JIT is applied each time you modify the dimension of your processes. I guess that it is related to @rsenne's comment too.

I think it would be possible to pass N as a field of the struct instead. Don't you think?

That makes sense, but isn't that true for History too?
Of course we do not need to pass it as a type parameter, passing as a field is also fine. The only thing is that this enforces the history and the process to have the same dimensionality at compile time, instead of at runtime with manual checks.
And wouldn't this be an issue only for applications where there are many processes with different number of dimensions? I don't think this is usually the case?

[gdalle]

The only thing is that this enforces the history and the process to have the same dimensionality at compile time, instead of at runtime with manual checks.

You don't need to perform any manual checks at runtime. If you try to compute an intensity for a 3-dimensional Poisson process on a 2-dimensional history, you'll trigger a DimensionMismatch error, which is much more informative to users than the MethodError you'd get from a type parameter.

In general, specializing on the length of an array is a bad idea because it gives rise to as many methods as you have array lengths and creates type instabilities, mostly without any benefits. At least, with StaticArrays.jl, you get the benefits of unrolled operations making things faster, but even then, it is only recommended for very small arrays (length < 100).

[JoseKling]

Ok, then. I see.
I guess we can add a field N in History instead of adding it as a type parameter. This will be needed in case the last marginal history is empty and to avoid having to compute minimum(h.dims) and maximum(h.dims) to get the dimensionality.
For the processes we don't actually even need this N, then. We can just add the method n_dims (or something like this) that returns the length of the vector of parameters of the process.
Thanks for the feedback!