Fixes for PR

Author: JoseKling

docs/examples/CustomProcesses.jl

@@ -35,6 +35,8 @@ using Plots
 using StatsAPI
 using Optim
 
+import PointProcesses: NoMarks, AbstractMarkDistribution, AbstractUnivariateProcess
+
 struct TwoStateModel{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess
     λ::R
     Δ::R
@@ -47,7 +49,7 @@ end
 # `PointProcessMarkDistribution` encompasses `NoMarks` and any distribution in package
 # `Distributions.jl`. In the next chapter we will discuss custom distributions as well.
 
-# > Even if you plan on only using a non-marked process, you should keep the `mark_dist` field. In you really want you can define an inner constructor like `TwoStateModel(λ, Δ, τ) → new(λ, Δ, τ, Nomarks())`.
+# > Even if you plan on only using a non-marked process, you should keep the `mark_dist` field. If you really want you can define an inner constructor like `TwoStateModel(λ, Δ, τ) → new(λ, Δ, τ, Nomarks())`.
 
 # Right now you cannot do much more than accessing the fields. The only methods you have
 # defined are `ndims`, `DensityKind` and `mark_distribution`. Not too interesting.  
@@ -243,10 +245,10 @@ imp = IndependentMultivariateProcess([
 
 sim = simulate(imp, 0.0, 10.0)
 
-print("ground_intensity(imp, 0.0, sim) → $(ground_intensity(imp, 0.0, sim))")
-print("ground_intensity(imp, 0.0, sim, 1) → $(ground_intensity(imp, 0.0, sim, 1))")
-print("intensity(imp, 0.0, 0.0, sim) → $(intensity(imp, 0.0, 0.0, sim))")
-print("intensity(imp, 0.0, nothing, sim, 3) → $(intensity(imp, 0.0, nothing, sim, 3))")
+println("ground_intensity(imp, 0.0, sim) → $(ground_intensity(imp, 0.0, sim))")
+println("ground_intensity(imp, 0.0, sim, 1) → $(ground_intensity(imp, 0.0, sim, 1))")
+println("intensity(imp, 0.0, 0.0, sim) → $(intensity(imp, 0.0, 0.0, sim))")
+println("intensity(imp, 0.0, nothing, sim, 3) → $(intensity(imp, nothing, 0.0, sim, 3))")
 
 # Lets end this tutorial with a nice plot of this process. 
 

src/PointProcesses.jl

@@ -45,7 +45,6 @@ export sample_mark
 
 ## Point processes
 
-export PointProcessMarkDistribution, AbstractMarkDistribution, NoMarks
 export AbstractPointProcess, AbstractUnivariateProcess, AbstractMultivariateProcess
 export BoundedPointProcess
 export ground_intensity, mark_distribution

src/history.jl

@@ -169,7 +169,7 @@ end
 
 function History(h::History, d::Int)
     times = event_times(h, d)
-    marks = event_times(h, d)
+    marks = event_marks(h, d)
     return History(times, min_time(h), max_time(h), marks)
 end
 
@@ -191,8 +191,8 @@ event_times(h::History) = h.times
 
 Return the sorted vector of event times for `h` in dimension `d`.
 """
-function event_times(h::History, d::Int)
-    h.N == 1 && d == 1 ? h.times : (@view h.times[h.dims .== d])
+function event_times(h::History, d::Union{Int, Nothing})
+    h.N == 1 ? h.times : (@view h.times[h.dims .== d])
 end
 
 """
@@ -230,7 +230,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::Int)
+function event_marks(h::History, d::Union{Int, Nothing})
     h.N == 1 && d == 1 ? h.marks : (@view h.marks[h.dims .== d])
 end
 
@@ -369,17 +369,37 @@ duration(h::History) = max_time(h) - min_time(h)
 
 Add event `(t, m)` inside the interval `[h.tmin, h.tmax)` at the end of history `h`.
 
-With `check_args=true` (the default), the event must satisfy
-`h.tmin <= t < h.tmax`, occur at or after the last existing event time, and
-(if `h` is multivariate) lie in dimension `d ∈ 1:h.N`. Violations throw an
-`AssertionError`. Pass `check_args=false` to skip these checks in trusted
-inner loops.
+With `check_args=true` (the default), the method checks if the event `t` satisfies
+`h.tmin <= t < h.tmax`, if it occurs after the last existing event time for that
+dimention, and (if `h` is multivariate) if the dimension lie in `d ∈ 1:h.N`.
+Pass `check_args=false` to skip these checks in trusted inner loops.
 """
 function Base.push!(h::History, t::Real, m=nothing, d=nothing; check_args=true)
     if check_args
-        @assert h.tmin <= t < h.tmax
-        @assert (length(h) == 0) || (h.times[end] < t)
-        @assert (d === nothing && h.N == 1) || 1 <= d <= h.N
+        if !(h.tmin <= t < h.tmax)
+            throw(
+                DomainError(
+                    (t, h.tmin, h.tmax),
+                    "Event time must lie in the half-open interval [tmin, tmax).",
+                ),
+            )
+        end
+        if !((d === nothing && h.N == 1) || (d isa Integer && 1 <= d <= h.N))
+            throw(
+                DomainError(
+                    d,
+                    "Event dimension must be between 1 and h.N for a multivariate history, or omitted for a univariate history.",
+                ),
+            )
+        end
+        if !(isempty(h) || event_times(h, d)[end] < t)
+            throw(
+                DomainError(
+                    (t, event_times(h, d)[end]),
+                    "Event time must be strictly greater than the previous event time in the same dimension.",
+                ),
+            )
+        end
     end
     push!(h.times, t)
     push!(h.marks, m)

src/mark_distributions.jl

@@ -1,4 +1,3 @@
-
 # Type definition
 "Abstract type for defining mark distributions not in `Distributions.jl`"
 abstract type AbstractMarkDistribution end
@@ -43,8 +42,6 @@ function DensityInterface.densityof(md::PointProcessMarkDistribution, t, h::Hist
     densityof(mark_distribution(md, t, h), m)
 end
 
-StatsAPI.fit
-
 # Support for `Distributions.jl`
 mark_distribution(d::Distribution, t, h::History) = d
 

src/simulation.jl

@@ -20,7 +20,7 @@ Simulate a temporal point process `pp` on interval `[tmin, tmax)` using Ogata's
 To infer the type of the marks, the implementation assumes that there is method of `mark_distribution` without the argument `h` such that it corresponds to the distribution of marks in case the history is empty.
 """
 function simulate_ogata(
-    rng::AbstractRNG, pp::AbstractPointProcess, tmin::T, tmax::T
+    rng::AbstractRNG, pp::AbstractUnivariateProcess, tmin::T, tmax::T
 ) where {T<:Real}
     M = eltype(pp.mark_dist)
     h = History(; times=T[], marks=M[], tmin=tmin, tmax=tmax)

src/univariate/poisson/simulation.jl

@@ -1,18 +1,11 @@
 function simulate(rng::AbstractRNG, pp::PoissonProcess, tmin::T, tmax::T) where {T<:Real}
-    times = simulate_poisson_times(rng, pp.λ, tmin, tmax)
-    h_temp = History(times, tmin, tmax)
-    marks = [sample_mark(pp.mark_dist, t, h_temp) for t in times]
-    return History(; times=times, marks=marks, tmin=tmin, tmax=tmax)
+    h = History(T[], tmin, tmax, eltype(pp.mark_dist)[])
+    inter_dist = Exponential(inv(pp.λ))
+    t = rand(rng, inter_dist)
+    while t < tmax
+        m = sample_mark(pp.mark_dist, t, h)
+        push!(h, t, m; check_args=false)
+        t += rand(rng, inter_dist)
+    end
+    return h
 end
-
-# function simulate(rng::AbstractRNG, pp::PoissonProcess, tmin::T, tmax::T) where {T<:Real}
-#     h = History(T[], tmin, tmax, eltype(pp.mark_dist)[])
-#     inter_dist = Exponential(inv(pp.λ))
-#     t = rand(rng, inter_dist)
-#     while t < tmax
-#         m = sample_mark(pp.mark_dist, t, h)
-#         push!(h, t, m; check_args=False)
-#         t = rand(rng, inter_dist)
-#     end
-#     return h
-# end

test/history.jl

@@ -1,16 +1,20 @@
 @testset "Univariate History" begin
+    # Constructors
     h_empty1 = History(0.0, 1.0)
     h_empty2 = History(0.0, 1.0, 2)
 
-    @test h_empty1 isa History{Float64}
-    @test h_empty2 isa History{Float64}
+    @test h_empty1 isa History{Float64,Any}
+    @test h_empty2 isa History{Float64,Any}
     @test ndims(h_empty1) == 1
     @test ndims(h_empty2) == 2
     @test isempty(h_empty1)
     @test isempty(h_empty2)
 
     h = History([0.2, 0.8, 1.1], 0.0, 2.0, ["a", "b", "c"]);
 
+    @test h isa History{Float64,String}
+
+    # Access
     @test duration(h) == 2.0
     @test nb_events(h) == 3
     @test nb_events(h, 1.0, 2.0) == 1
@@ -21,20 +25,20 @@
     @test event_times(h) == [0.2, 0.8, 1.1]
     @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
     @test event_marks(h) == ["a", "b", "c"]
     @test event_marks(h, 0.2, 0.8) == ["a"]
     @test event_marks(h, 0.8, 0.2) == []
     @test ndims(h) == 1
     @test event_dims(h) == fill(nothing, 3)
 
+    # Interface
     push!(h, 1.7, "d")
 
     @test has_events(h, 1.5, 2.0)
 
     h2 = History(; times=[2.3], marks=["e"], tmin=2.0, tmax=2.5)
-    h3 = History(; times=[[1], [2, 2.5]], tmin=0, tmax=3)
 
-    @test (event_marks(h3) == fill(nothing, 3)) && (event_dims(h3) == [1, 2, 2])
     @test string(h2) == "History{Float64,String} with 1 events on interval [2.0, 2.5)"
 
     h_cat = cat(h, h2)
@@ -75,9 +79,20 @@ end
     @test event_marks(h_multi) == ["a", "c", "b", "d"]
     @test event_dims(h_multi) == [1, 2, 1, 2]
 
+    h_multi1 = History(h_multi, 1)
+    h_multi2 = History(h_multi, 2)
+
+    @test event_times(h_multi1) == times1
+    @test event_times(h_multi2) == times2
+    @test event_marks(h_multi1) == marks1
+    @test event_marks(h_multi2) == marks2
+
     @test_throws DomainError History(rand(3), 0, 1, rand(3), [1, 2, 3], 2)
     @test event_dims(History([[0.5]], 0, 1)) == [nothing]
 
+    @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_throws DomainError History([1.0, 1.0, 1.0, 2.0, 3.0], 0.0, 5.0, fill(nothing, 5), [1, 2, 1, 2, 1], 2)
+
     # Test dimension-specific methods
     @test event_times(h_multi, 1) == [0.1, 0.5]
     @test event_times(h_multi, 2) == [0.2, 0.8]

test/hypothesis_tests.jl

@@ -1,6 +1,6 @@
-h1 = History([1, 2, 3, 4], 0, 5)
-h_empty = History(Float64[], 0, 2)
-PP = PoissonProcess{Float32,NoMarks}
+h1 = History([1.0, 2.0, 3.0, 4.0], 0.0, 5.0)
+h_empty = History(Float64[], 0.0, 2.0)
+PP = PoissonProcess{Float64,NoMarks}
 pp = PoissonProcess()
 
 @testset "Statistics" begin
@@ -10,7 +10,7 @@ pp = PoissonProcess()
     @test statistic(KSDistance{Exponential}, pp, h_empty) ≈ 1
 end
 
-h2 = History(collect(0:999), 0, 1000)
+h2 = History(collect(0.0:999.0), 0.0, 1000.0)
 
 @testset "BootstrapTest" begin
     @test_throws ArgumentError BootstrapTest(KSDistance{Uniform}, PP, h_empty)

test/mark_distributions.jl

@@ -42,4 +42,4 @@ end
     @test mark_distribution(md, 0.0, h) == Normal(1.0)
     @test eltype(md) == Float64
     @test densityof(md, 0.0, h, 0.0) == densityof(Normal(1.0), 0.0)
-end
+end

test/multivariate_poisson_process.jl

@@ -47,7 +47,7 @@ end
 
 @testset "Simulation" begin
     pp0 = PoissonProcess([0.0, 1.0, 0.0])
-    bpp = BoundedPointProcess(pp1, 0, 1000)
+    bpp = BoundedPointProcess(pp1, 0.0, 1000.0)
     h1 = simulate(rng, pp1, 0.0, 1000.0)
     h2 = simulate(rng, pp0, 0.0, 1000.0)
     h3 = simulate(rng, bpp)