Review Changes

.gitignore

@@ -3,4 +3,5 @@
 *.jl.mem
 /docs/build/
 **/Manifest.toml
-.vscode
+.vscode
+dev/

src/HypothesisTests/PPTests/monte_carlo_test.jl

@@ -110,6 +110,6 @@ function MonteCarloTest(
         throw(ArgumentError("Test is not valid for empty event history."))
     end
 
-    pp_est = fit(PP, h)::AbstractPointProcess # JET complains if not specified
+    pp_est = fit(PP, h)
     return MonteCarloTest(S, pp_est, h; kwargs...)
 end

src/abstract_point_process.jl

@@ -16,12 +16,13 @@ Base.ndims(::AbstractPointProcess)
 
 ## Intensity functions
 """
-    ground_intensity(pp, h, t)
+    ground_intensity(pp, t, h)
+    ground_intensity(pp, t, h, d)
 
 Compute the ground intensity for a temporal point process `pp` applied to history `h` at time `t`.
 For multivariate processes, it returns a vector of ground intensities for each dimension.
 
-ground_intensity(pp, h, t, d) computes the ground intensity for a multivariate process `pp` at dimension `d`.
+`ground_intensity(pp, t, h, d)` computes the ground intensity for a multivariate process `pp` at dimension `d`.
 
 The ground intensity quantifies the instantaneous risk of an event with any mark occurring at time `t` after history `h`:
 ```
@@ -32,47 +33,56 @@ function ground_intensity end
 
 """
     mark_distribution(pp, t, h)
+    mark_distribution(pp, t, h, d)
+    mark_distribution(md, t, h)
 
 Compute the distribution of marks for a temporal point process `pp` knowing that an event takes place at time `t` after history `h`.
 For multivariate processes, it returns a vector of mark distributions for each dimension.
 
-mark_distribution(pp, h, t, d) computes the mark distribution for a multivariate process `pp` at dimension `d`.
+`mark_distribution(pp, t, h, d)` computes the mark distribution for a multivariate process `pp` at dimension `d`.
+
+Given a mark distribution `md` instead of a process, `mark_distribution(md, t, h)` returns the
+distribution of marks at time `t` after history `h`. This is the method that a custom
+`AbstractMarkDistribution` must implement, and it must work for empty histories.
 """
 function mark_distribution end
 
 """
     intensity(pp, m, t, h)
+    intensity(pp, m, t, h, d)
 
 Compute the conditional intensity for a temporal point process `pp` applied to history `h` and event `(t, m)`.
 For multivariate processes, it returns a vector of intensities for each dimension.
 
-intensity(pp, h, t, d) computes the intensity for a multivariate process `pp` at dimension `d`.
+`intensity(pp, m, t, h, d)` computes the intensity for a multivariate process `pp` at dimension `d`.
 
 The conditional intensity function `λ(t,m|h)` quantifies the instantaneous risk of an event with mark `m` occurring at time `t` after history `h`.
 """
 function intensity end
 
 """
     log_intensity(pp, m, t, h)
+    log_intensity(pp, m, t, h, d)
 
 Compute the logarithm of the conditional intensity for a temporal point process `pp` applied to history `h` and event `(t, m)`.
 For multivariate processes, it returns a vector of log intensities for each dimension.
 
-log_intensity(pp, h, t, d) computes the log intensity for a multivariate process `pp` at dimension `d`.
+`log_intensity(pp, m, t, h, d)` computes the log intensity for a multivariate process `pp` at dimension `d`.
 """
 function log_intensity end
 
 ## Simulation
 
 """
     ground_intensity_bound(pp, t, h)
+    ground_intensity_bound(pp, t, h, d)
 
 Compute a local upper bound on the ground intensity for a temporal point process `pp` applied to history `h` at time `t`.
 
 Return a tuple of the form `(B, L)` satisfying `λg(t|h) ≤ B` for all `u ∈ [t, t+L)`.
 For multivariate processes, it returns a list of tuples [(B₁, L₁), (B₂, L₂), ...] for each dimension.
 
-ground_intensity_bound(pp, h, t, d) computes a local upper bound for a multivariate process `pp` at dimension `d`.
+`ground_intensity_bound(pp, t, h, d)` computes a local upper bound for a multivariate process `pp` at dimension `d`.
 """
 function ground_intensity_bound end
 
@@ -88,7 +98,7 @@ Compute the integrated ground intensity (or compensator) `Λ(t|h)` for a tempora
 
 For multivariate processes, it returns a vector of integrated ground intensities for each dimension.
 
-integrated_ground_intensity(pp, h, a, b, d) computes the integrated ground intensity for a multivariate process `pp` at dimension `d`.
+`integrated_ground_intensity(pp, h, a, b, d)` computes the integrated ground intensity for a multivariate process `pp` at dimension `d`.
 """
 function integrated_ground_intensity end
 

src/history.jl

@@ -14,6 +14,18 @@ Linear event histories with temporal locations of type `T` and marks of type `M`
 
 # Construction
 
+```julia
+History(times, tmin, tmax[, marks])                  # univariate, one vector of times
+History([times_1, ..., times_N], tmin, tmax[, marks])  # multivariate, one vector per dimension
+History(times, tmin, tmax, marks, dims, N)             # generic, everything explicit
+History(tmin, tmax, M[, N])                            # empty history with marks of type `M`
+History(; times, tmin, tmax, marks=nothing)            # keyword form of the first two
+```
+
+A plain vector of `times` always builds a univariate history (`N == 1`); the number
+of dimensions is never deduced from `dims`. To build a multivariate history, pass
+either one vector of times per dimension or `dims` together with an explicit `N`.
+
 The constructor validates its inputs and throws on any violation rather than
 silently coercing them. With `check_args=true` (the default), it requires:
 
@@ -99,7 +111,14 @@ struct History{T<:Real,M,D}
                 end
             end
             if N == 1
-                dims .= nothing
+                if !all(isnothing, dims)
+                    throw(
+                        DomainError(
+                            dims,
+                            "A univariate history (`N == 1`) must have `dims` filled with `nothing`.",
+                        ),
+                    )
+                end
             else
                 if any(d -> d < 1 || d > N, dims)
                     throw(
@@ -194,6 +213,13 @@ function History(tmin::R1, tmax::R2, M::Type, N::Int) where {R1,R2<:Real}
 end
 
 function History(h::History, d::Int)
+    if d < 1 || d > ndims(h)
+        throw(
+            DomainError(
+                d, "Dimension `d` must be between 1 and $(ndims(h)) for this history."
+            ),
+        )
+    end
     times = event_times(h, d)
     marks = event_marks(h, d)
     dims = fill(nothing, length(times))
@@ -227,7 +253,10 @@ event_times(h::History) = h.times
 Return the sorted vector of event times for `h` in dimension `d`.
 """
 function event_times(h::History, d::Union{Int,Nothing})
-    return h.N == 1 ? h.times : (@view h.times[h.dims .== d])
+    # For univariate histories `dims` is filled with `nothing`, so only `d == 1` (or
+    # an omitted `d`) refers to the events of this history. This condition must match
+    # the one in `event_marks(h, d)`, otherwise times and marks disagree in length.
+    return h.N == 1 && (isnothing(d) || d == 1) ? h.times : (@view h.times[h.dims .== d])
 end
 
 """
@@ -266,7 +295,7 @@ event_marks(h::History) = h.marks
 Return the vector of event marks in dimension `d` of `h`, sorted according to their event times.
 """
 function event_marks(h::History, d::Union{Int,Nothing})
-    return h.N == 1 && d == 1 ? h.marks : (@view h.marks[h.dims .== d])
+    return h.N == 1 && (isnothing(d) || d == 1) ? h.marks : (@view h.marks[h.dims .== d])
 end
 
 """
@@ -502,6 +531,11 @@ function Base.cat(h1::History, h2::History)
     return History(times, h1.tmin, h2.tmax, marks, dims, ndims(h1); check_args=false)
 end
 
+"""
+    time_change(h, Λ)
+
+Apply the time rescaling `t -> Λ(t)` to history `h`.
+"""
 function time_change(h::History{T}, Λ) where {T}
     new_times = Λ.(event_times(h))
     new_marks = copy(event_marks(h))

src/mark_distributions.jl

@@ -5,15 +5,8 @@ abstract type AbstractMarkDistribution end
 const PointProcessMarkDistribution = Union{Distribution,AbstractMarkDistribution}
 
 ## Standard implementations. Override as needed
-"""
-    mark_distribution(md, t, h)
-
-Compute the distribution of marks at time `t` after history `h`.
-
-Remark: This method must work for empty histories.
-"""
-function mark_distribution end
-
+# The docstring for `mark_distribution` lives in `abstract_point_process.jl`, which
+# documents both the process and the mark distribution methods of this function.
 function mark_distribution(md::AbstractMarkDistribution, t, h::History)
     return error(
         "Type $(typeof(md)) subtypes `AbstractMarkDistribution` but has " *
@@ -49,6 +42,13 @@ mark_distribution(d::Distribution, t, h::History) = d
 
 StatsAPI.fit(D::Type{<:Distribution}, h::History) = fit(D, h.marks)
 
+# `Type{<:Distribution}` and `Type{<:AbstractPointProcess}` intersect at `Type{Union{}}`,
+# because `Union{}` is a subtype of every type. Without this method, inference of
+# `fit(PP, h)` for an abstract `PP::Type{<:AbstractPointProcess}` therefore includes the
+# return type of the method above (some `Distribution`) in its union, which makes JET
+# report a missing method for every call that passes the result on to a process method.
+StatsAPI.fit(::Type{Union{}}, ::History) = error("unreachable")
+
 # Struct for non-marked processes
 "Mark distribution for non-marked processes. Always return the mark `nothing`."
 struct NoMarks <: AbstractMarkDistribution end

src/multivariate_process.jl

@@ -4,16 +4,19 @@
 Abstract type for multivariate temporal point processes.
 
 To implement a multivariate process, one must subtype `AbstractMultivariateProcess` and provide
-implementations for the methods below. 
+implementations for the methods below.
 - `mark_distribution(pp, t, h, d)`
 - `ground_intensity(pp, t, h, d)`
-- `integrated_ground_intensity(pp, t, h, d)`
+- `integrated_ground_intensity(pp, h, a, b, d)`
 - `ground_intensity_bound(pp, t, h, d)`
-- `fit(Type{pp}, h)`
-- `simulate(pp, h)`
 In all the cases, `pp` is the point process being implemented, `t` is the instant in
-which the function will be evaluated, `h` is the history up to `t` and `d` is the dimension.
-Other methods should be implemented if permance is an issue.
+which the function will be evaluated, `h` is the history up to `t`, `[a, b)` is the interval
+over which the intensity is integrated and `d` is the dimension.
+Given those, this file provides the vector-valued versions of the methods above, plus
+`intensity`, `log_intensity`, `sample_mark` and `time_change`.
+There are no generic fallbacks for `fit(Type{pp}, h)`, `simulate(rng, pp, tmin, tmax)` and
+`logdensityof(pp, h)`, so those must be implemented for each process.
+Other methods should be implemented if performance is an issue.
 The process is expected to have a field `mark_dist::Vector{<:AbstractMarkDistribution}`. If this
 field is not present, `Base.ndims(pp)` must also be implemented.
 """
@@ -64,9 +67,16 @@ function log_intensity(pp::AbstractMultivariateProcess, m, t, h::History)
 end
 
 function time_change(h::History{T}, pp::AbstractMultivariateProcess) where {T}
+    if ndims(h) != ndims(pp)
+        throw(
+            DimensionMismatch(
+                "The history has $(ndims(h)) dimensions but the process has $(ndims(pp))."
+            ),
+        )
+    end
     tmin = typemax(T)
     tmax = typemin(T)
-    transformed_times = [zeros(T, nb_events(h, d)) for d in 1:ndims(h)]
+    transformed_times = [zeros(T, nb_events(h, d)) for d in 1:ndims(pp)]
     for d in 1:ndims(pp)
         Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t, d)
         transformed_times[d] .= Λ.(event_times(h, d))
@@ -75,6 +85,6 @@ function time_change(h::History{T}, pp::AbstractMultivariateProcess) where {T}
         new_tmax = Λ(max_time(h))
         tmax = new_tmax > tmax ? new_tmax : tmax
     end
-    transformed_marks = [collect(event_marks(h, d)) for d in 1:ndims(h)]
+    transformed_marks = [collect(event_marks(h, d)) for d in 1:ndims(pp)]
     return History(transformed_times, tmin, tmax, transformed_marks)
 end

src/simulation.jl

@@ -32,7 +32,7 @@ function simulate_ogata(
             U = rand(rng, typeof(U_max))
             if U < U_max
                 if t + τ < tmax
-                    m = sample_mark(pp.mark_dist, t + τ, h)
+                    m = sample_mark(rng, pp, t + τ, h)
                     push!(h, t + τ, m; check_args=false)
                 end
             end

src/univariate/hawkes/fit.jl

@@ -1,5 +1,5 @@
 """
-    StatsAPI.fit(::Type{HawkesProcess{T}}, h::History; step_tol::Float64 = 1e-6, max_iter::Int = 1000, rng::AbstractRNG=default_rng()) where {T<:Real}
+    StatsAPI.fit(::Type{HawkesProcess{T,MD}}, h::History; step_tol::Float64 = 1e-6, max_iter::Int = 1000, rng::AbstractRNG=default_rng()) where {T<:Real,MD<:PointProcessMarkDistribution}
 
 Expectation-Maximization algorithm from [Lewis2011](@cite).
 The relevant calculations are in page 4, equations 6-13.

src/univariate/poisson/simulation.jl

@@ -1,9 +1,9 @@
 function simulate(rng::AbstractRNG, pp::PoissonProcess, tmin::T, tmax::T) where {T<:Real}
     h = History(tmin, tmax, eltype(pp.mark_dist))
     inter_dist = Exponential(inv(pp.λ))
-    t = T(rand(rng, inter_dist))
+    t = tmin + T(rand(rng, inter_dist))
     while t < tmax
-        m = sample_mark(pp, t, h)
+        m = sample_mark(rng, pp, t, h)
         push!(h, t, m; check_args=false)
         t += rand(rng, inter_dist)
     end

test/custom_processes.jl

@@ -69,4 +69,10 @@ end
     h_transf = time_change(h, fmp)
     @test event_times(h_transf, 1) == event_times(h, 1)
     @test event_times(h_transf, 2) == event_times(h, 2) .* 2
+
+    # A history and a process with different numbers of dimensions must be rejected
+    # explicitly, instead of failing with a `BoundsError` deeper down
+    fmp3 = FakeMultivariatePoisson([1.0, 2.0, 3.0], fill(NoMarks(), 3))
+    @test_throws DimensionMismatch time_change(h, fmp3)
+    @test_throws DimensionMismatch time_change(History([sort(rand(10))], 0.0, 1.0), fmp)
 end

test/history.jl

@@ -29,6 +29,10 @@
     @test event_times(h, 0.2, 0.8) == [0.2]
     @test event_times(h, 0.8, 0.2) == []
     @test event_times(h, nothing) == h.times
+    # `event_times` and `event_marks` must agree on which dimensions exist
+    @test event_times(h, 1) == h.times
+    @test event_marks(h, 1) == h.marks
+    @test length(event_times(h, 2)) == length(event_marks(h, 2)) == 0
     @test event_marks(h) == ["a", "b", "c"]
     @test event_marks(h, 0.2, 0.8) == ["a"]
     @test event_marks(h, 0.8, 0.2) == []
@@ -98,6 +102,15 @@ end
     @test_throws DomainError History(sort(rand(3)), 0, 1, rand(3), [1, 2, 3], 2)
     @test event_dims(History([[0.5]], 0, 1)) == [nothing]
 
+    # A univariate history must be built with `dims` filled with `nothing`
+    @test_throws DomainError History([0.1, 0.5], 0.0, 1.0, marks1, [1, 1], 1)
+
+    # `History(h, d)` must reject dimensions the history does not have, instead of
+    # returning times and marks of different lengths
+    @test_throws DomainError History(h_multi, 3)
+    @test_throws DomainError History(h_multi, 0)
+    @test_throws DomainError History(History([0.2, 0.8], 0.0, 1.0, marks1), 2)
+
     @test_throws DomainError History(
         [1.0, 1.0, 2.0, 3.0, 4.0], 0.0, 5.0, fill(nothing, 5), [1, 1, 1, 2, 1], 2
     )

test/poisson_process.jl

@@ -54,6 +54,22 @@ ppvec = PoissonProcess(1.0, MvNormal(Matrix(I, 2, 2)))
     @test nb_events(h2) > 0
     @test !has_events(h3)
     @test !has_events(h4)
+
+    # Events must lie inside [tmin, tmax), also for tmin != 0
+    h5 = simulate(rng, pp, 5.0, 10.0)
+    @test all(5.0 .<= event_times(h5) .< 10.0)
+    @test logdensityof(pp, h5) < 0
+
+    # Both the times and the marks must be reproducible from a seeded rng
+    h6 = simulate(Xoshiro(1), pp, 0.0, 100.0)
+    h7 = simulate(Xoshiro(1), pp, 0.0, 100.0)
+    @test event_times(h6) == event_times(h7)
+    @test event_marks(h6) == event_marks(h7)
+
+    h8 = simulate_ogata(Xoshiro(1), pp, 0.0, 100.0)
+    h9 = simulate_ogata(Xoshiro(1), pp, 0.0, 100.0)
+    @test event_times(h8) == event_times(h9)
+    @test event_marks(h8) == event_marks(h9)
 end
 
 @testset "Fitting" begin