From 821dccf9eedd71c5ebf1aabfb22cd9f32407877d Mon Sep 17 00:00:00 2001 From: jkling Date: Tue, 19 May 2026 11:26:43 +0200 Subject: [PATCH 1/6] Implemented interface for mark mark distribution --- src/PointProcesses.jl | 4 +- src/abstract_point_process.jl | 2 +- src/history.jl | 7 +- src/mark_distributions.jl | 54 ++++++++++++ .../multivariate_poisson_process.jl | 2 +- src/simulation.jl | 4 +- src/univariate/poisson/fit.jl | 11 --- src/univariate/poisson/inhomogeneous/fit.jl | 8 +- .../inhomogeneous_poisson_process.jl | 15 +--- src/univariate/poisson/poisson_process.jl | 31 ++----- src/univariate/poisson/simulation.jl | 17 +++- test/bounded_point_process.jl | 2 +- test/hypothesis_tests.jl | 3 +- test/inhomogeneous_poisson_process.jl | 86 +++++++++---------- test/poisson_process.jl | 6 +- 15 files changed, 142 insertions(+), 110 deletions(-) create mode 100644 src/mark_distributions.jl diff --git a/src/PointProcesses.jl b/src/PointProcesses.jl index e37f5da..bbfeb02 100644 --- a/src/PointProcesses.jl +++ b/src/PointProcesses.jl @@ -8,7 +8,7 @@ module PointProcesses # Imports using DensityInterface: DensityInterface, HasDensity, densityof, logdensityof -using Distributions: Distributions, UnivariateDistribution, MultivariateDistribution +using Distributions: Distributions, Distribution, UnivariateDistribution, MultivariateDistribution using Distributions: Categorical, Exponential, Poisson, Uniform, Dirac, Gamma using Distributions: fit, suffstats, probs using Integrals: Integrals, IntegralProblem, solve, QuadGKJL @@ -40,6 +40,7 @@ export time_change, split_into_chunks ## Point processes +export PointProcessMarkDistribution, AbstractMarkDistribution, NoMarks export AbstractPointProcess, AbstractUnivariateProcess, AbstractMultivariateProcess export BoundedPointProcess export ground_intensity, mark_distribution @@ -79,6 +80,7 @@ export PointProcessTest, BootstrapTest, MonteCarloTest ## General include("history.jl") +include("mark_distributions.jl") include("abstract_point_process.jl") include("simulation.jl") include("bounded_point_process.jl") diff --git a/src/abstract_point_process.jl b/src/abstract_point_process.jl index 787b7b5..9209ff7 100644 --- a/src/abstract_point_process.jl +++ b/src/abstract_point_process.jl @@ -55,7 +55,7 @@ For multivariate processes, it returns a vector of mark distributions for each d mark_distribution(pp, h, t, d) computes the mark distribution for a multivariate process `pp` at dimension `d`. """ -function mark_distribution end +mark_distribution(pp::AbstractPointProcess, t, h) = mark_distribution(pp.mark_dist, t, h) """ intensity(pp, m, t, h) diff --git a/src/history.jl b/src/history.jl index 3502d96..8816694 100644 --- a/src/history.jl +++ b/src/history.jl @@ -125,6 +125,11 @@ function Base.show(io::IO, h::History{T,M}) where {T,M} ) end +function History(tmin::R1, tmax::R2, N::Int=1) where {R1<:Real, R2<:Real} + R = promote_type(R1, R2) + History(R[], tmin, tmax, [], [], N) +end + """ event_times(h) @@ -318,7 +323,7 @@ Add event `(t, m)` inside the interval `[h.tmin, h.tmax)` at the end of history 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 (length(h) == 0) || (h.times[end] < t) @assert (d === nothing && h.N == 1) || 1 <= d <= h.N end push!(h.times, t) diff --git a/src/mark_distributions.jl b/src/mark_distributions.jl new file mode 100644 index 0000000..392902f --- /dev/null +++ b/src/mark_distributions.jl @@ -0,0 +1,54 @@ + +# Type definition +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`. +""" +function mark_distribution end + +function mark_distribution(md::AbstractMarkDistribution, t, h::History) + error("Type $(typeof(md)) subtypes `AbstractMarkDistribution` but has " * + "not implemented the required `mark_distribution(md, t, h)` method.") +end + +""" + sample_mark(rng, md, t, h) + +Return one sample from the distribution of marks at time `t` after history `h`, using the random number generator `rng`. +""" +sample_mark(rng::AbstractRNG, md::PointProcessMarkDistribution, t, h::History) = rand(rng, mark_distribution(md, t, h)) + +sample_mark(md::PointProcessMarkDistribution, t, h::History) = sample_mark(default_rng(), md, t, h) + +Base.eltype(md::AbstractMarkDistribution) = eltype(mark_distribution(md, 0.0, History(0.0, 1.0))) + +DensityInterface.densityof(md::PointProcessMarkDistribution, t, h::History, m) = + densityof(mark_distribution(md, t, h), m) + +# Support for `Distributions.jl` +mark_distribution(d::Distribution, t, h::History) = d + +StatsAPI.fit(D::Type{<:Distribution}, h::History) = fit(D, h.marks) + +# Struct for non-marked processes +struct NoMarks <: AbstractMarkDistribution end + +mark_distribution(::NoMarks, t, h::History) = Dirac(nothing) + +sample_mark(::AbstractRNG, ::NoMarks, t, ::History) = nothing + +Base.eltype(::NoMarks) = Nothing + +StatsAPI.fit(::Type{NoMarks}, h) = NoMarks() + +StatsAPI.fit(::Type{NoMarks}, marks, weights) = NoMarks() + +DensityInterface.densityof(::NoMarks, t, h, ::Nothing) = 1.0 + +DensityInterface.densityof(::NoMarks, t, h, m) = 0.0 diff --git a/src/multivariate/multivariate_poisson_process.jl b/src/multivariate/multivariate_poisson_process.jl index 43a2190..052d3db 100644 --- a/src/multivariate/multivariate_poisson_process.jl +++ b/src/multivariate/multivariate_poisson_process.jl @@ -16,7 +16,7 @@ function PoissonProcess( end function PoissonProcess(λ::AbstractVector{R}; check_args::Bool=true) where {R<:Real} - return PoissonProcess(λ, [Dirac(nothing) for _ in eachindex(λ)]; check_args=check_args) + return PoissonProcess(λ, [NoMarks() for _ in eachindex(λ)]; check_args=check_args) end function PoissonProcess( diff --git a/src/simulation.jl b/src/simulation.jl index 08d4fd5..928c265 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -22,7 +22,7 @@ To infer the type of the marks, the implementation assumes that there is method function simulate_ogata( rng::AbstractRNG, pp::AbstractPointProcess, tmin::T, tmax::T ) where {T<:Real} - M = typeof(rand(mark_distribution(pp, tmin))) + M = eltype(pp.mark_dist) h = History(; times=T[], marks=M[], tmin=tmin, tmax=tmax) t = tmin while t < tmax @@ -34,8 +34,8 @@ function simulate_ogata( U_max = ground_intensity(pp, t + τ, h) / B U = rand(rng, typeof(U_max)) if U < U_max - m = rand(rng, mark_distribution(pp, t + τ, h)) if t + τ < tmax + m = sample_mark(pp.mark_dist, t + τ, h) push!(h, t + τ, m; check_args=false) end end diff --git a/src/univariate/poisson/fit.jl b/src/univariate/poisson/fit.jl index 0c34b0e..52ade53 100644 --- a/src/univariate/poisson/fit.jl +++ b/src/univariate/poisson/fit.jl @@ -1,16 +1,5 @@ ## Fit MLE -#= -A separate `fit` method for unmarked Poisson processes is needed, because -`Distributions.jl` does not provide a `fit` method for the `Dirac` distribution -=# -function StatsAPI.fit( - ::Type{PoissonProcess{R,Dirac{Nothing}}}, ss::PoissonProcessStats{R1,R2}; kwargs... -) where {R<:Real,R1<:Real,R2<:Real} - λ = convert(R, ss.nb_events / ss.duration) - return PoissonProcess(λ, Dirac(nothing)) -end - function StatsAPI.fit( ::Type{PoissonProcess{R,D}}, ss::PoissonProcessStats{R1,R2}; kwargs... ) where {R<:Real,D,R1<:Real,R2<:Real} diff --git a/src/univariate/poisson/inhomogeneous/fit.jl b/src/univariate/poisson/inhomogeneous/fit.jl index abba084..a38c1a2 100644 --- a/src/univariate/poisson/inhomogeneous/fit.jl +++ b/src/univariate/poisson/inhomogeneous/fit.jl @@ -81,7 +81,7 @@ function StatsAPI.fit( end function StatsAPI.fit( - ::Type{<:InhomogeneousPoissonProcess{F,Dirac{Nothing}}}, + ::Type{<:InhomogeneousPoissonProcess{F,NoMarks}}, h::History, init_params; integration_config=IntegrationConfig(), @@ -89,7 +89,7 @@ function StatsAPI.fit( ) where {F<:ParametricIntensity} # For unmarked processes, skip fitting the mark distribution intensity = fit(F, h, init_params; integration_config=integration_config, kwargs...) - return InhomogeneousPoissonProcess(intensity, Dirac(nothing)) + return InhomogeneousPoissonProcess(intensity, NoMarks()) end function StatsAPI.fit( @@ -127,7 +127,7 @@ function StatsAPI.fit( end function StatsAPI.fit( - ::Type{<:InhomogeneousPoissonProcess{PiecewiseConstantIntensity{R},Dirac{Nothing}}}, + ::Type{<:InhomogeneousPoissonProcess{PiecewiseConstantIntensity{R},NoMarks}}, h::History, breakpoints::Vector; kwargs..., @@ -149,7 +149,7 @@ function StatsAPI.fit( end intensity = PiecewiseConstantIntensity(breakpoints, rates) - return InhomogeneousPoissonProcess(intensity, Dirac(nothing)) + return InhomogeneousPoissonProcess(intensity, NoMarks()) end function StatsAPI.fit( diff --git a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl index 0ce8e44..29e027a 100644 --- a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl +++ b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl @@ -50,8 +50,8 @@ end function InhomogeneousPoissonProcess( f::F; integration_config::C=IntegrationConfig() ) where {F,C} - return InhomogeneousPoissonProcess{F,Dirac{Nothing},C}( - f, Dirac(nothing), integration_config + return InhomogeneousPoissonProcess{F,NoMarks,C}( + f, NoMarks(), integration_config ) end @@ -61,21 +61,12 @@ function Base.show(io::IO, pp::InhomogeneousPoissonProcess) ) end -## Access methods - -""" -Mark distribution is time-independent for this implementation. -""" -mark_distribution(pp::InhomogeneousPoissonProcess) = pp.mark_dist - ## AbstractPointProcess interface ground_intensity(pp::InhomogeneousPoissonProcess, t, h) = pp.intensity_function(t) -mark_distribution(pp::InhomogeneousPoissonProcess, t, h) = pp.mark_dist -mark_distribution(pp::InhomogeneousPoissonProcess, t) = pp.mark_dist function intensity(pp::InhomogeneousPoissonProcess, m, t, h) - return ground_intensity(pp, t, h) * densityof(mark_distribution(pp, t, h), m) + return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m) end """ diff --git a/src/univariate/poisson/poisson_process.jl b/src/univariate/poisson/poisson_process.jl index 4b13ad9..1f372fc 100644 --- a/src/univariate/poisson/poisson_process.jl +++ b/src/univariate/poisson/poisson_process.jl @@ -12,7 +12,7 @@ Homogeneous temporal Poisson process with arbitrary mark distribution. PoissonProcess(λ, mark_dist) """ -struct PoissonProcess{R<:Real,D} <: AbstractUnivariateProcess +struct PoissonProcess{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess λ::R mark_dist::D @@ -29,7 +29,7 @@ struct PoissonProcess{R<:Real,D} <: AbstractUnivariateProcess end function PoissonProcess(λ::R; check_args::Bool=true) where {R<:Real} - return PoissonProcess(λ, Dirac(nothing); check_args=check_args) + return PoissonProcess(λ, NoMarks(); check_args=check_args) end PoissonProcess() = PoissonProcess(1.0) @@ -39,39 +39,24 @@ function Base.show(io::IO, pp::PoissonProcess) end ## Access -ground_intensity(pp::PoissonProcess) = pp.λ -mark_distribution(pp::PoissonProcess) = pp.mark_dist +ground_intensity(pp::PoissonProcess, t, h) = pp.λ -## Intensity functions -function intensity(pp::PoissonProcess, m) - return ground_intensity(pp) * densityof(mark_distribution(pp), m) -end - -function log_intensity(pp::PoissonProcess, m) - return log(ground_intensity(pp)) + logdensityof(mark_distribution(pp), m) -end +intensity(pp::PoissonProcess, m, t, h) = pp.λ * densityof(pp.mark_dist, t, h, m) ## Time change function time_change(h::History, pp::PoissonProcess) - times = (h.times .- h.tmin) .* ground_intensity(pp) - tmax = (h.tmax - h.tmin) * ground_intensity(pp) + times = (h.times .- h.tmin) .* pp.λ + tmax = (h.tmax - h.tmin) * pp.λ return History(; times=times, tmin=0, tmax=tmax, marks=h.marks) end ## Implementing AbstractPointProcess interface - -ground_intensity(pp::PoissonProcess, t, h) = ground_intensity(pp) -mark_distribution(pp::PoissonProcess, t, h) = mark_distribution(pp) -mark_distribution(pp::PoissonProcess, t) = mark_distribution(pp) # For simulate_ogata -intensity(pp::PoissonProcess, m, t, h) = intensity(pp, m) -log_intensity(pp::PoissonProcess, m, t, h) = log_intensity(pp, m) - function ground_intensity_bound(pp::PoissonProcess, t::T, h) where {T<:Real} - B = ground_intensity(pp) + B = pp.λ L = typemax(T) return (B, L) end function integrated_ground_intensity(pp::PoissonProcess, h, a, b) - return ground_intensity(pp) * (b - a) + return pp.λ * (b - a) end diff --git a/src/univariate/poisson/simulation.jl b/src/univariate/poisson/simulation.jl index aa22d8d..63d58c4 100644 --- a/src/univariate/poisson/simulation.jl +++ b/src/univariate/poisson/simulation.jl @@ -1,5 +1,18 @@ function simulate(rng::AbstractRNG, pp::PoissonProcess, tmin::T, tmax::T) where {T<:Real} - times = simulate_poisson_times(rng, ground_intensity(pp), tmin, tmax) - marks = [rand(rng, mark_distribution(pp)) for _ in 1:length(times)] + 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) 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 \ No newline at end of file diff --git a/test/bounded_point_process.jl b/test/bounded_point_process.jl index 494d756..e5ece17 100644 --- a/test/bounded_point_process.jl +++ b/test/bounded_point_process.jl @@ -8,7 +8,7 @@ h = simulate(rng, bpp) @test max_time(bpp) == 1000.0 @test ground_intensity(bpp, 0, h) ≈ intensities @test mark_distribution(bpp, 100.0, h, 1) == Normal() -@test mark_distribution(bpp, 0.0) == fill(Normal(), length(intensities)) +@test mark_distribution(bpp, 0.0, h) == fill(Normal(), length(intensities)) @test intensity(bpp, 0.0, 0.0, h, 1) ≈ pdf(Normal(), 0.0) * intensities[1] @test log_intensity(bpp, 1.0, 1.0, h, 2) ≈ log(pdf(Normal(), 1.0)) + log(intensities[2]) diff --git a/test/hypothesis_tests.jl b/test/hypothesis_tests.jl index 76cfa9e..e75a203 100644 --- a/test/hypothesis_tests.jl +++ b/test/hypothesis_tests.jl @@ -1,12 +1,11 @@ h1 = History([1, 2, 3, 4], 0, 5) h_empty = History(Float64[], 0, 2) -PP = PoissonProcess{Float32,Dirac{Nothing}} +PP = PoissonProcess{Float32,NoMarks} pp = PoissonProcess() @testset "Statistics" begin @test statistic(KSDistance{Uniform}, pp, h1) ≈ 0.2 @test statistic(KSDistance{Exponential}, pp, h1) ≈ 1 - exp(-1) - @test statistic(KSDistance{Uniform}, pp, h_empty) ≈ 1 @test statistic(KSDistance{Exponential}, pp, h_empty) ≈ 1 end diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl index 9d18257..0ce768b 100644 --- a/test/inhomogeneous_poisson_process.jl +++ b/test/inhomogeneous_poisson_process.jl @@ -39,7 +39,6 @@ rng = Random.seed!(12345) @test ground_intensity(pp, 0.0, h_empty) ≈ 1.0 @test ground_intensity(pp, 2.0, h_empty) ≈ 2.0 @test mark_distribution(pp, 0.0, h_empty) isa Normal - @test mark_distribution(pp) isa Normal end @testset "Simulation" begin @@ -791,8 +790,7 @@ end @test pp isa InhomogeneousPoissonProcess @test pp.intensity_function === intensity - @test pp.mark_dist isa Dirac{Nothing} - @test pp.mark_dist.value === nothing + @test pp.mark_dist isa NoMarks @test pp.integration_config isa IntegrationConfig # Test that it simulates correctly @@ -808,7 +806,7 @@ end @test pp isa InhomogeneousPoissonProcess @test pp.intensity_function === intensity - @test pp.mark_dist isa Dirac{Nothing} + @test pp.mark_dist isa NoMarks @test pp.integration_config === custom_config @test pp.integration_config.abstol == 1e-10 @test pp.integration_config.reltol == 1e-10 @@ -837,27 +835,26 @@ end end end -@testset "Fitting - PiecewiseConstant with Dirac (unmarked)" begin +@testset "Fitting - PiecewiseConstant with no marks" begin @testset "Fit with explicit breakpoints vector" begin # Generate data from known piecewise constant process breakpoints_true = [0.0, 3.0, 7.0, 10.0] rates_true = [1.5, 4.0, 2.5] intensity_true = PiecewiseConstantIntensity(breakpoints_true, rates_true) - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 10.0) # Fit using same breakpoints pp_est = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h, breakpoints_true, ) @test pp_est isa InhomogeneousPoissonProcess @test pp_est.intensity_function isa PiecewiseConstantIntensity - @test pp_est.mark_dist isa Dirac{Nothing} - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks @test length(pp_est.intensity_function.rates) == 3 # Check breakpoints are preserved @@ -877,13 +874,13 @@ end # Create uneven breakpoints custom_breakpoints = [0.0, 2.0, 5.0, 12.0] intensity_true = PiecewiseConstantIntensity([0.0, 6.0, 12.0], [3.0, 2.0]) - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 12.0) # Fit with different breakpoints than data generation pp_est = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h, custom_breakpoints, ) @@ -891,7 +888,7 @@ end @test pp_est isa InhomogeneousPoissonProcess @test pp_est.intensity_function.breakpoints == custom_breakpoints @test length(pp_est.intensity_function.rates) == 3 - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks # All rates should be positive @test all(r >= 0 for r in pp_est.intensity_function.rates) @@ -904,7 +901,7 @@ end h = History([4.5, 4.8, 5.2, 5.5, 5.9], 0.0, 10.0) pp_est = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h, breakpoints, ) @@ -924,19 +921,19 @@ end # Edge case: single bin (constant rate) breakpoints = [0.0, 10.0] intensity_true = PolynomialIntensity([3.0]) # constant - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 10.0) pp_est = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h, breakpoints, ) @test pp_est isa InhomogeneousPoissonProcess @test length(pp_est.intensity_function.rates) == 1 - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks # Rate should be approximately count / duration expected_rate = length(h.times) / 10.0 @@ -949,7 +946,7 @@ end h = History([1.0, 2.0, 6.0, 7.0, 8.0], 0.0, 10.0) pp_est = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float32},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float32},NoMarks}, h, breakpoints_f32, ) @@ -959,25 +956,24 @@ end end end -@testset "Fitting - ParametricIntensity with Dirac (unmarked)" begin - @testset "Polynomial intensity with Dirac marks" begin +@testset "Fitting - ParametricIntensity with no marks" begin + @testset "Polynomial intensity with no marks" begin # Generate data from linear process intensity_true = PolynomialIntensity([2.0, 0.3]) - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 20.0) # Fit polynomial intensity pp_est = fit( - InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks}, h, [2.0, 0.3], ) @test pp_est isa InhomogeneousPoissonProcess @test pp_est.intensity_function isa PolynomialIntensity - @test pp_est.mark_dist isa Dirac{Nothing} - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks # Check parameters are reasonable @test abs(pp_est.intensity_function.coefficients[1] - 2.0) < 1.0 @@ -992,13 +988,13 @@ end @testset "Polynomial intensity with log link" begin # Generate data with log link (ensures positivity) intensity_true = PolynomialIntensity([1.5, 0.05]; link=:log) - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 30.0) # Fit with log link pp_est = fit( - InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks}, h, [1.5, 0.05]; link=:log, @@ -1006,7 +1002,7 @@ end @test pp_est isa InhomogeneousPoissonProcess @test pp_est.intensity_function.link === :log - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks # Verify positivity over range @test all(pp_est.intensity_function(t) > 0 for t in 0.0:0.5:30.0) @@ -1016,24 +1012,23 @@ end @test abs(pp_est.intensity_function.coefficients[2] - 0.05) < 0.1 end - @testset "Exponential intensity with Dirac marks" begin + @testset "Exponential intensity with no marks" begin # Generate data from exponential process intensity_true = ExponentialIntensity(3.0, 0.08) - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 15.0) # Fit exponential intensity pp_est = fit( - InhomogeneousPoissonProcess{ExponentialIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{ExponentialIntensity{Float64},NoMarks}, h, [log(3.0), 0.08], ) @test pp_est isa InhomogeneousPoissonProcess @test pp_est.intensity_function isa ExponentialIntensity - @test pp_est.mark_dist isa Dirac{Nothing} - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks # Check parameters are in reasonable range @test 0.5 * intensity_true.a <= @@ -1045,16 +1040,16 @@ end @test all(pp_est.intensity_function(t) > 0 for t in 0.0:0.5:15.0) end - @testset "Sinusoidal intensity with Dirac marks" begin + @testset "Sinusoidal intensity with no marks" begin # Generate data from sinusoidal process intensity_true = SinusoidalIntensity(6.0, 2.5, 2π, 0.0) - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 8.0) # Fit sinusoidal intensity pp_est = fit( - InhomogeneousPoissonProcess{SinusoidalIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{SinusoidalIntensity{Float64},NoMarks}, h, [log(6.0), 0.4, 0.0]; ω=2π, @@ -1062,8 +1057,7 @@ end @test pp_est isa InhomogeneousPoissonProcess @test pp_est.intensity_function isa SinusoidalIntensity - @test pp_est.mark_dist isa Dirac{Nothing} - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks @test pp_est.intensity_function.ω ≈ 2π # Check constraint a >= |b| is satisfied @@ -1080,44 +1074,44 @@ end @testset "Custom integration config" begin # Test that integration_config is properly passed through intensity_true = PolynomialIntensity([2.5, 0.2]) - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 20.0) # Fit with custom integration config custom_config = IntegrationConfig(abstol=1e-10, reltol=1e-10, maxiters=5000) pp_est = fit( - InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks}, h, [2.5, 0.2]; integration_config=custom_config, ) @test pp_est isa InhomogeneousPoissonProcess - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks # Verify fitting succeeded with custom config @test pp_est.intensity_function isa PolynomialIntensity @test abs(pp_est.intensity_function.coefficients[1] - 2.5) < 1.5 end - @testset "Quadratic polynomial with Dirac marks" begin + @testset "Quadratic polynomial with no marks" begin # Test higher-degree polynomial intensity_true = PolynomialIntensity([1.0, 0.5, 0.05]) - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 15.0) # Fit quadratic pp_est = fit( - InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks}, h, [1.0, 0.5, 0.05], ) @test pp_est isa InhomogeneousPoissonProcess @test length(pp_est.intensity_function.coefficients) == 3 - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks # Check all coefficients are in reasonable range for i in 1:3 @@ -1130,13 +1124,13 @@ end @testset "Pass through additional kwargs" begin # Test that optimizer kwargs are properly passed through intensity_true = PolynomialIntensity([2.0, 0.1]) - pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing)) + pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks()) h = simulate(rng, pp_true, 0.0, 25.0) # Fit with custom optimizer settings pp_est = fit( - InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks}, h, [2.0, 0.1]; optimizer=LBFGS(), @@ -1145,6 +1139,6 @@ end @test pp_est isa InhomogeneousPoissonProcess @test pp_est.intensity_function isa PolynomialIntensity - @test pp_est.mark_dist.value === nothing + @test pp_est.mark_dist isa NoMarks end end diff --git a/test/poisson_process.jl b/test/poisson_process.jl index 23ea80b..00ea0f0 100644 --- a/test/poisson_process.jl +++ b/test/poisson_process.jl @@ -13,8 +13,8 @@ rng = Random.seed!(63) @test pp2.λ == 1.0 @test pp3.λ == 1.0 @test pp1.mark_dist == Normal() - @test pp2.mark_dist == Dirac(nothing) - @test pp3.mark_dist == Dirac(nothing) + @test pp2.mark_dist == NoMarks() + @test pp3.mark_dist == NoMarks() @test_throws DomainError PoissonProcess(-1.0) end @@ -27,7 +27,7 @@ pp = PoissonProcess(1.0, Normal()) h = History([0.1, 0.5], 0, 1, [0.0, 0.0]) @test ground_intensity(pp, 0, h) == 1.0 - @test mark_distribution(pp) == Normal() + @test mark_distribution(pp, 0.0, h) == Normal() @test intensity(pp, 0.5, 0, h) == 1.0 * pdf(Normal(), 0.5) @test log_intensity(pp, 0.5, 0, h) == log(1.0) + logpdf(Normal(), 0.5) @test integrated_ground_intensity(pp, h, 0.0, 1.0) == 1.0 From 2f317693b9a59d6d41c4973eca1a600bd934b55e Mon Sep 17 00:00:00 2001 From: jkling Date: Tue, 26 May 2026 09:11:05 +0200 Subject: [PATCH 2/6] Fixed last bugs and added `CustomProcesses` to docs --- docs/Project.toml | 1 + docs/examples/Basics.jl | 2 +- docs/examples/CustomProcesses.jl | 275 +++++++++++++++++ docs/examples/Hawkes.jl | 2 +- docs/examples/Inhomogeneous.jl | 2 +- docs/src/api.md | 21 +- docs/src/assets/selva2022.png | Bin 0 -> 70137 bytes docs/src/examples/Basics.md | 2 +- docs/src/examples/CustomProcesses.md | 302 +++++++++++++++++++ docs/src/examples/Hawkes.md | 64 ++-- docs/src/examples/Inhomogeneous.md | 2 +- docs/src/refs.bib | 12 + src/PointProcesses.jl | 4 + src/abstract_point_process.jl | 11 +- src/history.jl | 13 +- src/mark_distributions.jl | 10 +- src/multivariate/independent_multivariate.jl | 18 +- src/univariate/poisson/inhomogeneous/fit.jl | 42 +-- test/mark_distributions.jl | 26 ++ test/multivariate_poisson_process.jl | 14 +- test/runtests.jl | 3 + 21 files changed, 733 insertions(+), 93 deletions(-) create mode 100644 docs/examples/CustomProcesses.jl create mode 100644 docs/src/assets/selva2022.png create mode 100644 docs/src/examples/CustomProcesses.md create mode 100644 test/mark_distributions.jl diff --git a/docs/Project.toml b/docs/Project.toml index 9ca9240..f8799cd 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,6 +8,7 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PointProcesses = "af0b7596-9bb0-472a-a012-63904f2b4c55" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/docs/examples/Basics.jl b/docs/examples/Basics.jl index 6f10d9f..96226c7 100644 --- a/docs/examples/Basics.jl +++ b/docs/examples/Basics.jl @@ -157,7 +157,7 @@ p_waiting # Let's fit the actual model now using PointProcesses.jl and run some formal statistical tests. # ### Fitting a Homogeneous Poisson Process -pp_model = fit(PoissonProcess{Float64,Dirac{Nothing}}, h) +pp_model = fit(PoissonProcess{Float64,NoMarks}, h) println("Estimated rate λ̂ = ", pp_model.λ) # hide diff --git a/docs/examples/CustomProcesses.jl b/docs/examples/CustomProcesses.jl new file mode 100644 index 0000000..92214ac --- /dev/null +++ b/docs/examples/CustomProcesses.jl @@ -0,0 +1,275 @@ +# # Defining Custom Processes + +# Although several processes are already implemented and more are to come, we might need +# to define our own custom processes. +# Here we will explain the interface for defining new processes, starting with a simple case +# and building up to more complex cases (which should not be too complicated). + +# ## Basic interface + +# For this example we will implement the two-state model from [Selva2022](@cite). +# This model takes three parameters, λₗ, λₕ and τ. The process starts +# as a Poisson model with intensity λ = λₗ. Whenever an event occurs, the process +# immediately switches to a Poisson process with intensity λ = λₕ for τ +# units of time. If an event occurs during this high intensity phase, the "counter" τ simply +# resets. You can visualize it below. + +# ![Adapted from Selva (2022)](../assets/selva2022.png) + +# If the history ℋₜ up to time t contains the events t₁ < t₂ < ... < tₙ, then the equation +# for the conditional intensity is + +# λ(t | ℋₜ) = λₗ + 𝟙(tₙ > t - τ) (λₕ - λₗ), + +# where 𝟙 is the indicator function. Notice that tₙ < t by definition of ℋₜ. + +# We will use a different parametrization, namely λ = λₗ and Δ = λₕ - λₗ, so the constraint +# λₕ > λₗ becomes Δ > 0. + +# The first step in our implementation is to define this new struct by subtyping +# the `AbstractUnivariateProcess` type (we will get to multivariate processes). + +using PointProcesses +using Distributions +using Plots +using StatsAPI +using Optim + +struct TwoStateModel{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess + λ::R + Δ::R + τ::R + mark_dist::D +end + +# > Note: Subtype `AbstractUnivariateProcess`, not just `AbstractPointProcess`, this way we can use `TwoStateModel` to build multivariate processes with `IndependentMultivariateProcess`. + +# `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())`. + +# 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. + +# We will begin by implementing the methods `ground_intensity`, `integrated_ground_intensity` and +# `ground_intensity_bound`. They are sufficient for almost all methods to function properly. + +function PointProcesses.ground_intensity(pp::TwoStateModel, t, h::History) + if isempty(h) || t <= h.times[1] + return pp.λ + else + t_n = h.times[searchsortedfirst(h.times, t) - 1] # Last event time before `t` + in_high_phase = t < t_n + pp.τ + end + return pp.λ + (in_high_phase * pp.Δ) +end + +function PointProcesses.ground_intensity_bound(pp::TwoStateModel, t, h::History) + return (pp.λ + pp.Δ, Inf) # Also not the most efficient implementation, but works +end + +function PointProcesses.integrated_ground_intensity(pp::TwoStateModel, h, a, b) + integral = pp.λ * (b - a) + interval_events = event_times(h, a, b) + if length(interval_events) > 0 + for i in 1:length(interval_events) - 1 + integral += pp.Δ * max(interval_events[i + 1] - interval_events[i], pp.τ) + end + integral += pp.Δ * max(b - interval_events[end], pp.τ) + end + return integral +end + +# > Notice that in `ground_intensity` we searched for the last event that occurred before time `t`. We need to assume that `h` might be the full history, and not only the history up to time `t`. + +# With this we can now call many other methods: `intensity`, `log_intensity`, `logdensityof`, `time_change` +# and `simulate`. Here are the dependencies: +# - `intensity` → `ground_intensity` + `densityof` from `mark_dist` +# - `log_intensity` → `intensity` +# - `logdensityof` → `integrated_ground_intensity` + `logdensityof` +# - `time_change` → `integrated_ground_intensity` +# - `simulate` → `ground_intensity` + `ground_intensity_bound` + +pp = TwoStateModel(1.0, 2.0, 0.5, Normal()) +h = simulate(pp, 0.0, 20.0) + +println("intensity(pp, m, t, h) = $(intensity(pp, 0.0, 0.0, h))") +println("log_intensity(pp, m, t, h) = $(log_intensity(pp, 0.0, 0.0, h))") +println("logdensityof(pp, h) = $(logdensityof(pp, h))") +println("time_change(pp, h) = $(time_change(h, pp))") +println("simulate(pp, 0.0, 10.0) = $(simulate(pp, 0.0, 10.0))") + +# Let's plot our simulated process `h` along with its ground intensity. + +sim_plot = plot( + event_times(h), + fill(1.0, nb_events(h)) + ; line=:stem, + grid=false, + label=false, + xlabel="Time", + ylabel="Events", + yaxis=false, + xlim=(min_time(h), max_time(h)) +) + +xs = LinRange(min_time(h), max_time(h), 1000) +intensity_plot = plot( + xs, + [ground_intensity(pp, x, h) for x in xs], + ; xaxis=false, + xgrid=false, + label=false, + ylabel="Ground Intensity", + xlim=(min_time(h), max_time(h)) +) + +plot( + intensity_plot, + sim_plot, + title=["Two State Model Simulation" ""], + layout=grid(2, 1, heights=[0.7, 0.3]) +) + +# Lastly, we can define the `fit` method for estimating the parameters from an observed +# history. `logdensityof` is simply the log-likelihood of our process, so for a quick +# and simple example we keep the mark distribution fixed and only optimize the core +# process parameters. + +function StatsAPI.fit( + ::Type{TwoStateModel{R1,NoMarks}}, + h::History, + init_params::Vector{R2}, +) where {R1<:Real,R2<:Real} + objective(params) = -logdensityof(TwoStateModel(params..., NoMarks()), h) + + lower_bound = [0.0, 0.0, 0.0] + upper_bound = [Inf, Inf, Inf] + result = optimize(objective, lower_bound, upper_bound, init_params) + optimal = Optim.minimizer(result) + + return TwoStateModel(optimal..., NoMarks()) +end + +# > It is tricky to properly estimate the parameters from this process. This implementation is conceptually correct, but it should be improved in real applications. + +# With the `fit`, `simulate` and `time_change` methods, we can perform +# goodness-of-fit tests for our model. + +function random_params() + λ = 10.0 + rand() * 10 + Δ = λ * rand() + τ = inv(λ + Δ) * (0.5 + rand()) + return [λ, Δ, τ] +end + +true_params = random_params() +h_unknown = simulate(TwoStateModel(true_params..., NoMarks()), 0.0, 100.0) + +init_params = random_params() +pp_estimated = fit( + TwoStateModel{Float64,NoMarks}, + h_unknown, + init_params, +) +estimated_params = [pp_estimated.λ, pp_estimated.Δ, pp_estimated.τ] # hide + +println("True and estimated parameters:") # hide +println("λ: $(true_params[1]) - $(estimated_params[1])") #hide +println("Δ: $(true_params[2]) - $(estimated_params[2])") #hide +println("τ : $(true_params[3]) - $(estimated_params[3])") # hide + +pp_to_test = TwoStateModel(20.0, 30.0, 1.0/30.0, NoMarks()) +test_result = MonteCarloTest(KSDistance{Exponential}, pp_to_test, h_unknown) + +println("p-value for hypothesis test: $(pvalue(test_result))") # hide + +# > The statistics for goodness of fit tests currently implemented are suited only for testing the distribution of event times. They will run for marked processes, but the test results are only related to the event times. + +# ## Custom Mark Distributions + +# Mark distributions can also be customized. Analogously to processes, we need +# to subtype the `AbstractMarkDistribution` type. Lets make a mark distribution +# that changes with times. + +struct LinearTimeNormal{R<:Real} <: AbstractMarkDistribution + L::R + H::R +end + +# There are five methods that can be used with an `AbstractMarkDistribution`: `mark_distribution`, +# `sample_mark`, `DensityInterface.densityof`, `Base.eltype` and `StatsAPI.fit`. The number of +# methods we need to implement will depend on what `mark_distribution` returns. + +# `mark_distribution` is used to return the distribution of the marks at a specific time `t` after +# some history `h`. Say `mark_distribution(md::AbstractMarkDistribution, t, h::History)` is of +# type `T`. +# - `rand(::AbstractRNG, ::T)` implemented ⇒ `sample_mark(md, t, h)` works +# - `eltype(::T)` implemented ⇒ `eltype(md)` works +# - `densityof(::T, m)` implemented ⇒ `densityof(md, t, h, m)` works +# - `fit` must always be implemented + +function PointProcesses.mark_distribution(md::LinearTimeNormal, t, h::History) + time_ratio = (t - min_time(h)) / duration(h) + return Normal(md.L + (md.H - md.L) * time_ratio, 1) +end + +# In our case `mark_distribution` returns a `Distribution`, so `sample_mark`, `eltype` and +# `densityof` are already take care of. To have all the functionality from the standard +# mark distributions, we would only need to implement `fit`. + +pp = TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)) + +println("TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)) = $(TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)))") +println("PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)) = $(PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)))") +println("log_intensity(pp, m, t, h) = $(log_intensity(pp, 0.0, 0.0, h))") +println("logdensityof(pp, h) = $(logdensityof(pp, h))") +println("simulate(pp, 0.0, 10.0) = $(simulate(pp, 0.0, 10.0))") + +# ## Multivariate Processes + +# Lastly we can define a multivariate process composed of multiple independent +# univariate processes. If a method is available for all individual processes, +# it will be also available for the multivariate processes. +# All methods now return a vector, one for each dimension (each univariate process), +# and we can add an integer argument `d` for accessing the value for a specific +# dimension. + +imp = IndependentMultivariateProcess([ + TwoStateModel(1.2, 2.0, 0.3, Normal()), + TwoStateModel(0.8, 2.0, 0.4, NoMarks()), + PoissonProcess(2.0), # Any type of process can be added +]) + +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))") + +# Lets end this tutorial with a nice plot of this process. + +λ(t) = ground_intensity(imp, t, sim) # Returns a vector with the 3 individual intensities + +events_plot = scatter() +for d in 1:ndims(sim) + scatter!(events_plot, event_times(sim, d), fill(d, nb_events(sim, d)), label=nothing, markersize=3) +end +plot!(events_plot, yaxis=false, xlabel="Time", xlim=(min_time(sim), max_time(sim)), ylim=(0.0, 3.2)) + +xs = LinRange(min_time(sim), max_time(sim), 1000) +intensities_plot = plot() +for d in 1:ndims(sim) + plot!(intensities_plot, xs, getindex.(λ.(xs), d), label="Dimension $d") +end +plot!(intensities_plot, xs, sum.(λ.(xs)), label="Total ground intensity") +plot!(intensities_plot, xaxis=false, xlim=(min_time(sim), max_time(sim)), legend=:topleft) + +plot( + intensities_plot, + events_plot, + title="Event times and ground intensities", + layout=grid(2, 1, heights=[0.7, 0.3]), +) diff --git a/docs/examples/Hawkes.jl b/docs/examples/Hawkes.jl index f5f90e7..649e910 100644 --- a/docs/examples/Hawkes.jl +++ b/docs/examples/Hawkes.jl @@ -100,7 +100,7 @@ n_days = length(unique(day)) h_day = History(sort(Float64.(tod)), 0.0, 1440.0) # build a "history" on [0, 1440] minutes nbins = 96 # 96 bins = 15-minute bins pp_day = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h_day, nbins, ) diff --git a/docs/examples/Inhomogeneous.jl b/docs/examples/Inhomogeneous.jl index 256e763..ccfc45b 100644 --- a/docs/examples/Inhomogeneous.jl +++ b/docs/examples/Inhomogeneous.jl @@ -70,7 +70,7 @@ scatter!( # The simplest approach is to bin the data and estimate a constant rate in each bin: pp_piecewise = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h, 20, # number of bins ) diff --git a/docs/src/api.md b/docs/src/api.md index 27d1f26..ad2454c 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -38,6 +38,12 @@ time_change split_into_chunks ``` + +### Interface + +```@docs +``` + ## Point processes ```@docs @@ -47,20 +53,24 @@ AbstractMultivariateProcess BoundedPointProcess IndependentMultivariateProcess Base.ndims(::AbstractPointProcess) +eltype ``` -### Intensity +### Mark Distributions ```@docs -intensity -ground_intensity -log_intensity +AbstractMarkDistribution +NoMarks ``` -### Marks +### Intensity ```@docs +intensity +ground_intensity +log_intensity mark_distribution +densityof ``` ### Simulation @@ -68,6 +78,7 @@ mark_distribution ```@docs simulate_ogata simulate +sample_mark ``` ### Inference diff --git a/docs/src/assets/selva2022.png b/docs/src/assets/selva2022.png new file mode 100644 index 0000000000000000000000000000000000000000..7d6391c85cfd145713eb80dc434b73d9b4507ddc GIT binary patch literal 70137 zcmeFZ^L7scA6dUsv4ev{YlZ#0*I|{d4}QVl z>W+7Y$Ch$`c5cm<%_Fk~;8@1Gb}m z1ZWrtf=gXbOf#jE#HFO9Y)k*V*J(*}905n7q~Br{lL@ciC9L(|>n`(B(!q<;DGMY+rnwW6Xz>PUS&1N<6e z28G=b5js@5xr9jSx-9NnU2v5;qtSbtt~Mg zBnQrJk5bYU%aTq~F40Ah&yvp2YVuUK>t#!H-OTCA!}@chd%zPC;4XcMTy!d>pLJUO z2CJ;JN=r*)nGKK`>-X2Q6zjBd9(GbhC{g|~>wU11Xk|JP z2M0EiQ2&`3^?FUl6xjeJEsv{g97YiVS?~Jd?1QwlWpduA<7XWi|`=WRF zg(H!((ElrnHjdYYc{ocZtFVxP>!LWTvF7dnuF$v2_^@x(6;8>!|Kq_|VxDL<%NpOT zWY;hFEQYe0n#qX@cb{|mS{%U(>_tIAp$8%WX}Z1Mn{nMLfc;L1UI(+yk}&_}Zwhg7 z@o+q|lacj1pxmbIKR;ORpzi1Wz|73t+12%omKNnB2n;z@)n7|(K^N!Cp>ZrmXh&5N?1QtFua?dt_ z-GXx|fCBdKDYyL#a6ot|L_?laP*Auy&`UqZGC{t9*$4UN&Mlf^7YqHBEsL+V#*4DD zvja?qGY_WBib`BhHg^u48#G}nG|5V8NB24x5(y!y(QuV0v)Eq;9$Ac(mlo^$6JO+| z5;UWWM;w7|x}RHh?t8;|zy4^<8lUJ3@TlJ9D3;TPT(2W!HOXmi(uDu+fV2w!8i(c- z7FurZJI~d-&}r00u6D=Bzkh%4qs8>`-aO>q1H^$MErF?Wlh&aO376aJvz)xV<)7Fx zD}mI>+CKv@yYi$>+S=O6t>#IJ#=G;18h@pWMV$9qd8Wk++`*!f9P_!nre z^O4a}wpDy!Z&6~0Cmlp~kHNt)`tkAJ@$bP0kbY4c8>Uc!8_s^dGa@!NHt5Y+6Kt!` z=;28ZTc`hktkf#!_J%_Ad#;JlZ=(!^ChMt21uV_5sUl1)y|FEBB)SnDJ&JtKn@h8y*D}x)7=OxPNyON2jCFEy(xsvr{A4p(><<{I21Fv?gS;7wMC$$xq~hn)dyy8? zl1bcAq9iV#o1m+V=2PV~ToOT9gN=WFOZj2QW zR$9m_DT!)mXxK(83ikf{OTftf=K)&y`1oCAPPG5{K*Ga=FE}{3tBXF_^zVUh6wg_Z z5fMQo2}KI+C{w!k{nsyXzyUxqu;xrmoJpnhb5Z)wM*Zo|<<-@DVzfNFdg53vKvs_9 zup$B2L;3VE7Mj0Gy)%>Ddf)0+QPUbYep5Q}AfU?(dvyVRx$ie4{nTYDEi3El>5;t2 zRup2p%MbciD!b3Vyx_f9iQv739w|rrT1=n2^Z)wwOS#IDxNvJIONMl+%1VX~ciLp2 zwNh#(x4JrZXh@EzW07b60!@bZ0RH;&M51p_w@sHRQ!=gvcI?|2J)7fYh+a~@|JR8sLjWxOus|1MI$h^Ta-zt8=A(SOEyIT-$v;LS;qS5Bb3-2OU1 z?bYqFe(C=}vcGwGQPLy$`@(7OQwLu(Yro?9o2X9ZK;j%48rpZ}1%`L$*Rnm-p&uNP zx(h-p5Kl&(-=56??1G+FH5m+i7wqctr|d2*rlWPKyh|S>Vj!AAN%+nu`#szHJx^Hx ztRyGikvVcVnhc*iASN&MqB$SF|2*C^r~mIkQUeERxt0NN)yT+5Sp|i+z}S65OaHk- z33k*1g`GDtVxICi(!tP~ki$~4Z^$%z?Qc3Sz$Jx()YB8o^a(_2oYTBV4B(n8bOk=9 zw|7fzG(Ko1{rL6k*UZ~#$GmfKo{WEUM3t(0`g0Z*Y)F3$Z{6rS`_sD}lrqwtVSJ#D zKosYyx7erepK5m!*&T*eYu%eb9-|I zJMDIBJ&S!jd-^xf3w?ckcRd#B=}H?Se}K*egoJ7i`$j|QV$gkmft4A%eyRU>9;sn+ z{PJW=r9}7b@I0+z{`d5B_>bDJ-o1T`fP;erFeflPJbYFfSTa35{d0SJy9IyB$@^i( zJ1VDy#CJZNA%R!xxHEO^)St|^0kX@_PecGUKBk3me*tGbp7i>+sJcxX4o^~HC;a8g-8LzXx=(2D z^YHMv*!P9hTPB?1j~ykC=X<+MrIpPH*gUqbJaPD3We_%8AdgA>b)fA8_itkU?zseRhnHaYUNyi{H48o{;SqP( z=17xaded2$GynG};UV9YljWx0F@86RW(Y3X;9}Wb{Lh(&HA5=Q#t(tf56_xzlVIdO zH8pkRF4oDzzcfSW)hb`P9O~oJDSu0Z{{bMs2Bc%y!P9I1Y7#%^vQC^jXl#3``1aq_ zW0j%4Z}cmY<|`pb2q4vE?iYtUL*V#19}kMUT&%`S*1NEP^dgnOesDT%WjF6}#^iMa z-8rnHGCjRp?G-ry3c%zf*1j-$YNZJpRs>~wa8C_ z3Oj+@kstQ`a27h;Cq&ex5qTMe@8#u@z<#Xs!>Z3kuH^pBqH2I50vHBx4T43;AC0(D zg7)RhR6w`u>v<8;NhDQORpnLJ`;tHqgNtCn*L^1`5cSEuJdgN2!ig1Vj`tH(qOe-? z%NJyDx<%(3J*rfGr9YUcoPEh7BJEB+`ft?1NF2C3i2v-uLc_K%I8sMS`O3{-@mb)W za$1vt^l2uJ0)IbWqb}N*Rou~$JuWUz%H)pD`_}>P^7!eikD3jxm;imRPsUvqeeq(M z^bkQVeU5@lN5;dG5Jn+tIgU5b;)`+tf(q{UU>Xvi>+ytPaDo57@$V4^R?_?liN*EU z7*uL7wy;~bX9%zw3Zk%W@62u`OqAL6_;)7MY%G5fWN@90p$t%7+XB`rr>tB(%dHkP zHj#Ljn$YUO(pFYi4fyMI;rlw=nJ^w&e$YD_P;4eDH#N9s7$iU3y0qW-k`dEudMv;>U3vgwvK~3$sUVi z#`*!Crt_ks$7Om(#^A)nyE?}m7`NmEVO6Q7WLDyN?HWy|@&*N-TvtB@g{moyR7GHYt{m$)(NwqpOeL&k@) zB7LAX0B&p*<4yWeI*FS}`^QH`!gDcDl;4#&jB?!PclbnpONij53fEoAI#T!N=x;|v z?jCjQoyn`JiUR=vs5#S>P6Tj(EL(F*f*Qc%>XZ@yB4DBH7E`QXs|bCZ%a0lAwo6y^ z?@&xbylX*5MFk(01zC*x^q1R%M@B~6_4cWsAlr<*7q!_csIufBayR?_SApEli$uA3 zf`h{Zs&GyBgSUX)g-1ATby9?`04IaPBmGX=bnwaFkoACtY;7$!iR2U54-;zO$^`VMjzGiUovK3h*l6A8(LS>;T$OQwQB`8%A*& zZ!XO|Px7>w+)mAzSy;H1zp`+i&Dg=5YcW;cynL|DdLChX8U=yKG`ja$Ae-OJ;RdH~H|3E0H{EeL68*hcD>IRWW% z=Oj3=C@(;Hx|1!Nza2ItzsaTtfccXn@_w-_PjKa z=D%23NyoJSRdxOGkVK+&%OiD{RZt7vp;TQYh!~(e+JKx=V!t8xzf`QYXkDGN*1ypF z2@rhNf4MWPuF9&cj0Qs8^4V2_kAPdLOQ4}vK}A35|9z3Dqx8StRDyYP8%S9^JHJ)` zkeZc*^6CYd=U9_E*-m%A~c%#6%*X*TL;yxGP7GMtU7A z{gYqEwY-vK1g`B9*eXZ@_wxS3f`KcrA0x)GnqYZBmq}sa-s*?C45;fms6w#QL_;lq zRaxvCt#pQceahqm=o#P#leG?J^6_^hW03)z7odBs0AVwwwl*FncOsEEf8^P#j{sO5 zFxTM5u{vFe*fWqKgwLc44`()BXaK0s|C4#>)GEI!IJg>ME$Z5gg9%NX$P(0!TLioCa#_jDR{E4agaYbf9HHJ#QV8 z1n=O1Xcuf(Z<@W*0z0V9nuwaqC1!JL3-CQqibX_6M<>oKsYL&$MvwuQ4;O0ifaKp1 zO|KR!;B`ZSz>}MorveoFLE}gO&4WOy21Q3pWs>^cr)=lI5-BlvB*=gu>R`6rYrJBB z1LY+693l7YGf*4D0l^I9ECHD8c#jjPI{2K{3+!g*8mG^nQr&@k0k{!NA^~SI7LAy@ z7%0fQwOb4@q=Z$#fDhbzK})*~-d+Gw)KaqB?&3iQX%?XO!02`$dCSSleyg9tWqiLS z`fqrvMJMQD`T{-yJ8)FeiSHB?aQGpi$$V~d8X8gt29!I~Rav5;B*31sA3ngp2J%() z_389QmyG8dX=&-hmumGuC_OyvW@NGdEnC=dI`!Y20pu1qaD(-jCMYFHvg<%HY_=g9 z%U2NscmDZF2`W~~uZW4F!>jsPyb z>Q7*QV`JkVHbl0zw!eomP=J{7+$hO#pB{fPz z1C?lQ_djKRSi>0kUo;~xwR!?I5&SO~4MM*WCdY&k=k{R2j-!478Q4Q1pxvZ-{GoWA zjqiXQ4*VDQO6NYT*a44ROl%LORifM0AzpnHp}Z}C(HH?f4EZZy}SkP0Y>HK`!O$zFFkNc=+0Zjg|E^&|{ZBzW49Hi>M+H(GR|X zCmN29j)&W`zxw*xw@!uz23~u3c-S7^yjMdZ^?cZLj-c+vH{(KGe@61*!}Ik!ax8`q zvOrg7=Z~!u-D!51?9^5F;hzGY5-Sh`**84mEc6|`xWwYH%;HYfV$-d!KdGORmTzh4 z$_&!yf4wqy_?1+^)BUjLQwpHZ0GrE9hN;9QBml0^0-C7}F#Y`66Y=rW*NTd`K;G5` z%mdV~CUrUB^6+r1M(Cl$%@O1lM5i!d(6A~A$my>SJFfSLIwVz81_}~#4$5>9SlQp4 z?G6-?5`CABlaP>5X>{j8KtKREM*de6(T@Sez5?~dc+*XrybPbOfeJ{E<69X3$6!y< zg4%uq)f~cmPjcH2%Z)h=IIg<#nazc~&=1**|A= zad9!@aW)G*jA@$5qtd%F@4~8HCy8Z7kmi2W9_l}s7Gbx82lTcjkc-ot_U8Tw`{g1m z0-q63QtXU4PvL(L3UK)cbn^26HOdQ=oY0GYFU#DE<)e%2+}x?cH5p(w0DC+I7ZKl* zlchm8(ZS(_*r=R^ve7vBotO#|x+)r87KC)hm!+U{81a06^1>_%M zU>i`qP`l9n<1k(Y={w!?oDPR)Pi#%hS#J@G>G0HpG-77Z3ubxc(*3yhnkkz6!#Y|i z1$>h;*@S*soiUV;>_|IIPJWAfL=sw9QR)O_yDZ>cR6=(1@olry?MdF_0U>}u)_`Ol z?pS%bWi>;77OzJeupxNW+HaIQf$p4{4WAGI8AIK$!9r7BV{OsD2p#9NcJbR zt3F!H`9@dI((5G9qwLlXFHM84OB3g5DpOB-$#Hxc7Bgp53&BX@zZ>h$5b*A zLm-O+a01sibv%);AIthmH|B2%8KlDyj~1F%QGVNUm-8UIt8~?*c@ z;p7}VJjK_>vZa|(o`U{U-6eQ9>3N!#F5eRr0Dl>VAFaMoXV6KVcSd^7Ts65m|REid1?nV95W-1 zNWQpr&dkv;1D#2#%T1~|=?Qj4i8?__plbF5l2NnhN0^ZbwRt#aV zPlsM~1zj0KM|0^;<_XVY@da=PGeKg z-S%l+WfK!AgqGNOtFgT!9Ti;Z=RSmi+8G%-9ZjFGAG5!GWL`9GF3M@%XkfH6{XH{ zh|fTBa9>Z$*{^@hBS@Q7uZZ*#``!~(<9+G_uW`1wKMu}R)g35TTp#=-AlP47Tk7gV zBvV!J-QsjEEh}MRv1!rK&D8!_^&PHzi*jt>dTFV{NtJqsfbAqN4a z-21_DR7Y2CX|mVlmAFLy<$X!bm4Ra;7bfZyStpP_h!^zM z`Jh|ExLKiASngZLgEMH~(iFAu*o#56`LnO!Ei+UTpKkDPD_7s?;5Dgz}5*WO6&83!U@9U#F%P)Jk zH4?ZG&&JvpkMO$LAV$ZnVZ`I`kv9%)C*OrWvp*_@qWmYXIbMSG{;cazudh`WYP!L} zBaKIB)eKzJkgr9>-}c*sKk%LrS@iM6^yst}GqacgY`fnb6T#=s$wi&~@byW871m^J z)XcCnh26sI7PZO+P;e+K;OrFEW5M0F_>OoTAmDpkg}YsDU<6?*SA5xw(M$#BwTpM@ z;l6V=ggsdHJtkLHviRC}!RizYHBzs{uig*3_!r@axSgnqbrBqfGZ(ie!Yc7J;9p-G zK`@nT!t-@nJGmC&$yZlLJT7TcM)NkTuTRzpl3rj}yuMR$vVlkg6r4_QW`BXwRZ>t; z&^CGmBRF1XE=%eYyTx}cf#Lm{LbMQGEZ38%B9Gh1=;%`2Hsr&jBQDE{S3X-~i|r^y z%D5%*1T7gEncri=ys?9^RGX5O^(dmpgLaNEQQ#lq;bEJZ{{(!g)Pg)JO1S%%b*m2V zrff z@GM&s*fHrDIea(2K2Efb{qY578b_TXD(W+bMa#j}!P0i9^?Qa)(1D4Cm8tYSSv3!1gzTV$~uzhS_c2l7C2wd~|jq}kV+(SAo>GIrP zPn?H`1II=&!A?4(>CdOCkjPZ!g{*f}6Wc|$6NzSwBlp*O;c*$;B;e4Wo12o&wDiI6-XZ|Y2YjY*Z{UT5X4ZaCQu!=|IhA4cyrWRtPJ!RE zNs3H+F5Lu3^}{LgK7vnQbuLvms+m-6*b%HTkfmt}&=tjiZH{>E-vdx$GLA<5l&aXB zl1)Byq3lbEpQg`NSBB+g2Utr>U(e5o7;q-SzcgK=s^iH(o0*H9nNA=cBKeik$0U4k z#^;9;n^{t=;xZl6jeLVss!$Z%f&Eb>{l`QPqu&}|yPD2>L6qMj%z>j>8UP+K9h@~K z^qUPMb}y&c?5b~O$#ev_?K^Dv5C*-!fF22pi|?zBOBls$CBf6qA-Mn_Oq|&xtzk=>nsSCa`{86j*yAT$?C;d ztGDM!(t2r=A?2)4(3W#FrZxeJJ{;WpvPW}n=j*+#00+lQC>k+K?G#|Q(ZOz0f3q4>fYW(8We1E9c zfdw7)V7^1ZDw(F(t%0>C?$d1phl*BwvylYPV6_YCJjo&MbzMW zgUVapfO8;wJ!Q7$M8f@qhvHMaN?A5|x|tdk1>6*XjW*2b#DLS!A5>&#RKGKgomZH- zyBwNxdo0iFkJc6u^;ZAL@MW?`w6OwCqjXYnS?R+<-K1cwRfyU(`4w z*Oh9Qw7KoQ@yiO!6S?sp_Uk9ZJ!ST4Wirp-4GGxRibhE5E9c6cR^0a0TF|PAuYaFb zJT&w$CQ|nLa!VQb^K2)Gnl@;l*5O0LIXOO`TlicSIeKUQ!K z-)I4(9uxEQ(}uGE#Sh=+ib&^t?zbeczrUEV(?=r^Rt+TqshidFnl0FFeQ9gJyGQU7 zT=w<8(3?D^2yperY%L}|ZEKVEy~QQn!=*G+q-grI#j^flJKIhI?k60$6rgSzr7imA z(Zh6(j3R}N%pFj!lwH^!7E`$z~fM(O8VJWaA{NBW=h_Z6E?2^{3J zK*GUMk1z?8x*3N#z4)bH1rR>QQf^HAVg#Z5@IikrzFl`j1|r?7W`(IWy_WbcsEkunl-r02s;R!4^z1|a+718b)M-)ybxSuDscX>bFCLn zxY9gjw-r6d<&0cT}U*H@tC*?n#STA$Z9BpEMj|DFTn!Ow9oJgcDoxLcep_BkaO;C`jz?!a z<<~h@z}golWW1wlO_ZpOLH+c3xgn|p#M5-ww=igU;+=IHGBgO6ZI*P_I1(VpLmZqq zw}NF18u|>IvHZ*{_}2id?@PqQ<30a~{SZUFCcNDjrPt{`xq$(EG~F}c=3T9lr-P-; z1cv?W-EF$n&1B$2+?R-6Xy4nhB>U*qg>~ffoS#uq$sge3(Uox8LhWY`%jea$lmMa2 z-iJMG9?g5ipk=w6@^N_JHAtL3v`X<)>aWY|6ZgMlWJU)NW(Zm0Kj6jxD0v6t!gZw-i%H zM-5ZXneb?6nuEs&$O$MF<%R6G6w9=hx4uNiU>D?=OpRN<5Z3R;mDl{DB>x3a&4PS$ zZp(nS00jliyj)Z6`C03Uk>i{k6V<+hE2)B%a~JbzToQEsO9K2!H{1=yfX1EiAtonU zewinbesYRWiUJ0P@y%Xr#OE|TgtoS-w~@(>N2=dQ(4Fb*MOP{sVy^0?TgQBPhV-Mam-6Z}w(-&Rp5FS(wtH4N8O`W6?JDKw5qRp!d-`-_ zytSaJrK-uo^giF|(zek)7iELjl1N3!A4|2&Z*yFW zy4;gbqdkf;n|ZFupm$oC6&6OEriovPIA&UGI+m#y%)h||Ki4=JKN>jA(f5I{UU+Yn_}gslND9kT|7tFKh}CnW27QWK5;J+f7~3AAfY|U3{ptgScZgrbDqtq zxzF3T`k?s~|N06$uZX(bk$!PVZVI#nX-|Ka5Z~{IgjQ+rl!kC4cdt{ZDvI>-?0xfC z3e-Q|2}%mnxlc{V{%x?117FDR=B9h}yi`-9^M5`PKPZmCYEPuBX68gImO0S(eIpIOR_Q!Y;1WC%kcS@Z$v(^_IlR z>=vE;hwq3wp2^AicrCsGi~luZi0uCe2M1I()NVf?^18Ty^1&DoLYe7G{o{>KNfAQj zsrPT|UCwLmpXn{UHiw|!4rPFA-fhh#9q-QC_y__fzaT=ck4ETMB}FK2=P=7XV;!3a zLf4)D&DQv``OBrU6Azc&ipNt5(1a_-+R*xZS8`0!d9O=@%p$I&cG*)+j%!o;dd>44 z=OBpY)=-Eoc>>4s)>v?P2ywH+BV3$C5%P1O-t<0)dIClzNGIYt-~Y+&>AQbI|2gzzW+BN@n%TFBIm+c>sb>z zi%|X@h3H#Z)YnQVPc4 z7i2mJy47^0vL+*tCWiBA_=o|J0nP++KF3p51ixc@P}(SKB+;aoTUo{PAbgTSz>;yl zN_zd{3F1_x#U}bQIgMb#V4&!yz)MfeB!4&j11yb z*_+8qUVkii+##Fue*;Ka_i}SYlEo!gfN;PEg~hHelg%r#v55d|WSy)GHVZRQ4ztp9 z|7wwzXno{w1NBtzXdbUDxZpoc0N3m9A@t@$@Pdbrb5Cf$YYV( z9~zcYyRv@D@$@@?Y&+y8zE97~K0D{Y+Ww5|>s23xJ>$er(z%9DGGmFM7PRLSrA_5O zTNayfg9JJ2^$bQ&cNQ|4aBugN*by$~_l0VoWvjq5Cc%l$BYg5g<64tSDAy!)e++3j zr@X!|vmY_6$-2i-Z#^A^zTW~h%RY8$sg%zOa3L#1Dv32YX)0UHw`!bR@*^}4;aq}q zw?ow{un=f--exLG-YW7BVl?avr>wuZ+_K-9iXvRD&`V#< z$+5XOieL!68BED+UkgvuTwAk)eAm#VWbAG9W+b#C7Z~R=pO5^qZB^OejsOIU{m-#2c0>xFMfh)fv#Dn z`Cc$V5S6&rsCzrs$7`vlg!P*sP#m;?B#c=NJ7cLP*Y+Cg)>SWE@p4 z#EZo_Cw6^Y)|VUBN4?g^gVUa04$?@Q^Hr?-`K&KidsK>FH$xP^XST1B-dbN=?>(DX zSf4v;KU$W*4PLF7SgA>b9+6}%!2P6uE9?lGyERDMdUQ5NcGj=lt%s#>m)my)_eIep zP}Y%5&Nusy^F2RXT2#fK;qu2eo^yV26X+dyFFKs3w8ig;%jkGyISKes>7?Y26GxzD))&c)UTiGi3(Yu?2TcOE9`GA!%hu-JYM7*YxoP zz&Th5YTX>h^j>7V%_(pIO}!esb^WOgR`}?id(w?=D?C5d(Y>LA*6ji&!x?6Wy5!_Y zJ?oGc>R#;eo{YYL>wGvxXv@<%%0s@k7T##CwR{%JsYtS2vnMhPoqD6G{CwC36Najjc}u{BuB z^6o1zdgGD^Vj|IDOa(NuGr3e;iAtgN#Ohow)6jj6keO?gA%&oGJNd~=GhJWwBtAM= z97qsO%UqLTMdovTU7vM8Rr1t)7CC_ZVKpQrQ68kD; zc92CYR$ zxqc$Qf7yFV^#zSdGQ0UJ2zbctaTR~U8tnR=Ox8@>O_JKrgC{@3LYi$tdwD7FL%8CJ zLVXjYwf;E`IeAX6L+8VMGKz9fe;@t$@sxtHkCIM#yGH)r_PL#QV+$3P^W!8ar4w%| z>ul;S$jNPgKc6o{&jE8~rCf7cu&%eRPw4!epqr{$j5_J(5AO#X<3*$o5FbzIY|b|c zfIiYxr6NszeWH(!Z4M4V1E*7k-`L1BQ5rAgIVFh=(-FDhUV4bV5GOD{a%`g^?^RSe z#!8IhyO0*TP~^qd+D@G*aEqVp)Z)wZt!9sIgC=nicz~HBWZEQXcwJJ+b(*9@Wgr`FkdlGxB;L~WHR50Wne5JW}6m5#0DE5;}sGa<<(e%bqe{ZeZnO}YVip)3M9t3q1 zZvB{nqSlv(!)DD)&7_HBBZZ&Z<*XInVPS?n87sy#ReY0e+Yz3Q%4+F-srzu8no>Lu z3-xs?Rkm8y(K4AbT7}s&+2qI<SLjkF-+WtyDq-13!VtS||mo1fa~ z4DyU#?lEV>82!t=HcXU*f%{ODvBBW5>ujDhRAg|jo&RS9qe4V4nnOlD%~G3v0Y<@jeJ)N z-TvW_Vu=I)GF<*K?M-HyE8B?Gz;%9wt%m%S*%mL$H?hMC@&hN2N0I$QHe}@C_hJ1{ zAnSu3AJxXQhP?yubY&kj*OpSbxYtF$#P}F|AiW0^#CgH#!z(HI;&u^iXjjh+zlJJr zw=GTQ0ZVqbcqv~n&*`U)OP$$-h56h}EpqbkRHAH^2aB`y-^;vzxyaUFq#7P8npmh0 ze`wHI-i1(+k7jPfq^b4-hfeJ>j3a5#;iiF}F?ElVQ80WOEpWXd0~_{)jaR`YKhQzJ z3YxGE7yU_hUrB*xaxgKPE^U#b-%Msc1#dAOpQ*t=AF>g+AX{abdGy<&x-bt7qQ(iK ze$y1PwV_!YXfI=YC5&>ETi*CIS1C4k2w=OnC4H$7||i-Jc)_# zOQyJ)B}I~cd_PLIjGPM$&Z3oqMl$o+TF{~A;tT5q1dVnUF^XXD6--DfSDNE@$1u)Y zs)84=-aODaJ~E!Q(1}%O=EFBz+l+6tJphrBCxM_#^I~2ePd+<^E{*vs}SR%>YOQ(A*j$aduD@b&o`ho#9V=)+Lez#PH( zw7h%6Y3(h`_kmtU8V57%6CG$F1br{C7Q+Z>?m+OAGZgmE4h&_2p<~#qJ}&>VKJK6p zVjilZ(}6NN?DLDseo3T($;mD*$*h4a2%27dtjWa0;I{Ics$jM{{K?z3Ju7Zp^59JK zWo~~L2O7*Up>A9G7pVgi6P=(T2sD*H0Rj*SA}?JxIYDqZ=%WISG+o99QkiRgN!?%~ zPP!w6xCS(eheu>kJ{Uf>q^Hj=oHftCgTAqG;V`M{u`uXOh9_W*lt%0rUO!_ybquoQs~FTTzK= zZ3!$Lt6kK@phKCsJC|ct-)1?OU6iP#jx@(Nm(j1t(!Zi8{>R*|CMC~D>Id8;59VDS z3{Cqa@iNp0n0MQNUKwH%_W_eloac*QSr(Gr%teByW&^29?PwfE^o@-EY0tM^d3SlX z_YY9L(*D6A4Q5o0z(g46w$B?y@C-nz8rjbC4Pt%0${fJmMbk_7;{m*>cDJNMI4o;Y z;QD>Jz%UPJDmX`E$9NR~Y`)uwI!jzav9A5`A_?jqKis=|b8M717MPK}?1I7|^FIez zLs6GG*r@tyoUIwMxtdvDU#2QBZR`ET7h@uIbCBHj^#&G&+fRpQoBl&4RNq@UW;!ZY z=XG~YxF=LYYb!7ULJO~?It0o22H*(iyp3OX_ylx5OI--o(oAd5VAmaBDc!IH=NqnOR-pturR+q?NfS)Up996dB2{fwrdTi$r zN&x$tWs%^^QPt08Qd#qhn+-f@GY>}!;cpeKNA6Vzbfy$jho&O^gbh?IliqfMkOvdV zHnxW(K_Y3eSr-^^2;3Zmf~gI5E{VMF;OG_0LJN$H(M3f|UhVwe;)o^t7SMTbc~|9@ z-qsm7WRy8zclAZ2qI={kuPBqt8DwHd`dD_adu(RCK?Qs?!>l&Xlu;`3B(dLs#r}rP zOM(fyWtsT5)Z&g$Yv((XR@o6_-45O<#3ObXOT76+-gOxv(YNS*H9it@zi`HLhS z?t4xnSt$XX&6Ro6=l7~CSTc=AWVd#&1kVHug}V zd-=Ly+hCYU4Aqiuj#LMelZX96oRv&NX*A+e7A>SMBFG;pTirOmWveYZpi?l$olxL6 z2&R8#2f_H!sYxxYzZOhg3*N&LdwWAKTZfjKLX0dv-)WPQ@?`U@L52=9%qsQE(^9pR zm)RPqi^Nw2K??>R4l(mpuT!86Y}bkO6}#uo^X# z_wXSY{8`bI0;AeP#zDe$@}>FD8+ zxCcTz;@O&L!b?&YaFEaf6KduHN1`KklGY*-fo3}njqW|h+Jvv=7?7yjNE~&9nwm9E zvYVCnRcM>+ubnxFMS}u5J*y(9`#$A?DNA*SQ3cpvm;m#rd~R&B;IE#*UV#Cbmo60r zUkp?VYQ@y@6l$LS!b`~A?J|l-YE<>QkM6Ng*e+8OgB}fO&_g9d2@v^S`>@h83h?-$ zHy?WZFHgGcGAQ@0n9w%{mF*&upu4ZOsdV|8MClqZx_u-e5NNsT5U&o=xDaa<_=6g^ zNY{lQzlIzjd5zDYn<@S@yDfJwbW<+ZW=V5Xp!c_~=Gy;+pdmxeeQXHb=_7yAMJABu z*XN>U9HS73aL$5+z3oK*xRHd{q}Th0HhpBTZ*A&|3@AVzytoKE{w{W`v&@N>GSJv` z?~5?Gzmz=V$<=VnNFlRTmAcVYyntqGf^%Zq*&kHxp{#n5T!Q8(ea?y{34+zTds z5#6}g_#YUUwlFbGOr}Uti9;olY3T*Y=B$z)aFTTHJQAJ1Kg;mr*Ka=qmwPDV70nga zCHGLKPjyyNeLR&~K6XFL3v$^H<<6XswYE5<7JCUW{)E|%=;tf%#NBzHd*g5X2fvLW z!r`H)TWfSfrz3o$IN)715k=J>RgOnm=`zU&y%W>{EDJWq3IkaSLP-Caz{xGxQ}g!? z_y0rIS%yXReQzH@LK;Cp8cCJzZcwCDq?B%?yM}IQLAq5?knS8(kRe23Kw@a589D}d z_I!W;ch3thxddj;?6c3>d#|PZztTigN zkFAG#!JUkQBtu@gDrs{Pu6UNxnexzti2>uoM8>*B^5w3Kf|Al|aE!7#d3j{~qH2@~ zmg#x?@BP;UTTTm>J0A$rKH)e<^^obw->NZv?wo^*1Q8QuGf`?9zj&RzBKJ*XVU?gN z;I}6Q;U&L^I98iuwPC~Nm^jk8%CPDO;Q;O3uppdk-iSJ>6l3Te8qv%Zb+DuTph$^& zDEyuR9i_bITyYg62~{;0ECA^Sod%h)N=LTu2PAfPccT^ISpoHYUHG<~kuLF%v6>h4 zWwyp6k%_z&8}n=%$@#`G5-lLQSGkVA$(^eDC4&_?N0_UHp_Sm#C=-1s@X0}+LG;?of)oZYIaX8%TS$8L+{JuLgH+iXE_DST4!e?cd z?qVFbT^hk~1NgW)Q^{;*(hRfs*~mpr;dQmkYGq-+c+w1I=o)*~x3{L?&@s$Ra6?Ed z5@&LI2{!~tpyQnxL2=(i-*)2>*;@z@Y7798#|uP9KuT>Le#Lg_jGo$SC^3zQD0Jr7 zaDN>xr@29vnaKCq)|8*cG3!qIa1X<@yPp4L>D7{1iAIw{peur@UOI05Z|d5>*g@U5b%?)^+YNsmDViyNh4qXND!@i~rsgZf}x)6S!t6)lvkXzDVb8 z&|As))5&fn#i|Q6+d12tk!|wI!@vg&u)e1M;!3!OhVRJ`cfTFSnlNYxQem!0QdtrH;dnA z_r;dq56lFNL;CCQJ^K@l(;jYN2S_6Cy@d$deonB?7zN7_w?7Yi7+$pDU2qoh>$~Lh zqi|*-NC#d6$aTUc;(~AVGH-F8q7Rv?7}N8$-@Z_vT#PHc;SPy>`_JPMop0D@6#hqzqsg-PTI4Q!H`H1dk;yp9ni_)8R3VUrcK|1nlJ zr;C;V+8@T7xWXtHIMgzjXDZ$*Ln(#qD3HjnJG;gOTMn6-Ul!k@cd~tiH%i%l7vmCM zCVy{Fnf(4Z%5i_S@CiQ!#O(l=shr3|92@7B`Y%anJY&{kB4!BK?i5TLlYf~Bd)59v z>|zKZ7D|$q1G^jyRLV8Nj%`>{#Y!aMdZ4JCwJWGtYeT_W;?Q4X*j+7TT_Kry3%OyD z(Rbr{w=LP9-MnAiwhRkB#yqfO1@nL4PC=^^K_S?vWAyW1v-;GkFP5A&CmBt@X{!IR zzDN=rAGYRNB~46Oy-u9pXk}1rI2Du+(SvDq9tw(K14%qE+jOHC$K4-9wRK*U#J&$tJVp>NW_O;9UkJYc%Jb;s3m^cV zG2wz8n^naI5^@nppdOf-Sw3DfvX8l^@&|Ll@&-jPGk&r63 zIDrUEMw zs$b$?Cu;{@5V&0M^CDZH?yp4(<`p_$6`fl|4pxE_t^&{&g#8ddx44sc{zV()Ro^n; zX4?)Cn4LO&;kBy=^8`TG6JiDhT^{ONJQ%RfiI{Lt}*gtZ;HKQ`#!}qwvlLM;4xIbn6Ys5a3+W{!Z#aFnHxJ$iPehr zGJ|U7@q8u7g1n^wf^j^elLZOZtboj*R{e(9U4DyVv4r>6Ag2i{DOg%kR;D6LPLoGp z#<%8mt)eQ+5c{PTO3nY|HFkD}T(CQ{ecg8W*mT5qn1m-6bj@$GgR-HMCyCL|pS(h| zlWHrcx+YN7UYr!Ia{v1)Gea≧KANy@r*07E%2D^)V90bZ}iyCXnj=vZzH+WNpIz zuIM&SX*l1p*)%J_Q5JCdSd}yX(hJ`})+=NqHwiQ?9F!H$plL4G%`wE=zLxUw~z&&|(Xz}d}Cy1PWK2k__ z(-lX_{LG2iVMt0J>Utw*pVoA||@}k@Hpg9nhen&uN>Ro2YeBsH$cfczT)dhwKN* z{hc@xo%2tyG9Y~)-DP>GM8``ay}C5`lpBbK=KDu}81MvLrt}2e z)LOi{T+ilqY8$ie44%5C{b96yI2t4&br76=tU0rUI^uR}?lV8OZMUbKp-sRVmDvXA zV`BR$hkUP}4?$s4QD^~v4p3G7Qsw>& z@*}=0@uC9*11}Fsu-%QKRhEik&DFEl?-O7TjQ;4c=*MflAQYE~@V{H062{QCd;URDjDd~ANLTRINdNt zTULLqNcx!NG!f!1?@sV1N7B!mN-q(8Qqu zZ-fs8?qUE1Z8RB6)9GLP=Ix5sCbW_k{fo|lN3<0!TTc>XE^QrV*(qcB-%oI&i)J91 zmyb$*^w{F6aAE%ZIkmM=-1TqF*I?sNaOnGJ)n;-`1e&!|r?S)Ikm224$1mj-m**(T zS9SK{((A3FlVAJgH3*jcXi5;ERs{k9CR*?A3Z^PLO)I1UvfnE}%v@4f7&4M8=?-*z z0zk=*)+z(`hYVV`yw`E%bOXpm2h$db{NQIm9s3MkVf4MUv451A|t5C1MdHZ%dQr;Mj;2KPJjsD8o64Hf6FtjS>^h z$}Uc^9k?*MuV26ZhoJ*ku@1z@g#`tst&lXUF;G5NiKpq4KAYGc1}-gmr#6c%00 zC!pf71cG5Skckzr>BM2gXygO{grNbCM;oIyK7Oq{d?QlhF#+N053F`GbKe4P~8B-2|e5Tgsy51 ze9SFNE`Dj0_1ip)?nBcwk|`~z&bVMNCfq*tq_EfpW$<`#$ z>I&J6&`nhssrmUy7*~plH%|hd2ne@)IsXo3`UF0IZK5Rpzd;7qq@;9O3oZv;8Na!{ z0(w&#v@+0oG%q*DN)o_L&=f+BJfBb?83Y8gO`wHy6G=|}BlTBLkSbby%ST7)wc;cF zTuAwS)#q~+T!s#kI1E&?A2zvvn!)-kl)JDR=zA0tGeB)I3%m{J|H=l_DE;`@11*sR z_=g@v66xO7{SH_tHi3S-pD-X=p`ohkc~Y5ZvHdpr3t*0?NrZ&Ps7>50KgTV;C&V8`I<*luuXMA~yi$cbm94U`I4llCP zpISoxuI&82Io{vRMy+b|cWVr_u~L>%vxxb^zBTykUIADQR0ow`w0Lg+MDwsf6rh&! z;Y(PldO!K3SSug!!P|8Ob>urWx&zl!XZ!GCM>ecy)6)q}4qa!Rb9zB(M(yaU-OD?onA@mY@fPtMOMlLI&ktTyC! zWJB-xqOh|QR+k!`*Y@=GtW?JAzW7>*C6qpo1(7>B+jumIhd}%!+RnqimLc)$DJ+U5 zo7%B?DA!F*_2_O2dZ*5eMgf?Ez1>uo$NA3dHD}ne#xo!$!5iEnsOHuK2p3i> z+i%`8^zeEYqu4#{?55@3Teqt|64&>L(dTE8EAtqzrOr5N+vHD_t3{uU*Fs=>tCz2W z3bn3IiuDxrKoR44>GO+4fF%4{Q}YqctJT$QfY||zLTY+C08*ytF@Jn5#XtbVgxAN2 z*@atdidEpkbpR>I?O-#D7aja_)OAlPw$BcALw}b=e-)SVCorou z2C;(xWz^?%#^xbFjG^Z%dZvIb7Bv4B9aI3Ec%M%Ld_!ddNk#+90Jd&b=ff9@@f}!^ zqwa&A7VhPL#_sx1FvQg z0szp03Y5hFlnbyljRS@=c(5AAzXSN3$NoYt8e%%xKm!Wso(DdiH4X_oEg2^HlF10; zV}H#}!7!I;MXusgkRSi!_Ade8%l-(UNr_;+YXgFLje2jevW@~({3yUS)WDV&IMg$K zr;0euiGd!y=+#MZ<=RKMhXPyw7R*@M@kM`&wBo-6muH&B|Ap9p)nzKA~3^LJ9e2#aD8RhI&VSs8&T023+uYd-3Ry{4fX6?S$UekL-G0bs-p9QAB4&j@zi6PGCrqkxE zG2()~xlo-Pbb9pY>eHv!eIzru752q1xpKKZBYr&HS9C1E=%4J6Mm2%WJFJ0lDU1za zHWiut=86D=st*H<|FakKh4%LLzF?!KUuAj+ps~yV_Zqe54Yk;rt_IM!;?3J%<6yNz zKd0cqcNw=e!;RYw8S+|$VFqFU2%w^r!J?HF>S9iI21m8#UOJMl=QuadDf#>DB_LDG zfU3)#P&>>TRk@ppHHkq^rpgT_5gMBeYDuddS9>dZ@9w{ zES!Ayea|}#eeNqiXup%L{|4X+S%v>~X?StCR@QJw5juzO95cxzezrbvjY3?kEY}-> zgaIUvI&Dp4l9(LqnUcVsiDO~b5$un^!toU_^8KYe0RR@@RT#qym$(OP34Eo^7KMLr z&&u6x;0h<&XS1Q433*J|!o?;25f{;#a_he$Z%q-8L-*h#2!+_y|K^?OMRWN>ptt95 zM;?krs5^6tr_}jTZZ7bZoLr@DKa=hAj)KaYEv)-gz66~yPe2hIBuxj<(S-5#{26-v z15i;wiRl*AjDeuH*wHKr&%tP<_{3YxSP$CbVuvxc;tHkq6n`I%Up!+0$;_5>$U2{! zVeHEZHNIt*^^EH>8L$d@vbG#Qn0e?9p|FIzc)TU>2)&Ksw|S%2;3Ha zd!J86Xu0*pr2Ly{pst$O24{JEH12h_!Lj*YQSA8|Z_qDv>_LCAr5hf`|NoK5e{tH_ zR+dmA@oovj;6_$5)KIz$R&fD|(abMDe{9PR*Fi56Zz*YUK?hF5Xkz~ddR>j&)1%f7 z%HbNLnsM6AByCM&6?$XIQjD7D=O55-AG`BAhvJX2~Tgbo^U7JZy~Kh6Mcv!uyd?THm5=Hl{5*& zXLRS57+o^K}w>xGP zlKh*&y)ikrz1tW7|9xUrj}OwP4I2Aiv`vc{j}5(+VT3$}k=uWJ&7Y{Fw;#oC(V)NY z93Wn2i#o4wmgJ-1iQt)x1c>{@7D>>7ED&(%fhOiJ z`0uw%=b@=Z=nW=Vc$EQR|L{EtAqRdoG`yUk3x201A(H1aGmgW}lZ9jZTd806b(K?W zV;iF*-@s}^w_rGZ$TrkQpC;439sBhQsdv<4ez8(Ir~9v%l*2$ceg(oYI69+2^=O7> z7%s`erjX?TDS-Mt2NNti4F=HT;AptkO~CsF;HGEGkAlA=@Umq>w$_qn$uF-lE^jd) zau-~n@ecu-e737v7ZU+hXuyp3v$SPh#E**ay6VYB2h;oqq5Xeo<9ae9$}|JP8|%OL zR0^0e-SMIe=-nimOi#|PKD{XfKuYK*3>^QS0fz?w^g^@CvYlFDz=?+gjO_HK@SAdM zNg@ftaGsJ3GsIDgoQ3p{AgZk?cHU2{$TWvh^u~{w+i?T1NnY;%9e(r?-^(jk;(%*l zlURB*?DQFZX?&7&Qj-q$qlXTqdoNz8@;;%NC3k;aYWo;BmVm6EfGoOCO6&;?>{)VB zTNG`qtXu(U7zqyY>9wOW+0Fsm?iWM~-eA8;mCOa3{oobmh;jQWMY6N?Y-@*e zEQB29xHPE>Qd{V2{||2qwoVu5^LAihpj*RrnFkKURq_o_I7>zNXXg+Rtf7A4jOCvdopR<@#>frG_GyzHX?rRp9uz{d@u%KMwHxc zKxcV)xERgCg+s`v7s{q)s9Wxayeib=rIwn0#7pwziD9uumxZ zf30Q&w2xK?ScWfaYPm(6JPaC!z80C^QkoOPnWC(^20!tP%Ge7^-E#lTlTA(l$3wfp z+m*>RI5JAM37r2ZeSKR>V+TSW#CJ>Xi=hy=K$_4JkOe2O$tSYOvr9>vl9C-uhlxIr zI(wGae(}R-KL0*jj`U-$ebu~hqX6UXJg2KYnQ>3=+xdoij*C;F(S}>ZLzqBFymG9e zv2Wj6S*4)^tgpMzb(qR;`ChgfZ~easJJF}_u{%(I)?zs5j3NB5QX8^D>EXuq$R5VG z$N#Fg5vE!R$pM~BQ?0ct7K+s{HC}P4qf|xfx6y>oNzaF(WIuK_Z(N+8RC22%-b7+) zDK)H!9V7;cNu3ucz8X6P;M#6wDg#0qOGo%yrcZh0hH$BTuTAqkDp*?P)#yDAz=oqV zu6z)5>TxfkprA8Vw9Objm{pxtsy(-0d%2f<{Du-giUUj zZ_<^IxMbtyZVzeQ_i8?U8;gP%4JfU4;zup@QPDfKPZs`&I}fv`r~=!=P)yrduPgi8 zyt|~+JyCB{Y~Yz*rv=S|t+U0__ilDn^DJ90+Z=a)IahIP&!p`BO&R^04-$~yes~x( z5)Z;+^jI60Lr9`x$Tc#*xNkV?8xWk$*zk~ZNVMJf6-CBA$REP`TOhx_HoFHI-*2xx zwKUl)bA-QpNAOg)iA!JZbr@z80Z$dW90T-4_im()lN4ViH&=vGx%`5KIn7V+DOYb7 z>o$^A5J>xVy7tRkvjc$dy0O>UJY`9ad2*&5aygdx_NGxy5$kvpE9kl`Dvprzi5ptk zzzYrURsh**fQgsFnDcjz+Z-?KVBl%K8;G)KorGFWBiY=EZSbu{e!V6h7{_=lk5 z-LIL`M12a(cT4lJ1(1iw;v@-|Y!rDd0Ua}E-Jwu>T6T>_!upiPcd$)W;X)>duL}g!4dKekoSwQS!F0 z9v<`%yKk^Ne?mDsXzrdpOwT(MaICTS**%KQFh&Igy0^kk;4|#!Ty~FJs1Xzm+*{^uFNo6$gLi=Cp=`qWSDf5J#X5-H|{0slB;m z!+fQrv`T}4p1HAk^Ld!jxUY?G7BrDOJ}40K0=1_<+}}tZn=d=N9)qCG8iAiTn6o*C zk(+3l(=LC*bTLyjT7!YD;+rOl2NxmnaB;AV2a)Mf^AZ7Cx5I3NSld~i+GyUv3gf#tM(06ZJ(Cr_HOHtdnx-3^SV^bG&4p{1 z_Nl$>G}q%D6_5i*`^#Veh(97^tDw`9nKD>|QC@@0vOZ+5L~*w8mw4fU(7-O?iE@<{ z<%BcekrDiu>P_Pvfxq&8*^c;6TCLfP5&CjIph zH?W4KBclT0KWZVyM=d#gLTR1FiNkK-=ION346gTX^bu^Yxfcp_7-4k1_@<9hZ=cmz z=#)DX`O_!gpqg}oQvMeQsa^D~h>o@NBrFSyr11T;RL?PdtvAd*K3g4kFT|}DG!O@# zZG9!Ju)mRc2Nl_N>t;9UQK`L&)IhFvv$s2V9XZp04*0kBHPnB z)0Id$#S?+?Iw9mN02Quke=~y91>1A+Z&$=JvY7s91eNvoaN*CyQt6r$ZvIOSumxMa8ozR)Nw+>H@S8PNML|*z*%Qg)zyiYOcsg|}u?walY zDcgP9md9z=RqvyiJuJT3JH-FCd4O)V!ot{TqC!M_Ig6ek%U*R};>dYvlHx@Suml}R+Md$Sk!W%lnFH%Hlgf`Tc|I(HwKv`q5wiQnf*X?7rv$o&u($tol zj$CC;4#V#U|9t)JUC0BTp6(C0S&bN4m~2QW2`|G_m9tRugstk}KzhBD+^OJx0tU&u zA+Ozc)#kz!gxB-Ns;V#6UVWRIgJZ1?pB4NZw!XbPxN^l|T@jDo2=SJL+Yh&FH;AiX zu36{F1~SHBubkCJtRtLId2d^%)MNZ#e(5sIC!bwZ8WMTK?I^Qg*iXj@90)3JlX2X& zY&q>Za4#BvHk!qPwIFMy_`MwpuqhAg?T3Anns6h)YtxAKH!@pStbB@H$xjyx`Y=?Z zv7t4v`rX7i!FL&c3qIQ@v+0dd?I3D_btIq2y7m59!E=ivH}G2#kDbHcFK=iK&DoUc zW%20pPZl>F-FJq~n9>uTwl)Y=zF!(X&kRD`zgl2DQ*jDVYj3!-}@Y7*&I)Xm9j{Rwn zg(_->%XlF!V)?4r>T-RcM}8Ghou&OnoZC+*83WIkcuIi&NWgRJ4Vt5m#<(KEfgSVZ z%a?M%MI}YAHHeh~Bkh(Pb@JQw%flGfb|!HS54t<_Q6afRpPp(g$-0i?>t@92syuss z71`1`3z}S*n(k00M%9cm|HYebinnUYmJMsc{zcGN*++c$h3Wiu(~^q&`CAdEO%MGe z?{6qbHcQvfa)Vgf@0djza5tbrenKO`=Y8!uWauuyR&P{|m$Ggrr1%vx6+)<%Cvb;| zPc(tZxu*?xn22mh!}~(kbOw=4KsTvocB5EjSK5sCGj#n<;Aw!hxafOKmZvv&p5{C9 zP~3fDI%QO~j(y}D{K2>KQv4-gmts^i|_1aoxobLdggSJSFYfX zp`mEdy?lfNzHoJccwnq?R0O%YS;Dz*!}&dZix&&Z<1u>y_c6fQHw5 zr(9c=(5;#NNxqL6CV8~MsUByp^-Y3M_lKZthz9II_6b)OwX6w0iRjUdX{f7HufwxV zk6NP`OcgG|r}?5#vq^!Y02(#ZpBn4_?D&Q6I)eKW6_Nhy&cQ_5rd>V7r12I`!9Y{; zRZAG|_T3PlfHABHt}g$SQV1LR5GnwD_N8miHWm?vov4s`uZoiW^=?fXz|0PU8q&BkaZ+|C`gq&XECBSU%MjkmHs~O!u&? zucC1`(qypIZYfE`n%Sd;&jR(&bln%-#^2~4a~uUNW6zy9_sk*spP73MAA%Q|O6Q!` z;TFzsL?WH|zs)*R04ra{2Y{Tt6dqwcrbgY4Brb<60Ger>3Hem^n(obPYV#wasuTB&z+7>HJglU8iKAJ zq&1X1)Ky1KL)`qLb~syGqGzXHgE3pdpBwDh#$ z#R~QA%kO8>Tx-B*q7AG!$B|d=%M2cjKm6M-80uV0#c^U~`GPq1%|08Q7?GLI4o5vZ*=Y+g66@tZx` zkX>(7oztdOQx~_Ja?Up4V%w7kGaY?vYr1Siv#1EG3w)W4k0rQ#dWW!TcZ&Lly-ST8 zt--I&=E>47gkc9c@UVp2tIw>}(hdiDT3BCL1qrT|>JLEn|4v%pI$ZXma5N{{C^G$_ zAlUOh<@F)dm70YFIYp;0fp6$ipX54?9K|zcLRpB`8?*g`+?lEir-!m{zda@i_H`B5 zC$YZJpH8h+E9n`MT;|1F;jPfs5`Z_BvRRv4HD;_h_fQ-KbfJo4^TU_DhY_T_-qM5* z+=m4-E?zG7oJ>LcU57h)x8*6JdN&jK$E2=Fd_}N7Zi`d^5VoK;p#p;hhO#cCj!v z5BX@m#%N(RhQkG({kHm?i*&_qSGiGG>!C^2*8cE@*6B>MMoPwVzUDTNM90=}rFtO)Yn*xSq#EYd&FV+`!yP@c&Y? z0w7K6oC7%g;=oq`Sm0Sse0v4VOjyv$`)@cR8c{^~cbGY)$lO{vHU$2-754ik&5qL9 zt{n|U2q23^C?j?hQk+Q)niiXq#%LoE_vp3mIukE0T6?hm-Lo8^C1yP-Suw{y%`|fH z&ryr!!YDGsE$LDHz5nd@-rZ=nM^eONKa3Y2W3-8lEX~S}hHPZ5bMZ@3J?*-79d92dtJD>o6>B+brRhmSN8HRQqL z=x2o2-eI4h%r(Fn47$kA+_S>69C0!{y(OhPhUF673lLwE;4oC~7+Ficqm1C# zF;`iwLB*9pQl9f=?|WIoNQ6UL?RaO3S(f@~E8*c)hNpqZaWl1F9>I*!P}`ib^*hPj zi8i2%i0YtHZ1ov#tuy=wc|>^qz|wOhxWwJaFi{wa5?gU- zKL-JL!YO38BQ^@U&uS%WBr^BscH`EI9d{!Y6P*=X&DNcrBUShPe~G|6g`h_J1jBw^ z5KSm-XydJME&Xm}XYkWu?}m2k!1PzYuV<5Hdm#C7UVP^zzmfXmaTEmjIP#-Mg+2V% z^QM3NgGdgqGSwh&J5=`rHQmW8+xD!RHCZ$S#U+PPB4L5%r_G{05bLYY1@Ma`1Ok%R zj~l=n+O5s>x!gd{ORqLJSIcQ3Qsh}ifC=C4wiZ=bm?nn;LXBB^^vBa$77uzNBdw@9 zhg!!>`eEmM#3;%*k5Id{X5(U*er2{noTMe;y^NQjWL9C9B;EX6UAbw}OcWRaM4r%5 zfM1^TVVJlkfmTf7A8IJQHC5L2TfSaM$FBWNyMScIt~xBVJ7qX^Y{5a%uVQAD{+Wo$ ztb1#=XCuPvHjeK~q^$Nr@^k!*jdTaSv!8Wl!2F3&iYRygP9jQXoqA}?{LuH$5E5}L zX&SiQJw9}iJT@`}+p-*q@7>3+hI+_6aJwRnSKtKocU`JGY9Te0Ctk0Jxf36-^iZ*K z3&xF#@M1vd6FNt>-|xvb^yhHUCkmfo!9?J>q4L%*Vok`~9xAR;xB7iP(evy(H@PR+Bi&il3szcqqYTP0KGF#HUb`ZH$zMX0`TzU zZY?6t>V6%~HkNh#jCLiNql5eQHyWtCTqA1e{?BX{e?exbH%csCw0gwWe(-IF(rr$| z1?mSVog>{8HND;2*!S2r9fw-T`|Zq)n&($$oieY2bT`fiB9VxnDQ0Q*e{o=RE8<4k z>9tbpjlcf=o|~9I?Sr~sty>RuaA^~(&XIIFC|dWgBStj^k7@|AvvvLc3|;elc0@HI zGxl$9XC(Izo_weexjjwIIaa|e4FC8$n?m`+BRq4(s;8XhUB*b`%H7tjdXuY3do2Hb zI3{N(K4vJXpXzhwFK2-tb6SqB`8RS{ql8F>!Y`ypWyEKEy}i7>yg&;Su$dJyUjP2k zQe6L$rSZDTQ*K-pn7@d=lrQ#RK$=P{@^1K`y@0RB!(qI2zENQx4+F(uE=9HQjTWgT z+)#6A`0sBHn5{zp^t0si5ckIpeh^tV^L8Fx-!V81NoEKbbT^dCkSb|J5wiOt(heed zYR4WO@%!aRtT*=$L>?8YiffzpT-)OK@>J%}bdrWV&oGdNoVs$+k*o z_t#EmKf_gluXW!)qyD|M-p#_@BPZ5HXb55H?cbSD#B=AaR-C?jHx%yDO!G6>m<muT~h%G}6<18(a5`*+ia z9~gRmx+6qh*-?756eST6;WA}FLqKUuNX<0eJpOAEo>4Vx>2VC@_Ha8?4ZSu{k-455 zpmoNRTgkHbST}O|V#6EH7)Z)x{G&F&|MU0D4~jqEMb^;UN56jJA#d3<;@D_6lvhqh zX?@U(LiU>S;|wIm??Xp1bMj9!nJ+S<;|do?3w@13@|w(K6#K&)0`=2xUgW-+L_iIDHwJ`1VD4r%$& z5D)(rYWUb{Z2WnQymr7b{mkiL=V+nU1?G0ju#szQeVE#QCP8#3Z(}9KeZ|YivaLFi z=2Pu@jQwiO3%0H+NE?F&jzy1ji$mVI{DS{+YW!ZB6$#Nvm%r?S#B%A8VQZs9+uBaM z?&!aauGJYEhbG41-$_vK)C%MgQby1{`b|+={1YDv+uARsO;I`N76$t<)zn@Hs&dvSI#f6S$By}vL{Ruso$v{&MPc5xT?tZUb4_7ag6W&3W^$MS%& z%{Jt8by8Pmx4pQ*(ze9Vzw$d0ktSxivBGrtURA}^)MFIU954vY6g$?-yI41tnfBu5 zZFexg^k89Tj^B39h1?2{djA49c2@p=3!RZ@hV%)XbGxGP#!u4D7#?6z0D0E#7`yi| zSHSi>%@#+iwd9xgn88zo;C>_gh-Y9w?~GQDrE=zK*j19oC~QO8YckBG=Fmovx6KmR z{Ru6aboO^jE#3wz`6gIk#{XNyXOP*JSfYi1 z<_opK(e61jIrHhxQl#x@-}5uL`PqqM)ctDnp%}7G$0;hcXWX~{2r@^UnaE#NXawIS zUA))KWJrW?&0kpjxPVz)5P8(>NYoOcM|I>ELfyXm5q5Na<~2veQ2DlGha)4DwxSD6 zFk3#+I&M4zkt~jd&*NMC3_h@*@4Th^3q{rK1aM$(q%B;uiosI(XV{!$(`mn)YF2 zYJx%&WMKv!@b2+Q|5@{}?}Z6?-eF7Nv4A_{r=CUIye+J)u$vK+);83TuW<;sdBvoi zwBc`YY#8HT_ZYtWQ?s^>o|S5%KH2H|k&Ro0o{VS_vq@?3Nns9-ms| z$E1;Jfwmh?8C2%29r`rc{Zv6vYTXXnL!|IGU;ae=z9|9r%8z~eX#^hvD=HqWa!WT`O9gMXS=^snJOUma>?~08{CR+*% z6aJywwKrcgk5BaJ)%52gp{;FfW#y!e%)YT&j@~}IyEubuOYZJgZ%)#jbjAW9)fqyG zWq!)|O%Bojx~M65R*m_%zi=FRy_@{GwF)cskPxOkdU5lD^~pv3vIj*+>b ziCV;AU=hH=)>E_?6vyM{yb=BCSspZIW{zi!xA=7`WrR4`@=9L_2lG={uH>+A6+G%9 z$L?Zdht4^;w5)#IV^MieEmSqL`r_k*2c*<^My(3_3j*WEgq%+wG|iJzaa=(|;^{^y z9LEdeWHn1tfXD^(zw>1e?+V;_0CMM|>FD#JWFPz(*xe#iT@10!o z(=I*{nWrdZ{6!HNH73M#OErDXJfX)FX@YspXtpuh zbYo)nJc8VjloTY!wX}M0*hiI*47&B^XxU?zk^&s^fQawfqpYl}rl8OjKT9sK;WEw; zX8w*kXHBViUdv@GN^Dv4%lqG0)PCyhT*pXv|K zl5%F}5@@!CZ{FzYf-3x+xC7A~w-@Wx}1;uTOzhoju?nsbNsZ z(|FH&a)5RsZU8+yq}p$r9yA@lu3A(4C&=*gn%qzQH2pOv#qJm2(l{!6n|if>D|>jf z_Q4b=LY6=pbqTm#_DsweJ$>*jocO^(v-V1Zyaw@v9?W-pfh|={#J_<%&x$fP2y(9M znrdphs(iYUQ1O%kA3IM+L{Lx=II-3`n`med0kKYMnxEeFYtFyF_`h|ee$q+$xNPu{ z(4`d{IE=B#LaPi&PBk??Pk+acSWbNP`Zbt^64=-g9RZ*X&;p3ACno8+%K{%u&^(Jk zqAeu>4kzW~M}-$J@PS-lV#+T~L*}iEWT!6XomXy~T>&P))_!O&q`~@2+ORSaK9~ii zy5NBVK_CzKjh^n!tEtA=*g4AK$%AG?U{%;}yrJWIfHzv_Gkw1;8**uSY;Qm(fJJlXoXXp-_whh%SB&4!-|Gtf#9BI;0fx?m? zQ^^YGYmABo&8UCIY1=5<(6=5uTMX}V_4*2zjL0~d25XnKX2)kF34kC4_^trI##zfC z;A;xBYQPFpbOv!T*2$FqT_&p&jx}?ezR*J|3guRG6H-Ak4NF8Z{ED2TdWU$);w2luScH96j5ZL z;3v22KWdc=rudQbP8?64ZcqCK#aTx);J0(g?bzB8z5zoK zUa<|(9dxGpKgVVC-n`)~!RTWiTlzr=T6f!k^h)ZW6AN{@T{-O_BBa40Ba^q(6D0_| zOqW+y_%qeP{{p=Y3=IGTCh!V6I-i2Pz(WPbIIC;dOB>iQynUM$bibt{1^B3HY9<7I zR2WNL-23fJ8P`H@-XGkx%sU3u+rR9l+=-c_rLV#D4Vm#--pRQjj%#HsiaJ;0r1;?V z5g6HmX$U&PZd8H6B#kI`%}&{B48GxPZO-|9yh%bdNmEb>%L3BNDP!G(`ug2xEW z5nEn1zhyn>r3Iecr-<^GLVc!pwMp-Whw>K~m1TJ0IEc`KIX#>(gyn$f$OZO!bOHQN zNFEF`Q3>6P3k~Mu%eZEvz*m->=?IYYKi#pL>I6Q;gX_cDGteL?dJxz_xB|n|3p(Jn zKxdv&i^o#Sb_JSsxq+5(g%a?X0w3NSufKL^19%je-z#VqW^p#?hJF+p8j1P zmIXE7EBE)B-ch5z8yUtwAgcq0K7kW-%v%6`-SCripTDc1GkfiA#Y-VYP$~zX4ur!? zju4x#4EuQP;6~L^8NpSxw47ZJ`qtd4fkF+ayRK6xK8T?ACUbL!zz`n06AgAmQ(FPt z@wyL!w$Pq+YK;dqp`vXG|FbKApwB%T3%j8 zDHAjC@FW2B5}8!L_tc`_`z<{V)*mFSIapY{k(|P!#)U6YzMw1gBT${m%ZCJgBq^JG z0cNT5OO-;~v=)VF=-q~PVzzg3-sfr_5_3+XV0tr1S}6UMl8r=ZUXKWG;~vYBQXA z*pxUT*ySMy+WPf|75v2jnEyJi6P`e(K?efe-+A_@3G!R`rijKWPqI`83HCVI$;nBF z*@p7r3|bHccMwo$-}g_s_LDT|79|R4px80d}it-bn?K?m-vZfw+or5__SENt^{eHk?Rzu<^ZqJRe` z!LQ+AWw__J5S3t!5pTgheFWM^(UXRMAu%1;RRYb@+uV*ZJ5mX; zCxs0ZTd>4?ECmX#O|H#`;PGhJ_?+PFVhvW9_bw`~FFiO{v>*^Ew4w?Fcp)wVBCr|A zF>~u~fy4t6bjDTzIw~O!7U-U3XJ+m)+4=U&_uAvNM>^eO_oDKgS+Y@54#DpT@%lTibs;YK+86lUKz=&|B2YB0Fr${hn;oy){CkILD zQYjXSJwdz50>z0P|FpXW$WKhy*=fhR_EUUV|M${l{og;^5t%Cov~RSsGB&yh3Ph0q zHT0uRXn}rr75xPsh_w)0hDE?XTm(*)J@4y*OOP{9%fClpJg;kYI3?KG-kGgfRGm5v zyj#mJMX!zF;K^-poSmxBXXQzc1>@%R@B_sMiwNLWUcY&=@)dX&0CC&ZzAgxRSkF~e zftC{^>a*i$GjP54+M81i4D3Khk4IqQV}X|WyXYQ3pq7pw-&_Te4Xrz{1m@Km4Zu4# z95}Gw%N|xa+na9&l(kF|$1D(=mHWR|U>v3}*<-yaYyN*cU3WOu|NBk}2~jqgg=~_L zl|4d6Av2@M-kW4*%ieosM)ux&Z$kFoBirvepYL`3{^l#2&&1y;L5w-x|dT>KtoIG0|$|wp8l5a_wVyt zCH*T1bp$xCKc%DVU@7nO$TVb>vru7{=sC3F$Q9yT~M)r4n{hcCGvDUoUuq4v44JJozc_d z`pJ4@e+vi*@L)&K$csPslR^3YLBmDFB9BTC4LflzAs%W~`L}NecI@Cp!e<;Ak%S~d zRaKy1b;{mQic|de@88vJC-r&3ClhMO!y}92I!zh++>C_stT1PF@~w zwqZ89s1N5g2Ms;FUx7yT5L6Rz9)fc4ih?qn-oYDa4b_Ezg(z{9l~OnqrxKpmA5**4gr-6w%HQRA_d~B5 zNjF#^>U~JAt4oAm0ng>Y<6t2tHy81ij3Rt!SAdsn=YCLjI2U(c@)lL9x2Dy!^3z`YKE< z1M1kNB|WIduzu>=DSgE9_O$xgB%s|9fj!%=1?P;z1N$4q30Lqc^Efzfv@({IyqBw7 zcv5|R!ABNvT>_!*4-4U`eE z_3b7Bv{umxuk6nK>2MQJ=TXYBnEQDhpaA0%LZgHjNB#8aa&lJQ-{r{xsP?t!bB?MGCzNR z)k@2!P*+etd-ffIkz1%P{_@%Iy?g6-YiVCEJj3R0_^#xrlz>VI=SWv-QXaoJ@4g16yakwcUimlDu#BfPr<#9%9 zG+%5cs0bT>8xDg8-oUn@kK}oFa&oej?Z6+bY$>11xnyGkuAi4s*?Ez}6<1Jbf*O7! zX^EPVF<>xFJS;lew{yf<&?T<74V3?{R?9^#R1jqp95DD$ti*vB#2mhmT-#$Mkk4NJm3k`=X*~%$dNxO|Kz`Nf2!$~zCNUJ)Zh8;?(T@) z@VT(l2chtQ`a34Kpvn&xL%%e-MbE-Q`{G3eXzFjCwEptTCClS}8;MfJj@?Tl5%k^@ z8x|rLIsG%ed2$QVTCXU9M$6^dA;Ja?kyxNArr_tF)xy;nhqm|myLSl#3}SXdqU?G% zw_D>HHP;cr(T0++nTbh&Qe7~&m9iByLr>o2`VhQ(Zq7!qJrsFfXK*9XyXg}5(dIjKeh6*ZELMnJb|6yWUQZ93sIZ{EIc*m>T?&(f z)bRjXU%0H$>lLUHU;13WDN( zJ$5$pzk-MYh`IjKvggxtgtYhHtP=}FY`x3VeMXg{uZF{!XS|lR|Ld0vP@dqnBf|-U zc>TXN8;TWfNQ9NJVsCA>_7x)XR={sKg9dbcFbx%#mMB4KmoZJTPIS8c zy2Ffn;J<=}e)bsyw0is7&tu+zg|PqG;i?t*xeV1h*#T7owvD8vrrypS^_$*Cm~_m; zJZ<$i?3ra{9k6qTI|y|>SAy&uzgKF>LPF$xU{r=7F;wyh0W;#cf#=_XZS}7Y`EWtr zDk?&_SZR4Aw}kMMxjgKr<>mbol>{w5d|rqmVxOTGXeWjW6Cf7wLv*>gedi8BmsCFM z$bisVeJY8GjU|SK^oX}rP!I%X%VNf6++sOGc_>ScyysKN?&;Ff(w!P-Jij(X!49n> zTnElM*Nc1>m3NR6!K~shj7#S;qBIu}f9h|p%)++e|8xdZMG>yP(4@=>sKH%Hebu-g z5(nFUp5gj62M@2OS5#DdU+w>jbO+(+5-s57`lt|k(bVzb5p7mh5S1g$R6ruW62?GB zrwt(m?3yo5p5Z1$z=wg1v2!5+feO)oB6Msi!anCGBp3=Ardv6DIxS!?e`RZ0d;I?P z?c2^5ii(P(P>?MU(n1jm#sYvN>`uSb@4`iVlbLT~hFWsl4LCT@K#KUy;c&%b&gCDQ z!2reFHPql?NlC#ovuE!$lt{^V!uPZC@&W)~`UW>M`#Y+fgaq3B{Jf{PcUi?rgy}s||`n7AXG!a(RS{+X7Xt!$Bu4B5LCMxBg?=pEFo;w?{0_fI$9SyW!Se z3-$+4b8>Uj$9twaz^4-*=1AKA1Uw#uuPMS5f%ocY7{U5sGVSNTVC?5T!~l$w09vw> z{<7g@^-~_A*fDEGXi}QrsHzeH1keiK9sFj584kk5LPtjj(%tm}KQTizzQ5ON50tY= z_uznlC(oaU)zsAV5&V7+VBpINjlZ|839PB!AsQd3-yv5ae23s6(?R)@31L(w^8-iQ z#g08&T<_85c(74HrMJ^grH`*K+R65`I6#XL5edekjPocnhvm=TA$u*bF4|Uic6U?5 z3&9T|pdiG1Lct^!dHXhG_tnscxqbtVhmRlkB6caK1@m{ocfm%8fq-$~e7RY8-G(fH z+Rq++3-$v1*UA)B?<};Dm?2u0Wr#y&w$)(uQ-NryUYe2dIh`E!A7p3!%20;fgX?$l z7hb-rO9wn@5J?cAzd{Xfa<>|kWWdc1%j>TK!My*kIafzlS88Tvo0hso-bW1YB2DZB zG1pdEV3?s#r1SV1rD@0}D`q{8@whQ_k@#W<_oj=GSL5ZIz@ zBDre@n{u$;71lJGd`l)}K?AupH=Ka#M^%BSx~tKyg_{5_--V+M?N>Jh z-u2Pkamvjsh&<5QlxR2J2BV#AwdJhgJOv{RPtB~xwK$6X+o~*pnw3E4fehOXy6{SP z)QAxn8*w12{^`CB2l015KOIN`f6SXP5zak;`1Gx+W^*wjqLi^QBV-83qdOnL$L6}T zoj1So8fr^ueWEuuzA)n7NBFlOJY7IhT9^F$daWvGcwGc-emb^Uo7KGQSr#jFhKREQ zcX$}_$HwXuLvHovT`DG^2b779Nqf~@$pyslMvSNSE9~nl+O1f)~aR~F7B{? z>LEmUQ?&&>_V3i&;LTl$`Ifwz!@5Bu!9^|v7cIpAOJf}TpmhE=xc9(v@n?#ezj^4= zHXXw_tGEI6fT*JGyx)ol`y!&EIZx;CdWUN8PYH1F@LCam7=UXvj^rp_h2i-@bC>MI znZhFnNqFxeDQOSlc|igAJt*iz&h*=)F%6)#EiEEV&CP1>L}X|lA(k~j%1HyfpNT#V zy>&2PdixVtJH^9|0u%Nzfy#SOd_sMCe!M9T#lciOKEmXKxpL{@9Q%I4?Ouv3r@L)m zpORwaM*;MmN9CCewLLGvM} zq7VCVo;-PiO-j1O;&@9FZ{e-_C|&|GrJlJKC9gk$b&x1V4Vej(L-&i!%uH{J4K366-y@{)Ck^1WxTEBBFH`OQppZ*IU!&!JiWy+5}WBek*99zi>!? zMVW`{d9)YxzlL3R(e_&|zV)l<<_2uEpSqS?l*q|WyfK(VkFp;ps38jGONYzlwj{)p zI9nls``LH2K717>DJ_i%Ioi-?@h{3*MJhyfivm{W26a#jmqn`EfqWp$3>7@Fsy=fo z^`d!4W~m|o0--_5XVZ5(cZFRa2GpUP_uaP2t`AQwSu7G-Q4YC!ZcTLis?Q_a{aOKj zH0AXEnY7O**bMJ*x!+z1DcByq1SsH62vXqDL62GjkB-xLqG7lGrq3Ib1Zx61Z@|YM zK|gq{YKF_)FmGE8=!YKY7XaGQDW8tRg!KD6s=ntm zEHtzc(8++_Pr4IDAOAr82u|N!T6vv6O2AHM8S=)$E%H z%0p=IbqcAG01qM{h6zL5vU*(;0Tl-v2YA3!>b31 zP|J?47yYxEIOY#G{U-;2kDovwRELOz$$xdE2p}ai*4I?%q70vN-@xw%)OxPBf@r_f09E!7`N=$20-{2m3aPpM(9 zxC(V^vRo4XwazVQX$}BM4F#JsF#B^(T#T4+ewp#u*Ir;WI z3Nry0X^|D}(vpOO|J|QIwRbh?=t9n-qJ*W;m6SR*(J{iu7p=RT%gmzOgoPW8O{<1{ zx}WyP{h%sr)CrwV(JAlVl#{uei}&r@(yq%za@B{GZQE{u%ze%s(e492qIWj2jk!Eb zTw02`nv);lkeMa?K9X^F;H|1=XJ>JMjYgu;jjMxclyu@h^Y-l_n~wgd{;Z<6*Sqz< zi-zxsiWbY_;&*Mz$V9Xg1s3;QDlg|`hd8V(bvWhchPC73bKa;Y$kIIvU5<^LE^^N! z7pu2_fA?-kBuzIi#gnB5qrC;Di{BioNoi?unfEC-^(>2dvwCT|r=0o@EYIGsv+Fs! z*3}f&@H!Z4C+9KDE>cr3mrPeH^9y7>`sUH7_r-A*fsDeOm;Z84IPX2fXIMrIyA-LP zJZaHVN(o_5dY}2)#@l#~^xcq0sgA%AH<3=eMCDJ6>{(jfMr2Q$FLEL9Z35+`1p>@} z6c;WFD(H5ubrMjsFn^wOI%*dp_hQ6;tw~BYClERDuDE9zP^rT*(cMUhp2;~ieUiM1 zId{DAWs=6h28Ts9%Xn2$F{N0c?C_fX8jbRBwpNX7(^`tJ$=Fl^b&`K;nsyo^DmMkQ zWP176Y?%P9_y*6X5@TU{dc9XMVbsAkJo-*yVbq5uAC}MlZG~>9PclxgZG3SkjEPZ@ zyN~I$pp)Yh9@%xm%%wu6raW#QluG50IHcge4) zK#h}7<$UuH5hD>RiY}o=M7RbbZf->d7a%}VQc@~)&fGA>LrBR1bU62atFGPy$bVq- zy!M|8ivV(gS*@jPK$3nyHLZWi{lS7K*`C(enRKZBmX$KvHd1Nj+tXKT#{>FwN>M^X zcS7QjlP=MdoKY*!2KcZR7Kn-Ysvp14`#P@Wvs)j9-dXff;g>}<$ol6K?DVoKr3u!? zWH_(MUD3cK!1ySN*-O8Bsmw4;Ny^Xfuhz{E2YiC=puM)hy@T$ej~+_Q@SBviB8soh zJAqdnVm6I(^Uv)ddbT#jXdI>4uls2rs#!f#VRa_T`H9!WgtsM+%FmyFX`3`)I=c1Ux+A+x(138wL_qrGq`c;dO z-xq}QNI2IRZ`y8uyV&(HtM^KB{b8c&QXB4aVZzP&UIVkA5A~I;U=qKA&vb=++2Nck z!*rhe9et-#gHYp}o=*}IA-k#;zl{Cs&!sZ+)t?=Q!WXW&Vl6eZyd3p$^L-AjL-5t} zXN=+cr(h2}E(ob(%EPfcn<6)b%cScY$AU{GOu0`Y)bPIYN zNKPeyKf;*7w46j{m?wjBjir(#ZP0g8SJIA)6vja3{-3gPnR|QVpm3Ll>V0sd5_hJ_(6g@-7PsAO)F$~pkjQl$!2w4S*qw6je=*2%LG4&r0iWd*xz;^`3!S3K} za|RP0G6C!CKjZknejV^${Oe|ylnz^ErVClV05;X%GV2^?lq%-sx-&6)PY+abBCl5C z$mNNwU}sI1dd9%VSmGX(t5I6%2fxjdIEs6sk#yTSoc4&hYTR8knAdR*w@UU7`-iyf zF87C`E&o_3gw|Zy#zYNICexWsJ?XDCBkzfQiue1{Sv9HdREQYA);noLN%ejGB^9Qe z2Wt7w6&?X?Qdp)WsFiO_I^@s{11{!K9coElN{PB%*s5xIg>^cD`gltT5BT97{Aqt|ms}nSt&A>KL|gp1tA)f6FDSBH_|Fb|#weqG;+HUqluKkX{pQ56j`9GzIxfB$Z&Zwv}Q2XM{q?WCl zknQ(gSGx5O8pcnS$wW#9lNI-cM>HnbNoAN%aj2Gy1%>}3!=c3cW3##B;DZ%$taton zJt8cgF=EJRCseF$%Qz<>D!k)>jPnisfZfc~q|x7o0GJF_jC`0R3{;Q|`0l>ky?jO^ z^W$R6)^O?K@0~MJ@u2Y4yn@K_fzn*(>p3@;?$*B+PW8L|IqEx;)WbF!#Q8f_(vH|y zdmc({W_7Ck!fC&@c`U=zc-?SOWmE)ME0UcFa7|EZ@l#O*=Y}b7lqU|KinKLbj0;_^ zYoFfdiFTuEk6(QtSkWG>7)oP9VuXc^Ec~umy`?Gkek3a)FG1TZ3X;3euG&;?l(C1E z)pGft^>7nM!LP=TW*x+_ZucbIkJ6fD6&xc&!Te$J<;zcwB@vO1GFTyu_4Vr9$m6Af zW+EvUY@$wd>c*dZKRVKJvQn|AH?5@ltk4Kjdi*GTp#bIUvS5m7<0znU-zun<#Zq@Z zO&O^@KMK7gLoL=(#>n_)1MIYmH45@0Omq#h8T}zQd5|{yH21_DR9rg4?r(>{C{6Uk z)o%Biv8lh*UX2umbMql>Md7#(+YsS4A=R&!{;7Vc?vy`5bd?oCPHoTUhG?ah&#aJv zFb6Ad!eCx|0M^Eb(`NO*ESIIm^3`A=9veU}dwb@9TL+yN-$(}tpm?|g{P8|dzfxvH@SX znOB!yb~ln;ASqKI@jXQ?&kla=d&cI2Gga^CxhN@%b*^+1DTs=Y>T_BpF7W8AF6I zv=RJQC?l&qfM|$rFmX2!!g+75vFV?Qx11D#hV%17d*NZ!N%GFTf9%sZERpY? zt=gZm5-wfYZPH*elW=m_Tf3a5C@Ij&de+})5E}=G8TM?9N)s=P6-K#vcw`ixuEZ&v zfh{}NJaUfnmssZlmy0gPBG=lQtf>#&I_0+cN13ppZrKK)2cEDvhD||>U*T8tsOu>& zEYmsI6YO=S_0tTC6CDE0x+6Oeg1(2&3Tru27VeL!Wl&mm;~NhAnoOTxpElDcI@)3C zul^|9Z2i*gPbR;M-s8^xs9C|8>gL*iUH=+hJuL`ZG`o~%bHg=BC8ErGS+SYo{L%o4 z^)KmJov{;UKtl&f${e@lURqeacwnFaxL|Hxwvq34X;RKmQs8#kz7}OyaTGK+%kI7N z@vG;3?bCPIgNg;WZ7cRB42x?2bH)~hv(#kj(=j8 zG+UAr@EVr$i8K9yl2WEaM)s`cTrHSPq4505PL+%{o`TO|&Q)!86HmEeYx00USRRgS zf8(px?s@Y)xDE0^rPI)O3@il1R2D%HL1Y@ZE(#hNGAb&BDJd!9%T}jwIc~xKLRKm< z9Hs%9Fk%=C{|H870Db2a7PiluzXd`AtobV~Ek)4R5)x!AE!n}>S!aKtwFXuZ!oU#+ zCMAf0Ah>?AlJ$sR(9^e zSo)sz5T6E<9|Fa2l~vZ3{rIS4#;28@%-NTL_pl<{0{zepttaa=Oz++`pEgm5>F5c6 zyu75B!*Lf?8vcy}J-uw!5eGeEIC{_DM$4sVY`+}58goOex9LfaXR4RE%Z)=HpZwBY zcK)YXOMIJ-42u!pK>!KKQYk+@z47m|xx`&YK8>J(-k3Yen$P&evD!u_4Z^))pr+@L zncl(P*=-d4Ls^+NDRhz%Qy?qpx?P93l%PbF-d$dPh~h^mWuYz<$wz&0s96bPBluStb=Pp49XZyPGp44omVi}9F#;A6P zlysQG&IIN3&pq??oUlaZ>7qaqc#-R~A+(2&hL3n;Xez-b%pqTkCHUf4nY+TW@$SKC z`V&>oxR>2!XV5aAN5ABp)$+>hH1Z7yubSEzJ(zX7R>jinOSn&g-WF&VabTs|FGu(>^@&KxB-yKVcjx>CKs z)=V~QxV%!G8LBk~<%0C0B5aLu3!VfZ6~RgGfJqt5jMLMf*lta{1{N-79A;Q2jLHV& zb=?M98hLx0dtf$M0W%dq9xBbWlK5Q^6FQim!0=OtJ;oWY0AO55RES>&GcA~L-zu5x z(T~z-ySbXX35||k>&`7Ha=USjWYLoS@`WDeiFt0t3kGTQq)O&usn~yT}L&n zxe2PIeYLayMqH>ushe_)+qSFHMlrRCq8tflOC4wW@y;^YyJkEz{f zm?M&2?S#zI$YxnLmrbUpP_e$H7Y^p5;nV)$|9luRb&$BbTq1KxLYuR<2hJuCdw+_wj;Fg+Z4k z93ttIqPKmM^&9jRr#gnl#%~OBjug~>OVBrb{HU;AyUENw1(XrBiR~d3hi1%!n2jpl zU6XRjr|M|jU1wV*mLseJT2;TbT)EY)0`u+wxFK$w`%~}iXya+1=PTQeBWY&$`YRb? zt{aKR++71U2>`+NaV=s%NOS;xBVxSFcAXF?;C_<*WIJ?!;T`LAcVkzXa@{X=wRTde zI+ITFo-@NNjITSzX?Qnc+6N5Cc-OaY=hc4FA$NTxf?+alU3e0oo6o*4cX(6(bV||n z+VJ}3&eFYrUcM|09n{(3lD;!XNe zq?&mUIm>i{?GgjAjnKY55?Yl1U+5ZEJRYnS>!LFIdtKsCXFIP^#%olkl3+xa8B?w- z2Cr9s!fet@<$dD$r+T;Tg()d8?!RB0>Q)kIV#S>p{?ipwm+wVA~?Nm`tSMS z%5;3AQ2F$|H{w2Fr^x zF0uwoQZTtM3r4kB->Ry1A`nMYQ`6jWJbZjH7!d#&x7%3LKq^!>JLl8X7@xa8H)dd#c-sAi__f#_i%D|r&P%l3Soc@O z0;N~Xv;wBgRnwE^k7K<>(Nb-9X;xjk(H46gr;qJ3n=ra|VLb%`ossr#BrwWZtk_HsfOx1}Yya^`ntKxw#{6&hvw zL7Jh+$OtPeMH!4_5x3&eS5K_snW(_QdHC%r;JHNoBG#9w!J><#xqnHWEe(F?spL)= zhnX@%f-8?Bq-0T&FOzPb=ng&E_{Jw)$k`+?7iqA0UUS+dHH>CiO;l>b@cQQZrMDk? zOSU|&x-*;VV!N(Q^&s&w14p2shkIQd=iugHKbz^~w`yqudNP-AjjQvIFA1@){roa< zwAtA{ElV{4psA}#7StKmS$(QAgo@Fyj!C-RDKMd(Ja@8vbbcWHhUDc!kZPJ4ioP0HLKWe-AW5>l4hoKIWd+BZF`nB*;7$!w0>}X5p8Mp zE0|QWZ)-9=-ku7jJdW+@)0RJO0+|asf2hThX=zgNaJ8q3kob7xV>b&w58tEJei0Yk zVfUk6u&mOh7kHn2z<`r`_YM!XTd<{x1yXTv{GMEln;s|KjyN0|g5^q@+`PrxE=gRz zb`R>yewS*;LqV7*D`J9mzM_+mV>YCei}7svBH_UzW*37N`zko&1OHSe?B zvqja?tWxJ%hpO6QlsnOHjwGdZ|7?t6ypeAH#93n4fz6rd#$QVF$duBTiBWHBxBj_l zH@-Cm1N#qs0ZqM~okj+JGZ7+dZ+$AmOQK9X&iDrNA-pT(T%S6NLm)2%5ssl87sLS(-w z?bId6ao}m_j5({T*2JJ&GMSx|gG513kFNOK<8xsB&35N_oSjke^}jjI$+0Daw7kpx zC=;gLm9B;xvP;AGB%tmV6!a0SRhC*I!GfWITxAv@WV9gII|;3l(OWA-!-oHmcmF4is&;cs`A?~A|v z7ch`aM2M%^9`0M}`DzG+1nhQcd!zl3$0v z>t5Q|J?4)O-!tEMvE6{8E}<(}G2EGcb;14d!uAa}_Zkus;lje{+(^!o@Nib`VHy~! zVK+8LQ&SB%9XIV=>~b;T_f!QTnHv%}ZMTsTOLa-o?yDy8SCG6Q4(yg5=5PoWVx9V; z-_;gq-QD@xuF~Ao`-kwymwt>ji*xByQo6^fuFFjMfaoA zow-{rkx^DwABZ)&yV)w$Cz_ilM@ne5T3h9D@Vy7EW_Yl*R6ExvCj(gM8kDqx)zJJN zJ+!Y=X7S_GAGTC3O(C&>m<0cIWm0Db|UIcnKq=JJaBxff#jYYwMrL z1Q{V!A|%t7Hm_gj1qUv7S8%-%qb|qAA7Xh$k~SGItsR11RG@3gGdN?~Rb{i7LB`kN z_q|b{osxHs!;H`1-tNEP_lfUxr_>%0yXk%iJH8=%5oxkbB1)~AH7p5@2@J&%kj|Yu zadZmT4rP&ryS0Ws@&C-!iaCCmNn_P-n&js1C)V?EC}YLaz892U-CJSeP1qdMGgK15 zKoKrUc*+u3I2mUq&B3e{);Xe+$YG_$L=<7&V8EhlKcLws(?iYAdCN1wWKyiIljL>E z8fqp5aZG>XZrv2UT_NMP{f<}~PTJr^LDf=iVVNF@z;v0O1~Wly=0_;pH8+_$U)Hfq zRlMANVXuEV?Gw4zWwgE4;;65dYoYw>qgM09WTJ*yj@05Gy-bzRGmA=Xow=$|vv=)ylgOBp`kTiy2#jY-o#6@Q zxA%?4SZ&9d&p$9u1d0u9Zm96-qU=b9v2M?F@8K~JYHG+RnX_h;mKt*qW*4h))#e3c zcSu<0M)1wg8z=C_#CC?onXq_adp@BZc3W3yBOuV-84GNx;dLe9i~UjOvR-1xk~TEt zi}|4agIaeX4i4_Dxko(nn3RKk_s@*ae_P{bXN!;fl?{k-@gMFV_6I<@cz(PSSY7<% zG0L470(HTw@3cAwz0EViM#mX!CZqXn_H$k*9S5jg$)%UIt1k&L>Rbzk*Nsm|{@%zK z26}`VnGx9+mhh2E5&iaQGtpkuo#Lt6N3fqQfcsGbBCt>0+(2&4&ZYwr1QEY87a&2p zKvqXkpUBB~S!OBMvD?U5gwTNG8%yX0lXoDCCnr;Yc*(gf3CMg8FFv;nuiXM}$N^+w zL?8W<#G`TA{K8H}E%n!YMydC#z}%STm%==^#hHhwy7~WBrl(5Jhu_y})>xuaPm{$) z;PH$oo4cMm%D(Co73^XC<1N7$WbAiva=0I13(8XDFtZDBD%BO`Pw^L?M61j*pypVi7( zethaCauDhFTxkNvYTI87FQZAGFkGfg%$|p}*9!efNOQ#xfIom@5}bguua0(OXr7? zOuD9$_jDFS)bHJk4*M~Y6xwKZ2lWS!32~bkoAc36`tt&n#`FM;hbWj$Y_^P%0T@T1 z*b95@(ZPLXtWdfAnKI4TB>sjMwkNyO?LX(6pK+~{_9T9O z`sE8@k@Yj$C$x{|%302RaR^Vgn=5S?k{%EW&U3F)EG$epm)(kO{F+h_l$`!V_I!l| z#dZ_rbK-YO2LG;nwMiVJ&LE=a*fZ6`o6Au#^6x|^PWK1$NpyBKE!m$>;^KLc5qCK} zdq()Vr(q_aq;1Z%IXW?D7GA2tg79Z9*x9R07Lu!4tuCxh4rdvcA7c;?C-O+A@lI3w zCn*f$n-CGyoN9R}S00t&1hSE)>rBobmX{mHg&A#=ULP=ZIu@p6Si99YI4+aS2|FDa zm~oT6soQ1}H{+k>o$iX_-1uQW&9=fw$-DTN3bQoCoM%F6;JIshUN`QRN1@l6{Hfd{ zQ!}E00As81h0#|Zo-3UVPr1BkyCM2*8_=k|=Q4dYNA5PIcKw(v3-e&CF(`_{?&(qV zJt4tP`Zab^?>V%qKF_arY{mrwUdr|`6m}Tyr-^kw>zi|?zC#XsbJ~1uhpvV&Zw~1<2gG8I~guqFq?2Rlc`K$M`)mp{i zYT;I|9`HA}a0z#jQdp0-ib+i`AhAP&D9db|nds$i&+}!fdv{PdU81Zt+Fm;E@>RCB zcIF)K)-h!5w6}F}1-rYmbCue?%a79hu{`9Bv?d4O6Pv51vwRyv)s-z#n!)x5qD_Ow z+eD~2j~?tX@Zcj9v>*yc@Qv$SfLsfDC*WQa^4b^muRVu-j1T6r%oPYKc`5qcPeTfs%C=_u(tq(FJuzifH$ms2_RY7&cxeZ>GvrO4WJ z**%w0dWZN%3eQ2hBZflUog_E>opCJ!DxMF%A=s2FJva-Qah_|V(r7Y|befPU$(-K! z0!Xx;CVoFOmcD;JxZu9L7R$2JKYX1w49dD2O^6|sWUeZ{MXZ%9&r8l0AKZKHp`95%csHtTHOev|Th z`=^USZ(^vrpEJw)<`&|ES^_WC$*!{YI1>B)rhEUYdFrMRE-u=B7lV+e&wU(H#unNn zGC9Yim8d(kmR_JgLAQkL(XZ^EmLyy0gpLm5Kx)<3Xks_t3z z^9kvFzP~?ZBm}_P=YvL6q421SIhn%>b)l^Z$`&sI^pEWKEOu)T&VXvs!o7;dW)Sk{ z-+C7NO(Ng&Au}EJ!gt!G8^nWX>2SQBHWP9_j=+9I+jY+B3XI{P}5 zettn)C5m)GpZ+jYc9o(3I$)=q4ZvH7+VnJ6@yhOexG-!X#2AFDGaV*0g4cB-H+VKM zBeZ1tNT>5yh?ktArJCB4uK!Ug#XH}AIf4alf9$0GhdOigs#;jIrsX4|1-0*e8w4k3 zpFF)MuSngpp4DAuyrWJzmTzQmG35`RdSg`RO#=P&=w=}S70ks+Qc{W^SDk;Ml;u~= zOtOu|&V7F`G>)NVzrEXb$RONpwFgIMQLGY^h-y0CVo^_x!@g+bQE{=S{Sr?zw6e`X}Tmq3A~D(`Wx3AkFUn$Np6~D1$79%-Gs?s*AKm@t~)cHp)=SMhV zY4XS@D1IFs5^UCMZ!LypUVaDN@pqum&GScJ3!^TvsBKFl*re1H6ki7h25e1$)Ph=i z7ZtU<5=i(sgoG12y*j^sSAtgMx6$*sUZAiiG1xe97GBt?#n~whzCs@#vt{K7Emsyk zAy#ohRf=ltcr}0$??IVK9{Uud_XB2S{p&%oH7k1D@n(OPPVN3jv`sU1x`e!h`EulR z^g@RPE4X-1Jzl&#L8qhtA^p8T{h4@EL8Vy1EH3ufGRq8%5fVZ)0!`mfrKzp3yMagVL%sqe7bnxw)<03qpAas=42~0 zSN<+)w|u3qZ)Q6=`rW&h>qD*n`g^l>a}o8*`CrhEqGO)z!WOBep`19*&ZjyTN1+yh zQ^Jlsyq86!$8j&DXK#KywY-bD>&KjvV?OI(IbzyRrWz?PU%Z9YWVD!EG_C*GVdEM1 zx07JEj^tK&)e0GXeM>zK0WIdwWF!pAj~z}wa{uJ|jHT{YONtZ_6ggh`CE|^=L`d2X z*{r-rc5_VPJX(&D%-cII1!TcQKf;zbNULo-le_9~us+8ljcYO}OKwgowrNk*s1`_Y z+86cIGaVf0{oQ*x8V-vMsx6m5f_~n1EJXfPo<+pq0u) zhTJY5!+1x7DX;2->GMw=o9d-fUH>*|Zy zrJF#08rS#`Wl*@md<*%lz^Ck7&m}v`WhS%@th603*GH63h%9_V;{K-%0mZ2VMiGY51Idn*0XI_5FvxL9qtAbdLCAcvUyvripiR zK4n*c4QHIMZOe@13hF+TaW^(9?6i|-L=AY}JD7N;nBIr?7!zm()4!QDmTupN*nJ#r zRv%$P+i{m5^G#^`tAyg-p5BO6a<R6lypC>Ky34e?zKm>?tABkX5-mV((uTn;)2S^; z!Y)Ej>vSS9$O4I#}t+v48Vb_fu1@6gOHDVdZf8eh-GAb$DkzX#Av6UH@{wo7W zl)?(uzd5drXBjac?e7=M3n+d0Lj0$g2|wrb2`z#z8IitIAxw{9XgG9y?)gZ2g4M8N zqP(#2*2vgT*~-ej`kPL?^wiWP2O#AYbSh?cIBicV5|hKMcv5R+EGt8ZZ&`^c>-rk6|j-T|EsqYnj)FkE>Mz%@Ae6qgT zTe9`}Q&Hv=5)y_@dFg^xFE94`D6Pm<=aWhA(J7ekI~RY*NZqvUwz^s*KbRyAwB-H= zPSYqc5VS+8(wKf%GM2k`f?n_Nrcc~ zpa%XuH6Swpw42nUQ4prf`2HC)!v?tR*jNl(jzkLeUVGQ9I!H(H-&mu~c{E0db$_Y) z|5~M8%w}$UOz6S(B{?8%=pL}6pk69+KC{GU(F_><@|2zZ2S~>IF0Lo?idJ1Du2sCP ziS2O$RBx?NOpbZAy5|#j4J8{tsQ8lHF2bFE6}({mO*v`M#!xUuKx^YsEKMM474!62 zt`n})zQ{ddka>^!zWdImOzHWw!S z7M1)yGQ)qP&R{> zaB46O*o6Re2Rh7RgFz}#7D+fe*RarNh1Gd4Jos2LiO?W}3Kk)M`eWJ6)RV{;_p#Y~ z5K5OMerrMUT{;mWH&DvUXb>%#dH%p0&}R*(CaJ=@TSDub(1`Im6npB?7llTR`%bwO z6!@Q#SYBlZ)*1^5FTa2Hlw=^8Vm$>L_sOeD;L7o7uJnZZNxs3@nq!H;(unGM1N^EI zV~!@E<=q?9i>Iuz2|G)$@K)m8*wC9=8?=1#v*<0^D}xXvQVto-_MwboMYld~^GSUL z3${4@s!7H|YR>n>=ol#6ZZp*iV=2J1f3?jS7r)ByVClhbzZ`Kei;5LqTv<6J0Sa4% zlOm+8FM#e!Mt;~;hWhIHHW)nqa zfp1X3P)+HU-S!VMw(>dVeqGiM_7^efzo_XXnD1z*KErF~b4sV;GD~z$oV7JnfNtyLCatFF~&Abl0F^>XQ=sSvf!-M?rd>hmYmt!4osA@YPFX! zcKR8Y3DFalZ@l5R(FS!t#gto(r`E6ir=t37)7Ih_2L~1r72a_XD|b_+MH+_)BxD{f zW!kwoPvl83b-%5O%0RD|K=|DICL( z>PL0DlbUg(B6ryStkW!e5lb^2>w1M?JNMF{in*S?Ua@G!W3|dO7JnhLfj>HGxko`E zxba-KU8-3@JDH;om}SP|K&o?tXiO_AA`zVOn`5#~Ohlr8WXqVF>k%2VO_@bO{&19cptnn(uV(%w3VYjDy8U?6*RzL?! z&b^9(&_8YKFG{@I$EO_WPT()9nT7rK?N_z&&R|%KBPY;XKIhzAoFT1RW=y22stTII z!Tc9|J9L?#XNJu#Kr?qsMEIdit%6z>K&N|Ox|4GZUjEMr1Ei zVDU@%uJe03o~O19O8T#5s;s{UlRU_5$P{}p!<6NXca36-x`Y&@EkSte>YQwiNM7|> zt-7Coh7PNBHU}CWu7vK_X(mjhTE|XZX=#1kO-sqi1Ix<(^>9Cle>)09P&%H8i%IDV zjb{kymP+Z(jXdAp)Hl9L5PyJO46dCy8FE^-{*P$f}i7X9O!o6 z*L_~+^*KND^}aYD!I7V>vGM8J(IW}f70;C%FHXGE^1aE$@UAhpN~iVB>qL696Dh~S zE1YijuyXIZnekL%BW;_DQtptL;>m9ZE~Y3uySj$FePkwiQuY3QlZ2_S7C+ls6V839 zsZuW5zQ>f8T~>C!wpwB|M~7eG@%^d!+dn=Fwe8xe8B5jtQg=*BB=5d*23fSFWJZCQ z&%oH&{K@S>)pbuEUO0E|UC!$ZVho#;I+kHP-it-%P*${XCayT2Rf32!HL0h?b9! zwHzAX8gs` z`;IB~u51lrMoz>N|T3t9oEh+2c_Gt-8QVM*BULa@W4rJXYaWHRo}v?nLv@w5>_v z#i`KvO930C$J63_lBoI0PJXxb@OXC(&4GQr(F~t0ho03I6i~0XUcA-QQ$85Ay?TQ@)f`10N z=L;M-a0lA-N%{z+?V z{$x9Zgl@}T*c!ztDNF6ESq2^cIw6X^2g8! zCq5Ztg5**>2B7|hDKCB5pI(tO4;66$Oyf$>D>_hng_36e{ZS+7%irE~IO4NOn($>{ z=Egwo$b0bMb>)iz@M{Qw1waoK8-nQoSTPWQ;S8T|w?$SneZTMe){)-lSa&~CNN;?> z#*Dc7VHP?j`t$EjUjEw`#LJ|^_N?P-zf(SNQu+FRbJRc0#CoINI_#8_tGFJKR#b%n zW&ugb;_DafF=Y{(^TuyM& z6w%3d>lRU}OWcZH1{b~VIrd09DKyni_a!6omHAeKmgr$ZKiK)?th>tjMgYhGBcop%}_fK$vN0cgw&sppbCk%r7ws2|FYNFG;|eB zZbxY%Xi@>E*hIeb0~DoeXvDBJI5LrQm@TWE{|HI$ZjjDV!nHZ zDd475*j*0T1Q?%h*bTresMN0GJm=-rx3(4#Nowt}`AxcHMDu>nj4TNi_W1btdT53# ztc)g2cE5=f1xJVIz6H#w)N95vbFRgy5?9riYwjDPW?4+*tn<;XY$aMZhG0xRMJp}0 zitiH>2@aa}ce#4cri!Wn2V9))mWc@t4i+hOsh-jBi&0C>TlvagJeO})Y{?t6Hb)2; z0|~O7#r)ax=cxqQ8pubY4PHwttaEuJAgDKo+HZinnDxV$+mMSwL1AIzq@gm;A8Dri&|#w*ttg%rtZblvnwF@e3wnY&w<&T++Rk)hR}*E1`e}nsO)dXS+>|P!>9N5w%p6%Fss7B3FC>|ynr~>2( z(T@jUzz#5}&H*&xklHL+xAp=zi0G6?;EN3=?-O*&uXvp`!*Tl~R#2>PH3UuhwrYS9 z22=4#I27h;L%`p$1<~iBMO8)n3AKn`RScDb?1YPX+0-Nr;10BnmF+sxt~Xd)x5RG8 z?E3q0#nVTQ3~ke$E2HJjats6how}$-486?9ux}$Z!+3X3raXE!$8QG%~{e zz68V#{S({3&;m*Hfq_O>_8DAoRFo|yYL!8lM}2ngtKZBV?5pvqQ@+7oeM9xuQY`H3 zn(pp0XgnUdmt=Sbo0Bj)DBYIx#HqVd5r+?`i4RX##sc+6FwLy0Gee{&DmamC0E?CA zMjCn%2Bgb840jn42!L9~eUcYP3KwS?$fd2gB>~e2(Q5iswi^H(^4HbaPSu(3Fy*`K zL$~(p5yR|7uS1zsPgk6tPAF}JAUf9e@Q?)n@kehtF^^A}!=pvjS9bPIT^UJt8^g0y zT#Qpa6{s8HOcr}cNy9Wc93}Y;U4kR{y|C*ROa%^C7+K&u1ORW(vgL6|4R1ZW0~dxT z0tV<>WaM3)I6h*=3l%1;_4{jghf4s5)G;;P3R*5f$G~_P)-ttP?dy^G3YLufrnr+Z zm>^S%?b?+M1kmB}wlk$rB9~9ReY+JMT4K9*Qy8r?!)pNb2jRw+EnBWOotdjqn@;!q zxV%nBPY+mlBIcm+0}uhI5bym~(+2_64(LlHxp#_+60~xlsbIiry$=`~&3;?-9=xV~ zWpLW_0H_5x@aXC3Q}>qfadEwZi=?ZsKYGWxR&dXrwq!fGJ8;i7!_xn5nm>@@ZD$!L zRQ6dPdJKRU5hXQSTVex4!<%87oc<;sQ{Jrv%f)!@_tSpsM^{HL%}(Aa84~3lJxMFY zV>3!OJJ=@W9=(wONMy^H+lo>efl)Kv<5p=qrdaGY2_y6R9LTbcj*dS*m4>{3f6UFz zjYZvBhG!C!r9eS|GoItf19yO62p!4MCb(+Y$+Vj{9}plJpZ30~7uBQ=YYuo+ZT+p? zC1u09*?w2EO-)Trfar{p9V1L=`1g$90*U0d9>unf4m;dX0(A};G=ga|d{eYhQQv9D z7(c!rv8IY{J9Bj6^7N7=hz>crPha^$@_UxqxVR?zaL)mDx#TtI3CyB54k$c=%Ujj4 zmz?)0bZ@n6jO&|TWPC#;!hn__I>|w42F_fS(hWuoHV)nod=OBzmdZPEK~(o}kjt zq+_n77R)B;z!MY4N)g`)43k)*t-MwDE*M1)hN>%Pi)#OzH&mM330U@x9RA+vkM}sNzNyie5_09;zCo+gs%B(pA zZ=QI|oY#4kO=Ok-kumQFB$^4v8Z|RtblcdO_EcGecVP1*i5Lv?Z!FUg79V$YdVcnZ zSnMWPS}dA?ZiWr=cBcUWRNbMhoQ7l>E3i?R5hKvd ziG1ALt+GAljlii)0vH`Rm#c8z_|!cBU=~;@+_zV(Q%kB)6|2$-=d$RB4C&D9)#g&5 z6=|RCL?4k1fmXWQp%?+@K;-8TaSijMt!ZwrcNvUw&}&fwBu)G!dYgPvE9iw>O73D< zO8=yFfFgVf-2Xb!ws0E!>9gU;_uU3e z7$c;4t${$=I7p)2h4aX!ROYxzmDfyBt$2`hYF?@Hfa*qnq@3};pO=_UUaZ1YtVYu( zf5nM1B>t;v4C}$T+IR&j{qraFw8pi%Wumo(v~xjs;$69&OO9=1wfzl-3$ES< zrcP03RdLFBZiRm^iPm(ckb*&A$U}gy|h@>ix-QN%~EO(cDXn4*)w7&5|vP~eU zu4Q)Q_O}mQ+oc>_hS^h8w>@|N&1Nj7#BA9!hOiyXn6}xW#{V|^;K%_9sAYo%3rU;Y z76!i*SVf@CwWGNjucb3+!jOQRwI5crA6jl9(0Rbu+5lrV^USlOrID)u| zB?pymGU(c{zP`Rm%lk@8OW^>-5F2Yga&hjh4nN4w-|MO@ zhXA^x#U3R-BFCy-WmUh%g5jm8w*K3$4WXR0HS4Frhf``q5<^T(lG zUzH=6Ek$};Q8i_=!5r3>Kq3(|HXR}0Ft@;sUl{0@u=Jf2c6v6pb)#Q(Ek=!idIJ)9 zob9$P_l8q4Ue3$6yUlu)dU;LUdqnHgcj+~<6w3N7RN6^t!ymZhS=a?ARhGCXk_-O# z_OYgbkHHM=6ca0G{rE|*`=T7%mM!mgs=x&+DfwtmNifP7xc7w;Hc4&W-R_O?9%^}Q zAHs2Oj4ILycoaReQ0HYc;%qz>Ze7Ok<)XAkWzpc{Lx)t4-3lpvYyXPr*YsvD6}s+@ zr0^a8ej~y@y-&J%Q_H}>Ac@0J)__n1IP>AbaWdt#*+zmIW^nb=)N@0k&p$G9e)>Wz z*s~YqPTz5B+40w`GD%kWF0XT*n3o=%j+o$@;A8IhJ*aOq!Od*9*?(Q);mQ$T+Rqox zT@(mb8vJp2=ug$q(>c;J$c_B_`H_arKI7>vAdWEBV}!L60me_syUQW}-kj=B-?Q(L zMQhVCTp`=f6`jRhKlu5pA9<}W%^#*4&DqXsRCf&c_^q$dScN;53I4O`IYF2H`^HKK zz(^d19!*FH?fZ+BefPo59$@5izB@ZRYZ^l8>FGIYWd-ap@O5M*t|Qu@B&JB&P;8Cb z6pV8n_B=e#5xC#H_DOp(I0I{NkP%$pb>soh!5+^P!5@nvIToaVCJtDd^_8cVUFD*h z9#^_(X_>;ZtV^g~7!A}XcG?WOx=b>@s}dD3ZvInBsSz7;rBTpaJoEnl?t1Dt+`xN5 z8>Dt8JaU&VA1~#!^u7Y4yB~~IOzSp)S4+6wNEz+M!=|cE8%SEQZ;s>d1h(zAlE-dX zdSY6QbaDl8%}D0PnkZbsVWJ`lz~H%ttDA`ED9KRiI)Z{_{`q_qB#CQ2BlSu_MTE`! z)|XB7lc8S%X20%%)=++i5ThFFTb8YYF=|RuXLvQeZZ>4$B>KDk3351)ANLCi96a|( zNHoI};eQtj9}XQo@WO@&upvPbnf4S4wMH7Az>t6OUugf)U&To#NCm#?*{{_;gkFTb z4!T*1Z$pq7F^f!w^|aa78rTjxo5D=>!ZKXg_1Lw`e=V)+sH=j2J%2UHxnOt1eu>Z# zqg&dUS=QGkD4tWZ9&%Yb%u+=0oTOMXjc?}S) z0_uPv?g&`wZ_$&i(bSJ-P?;mB$Tx1V6VEb{*8(S)@gU04kW`78IUybjOLEQkw$&ujzD$!c$DkQ(8lf5rC}Y1!cZ_vw24{-mJKHvbJp1|FvIH(@_(&K`7#cMX1V zw|)ylz6F`Z@aPTc|K50o^H~#r%Q-8+t-!7{fT?dXdQ_Y{D z3oUS%;2~$r4;d?t3NHTd^TErYPHQf~iDvVy?^J#17%!#xbS*2&S$-~@3k=C_xM6bS zNCT*QWbnSv;<>Dq2_lphW?ybVZsNC*>VH0QA1(Q&lO=0V)&tTV0AXDtBT*HVbU4BZu;ge-P+R}S2K^{m#7>~1Z|vo0A%j0rd*wOW z@V@#l7;F_k^>Vl=>dwXD>iNh{4P9Cw0j^sa>TET5a<&%QQjuyeyk;6(q6Bs_I6ddr zj&;7a(4BIskb(A+)Il%qLsr=n7r|9W0f3nB!sQ_E2|ylzqDYkh7>jXF7Jy+b~&=~Doh*C1g*tgLGaHYUwX%@Q9$-p8-&f?$38kqc=UT%Ec!AS5KT0>!M!p5jdB zlg>>I(icANCB$jSaOLeH4LAFuz@`EMl30WK`hx`Y3;qsdWMJ|Rn0tIUPG;^mU0(ZC zg{CS9)GM%2D1#B^4Lm9Ss4xD!=rR~BQUQP5H|`0A|Gg89ogLsy9W}OK;xuQnx1cyA zNHW<{xy`*&hQxK;_Wp!Ff)k+WULv+5_@$s064@Qks({FQt!I@EX8+>b1~0iArLkT4 zX~UW2cvL9!yRWy~Y4ms?Tv$tj`J{+wJOq=%tUODG)MNgWC+rc8t_NEzx}S(_%NVRh zL5u(r0wNKN^GPn5(@duG?Zq`D>(1bdVyTHCOagJ)$*r{Yd@sBUg=+TEdQN8zM{<2- z4w3YakB4aP$;T@mLrTAZ%U-1br+B4=I% zvy4bW5lq~NWgWwWI|xZ<5wU3`~HeAptWK%4V2_Q>0Y>=0m;tI=o8OG*=h*S%D|xmu8R0fk zTwOBlOJ=l(@I!1P?B9KXYu+0VIq<=tzV5^m7$sSoa-B`o4Mq>1Co>t=ny8Z^o|GUX z@-sM$k)J`~la`K-;MMS{__81#PA(g^R^-i4!*~ZWy(`K=cxY1he)wsmm!V8k#Sd1{ zGk?_0$cWRk9iqkt$YoHAwB^ADf<76M z{EAaqlINxWO4Vd?FVdbM%(XqpwRKR6trf0U2i09O^`*{XxORwci8UU-dQ!0WK?J$pnEI z8I{!0LYQlrKc0ll83Vs+_4U$vRJ1HCrdQRJmCZL@{%KteJq3cmj?{K6BB=}lK@G1m zq=m%|Y;OYb?gn&vDc%n`LY0lZOxp;h0i;E#1V9~!1i)qd7@fVvv-r1JROx^K-~@tA zbQU$e*@>LEyHn4Dpu8&EH|>ez2ay@FbZ;GAo+&LjK~AqFRNY;HZFvpq&4fk~^Z-WHHddGk zi$=smD83kUUyF_XhSTmF;e)~-Q6NOpL%i)h_tB3kN3@vLaEPc z7T#?sF*T0!46}E=qsz_Xb4;(czYdfnHZ5`fk|9gi6vcTNLEmL$H_&vvH2=ZOYUTSe zB+ANscqLL7bdUaxueC_Z17pg)=R6qq3)`GMv5!-ID5JwT%X2#O|F<` zyKRox892!ApejrB;}Y}gsO@$*=BJ|`upeFsRVTeN65i&Q+hxrISH{e}hr;yl$hC7O zLtRLWX%h|tD3r`>;1R^AG#hu*kPu{Y5YZtuLC3@d(DK&e9(O-^XJt{zhK7dPC(U2) z(NqoSA$3SRYjE{;5s`qPAY;^5Dk>^soDnnWEs4T<-Wd6zqos9lTT9R;dMel>VjBe+ z)r!^VBxtH<>D*^1McBlcsl)DS!7TZFzn|~@`-2&n$S~n$s1+mgz2UBPFngkaR`6A4 zr>MMqGK!j8I5@`i@M~lMQ^rJ$3r7-pq0%!=^tj;T<3z%AD>5oAa6NI@=#hDkCg?o< z{-7}RCAkT&Dg@;^S4DU3yl3&p7)cG!s{SCas1ti{iG&MBG!8AibWHOsSv~!tiXxZg zDnv<;Hl$@{UgBKsHqIb|l&;T<+x(YU6@IvUP@i~fhQlq=bMW`tTaWFiexWe4Y14fL zypA6^0OTo#gIP_Hfb4hDp}m5#+?i6DsdM>w8Y3dRCtsc~|3hM?r#ERH6VipO6iTva zT`6yy4QXL(@zeK{3$ishc<^AEBTZGE@D?_<0KP&P>_2*pYX7>dXCHUlG?~2m!knK$ zDcWuv&L4J{@Ci{CdtS7vR0N|zq-x&M8g!#VC!nCPOx9S?@3bV2E=!_vjlEsy1$)T5 zmEa67QbOTzUK`ndG9zs#nXipo-6cKa?ECVgvqp1$K{9)sbJ@Ool2vY7h1u4j zK$7r9$vUQr#j5?JT{qb4=Zgi-RhfU7db3$@_eZOLkjV;e95UW+1`ILS#Qxos(b_pY z=mTJn-~b#PS=BRZ*juIp+)V`US&Ve-M`pnF#oL1?mi}Q`oK}}y%e+Uwd{%^F zq;F`u_Bxx9?OJEyr%eB?N33Qf=X-|S>r?rM#%BcNs5`3Iw}db2!OWTh}wOMQM9yEr^S|H6ARfe|;zAs)!f1n(W$(%NMS>X0F_k+{YoqA1epHcPLiqs^VY z67igjjJ5x;=AmOj>_NLvoe*>CH+c3l&z$x(BhOf>gJ6G4Nsvn>#gi#IgZ4bTX10PDY)A3Pc`P|Sx({4A)QyiEIb>s#4Hp+@ zQP-!1mQ|+vtdKKuecI(wb3M7J=hYYQUNd+=IJ4kSmcLDapqUoh0)yySQ_0(!I9cZU zIW_%%JC>cEP3C7Kt>|#kF#b|ApEJvvvVIh9#vgH4IYm*YdF9CK3U${mcxvtz5!&6g zJYMD=7WQy%2nI7;2w5eQY~w-Bg_gcu-7*IMW^BS)?jAmUsq_Tq?vSkqTFMzj6rY?N zRAGq_fwUz$HN4bIx!FcTI##vlbGwWkd*0iY?64=>57(4sZ=$)9mn~Qm{A5rdC?St` za#B$J;3Y=?b*6`>7?0XDJUPkow#f<~%DWTXFVtB?(oNylz%4}R^7_u9YzOL@^HwsE zE8Sc(?+?p{5wBzlWx*F1Skfe%wE$>)bs7p(+a1S@<5@d=1*B@li9}q#JFwX zALTD;v$M;!XN4sgF8x$dP+g~3_2WO|w?xTnhbGP3qV>6o^7MN8xe0# zH=QZ%%GKk1?m11LG1-f;q^JFl)sC2&He(PHF<=^+GXha3*v^La;Psd|sdy~iSEM2- z|ES-*ds^M4Sx)l$a>xB2Tlc$UD5rmqk_>sC*L=HK!L+!+la%}3VjuIuItn`z* z0iEd|Q3a=$nSVi4NQ10v&+hxFu1xx>F}WJRU9r-^WOt}hQ(tctlSX2}|z0G=# z(`8YZ(7qZLz!8m}8G7aP_aLfglE0evEJL>t){r2mpE5PY{;}u@(B3}tq_`RVDy4zz zzG%SxbhZwOq9-&4-MYHEgHWX*kIV1GsXz#jOioS`{_Eb8k3%6o!nUo@wE8&jcCi=} z*Y`HMKO#!)y{Ej+U(Bc-#}nV*Ms0OC^l94IBBZK>WvRpB#P}Q0(fs3@ zN)}X(JcmH6uoH6l?Dnu{;HrZNzqu0yAQ=WiEVRE2SLVCC)lxq`;7q!(VESXBNW`Jv zSKVA9n&&lFT%7xV26->M6s>P(u7BxB+JyyUxrXvCxbl!P zB||$Z&Q!ufUWT*kgH!e-q|V?IMDr?|G)Y0tGyBo9DqIl6P=oGDr*9IvC*QhcV&H?J z^4*H&353!_Cms!-ckqFk>LhsAi73JoLMV{Z*%*Sp-pl`x(2r&-m?+#5Z%IoHCui(w zWs}bIzFMby$AOyiZUI;0%r;kYrpnE2!km0kzgri-1WAo@4r*Ih51oo`WHOKa@{HtV zUb%OtanJpC9xvD0=&Gj-_nNhb(21UDYi&GXu5TPDC`a`zNkmT3l8vHeqqgaUX)H$d zYtAzceSUVrKI2MYllhmUA|qX!!q-o}|E|Q^l>2GN3u%6}yu+&D5^tXxaz)uJ>PL@f zgxK^ROS|xqb;5m!iB52KvPIg`T{Zl2hq~H`jJ~MR3l{MjY91l?^xnFDgHN@4zlYUJ zd$aUjl0Q@H&(wUYxQjuvk8)@f{Oql;o)6DC)TqGw{+c*MWWfJAsp$V{q7zN4KWU%dghWQ^iM^zV-7t~{ z1!CTjYOGV@o9X^IanqJqxRKx5^X+O**eA8;v4VMkGcq7YVd7<2%#}QujXiE zF4wT%%Ddk`cpuok+4;0X+rojCbk^fN*|rt3Fk^*15=0)M%N2CcXj7M`W6yz~0)2;HNHd2&0_UOn1Y|Ne0`r;23e14nySQ!^urfv=%~ zB*AdBXgy|*pzMiLvq|%PEzPZeI(tHF;<+J&X`v(}2+cmow4wxyUK=|&zS!A`@$oDh zjgS#1awSlAX+gUn(UnI;2M{{ifg&VyXy(zN-3=SNb?=FKBGdZ=*C(I<8kFWa=v+%5et~E4wms|oyE4*Y+NzN6WgWs& z;*>IjlcFW&{q(I`TZOYdET4gUNmzNzi;~ zd1dQoHx{#+2aykw@>ordYgl-5v0ex|fBb8dUyH^ zA@`yOZXhOU(GNW~on!pu5K`YnMu3p(Mu!hlW|tw8wFT8x_|fQ}BO(SuqxtnULgAG@ zEOI4hxv=XB3;A0m%^caP+cK1Q@4nbDCX>-4`?f2o)PaZG;i;Squ=3dF`ZDKAAVE+V$M~dGazHQ#-3#6wYuh z4qD5|XnHXIeE1_Tf?tGT{ZE#b_VAXhO9R)}KQs-xtwh1yV+uvr2BfWsV0qiNLLyUn z;97JvA$LLq=xN@=*YRx@Amu7%?HzA!$$qcy=u;^t zeuKx7s*mE*rb0Tw>>Ywx!S_Lt~dBkUBL60cIxp14Y_QF86AX3ovZ zTk2Odqo#lS==jn0Muk+`cgwQZ;DX16h#Aq_7u{wSzTEB)vRi0V?DR=S9xzTTYAg4e zZlX-{O|?|xK2s8djzfmqmX%ElWx{YY;C)Hd;@B!Q=^%9ms1&6muj zRUOkL)!|%xNf4*q)$n-M?Sz@$9 z+@8Hb!hd3VT3Ehg@%-UL)r$;4GT9YjYKjygJpm!HdnDR!{V4t9Q1W~KxXIUF62HmD zB4(gK#3xi2kkK5ew=1H)Ug31+dCrGfPZovGOP_N~+_oeWZ>QgyWUTV7;O!-gK2kRC z=Ah@6Rz~7fy9><{G7pIMsh+Alt7uVNqZ-(s5f3Dcx?AA&4Yy1AB-m5F`zjgMF6}rA| zmHE~w9JJ=z)w?b|HIkZ@CVS%__rP)L=cdJX=f?EzZO z)I}sE$no?EA6BvsWzkK9Y|O#Z6ObH^-+8?uZs%d{a5beY4s@t#NN>6F@za;g1Wv0R z!X}JYUU281J6*B)CjAPFSop#5*TwJ3zK{P>b{}xE+Rogsvl1cSNWL#lNX=d@y~j## zDJiaN{Z$CVfVT?exP+bFGkvPpacb@LF#~NChZA=mPV9{yieqvJ+q7`s=zH!anT?9r zZ<{Y1w(^Ogl^W5V7tM&j^jX5|$J?C;?P>jFf)ztcY*Q9A`y!SYbGJNGEZb!!Go?G~8)&H`%!Yhj$I%dh~fK@!lreq}Se(%9yV-%MA_2nJD{wooKZ@UHIO)F62aE zf{?M&g1gvgnRp+Ep<2xOR0^s2Pa{&wY2u&jnk%k1aV>p!{@`mXBowETP*ZK;^ih9a zv{Q(DxSs5hfwR3Hu6*a=OD|`2Rt=UV(_b7>SYSxiBnKQt{9maz)F`EADX9=G8AW4cV4sp1L5hnKS&t-cHG_3Z7jj(%d- zzSC=XPD(jZysv1#cB)=y_S0W!KORlLtCTxpQ5~%wl+L;%H(pBSx3P%ID!aSq**33t zzAGHmx2$x1vJCfYhsUXfuMLOoJjx`(#FS7pQ2%N{vTL>Z&anl${zs_)6sS z1Evr6%RSoq;ha_0T8{_Amcwx=hOSLvp+ypBY~zQ+V`Rny6Nkq;2ku`Q{cOc1QIkmJ zpzWZqmiFRBnQmSX?!_eCwLEuG#XoA5(Lee7Q^oso zhX#~eUX}M8CAEB;?n%{j3w_&fYMWVXIhpOUD$c$&qddpQ;o76k$!=Yz?-bXRZcm%Q0@Y?g<4J2mnV&7aq<7&k;(u^^{}ser)oG-YdFr|iq$4tW={WKw6X7{Gt~&Bv5|>_Y2~ z|CArQyYVT9yZX-__Mny7jt1tO%TbA~N!*iux~yvFkoji^)-Lwv7Uo@S)%)Q z8yPpc@)Q$;qYrX-ge3X4wx=arBG0m2Sz+g|jcZPpxVt!a|JdeD6%?nhA3LDLrQj=~ zv8VW6^>3#G6-Y#VbRt>K+*@ib;)QX1M(S*_87m=9p4yzB2U4_%`G4 m&zlwfKY#w;k2r6;LH&sN)52vdZ+9~MI;>@=S)^`#_5T1aQ=sPn literal 0 HcmV?d00001 diff --git a/docs/src/examples/Basics.md b/docs/src/examples/Basics.md index 13a5eb8..3e5ad31 100644 --- a/docs/src/examples/Basics.md +++ b/docs/src/examples/Basics.md @@ -185,7 +185,7 @@ Let's fit the actual model now using PointProcesses.jl and run some formal stati ### Fitting a Homogeneous Poisson Process ````@example Basics -pp_model = fit(PoissonProcess{Float64,Dirac{Nothing}}, h) +pp_model = fit(PoissonProcess{Float64,NoMarks}, h) println("Estimated rate λ̂ = ", pp_model.λ) # hide ```` diff --git a/docs/src/examples/CustomProcesses.md b/docs/src/examples/CustomProcesses.md new file mode 100644 index 0000000..1857cc1 --- /dev/null +++ b/docs/src/examples/CustomProcesses.md @@ -0,0 +1,302 @@ +```@meta +EditURL = "../../examples/CustomProcesses.jl" +``` + +# Defining Custom Processes + +Although several processes are already implemented and more are to come, we might need +to define our own custom processes. +Here we will explain the interface for defining new processes, starting with a simple case +and building up to more complex cases (which should not be too complicated). + +## Basic interface + +For this example we will implement the two-state model from [Selva2022](@cite). +This model takes three parameters, λₗ, λₕ and τ. The process starts +as a Poisson model with intensity λ = λₗ. Whenever an event occurs, the process +immediately switches to a Poisson process with intensity λ = λₕ for τ +units of time. If an event occurs during this high intensity phase, the "counter" τ simply +resets. You can visualize it below. + +![Adapted from Selva (2022)](../assets/selva2022.png) + +If the history ℋₜ up to time t contains the events t₁ < t₂ < ... < tₙ, then the equation +for the conditional intensity is + + λ(t | ℋₜ) = λₗ + 𝟙(tₙ > t - τ) (λₕ - λₗ), + +where 𝟙 is the indicator function. Notice that tₙ < t by definition of ℋₜ. + +We will use a different parametrization, namely λ = λₗ and Δ = λₕ - λₗ, so the constraint +λₕ > λₗ becomes Δ > 0. + +The first step in our implementation is to define this new struct by subtyping +the `AbstractUnivariateProcess` type (we will get to multivariate processes). + +````@example CustomProcesses +using PointProcesses +using Distributions +using Plots +using StatsAPI +using Optim + +struct TwoStateModel{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess + λ::R + Δ::R + τ::R + mark_dist::D +end +```` + +> Note: Subtype `AbstractUnivariateProcess`, not just `AbstractPointProcess`, this way we can use `TwoStateModel` to build multivariate processes with `IndependentMultivariateProcess`. + +`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())`. + +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. + +We will begin by implementing the methods `ground_intensity`, `integrated_ground_intensity` and +`ground_intensity_bound`. They are sufficient for almost all methods to function properly. + +````@example CustomProcesses +function PointProcesses.ground_intensity(pp::TwoStateModel, t, h::History) + if isempty(h) || t <= h.times[1] + return pp.λ + else + t_n = h.times[searchsortedfirst(h.times, t) - 1] # Last event time before `t` + in_high_phase = t < t_n + pp.τ + end + return pp.λ + (in_high_phase * pp.Δ) +end + +function PointProcesses.ground_intensity_bound(pp::TwoStateModel, t, h::History) + return (pp.λ + pp.Δ, Inf) # Also not the most efficient implementation, but works +end + +function PointProcesses.integrated_ground_intensity(pp::TwoStateModel, h, a, b) + integral = pp.λ * (b - a) + interval_events = event_times(h, a, b) + if length(interval_events) > 0 + for i in 1:length(interval_events) - 1 + integral += pp.Δ * max(interval_events[i + 1] - interval_events[i], pp.τ) + end + integral += pp.Δ * max(b - interval_events[end], pp.τ) + end + return integral +end +```` + +> Notice that in `ground_intensity` we searched for the last event that occurred before time `t`. We need to assume that `h` might be the full history, and not only the history up to time `t`. + +With this we can now call many other methods: `intensity`, `log_intensity`, `logdensityof`, `time_change` +and `simulate`. Here are the dependencies: +- `intensity` → `ground_intensity` + `densityof` from `mark_dist` +- `log_intensity` → `intensity` +- `logdensityof` → `integrated_ground_intensity` + `logdensityof` +- `time_change` → `integrated_ground_intensity` +- `simulate` → `ground_intensity` + `ground_intensity_bound` + +````@example CustomProcesses +pp = TwoStateModel(1.0, 2.0, 0.5, Normal()) +h = simulate(pp, 0.0, 20.0) + +println("intensity(pp, m, t, h) = $(intensity(pp, 0.0, 0.0, h))") +println("log_intensity(pp, m, t, h) = $(log_intensity(pp, 0.0, 0.0, h))") +println("logdensityof(pp, h) = $(logdensityof(pp, h))") +println("time_change(pp, h) = $(time_change(h, pp))") +println("simulate(pp, 0.0, 10.0) = $(simulate(pp, 0.0, 10.0))") +```` + +Let's plot our simulated process `h` along with its ground intensity. + +````@example CustomProcesses +sim_plot = plot( + event_times(h), + fill(1.0, nb_events(h)) + ; line=:stem, + grid=false, + label=false, + xlabel="Time", + ylabel="Events", + yaxis=false, + xlim=(min_time(h), max_time(h)) +) + +xs = LinRange(min_time(h), max_time(h), 1000) +intensity_plot = plot( + xs, + [ground_intensity(pp, x, h) for x in xs], + ; xaxis=false, + xgrid=false, + label=false, + ylabel="Ground Intensity", + xlim=(min_time(h), max_time(h)) +) + +plot( + intensity_plot, + sim_plot, + title=["Two State Model Simulation" ""], + layout=grid(2, 1, heights=[0.7, 0.3]) +) +```` + +Lastly, we can define the `fit` method for estimating the parameters from an observed +history. `logdensityof` is simply the log-likelihood of our process, so for a quick +and simple example we keep the mark distribution fixed and only optimize the core +process parameters. + +````@example CustomProcesses +function StatsAPI.fit( + ::Type{TwoStateModel{R1,NoMarks}}, + h::History, + init_params::Vector{R2}, +) where {R1<:Real,R2<:Real} + objective(params) = -logdensityof(TwoStateModel(params..., NoMarks()), h) + + lower_bound = [0.0, 0.0, 0.0] + upper_bound = [Inf, Inf, Inf] + result = optimize(objective, lower_bound, upper_bound, init_params) + optimal = Optim.minimizer(result) + + return TwoStateModel(optimal..., NoMarks()) +end +```` + +> It is tricky to properly estimate the parameters from this process. This implementation is conceptually correct, but it should be improved in real applications. + +With the `fit`, `simulate` and `time_change` methods, we can perform +goodness-of-fit tests for our model. + +````@example CustomProcesses +function random_params() + λ = 10.0 + rand() * 10 + Δ = λ * rand() + τ = inv(λ + Δ) * (0.5 + rand()) + return [λ, Δ, τ] +end + +true_params = random_params() +h_unknown = simulate(TwoStateModel(true_params..., NoMarks()), 0.0, 100.0) + +init_params = random_params() +pp_estimated = fit( + TwoStateModel{Float64,NoMarks}, + h_unknown, + init_params, +) +estimated_params = [pp_estimated.λ, pp_estimated.Δ, pp_estimated.τ] # hide + +println("True and estimated parameters:") # hide +println("λ: $(true_params[1]) - $(estimated_params[1])") #hide +println("Δ: $(true_params[2]) - $(estimated_params[2])") #hide +println("τ : $(true_params[3]) - $(estimated_params[3])") # hide + +pp_to_test = TwoStateModel(20.0, 30.0, 1.0/30.0, NoMarks()) +test_result = MonteCarloTest(KSDistance{Exponential}, pp_to_test, h_unknown) + +println("p-value for hypothesis test: $(pvalue(test_result))") # hide +```` + +> The statistics for goodness of fit tests currently implemented are suited only for testing the distribution of event times. They will run for marked processes, but the test results are only related to the event times. + +## Custom Mark Distributions + +Mark distributions can also be customized. Analogously to processes, we need +to subtype the `AbstractMarkDistribution` type. Lets make a mark distribution +that changes with times. + +````@example CustomProcesses +struct LinearTimeNormal{R<:Real} <: AbstractMarkDistribution + L::R + H::R +end +```` + +There are five methods that can be used with an `AbstractMarkDistribution`: `mark_distribution`, +`sample_mark`, `DensityInterface.densityof`, `Base.eltype` and `StatsAPI.fit`. The number of +methods we need to implement will depend on what `mark_distribution` returns. + +`mark_distribution` is used to return the distribution of the marks at a specific time `t` after +some history `h`. Say `mark_distribution(md::AbstractMarkDistribution, t, h::History)` is of +type `T`. +- `rand(::AbstractRNG, ::T)` implemented ⇒ `sample_mark(md, t, h)` works +- `eltype(::T)` implemented ⇒ `eltype(md)` works +- `densityof(::T, m)` implemented ⇒ `densityof(md, t, h, m)` works +- `fit` must always be implemented + +````@example CustomProcesses +function PointProcesses.mark_distribution(md::LinearTimeNormal, t, h::History) + time_ratio = (t - min_time(h)) / duration(h) + return Normal(md.L + (md.H - md.L) * time_ratio, 1) +end +```` + +In our case `mark_distribution` returns a `Distribution`, so `sample_mark`, `eltype` and +`densityof` are already take care of. To have all the functionality from the standard +mark distributions, we would only need to implement `fit`. + +````@example CustomProcesses +pp = TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)) + +println("TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)) = $(TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)))") +println("PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)) = $(PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)))") +println("log_intensity(pp, m, t, h) = $(log_intensity(pp, 0.0, 0.0, h))") +println("logdensityof(pp, h) = $(logdensityof(pp, h))") +println("simulate(pp, 0.0, 10.0) = $(simulate(pp, 0.0, 10.0))") +```` + +## Multivariate Processes + +Lastly we can define a multivariate process composed of multiple independent +univariate processes. If a method is available for all individual processes, +it will be also available for the multivariate processes. +All methods now return a vector, one for each dimension (each univariate process), +and we can add an integer argument `d` for accessing the value for a specific +dimension. + +````@example CustomProcesses +imp = IndependentMultivariateProcess([ + TwoStateModel(1.2, 2.0, 0.3, Normal()), + TwoStateModel(0.8, 2.0, 0.4, NoMarks()), + PoissonProcess(2.0), # Any type of process can be added +]) + +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))") +```` + +Lets end this tutorial with a nice plot of this process. + +````@example CustomProcesses +λ(t) = ground_intensity(imp, t, sim) # Returns a vector with the 3 individual intensities + +events_plot = scatter() +for d in 1:ndims(sim) + scatter!(events_plot, event_times(sim, d), fill(d, nb_events(sim, d)), label=nothing, markersize=3) +end +plot!(events_plot, yaxis=false, xlabel="Time", xlim=(min_time(sim), max_time(sim)), ylim=(0.0, 3.2)) + +xs = LinRange(min_time(sim), max_time(sim), 1000) +intensities_plot = plot() +for d in 1:ndims(sim) + plot!(intensities_plot, xs, getindex.(λ.(xs), d), label="Dimension $d") +end +plot!(intensities_plot, xs, sum.(λ.(xs)), label="Total ground intensity") +plot!(intensities_plot, xaxis=false, xlim=(min_time(sim), max_time(sim)), legend=:topleft) + +plot( + intensities_plot, + events_plot, + title="Event times and ground intensities", + layout=grid(2, 1, heights=[0.7, 0.3]), +) +```` + diff --git a/docs/src/examples/Hawkes.md b/docs/src/examples/Hawkes.md index d680a98..b3715f8 100644 --- a/docs/src/examples/Hawkes.md +++ b/docs/src/examples/Hawkes.md @@ -15,9 +15,9 @@ where μ(t) is the baseline intensity, which accounts for spontaneous events, `` which controls how past events influence future events. One of the most common choices for this kernel is the exponential function: ``\phi(s) = \alpha\exp{-\beta s}``. In this model, α is the strength of self-exciting, and β controls the rate of excitation decay. -Given these, parameters, the intensity function is: ``λ(t|\mathcal H_t) = μ + \sum_{t_i < t} \alpha \exp{-\beta (t - t_i)}``. We assume that μ is constant for simplicity. But, we could generalzie this to be inhomogeneous. +Given these parameters, the intensity function is: ``λ(t|\mathcal H_t) = μ + \sum_{t_i < t} \alpha \exp{-\beta (t - t_i)}``. We assume that μ is constant for simplicity. But, we could generalize this to be inhomogeneous. -An important statistic of the Hawkes proces, is the branching ratio which is a measure of the expected number of "daughter" events, a "parent" event will create. For the exponential kernel, this is given by the following equation: +An important statistic of the Hawkes process is the branching ratio which is a measure of the expected number of "daughter" events, a "parent" event will create. For the exponential kernel, this is given by the following equation: ``n = \frac{\alpha}{\beta}``. In the event that n < 1, the process is subcritical, meaning that events will eventually die out. If n = 1, the process is critical, and if n > 1, the process is supercritical, meaning that events can lead to an infinite cascade of events. Let's now apply this theory, and develop more through the process, by fitting a Hawkes process to some data. @@ -36,7 +36,8 @@ First let's open our data. This data records litter box entries taken from three ````@example Hawkes root = pkgdir(PointProcesses) -data = CSV.read(joinpath(root, "docs", "examples", "data", "cats.csv"), DataFrame) +data = CSV.read(joinpath(root, "docs", "examples", "data", "cats.csv"), DataFrame); +nothing #hide ```` This dataset has three cats in it, so let's first separate out each cat, using their weight to infer their identity @@ -46,7 +47,8 @@ cat_weights = [parse(Float64, split(i)[1]) for i in data.Value] clusters = kmeans(reshape(cat_weights, 1, :), 3; maxiter=100, display=:none) data.CatWeight = cat_weights -data.CatID = clusters.assignments +data.CatID = clusters.assignments; +nothing #hide ```` Now that we have set this up we need to get the timestamps in a useable order for plotting. We can use the Dates package for this. @@ -70,7 +72,9 @@ Since this dataset only has resolution up to the nearest minute, we need to chec ````@example Hawkes if any(diff(sort(data.t)) .== 0) - println("Data contains tied event times. Please add small jitter to event times to make them unique.") + println( + "Data contains tied event times. Please add small jitter to event times to make them unique.", + ) else println("Data contains no tied event times. Proceed with fitting. Yay!") end @@ -79,8 +83,23 @@ end Great! We can now move forward. First, let's visualize the data using a simple event plot. ````@example Hawkes -function eventplot(event_times::Vector{Float64}; title="Event Plot", xlabel="Time (minutes)", ylabel="Events") - scatter(event_times, ones(length(event_times)); markershape=:vline, markersize=10, label="", title=title, xlabel=xlabel, ylabel=ylabel, yticks=false) +function eventplot( + event_times::Vector{Float64}; + title="Event Plot", + xlabel="Time (minutes)", + ylabel="Events", +) + scatter( + event_times, + ones(length(event_times)); + markershape=:vline, + markersize=10, + label="", + title=title, + xlabel=xlabel, + ylabel=ylabel, + yticks=false, + ) end cat1_times = data.t[data.CatID .== 1] @@ -106,7 +125,7 @@ n_days = length(unique(day)) h_day = History(sort(Float64.(tod)), 0.0, 1440.0) # build a "history" on [0, 1440] minutes nbins = 96 # 96 bins = 15-minute bins pp_day = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h_day, nbins, ) @@ -131,10 +150,7 @@ does that increase the likelihood of another cat using it soon after? To do this ````@example Hawkes full_history = History(data.t, 0.0, maximum(data.t) + 1.0) -hawkes_model = fit( - HawkesProcess, - full_history, -) +hawkes_model = fit(HawkesProcess, full_history) println("Fitted Hawkes Process Parameters:") # hide println("Base intensity (μ): ", hawkes_model.μ) # hide @@ -153,17 +169,21 @@ Finally, we can visualize the fitted intensity function over time. ````@example Hawkes ts = sort(data.t) -λ_hawkes(t::Real) = +function λ_hawkes(t::Real) hawkes_model.μ + sum((hawkes_model.α * exp(-hawkes_model.ω * (t - ti)) for ti in ts if ti < t); init=0.0) +end u = range(0.0, maximum(ts) + 1.0; length=2000) -plot(u, λ_hawkes.(u); - xlabel="Time (minutes)", - ylabel="Fitted Hawkes intensity (events/min)", - title="Fitted Hawkes Process Intensity Function", - legend=false) +plot( + u, + λ_hawkes.(u); + xlabel="Time (minutes)", + ylabel="Fitted Hawkes intensity (events/min)", + title="Fitted Hawkes Process Intensity Function", + legend=false, +) ```` From the plot, we can observe how the intensity function varies over time, capturing the self-exciting nature of the litter box entries. @@ -176,11 +196,15 @@ test = MonteCarloTest(KSDistance{Exponential}, hawkes_model, full_history; n_sim p = pvalue(test) # hide println("Monte Carlo Test p-value for Hawkes Process fit: ", p) # hide -println("Assuming a significance level of 0.05, we " * (p < 0.05 ? "reject" : "fail to reject") * "",) # hide +println( + "Assuming a significance level of 0.05, we " * + (p < 0.05 ? "reject" : "fail to reject") * + "", +) # hide println("the null hypothesis that the Hawkes process is a good fit to the data.") # hide ```` -From this analysis, it seems that we can can say that litter-box usage is indeed a self-exciting process, as the fitted Hawkes process provides a good fit to the data. +From this analysis, it seems that we can say that litter-box usage is indeed a self-exciting process, as the fitted Hawkes process provides a good fit to the data. It's worth considering that we have ignored the identities of the cats in this analysis. A more thorough analysis could involve fitting a multivariate Hawkes process, where each cat is represented as a separate dimension. This would allow us to capture the interactions between the cats more accurately. This will be the subject of a future tutorial. diff --git a/docs/src/examples/Inhomogeneous.md b/docs/src/examples/Inhomogeneous.md index c0b8a08..93d3dcd 100644 --- a/docs/src/examples/Inhomogeneous.md +++ b/docs/src/examples/Inhomogeneous.md @@ -92,7 +92,7 @@ The simplest approach is to bin the data and estimate a constant rate in each bi ````@example Inhomogeneous pp_piecewise = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}}, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h, 20, # number of bins ) diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 8000823..9e5b630 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -110,4 +110,16 @@ @Book{Cox2018 isbn = {9780203743034}, month = dec, doi = {10.1201/9780203743034}, +} + +@Article{Selva2022, + author = {Jacopo Selva and Laura Sandri and Matteo Taroni and Roberto Sulpizio and Pablo Tierz and Antonio Costa}, + journal = {Science Advances}, + title = {A simple two-state model interprets temporal modulations in eruptive activity and enhances multivolcano hazard quantification}, + year = {2022}, + month = {nov}, + number = {44}, + volume = {8}, + doi = {10.1126/sciadv.abq4415}, + publisher = {American Association for the Advancement of Science ({AAAS})}, } \ No newline at end of file diff --git a/src/PointProcesses.jl b/src/PointProcesses.jl index bbfeb02..5bbed1d 100644 --- a/src/PointProcesses.jl +++ b/src/PointProcesses.jl @@ -38,6 +38,10 @@ export min_time, max_time, min_mark, max_mark export nb_events, has_events, duration, ndims export time_change, split_into_chunks +## Mark Distributions +export AbstractMarkDistribution, PointProcessMarkDistribution, NoMarks +export sample_mark + ## Point processes export PointProcessMarkDistribution, AbstractMarkDistribution, NoMarks diff --git a/src/abstract_point_process.jl b/src/abstract_point_process.jl index 9209ff7..3d4b67f 100644 --- a/src/abstract_point_process.jl +++ b/src/abstract_point_process.jl @@ -67,7 +67,9 @@ intensity(pp, h, t, d) computes the intensity for a multivariate process `pp` at 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 +function intensity(pp::AbstractPointProcess, m, t, h) + return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m) +end """ log_intensity(pp, m, t, h) @@ -78,7 +80,7 @@ For multivariate processes, it returns a vector of log intensities for each dime log_intensity(pp, h, t, d) computes the log intensity for a multivariate process `pp` at dimension `d`. """ function log_intensity(pp::AbstractPointProcess, m, t, h) - return log(ground_intensity(pp, t, h)) + logdensityof(mark_distribution(pp, t, h), m) + return log(intensity(pp, m, t, h)) end ## Simulation @@ -147,3 +149,8 @@ Fit a point process of type `PP` to one or several histories using maximum a pos Not implemented by default. """ function fit_map end + +function time_change(h::History, pp::AbstractPointProcess) + Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t) + return time_change(h, Λ) +end \ No newline at end of file diff --git a/src/history.jl b/src/history.jl index 8816694..c09b105 100644 --- a/src/history.jl +++ b/src/history.jl @@ -78,8 +78,9 @@ struct History{T<:Real,M,D} end function History( - times::Vector{Vector{R1}}, tmin::R2, tmax::R3, marks::Vector{Vector{M}}; check_args=true -) where {R1<:Real,R2<:Real,R3<:Real,M} + times::AbstractVector{<:AbstractVector{R1}}, tmin::R2, tmax::R3, + marks::AbstractVector{<:AbstractVector}; check_args=true +) where {R1<:Real,R2<:Real,R3<:Real} vec_times = vcat(times...) vec_marks = vcat(marks...) perm = sortperm(vec_times) @@ -101,7 +102,7 @@ function History( end function History( - times::Vector{Vector{R1}}, tmin::R2, tmax::R3; check_args=true + times::AbstractVector{<:AbstractVector{R1}}, tmin::R2, tmax::R3; check_args=true ) where {R1<:Real,R2<:Real,R3<:Real} nots = [fill(nothing, length(times[i])) for i in 1:length(times)] return History(times, tmin, tmax, nots; check_args=check_args) @@ -130,6 +131,12 @@ function History(tmin::R1, tmax::R2, N::Int=1) where {R1<:Real, R2<:Real} History(R[], tmin, tmax, [], [], N) end +function History(h::History, d::Int) + times = event_times(h, d) + marks = event_times(h, d) + return History(times, min_time(h), max_time(h), marks) +end + """ event_times(h) diff --git a/src/mark_distributions.jl b/src/mark_distributions.jl index 392902f..0ed809d 100644 --- a/src/mark_distributions.jl +++ b/src/mark_distributions.jl @@ -1,5 +1,6 @@ # Type definition +"Abstract type for defining mark distributions not in `Distributions.jl`" abstract type AbstractMarkDistribution end const PointProcessMarkDistribution = Union{Distribution, AbstractMarkDistribution} @@ -26,17 +27,22 @@ sample_mark(rng::AbstractRNG, md::PointProcessMarkDistribution, t, h::History) = sample_mark(md::PointProcessMarkDistribution, t, h::History) = sample_mark(default_rng(), md, t, h) +"The type of the marks returned by the mark distribution" Base.eltype(md::AbstractMarkDistribution) = eltype(mark_distribution(md, 0.0, History(0.0, 1.0))) +"The likelihood of a mark `m` occurring in an event at time `t` after history `h`" DensityInterface.densityof(md::PointProcessMarkDistribution, t, h::History, m) = densityof(mark_distribution(md, t, h), m) +StatsAPI.fit + # Support for `Distributions.jl` mark_distribution(d::Distribution, t, h::History) = d StatsAPI.fit(D::Type{<:Distribution}, h::History) = fit(D, h.marks) # Struct for non-marked processes +"Mark distribution for non-marked processes. Always return the mark `nothing`." struct NoMarks <: AbstractMarkDistribution end mark_distribution(::NoMarks, t, h::History) = Dirac(nothing) @@ -49,6 +55,6 @@ StatsAPI.fit(::Type{NoMarks}, h) = NoMarks() StatsAPI.fit(::Type{NoMarks}, marks, weights) = NoMarks() -DensityInterface.densityof(::NoMarks, t, h, ::Nothing) = 1.0 +DensityInterface.densityof(::NoMarks, t, ::History, ::Nothing) = 1.0 -DensityInterface.densityof(::NoMarks, t, h, m) = 0.0 +DensityInterface.densityof(::NoMarks, t, ::History, m) = 0.0 diff --git a/src/multivariate/independent_multivariate.jl b/src/multivariate/independent_multivariate.jl index 193c288..affebc0 100644 --- a/src/multivariate/independent_multivariate.jl +++ b/src/multivariate/independent_multivariate.jl @@ -6,8 +6,7 @@ A multivariate point process where each dimension is an independent univariate p # Fields - `processes::Vector{P}`: vector of univariate point processes, one for each dimension. """ -struct IndependentMultivariateProcess{P<:AbstractUnivariateProcess} <: - AbstractMultivariateProcess +struct IndependentMultivariateProcess{P<:AbstractUnivariateProcess} <: AbstractMultivariateProcess processes::Vector{P} end @@ -19,7 +18,7 @@ Base.ndims(pp::IndependentMultivariateProcess) = length(pp.processes) ## AbstractPointProcess interface function ground_intensity(pp::IndependentMultivariateProcess, t, h, d) - return ground_intensity(pp.processes[d], t, h) + return ground_intensity(pp.processes[d], t, History(h, d)) end function ground_intensity(pp::IndependentMultivariateProcess, t, h) @@ -27,7 +26,7 @@ function ground_intensity(pp::IndependentMultivariateProcess, t, h) end function mark_distribution(pp::IndependentMultivariateProcess, t, h, d) - return mark_distribution(pp.processes[d], t, h) + return mark_distribution(pp.processes[d], t, History(h, d)) end function mark_distribution(pp::IndependentMultivariateProcess, t, h) @@ -39,7 +38,7 @@ function mark_distribution(pp::IndependentMultivariateProcess, t) end function intensity(pp::IndependentMultivariateProcess, m, t, h, d) - return intensity(pp.processes[d], m, t, h) + return intensity(pp.processes[d], m, t, History(h, d)) end function intensity(pp::IndependentMultivariateProcess, m, t, h) @@ -49,10 +48,11 @@ end function log_intensity(pp::IndependentMultivariateProcess, m, t, h, d) log(intensity(pp, m, t, h, d)) end + log_intensity(pp::IndependentMultivariateProcess, m, t, h) = log.(intensity(pp, m, t, h)) function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h, d) - return ground_intensity_bound(pp.processes[d], t, h) + return ground_intensity_bound(pp.processes[d], t, History(h, d)) end function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h) @@ -60,7 +60,7 @@ function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h) end function integrated_ground_intensity(pp::IndependentMultivariateProcess, h, a, b, d) - return integrated_ground_intensity(pp.processes[d], h, a, b) + return integrated_ground_intensity(pp.processes[d], History(h, d), a, b) end function integrated_ground_intensity(pp::IndependentMultivariateProcess, h, a, b) @@ -77,7 +77,7 @@ end function time_change(h::History{R,M}, pp::IndependentMultivariateProcess) where {R<:Real,M} histories = [ - time_change(History(event_times(h, d), h.tmin, h.tmax), pp.processes[d]) for + time_change(History(h, d), pp.processes[d]) for d in 1:ndims(pp) ] tmax = maximum(histories[d].tmax for d in 1:ndims(pp)) @@ -103,7 +103,7 @@ function StatsAPI.fit(proc_types::Vector, h::History{R,M}; kwargs...) where {R<: processes = [ fit( proc_types[d], - History(event_times(h, d), h.tmin, h.tmax, event_marks(h, d)); + History(h, d); kwargs..., ) for d in eachindex(proc_types) ] diff --git a/src/univariate/poisson/inhomogeneous/fit.jl b/src/univariate/poisson/inhomogeneous/fit.jl index a38c1a2..64f2ebb 100644 --- a/src/univariate/poisson/inhomogeneous/fit.jl +++ b/src/univariate/poisson/inhomogeneous/fit.jl @@ -80,18 +80,6 @@ function StatsAPI.fit( return from_params(F, optimal_params; intensity_kwargs...) end -function StatsAPI.fit( - ::Type{<:InhomogeneousPoissonProcess{F,NoMarks}}, - h::History, - init_params; - integration_config=IntegrationConfig(), - kwargs..., -) where {F<:ParametricIntensity} - # For unmarked processes, skip fitting the mark distribution - intensity = fit(F, h, init_params; integration_config=integration_config, kwargs...) - return InhomogeneousPoissonProcess(intensity, NoMarks()) -end - function StatsAPI.fit( ::Type{<:InhomogeneousPoissonProcess{F,D}}, h::History, @@ -100,7 +88,7 @@ function StatsAPI.fit( kwargs..., ) where {F<:ParametricIntensity,D} # Fit mark distribution independently - mark_dist = fit(D, h.marks) + mark_dist = fit(D, h) intensity = fit(F, h, init_params; integration_config=integration_config, kwargs...) return InhomogeneousPoissonProcess(intensity, mark_dist) @@ -126,32 +114,6 @@ function StatsAPI.fit( return fit(pptype, combined, args...) end -function StatsAPI.fit( - ::Type{<:InhomogeneousPoissonProcess{PiecewiseConstantIntensity{R},NoMarks}}, - h::History, - breakpoints::Vector; - kwargs..., -) where {R} - # Count events in each bin - n_bins = length(breakpoints) - 1 - rates = Vector{R}(undef, n_bins) - for i in 1:n_bins - # Count events in bin [breakpoints[i], breakpoints[i+1]) - bin_start = breakpoints[i] - bin_end = breakpoints[i + 1] - bin_width = bin_end - bin_start - - # Count events in this bin - count = sum(bin_start .<= h.times .< bin_end) - - # Estimate rate as count / duration - rates[i] = count / bin_width - end - - intensity = PiecewiseConstantIntensity(breakpoints, rates) - return InhomogeneousPoissonProcess(intensity, NoMarks()) -end - function StatsAPI.fit( ::Type{<:InhomogeneousPoissonProcess{PiecewiseConstantIntensity{R},D}}, h::History, @@ -159,7 +121,7 @@ function StatsAPI.fit( kwargs..., ) where {R,D} # Fit mark distribution independently - mark_dist = fit(D, h.marks) + mark_dist = fit(D, h) # Count events in each bin n_bins = length(breakpoints) - 1 diff --git a/test/mark_distributions.jl b/test/mark_distributions.jl new file mode 100644 index 0000000..2f70b02 --- /dev/null +++ b/test/mark_distributions.jl @@ -0,0 +1,26 @@ +@testset "Distributions.jl" begin + h = History(rand(10), 0, 1, rand(10)) + md = Normal() + @test md isa PointProcessMarkDistribution + @test mark_distribution(md, 0.0, h) == Normal() + @test fit(Normal, h) isa Normal + @test fit(Exponential, h) isa Exponential + @test densityof(md, 0.0, h, 0.0) == densityof(md, 0.0) + sample1 = sample_mark(Random.seed!(1), md, 0.0, h) + sample2 = rand(Random.seed!(1), Normal()) + @test sample1 == sample2 +end + +@testset "NoMarks" begin + h = History(rand(10), 0, 1, rand(10)) + md = NoMarks() + @test md isa AbstractMarkDistribution + @test md isa PointProcessMarkDistribution + @test mark_distribution(md, 0.0, h) isa Dirac{Nothing} + @test fit(NoMarks, h) == NoMarks() + @test fit(NoMarks, rand(10), rand(10)) == NoMarks() + @test densityof(md, 0.0, h, 0.0) == 0.0 + @test densityof(md, 0.0, h, nothing) == 1.0 + @test sample_mark(md, 0.0, h) === nothing + @test eltype(NoMarks()) == Nothing +end \ No newline at end of file diff --git a/test/multivariate_poisson_process.jl b/test/multivariate_poisson_process.jl index c0c6a58..c59cbca 100644 --- a/test/multivariate_poisson_process.jl +++ b/test/multivariate_poisson_process.jl @@ -37,10 +37,10 @@ end @test integrated_ground_intensity(pp1, h, 0, 1) == [1.0, 2.0, 3.0] @test integrated_ground_intensity(pp1, h, 0, 1, 1) == 1.0 - h1 = simulate(rng, pp, 0.0, 1000.0) + h1 = simulate(rng, pp2, 0.0, 1000.0) f1(λ) = logdensityof(PoissonProcess(λ), h1) - gf = ForwardDiff.gradient(f1, 3 * ones(10)) + gf = ForwardDiff.gradient(f1, 10 * ones(10)) @test all(gf .< 0) end @@ -67,18 +67,18 @@ end end @testset "Fitting" begin - h1 = simulate(rng, pp1, 0.0, 1000.0) - h2 = History(fill(rand(100) .* 1000, 3), 0, 1000, fill(rand(100), 3)) + h1 = simulate(rng, pp2, 0.0, 1000.0) + h2 = History(fill(rand(100) .* 1000, 3), 0.0, 1000.0) - pp_est1 = fit(fill(PoissonProcess{Float64,Normal}, 3), h1) - pp_est2 = fit(fill(PoissonProcess{Float64,Normal}, 3), h2) + pp_est1 = fit(fill(PoissonProcess{Float64,NoMarks}, 3), h1) + pp_est2 = fit(fill(PoissonProcess{Float64,NoMarks}, 3), h2) λ_est1 = [pp_est1.processes[d].λ for d in 1:ndims(pp_est1)] λ_est2 = [pp_est2.processes[d].λ for d in 1:ndims(pp_est2)] λ_error1 = mean(abs, λ_est1 - λ) λ_error2 = mean(abs, λ_est2 - λ) - l = logdensityof(pp, h1) + l = logdensityof(pp2, h1) l_est = logdensityof(pp_est1, h1) @test λ_error1 < λ_error2 diff --git a/test/runtests.jl b/test/runtests.jl index 762da97..7d17495 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,6 +31,9 @@ DocMeta.setdocmeta!(PointProcesses, :DocTestSetup, :(using PointProcesses); recu @testset verbose = true "History" begin include("history.jl") end + @testset verbose = false "MarkDistributions" begin + include("mark_distributions.jl") + end @testset verbose = true "Bounded" begin include("bounded_point_process.jl") end From de48b7541ebf48ebb9597cb1bbf23e36f1f27f80 Mon Sep 17 00:00:00 2001 From: jkling Date: Mon, 8 Jun 2026 09:53:34 +0200 Subject: [PATCH 3/6] Fixes for PR --- docs/examples/CustomProcesses.jl | 12 ++++---- src/PointProcesses.jl | 1 - src/history.jl | 44 ++++++++++++++++++++-------- src/mark_distributions.jl | 3 -- src/simulation.jl | 2 +- src/univariate/poisson/simulation.jl | 25 ++++++---------- test/history.jl | 23 ++++++++++++--- test/hypothesis_tests.jl | 8 ++--- test/mark_distributions.jl | 2 +- test/multivariate_poisson_process.jl | 2 +- 10 files changed, 74 insertions(+), 48 deletions(-) diff --git a/docs/examples/CustomProcesses.jl b/docs/examples/CustomProcesses.jl index f458357..fe21b52 100644 --- a/docs/examples/CustomProcesses.jl +++ b/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. diff --git a/src/PointProcesses.jl b/src/PointProcesses.jl index 5a7d717..9315f5b 100644 --- a/src/PointProcesses.jl +++ b/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 diff --git a/src/history.jl b/src/history.jl index be761f0..8a1997e 100644 --- a/src/history.jl +++ b/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) diff --git a/src/mark_distributions.jl b/src/mark_distributions.jl index 224e5b9..b38363b 100644 --- a/src/mark_distributions.jl +++ b/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 diff --git a/src/simulation.jl b/src/simulation.jl index 928c265..69bff15 100644 --- a/src/simulation.jl +++ b/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) diff --git a/src/univariate/poisson/simulation.jl b/src/univariate/poisson/simulation.jl index cbcc763..60f565a 100644 --- a/src/univariate/poisson/simulation.jl +++ b/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 diff --git a/test/history.jl b/test/history.jl index 0ff55fb..a164e72 100644 --- a/test/history.jl +++ b/test/history.jl @@ -1,9 +1,10 @@ @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) @@ -11,6 +12,9 @@ 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] diff --git a/test/hypothesis_tests.jl b/test/hypothesis_tests.jl index d71dbca..307a51d 100644 --- a/test/hypothesis_tests.jl +++ b/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) diff --git a/test/mark_distributions.jl b/test/mark_distributions.jl index 9b03ded..f0abd21 100644 --- a/test/mark_distributions.jl +++ b/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 \ No newline at end of file +end diff --git a/test/multivariate_poisson_process.jl b/test/multivariate_poisson_process.jl index c59cbca..151d713 100644 --- a/test/multivariate_poisson_process.jl +++ b/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) From 257f73a2a97c035a445b6518483b63eb3ed23a44 Mon Sep 17 00:00:00 2001 From: jkling Date: Tue, 9 Jun 2026 14:22:47 +0200 Subject: [PATCH 4/6] Fixed underflow problem in simulation --- src/abstract_point_process.jl | 12 ++++++++++++ src/history.jl | 4 ++-- src/simulation.jl | 14 +------------- src/univariate/poisson/simulation.jl | 2 +- test/history.jl | 8 ++++++-- test/hypothesis_tests.jl | 2 +- 6 files changed, 23 insertions(+), 19 deletions(-) diff --git a/src/abstract_point_process.jl b/src/abstract_point_process.jl index 263d154..72810f2 100644 --- a/src/abstract_point_process.jl +++ b/src/abstract_point_process.jl @@ -154,3 +154,15 @@ function time_change(h::History, pp::AbstractPointProcess) Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t) return time_change(h, Λ) end + +""" + simulate([rng,] pp, tmin, tmax) + +Generate an event record from the point process `pp` on the interval [tmin, tmax). + +# Technical Remark +Simulating events using types other than `Float64` may cause underflow issues, causing an error to be raised because of repeated event times. This is because the `rand` method from `Distributions.jl` always returns `Float64` values. +""" +function simulate(rng::AbstractRNG, pp::AbstractPointProcess, tmin, tmax) + return simulate_ogata(rng, pp, tmin, tmax) +end diff --git a/src/history.jl b/src/history.jl index 8a1997e..b125f5f 100644 --- a/src/history.jl +++ b/src/history.jl @@ -191,7 +191,7 @@ 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}) +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::Union{Int, Nothing}) +function event_marks(h::History, d::Union{Int,Nothing}) h.N == 1 && d == 1 ? h.marks : (@view h.marks[h.dims .== d]) end diff --git a/src/simulation.jl b/src/simulation.jl index 69bff15..dfed8ca 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -15,9 +15,6 @@ end simulate_ogata(rng, pp, tmin, tmax) Simulate a temporal point process `pp` on interval `[tmin, tmax)` using Ogata's algorithm. - -# Technical Remark -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::AbstractUnivariateProcess, tmin::T, tmax::T @@ -27,7 +24,7 @@ function simulate_ogata( t = tmin while t < tmax B, L = ground_intensity_bound(pp, t + eps(t), h) - τ = B > 0 ? rand(rng, Exponential(inv(B))) : typemax(inv(B)) + τ = B > 0 ? T(rand(rng, Exponential(inv(B)))) : typemax(inv(B)) if τ > L t = t + L elseif τ <= L @@ -45,15 +42,6 @@ function simulate_ogata( return h end -""" - simulate([rng,] pp, tmin, tmax) - -Alias for `simulate_ogata`. -""" -function simulate(rng::AbstractRNG, pp::AbstractPointProcess, tmin, tmax) - return simulate_ogata(rng, pp, tmin, tmax) -end - function simulate(pp::AbstractPointProcess, args...; kwargs...) return simulate(default_rng(), pp, args...; kwargs...) end diff --git a/src/univariate/poisson/simulation.jl b/src/univariate/poisson/simulation.jl index 60f565a..bf4872e 100644 --- a/src/univariate/poisson/simulation.jl +++ b/src/univariate/poisson/simulation.jl @@ -1,7 +1,7 @@ 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) + t = T(rand(rng, inter_dist)) while t < tmax m = sample_mark(pp.mark_dist, t, h) push!(h, t, m; check_args=false) diff --git a/test/history.jl b/test/history.jl index a164e72..753057e 100644 --- a/test/history.jl +++ b/test/history.jl @@ -90,8 +90,12 @@ end @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_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] diff --git a/test/hypothesis_tests.jl b/test/hypothesis_tests.jl index 307a51d..32f5795 100644 --- a/test/hypothesis_tests.jl +++ b/test/hypothesis_tests.jl @@ -1,6 +1,6 @@ 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{Float32,NoMarks} pp = PoissonProcess() @testset "Statistics" begin From 2da3339f3ca36b5cf61000c716b07a83114f2014 Mon Sep 17 00:00:00 2001 From: Ryan Senne <50930199+rsenne@users.noreply.github.com> Date: Tue, 16 Jun 2026 16:23:45 -0400 Subject: [PATCH 5/6] Format --- docs/examples/CustomProcesses.jl | 2 +- src/mark_distributions.jl | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/examples/CustomProcesses.jl b/docs/examples/CustomProcesses.jl index fe21b52..4c81323 100644 --- a/docs/examples/CustomProcesses.jl +++ b/docs/examples/CustomProcesses.jl @@ -177,7 +177,7 @@ println("λ: $(true_params[1]) - $(estimated_params[1])") #hide println("Δ: $(true_params[2]) - $(estimated_params[2])") #hide println("τ : $(true_params[3]) - $(estimated_params[3])") # hide -pp_to_test = TwoStateModel(20.0, 30.0, 1.0/30.0, NoMarks()) +pp_to_test = TwoStateModel(20.0, 30.0, 1.0 / 30.0, NoMarks()) test_result = MonteCarloTest(KSDistance{Exponential}, pp_to_test, h_unknown) println("p-value for hypothesis test: $(pvalue(test_result))") # hide diff --git a/src/mark_distributions.jl b/src/mark_distributions.jl index b38363b..9a1612b 100644 --- a/src/mark_distributions.jl +++ b/src/mark_distributions.jl @@ -13,7 +13,7 @@ Compute the distribution of marks at time `t` after history `h`. function mark_distribution end function mark_distribution(md::AbstractMarkDistribution, t, h::History) - error( + return error( "Type $(typeof(md)) subtypes `AbstractMarkDistribution` but has " * "not implemented the required `mark_distribution(md, t, h)` method.", ) @@ -25,21 +25,21 @@ end Return one sample from the distribution of marks at time `t` after history `h`, using the random number generator `rng`. """ function sample_mark(rng::AbstractRNG, md::PointProcessMarkDistribution, t, h::History) - rand(rng, mark_distribution(md, t, h)) + return rand(rng, mark_distribution(md, t, h)) end function sample_mark(md::PointProcessMarkDistribution, t, h::History) - sample_mark(default_rng(), md, t, h) + return sample_mark(default_rng(), md, t, h) end "The type of the marks returned by the mark distribution" function Base.eltype(md::AbstractMarkDistribution) - eltype(mark_distribution(md, 0.0, History(0.0, 1.0))) + return eltype(mark_distribution(md, 0.0, History(0.0, 1.0))) end "The likelihood of a mark `m` occurring in an event at time `t` after history `h`" function DensityInterface.densityof(md::PointProcessMarkDistribution, t, h::History, m) - densityof(mark_distribution(md, t, h), m) + return densityof(mark_distribution(md, t, h), m) end # Support for `Distributions.jl` From 9544cf53bc5a7928762777f9d89a8791959f7939 Mon Sep 17 00:00:00 2001 From: Jose Kling Date: Fri, 19 Jun 2026 08:52:08 -0300 Subject: [PATCH 6/6] PR Fixes - Added tests - Fixed `CustomProcesses` example in docs --- docs/examples/CustomProcesses.jl | 9 ++- docs/src/examples/Basics.md | 2 +- docs/src/examples/CustomProcesses.md | 92 +++++++++++++++++----------- docs/src/examples/Hawkes.md | 13 ++-- docs/src/examples/Inhomogeneous.md | 2 +- src/mark_distributions.jl | 4 +- test/bounded_point_process.jl | 17 +++++ test/history.jl | 4 ++ test/multivariate_poisson_process.jl | 7 +++ 9 files changed, 101 insertions(+), 49 deletions(-) diff --git a/docs/examples/CustomProcesses.jl b/docs/examples/CustomProcesses.jl index 4c81323..c8119e5 100644 --- a/docs/examples/CustomProcesses.jl +++ b/docs/examples/CustomProcesses.jl @@ -76,9 +76,9 @@ function PointProcesses.integrated_ground_intensity(pp::TwoStateModel, h, a, b) interval_events = event_times(h, a, b) if length(interval_events) > 0 for i in 1:(length(interval_events) - 1) - integral += pp.Δ * max(interval_events[i + 1] - interval_events[i], pp.τ) + integral += pp.Δ * min(interval_events[i + 1] - interval_events[i], pp.τ) end - integral += pp.Δ * max(b - interval_events[end], pp.τ) + integral += pp.Δ * min(b - interval_events[end], pp.τ) end return integral end @@ -178,7 +178,8 @@ println("Δ: $(true_params[2]) - $(estimated_params[2])") #hide println("τ : $(true_params[3]) - $(estimated_params[3])") # hide pp_to_test = TwoStateModel(20.0, 30.0, 1.0 / 30.0, NoMarks()) -test_result = MonteCarloTest(KSDistance{Exponential}, pp_to_test, h_unknown) +# n_sims=100 for the sake of this example. In real applications, this should be higher +test_result = MonteCarloTest(KSDistance{Exponential}, pp_to_test, h_unknown; n_sims=100) println("p-value for hypothesis test: $(pvalue(test_result))") # hide @@ -212,6 +213,8 @@ function PointProcesses.mark_distribution(md::LinearTimeNormal, t, h::History) return Normal(md.L + (md.H - md.L) * time_ratio, 1) end +# > It is important to think of the case where the history `h` is empty. If your custom distribution depends on past events, there should be a specific branch covering the case `isempty(h)`. + # In our case `mark_distribution` returns a `Distribution`, so `sample_mark`, `eltype` and # `densityof` are already take care of. To have all the functionality from the standard # mark distributions, we would only need to implement `fit`. diff --git a/docs/src/examples/Basics.md b/docs/src/examples/Basics.md index 3e5ad31..c565a5a 100644 --- a/docs/src/examples/Basics.md +++ b/docs/src/examples/Basics.md @@ -98,7 +98,7 @@ counts to be roughly Poisson distributed with rate λ * Δt in each bin of width ````@example Basics bin_width = 1.0 -bins = collect(h.tmin:bin_width:h.tmax) # bin edges +bins = collect((h.tmin):bin_width:(h.tmax)) # bin edges counts = fit(Histogram, h.times, bins).weights; # counts per bin nothing #hide ```` diff --git a/docs/src/examples/CustomProcesses.md b/docs/src/examples/CustomProcesses.md index 1857cc1..8104634 100644 --- a/docs/src/examples/CustomProcesses.md +++ b/docs/src/examples/CustomProcesses.md @@ -40,6 +40,8 @@ using Plots using StatsAPI using Optim +import PointProcesses: NoMarks, AbstractMarkDistribution, AbstractUnivariateProcess + struct TwoStateModel{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess λ::R Δ::R @@ -53,7 +55,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. @@ -67,7 +69,7 @@ function PointProcesses.ground_intensity(pp::TwoStateModel, t, h::History) return pp.λ else t_n = h.times[searchsortedfirst(h.times, t) - 1] # Last event time before `t` - in_high_phase = t < t_n + pp.τ + in_high_phase = t < t_n + pp.τ end return pp.λ + (in_high_phase * pp.Δ) end @@ -80,10 +82,10 @@ function PointProcesses.integrated_ground_intensity(pp::TwoStateModel, h, a, b) integral = pp.λ * (b - a) interval_events = event_times(h, a, b) if length(interval_events) > 0 - for i in 1:length(interval_events) - 1 - integral += pp.Δ * max(interval_events[i + 1] - interval_events[i], pp.τ) + for i in 1:(length(interval_events) - 1) + integral += pp.Δ * min(interval_events[i + 1] - interval_events[i], pp.τ) end - integral += pp.Δ * max(b - interval_events[end], pp.τ) + integral += pp.Δ * min(b - interval_events[end], pp.τ) end return integral end @@ -115,32 +117,33 @@ Let's plot our simulated process `h` along with its ground intensity. ````@example CustomProcesses sim_plot = plot( event_times(h), - fill(1.0, nb_events(h)) - ; line=:stem, + fill(1.0, nb_events(h)); + line=:stem, grid=false, label=false, xlabel="Time", ylabel="Events", yaxis=false, - xlim=(min_time(h), max_time(h)) + xlim=(min_time(h), max_time(h)), ) xs = LinRange(min_time(h), max_time(h), 1000) intensity_plot = plot( xs, [ground_intensity(pp, x, h) for x in xs], - ; xaxis=false, + ; + xaxis=false, xgrid=false, label=false, ylabel="Ground Intensity", - xlim=(min_time(h), max_time(h)) + xlim=(min_time(h), max_time(h)), ) plot( intensity_plot, - sim_plot, + sim_plot; title=["Two State Model Simulation" ""], - layout=grid(2, 1, heights=[0.7, 0.3]) + layout=grid(2, 1; heights=[0.7, 0.3]), ) ```` @@ -151,9 +154,7 @@ process parameters. ````@example CustomProcesses function StatsAPI.fit( - ::Type{TwoStateModel{R1,NoMarks}}, - h::History, - init_params::Vector{R2}, + ::Type{TwoStateModel{R1,NoMarks}}, h::History, init_params::Vector{R2} ) where {R1<:Real,R2<:Real} objective(params) = -logdensityof(TwoStateModel(params..., NoMarks()), h) @@ -183,11 +184,7 @@ true_params = random_params() h_unknown = simulate(TwoStateModel(true_params..., NoMarks()), 0.0, 100.0) init_params = random_params() -pp_estimated = fit( - TwoStateModel{Float64,NoMarks}, - h_unknown, - init_params, -) +pp_estimated = fit(TwoStateModel{Float64,NoMarks}, h_unknown, init_params) estimated_params = [pp_estimated.λ, pp_estimated.Δ, pp_estimated.τ] # hide println("True and estimated parameters:") # hide @@ -195,8 +192,13 @@ println("λ: $(true_params[1]) - $(estimated_params[1])") #hide println("Δ: $(true_params[2]) - $(estimated_params[2])") #hide println("τ : $(true_params[3]) - $(estimated_params[3])") # hide -pp_to_test = TwoStateModel(20.0, 30.0, 1.0/30.0, NoMarks()) -test_result = MonteCarloTest(KSDistance{Exponential}, pp_to_test, h_unknown) +pp_to_test = TwoStateModel(20.0, 30.0, 1.0 / 30.0, NoMarks()) +```` + +n_sims=100 for the sake of this example. In real applications, this should be higher + +````@example CustomProcesses +test_result = MonteCarloTest(KSDistance{Exponential}, pp_to_test, h_unknown; n_sims=100) println("p-value for hypothesis test: $(pvalue(test_result))") # hide ```` @@ -235,6 +237,8 @@ function PointProcesses.mark_distribution(md::LinearTimeNormal, t, h::History) end ```` +> It is important to think of the case where the history `h` is empty. If your custom distribution depends on past events, there should be a specific branch covering the case `isempty(h)`. + In our case `mark_distribution` returns a `Distribution`, so `sample_mark`, `eltype` and `densityof` are already take care of. To have all the functionality from the standard mark distributions, we would only need to implement `fit`. @@ -242,8 +246,12 @@ mark distributions, we would only need to implement `fit`. ````@example CustomProcesses pp = TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)) -println("TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)) = $(TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)))") -println("PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)) = $(PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)))") +println( + "TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)) = $(TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)))", +) +println( + "PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)) = $(PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)))", +) println("log_intensity(pp, m, t, h) = $(log_intensity(pp, 0.0, 0.0, h))") println("logdensityof(pp, h) = $(logdensityof(pp, h))") println("simulate(pp, 0.0, 10.0) = $(simulate(pp, 0.0, 10.0))") @@ -260,17 +268,17 @@ dimension. ````@example CustomProcesses imp = IndependentMultivariateProcess([ - TwoStateModel(1.2, 2.0, 0.3, Normal()), + TwoStateModel(1.2, 2.0, 0.3, Normal()), TwoStateModel(0.8, 2.0, 0.4, NoMarks()), PoissonProcess(2.0), # Any type of process can be added ]) 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. @@ -280,23 +288,35 @@ Lets end this tutorial with a nice plot of this process. events_plot = scatter() for d in 1:ndims(sim) - scatter!(events_plot, event_times(sim, d), fill(d, nb_events(sim, d)), label=nothing, markersize=3) + scatter!( + events_plot, + event_times(sim, d), + fill(d, nb_events(sim, d)); + label=nothing, + markersize=3, + ) end -plot!(events_plot, yaxis=false, xlabel="Time", xlim=(min_time(sim), max_time(sim)), ylim=(0.0, 3.2)) +plot!( + events_plot; + yaxis=false, + xlabel="Time", + xlim=(min_time(sim), max_time(sim)), + ylim=(0.0, 3.2), +) xs = LinRange(min_time(sim), max_time(sim), 1000) intensities_plot = plot() for d in 1:ndims(sim) - plot!(intensities_plot, xs, getindex.(λ.(xs), d), label="Dimension $d") + plot!(intensities_plot, xs, getindex.(λ.(xs), d); label="Dimension $d") end -plot!(intensities_plot, xs, sum.(λ.(xs)), label="Total ground intensity") -plot!(intensities_plot, xaxis=false, xlim=(min_time(sim), max_time(sim)), legend=:topleft) +plot!(intensities_plot, xs, sum.(λ.(xs)); label="Total ground intensity") +plot!(intensities_plot; xaxis=false, xlim=(min_time(sim), max_time(sim)), legend=:topleft) plot( intensities_plot, - events_plot, + events_plot; title="Event times and ground intensities", - layout=grid(2, 1, heights=[0.7, 0.3]), + layout=grid(2, 1; heights=[0.7, 0.3]), ) ```` diff --git a/docs/src/examples/Hawkes.md b/docs/src/examples/Hawkes.md index f9a134a..05c7237 100644 --- a/docs/src/examples/Hawkes.md +++ b/docs/src/examples/Hawkes.md @@ -59,7 +59,7 @@ dt = DateTime("12/22 3:09 pm", fmt) function parse_timestamp_min(s::AbstractString; year=2024) dt = DateTime(s, dateformat"m/d I:M p") - DateTime(year, month(dt), day(dt), hour(dt), minute(dt)) + return DateTime(year, month(dt), day(dt), hour(dt), minute(dt)) end data.TimestampDT = parse_timestamp_min.(String.(data.Timestamp)) @@ -89,7 +89,7 @@ function eventplot( xlabel="Time (minutes)", ylabel="Events", ) - scatter( + return scatter( event_times, ones(length(event_times)); markershape=:vline, @@ -135,9 +135,7 @@ tod = sort(tod_raw .+ (1:length(tod_raw)) .* 1e-6) h_day = History(tod, 0.0, 1440.0) # build a "history" on [0, 1440] minutes nbins = 96 # 96 bins = 15-minute bins pp_day = fit( - InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, - h_day, - nbins, + InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h_day, nbins ) λ_avg(u) = pp_day.intensity_function(u) / n_days @@ -180,8 +178,9 @@ Finally, we can visualize the fitted intensity function over time. ts = sort(data.t) function λ_hawkes(t::Real) - hawkes_model.μ + - sum((hawkes_model.α * exp(-hawkes_model.ω * (t - ti)) for ti in ts if ti < t); init=0.0) + return hawkes_model.μ + sum( + (hawkes_model.α * exp(-hawkes_model.ω * (t - ti)) for ti in ts if ti < t); init=0.0 + ) end u = range(0.0, maximum(ts) + 1.0; length=2000) diff --git a/docs/src/examples/Inhomogeneous.md b/docs/src/examples/Inhomogeneous.md index 489984f..c625c37 100644 --- a/docs/src/examples/Inhomogeneous.md +++ b/docs/src/examples/Inhomogeneous.md @@ -307,7 +307,7 @@ models = [ ] println("\nModel Comparison (Negative Log-Likelihood):") # hide -println("-" ^ 50) # hide +println("-"^50) # hide for (name, model) in models # hide nll = compute_nll(model, h) # hide println(" $name: ", round(nll; digits=2)) # hide diff --git a/src/mark_distributions.jl b/src/mark_distributions.jl index 9a1612b..8174fcd 100644 --- a/src/mark_distributions.jl +++ b/src/mark_distributions.jl @@ -9,6 +9,8 @@ const PointProcessMarkDistribution = Union{Distribution,AbstractMarkDistribution 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 @@ -34,7 +36,7 @@ end "The type of the marks returned by the mark distribution" function Base.eltype(md::AbstractMarkDistribution) - return eltype(mark_distribution(md, 0.0, History(0.0, 1.0))) + return typeof(sample_mark(md, 0.0, History(0.0, 1.0))) end "The likelihood of a mark `m` occurring in an event at time `t` after history `h`" diff --git a/test/bounded_point_process.jl b/test/bounded_point_process.jl index e5ece17..718b9ae 100644 --- a/test/bounded_point_process.jl +++ b/test/bounded_point_process.jl @@ -14,3 +14,20 @@ h = simulate(rng, bpp) @test all(ground_intensity_bound(bpp, 243.0, h, 1) .≈ (intensities[1], Inf)) @test integrated_ground_intensity(bpp, h, 342, 598) ≈ intensities * (598 - 342) + +struct FakePoisson <: AbstractUnivariateProcess + λ::Float64 + mark_dist::NoMarks +end + +PointProcesses.ground_intensity(fp::FakePoisson, t, h::History) = fp.λ +function PointProcesses.integrated_ground_intensity(fp::FakePoisson, h::History, a, b) + return fp.λ * (b - a) +end +PointProcesses.ground_intensity_bound(fp::FakePoisson, t, h::History) = (fp.λ, Inf) + +pp = BoundedPointProcess(FakePoisson(1.0, NoMarks()), 0.0, 10.0) +h = History([1.0, 2.0], 0.0, 10.0) + +@test simulate(pp) isa History +@test time_change(h, pp).times == h.times diff --git a/test/history.jl b/test/history.jl index 01ac65c..339dc5e 100644 --- a/test/history.jl +++ b/test/history.jl @@ -33,6 +33,9 @@ @test event_dims(h) == fill(nothing, 3) # Interface + @test_throws DomainError push!(h, 100.0, "d") + @test_throws DomainError push!(h, 100.0, "d", 2) + @test_throws DomainError push!(h, 1.0, "d", 2) push!(h, 1.7, "d") @test has_events(h, 1.5, 2.0) @@ -112,6 +115,7 @@ end @test event_dims(h_multi, 0.0, 0.3) == [1, 2] # Test push! with dimension + @test_throws DomainError push!(h_multi, 1.0, "d", nothing) push!(h_multi, 0.85, "e", 1) @test nb_events(h_multi, 1) == 3 @test event_marks(h_multi, 1) == ["a", "b", "e"] diff --git a/test/multivariate_poisson_process.jl b/test/multivariate_poisson_process.jl index 151d713..416843d 100644 --- a/test/multivariate_poisson_process.jl +++ b/test/multivariate_poisson_process.jl @@ -83,6 +83,13 @@ end @test λ_error1 < λ_error2 @test l_est > l + + pp_marks = PoissonProcess([1, 1], Normal()) + h_marks = simulate(pp_marks, 0.0, 1000.0) + pp_est_marks = fit(fill(PoissonProcess{Float64,Normal}, 2), h_marks) + @test pp_est_marks isa MultivariatePoissonProcess + @test pp_est_marks.processes[1].mark_dist isa Normal + @test logdensityof(pp_est_marks, h_marks) >= logdensityof(pp_marks, h_marks) end @testset "Time change" begin