From b73914becfee96f367ae0949d5286e1a3aac124f Mon Sep 17 00:00:00 2001 From: Ryan Senne <50930199+rsenne@users.noreply.github.com> Date: Fri, 7 Aug 2026 21:05:51 -0400 Subject: [PATCH] Review Changes --- .gitignore | 3 +- .../PPTests/monte_carlo_test.jl | 2 +- src/abstract_point_process.jl | 24 +++++++---- src/history.jl | 40 +++++++++++++++++-- src/mark_distributions.jl | 18 ++++----- src/multivariate_process.jl | 26 ++++++++---- src/simulation.jl | 2 +- src/univariate/hawkes/fit.jl | 2 +- src/univariate/poisson/simulation.jl | 4 +- test/custom_processes.jl | 6 +++ test/history.jl | 13 ++++++ test/poisson_process.jl | 16 ++++++++ 12 files changed, 123 insertions(+), 33 deletions(-) diff --git a/.gitignore b/.gitignore index 81d6668..5d32027 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ *.jl.mem /docs/build/ **/Manifest.toml -.vscode \ No newline at end of file +.vscode +dev/ \ No newline at end of file diff --git a/src/HypothesisTests/PPTests/monte_carlo_test.jl b/src/HypothesisTests/PPTests/monte_carlo_test.jl index 26e27fb..d1fc8d7 100644 --- a/src/HypothesisTests/PPTests/monte_carlo_test.jl +++ b/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 diff --git a/src/abstract_point_process.jl b/src/abstract_point_process.jl index 8003676..e02a260 100644 --- a/src/abstract_point_process.jl +++ b/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,21 +33,28 @@ 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`. """ @@ -54,11 +62,12 @@ 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 @@ -66,13 +75,14 @@ function log_intensity end """ 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 diff --git a/src/history.jl b/src/history.jl index 297ab9f..b55bd73 100644 --- a/src/history.jl +++ b/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)) diff --git a/src/mark_distributions.jl b/src/mark_distributions.jl index 74aa06d..814d54c 100644 --- a/src/mark_distributions.jl +++ b/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 diff --git a/src/multivariate_process.jl b/src/multivariate_process.jl index 9d9b6c6..80ede22 100644 --- a/src/multivariate_process.jl +++ b/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 diff --git a/src/simulation.jl b/src/simulation.jl index 144b308..7cbd9e7 100644 --- a/src/simulation.jl +++ b/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 diff --git a/src/univariate/hawkes/fit.jl b/src/univariate/hawkes/fit.jl index bd67280..955d0ff 100644 --- a/src/univariate/hawkes/fit.jl +++ b/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. diff --git a/src/univariate/poisson/simulation.jl b/src/univariate/poisson/simulation.jl index a7a33ce..cb2448e 100644 --- a/src/univariate/poisson/simulation.jl +++ b/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 diff --git a/test/custom_processes.jl b/test/custom_processes.jl index d6813d8..97ab3ba 100644 --- a/test/custom_processes.jl +++ b/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 diff --git a/test/history.jl b/test/history.jl index b4e49fb..5e1c706 100644 --- a/test/history.jl +++ b/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 ) diff --git a/test/poisson_process.jl b/test/poisson_process.jl index 4b994b6..4ffdb1c 100644 --- a/test/poisson_process.jl +++ b/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