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 99.40476% with 1 line in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/bounded_point_process.jl	80.00%	1 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:14:31Z

-    end
-    return l
-end
+DensityInterface.logdensityof


rsenne · 2026-06-27T21:18:53Z

+HawkesProcess(μ, α, ω) = HawkesProcess(μ, α, ω, NoMarks())

-∑_{i=1:n} pᵢᵢ = ∑_{i=1:n} (1 - ∑_{j < i} Dᵢⱼ) = N - D.
+function ground_intensity(hp::HawkesProcess, h::History, t)


This will throw on a Univariate hawkes no? as in univariate_process.jl the order is t, h?

if yes we should also add a test for a univariate hawkes intensity

The test is already there, but I switched the order in the tests too : P

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-27T21:23:32Z

+function time_change(h::History, pp::AbstractUnivariateProcess)
+    Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t)
+    return time_change(h, Λ)
+end


redundant with abstract_point_process.jl:126?

rsenne · 2026-06-27T21:24:21Z

+    h_temp = History(sim, tmin, tmax)
+    marks = [sample_mark(hp.mark_dist, t, h_temp) for t in event_times(h_temp)]
+    return History(; times=sim, tmin=tmin, tmax=tmax, marks=marks)


why build history twice?

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:01:51Z

+```
+The default method uses a loop over events combined with `integrated_ground_intensity`, but it should be reimplemented for specific processes if faster computation is possible.
+"""
+function DensityInterface.logdensityof(pp::AbstractUnivariateProcess, h::History)


With the new type hierarchy, BoundedPointProcess now has no logdensityofmethod i.e., this would fail:

using PointProcesses h = History([1.0, 2.0, 3.0], 0.0, 10.0) logdensityof(BoundedPointProcess(PoissonProcess(1.0), 0.0, 10.0), 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.

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?

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

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

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

		@@ -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!

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!

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!

Reviewers

Assignees

Labels

Projects

Milestone

Development

codecov Bot commented Jun 2, 2026 •

edited

Loading