Interface for Multivariate processes by JoseKling · Pull Request #92 · JoseKling/PointProcesses.jl

JoseKling · 2026-06-02T20:50:45Z

This PR has the objective of implementing multivariate Hawkes processes.
In the process, some other modifications were carried out.

Split implementation of univariate Hawkes process into three files and added mark distributions
Defined separate common interfaces for AbstractUnivariateProcess and AbstractMultivariateProcess

- Added mark distribution to the process

- Added `univariate_process.jl` and `multivariate_process.jl` files defining the common interfaces - Updated the multivariate processes to comply with the interface - Split univariate Hawkes implementation into thre files

codecov · 2026-06-02T20:54:29Z

Codecov Report

❌ Patch coverage is 94.90741% with 11 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/bounded_point_process.jl	44.44%	5 Missing ⚠️
src/multivariate_process.jl	88.57%	4 Missing ⚠️
src/history.jl	94.28%	2 Missing ⚠️

📢 Thoughts on this report? Let us know!

JoseKling · 2026-06-21T16:33:13Z

I figured it would be a good idea to split this PR in two.
This part only fixes the univariate Hawkes process and defines the interface for multivariate processes.
The next PR will be the implementation of multivariate Hawkes processes.

rsenne · 2026-06-27T01:25:41Z

Is this ready for me to review?

JoseKling · 2026-06-27T06:15:20Z

Sorry. I guess I forgot to request the review. But it should be ready, yes

rsenne · 2026-06-27T21:22:07Z

+    descendants = T[]
+    next_gen = immigrants
+    while !isempty(next_gen)
+        # OPTIMIZE: Can this be improved by avoiding allocations of `curr_gen` and `next_gen`? Or does the compiler take care of that?


probably. we would want to make a workspace type that we use inplace operations on

rsenne · 2026-06-27T23:16:06Z

+"""
+abstract type AbstractMultivariateProcess <: AbstractPointProcess end
+
+Base.ndims(pp::AbstractMultivariateProcess) = length(pp.mark_dist)


This assumes there is a mark_dist field but IndependentMultivariate does not have one?

Yes, IndependentMultivariateProcess is a bit different, since we are not really building a new process, but only gluing a bunch of processes together.
For new multivariate processes, I would expect we use the interface for mark distributions. The idea is that the mark distributions at time t are independent conditional on the history up to t, so we can model them separetely (well, technically not really true, only for processes in which the probability of simultaneous events is 0, which are all interesting processes I am aware of)

rsenne · 2026-06-27T23:21:16Z

@@ -0,0 +1,102 @@
+"""
+    StatsAPI.fit(rng, ::Type{HawkesProcess{T}}, h::History; step_tol::Float64 = 1e-6, max_iter::Int = 1000) where {T<:Real}


stale docstring

rsenne · 2026-06-27T23:22:20Z

    end

-    pp_est = fit(PP, h)
+    pp_est = fit(PP, h)::AbstractPointProcess # JET complains if not specified


this is bein caused by the type stability issue i'm p sure

This is the message I get.

┌ MonteCarloTest(S::Type{<:Statistic}, PP::Type{<:AbstractPointProcess}, h::History; kwargs...) @ PointProcesses /home/jkling/projects/PointProcesses.jl/src/HypothesisTests/PPTests/monte_carlo_test.jl:114
│ no matching method found kwcall(::NamedTuple, ::Type{MonteCarloTest}, ::Type{<:Statistic}, ::Rician, ::History) (1/4 union split): Core.kwcall(merge(NamedTuple()::@NamedTuple{}, kwargs::@kwargs{…})::NamedTuple, MonteCarloTest, S::Type{<:Statistic}, pp_est::Union{Rician, HawkesProcess, IndependentMultivariateProcess, PoissonProcess}, h::History)
└────────────────────

Apparently, the static analysis thinks that Rician is a possible return type of fit (Rician is a Distributions.Distribution). Not sure how to deal with this if not by type annotating this line.

rsenne · 2026-06-27T23:22:49Z


    # Estimate process and calculate statistic from data
-    pp_est = fit(PP, h; rng=rng)
+    pp_est = fit(PP, h; rng=rng)::AbstractPointProcess # JET complains if not specified


same as my other comment

Added tests for custom custom processes

rsenne · 2026-07-06T23:53:47Z

+end

-"""
-function StatsAPI.fit(


These fit methods e.g., fit(HawkesProcess, h) were delted yet still referenced in docs (see index.md, hawkes.md, and montecarlo.md)

rsenne · 2026-07-06T23:56:04Z

+    sort!(sim)
+    h = History(eltype(sim)[], tmin, tmax, eltype(hp.mark_dist)[], Nothing[], 1)
+    for event in sim
+        push!(h, event, sample_mark(hp, event, h))


Currently ignores the RNG

Suggested change

push!(h, event, sample_mark(hp, event, h))

push!(h, event, sample_mark(rng, hp, event, h))

rsenne · 2026-07-07T00:04:38Z

+    transformed_times = [zeros(T, nb_events(h, d)) for d in 1:ndims(h)]
+    for d in 1:ndims(pp)


A 3-dim process with a history whose keyword constructor computed N=length(unique(dims)) = 2, which split_into_chunks would do, would cause a BBounds error. I think changiing transformed tiems to index by 1:ndims(pp) would fix.

but should we allow for time_change to be applied with process and history with different dimensions? I think it is better for split_into_chunks to guarantee all the chunks have the same number of dimensions as the original history.

sure, that makes sense. In that case we probably want to thread N through:

chunk = History(times, a, b, marks, dims, h.N; check_args=false)

rsenne · 2026-07-07T00:06:17Z

+    return sum(log.(hp.μ .+ (hp.α .* A))) - # Value of intensity at each event
+           (hp.μ * duration(h)) - # Integral of base rate
+           ((hp.α / hp.ω) * sum(1 .- exp.(-hp.ω .* (duration(h) .- h.times)))) + # Integral of each kernel
+           sum(log.([densityof(hp.mark_dist, t, h, m) for (t, m) in zip(h.times, h.marks)])) # Loglikelihood of marks


Suggested change

sum(log.([densityof(hp.mark_dist, t, h, m) for (t, m) in zip(h.times, h.marks)])) # Loglikelihood of marks

sum(log(densityof(hp.mark_dist, t, h, m)) for (t, m) in zip(h.times, h.marks); init=0.0) # Loglikelihood of marks

rsenne · 2026-07-07T00:14:18Z

    return PoissonProcess(λ, mark_dist)
 end

-function StatsAPI.fit(pptype::Type{PoissonProcess{R,D}}, args...; kwargs...) where {R,D}


this deletion means we no longer support weighted fitting. is that intentional?

I think I was testing to see if this was the origin of the problem with JET and forgot to add back. Anyway, good catch

Even though this is not really in tune with the rest of the package anymore. Not sure if we should keep it, or try to implement something similar for the other processes, or get rid of it...

rsenne · 2026-07-07T00:16:27Z

+    end
+    return sum(log.(hp.μ .+ (hp.α .* A))) - # Value of intensity at each event
+           (hp.μ * duration(h)) - # Integral of base rate
+           ((hp.α / hp.ω) * sum(1 .- exp.(-hp.ω .* (duration(h) .- h.times)))) + # Integral of each kernel


should be h.tmax .- h.times ?

rsenne · 2026-07-07T00:19:12Z

+function ground_intensity(hp::HawkesProcess, t, h::History)
+    activation = sum(exp.(hp.ω .* (@view h.times[1:(searchsortedfirst(h.times, t) - 1)])))
+    return hp.μ + (hp.α * activation / exp(hp.ω * t))
+end


need the shift invariant form

rsenne · 2026-07-07T00:20:24Z

+            dims = fill(nothing, length(times))
+            N = 1
+        else
+            N = length(unique(dims))


This counts distinct dim values present, not the actual dimension count

yes. In the case we are passing vectors, instead of one vector per dimension, this is the only way to deduce the number of dimensions in the history.
I agree this has some weird edge cases. We should probably think in a more structured way about the interface for the constructor. Perhaps this is trying to be too general.

some ideas:

Deduce with N = maximum(dims) instead of length(unique(dims))

Additionally accept an explicit N/ndims keyword so callers like split_into_chunks can preserve the parent's dimension count

In the split_into_chunks I planned to do exactly this.
About deducing the dimensions with max(dims), this also has edge cases in which there are no events in the dimensions with higher values. Perhaps for multivariate histories the dimension must be passed, if the number of dimensions is nos passed we assume it is univariate?

rsenne · 2026-07-07T00:21:11Z

        any((μ, α, ω) .< 0) &&
            throw(DomainError((μ, α, ω), "All parameters must be non-negative."))
        (α > 0 && α >= ω) &&
            throw(DomainError((α, ω), "Parameter ω must be strictly smaller than α"))


error message inverted?

rsenne · 2026-07-07T00:23:18Z

found a few more bugs. i'll see if i can figure out the JET thing

- Fixed constructor for empty history - History does not deduce number of dimensions - Small fixes for Hawkes processes

JoseKling · 2026-07-18T10:42:32Z

I think I got everything, except for the JET warning.
I changed History slightly. Instead of guessing the number of dimensions, if a Vector is passed to the constructor, we assume it is a univariate history. The inner constructor does not have any default parameters, so everything must be passed explicitly.
Also, I changed the constructor for the empty history to expect a type for the marks and the number of dimensions (not necessary for univariate history)

rsenne · 2026-07-21T19:25:03Z

i'll try to look at this soon

JoseKling added 2 commits May 28, 2026 19:44

Univariate Hawkes process

e98f3e6

- Added mark distribution to the process

Interface for Multivariate PointProcesses

7f45b10

- Added `univariate_process.jl` and `multivariate_process.jl` files defining the common interfaces - Updated the multivariate processes to comply with the interface - Split univariate Hawkes implementation into thre files

JoseKling mentioned this pull request Jun 2, 2026

Marked hawkes #64

Closed

github-actions Bot reviewed Jun 2, 2026

View reviewed changes

JoseKling added 3 commits June 20, 2026 10:48

Merge branch 'main' into MultivariateHawkes

28c32f6

Automatic formatting

c2d4ba1

Fixed convergence in example CustomProcesses

af47ed7

JoseKling changed the title ~~Multivariate hawkes~~ Interface for Multivariate processes Jun 21, 2026

JoseKling marked this pull request as ready for review June 21, 2026 16:33

rsenne self-requested a review June 27, 2026 20:59

rsenne requested changes Jun 27, 2026

View reviewed changes

JoseKling force-pushed the MultivariateHawkes branch from c3eec46 to 774c595 Compare July 2, 2026 16:06

PR fixes

4052588

Added tests for custom custom processes

JoseKling force-pushed the MultivariateHawkes branch from 774c595 to 4052588 Compare July 5, 2026 16:00

JoseKling requested a review from rsenne July 5, 2026 16:08

rsenne requested changes Jul 7, 2026

View reviewed changes

Fix for PR

5d08bc9

- Fixed constructor for empty history - History does not deduce number of dimensions - Small fixes for Hawkes processes

JoseKling force-pushed the MultivariateHawkes branch from 1328b5e to 5d08bc9 Compare July 18, 2026 10:33

JoseKling requested a review from rsenne July 18, 2026 10:38

		@@ -0,0 +1,102 @@
		"""
		StatsAPI.fit(rng, ::Type{HawkesProcess{T}}, h::History; step_tol::Float64 = 1e-6, max_iter::Int = 1000) where {T<:Real}

	push!(h, event, sample_mark(hp, event, h))
	push!(h, event, sample_mark(rng, hp, event, h))

		transformed_times = [zeros(T, nb_events(h, d)) for d in 1:ndims(h)]
		for d in 1:ndims(pp)

	sum(log.([densityof(hp.mark_dist, t, h, m) for (t, m) in zip(h.times, h.marks)])) # Loglikelihood of marks
	sum(log(densityof(hp.mark_dist, t, h, m)) for (t, m) in zip(h.times, h.marks); init=0.0) # Loglikelihood of marks

Conversation

JoseKling commented Jun 2, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

codecov Bot commented Jun 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

JoseKling commented Jun 21, 2026

Uh oh!

rsenne commented Jun 27, 2026

Uh oh!

JoseKling commented Jun 27, 2026

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rsenne commented Jul 7, 2026

Uh oh!

JoseKling commented Jul 18, 2026

codecov Bot commented Jun 2, 2026 •

edited

Loading